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.
99,790 characters · 10 sections · 171 citation commands
Inference in Unbalanced Panel Data Models with Interactive Fixed Effects
\onehalfspacing
Economists are often concerned that unobserved heterogeneity is correlated with some regressors, leading to inconsistent estimates of the parameters of interest. When panel data are available, fixed effects models are frequently used to address this issue. A critical assumption of these models is that unobserved heterogeneity enters additively. If this fails, for example, because an unobserved financial crisis shock affects each country's output differently, fixed effects models are no longer appropriate (see b2009 for additional motivating examples). This concern motivates interactive fixed effects (IFE) estimators, which model unobserved heterogeneity as a low-rank factor structure $\boldsymbol{\lambda}_{i}^{\prime} \mathbf{f}_{t}^{\phantom{\prime}}$, where $\boldsymbol{\lambda}_{i}$ and $\mathbf{f}_{t}$ are unit- and time-specific effects, respectively (see, among others, hnr1988, p2006, and b2009).\footnote{bm2015 suggest a related but different approach. Instead of imposing rank restrictions on the time-varying unobserved heterogeneity, they use a clustering approach to assign each cross-sectional unit to a specific group, where the corresponding group-specific heterogeneity is allowed to vary over time.} Throughout this article, we refer to $\boldsymbol{\lambda}_{i}$ as factor loadings and $\mathbf{f}_{t}$ as common factors.
Inspired by ah1982, hnr1988 propose a quasi-differencing approach for panels with large $N$ but small $T$. They first remove factor loadings from the estimation equation and then estimate the remaining common factors and parameters using lagged regressors as instruments. While this estimator is consistent under asymptotic sequences in which $T$ is fixed, it is well known that for large $T$, the number of instruments and parameters causes bias (see ns2004). More recent work has considered estimators that require both $N$ and $T$ to be large. p2006 proposes a common correlated effects (CCE) estimator in the spirit of \textcites{m1978}{c1982}{c1984}, which uses cross-sectional averages of the dependent variable and the regressors to proxy for the unobserved common factors. p2006's estimator is at least $\sqrt{N}$-consistent without requiring knowledge of the true rank of the factor structure or strong factor assumptions as in \textcites{b2009}{mw2015}{mw2017}. However, it requires additional parametric assumptions on the joint distribution of the dependent variable and the regressors in order to use cross-sectional averages as valid proxy variables. b2009 proposes a different estimator that treats the common factors and factor loadings as additional parameters.\footnote{For a detailed discussion of the different interactive fixed effects estimators, we refer the reader to \textcites{b2009}{mw2015}{mw2017}.} This estimator is closely related to b2003's principal components estimator for pure factor models and has the advantage of not requiring distributional assumptions about the unobserved heterogeneity. Under the assumption that the true number of factors is known, b2009 establishes $\sqrt{NT}$-consistency irrespective of cross-sectional and/or time-serial dependence in the idiosyncratic error term. Such dependence does, however, induce an asymptotic bias in the limiting distribution, which can be corrected (see {\citeresetb2009}). mw2017 derive an additional correction for the feedback bias (essentially a n1981-type bias) that arises from the inclusion of predetermined regressors such as lagged dependent variables. Because the true number of factors is generally unknown, mw2015 show that as long as the number of factors used to estimate $\boldsymbol{\beta}$ exceeds the true number, the estimator remains at least $\sqrt{\min(N, T)}$-consistent, at the potential cost of some efficiency loss from including irrelevant factors. Given a consistent estimator of $\boldsymbol{\beta}$, the number of factors can then be estimated using estimators for pure factor models (see, among others, be1992, bn2002, hl2007, abc2010, o2010, ah2013, and do2019). A recent comparison of popular estimators for pure factor models is given in cj2019.
In applied work, observations are often missing. A frequent cause is attrition: individuals may drop out of a panel because they move or leave the participating household, and in some cases they are replaced by new survey participants. In macroeconomic panels, countries are sometimes divided into several independent states. Non-response can also lead to the replacement of survey participants. These cases give rise to very different missing data patterns that, in the absence of sample selection, generally do not affect the properties of estimators (see fw2018). In the presence of missing data, the principal component estimator of b2009 requires an additional imputation step based on the EM algorithm of \textcites{sw1998}{sw2002} (see the appendix of b2009 and bly2015). bly2015 demonstrate consistency of the EM-type principal component estimator through simulation studies, but provide no guidance on inference. The asymptotic properties of the EM algorithm for factor models were recently studied by jms2021.
We make the following contributions. First, we extend mw2017 using the insights of fw2018 to derive the asymptotic distribution of the IFE estimator in unbalanced panels under the assumption that attrition is conditionally random. Second, we propose a novel method of alternating projections algorithm to compute the residuals required for inference. The algorithm relies on straightforward scalar expressions and is particularly suited to settings with missing data, though it can also be applied to balanced panels. We also propose an alternative estimation procedure to those of \textcites{b2009}{bly2015}. Specifically, we combine the profile-objective-function reformulation of \textcites{mw2015}{mw2017} with matrix completion methods such as the EM algorithm. This procedure eliminates the need to optimize explicitly over the high-dimensional nuisance parameters $(\boldsymbol{\lambda}_{1}^{\prime}, \ldots, \boldsymbol{\lambda}_{N}^{\prime})^{\prime}$ and $(\mathbf{f}_{1}^{\prime}, \ldots, \mathbf{f}_{T}^{\prime})^{\prime}$, and typically converges in few iterations. We also present a regularization-based matrix completion approach as an alternative to the EM algorithm, which is particularly advantageous for larger-scale panels, as noted by fll2021. An R package implementing all proposed methods is available at \url{https://github.com/dczarnowske/InteractiveEffects}. Third, we analyze the finite-sample properties of the IFE estimator for a dynamic model through simulation experiments that explore different shares of missing data, confirming that our asymptotic results provide a reliable approximation to finite-sample behavior. Additional Monte Carlo results for static models, covering various error term configurations and missing data patterns, are reported in Section (ref) of the Online Supplement. Fourth, given that our results assume the true number of factors is known, we also examine the performance of various estimators for this quantity. For sufficiently long panels, all estimators perform similarly regardless of the share of missing observations. In configurations with high persistence, only a few estimators achieve reliable predictions; without high persistence, all estimators predict the correct number of factors almost perfectly. Fifth, we reassess the baseline analysis of anrr2019 using the IFE estimator. Our findings qualitatively confirm their main results. However, in their preferred specification, the estimated short-run and long-run effects are roughly halved relative to those they report. Sixth, our findings and algorithms extend to several related estimators, including the minimum distance estimator of \textcites{mw2017}{msw2018} for endogenous regressors, the nuclear norm regularized estimators of mw2026, and the estimator for nonlinear factor models of cfw2021.
Related work to ours is sww2026.\footnote{Some of the ideas in sww2026 build on our earlier work. The bias correction formulas and algorithms presented in our paper first appeared in an arXiv preprint circulated in April 2020 (see cs2020).} Their approach extends the asymptotic expansion of fw2016, while ours builds on the expansion of mw2017. Both expansions eliminate the effects of high-dimensional nuisance parameters through projections, resulting in related inference procedures.
The paper is organized as follows. Section (ref) introduces the model and presents estimation and inference procedures. Section (ref) briefly reviews estimators for the number of factors. Section (ref) presents simulation results. Section (ref) reassesses anrr2019 using the IFE estimator. Section (ref) presents algorithms for related estimators. Section (ref) concludes.
Throughout this article, we follow standard notation: scalars are in roman type, vectors and matrices in boldface, and all vectors are column vectors. Let $\mathbf{A}$ be an $M \times N$ matrix. We write $[\mathbf{A}]_{ij}$ for the $(i,j)$-th element of $\mathbf{A}$, where $i$ is a row index and $j$ is a column index. $\operatorname{\mathbb{I}}_{M}$ denotes the $M \times M$ identity matrix.
We consider the following unobserved effects model:
where $i$ and $t$ denote individual and time indices, respectively. Here, $\mathbf{x}_{it} \coloneqq (x_{it, 1}, \ldots, x_{it, K})^{\prime}$ is a vector of $K$ regressors, $\boldsymbol{\beta}$ is the corresponding parameter vector, and $e_{it}$ is the idiosyncratic error term. To accommodate unbalanced panels, let $\mathcal{D} \subseteq \{(i, t) \colon i \in \{1, \ldots, N\} \times t \in \{1, \ldots, T\}\}$ denote the set of observed index pairs, where $N$ and $T$ denote the number of individuals and time periods, respectively, and $n \coloneqq \lvert\mathcal{D}\rvert$ denotes the total sample size. The unobserved effects in (ref) are modeled as a factor structure, with $\boldsymbol{\lambda}_{i} \coloneqq (\lambda_{i1}, \ldots, \lambda_{iR})^{\prime}$ a vector of factor loadings and $\mathbf{f}_{t} \coloneqq (f_{t1}, \ldots, f_{tR})^{\prime}$ a vector of common factors. The factor structure is assumed to be of low rank, $R \ll \min(N, T)$.
Given the number of factors $R$, the estimator of the common parameters is defined as:
where
is the profile objective function. Here, $\boldsymbol{\Lambda} \coloneqq (\boldsymbol{\lambda}_{1}, \ldots, \boldsymbol{\lambda}_{N})^{\prime}$ is an $N \times R$ matrix of factor loadings and $\mathbf{F} \coloneqq (\mathbf{f}_{1}, \ldots, \mathbf{f}_{T})^{\prime}$ is a $T \times R$ matrix of common factors.
Let $\mathcal{C}$ be a conditioning set containing the sigma-algebra generated by the true factor loadings and common factors, and let $\mathcal{Z}_{i}^{t} \coloneqq \sigma(\{(\mathbf{x}_{is}, e_{i(s - 1)}) \colon s \leq t\})$ hold for all $i, t, N, T$. For balanced panels, mw2017 derived the asymptotic distribution of the interactive fixed effects estimator (ref) under an asymptotic framework in which $N, T \rightarrow \infty$ at the constant rate $N / T \rightarrow \kappa^{2}$ with $0 < \kappa < \infty$. Their assumptions further require that the true number of factors is known, that $\{(\mathbf{x}_{it}, e_{it}) \colon t = 1, \ldots, T\}$ is independent across $i$ (conditional on $\mathcal{C}$), that $\operatorname{\mathbb{E}}[e_{it} \mid \mathcal{C} \vee \mathcal{Z}_{i}^{t}] = 0$ holds for all $i, t, N, T$, and that the regressors are not fully absorbed by the factor structure (a non-collinearity condition).
As argued by fw2018 in Section 4.1, missing observations do not pose major theoretical challenges if the attrition process is deterministic or conditionally random. Let $\mathcal{I}_{t} \coloneqq \{i \colon (i, t) \in \mathcal{D}\}$, $\mathcal{T}_{i} \coloneqq \{t \colon (i, t) \in \mathcal{D}\}$, and $\delta_{it}$ be an attrition indicator for all $i, t, N, T$. To derive the asymptotic distribution of (ref) for unbalanced panels, we augment the assumptions of mw2017 with one of the following assumptions.
Let $\bar{p}_{itt^{\prime}}^{f} \coloneqq \mathbf{f}_{t}^{\prime} (\sum_{t^{\prime} = 1}^{T} \operatorname{\mathbb{E}}[\delta_{it^{\prime}} \mid \mathcal{C}] \, \mathbf{f}_{t^{\prime}}^{\phantom{\prime}} \mathbf{f}_{t^{\prime}}^{\prime})^{- 1} \mathbf{f}_{t^{\prime}}$, $\bar{\xi}_{it}^{\dagger} \coloneqq \boldsymbol{\lambda}_{i}^{\prime} (\boldsymbol{\Lambda}^{\prime} \overline{\boldsymbol{\Phi}} \boldsymbol{\Lambda})^{- 1} (\boldsymbol{\Lambda}^{\prime} \boldsymbol{\Lambda}) \mathbf{f}_{t}^{\phantom{\prime}}$, $\overline{\boldsymbol{\Phi}} \coloneqq \operatorname{\mathbb{E}}[(\boldsymbol{\Delta} \odot \boldsymbol{\Lambda} \mathbf{F}^{\prime}) (\boldsymbol{\Delta} \odot \boldsymbol{\Lambda} \mathbf{F}^{\prime})^{\prime} \mid \mathcal{C}]$, and $\boldsymbol{\Delta}$ be an $N \times T$ matrix with $[\boldsymbol{\Delta}]_{it} = \delta_{it}$. Let $\odot$ denote the Hadamard (element-wise) product, $\mathbf{A}^{\cdot} \coloneqq (\mathbf{a}_{1}^{\cdot}, \ldots, \mathbf{a}_{T}^{\cdot})^{\prime}$, and $\mathbf{C}^{\cdot} \coloneqq (\mathbf{c}_{1}^{\cdot}, \ldots, \mathbf{c}_{N}^{\cdot})^{\prime}$, where the dot in the exponent is a placeholder. Under the assumptions of mw2017 augmented by Assumption (ref), the estimator in (ref) has the following asymptotic distribution when data are conditionally missing at random:
where
with
denoting residuals from population projections. When the attrition process is deterministic, we augment the assumptions of mw2017 by Assumption (ref), and the asymptotic distribution follows immediately by replacing $\operatorname{\mathbb{E}}[\delta_{it} \mid \mathcal{C}]$ with $\delta_{it}$ in (ref). The derivation is provided in Appendix (ref).
The bias term $\mathbf{B}_{1}$ represents feedback bias (a generalization of the n1981-bias) arising from potential feedback from past outcomes to future realizations of the regressors. Specifically, $\mathbf{x}_{it}$ may depend on $(e_{i(t - 1)}, e_{i(t - 2)}, \ldots)$, $\boldsymbol{\lambda}_{i}$, and $\mathbf{f}_{t}$ in an arbitrary nonlinear manner. Our framework thus naturally accommodates dynamic specifications such as $\mathbf{x}_{it} = y_{i(t - 1)}$. Feedback bias is ruled out by assumption in b2009 and was first introduced by mw2017.
The remaining bias terms, $\mathbf{B}_{2}$ and $\mathbf{B}_{3}$, are also present in b2009 and mw2017 for balanced panels. They arise when the idiosyncratic error term is heteroskedastic across individuals or over time, respectively. The reason is that $\sum_{i = 1}^{N} \sum_{t = 1}^{T} \operatorname{\mathbb{E}}[\ddot{\mathbf{x}}_{it}^{\lambda} \mid \mathcal{C}] \, \bar{\xi}_{it}^{\dagger} = 0$ and $\sum_{i = 1}^{N} \sum_{t = 1}^{T} \operatorname{\mathbb{E}}[\ddot{\mathbf{x}}_{it}^{f} \mid \mathcal{C}] \, \bar{\xi}_{it}^{\dagger} = 0$ follow from the definition of the population residuals. In unbalanced panels, the attrition process can induce a form of heteroskedasticity. Consequently, even when the idiosyncratic error term is homoskedastic, $\mathbf{B}_{2}$ and $\mathbf{B}_{3}$ are generally non-zero, since missing probabilities may also be heterogeneous. This finding is consistent with sww2026.
The covariance matrix $\mathbf{W}^{- 1} \boldsymbol{\Omega} \mathbf{W}^{- 1}$ allows for arbitrary heteroskedasticity. For balanced panels, b2009 and mw2017 discuss simplifications that arise under homoskedasticity in the cross-section and/or time dimension. However, as noted in the discussion of $\mathbf{B}_{2}$ and $\mathbf{B}_{3}$, the attrition process can induce heteroskedasticity, rendering such simplifications generally invalid.
The bias terms $\mathbf{B}_{1}$ and $\mathbf{B}_{3}$ are of order $\overline{T}^{- 1}$, while $\mathbf{B}_{2}$ is of order $\overline{N}^{- 1}$, where $\overline{T} \coloneqq N^{- 1} \sum_{i = 1}^{N} \sum_{t = 1}^{T} \operatorname{\mathbb{E}}[\delta_{it} \mid \mathcal{C}]$ and $\overline{N} \coloneqq T^{- 1} \sum_{t = 1}^{T} \sum_{i = 1}^{N} \operatorname{\mathbb{E}}[\delta_{it} \mid \mathcal{C}]$. The bias terms are therefore larger, to a degree that depends on the extent of missing data. This is consistent with the results of fw2018 on the asymptotic distribution of traditional fixed effects estimators in unbalanced panels.
For balanced panels, the asymptotic distribution reduces to that derived by mw2017, rendering Assumptions (ref) and (ref) redundant.
For balanced panels, \textcites{mw2015}{mw2017} showed that the profile objective function (ref) can be reformulated as
where $\boldsymbol{\Gamma}(\boldsymbol{\beta})$ is an $N \times T$ matrix with $[\boldsymbol{\Gamma}(\boldsymbol{\beta})]_{it} = y_{it} - \mathbf{x}_{it}^{\prime} \boldsymbol{\beta}$, and $\mu_{r}(\cdot)$ denotes the $r$-th largest eigenvalue. This reformulation is advantageous because it eliminates the need to optimize explicitly over the high-dimensional nuisance parameters $\boldsymbol{\Lambda}$ and $\mathbf{F}$. Moreover, since modern algorithms for symmetric eigenvalue problems are highly optimized, computing $\hat{\boldsymbol{\beta}}$ remains efficient even for large $T$. Estimates of $\boldsymbol{\Lambda}$ and $\mathbf{F}$ are subsequently recovered by decomposing $\widehat{\boldsymbol{\Gamma}} \coloneqq \boldsymbol{\Gamma}(\hat{\boldsymbol{\beta}})$. Specifically, under the normalizing restrictions $\mathbf{F}^{\prime} \mathbf{F}^{\phantom{\prime}} / T = \operatorname{\mathbb{I}}_{R}$ and $\boldsymbol{\Lambda}^{\prime}\boldsymbol{\Lambda}^{\phantom{\prime}}$ diagonal, $\widehat{\mathbf{F}}$ equals the first $R$ eigenvectors of $\widehat{\boldsymbol{\Gamma}}^{\prime} \widehat{\boldsymbol{\Gamma}}$ multiplied by $\sqrt{T}$, and $\widehat{\boldsymbol{\Lambda}} = \widehat{\boldsymbol{\Gamma}} \widehat{\mathbf{F}} / T$.\footnote{Other valid normalizing restrictions are discussed in bn2013. Moreover, if $T > N$, it is computationally more efficient to minimize $(NT)^{- 1} \sum_{r = R + 1}^{N} \mu_{r} \big(\boldsymbol{\Gamma}(\boldsymbol{\beta}) \boldsymbol{\Gamma}(\boldsymbol{\beta})^{\prime}\big)$ and estimate $\widehat{\boldsymbol{\Lambda}}$ as the first $R$ eigenvectors of $\widehat{\boldsymbol{\Gamma}} \widehat{\boldsymbol{\Gamma}}^{\prime}$ multiplied by $\sqrt{N}$ and $\widehat{\mathbf{F}} = \widehat{\boldsymbol{\Gamma}}^{\prime} \widehat{\boldsymbol{\Lambda}} / N$, imposing $\boldsymbol{\Lambda}^{\prime} \boldsymbol{\Lambda}^{\phantom{\prime}} / N = \operatorname{\mathbb{I}}_{R}$, where $\mathbf{F}^{\prime}\mathbf{F}^{\phantom{\prime}}$ is diagonal.}
The estimation procedure for unbalanced panels is motivated by the following decomposition of (ref):
where $\boldsymbol{\Gamma}(\boldsymbol{\beta})$ has missing entries corresponding to unobserved index pairs, i.e., entries are missing whenever $(i, t) \notin \mathcal{D}$. The key idea is that, for a given $\boldsymbol{\beta}$, the observed entries of $\boldsymbol{\Gamma}(\boldsymbol{\beta})$ can be used to estimate $\boldsymbol{\Lambda}$ and $\mathbf{F}$, which in turn permit imputation of the missing entries via $\boldsymbol{\lambda}_{i}^{\prime} \mathbf{f}_{t}^{\phantom{\prime}}$.
We introduce two matrix completion algorithms to accomplish this. Algorithm 1 is the classical Expectation-Maximization (EM) algorithm, originally developed by \textcites{sw1998}{sw2002} for pure factor models. Algorithm 2 follows fll2021, combining nuclear norm regularization with a debiasing step to mitigate regularization bias. The second approach is particularly advantageous for large-scale panel data, as it typically offers substantial computational speed gains over the EM algorithm. Numerical comparisons of both algorithms are provided in Appendix (ref). We then demonstrate how these matrix completion algorithms integrate into the reformulation of the profile objective function proposed by \textcites{mw2015}{mw2017}.
Before introducing the matrix completion algorithms, we adopt the notation of ccs2010 for handling observed and missing data. Let
denote the projection operator onto the subspace of matrices whose support is contained in $\mathcal{D}$, and let $\mathcal{P}_{\mathcal{D}}^{\perp}$ denote its orthogonal complement, defined analogously with $\mathcal{D}$ replaced by its complement. By construction, $\mathcal{P}_{\mathcal{D}}^{\phantom{\perp}}(\mathbf{M}) + \mathcal{P}_{\mathcal{D}}^{\perp}(\mathbf{M}) = \mathbf{M}$ for any $N \times T$ matrix $\mathbf{M}$.
We now present the first matrix completion algorithm.
Algorithm (ref) is the classical approach to missing data for pure factor models. Heuristically, Step 1 is the E-step, where missing entries are imputed using current parameter estimates, and Step 2 is the M-step, which applies eigenvalue decomposition to the completed data, motivated by the decomposition in (ref). Despite its longstanding use in empirical work, the formal asymptotic properties of this algorithm were established only recently by jms2021. As noted by fll2021, the EM algorithm can be computationally burdensome for large-scale panel data relative to modern regularization-based alternatives.
Hence, we next present the regularized matrix completion algorithm proposed as Algorithm 5 in fll2021.
Algorithm (ref) is a modern matrix completion approach combining nuclear norm regularization with a post-estimation debiasing step. Step 1 implements the SOFT-IMPUTE algorithm of mht2010 to solve the nuclear norm penalized optimization problem. Unlike Algorithm (ref), which imposes a “hard” rank constraint by retaining only the first $R$ singular values, Algorithm (ref) restricts the rank implicitly via the tuning parameter $\nu$. Since the nuclear norm is the convex envelope of the rank operator, Algorithm (ref) is expected to outperform Algorithm (ref) in terms of computational speed in many settings, particularly high-dimensional ones (see mht2010 for details). Steps 2 through 4 implement the two-step least squares debiasing procedure of \textcites{chlz2019}{chlz2023} to mitigate regularization bias.
For unbalanced panels, we adapt the profile objective function to accommodate missing observations by using the completed matrix:
where $\boldsymbol{\Gamma}^{\ast}(\boldsymbol{\beta})$ is the completed matrix obtained after convergence of Algorithm (ref) or (ref). After obtaining $\hat{\boldsymbol{\beta}}$ by minimizing (ref), $\widehat{\boldsymbol{\Lambda}}$ and $\widehat{\mathbf{F}}$ are recovered by decomposing $\widehat{\boldsymbol{\Gamma}}^{\ast} \coloneqq \boldsymbol{\Gamma}^{\ast}(\hat{\boldsymbol{\beta}})$.
To solve this minimization problem efficiently, we recommend a Quasi-Newton method (e.g., BFGS) with an analytical gradient. Let $\check{\boldsymbol{\beta}}$ denote a trial value, and let $\check{\boldsymbol{\Lambda}}$ and $\check{\mathbf{F}}$ be the estimates of $\boldsymbol{\Lambda}$ and $\mathbf{F}$ obtained from the decomposition of $\check{\boldsymbol{\Gamma}}^{\ast} \coloneqq \boldsymbol{\Gamma}^{\ast}(\check{\boldsymbol{\beta}})$. The analytical gradient is
In our simulations and empirical applications, BFGS typically converges in few iterations and substantially reduces the computational overhead introduced by the iterative matrix completion task.
Our estimation procedure is based on a reformulation of the objective function in b2009 for unbalanced panels. Convergence of his alternating estimation algorithm was recently established by sww2026, and we expect their results to carry over to our setting.
Because the rank constraint renders the optimization problem in (ref) non-convex (see mw2026), the choice of starting values is critical. Following the intuition in sw2016b, one could construct initial estimates from balanced sub-panels. However, this approach requires sufficiently large sub-panels and still necessitates testing multiple starting guesses.
To overcome these limitations, we recommend initializing the optimization with the nuclear norm minimizing estimator of mw2026:
where $\sigma_{j}(\cdot)$ denotes the $j$-th largest singular value. The key advantage of (ref) is that its objective function is convex. Although mw2026 show that this estimator is consistent only at rate $\sqrt{\min(N, T)}$ (rather than the rate $\sqrt{NT}$ of (ref)), it provides a reliable and computationally efficient starting guess for minimizing (ref). We recommend solving (ref) using a Quasi-Newton method with analytical gradient
where $\mathcal{P}_{\mathcal{D}}^{\phantom{\perp}}(\boldsymbol{\Gamma}(\check{\boldsymbol{\beta}})) = \check{\mathbf{U}} \check{\boldsymbol{\Sigma}} \check{\mathbf{V}}^{\prime}$ is the singular value decomposition of $\mathcal{P}_{\mathcal{D}}^{\phantom{\perp}}(\boldsymbol{\Gamma}(\check{\boldsymbol{\beta}}))$ with $\check{\boldsymbol{\Sigma}} = \operatorname{\text{diag}}(\check{\sigma}_{1}, \ldots, \check{\sigma}_{\min(N, T)})$.
We obtain estimators for $\mathbf{W}$, $\boldsymbol{\Omega}$, $\mathbf{B}_{1}$, $\mathbf{B}_{2}$, and $\mathbf{B}_{3}$ by forming sample analogues, i.e., by dropping expectations and substituting the corresponding estimators for $\boldsymbol{\beta}$, $\boldsymbol{\Lambda}$, and $\mathbf{F}$. Let $L$ denote a bandwidth parameter for the truncation kernel of nw1987, depending on the sample size. Then,
where $\hat{p}_{itt^{\prime}}^{f} \coloneqq \hat{\mathbf{f}}_{t}^{\prime} (\sum_{t^{\prime} \in \mathcal{T}_{i}} \hat{\mathbf{f}}_{t^{\prime}}^{\phantom{\prime}} \hat{\mathbf{f}}_{t^{\prime}}^{\prime})^{- 1} \hat{\mathbf{f}}_{t^{\prime}}$, $\hat{\xi}_{it}^{\dagger} \coloneqq \hat{\boldsymbol{\lambda}}_{i}^{\prime} (\widehat{\boldsymbol{\Lambda}}^{\prime} \widehat{\boldsymbol{\Phi}} \widehat{\boldsymbol{\Lambda}})^{- 1} (\widehat{\boldsymbol{\Lambda}}^{\prime} \widehat{\boldsymbol{\Lambda}}) \hat{\mathbf{f}}_{t}^{\phantom{\prime}}$, $\widehat{\boldsymbol{\Phi}} \coloneqq \mathcal{P}_{\mathcal{D}}(\widehat{\boldsymbol{\Lambda}} \widehat{\mathbf{F}}^{\prime}) (\mathcal{P}_{\mathcal{D}}(\widehat{\boldsymbol{\Lambda}} \widehat{\mathbf{F}}^{\prime}))^{\prime}$, $\widehat{\mathbf{A}}^{\cdot} \coloneqq (\hat{\mathbf{a}}_{1}^{\cdot}, \ldots, \hat{\mathbf{a}}_{T}^{\cdot})^{\prime}$, $\widehat{\mathbf{C}}^{\cdot} \coloneqq (\hat{\mathbf{c}}_{1}^{\cdot}, \ldots, \hat{\mathbf{c}}_{N}^{\cdot})^{\prime}$,
A debiased estimator for $\boldsymbol{\beta}$ is then constructed as
such that
The factor $\mathcal{T}_{i} / (\mathcal{T}_{i} - j)$ in $\widehat{\mathbf{B}}_{1}$ is a finite-sample adjustment proposed by fw2016. Following fw2018, we use $\sqrt{n}$ rather than $\sqrt{NT}$ as the normalizing factor in (ref) to improve finite-sample approximation. The uncorrected estimator $\hat{\boldsymbol{\beta}}$ is obtained using the algorithms of Section (ref).
To construct $\tilde{\boldsymbol{\beta}}$, we require a computationally feasible method for the residuals defined in (ref), (ref), and (ref). Consider an arbitrary $n$-dimensional vector $\mathbf{v}$. The minimization over $\mathbf{A}$ in (ref),
is separable across $t$. For each $t \in \{1, \ldots, T\}$, the solution reduces to a cross-sectional regression:
The corresponding residuals are
The same argument applies to the minimization problem in (ref).\footnote{The separability of both minimization problems is also exploited in the two-step least squares debiasing procedure of \textcites{chlz2019}{chlz2023} (see Steps 3 and 4 of Algorithm (ref)).} For each $i \in \{1, \ldots, N\}$, the solution reduces to a time-series regression, with residuals
For the minimization problem in (ref), no closed-form expressions analogous to (ref) and (ref) are available.
We propose a novel algorithm based on the Method of Alternating Projections (MAP, see \cites{vn1949}{vn1950}{h1962}) as a computationally feasible method to compute the residuals (ref).\footnote{Our algorithm adapts s2018, who introduced MAP as a powerful tool for demeaning variables in the optimization of fixed effects estimators for nonlinear models with multi-way fixed effects, such as binary choice models with individual and time effects. s2020 (Chapter 3) and cs2019 noted the usefulness of this approach for unbalanced panel data. MAP is particularly appealing because it computes residuals from complex regressions (including unbalanced, weighted, and multi-way fixed effects) by alternating between one-way fixed effects demeaning steps.} Let $\hat{\mathbf{v}}^{\lambda}$ be the $n$-dimensional vector with entries (ref), and $\hat{\mathbf{v}}^{f}$ the $n$-dimensional vector with entries (ref). We define two orthogonal projection operators, $\mathcal{M}_{\hat{\lambda}}(\mathbf{v})$ and $\mathcal{M}_{\hat{f}}(\mathbf{v})$, such that $\mathcal{M}_{\hat{\lambda}}(\mathbf{v}) = \hat{\mathbf{v}}^{\lambda}$ and $\mathcal{M}_{\hat{f}}(\mathbf{v}) = \hat{\mathbf{v}}^{f}$.
Algorithm (ref) iterates between orthogonal projections onto two closed subspaces and converges strongly to the projection $\hat{\mathbf{v}}^{\lambda f}$, i.e., to the residuals (ref), as established by \textcites{vn1949}{vn1950}. The linear rate of convergence was first proved by a1950. Acceleration techniques are discussed in er2011, among others.
We now summarize how the components of this section are combined to conduct inference on $\boldsymbol{\beta}$ using the debiased interactive fixed effects estimator for unbalanced panels.
If certain bias terms are not required, for example, when all regressors are strictly exogenous (so that $\mathbf{B}_{1} = \mathbf{0}_{K}$), the corresponding estimates, $\widehat{\mathbf{B}}_{1}$ in the example, can be omitted from Step 5. Since choosing an appropriate bandwidth $L$ for $\widehat{\mathbf{B}}_{1}$ is non-trivial, \textcites{fw2016}{fw2018} recommend a sensitivity analysis reporting estimates across different values of $L$.
Algorithm (ref) applies to balanced panels as well, by removing Step 2 and replacing the completed matrix $\boldsymbol{\Gamma}^{\ast}(\boldsymbol{\beta})$ with $\boldsymbol{\Gamma}(\boldsymbol{\beta})$.
\textcites{b2009}{mw2017} derived their results under the assumption that the number of factors is known. In practice, this assumption is often very unlikely unless economic theory provides a clear prediction about the number of factors. Even in that case, it may be necessary to support the theoretical prediction with additional empirical evidence. We therefore need a reliable method to estimate the number of factors. We denote the true number of factors by $R^{0}$.
For pure factor models, i.e., (ref) without additional regressors, there is an extensive literature on estimating the number of factors (see, among others, be1992, bn2002, hl2007, abc2010, o2010, ah2013, and do2019). As pointed out by b2009,
is essentially a pure factor model. Thus, given an estimator for $\boldsymbol{\beta}$ such that the estimation error $\mathbf{x}_{it} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})$ is asymptotically negligible, the number of factors can be estimated consistently using methods developed for pure factor models (see {\citeresetb2009} Remark 5 and the corresponding appendix). Since mw2015 show that the interactive fixed effects estimator is at least $\sqrt{\min(N, T)}$-consistent for any $R \geq R^{0}$, the initial estimate of $\boldsymbol{\beta}$ should be based on a sufficiently large value of $R$.
We consider the estimators of \textcites{bn2002}{o2010}{ah2013}{do2019}. Specifically, we apply them to $\widehat{\boldsymbol{\Gamma}}$, where $\boldsymbol{\beta}$ is estimated using $R = \overline{R}$ and $\overline{R}$ is a known upper bound on the number of factors. bn2002 proposes model selection criteria that minimize the sum of squared residuals plus a penalty for the number of estimated parameters. \textcites{o2010}{ah2013}{do2019} segment the eigenvalue spectrum of the sample covariance of $\widehat{\boldsymbol{\Gamma}}$ to identify a cut-off between the common factors and the noise from the idiosyncratic error term. o2010 proposes the edge distribution estimator (ED), based on differences of consecutive eigenvalues. ah2013 proposes using ratios (ER) and growth rates (GR) instead of differences. be1992 proposes a specific version of parallel analysis (PA), which compares eigenvalues to those obtained from independent data to identify a cut-off between common factors and noise. Independent data are constructed by permuting each column of $\widehat{\boldsymbol{\Gamma}}$, which preserves the marginal variances while destroying the correlation pattern induced by the common factors. Theoretical justification for PA was recently provided by d2020.
For unbalanced panels, we follow jms2021 and apply the estimators to $\mathcal{P}_{\mathcal{D}}^{\phantom{\perp}}(\widehat{\boldsymbol{\Gamma}}^{\ast}) / (1 - \psi) = \mathcal{P}_{\mathcal{D}}^{\phantom{\perp}}(\boldsymbol{\Gamma}(\hat{\boldsymbol{\beta}})) / (1 - \psi)$ rather than $\widehat{\boldsymbol{\Gamma}}^{\ast}$, where $\psi \coloneqq 1 - n / (NT)$ is the share of missing observations.
We use Monte Carlo simulations to analyze the finite-sample properties of the debiased estimator $\tilde{\beta}$, defined in (ref), in the presence of missing data. Specifically, we compare relative biases (Bias), average ratios of standard errors to standard deviations (Ratio), and empirical sizes of $z$-tests with a 5% nominal size (Size) across different shares of missing data ($\psi$) and relative to the balanced panel case. We use Algorithm (ref) as matrix completion procedure for unbalanced panels. Because the number of factors is typically unknown, we also compare different estimators for the number of factors. Specifically, we consider the estimators of \textcites{bn2002}{o2010}{ah2013}{do2019}. Of the information criteria introduced by bn2002, we focus on $\text{IC}_{2}$ and $\text{BIC}_{3}$, which are also used in o2010 and ah2013. Performance is assessed by comparing the average estimated number of factors.
We follow mw2017 and consider an AR(1) model with $R = 1$ factor,
The idiosyncratic error term $e_{it}$ is homoskedastic with fat tails. Specifically, $e_{it}$ is drawn independently and identically from the $t$-distribution with five degrees of freedom. The factor structure is constructed from $\lambda_{i} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(1, 1)$ and $f_{t} = \rho \, f_{t-1} + u_{t}$, where $u_{t} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(0, (1 - \rho^2) \sigma^2)$ and $\rho = \sigma = 0.5$. We discard the first $1{,}000$ time periods to ensure that the simulated data are drawn from the stationary distribution of the model. All random variables are redrawn in each replication, and all results are based on $1{,}000$ replications.
We consider three shares of missing data, $\psi \in \{0, 0.2, 0.4\}$, with $\psi = 0$ corresponding to a balanced panel. The total sample size satisfies $n = NT(1 - \psi)$. As implied by the results in Section (ref), the biases shrink with $\overline{N}$ and $\overline{T}$. To ensure comparability across values of $\psi$, we therefore select $N$ and $T$ so that both $\overline{N}$ and $\overline{T}$ remain constant, setting $N = \overline{N} / (1 - \psi)$ and $T = \overline{T} / (1 - \psi)$. We consider panels with $\overline{N} = 100$ and $\overline{T} \in \{5, 10, 20, 40, 80\}$, and AR(1) models with $\beta = 0.3$ and $\beta = 0.9$.
Figure (ref) illustrates the missing data pattern for $\overline{N} = 100$ and $\overline{T} = 20$ across different values of $\psi$.
The pattern is taken from cs2019. All units are divided into two types. Type 1 consists of $N_{1} = 2 \psi N$ units observed over $T_{1} = T / 2$ consecutive time periods. The remaining $N_{2} = N - N_{1}$ units are Type 2 and are observed over the entire time horizon, i.e., $T_{2} = T$. The initial period is drawn uniformly at random from $\{0, 1, \ldots, T - T_{1}\}$. All unbalanced data sets are generated from initially balanced panels. Whether unit $i$ is Type 1 or Type 2 is determined by the value of $\lambda_{i}$: units with the lowest values of $\lambda_{i}$ are assigned to Type 1. Observations are therefore not missing completely at random but are conditionally missing at random. Note also that the missing probabilities are homogeneous across $i$ but heterogeneous across $t$: they are lowest for time periods at the beginning or end of the time series and highest for time periods near $T / 2$.
Table (ref) presents the simulation results for $\tilde{\beta}$.
Although $e_{it}$ is homoskedastic and the missing probabilities are homogeneous across $i$, we do not exploit this information. Instead, we apply the debiased estimator and its covariance matrix estimator exactly as described in Section (ref), correcting for all three bias terms and using a covariance estimator that is robust to arbitrary heteroskedasticity. This approach yields a more realistic assessment of finite-sample performance in practice, where the true data-generating process is unknown and heteroskedasticity-robust inference is standard. The bandwidth parameter $L$ is taken from Table 1 of mw2017. For both $\beta = 0.3$ and $\beta = 0.9$, the biases, ratios, and sizes are similar to those in the balanced case, regardless of the share of missing data. Overall, the finite-sample performance of $\tilde{\beta}$ in unbalanced panels is well predicted by our theory.
Table (ref) presents the simulation results for the various estimators of the number of factors, $\widehat{R}$.
The initial estimator uses $R = \overline{R}$, with $\overline{R} = 2$ for $\overline{T} = 5$, $\overline{R} = 5$ for $\overline{T} = 10$, and $\overline{R} = 10$ for $\overline{T} \in \{20, 40, 80\}$.\footnote{Our choice of $\overline{R}$ differs from studies such as \textcites{bn2002}{o2010}{ah2013}, which hold $\overline{R}$ fixed regardless of the sample size.} For $\psi > 0$, we apply the estimators to $\mathcal{P}_{\mathcal{D}}^{\phantom{\perp}}(\boldsymbol{\Gamma}(\hat{\beta}_{\bar{R}})) / (1 - \psi)$ as suggested by jms2021, where $\hat{\beta}_{\bar{R}}$ denotes the initial estimator with $R = \overline{R}$. For ER and GR, we use the mock eigenvalue of ah2013 to allow for the possibility of selecting zero factors. For sufficiently large $\overline{T}$, all estimators perform similarly across different shares of missing data. However, performance differs between $\beta = 0.3$ and $\beta = 0.9$. In the low-persistence setting with sufficiently large $\overline{T}$, all estimators recover the correct number of factors, $R = 1$, nearly perfectly. In the high-persistence setting, only ER and GR achieve good performance. The results also suggest that the estimation error in $\boldsymbol{\beta}$ is asymptotically negligible, and they support the conjecture of mw2015 that their main results extend beyond the case of independent and identically normally distributed errors.
Tables (ref) and (ref) in the Online Supplement report the simulation results for $\tilde{\beta}$ and $\widehat{R}$ using Algorithm (ref) in place of Algorithm (ref) as the matrix completion procedure. The tuning parameter $\nu$ is selected via the plug-in approach of \textcites{chlz2019}{chlz2023}. The results are virtually identical to those reported here. Additional simulation results for a static panel data model with one regressor, two factors, and other missing data patterns can be found in the Online Supplement (ref).
The effect of democracy on economic growth remains a highly debated topic among economists. anrr2019 provide evidence that democratization has a substantial positive impact on GDP per capita. Using annual data from 175 countries observed between 1960 and 2010, their main findings suggest a long-run effect of about 20%. The dataset they construct is well suited for our purposes: it is naturally unbalanced, spans a long time horizon, and contains several unobserved common shocks triggered by technological progress and financial crises.\footnote{The data are part of the \href{https://www.journals.uchicago.edu/doi/suppl/10.1086/700936}{replication package} provided by the authors.}
The sample consists of $6{,}934$ observations, of which $3{,}558$ are classified as democratic. Of the 175 countries, 88 transition between democracy and non-democracy or vice versa. Average GDP, measured in year-2000 dollars, is $8{,}150$ for democratic and $2{,}074$ for non-democratic countries. A total of 71 countries are observed over the entire time horizon; on average, the dataset covers 136 countries and 40 years. The fraction and pattern of missing data are comparable to the configuration of our simulation study with $\psi = 0.2$.
Figure (ref) illustrates the evolution of GDP per capita around democratization for transitioning countries relative to persistently non-democratic ones.
The pre-transition dip and subsequent recovery suggest that the timing of democratization is endogenous to past GDP shocks. Failing to account for these dynamics by including sufficient lags of GDP could bias estimates, incorrectly attributing natural mean reversion following a crisis to the effect of democracy. To address this, anrr2019 adopt the following dynamic panel specification:
where $y_{it}$ is the natural logarithm of GDP per capita of country $i$ at time $t$, $D_{it}$ is a democracy indicator, $\alpha_{i}$ and $\delta_{t}$ denote country and year fixed effects, and $u_{it}$ is an idiosyncratic error term. $\boldsymbol{\beta} = (\theta, \boldsymbol{\gamma}^{\prime})^{\prime}$ are the parameters of interest. This specification also allows us to distinguish between the short-run effect, $\theta$, and the long-run effect of democratization, $\phi(\boldsymbol{\beta}) \coloneqq \theta / (1 - \sum_{j = 1}^{p} \gamma_{j})$.
In contrast to anrr2019, we further decompose the error term into a factor structure $\boldsymbol{\lambda}_{i}^{\prime} \mathbf{f}_{t}^{\phantom{\prime}}$ and a residual idiosyncratic component $e_{it}$, i.e., $u_{it} = \boldsymbol{\lambda}_{i}^{\prime} \mathbf{f}_{t}^{\phantom{\prime}} + e_{it}$. This decomposition captures unobserved common shocks ($\mathbf{f}_{t}$) that simultaneously affect GDP growth and democratization in heterogeneous ways ($\boldsymbol{\lambda}_{i}$). Following anrr2019, we report results for $p \in \{1, 2, 4\}$, noting that $p = 4$ is the authors' preferred specification for modeling the GDP dynamics that follow a transition. As in the simulation study of Section (ref), we use Algorithm (ref) as the matrix completion procedure.
To reduce the number of parameters during optimization, we project out the country and time fixed effects before estimating $\boldsymbol{\beta}$:
where a dot denotes variables after projecting out country and time fixed effects. For example,
The remaining residuals, $\dot{D}_{it}$, $\dot{y}_{i(t - 1)}$, $\dot{y}_{i(t - 2)}$, $\dot{y}_{i(t - 3)}$, and $\dot{y}_{i(t - 4)}$, are defined analogously. These residuals can be computed using the MAP algorithm (Algorithm (ref)) presented in Section (ref) with $R = 1$, $\hat{\lambda}_{i1} = 1$ for all $i \in \{1, \ldots, N\}$, and $\hat{f}_{t1} = 1$ for all $t \in \{1, \ldots, T\}$.\footnote{In this case, Algorithm (ref) reduces to the algorithm proposed by g2013_paper, which is also used in popular fixed effects estimation software such as $\textit{lfe}$ g2013_software and reghdfe c2016.}
For valid inference, the true number of factors must be known, or at least consistently overestimated. Since the true number is unknown, we proceed as follows. We estimate each specification with $R = 5$ to obtain $\mathcal{P}_{\mathcal{D}}^{\phantom{\perp}}(\boldsymbol{\Gamma}( \hat{\boldsymbol{\beta}})) / (1 - \psi)$, where $[\boldsymbol{\Gamma}( \hat{\boldsymbol{\beta}})]_{it} \coloneqq \hat{\theta} \, \dot{D}_{it} + \sum_{j = 1}^{p} \hat{\gamma}_{j} \, \dot{y}_{it - j}$. We then apply the estimators of \textcites{be1992}{bn2002}{o2010}{ah2013} to estimate the number of factors. For ER and GR, we use the mock eigenvalue of ah2013 to accommodate the possibility of zero factors.
Table (ref) summarizes the results.
The estimates are nearly identical for $p > 1$, i.e., for specifications with flexible dynamics. For $p = 2$, most estimators select one common factor; for $p = 4$, half do so. The ER and GR results for $p = 4$ should be interpreted with caution, as they depend on the definition of the mock eigenvalue, which has multiple formulations (see {\citeresetah2013}). Estimates for $p = 1$ vary substantially across estimators, ranging from one to five common factors, which may reflect insufficiently specified dynamics.
Figure (ref) displays the singular values of the pure factor models alongside those of the permuted versions.\footnote{More precisely, we randomly shuffle each column of $\mathcal{P}_{\mathcal{D}}^{\phantom{\perp}}(\boldsymbol{\Gamma}(\hat{\boldsymbol{\beta}}))$ and compute the maximum singular value across 199 randomized samples, multiplied by 1.05. This scaling factor is suggested by do2019, whose reasoning is that a factor only marginally exceeding what would be expected from pure noise should not be included in the model.}
We focus on the flexible dynamic specifications with $p > 1$. The gap between the first and second singular values explains why most estimators that decompose the eigenvalue spectrum select one common factor. Comparing the spectra with those of the permuted data, however, reveals that this common factor has explanatory power beyond what noise alone would generate, even if it accounts for only a small share of total variance. In light of the finding by mw2015 that overestimating the number of factors is preferable to underestimating it, $R = 1$ is our preferred choice for $p > 1$.
Table (ref) summarizes our results.
Following anrr2019, we report results for the fixed effects estimator (FE), the Arellano-Bond estimator (AB, see ab1991), and the Hahn-Hausman-Kuersteiner estimator (HHK, see hhk2004). However, rather than the uncorrected FE estimator used by anrr2019, we report results from a debiased estimator with bandwidth $L = 5$ to correct for feedback bias.\footnote{Let $\hat{\boldsymbol{\beta}}_{\text{FE}}$ denote the uncorrected FE estimator and let $\dot{\mathbf{x}}_{it} \coloneqq (\dot{D}_{it}, \dot{y}_{i(t-1)}, \ldots)^{\prime}$. Then, the debiased FE estimator is constructed as
A similar estimator was also used by ccf2019 for the same empirical illustration on a balanced subset of the data. Our debiased FE estimator follows the formulation of fw2018, augmented with the finite-sample adjustment described in Section 2.3. The key difference between ccf2019 and fw2018 is that the latter uses residualized regressors in place of raw regressors. We adopt the formulation of fw2018 as it is more closely aligned with the bias expressions underlying our debiased IFE estimator.} In addition, we report results for three debiased interactive fixed effects estimators (IFE) with $R \in \{1, 2, 3\}$. We correct for both feedback bias and biases induced by heteroskedasticity, with bandwidth $L = 5$. We report estimates and standard errors for the short- and long-run effects of democratization and the persistence of GDP processes.
All estimators indicate strong and significant GDP persistence across all specifications. The democracy coefficients from FE and IFE are significant at the 5% level throughout, whereas those from AB and HHK are significant only for $p = 4$. Focusing on the preferred specification $p = 4$, the estimators used by the authors imply short-run effects of democratization between 0.725% and 1.178%, and long-run effects between 16.448% and 25.032%. After controlling for additional time-varying unobserved heterogeneity, however, both effects are substantially smaller. Our preferred specification, IFE with $R = 1$, yields short- and long-run estimates of 0.519% and 12.334%, respectively.\footnote{Additional sensitivity checks are provided in Appendix (ref). In particular, all IFE estimates are remarkably stable across bandwidth choices $L \in \{1, \ldots, 8\}$, and across different values of $R$, with the exception of $p = 1$.}
In summary, we find further support for the “democracy does cause growth” hypothesis of anrr2019. Controlling for time-varying unobserved heterogeneity via the interactive fixed effects estimator yields results that are qualitatively similar to those of the original authors. In the preferred specification $p = 4$, comparing HHK to IFE with $R = 1$ shows that both the short-run and long-run effects of democratization are roughly halved.
Although our analysis focuses on the interactive fixed effects (IFE) estimator of b2009, we briefly discuss three related estimators for which our findings and algorithms may also prove useful. First, in the presence of endogenous regressors, \textcites{mw2017}{msw2018} propose a minimum distance estimator in the spirit of \textcites{ch2006}{ch2008}. Second, because the IFE objective function is generally nonconvex, mw2026 proposes an alternative estimator that replaces the potentially difficult nonconvex optimization problem with a convex one. Third, cfw2021 propose an estimator for nonlinear parametric single-index models with interactive effects, such as logit, probit, ordered probit, and Poisson models.
\noindentMinimum Distance Estimator. Suppose that $\mathbf{x}_{it}$ can be decomposed into $K_{1}$ endogenous and $K_{2}$ exogenous regressors, so that $K = K_{1} + K_{2}$. We use superscripts to distinguish between endogenous and exogenous regressors. Let $\mathbf{z}_{it} = (z_{1, it}, \ldots, z_{M, it})^{\prime}$ be a vector of excluded exogenous instruments, where $M \geq K_{1}$. mw2017 suggest the following minimum distance estimator. In the first step, an estimator for $\boldsymbol{\beta}^{\text{end}}$ is obtained by
where $\hat{\boldsymbol{\pi}}(\boldsymbol{\beta}^{\text{end}})$ is the IFE estimator of
and $\boldsymbol{\Sigma}$ is a positive definite $M \times M$ weighting matrix. At the true value of $\boldsymbol{\beta}^{\text{end}}$, the instrumental variable moment conditions imply $\boldsymbol{\pi} = \mathbf{0}_{M}$. In the second step, $\hat{\boldsymbol{\beta}}^{\text{exo}}$ is the IFE estimator of
The properties of the minimum distance estimator are studied in msw2018, where the authors extend the random coefficient demand model of blp1995 to include interactive fixed effects, thereby accounting for unobserved product-market-specific heterogeneity, such as perceived utility from advertising at the product-market level. Under assumptions similar to those of mw2017, the authors establish consistency and derive the asymptotic distribution of the estimator. lmw2012 apply the same estimator to address measurement error in the dependent variable of dynamic interactive fixed effects models.
\noindentNuclear Norm Regularized Estimator. mw2026 show that the rank constraint on the factor structure renders the optimization problem nonconvex. They propose two alternative estimators based on a convex relaxation of this constraint. The nuclear norm minimizing estimator was already presented in (ref); the second estimator uses nuclear norm regularization. mw2026 establish consistency for both estimators, but only at the rate $\sqrt{\min(N, T)}$. To recover the properties of the IFE estimator, they suggest estimating the number of factors from $\mathcal{P}_{\mathcal{D}}^{\phantom{\perp}}(\boldsymbol{\Gamma}(\hat{\boldsymbol{\beta}}^{\star}))$ and then applying an iterative post-estimation routine. After a finite number of iterations, the estimator attains the same limiting distribution as the IFE estimator.
\noindentEstimator for Nonlinear Factor Models. Suppose the outcome variable is generated by
where $f(\cdot)$ is a known density, such as the logistic density. To maximize the corresponding log-likelihood, cfw2021 propose an EM-type optimization algorithm.
The assumption that unobserved heterogeneity is constant over time is often too restrictive. In panels that span a long time horizon, such as macroeconomic country panels, it is implausible that a global shock affects all units equally. Interactive fixed effects estimators offer researchers a flexible way to accommodate this form of heterogeneity (see, among others, hnr1988, p2006, and b2009). These panels are, however, often naturally unbalanced. Although b2009 proposed an estimation algorithm for this case, the practical aspects of inference remained unclear. Drawing on insights from fw2018 and extending mw2017, we derive the asymptotic distribution of b2009's interactive fixed effects estimator for unbalanced panels, thereby establishing a foundation for inference in this practically relevant setting. We also develop a novel algorithm to compute the residualized variables required for estimating the bias terms and the covariance matrix.
Our findings and algorithms may further prove useful for related estimators, including the minimum distance estimator of \textcites{mw2017}{msw2018}, the nuclear norm estimator of mw2026, and the estimator for nonlinear factor models of cfw2021.
\printbibliography