Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
73,400 characters · 18 sections · 80 citation commands
Statistical Inference on Partially Linear Panel Model under Unobserved Linearity
\runtitle{Statistical Inference under Unobserved Linearity}
A new statistical procedure, based on a modified spline basis, is proposed to identify the linear components in the panel data model with fixed effects. Under some mild assumptions, the proposed procedure is shown to consistently estimate the underlying regression function, correctly select the linear components, and effectively conduct the statistical inference. When compared to existing methods for detection of linearity in the panel model, our approach is demonstrated to be theoretically justified as well as practically convenient. We provide a computational algorithm that implements the proposed procedure along with a path-based solution method for linearity detection, which avoids the burden of selecting the tuning parameter for the penalty term. Monte Carlo simulations are conducted to examine the finite sample performance of our proposed procedure with detailed findings that confirm our theoretical results in the paper. Applications to Aggregate Production and Environmental Kuznets Curve data also illustrate the necessity for detecting linearity in the partially linear panel model.
Keywords: Semiparametric model, Panel data, Fixed effects, Linearity detection, Penalized estimation, Partially linear regression
Keywords: C01, C14, C33
Panel models have attracted much attention from economists and econometricians, especially for their flexibility in modeling homogeneity while preserving individual-level heterogeneity. With the rapid increase in availability of panel data in the past two decades or so, panel models in both parametric and nonparametric frameworks have been well studied in the literature; see rwc00, hcl08, f17, sz15, Hsiao1997, lr15, ht08, Hsiao1997, and h14. Still, either framework, parametric or nonparametric, is not fully satisfactory in modeling panel data, as each has its own advantages and drawbacks. In light of its simplicity and interpretability, the parametric model becomes a prominent tool for panel data analysis; see bg83, dj00, kt04, s57, g64, bcs04. However, when compared to the nonparametric model, it appears to be more sensitive to model misspecification, which is often the case in empirical applications. Based on fewer model assumptions, the nonparametric model can lead to a more robust estimator, especially when dealing with relatively large panel data sets. On the other hand, with a larger dimension of input data, a purely nonparametric model is usually not preferred in empirical applications due to the infamous “Curse of Dimensionality” issue and the poor model interpretability. To address these noted drawbacks and make the best use of the apparent advantages, the partially linear panel model strikes a balance between parametric and nonparametric frameworks. For instance, hcl08 studied both nonparametric and partially linear panel models with fixed effects and proposed a kernel estimator with a corresponding linearity specification test. Combining the works by hcl08 and mbT09, ll15 proposed a two-step estimator in partially linear panel; bd02 considered the problem of estimating a partially linear fixed effects panel model with possible endogeneity and lagged dependent variables in the linear part; sz16 proposed estimation and specification testing procedures for partially linear dynamic panel model with fixed effects, with either exogenous or endogenous variables or both in the linear part and the lagged dependent variables, together with some other exogenous variables entering nonparametrically in the model. Following sz16, sz15 extended their work to the panel model with interactive fixed effects.
In practice, however, when considering the partially linear model, the researchers need to consider the following two questions: (a) which variables should be included in the model? (b) what is the functional form of each variable? Various statistical variable selection techniques, such as ll09, x09, hhw10, are available to address the first question. Nevertheless, in the context of economic modeling, one would prefer to select the dependent variables also by economic theory, as relying on purely statistical variable selection procedures may fit a model which is lacking in its economic justification and interpretability, see bvw18, bpmv15, kp09, d08. Even though the economic theory can explain which variables should be included in the model, it fails to specify the functional forms of the variables. Therefore, the second question is of more practical importance than the first one. Misspecification of the functional forms of the regressors can either (a) result in inconsistent estimation if fitting nonlinear functions by linear forms, (b) or reduce the model interpretability and estimation efficiency if the linear functions are estimated nonparametrically. Thus, correct specification of the linear components, if any, is essential to improve estimation and model explainability. However, to the best of our knowledge, all the linearity detection methods advocated in the literature on partially linear panel model, are all based on specification tests; see hcl08, sz16 and sz15. One primary drawback to this approach is that the test statistics is often difficult to construct and may be deficient in its power when the number of dependent variables is large, which may lead to incorrect model specification. Under cross-sectional data settings, zcl11 propose a smoothing-spline-type estimator which is able to estimate the underlying regression function and discover the linear regressors simultaneously. However, how to conduct valid statistical inference using the approach in zcl11 is still unknown.
The main purpose of this paper is to propose a unified statistical procedure capable of simultaneously estimating underlying regression function, detecting linear components, and conducting inference in the partially linear panel. This paper is organized as follows. In Section (ref), we mathematically formulate the linearity detection problem in the partially linear panel model. In Section (ref), we propose a penalized estimator for linearity detection, and provide {the corresponding computational} algorithm. The asymptotic properties of the proposed estimator and the corresponding linearity detection procedure are established in Section (ref) for both short and large panels. In Section (ref), we {discuss} how to determine the tuning parameters involved in the proposed procedure. Section (ref) carries out a set of Monte Carlo simulations to investigate the finite sample performance of our {proposed} method. Applications to two real-world datasets are provided in Section (ref). {Technical details and proofs of the main theorems and auxiliary results are deferred in the Appendix. Throughout this paper, we use the following notation}.
\noindentNotation: Define $\otimes$ as the tensor product operator. For positive real number $m$, let $\left\lfloor m \right\rfloor$ be the largest integer that is strictly less than $m$ and $\ceil{m}=\left\lfloor m \right\rfloor+1$. Denote $(x)_+=\max(x,0)$ for $x\in \mathbb{R}$.
Suppose that the observations $\{(Y_{it}, \mathbf{Z}_{it}), i=1,\ldots, N, t=1,\ldots, T\}$ are generated from the following model
where $Y_{it}$ is the response variable, $\mathbf{Z}_{it}=(Z_{it1}, Z_{it2}, \ldots, Z_{itp})^\top\in \mathcal{Z}:=[0, 1]^p$ are explanatory variables, both observed for individual $i$ at time period $t$, $\alpha_i^0\in \mathbb{R}$ are unobservable individual-level fixed effects, $\epsilon_{it}\in \mathbb{R}$ is unobservable errors. Assume that the unknown regression function $f_0:\mathcal{Z} \to \mathbb{R}$ has the following semiparametric expression:
where $J_{\textrm{lin}}$ is a (unknown) subset of $\{1,\ldots, p\}$ and $J_{\textrm{lin}}^c$ denotes its complement, $f_{j,0}$ for $j\in J_{\textrm{lin}}$ are linear functions and $f_{j,0}$ for $j\in J_{\textrm{lin}}^c$ are nonlinear. Our aim is to identify $J_{\textrm{lin}}$ as well as to conduct statistical inference about $f_0$ based on the observations. Without loss of generality, we may assume that $J_{\textrm{lin}}=\{1,2,\ldots, d\}$ for some nonnegative integer $d\leq p$, therefore, $f_{1,0},\ldots, f_{d,0}$ are linear and $f_{d+1,0},\ldots,f_{p,0}$ are nonlinear. For convenience, define $\{1,\ldots, d\}$ as the empty set when $d=0$.
In this section, we propose a penalized sieve estimator based on a modified spline basis, which can consistently estimate the underlying regression function $f_0$, effectively identify the linear components, and validly conduct statistical inference.
To estimate $f_0=\sum_{j=1}^pf_{j,0}$, we follow the idea of sieve estimation, i.e., estimating each $f_{j,0}$ by a linear combination of basis functions. The common basis function used in literature includes B-spline basis, wavelet basis, etc. (see c07 for an excellent review of sieve basis). However, for linearity detection purpose, the existing bases are not adequate. Thus, we will propose a modified spline space and the corresponding basis to address this issue. Given $M+1$ strictly increasing knots $\mathbf{t}_M=\{t_0, t_1, \ldots, t_M\}$ with $t_0=0, t_M=1$ and integer $r\geq 1$, define $r$-th degree Centralized Spline Space
with
being the corresponding Centralized Spline Basis. Expressed by centralized spline basis, any function in centralized spline space can be decomposed two orthogonal parts. To be more specific, for any $f=\sum_{k=1}^{r}c_k\psi_k+\sum_{k=1}^{M-1}\widetilde{c}_k\widetilde{\psi}_k\in \textrm{CSpl}(r, \mathbf{t}_M)$, we decompose $f=f_{-}+f_{\sim}$, with
which are corresponding to the linear {and nonlinear components}. It can be verified that
We begin by introduing the following function spaces
and
where $\Theta_{NT, j}=\textrm{CSpl}(r_j, \mathbf{t}_{j, M_j})$ for some integers $M_j, r_j\geq 1$ and knots $\mathbf{t}_{j,M_j}=(t_{j,0}, \ldots, t_{j, M_j})$ with $t_{j,0}=0, t_{j, M_j}=1$ for $j\in [p]$. Clearly, $\Theta_{NT}$ is a linear subspace of $\mathcal{H}_0$ and in the following it will be the sieve space to estimate the underlying regression function $f_0$. Moreover, for $g, f\in \mathcal{H}$, we {introduce the} following notation when the corresponding values exist,
{We can show that under mild assumptions}, $\langle \cdot , \cdot\rangle$ is a valid inner product on $\Theta_{NT}$ {(see Lemma} (ref) and Lemma (ref) for details). By above notation, we define a penalized objective function on $\Theta_{NT}$ as follows. For $f(\mathbf{z})=\sum_{j=1}^pf_j(z_j) \in \Theta_{NT}$ with $f_j\in \Theta_{NT,j}$, let
where $f_{j,\sim}$ is the nonlinear component of $f_j$ as defined in ((ref)), and $p_{\lambda_{NT}}$ is a given penalty function with tuning parameter $\lambda_{NT}$. The penalized estimator is defined as the minimizer of ((ref)), {namely},
There are {several possible} choices {for the functional form of the penalty term} $p_{\lambda_{NT}}$. To name a few, Ridge penalty for $p_{\lambda_{NT}}(z)=\lambda_{NT}z^2$, Lasso penalty t96 for $p_{\lambda_{NT}}(z)=\lambda_{NT}|z|$, and Smoothly Clipped Absolute Deviation (SCAD) penalty fl01 for $p_{\lambda_{NT}}$ with first order derivative
where $\kappa>2$ is some predetermined constant. In general, with larger $\lambda_{NT}$, the penalty function $p_{\lambda_{NT}}$ will be larger and thus ((ref)) will {tend} to shrink the nonlinear components $f_{j, \sim}$'s. When compared to other penalties, the solution via SCAD penalty simultaneously enjoys three desirable properties, i.e., unbiasedness, sparsity, and continuity, see fl01 for {a} detailed discussion. Therefore, throughout this paper, we will consider $p_{\lambda_{NT}}$ as SCAD penalty, and extension to other types of penalties are left as future work.
{In this section} we propose a local quadratic approximation algorithm to solve optimization problem in ((ref)). For each $j=1,\ldots, p$, let $\psi_{j, 1}, \psi_{j, 2}, \ldots, \psi_{j, r_j}, \widetilde{\psi}_{j, 1}, \ldots, \widetilde{\psi}_{j, M_j-1}$ be the centralized spline basis and for any $f_j \in \Theta_{NT,j}$, it follows that $f_j(z)=f_{j,-}(z)+f_{j, \sim}(z)$, with $f_{j,-}(z)=v_j\psi_{j,1}(z) \textrm{ and } f_{j,\sim}=u_j^\top\mathbf{\Psi}_{j,\sim}(z),$ for some $v_j \in \mathbb{R}$, $u_{j}\in \mathbb{R}^{M_j+r_j-2}$, and all $z\in [0,1]$. Here $\mathbf{\Psi}_{j,\sim}(z)=( \psi_{j, 2}(z), \ldots, \psi_{j, r_j}(z), \widetilde{\psi}_{j, 1}(z), \ldots, \widetilde{\psi}_{j, M_j-1}(z))^\top$ is a $(M_j+r_j-2)$-dimensional vector of functions. Furthermore, for each $j \in [p]$, we define vectors
and matrices
By {using the} above notation, it is not difficult to verify the following equalities,
and
Therefore, the optimization problem in ((ref)) is adapted to the optimization problem in ((ref)), which {is reduced to finding the} corresponding minimizer $v$ and $u_j$'s. {As in} fl01, {we will also use quadratic functions} to approximate the penalty terms in ((ref)). Note that
provided $u^\top \mathbf{B}_{j,\sim}^\top M_H\mathbf{B}_{j,\sim}u>0$. Therefore, if $u \approx u^0$, Taylor expansion leads to
with $D_j(u^0)=\sqrt{NT}p'_{\lambda}\left(\sqrt{\frac{1}{NT}\smash[b]{u^{0\top} \mathbf{B}_{j,\sim}^\top M_H\mathbf{B}_{j,\sim}u^0}}\right)\left(u^{0\top} \mathbf{B}_{j,\sim}^\top M_H\mathbf{B}_{j,\sim}u^0\right)^{-1/2}$ and provided $D_j(u^0)$ exists. As a consequence, if $u_j\approx u_j^0$ for all $j=1,\ldots, p$, ((ref)) can be locally approximated, up to a constant, by
From above equation, we summarize the proposed algorithm below.
In this section we present several asymptotic results concerning our proposed procedure for both short panel (fixed $T$) and large panel (diverging $T$). However, before proceeding further, we {remind the readers the Holder-smoothness notion of a function}. An univariate function $f:[0,1]\to \mathbb{R}$ is said to be $m$-smooth, if $m=r+\delta$, for some $0<\delta\leq 1$ and integer $r$ such that $f$ is $r$-times continuously differentiable and $|f^{(r)}(u)-f^{(r)}(v)|\leq c|u-v|^\delta$ for some $c>0$ and all $u, v\in [0,1]$. {Additionally, in the sequel, we use the following notation.} {We let} $q_i(\mathbf{w})$ be the density function of $\mathbf{W}_i=(\mathbf{Z}_{i1}, \ldots, \mathbf{Z}_{iT})$ and $\pi_i(\mathbf{z})$ be the density function of $\mathbf{Z}_{i1}$. For a function $g: [0,1]^k \to \mathbb{R}$, we define $\|g\|_2^2=\int g^2(\mathbf{u})d\mathbf{u}-[\int g(\mathbf{u})d\mathbf{u}]^2$ whenever the integrals exist. {Finally we set} $\mathbb{Z}=(\mathbf{W}_1, \mathbf{W}_2,\ldots, \mathbf{W}_N)$ and $\bm{\epsilon}=(\epsilon_{11},\epsilon_{12},\ldots, \epsilon_{it},\ldots, \epsilon_{NT})^\top \in \mathbb{R}^{NT}$.
The main results of this section show that the proposed penalized estimator is consistent in terms of both estimation and linearity detection. However, these results require some technical conditions, which are stated as follows.
For each $j=1,\ldots, p$, let $h_j$ be the maximal length between two successive points of knots $\mathbf{t}_{j, M_j}$, i.e., $h_j=\max_{1\leq i\leq M_j}(t_{j, i}-t_{j, i-1})$. Under Assumption (ref).(ref), it follows that $h_j \asymp M_j^{-1}$. Theorem (ref) below proves that $m_j$'s and $h_j$'s play critical roles in the rate of convergence of the proposed estimator $\widehat{f}$.
Theorem (ref) states that the rate of convergence $\widehat{f}$ consists of two parts, namely, estimation error $\sum_{j=1}^p(NTh_j)^{-1}$ and approximation error $\sum_{j=d+1}^ph_j^{2m_j}$, which coincides with standard result in h98 and h03. It should be noted that for linear components, namely $j=1,\ldots, d$, the approximation error does not involve in the $O_P$ term. On the other hand, for the nonparametric parts, the rate of convergence can benefit from balancing the estimation and the approximation errors. Specifically, if $h_j \asymp N^{-\frac{1}{2m_j+1}}$ for $j=d+1, \ldots, p$, the rate of convergence improves. It should be observed that the convergence still holds even if $h_j\asymp 1$ for $j=1,\ldots, d$ and by doing so, the rate of convergence can be further improved. The choice of $h_j$ with constant order means the number of knots $M_j$ is not diverging. Since the first $d$ components are linear, setting the corresponding $h_j$'s to be constant does not ruin the estimation consistency. However, this is usually infeasible in practice, as the prior information about the linearity of the explanatory variables is typically unavailable. Furthermore, Theorem (ref) directly shows that the global minimizer $\widehat{f}$ is consistent, while previous work about SCAD penalized regression only establishes the existence of a consistent local minimizer, e.g., see fl01 and x09.
Theorem (ref) only addresses the issue for estimation, which is not adequate to distinguish the linear components from the nonlinear ones. While with appropriate choice of tuning parameter $\lambda_{NT}$, Theorem (ref) below proves that the estimator $\widehat{f}$ will automatically recover the linearity in the underlying regression function $f_0$.
The tuning parameter $\lambda_{NT}$ in Theorem (ref) (unlike in Theorem (ref)) can neither be too large nor too small. With suitable choices of $\lambda_{NT}$ and $h_j$'s, the proposed estimator $\widehat{f}=\sum_{j=1}^p\widehat{f}_j$ will automatically and correctly specify the linear and nonlinear forms with probability approaching one. Since in Theorem (ref), the tuning parameters $h_j$'s and $\lambda_{NT}$ play important roles in selection consistency, a fundamental issue in practice is the choice of these parameters. The discussion of this issue is deferred to Section (ref).
{In this section we define the} solution path of $\widehat{f}$ and provide its theoretical properties and practical implications. {For fixed} knots $\mathbf{t}_{j, M_j}$'s and {the} tuning parameters $k_j$'s and $h_j$'s, one can obtain a sequence of estimators $\widehat{f}$ by using a sequence of increasing $\lambda_{NT}$'s and these estimators forms a solution path. For sufficiently large $\lambda_{NT}$, all the nonlinear components $\widehat{f}_{j,\sim}$'s will vanish and result in a model consisting of all linear components. On the other hand, when $\lambda_{NT}$ is close to zero, all the $\widehat{f}_j$'s will be nonlinear. {Consequently, we may} obtain $p+1$ different models in the solution path by increasing $\lambda_{NT}$ from zero to infinity. The following corollary is a direct consequence of Theorem (ref).
Corollary (ref) indicates that the solution path is consistent in the sense that, one in the $p+1$ models will correctly identify both the linear and the nonlinear parts. Notice that for linearity detection problem, {one essentially needs} to identify the correct model out of $2^p$ candidates. Another {immediate} implication from Corollary (ref) is that in practice, any model selection method, e.g., Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) criteria, based on these $p+1$ models is valid, reliable and is equivalent to that based on $2^p$ models, which is a significant reduction on model complexity.
In this section, we will present the limit distribution of proposed estimator $\widehat{f}$. To proceed further, recall that $\mathbf{\Psi}_{j,\sim}(z)=( \psi_{j, 2}(z), \ldots, \psi_{j, r_j}(z), \widetilde{\psi}_{j, 1}(z), \ldots, \widetilde{\psi}_{j, M_j-1}(z))^\top$ is the basis of $\Theta_{NT,j,\sim}$, for $j=1,\ldots, p$. We further define $\mathbf{\Psi}_{j,-}(z)=\psi_{j,1}(z)=z-1/2$, $\mathbf{\Psi}_j(z)=(\mathbf{\Psi}_{j,-}(z), \mathbf{\Psi}{j,\sim}^\top(z))^\top$, and $\mathbf{\Psi}^0(\mathbf{z})=(\mathbf{\Psi}_{1,-}(z_1),\ldots, \mathbf{\Psi}_{d,-}(z_d), \mathbf{\Psi}_{d+1}^\top(z_{d+1}),\ldots, \mathbf{\Psi}^\top_{p}(z_p))^\top$. By this definition, we know $\mathbf{\Psi}_{j}(z)$ is the basis of $\Theta_{NT,j}$. If we define the space of correctly specified model
then $\mathbf{\Psi}^0(z)$ will be its basis. By Theorems (ref), it follows that $\widehat{f}\in \Theta_{NT}^0$ with probability approaching one, and thus we have the following expression for the proposed estimator:
Therefore, it is natural for us to study the asymptotic distributions of $\widehat{\beta}_j$'s and $\widehat{f}_j(z_{j,0})$'s with $z_{d+1,0}, \ldots, z_{p, 0}\in [0, 1]$ being some fixed constants. We consider following elements in $\Theta_{NT}^0$:
with
It can be shown that, for any $f(\mathbf{z})=\sum_{j=1}^d\beta_j(z_j-1/2)+\sum_{j=d+1}^p{f}_{j}(z_j)\in \Theta_{NT}^0$ with $f_j\in \Theta_{NT,j}$, $j=d+1,\ldots, p$, the following equality holds:
If we define linear functionals from $\Theta_{NT}^0$ to $\mathbb{R}$ such that $\mathcal{L}_j(f)=\beta_{j}$ for $j=1,\ldots, d$ and $\mathcal{L}_j(f)=f_j(z_{j,0})$ for $j=d+1,\ldots, p$, then $v_{NT,j}^*$'s are essentially the Riesz representatives of $\mathcal{L}_j$'s.
In order to establish the asymptotic distribution, more regular assumptions on the error terms $\epsilon_{it}$'s are needed. Thus, in the following, we define the standard deviation inner product and norm in $\Theta_{NT}$, which contains the information of $\bm{\epsilon}$. For $g, f \in \Theta_{NT}$, we define
where $\mathbf{g}=(g(\mathbf{Z}_{11}), g(\mathbf{Z}_{12}),\ldots, g(\mathbf{Z}_{it}), g(\mathbf{Z}_{NT}))^\top, \mathbf{f}=(f(\mathbf{Z}_{11}), f(\mathbf{Z}_{12}),\ldots, f(\mathbf{Z}_{it}), f(\mathbf{Z}_{NT}))^\top \in \mathbb{R}^{NT}$. In addition, denoting $\bm{\epsilon}_i=(\epsilon_{i1},\ldots, \epsilon_{iT})$ for $i \in [N]$, we propose Assumption (ref) on the error terms $\bm{\epsilon}_i$'s and $v_{NT,j}^*$ for statistical inference.
For presentation purpose, we choose $h_1=h_2=\ldots=h_p=h$ and define $m_*=\min_{d+1\leq j \leq p}m_j$. Theorem (ref) below states that, with suitable choice of $h$ and $\lambda_{NT}$, we can obtain the asymptotic distribution of $(\widehat{\beta}_1, \ldots, \widehat{\beta}_d, \widehat{f}_{d+1}(z_{j,0}),\ldots, \widehat{f}_{p}(z_{j,0}))$.
Theorem (ref) establishes the joint asymptotic distribution of both the linear and nonlinear components of $\widehat{f}$, which includes estimators with different rate of convergence. sc13aos, cs15 and dl18 also established similar joint asymptotic results in partially linear model. However, compared with their results, Theorem (ref) does not require the prior knowledge of linearity. The constant $m_*$ is the smallest degree of smoothness among all the $f_{j,0}$'s, which represents the effective smoothness of $f_0$. From Theorem (ref), a necessary condition is $m_*>1.5$ for short panel and $m_*>2$ for large panel, which requires the underlying regression function needs to be enough smooth. If one is of more interest in the marginal distribution of each $\widehat{f}_j$, Theorem (ref) below establishes the limit distribution of $\widehat{f}_{j}(z_{j,0})$ without Assumption (ref), where $z_{j,0}\in [0,1]$ is fixed constant for $j=1,\ldots,p$.
The choice of homogeneous $h_j$'s in Theorems (ref) and (ref) is not only simple for presentation, but also it is practically convenient. As discussed in Section (ref), homogeneous $h_j$'s will reduce the complexity of tuning parameter selection. For theoretical interest, we include the case of heterogeneous $h_j$'s in Appendix.
{In this section we discuss} how to determine tuning parameters $r_j$'s, $h_j$'s and $\lambda_{NT}$. {Motivated by two different objectives,} we propose two distinct strategies to select $\lambda_{NT}$ for estimation and {for} linearity detection. For convenience, we simply choose each of the knots $\mathbf{t}_{j, M_j}$, to be an uniform partition of $[0, 1]$ in practice when $h_j$'s are determined.
Before proceeding further, we formally define $k$-fold cross validation procedure in the framework of panel data. Given positive integer $k\geq 2$, $N$ individuals are randomly separated into $k$ disjointed groups and let $I_1, \ldots, I_k$ be the corresponding sets of indexes with $N_1,\ldots, N_k$ elements, respectively. By this notation, it follows that $I_1,\ldots, I_k$ is a partition of $\{1,\ldots, N\}$. Moreover, we denote $I_s^c$ as the compliment of $I_s$ for $s=1,\ldots, k$ and $\theta=(r_1,\ldots, r_p, h_1,\ldots, h_p, \lambda_{NT})$ as the tuning parameters. Given $\theta$, we set $\widehat{f}_{I_s^c,\theta}$ to be the penalized estimator based on observations $\{(\mathbf{Y}_i, \mathbf{W}_i), i\in I_s^c\}$ and tuning parameter $\theta$. The cross validation value is defined as follows,
Based on ((ref)), the optimal tuning parameter ${\theta}_{\textrm{opt}}$ is defined as the minimizer of $\textrm{CV}(\theta)$ among several candidates, i.e.,
where the minimum is taken over some pre-specified values. The procedure in ((ref)) is called $k$-fold cross validation, which {provides a} powerful tool {for choosing the} tuning parameters with solid theoretical {justifications} , see a91, h14b and y07. Other methods for empirical {choices} of $k_j$'s and $h_j$'s in the framework of sieve estimator can be found in ho14 and cc18.
In sieve estimation, the choices of $k_j$'s and $h_j$'s {play} essential roles in the estimation accuracy. For example, Assumption (ref).(ref) specifies lower bounds on $r_j$'s, while Theorem (ref) implies that if $h_j\asymp (NT)^{-\frac{1}{2m_j+1}}$ for $j=d+1,\ldots, p$, the rate of convergence for $\widehat{f}$ will be improved and a more accurate estimation is obtained. The procedure in ((ref)) to determine $k_j$'s and $h_j$'s also needs to specify $\lambda_{NT}$ simultaneously, which is inconvenient in practice {as it involves} too many parameters. To address this concern, we use cross validation criterion based on non-penalized estimator for the choices of $k_j$'s and $h_j$'s. To be more specific, {their optimal choices,} $k_{j,\textrm{opt}}$'s and $h_{j,\textrm{opt}}$'s are defined as follows,
The cross validation procedure in ((ref)) is motivated by Theorem (ref), since the rate of convergence is the same regardless of the penalty.
After selecting $k_j$'s and $h_j$'s, we may choose $\lambda_{NT}$ in three distinct ways different purposes.
For estimation, a similar procedure as ((ref)) is recommended. Specifically, given pre-determined $k_j$'s and $h_j$'s, the optimal $\lambda_{NT, \textrm{opt}}$ is selected as follows,
with the minimum taken over some pre-determined candidates of $\lambda_{NT}$.
However, for linearity detection, we propose a practically convenient approach to select a model based on solution path without determining $\lambda_{NT}$. By Corollary (ref), the solution path will select $p+1$ models, {in which one} correctly identifies all the linear components. Therefore, it is natural for us to perform model selection among these $p+1$ candidates. For $\nu=1,\ldots, p+1$, let $J_\nu \subset \{1,\ldots, p\}$ be the set of indexes of linear components selected by $\nu$-th model along the solution path. We further define the function space
and non-penalized estimator for $k$-fold cross validation
Similar to ((ref)), we propose following procedure to identify linearity based on $k$-fold cross validation,
The procedure in ((ref)) is completely data-driven without {the need to choose} $\lambda_{NT}$. Based on the solution path, other information {criteria, such as such AIC or BIC,} also can be applied to conduct model selection, see h14b and b06.
To conduct valid statistical inference, $\lambda_{NT}$ is selected based on the solution path and $\widehat{J}_{CV}$ defined in ((ref)). First, based on the solution path, we find the values of $\lambda_{NT}$ resulting in the model with indexes of linear components being $\widehat{J}_{CV}$. Then the turning parameter $\lambda_{NT,\textrm{inf}}$ is chosen to be the smallest one among these values.
To evaluate the finite sample performance of {the proposed estimation and selection procedure}, we consider the following data generating process,
The {functional forms of the underlying} regression functions are specified as follows
with the first two $f_j$'s being linear and the last two being nonlinear functions whose degree of nonlinearity is controlled by a factor of $r$. The function $\beta_{6,9}(z)$ is density of the beta distribution with parameters $(6, 9)$. The fixed effect $\alpha_i$'s and the idiosyncratic error $\epsilon_{it}$'s are i.i.d standard normal random variables across $i$ and $t$. The regressors $z_{itj}, j=1,2,3,4$ are generated as follows, (a)$\{u_{itj}, i\in [N], t\in [T], j\in [p]\}$ are i.i.d uniform random variables on $[0, 1]$; (b) $z_{it1}=u_{it1}+\alpha_i$; (c) for $j=2,3$, $z_{itj}=u_{itj}+\delta_i$, with $\delta_i$'s being i.i.d standard normal random variables; (d) $z_{it4}=u_{it4}$. For the sample size and degree of nonlinearity, we consider all combinations of $(N, T, r)$ with $N=(50, 100, 200)$, $T=(3, 10, 50)$ and $r=(0.01, 0.1, 0.2, 0.5, 1)$, which include both short and large panel settings with weak and strong nonlinearity. The number of replication is set to be $R=500$. For convenience, we choose the degree of polynomial spline $k_j=3$ for $j\in [4]$ and set $h_1=h_2=h_3=h_4=h$ with $h^{-1}$ determined by $5$-fold cross validation among $\{\ceil[\big]{c(NT)^{1/4}}+2, c=0.3, 0.4,\ldots, 2\}$.
In the following, we {consider} three numerical experiments to study the finite sample performance of the proposed procedure.
Experiment 1: For the estimation, the estimator $\widehat{f}(\mathbf{z})=\sum_{j=1}^T\widehat{f}_j(z_j)$ is evaluated using the root mean squared error (RMSE) defined as
An integration term is added in above equation, since $\widehat{f}_j$ essentially is the estimator of $f_j-\int_0^1 f_j(z)dz$ (see Assumption (ref).(ref)). A sequence of $\lambda_{NT}$'s in $[0, 1]$ are used in the experiment to obtain the estimator. In particular, $\lambda_{NT}=0$ results in a non-penalized estimator.
Experiment 2: For linearity detection, we generate the solution path along a sequence of $\lambda_{NT}$'s with $\log(\lambda_{NT})=\{-6,-5.9,\ldots, 0.9, 1\}$. Four different proportions are calculated among 500 replications, namely, proportion of solution path containing the correct model and proportions of correct linearity detection from solution path based on $5$-fold cross validation score (CV), AIC and BIC, respectively.
Experiment 3: To study the asymptotic normality of proposal estimator, we consider the setting with $r=1$. For $z_0=(0, 0.25, 0.5, 0.75, 1)$, we construct the point wise confidence intervals for $f_{j,0}(z_0)$ based on Theorem (ref). We calculate the percentages of the ground truth $f_{j,0}(z_0)$ falling in the 95% confidence intervals.
Figure (ref) reports the RMSE of the proposed estimator with different sample sizes and degrees of nonlinearity. Some interesting findings can be observed in Figure (ref). Firstly when varying $\lambda_{NT}$, for cases that $r=0.5, 1$ with strong nonlinearity, RMSE decreases and then increases, while for cases $r=0.01, 0.1$ with weak nonlinearity, RMSE decreases and then stays the same regardless of the sample size. For the case with moderately strong linearity, namely $r=0.2$, with small sample size $N=50, T=3$, RMSE follows a similar pattern as that of weak linearity cases, while with other sample sizes, there is a decrease on RMSE when $\lambda_{NT}$ varying from $0$ to $0.17$ and followed by a slight growth when $\lambda_{NT}$ increasing from $0.17$ to $0.2$. With larger $\lambda_{NT}$, the RMSE remains the same. Secondly, with larger $\lambda_{NT}$, RMSE stays at the same level regardless of $\lambda_{NT}$ in each case except $r=1$. Thirdly, for $\lambda_{NT}\geq 0.4$, the RMSE increases as the degree of nonlinearity becomes larger. Finally, with an appropriate choice of $\lambda_{NT}$, a penalized estimator can outperform nonpenalized estimator in terms of RMSE. For linearity detection, Figure (ref) reveals that with stronger nonlinearity, all the procedures are more likely to perform a correct linearity detection. In particular, when $r=1$, the solution path will contain the correct model in all replications except when sample size is small, $N=50, T=3$, which confirms the validity of Corollary (ref). Moreover, among three criteria for model selection, BIC and CV score can effectively choose the true model when $r$ is large, while AIC and CV score work better for small $r$. Figure (ref) reports the coverage rates of the 95% confidence intervals for $f_{j,0}(z_0)$. It is worth mentioning that, for the linear components $f_{1,0}$ and $f_{2,0}$, the coverage rates are almost $100\%$ when $z_0=1/2$. This is due to ((ref)) that $\widehat{f}_j(1/2)=0$ if $\widehat{f}_j$ is estimated as a linear function. In general, the coverage rate will approach to 95% when $N$ becomes larger or both $N$ and $T$ become larger.
In this section, we apply our linearity detection procedure to Aggregate Production data, which is extracted from version 9.0 of the Penn World Table. We keep a balanced panel dataset for $48$ countries across the world for the period 1950-2014. Following gks16, we consider following regression model,
where $y_{it}, k_{it}, l_{it}$ are the real log GDP, capital stock, and number of people engaged of the $i$-th country at time $t$, respectively. Besides, $\textrm{pub}_{it}$ is government/public expenditure, defined as the government spending, and $\textrm{xm}_{it}$ is net trade openness, which equals exports minus imports of merchandise. Using the same criteria as in the simulation study, we choose $h_1=h_2=h_3=h_4=0.2$ by cross-validation. The tuning parameter $\lambda_{NT}$ is selected such that $\log(\lambda_{NT})$ increases from $-4$ to $-3$ with an increment $0.05$. Firstly, the solution path in Figure (ref) indicates five candidate models can be obtained when $\lambda_{NT}$ increasing with models being summarized in Table (ref). Moreover, according to Table (ref), among these five candidates, the model with all linear $f_j$'s is the preferable based on CV. The estimated coefficients of explanatory variables are provided in Table (ref), from which we can see the model is highly significant. Non-penalized estimators of $f_j$'s are constructed and corresponding fitted curves are provided in Figure (ref). The fitted curves of the non-penalized estimators preserve linear patterns if one only looks at the interval between two dashed vertical lines, which is the $2.5\%$ to $97.5\%$ percentile range of the corresponding regressor. This information contained in Figure (ref) coincides with the findings based on our proposed linearity detection procedure.
In the second application, we estimate the environmental Kuznets curve (EKC), which is also studied in a07, ap09 and lqs16. Following lqs16, we consider the following nonparametric model:
where $y_{it}$ represents the per capita $\textrm{CO}_2$ emission of country $i$ in year $t$, $\textrm{e}_{it}$ is per capita energy consumption, $\textrm{gdp}_{it}$ stands for the per capita GDP, and $\textrm{trade}_{it}$ is the per capta trade. All the variables are taken logarithm and all the explanatory variables are scaled to $[0, 1]$. The data is obtained from World Bank Development Indicators and we keep a balanced panel for $N=89$ and $T=40$ after eliminating missing values.
From the solution path in Figure (ref), we extract 4 submodels and calculate their CV scores, which are summarized in Table (ref). Based on the CV score, the selected linear explanatory variables are trade and e, while the variable gpd will be treated as nonlinear. Based on selected, model, we further estimate the linear and nonlinear components, and the results are reported in Table (ref) and Figure (ref). Table (ref) shows the coefficients of e and trade are both highly significant. Meanwhile, Figure (ref) indicates that as gdp increasing, its effect on $\textrm{CO}_2$ emission will increases first and then begin to fall, which coincides with common hypothesis that the relationship between income and the emission of chemicals like sulfur dioxide ($\textrm{SO}_2$) and carbon dioxide ($\textrm{CO}_2$) or the natural resource usage has an inverted U-shape, see lqs16. Finally, we present the nonpenalized estimation curves of the explanatory variables in Figure (ref). If only screening the fitted curves in $2.5\%$ to $97.5\%$ pencentile range of the regressors in Figure (ref), we can draw the same conclusion that e and trade are linear, while gdp is nonlinearly correlated with CO$_2$ emission.
\setcounter{subsection}{0} \setcounter{subsubsection}{0} \setcounter{equation}{0} \setcounter{lemma}{0} \setcounter{proposition}{0}