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,216 characters · 11 sections · 61 citation commands
A projection based approach for interactive fixed effects panel data models
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 {
} \fi
{\it Keywords:} cross-sectional dependence, semiparametric factor models, principal components, sieve approximation, large panels.
\spacingset{1.9}
In this paper, we consider the consistent estimation of the following panel data model with interactive fixed effects, often referred to as a factor structure model. These models have gained significant attention over the past decade due to their ability to control latent variables and account for cross-sectional dependence. The model we consider takes the following form,
where $y_{it}$ is the response variable of individual $i$ at time $t$, $\mathbf{X}_{it}$ is a $Q$-dimensional vector of covariates, and $\boldsymbol{\boldsymbol{\beta}}$ is the $Q$-dimensional vector of parameters to be estimated. The error term $v_{it}$ has two components: (i) a latent factor structure $\boldsymbol{\boldsymbol \lambda}_i^{\top}\mathbf{f}_t$, where $\boldsymbol{\boldsymbol \lambda}_i=(\lambda_{i1},\ldots,\lambda_{iK})^{\top}$ represents individual specific factor loadings and $\mathbf{f}_t=(f_{t1},\ldots,f_{tK})^{\top}$ denotes common unobserved factors; (ii) an idiosyncratic error term, $u_{it}$, which is assumed to be zero mean and independent of the covariates and the factor structure. The number of factors $K$ is finite and does not depend on the size of the cross-section $N$ or the time dimension $T$.
A critical challenge of model ((ref)) is the potential correlation between covariates $\mathbf{X}_{it}$ and latent factors, $\mathbf{f}_t$, or the individual factor loadings, $\boldsymbol{\boldsymbol \lambda}_i$. This correlation introduces cross-sectional dependence (CSD) in the structural model, and standard estimation techniques for panel data models can lead to inconsistent estimates and spurious inferences of the parameters of interest because they do not account for the latent structure embedded in the error term $v_{it}$. In this context, this paper aims to provide an innovative estimation procedure to deal with this endogeneity issue and obtain consistent estimates of the slope coefficient vector $\boldsymbol{\boldsymbol{\beta}}$.
Addressing this issue is crucial in many areas of economics, where some of the regressors are decision variables that are influenced by unobserved individual heterogeneities. For instance, in macroeconomics, cross-country growth studies often seek to understand the determinants of economic growth, $y_{it}$, using observable variables, $\mathbf{X}_{it}$, such as physical capital investment, population growth, and trade openness. However, global shocks, such as financial crises or fluctuations in world oil prices affect all countries, introduce common latent factors that affect all countries simultaneously through trade and financial links and must be taken into account to obtain valid estimates Chudik_Mohaddes_Pesaran_Raissi2017, lu2016shrinkage. Another relevant example is found in financial econometrics, where excess returns $y_{it}$ for asset $i$ at time $t$ are often modeled using observable factors such as Fama-French factors (e.g., small market capitalization and book-to-market ratios), dividend yields and payout ratios, among others. However, asset returns are also influenced by unobserved global financial conditions, such as the state of the global credit market, which induce a strong CSD across assets Bernanke-Booivin-Eliasz2005, Fan-Ke-Liao2021.
To address this source of bias in large panels (where both $N$ and $T$ are large), two main strands have emerged in the literature: one focuses on controlling for common factors, $\mathbf{f}_t$, while the other focuses on modeling factor loadings, $\boldsymbol{\boldsymbol \lambda}_i$. The literature on controlling for common factors is extensive. A widely used methodology is the Common Correlated Effects (CCE) estimator introduced by Pesaran2006, which approximates unobserved factors through linear combinations of cross-sectional averages of both dependent and explanatory variables. While this approach offers computational simplicity and broad applicability, its reliance on asymptotic properties may lead to biased estimates in finite samples (see Westerlund-Urbain2013 for a deeper discussion). An alternative methodology, developed by Bai2009 and further examined by Moon-Weidner2015, Moon-Weidner2017, adopts a fundamentally different approach. Their method employs the principal component (PC) approach to directly estimate the unobserved factors and later compute consistent estimates for $\boldsymbol{\beta}$ solving a non-convex optimization problem. Although this PC-based approach provides an elegant solution, it presents practical challenges due to its computational intensity and sensitivity to the number of unobserved factors. After these original papers, a large body of literature emerged by extending these estimation procedures to more general settings. See Sarafidis_Wansbeek2012, chudik_pesaran2015, or Bai-Wang2016 for excellent surveys, and Westerlund-Urbain2015 for a comparison analysis between the CCE and PC estimation procedures.
The second strand of the literature focuses on modeling the factor loadings $\boldsymbol{\boldsymbol \lambda}_i$ as a function of observed time-invariant covariates $\mathbf{Z}_{i}$. Connor-Linton2007, Connor-Hagmann-Linton2012, Zhang_Zhou_Wang2021, and Cheng_Dong_Gao_Linton2024 propose modeling the factor loadings as $\lambda_{ik}=g_k(\mathbf{Z}_{i})$ for some unknown function $g_k(\cdot)$ and $k=1,\ldots,K$. However, this modeling imposes a restrictive assumption that the entire variation in $\boldsymbol{\boldsymbol \lambda}_i$ must be explained by $\mathbf{Z}_{i}$, which increases the risk of model misspecification. To allow for partial flexibility and mitigate the risk of misspecification, Fan-Liao-Wang2016 propose a more flexible framework that decomposes the factor loadings into a systematic component and an idiosyncratic error:
where $\mathbf{Z}_i$ is a $D$-dimensional vector of additional covariates of individual characteristics, $\mathbf g(\mathbf{Z}_i)=(g_1(\mathbf{Z}_i),\ldots,g_K(\mathbf{Z}_{i}))^{\top}$ is a $K$-dimensional vector of unknown functions, and $\boldsymbol{\boldsymbol \gamma}_i=(\gamma_{i1},\ldots, \gamma_{iK})^{\top}$ is a $K$-dimensional vector of errors that reflects the part of $\boldsymbol{\boldsymbol \lambda}_i$ that cannot be explained by $\mathbf{Z}_{i}$. Throughout the paper, it is assumed that $\{\boldsymbol{\boldsymbol \gamma}_i\}_{i\leq N}$ is zero mean and independent of $\{\mathbf{Z}_{i}\}_{i\leq N}$.
Building on this literature, this paper introduces a novel estimation methodology, very easy to implement from the empirical point of view, that leads to consistent estimators of $\boldsymbol{\beta}$ in ((ref)) when $\mathbf{X}_{it}$ are correlated with $\boldsymbol{\boldsymbol \lambda}_i$ and/or $\mathbf{f}_{t}$. More precisely, inspired by Fan-Liao-Wang2016 we propose to introduce the relationship in ((ref)) into a panel data framework, as Eq. ((ref)). Hence, plugging ((ref)) into ((ref)) and rearranging terms, we get the following regression model
This new framework addresses the bias caused by the term $\mathbf{g}(\mathbf{Z}_{i})^{\top}\mathbf{f}_t + \boldsymbol{\boldsymbol \gamma}_i^{\top}\mathbf{f}_t$, which introduces endogeneity in the ordinary least squares (OLS) estimator. Hence, the innovation of this paper lies in projecting the data onto the subspace generated by sieve basis functions of the covariates $\mathbf{Z}_{i}$. This orthogonal projection intends to “project away” the unobserved factor loadings to eliminate the bias asymptotically and obtain consistent and asymptotically normal estimators of the $\boldsymbol{\boldsymbol{\beta}}$'s in ((ref)) without the need for computationally intensive procedures.
Our model setup is closely related to the one in Zhang_Zhou_Wang2021, however, we want to point out crucial distinctions. First and foremost, the main issue of interest in the above paper is efficiency and they propose a GLS-type estimator that under broadly general conditions is oracle efficient. It is important to note that their asymptotic results require the consistency of the pooled OLS estimator in a first step which is not the case in our model setup. Second, Zhang_Zhou_Wang2021 assume that the factor loadings are fully explained by $\mathbf{Z}_i$, i.e., $\boldsymbol{\gamma}_i = 0$. Unfortunately, in the presence of error factor loadings, i.e., $\boldsymbol{\gamma}_i\ne 0$, the statistical properties of standard estimators for $\boldsymbol{\boldsymbol{\beta}}$ remain unclear. Therefore, it is of interest to derive a new estimator to obtain consistency and asymptotic rates.
In this framework, the estimation procedure proposed in this paper offers several advantages. First, it is based on a simple OLS framework, which avoids the complexities of iterative procedures such as the PC method with unknown convergence properties that may be computationally intensive. Second, it does not require prior knowledge of the number of common factors and does not require knowledge or assumptions about them, making it robust to various specifications. Further, the underlying limiting distribution is centered at zero. Therefore, the proposed estimation techniques circumvent the omitted variable bias problem as in Bai2009. Finally, the proposed estimator reaches the semiparametric efficiency bound under certain conditions.
We want to highlight that the asymptotic results of the proposed estimator hold irrespectively of the variance of the error factor loadings being zero, close to zero, or much larger than zero. However, there exists a discontinuity in the limiting distribution when this variance is close to zero. Then, the usual “plug-in” approaches would lead to valid but overly conservative inference. Similarly, ignoring the idiosyncratic part leads to invalid inference in the case of persistent variance. To achieve uniformly valid but non-conservative inference, we resort to the cross-sectional bootstrap originally proposed by kapetanios2008bootstrap. In this way, by stacking the time observations we are able to mimic the asymptotic distribution and conduct uniformly valid inference even in the presence of this type of discontinuities as shown by liao2018uniform and FernndezVal2022DynamicHD. The issue of uniformity is an important topic for modeling panel data. lu2023uniform consider a model with two-dimensional heterogeneity of varying degrees in the slope parameters and are interested in uniformly valid inference. kock2016oracle and kock2019uniform study uniform inference in high-dimensional panel regression contexts. menzel2021bootstrap showed that uniform non-conservative inference is impossible under general dependence in more than one dimension. The novelty of our bootstrap procedure is that we resample cross-sectional units after projecting the data, i.e., partialing-out the modeled part of the factor loadings.
The rest of the paper is organized as follows. In Section 2, we derive our projection-based interactive fixed effects estimator. Section 3 states our assumptions and studies the asymptotic properties of the proposed estimators. In Section 4 we validate the theoretical results in a simulation study. In Section 5 we apply our method to the identification of the determinants of economic growth. lu2016shrinkage argued that the GDP growth rates per capita might not only be determined by observed factors but might also be influenced by latent factors or shocks. Our projection-based interactive fixed effects estimator is well suited for such a setting. All proofs of the asymptotic results and further Monte Carlo results are relegated to a Supplementary Material document.
To nonparametrically estimate the unknown function $g_k(\cdot)$ without curse of dimensionality, it will be assumed that for each $k$, where $k=1,\ldots,K$, $g_k(\cdot)$ is an additive function of the form
For each $k$ and $d$, the additive component $g_{kd}(\cdot)$ can be approximated by the sieve method. We define $\left\{\phi_1(Z_{i,d}), \ldots, \phi_{J_N}(Z_{i,d})\right\}$ as a set of basis functions (i.e., splines, Fourier series, wavelets), which spans a dense linear space of the functional space for $g_{kd}(\cdot)$. Then,
where, for $j=1,\ldots,J_N$, $\phi_{j}(\cdot)$'s are the sieve basis functions, $b_{j,kd}$'s are the sieve coefficients of the $d$th additive component of $g_k(\mathbf{Z}_{i})$ corresponding to the $k$th factor loading, $R_{kd}(\cdot)$ is a “remainder function” that represents the approximation error, and $J_N$ denotes the number of sieve terms which grows slowly as $N\rightarrow\infty$.
As it is well-known in the literature, under some regularity condition of the functional class, the approximation functions $\phi_{j}(\cdot)$ have the property that, as $J_N$ grows, there is a linear combination of $\phi_{j}(\cdot)$ that can approximate $g_k(\cdot)$ arbitrarily well in the sense that the approximation error can be made arbitrarily small. Therefore, for a given $d=1,\ldots,D$, the basic assumption for sieve approximation is that $\sup_z|R_{kd}(z)|\rightarrow0$, as $J_N\rightarrow\infty$. In practice, an optimal choice for the smoothing parameter $J_N$ can be based on cross-validation.
For the sake of simplicity, we take the same basis functions in ((ref)) and, for each $k\leq K$, $d \le D$ and $i\leq N$, let us define
so the above equation can be rewritten as
Let $\mathbf{Z} = \left(\mathbf{Z}^{\top}_{1},\ldots,\mathbf{Z}^{\top}_{N}\right)$ be an $N\times D$ matrix whose $i$th element is a $D$-dimensional vector of random variables as $ \mathbf{Z}_{i}=(Z_{i1},\ldots,Z_{iD})^{\top}$ and denote by $\mathbf{\mathcal{Z}}$ its support. Let also $\mathbf{G}(\mathbf{Z})$ be an $N\times K$ matrix of unknown functions, $g_k(\mathbf{Z}_{i})$, $\Phi(\mathbf{Z})=(\phi(\mathbf{Z}_{1}),\ldots,\phi(\mathbf{Z}_{N}))^{\top}$ be an $N\times J_ND$ matrix of basis functions, $\mathbf B=(\mathbf b_1,\ldots,\mathbf b_K)$ be a $J_ND\times K$ matrix of sieve coefficients, and $\mathbf R(\mathbf{Z})$ be an $N\times K$ matrix with the $(i,k)$th element $\sum^D_{d=1}R_{kd}(\mathbf{Z}_{i,d})$. By considering ((ref)) in matrix form, we obtain
and substituting ((ref)) into the matrix form of ((ref)) leads to
where the residual term consists of two parts: the sieve approximation error, $\mathbf R(\mathbf{Z})\mathbf{f}_t$, and the error term, $\mathbf v_t$. $\mathbf{v}_t$ is an $N\times 1$ vector such as $\mathbf{v}_t=\boldsymbol \Gamma \mathbf{f}_t+ \mathbf u_t$, where $\boldsymbol\Gamma=(\boldsymbol{\gamma}_1,\ldots,\boldsymbol{\gamma}_N)^{\top}$ is an $N\times K$ matrix of unknown loading coefficients and $\mathbf u_t = \left(u_{1t},\ldots,u_{Nt}\right)^{\top}$ is an $N\times 1$ vector of idiosyncratic errors. Finally, $\mathbf{X}_t$ is an $N \times Q$ matrix of covariates.
To obtain consistent estimators of $\boldsymbol{\beta}$ in ((ref)) we propose a transformation that removes $\Phi(\mathbf{Z})\mathbf B\mathbf{f}_t$ and accounts for the error term, $\boldsymbol \Gamma \mathbf{f}_t$. A natural choice to remove $\Phi(\mathbf{Z})\mathbf B\mathbf{f}_t$ in ((ref)) is to define the following projection matrix
Premultiplying both sides of ((ref)) by $\mathbf{P}_{\Phi}(\mathbf{Z})$ and assuming that $(NT)^{-1}\sum_{t=1}^T \mathbf{X}_t^{\top} \left[\mathbf{I}_N- \mathbf P_{\Phi}(\mathbf{Z})\right]\mathbf{X}_t$ is non-singular, the following estimator for $\boldsymbol{\boldsymbol{\boldsymbol{\beta}}}$ is obtained,
In this section, we analyze the main asymptotic properties of the estimator. Firstly, we introduce some notation, definitions, and assumptions that will be necessary to derive the main results of this paper. Later, we present the main large sample properties of these estimators. All proofs of these results are relegated to a Supplementary Material document.
Let $n = NT$. For two positive number sequences $(a_n)$ and $(b_n)$, we say $a_n={\mathcal{O}}(b_n)$ or $a_n\lesssim b_n$ (resp. $a_n\asymp b_n$) if there exists $C>0$ such that $a_n/b_n\le C$ (resp. $1/C\le a_n/b_n\le C$) for all large $n$, and say $a_n={\scriptstyle{\mathcal{O}}}(b_n)$ if $a_n/b_n\rightarrow0$ as $n\rightarrow\infty$. We set $(X_n)$ and $(Y_n)$ to be two sequences of random variables. Write $X_n={\mathcal{O}}_{p}(Y_n)$ if for $\forall \epsilon>0$, there exists $C>0$ such that $\mathrm{P}(|X_n/Y_n|\leq C)>1-\epsilon$ for all large $n$, and say $X_n={\scriptstyle{\mathcal{O}}}_{p}(Y_n)$ if $X_n/Y_n\rightarrow 0$ in probability as $n\rightarrow\infty$. We use $\text{plim}$ to denote the probability limit. Further, for a real matrix $\mathbf A$, let $\|\mathbf A\|_F={\rm tr}^{1/2}(\mathbf A^{\top}\mathbf A)$ and $\|\mathbf A\|_2=\lambda_{\max}^{1/2}(\mathbf A^{\top}\mathbf A)$ denote its Frobenius and spectral norms, respectively. Let $\lambda_{\min}(\cdot)$ and $\lambda_{\max}(\cdot)$ denote the minimum and maximum eigenvalues of a square matrix. For a vector $\textbf{v}$, let $\|\textbf{v}\|$ denote its Euclidean norm.
For any scalar or vector function $\varphi(z)$, we use the notation $\Pi_{\mathcal{G}}[\varphi(z)]$ to denote the projection of $\varphi(z)$ onto the class of functions $\mathcal{G}$. That is, $\Pi_{\mathcal{G}}[\varphi(z)]$ is an element that belongs to $\mathcal{G}$ and is the closest function to $\varphi(z)$ among all the functions in $\mathcal{G}$. More specifically, we have
where the infimum is in the sense that
for all $h\in\mathcal{G}$, where for square matrices $\mathbf A$ and $\mathbf B$, $\mathbf A\leq \mathbf B$ means that $\mathbf A-\mathbf B$ is negative semidefinite.
Denote $\theta(z)=\mathop{\mbox{\sf E}}[\mathbf{X}_t|\mathbf{Z}=z]$ and $m(z)$ is the projection of $\theta(z)$ onto $\mathcal{G}$, i.e., $m(z)=\mathop{\mbox{\sf E}}_{\mathcal{G}}[\theta(z)]$. For $t=1,\ldots,T$ we define $\boldsymbol{\xi}_t = \mathbf{X}_t -\boldsymbol{m}(\mathbf{Z})$, $\boldsymbol{\eta}(\mathbf{Z}) = \theta(\mathbf{Z})-\boldsymbol{m}(\mathbf{Z})$, and $\boldsymbol{\varepsilon}_t= \mathbf{X}_t -\theta(\mathbf{Z})$, where $\boldsymbol{\xi}_t$, $\boldsymbol{\eta}(\mathbf{Z})$, and $\boldsymbol{\varepsilon}_t$ are $N\times Q$ matrices. Also, the following conditions about the data generating process, basis functions, factor loadings, and sieve approximation are required to obtain the large sample properties of the proposed estimator, $\boldsymbol{\boldsymbol{\widehat{\boldsymbol{\beta}}}}$.
Assumption (ref) allows for correlation between $\mathbf{X}_{it}$ and $\mathbf{Z}_{i}$ through $\theta(\mathbf{Z})$, $m(\mathbf{Z})$ and $\eta(\mathbf{Z})$. This assumption is standard in semiparametric estimation techniques (see for example Assumption 2.1(ii) in AHMADLEELAHANONLI:2005). Also, Assumption (ref) is commonly used in the estimation of factor models and enables to identify separately the factors $\mathbf{F}$ (see condition PC1 in Bai-Ng2013). Assumption (ref) (iii) guarantees that $(NT)^{-1}\sum_{t=1}^T \mathbf{X}_t^{\top} \left[\mathbf{I}_N- \mathbf P_{\Phi}(\mathbf{Z})\right]\mathbf{X}_t$ is asymptotically nonsingular.
As already pointed out in Fan-Liao-Wang2016, $N^{-1}\Phi(\mathbf{Z})^{\top}\Phi(\mathbf{Z})=N^{-1}\sum^N_{i=1}\phi\left(\mathbf{Z}_i\right)^{\top}\phi\left(\mathbf{Z}_i\right)$ and $\phi\left(\mathbf{Z}_i\right)$ is of order $J_ND$ much smaller than $N$. Thus, condition (i) can follow from a strong law of large numbers. This condition can be satisfied through proper normalizations of commonly used basis functions.
The following set of conditions is concerned with the accuracy of the sieve approximation.
As it is remarked in Fan-Liao-Wang2016, Assumption (ref) (ii) is satisfied by the use of common basis functions such as polynomial basis or B-splines. In particular, Lorentz1986 and Chen2007 show that (i) implies (ii) in this particular case.
The next assumption refers to the error factor loadings $\mathbf{\gamma}_i$, for $i=1,\ldots,N$.
Note that in Assumption 3.5 (ii) we assume cross-sectional dependence of the error factor loadings. To show the consistency of the proposed estimator for simplicity we can assume the independence of the factor loadings $\gamma_{ik}$ from the random part of the covariates, $\textbf{Z}_i$, but we do not need to impose a restrictive i.i.d. assumption.
Through the paper, some regularity conditions about weak dependence and stationarity are assumed on the factors and the idiosyncratic terms. In particular, we impose strong mixing conditions. Let $\mathcal{F}_{-\infty}^0$ and $\mathcal{F}_{T}^{\infty}$ denote the $\sigma$-algebras generated by $\{(\boldsymbol\xi_t,\mathbf{f}_t,\mathbf u_t):t\leq 0\}$ and $\{(\boldsymbol\xi_t,\mathbf{f}_t,\mathbf u_t):t\geq T\}$, respectively. Define the mixing coefficient \[ \alpha(T)=\sup_{A\in\mathcal{F}_{-\infty}^0,B\in\mathcal{F}_T^{\infty}}|\mathrm{P}(A)\mathrm{P}(B)-\mathrm{P}(AB)|. \]
Assumption (ref) is standard in factor analysis Bai2003,Stock-Watson2002,Fan-Liao-Wang2016. Part (i) is standard in partially linear models AHMADLEELAHANONLI:2005, Hardle2000. The independence assumption between $\mathbf u_t$ and $\{\mathbf{Z}_i, \boldsymbol\xi_t\}$ can be relaxed by allowing for conditional independence. Part (iii) is a strong mixing condition for the weak temporal dependence of $\{\boldsymbol\xi_t,\mathbf u_t,\mathbf{f}_t\}$, whereas (iv) imposes weak cross-sectional dependence in $\{u_{it}\}_{i\leq N, t\leq T}$. This condition is usually satisfied when the covariance matrix of the error term $u_{it}$ is sufficiently sparse under the strong mixing condition and it is commonly imposed for high-dimensional factor analysis.
A very intuitive idea of the asymptotic behavior of our estimator can be obtained by plugging ((ref)) in ((ref)) that yields
where $\mathbf M_{\Phi}(\mathbf{Z})=\mathbf{I}_N-\mathbf{P}_{\Phi}(\mathbf{Z})$ and $\mathbf{\Lambda} = \mathbf{G(Z)}+\mathbf{\Gamma}$. As the reader can see in the above expression, there is a direct dependence of $\boldsymbol{\widehat{\boldsymbol{\beta}}}$ on the unobserved factor loadings through $(NT)^{-1}\sum_t\mathbf{X}_t^{\top}\mathbf M_{\Phi}(\mathbf{Z})\boldsymbol \Lambda \mathbf{f}_t$. Nevertheless, using ((ref)) and given that it can be proved that $(NT)^{-1}\sum_t\mathbf{X}_t^{\top}\mathbf M_{\Phi}(\mathbf{Z})\mathbf{G}(\mathbf{Z})\mathbf{f}_t={\scriptstyle\mathcal{O}}_p(1/\sqrt{NT})$ (see the proof of Theorem 3.1 in the Supplementary Material document), we have that
In this situation, we can conclude that the limiting distribution of $\boldsymbol{\boldsymbol{\widehat{\boldsymbol{\beta}}}}-\boldsymbol{\boldsymbol{\beta}}$ only depends on idiosyncratic terms (related to both the error term and the approximation error of the basis functions to the factor loadings). Under Assumptions (ref)--(ref) it is also possible to show (see the Appendix A of the Supplementary Material document for a related proof) that, as both $N$ and $T$ tend to infinity,
The interesting feature of this asymptotic bound is that the rate of convergence of $\widehat{\boldsymbol\beta}$ can be slower than $\sqrt{NT}$ depending on the behavior of $\nu_N$ (i.e. the variance of the error factor loadings). In order to clarify this, we further take a look at the asymptotic distribution of our estimator. The previous result on the consistency and the convergence rate is based on weak dependence in the error term and idiosyncratic factor loadings. To show asymptotic normality, we have to impose the stronger condition of cross-sectional independence while still allowing for weak dependence in the time dimension.
Then, the asymptotic distribution of the projection-based interactive fixed effects estimator $\boldsymbol{\boldsymbol{\widehat{\boldsymbol{\beta}}}}$ is provided in the following theorem.
The proof of Theorem (ref) is provided in Appendix B.1 in the Supplementary Material document. The key component of the proof is following the Frisch-Waugh Theorem to partial-out the effect of the latent factors and corresponding loadings. We want to highlight that the relative rate requirements of $N$ and $T$ crucially depend on the smoothness parameter, $\kappa$. In particular, if $\kappa=4$ we have the requirement that $T/N$ tends to zero regardless of the choice for the sieve dimension $\varrho$. The constraints in the rates of growth imposed on $N$ and $T$ are similar to other assumptions used in similar literature such as in AHMADLEELAHANONLI:2005. Note that Assumption (ref) is introduced for the sake of simplicity. It is indeed used for the application of the corresponding central limit theorems (CLT) in the proof but it could be relaxed at the cost of a much cumbersome proof.
The results of Theorem (ref) have important implications for conducting inference on the estimated regression parameters. In particular, since the variance of the idiosyncratic part of the factor loadings decides which of the two terms will be the leading one, it will ultimately determine the convergence rate of our estimator. As a consequence, the asymptotic distribution of $\boldsymbol{\widehat{\boldsymbol{\beta}}}-\boldsymbol{\beta}$ has a discontinuity when the variance of the factor loadings is close to the boundary. In financial econometrics studies, this issue is often circumvented by assuming weak heterogeneity as a default setting Connor-Linton2007,connor2012efficient. However, recent studies such as Fan-Liao-Wang2016 find empirical evidence for the case of strong heterogeneity. Further, usual plug-in approaches based on estimated asymptotic covariance matrices will lead to misleading conclusions if $\nu_N=\mathcal{O}(T^{-1})$, as they provide confidence intervals that are too wide because the asymptotic covariance matrix is over-estimated, and the coverage probabilities will be too conservative liao2018uniform,FernndezVal2022DynamicHD. Similarly, simply ignoring the cross-sectional term will lead to under-coverage in the strong heterogeneity case.
Uniformly valid inference for panel data models is an important topic beyond the specifics of our model setup. For instance, lu2023uniform observe a similar issue in a panel model with two-dimensional heterogeneity in the regression parameters. In their case, the issue is caused by the level of temporal and cross-sectional heterogeneity in the slope coefficients.
Fortunately, the uniformity issue can be solved by using the cross-sectional bootstrap proposed by kapetanios2008bootstrap. Besides achieving uniformly valid inference, the approach is both intuitive and easy to implement. The basic idea is to sample with replacement cross-sectional units while keeping the entire time series of the sampled individual units unchanged. By doing this, the resampling scheme directly mimics the cross-sectional variations in $\boldsymbol \Gamma$, regardless of the underlying level of heterogeneity. The consequence is a uniformly valid inference. This is in direct contrast to andrews2000inconsistency, who found that the usual bootstrap will lead to inconsistency when a parameter is on the boundary of the support. The reason why this problem does not occur in our case is that we do not explicitly model the variance of the idiosyncratic factor loadings as a parameter, i.e., it does not appear in the loss function of our least squares problem.
A crucial assumption for the bootstrap validity is that the data is cross-sectionally independent. In fact, menzel2021bootstrap showed that uniform non-conservative inference is impossible under general dependence in more than one dimension. Recently, de2024cross studied the theoretical properties of the cross-sectional bootstrap for the CCE approach of Pesaran2006 and proposed a bias-correction procedure in the asymptotic regime $N/T\to\rho<\infty$. The uniform validity of the bootstrap procedure in settings similar to ours was recently shown in liao2018uniform and FernndezVal2022DynamicHD. The specific aspect of our procedure is that we only resample cross-sectional units after projecting the data, i.e., removing the effect of $\boldsymbol{Z}_i$ on the factor loadings.
As a positive side effect, the cross-sectional bootstrap is able to keep the dependence in the time dimension. Therefore, the inference is also robust towards serial dependence in the idiosyncratic error term and in the latent factors. In the following, we summarize the steps of the cross-sectional bootstrap procedure.
For the bootstrap validity, we need to assume the existence of a consistent estimator of the variance of $v^{\top}\widehat{\boldsymbol{\beta}}$.
The following theorem provides the bootstrap validity, uniformly over settings with varying degrees of variability in the idiosyncratic factor loadings.
The proof of Theorem (ref) can be found in Appendix B.2 in the Supplementary Material document. An essential part of the proof is to show that the asymptotic expansion of the bootstrap version of the estimator is identical to that of the original estimator.
In this section, we evaluate the finite-sample performance of our estimator in a simulation study. We are interested both in the estimation accuracy of the parameter vector, $\boldsymbol{\beta}$, and the empirical coverage probabilities of the cross-sectional bootstrap procedure. Throughout the study, we fix the number of factors, $K=3$, the dimension of the time-invariant variable is set to $D=2$, and the dimension of covariates is set to $Q=2$. The true regression coefficients are $\boldsymbol{\beta}=(2,-1)^{\top}$. The time-invariant variables are generated by i.i.d. $Z_{id}\sim U[-1,1]$. The covariates are generated by setting $X_{itq}=\mathbf a_{iq}^\top \mathbf{f}_t+2(\sqrt{ g_1(\mathbf{Z}_{i})},\ldots,\sqrt{g_K(\mathbf{Z}_{i})})^{\top}\boldsymbol{b}_q+\pi_{itq}$, where $\pi_{itq}\sim N(0,1)$ i.i.d., $a_{iqk}\sim U[-0.5,0.5]$ and $b_{qk}\sim U[-1,1]$. We generate the latent factors, $(f_{k1},\ldots,f_{kT})$, as MA$(\infty)$ processes with algebraic decay and under independence across factors for all $k$. The factor loadings are set to $\lambda_{ik}=g_{k}(\mathbf{Z}_{i})+\gamma_{ik}$, where $g_1(z)=\sin(2z_1)^3+\cos(z_2^2)$, $g_2(z)=-\tan(z_1^2)+2\cos(z_2+1)$ and $g_3(z)=z_2^3-\sin(3z_1)$.
Finally, for the idiosyncratic error term, we consider the case of i.i.d. standard normal $u_{it}$ as well as the case of weak temporal dependence, in which $(u_{i1},\ldots,u_{iT})$ are generated from a MA$(\infty)$ process with algebraic decay parameter $5$. For the idiosyncratic part of the factor loadings we consider three settings. First, in the strong heterogeneity case $\boldsymbol\gamma_{i}\sim N(0,0.5)$ (i.e., $\nu_N=\mathcal{O}(1)$). Second, we consider the special case $\nu_N=0$. Third, we consider the weak factor case in which $\boldsymbol\lambda_i^{\top}\mathbf{f}_t=\mathcal{O}(T^{-1/2})$ (i.e., $\nu_N=\mathcal{O}(T^{-1})$). In this numerical study, we rely on B-spline basis functions and we select $J_N=\lceil N^{1/3}1.5\rceil$. For each setting $500$ Monte Carlo runs are conducted.
We compare the performance of our projection-based interactive fixed effects (P-IFE) estimator for $\boldsymbol{\beta}$ with the principal component-based interactive fixed effects (PC-IFE) estimator of Bai2009 and a bias-corrected version of the same estimator (bc-PC-IFE). For these comparisons, we rely on the R package phtt bada2014phtt. The number of factors is selected according to the PC1 criterion in Bai-Ng2002. As performance measures, we consider the root mean square error ($RMSE$). The simulation results under Gaussian disturbances for different values of $\nu_N$, $N$, and $T$ are reported in Table (ref). The $RMSE$ of our P-IFE can be effectively reduced with increasing sample size. For the strong heterogeneity case, we observe an advantage of the PC-IFE for small and medium samples. For $N=500$ this advantage is reversed and the P-IFE has a higher accuracy. In the other two settings for $\nu_N$, we can see that the P-IFE outperforms its competitors in almost all cases. The outperformance is best visible for settings with large sample sizes in the case of $\nu_N=0$. In particular, for $N=500$ the $RMSE$ of our P-IFE is less than a third of that of the PC-IFE.
The results for serially dependent error terms are displayed in Table (ref). Again, we can observe that the PC-IFE outperforms the P-IFE in $\nu_N=\mathcal{O}(1)$ case for small and medium sample sizes. Also similar to the i.i.d. case, the P-IFE has the lowest $RMSE$ in all settings for $\nu_N=0$ and $\nu_N=\mathcal{O}(T^{-1})$. Interestingly, the bias-corrected estimator performs worse than the estimator without bias correction. As a robustness check, we also consider $t$-distributed error terms in the Supplementary Material document. See Table (ref) for the results. We also consider a more complicated additional data-generating process with a larger number of factors, $K=10$. The results are displayed in Table (ref) for i.i.d. errors and Table (ref) for serially correlated errors. The most notable difference to the results of the first DGP is that the P-IFE outperforms the PC-IFE even in the strong heterogeneity case for settings with a small time dimension, $T=10$. Again, the P-IFE dominates in all settings for $\nu_N=0$ and $\nu_N=\mathcal{O}(T^{-1})$.
In the following, we look at the performance of the cross-sectional bootstrap procedure and show its validity in finite samples. As a comparison, we look at the empirical coverage of the PC-IFE estimator. By Corollary 1 in Bai2009, under the assumption of i.i.d. error terms, $\sqrt{NT}(\boldsymbol{\widehat{\boldsymbol{\beta}}}_{\text{PC-IFE}}-\boldsymbol{\beta})\stackrel{\mathcal{L}}{\rightarrow}N(0,\sigma^2\mathbf D^{-1})$, where $\mathbf D=\text{plim}(NT)^{-1}\sum_{i=1}^N \mathbf Z_i^\top \mathbf Z_i$, $\mathbf Z_i=\mathbf M_F\mathbf X_i-1/N \sum_{k=1}^N\mathbf M_F\mathbf X_ka_{ik}$, $\mathbf M_F=\mathbf{I}_N-\mathbf F\mathbf F^\top/T$ and $a_{ik}= \boldsymbol\lambda_i^\top(\boldsymbol \Lambda^\top\boldsymbol \Lambda)^{-1} \boldsymbol\lambda_k$. We construct confidence intervals based on the asymptotic distribution with an estimated covariance matrix based on estimated factors and factor loadings.
Table (ref) shows that the empirical coverage of our cross-sectional bootstrap procedure approaches the nominal coverage level as $N$ and $T$ increase. As is often the case, we can observe slight under-coverage in small samples. However, the issue becomes virtually absent in settings with the largest sample size. We want to highlight that these findings hold for all settings for the variance of the idiosyncratic factor loadings, $\nu_N$. We have thus provided evidence for the uniform validity of the bootstrap procedure in finite samples. In the Supplementary Material document, we show that the cross-sectional bootstrap procedure is also robust towards $t$-distributed errors and serially correlated error terms. See Tables (ref) and (ref). Our uniform bootstrap procedure naturally adapts to the data, particularly to potential serial dependence and a varying degree of heterogeneity in the factor loadings.
Looking at the coverage of the asymptotic distribution of the PC-IFE estimator, we can observe under-coverage in all settings. Moreover, the coverage does not improve with increasing sample size. On the contrary, the coverage is worst for the setting with $N=500$ and $T=100$.
The aim of this section is to show the performance of our estimator in empirical analysis. More precisely, we will apply our estimator in the analysis of the determinants of economic growth. We refer to durlauf2005growth for a comprehensive review of the growth literature. While many studies focus on a cross-sectional analysis (see for instance barro1991economic), there are also numerous studies employing a panel data approach with country-specific fixed effects acemoglu2019democracy,islam1995growth. However, lu2016shrinkage argue that economic growth rates might not be solely determined by observable regressors, but could also be influenced by latent factors or shocks. Our projection-based interactive fixed effect estimator is well suited as it is flexible enough to model such latent factors.
The yearly data on GDP growth rates and the country-specific characteristics are obtained from the Penn World Table (PWT) and the World Bank World Development Indicators (WDI). Our sample contains 129 countries in a period from 1991--2019, $N=129$ and $T=29$. Countries with incomplete data availability or which did not exist yet in 1991 are excluded from our analysis. Our dependent variable is the real GDP growth rate per capita. The set of regressors is identical to the regressors in lu2016shrinkage. Summary statistics of all dependent and independent variables can be found in Table (ref). Figure (ref) shows the time series of the mean growth rates, averaged over all countries in our sample. We also visualize the time series of the cross-sectional $5\%$ and $95\%$-quantiles of the growth rates in the same figure. For the time-invariant characteristics used for modeling the systematic part of the factor loadings we take the longitude and latitude of the respective country\footnote{Data obtained from developers.google.com/public-data/docs/canonical/countries_csv.}.
We first fit our projection-based interactive fixed effects model using the complete sample of $N=129$ countries. To be consistent with the simulation section, we use B-spline basis functions with $J_N=\lceil N^{1/3}1.5\rceil$. The estimation results can be found in Table (ref). We report the estimated coefficients and the $95\%$ confidence interval based on the cross-sectional bootstrap with $1000$ bootstrap iterations. As a comparison, we also report the estimated coefficients and confidence intervals following the PC-IFE approach of Bai2009. We obtain negative significant coefficients for consumption share and government consumption share and a positive significant coefficient for the age dependency ratio at the $5\%$ confidence level. These results are similar to the estimation based on the PC-IFE. However, investment share and fertility rate also become significant for the PC-IFE. The remaining variables are insignificant for both estimation procedures.
We now restrict our analysis to the subset of countries that are members of the OECD (Organisation for Economic Cooperation and Development). See Table (ref) for the estimation results. Similar to the previous results, both approaches find a negative significant effect on the consumption share. However, our P-IFE identifies a positive effect on investment share and a negative effect on the investment price level. Both variables are insignificant for the PC-IFE approach. Moreover, PC-IFE additionally finds significant effects on the fertility rate and life expectancy.
In this paper, a new estimator for the regression parameters in a panel data model with interactive fixed effects has been proposed. The main novelty of this approach is that factor loadings are approximated through nonparametric additive functions, and it is then possible to partial out the interactive effects. Therefore, the new estimator adopts the well-known partial least squares form, and there is no need to use iterative estimation techniques to compute it. It turns out that the limiting distribution of the estimator has a discontinuity when the variance of the idiosyncratic parameter of the factor loading approximation is near the boundaries. The discontinuity makes the usual “plug-in” inference based on the estimated asymptotic covariance matrix problematic since it can lead to either over- or under-coveraging probabilities. We show that non-conservative uniformly valid inference can be achieved by cross-sectional bootstrap. A Monte Carlo study indicates good performance in terms of mean squared error and bootstrap coverage. We apply our methodology to analyze the determinants of growth rates in OECD countries.
Juan M. Rodriguez-Poo and Alexandra Soberon acknowledge financial support from the I+D+i project Ref. PID2019-105986GB-C22 financed by MCIN/AEI/10.13039/ 501100011033. In addition, this work was also funded by the I+D+i project Ref. TED2021-131763A-I00 financed by MCIN/AEI/10.13039/501100011033 and by the European Union NextGenerationUE/PRTR. Georg Keilbar acknowledges gratefully the support from the Deutsche Forschungsgemeinschaft via the IRTG 1792 "High Dimensional Nonstationary Time Series". Weining Wang's research is partially supported by the ESRC (Grant Reference: ES/T01573X/1).
\spacingset{1.3}
{
\centerline{SUPPLEMENTARY MATERIAL} }
\spacingset{1.8}
This supplementary document is organized as follows. In Section (ref) we provide the proof of lemmas. The proofs for Theorems 3.1 and 3.2 are in Section (ref). We provide more details on the estimation of latent factors and factor loadings in Section (ref). Finally, we provide more details on the simulation study in Section (ref).