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.
162,920 characters · 22 sections · 103 citation commands
Large Dimensional Latent Factor Modeling with Missing Observations and Applications to Causal Inference
\onehalfspacing
Large dimensional panel data with missing entries are prevalent. In causal panel data, the main focus is to estimate the unobserved potential outcomes. In financial data, stock returns can be missing before a company is listed, after its bankruptcy, or because of illiquidity. In macroeconomic datasets, panel data might be collected at different frequencies or not for all geographical locations resulting in missing entries. In the famous Netflix challenge, a majority of users' ratings for films are missing. Estimating missing entries in panel data is a fundamental problem with applications in social science, statistics, and computer science.
This paper develops the inferential theory for latent factor models estimated from large dimensional panel data with missing observations. We propose a novel and easy-to-use approach to estimate a latent factor model by applying principal component analysis (PCA) to an adjusted covariance matrix, which is estimated from partially observed panel data. We derive the asymptotic normal distribution for the estimated factors, loadings, and imputed values. The key application is to estimate counterfactual outcomes for causal inference. The unobserved control group is modeled as missing values, which are inferred from the latent factor model. The inferential theory for the imputed values allows us to test for individual treatment effects at a particular time. This granular test is of practical relevance because we learn not only for whom but also when a treatment is effective.
The inferential theory for latent factor models with missing data is important for a number of reasons. First, we show how to consistently impute the missing observations in a large dimensional panel data set, which can then be used as an input for other applications. Our confidence intervals for the imputed values can serve as a decision criterion if the imputed data should be used. Second, the distribution of the missing observations can actually be the object of interest itself. For example, the imputed values serve as the synthetic control in causal inference for which we need an asymptotic distribution theory. The inferential theory is key for deriving test statistics for treatment effects. Last but not least, we provide the complete inferential theory for the latent factors themselves, which is relevant when the factors are the object of interest and are used as input for other applications.
Our method is very simple to adopt and works under general assumptions. We provide an “all-purpose” estimator that performs well under all empirically relevant missing patterns and only assumes a general approximate factor model. Our estimation consists of two simple steps, where we first apply PCA to a re-weighted covariance matrix to obtain the loadings and, in a second step, run a regression on these loadings using only the observed units to obtain the factors. The missing entries are estimated by the common components of the factor model. Importantly, our estimator does not require the estimation of the observation pattern itself. In some cases, we might have additional information about the missing pattern. We provide a modification of our estimator that can take advantage of a probabilistic model of the missing pattern and use an inverse probability weight in the second step regression to obtain the factors. It is inspired by the inverse propensity weighted regression from causal inference that enjoys the doubly-robust property, meaning the estimator is robust to some form of omitted variable bias. Our probability weighted estimator also has similar desirable robustness properties when we omit latent factors, but it is generally less efficient than our all-purpose estimator.
Our framework stands out by the very general patterns of missing observations that it can accommodate. We cover the common scenarios of missing at random or a simultaneous/staggered treatment adoption, where the treatment cannot be removed once implemented. Importantly, the missing pattern can depend in a general way on the unobserved factor loadings or unit-specific features. Hence, the observations can be missing because of how the units are exposed to the latent factors. Our simple all-purpose estimator does not require us to explicitly model this relationship, but takes it automatically into account. In the case of the propensity weighted estimator, we provide feasible estimators of the probability weights that result in the same distribution as the population weights.
Deriving the inferential theory under these general conditions is a challenging problem. The missing observations have a complex effect on the asymptotic covariance matrix of the imputed entries. In particular, the asymptotic variance has an additional variance correction term compared with the fully observed panel. This term results in a larger asymptotic variance than in the fully observed case. The variance correction term arises because, in a panel with missing observations, we take averages over a different number of time periods for the different entries in the estimated covariance matrix. The variance correction term is larger if the observation pattern has many missing entries, or if it deviates more from a missing at random scheme. The propensity weighted estimator has a similar asymptotic distribution structure as our all-purpose estimator but in general a larger variance.
Our work contributes to three distinct fields: large dimensional factor modeling, matrix completion, and causal inference. First, we extend the inferential theory of latent factors to large dimensional data with general patterns in missing entries. Second, matrix completion methods impute missing entries under the assumption of a low-rank structure, which is corrupted with noise. We provide confidence intervals for the imputed values. Lastly, the key question in causal inference is the estimation of counterfactual outcomes, i.e., what would have been the outcome if a unit had not been treated or if a unit had been treated. The unobserved counterfactual outcome can naturally be formulated as a missing observation problem. We are the first to provide a test for the point-wise treatment effect that can be heterogeneous and time-dependent under general adoption patterns where the units can be affected by unobserved factors:
This paper works under the framework of an approximate latent factor structure where both the cross-section dimension and time-series dimension are large. When the data is fully observed, bai2002determining show that the factor model can be estimated with PCA applied to the covariance matrix of the data. bai2003inferential and fan2013large derive the consistency and asymptotic normality of the estimated factors, loadings and common components. Extensions of latent factor models with fully observed data include adding observable factors in bai2009panel, sparse and interpretable latent factors in pelgerxiongsparse2020, time-varying loadings in fan2016projected and pelger2018state, high-frequency estimation in pelger2019large and including additional moments to estimate weak factors as in lettaupelger2019. When a panel has missing entries, a common approach is to estimate the factor model from a subset of the data for which a balanced panel is available. This approach has two drawbacks: First, it is, in general, less efficient as our approach makes use of all the data. Second, it can lead to a biased estimate if the data is not missing at random.
The inferential theory of large dimensional factor models with missing observations is an active area of research. Our paper is most closely related to the recent papers by jin2020factor, bai2019matrix, and cahan2021factor. The papers differ in the algorithms to impute the missing observations, the generality of the missing patterns, and the proportion of required observed entries relative to the missing entries. Our main results are derived under the assumption that entries are observed at the same rate as missing entries, but we show that this assumption can be considerably relaxed. Importantly, in contrast to the other papers, our framework allows the missing pattern to depend on unit-specific features and to test for an individual treatment effect at any time for any cross-section unit or a weighted treatment effect. {This is exactly what we need for the main application in causal inference.} jin2020factor provide the inferential theory for the estimated factor model with the expectation-maximization (EM) algorithm under the assumption of randomly missing values. This is a major advance in the literature on using the EM algorithm to impute missing values on cross-sectional data rubin1976inference,dempster1977maximum.\footnote{stock2002macroeconomic,banbura2014maximum,negahban2012restricted propose to use EM algorithms to estimate the factor model from panel data with missing observations. giannone2008nowcasting,doz2011two,jungbacker2011maximum,stock2016dynamic propose to use the state-space framework and Kalman Filtering to estimate the factor model with missing observations. gagliardini2019diagnostic propose a simple diagnostic criterion for an approximate factor structure in large (unbalanced) panel data sets.} bai2019matrix provide the inferential theory for the factor-based imputed values based on the innovative idea of shuffling rows and columns such that there exist fully observed TALL and WIDE blocks for estimating the factor model. Their TALL-WIDE algorithms involves two applications of principal components on the two fully observed blocks. cahan2021factor propose the TALL-PROJECT estimator that extends the TALL-WIDE estimator by first using only the fully observed TALL block for a PCA estimation of the factors and then obtains the loadings from a time-series regression that that uses all observed entries. They provide the inferential theory for this TALL-PROJECT estimator. Each of these estimators is designed for a specific observation pattern under which it performs particularly well, but might not generalize to other patterns. In contrast, we view our estimator as a simple all-purpose estimator that can reliably impute missing data and provide the correct confidence intervals for general missing patterns and factor structures, which makes it appealing for applied researchers in causal inference. In an extensive simulation study, we show that while our estimator has a similar performance as jin2020factor for data missing at random, and as bai2019matrix for missing with a block structure, our estimator can have a better performance for a staggered design or when the observation pattern depends on unit-specific features.
Our imputed values are point-wise consistent and have asymptotic normal distributions, which is relevant for the matrix completion literature that studies a similar problem. Both our paper and the matrix completion literature assume a low-rank structure in the panel data. In the matrix completion literature, the most popular method is to estimate the low-rank matrix from a convex optimization problem using a nuclear norm regularization mazumder2010spectral,negahban2011estimation,negahban2012restricted. The main results in the matrix completion literature are upper bounds for the mean-squared estimation error of the estimated matrix. However, point-wise consistency does not hold in general because the typically used nuclear norm regularization results in a bias in the estimated matrix. In their path-breaking work, chen2019inference propose de-biased estimators and provide an inferential theory under the assumption of i.i.d. sampling and i.i.d. noise. There is a trade-off in terms of the generality of the model and the required observations, where our work allows for the most general patterns in missing observations with a general approximate factor structure at the cost of observing entries at a higher rate than chen2019inference. Our paper contributes to the matrix completion literature by allowing general observation patterns and dependent error structures, which is particularly relevant for applications in social science.
Our paper allows for heterogeneous and time-dependent treatment effects of an intervention and more general intervention adoption patterns compared with the synthetic control methods in causal inference. Furthermore, our paper provides a flexible test for treatment effects. In comparative case studies, a key question is to estimate the counterfactual outcomes for treated units. A valid control unit is “close” to the treatment unit except for the treatment effect. Typically synthetic controls are weighted averages of untreated units where the weights depend on unit-specific features. A popular model assumption is that the potential outcome is linear in observed covariates and unobserved common factors. abadie2010synthetic,abadie2015comparative, doudchenko2016balancing, xu2017generalized, li2017estimation and li2017statistical propose to match each treated unit by weighted averages of all control units using the pretreatment observations. li2017estimation, li2017statistical and masini2018counterfactual show the inferential theory for the average treatment effect over time. These methods rely on the assumption that there is only one treated unit and the treatment effects are either constant or stationary. Another method is to regress the post-treatment outcomes for the control units on the pre-treatment outcomes and covariates and use the coefficients to predict the counterfactual outcome for the treated/control units. athey2018matrix propose to use matrix completion methods to impute the control panel data and allow for more general treatment adoption patterns: multiple treated units and staggered treatment adoption. However, they do not provide point-wise guarantees for the imputed values. In this paper, in addition to allowing for general treatment adoption patterns, we also provide the point-wise inferential theory for the imputed counterfactual outcomes. Furthermore, we can test for treatment effects even if they are heterogeneous and time-dependent. Our approach does not require a priori knowledge about which covariates describe if treated and control units are a good match. Instead, our latent loadings capture all unit-specific information in a data-driven way. The synthetic control, that we impute, is a weighted average of the untreated units, that takes all unit-specific information into account. In causal inference, we can either model the relationship between the covariates and the outcome, or model the probabilities of missingness to estimate causal effects. Doubly robust procedures, as discussed, for example, in kang2007 combine both by using a propensity weight in regressions to mitigate the selection bias. Our propensity weighted estimator builds on this intuition. Interestingly, we prove that using the estimated feasible propensity instead of the population weights does not affect the asymptotic distribution. This observation is aligned with the results for the classical inverse propensity weighted estimator in hirano2003efficient.
The rest of the paper is organized as follows. Section (ref) introduces the model and provides the simple all-purpose estimator for factors, loadings, and common components. Section (ref) states the necessary assumptions for the asymptotic distribution results that are presented in Section (ref). Sections (ref) and (ref) extend the results to the propensity weighted estimator. Section (ref) shows how to apply our model to test treatment effects. We discuss the feasible estimation in Section (ref) and how to relax the rate conditions in Section (ref). The extensive simulation in Section (ref) shows the good finite sample properties, the strong performance relative to other methods, and robustness results under misspecification. The Internet Appendix collects additional simulation results and all proofs.
Assume we partially observe a panel data set $Y$ with $T$ time periods and $N$ cross-sectional units. $Y \in \+R^{N \times T}$ has a factor structure with $r$ common factors. We denote by $F_t \in \+R^{r}$ the latent factors, $\Lambda_i \in \+R^{r }$ the factor loadings, $C_{it} = \Lambda_i^\top F_t$ the common component, and $e_{it}$ the idiosyncratic error:
\[Y_{it} = \Lambda_i^\top F_t + e_{it} \quad \text{for $i = 1, 2, \cdots, N$ and $t = 1, 2, \cdots, T$}\] or in vector notation, \[\underbrace{Y_t}_{N \times 1} = \underbrace{\Lambda}_{N \times r} \underbrace{F_t}_{r \times 1} + \underbrace{e_t}_{N \times 1} \qquad \text{for $ t = 1, 2, \cdots, T$}.\]
In an asymptotic setup where $N$ and $T$ are both large, we randomly observe some entries in $Y$. Let $W_{it} \in \{0, 1\}$ be a binary variable, where $W_{it}= 1$ indicates that the $(i,t)$-th entry is observed and $W_{it} = 0$ otherwise. In this paper, we will estimate the latent factors $F$ and loadings $\Lambda$ from the partially observed $Y$, impute the missing values, and provide the inferential theory for all estimators.
We allow for very general patterns in the missing observations. Figure (ref) shows three important examples widely seen in empirical applications. The first one is a randomly missing pattern, that is, whether an entry is observed or not does not depend on other entries or observable covariates. For example, the observational pattern of the Netflix challenge is usually modeled as entries missing at random. The second and third ones are the observation patterns for control panels in simultaneous and staggered treatment adoptions. Once a unit adopts the treatment, it stays treated afterwards, which will be modeled as missing values. These two patterns are widely assumed in the literature on causal inference in panel data.\footnote{See candes2009exact,zhou2008large for the Netflix challenge and athey2018matrix,athey2018design for missing patterns used in causal inference.}
$\mathcal{Q}_{ij} = \{t: W_{it} = 1 \text{ and } W_{jt} = 1 \}$ denotes the set of time periods $t$ when both units $i$ and $j$ are observed. $|\mathcal{Q}_{ij}|$ is the cardinality of the set $\mathcal{Q}_{ij}$. Assumption (ref) states the conditions on the observation pattern.
Assumption (ref) allows very general observation patterns that can vary over time and depend on unit-specific features. In particular, the observation pattern can depend on the factor loadings that capture cross-sectional information. For the purpose of identification, we assume that the observation pattern is independent of the factors. Note that the estimator of the common components is “symmetric” in $N$ and $T$, and therefore we could switch the roles of $N$ and $T$ in the above assumptions. In that case, the observation pattern would be independent of the loadings but can depend on the factors. The assumption that the observation pattern is independent of the errors is closely related to the unconfoundedness assumption in rosenbaum1983central. Assumption (ref) implicitly assumes that for any two units, the number of time periods when both are observed is proportional to $T$. This simplifies the presentation of our results and is sufficient for most empirically relevant cases, but we will also discuss how this assumption can be relaxed.
Our framework allows for the following important examples:
We provide an “all-purpose” estimator without the need to explicitly model the probability distribution $P(W_{it}=1)$. However, in some cases, we might have additional information about the missing pattern. In Section (ref), we provide a modification of our estimator that can take advantage of a model for $P(W_{it}=1)$. More specifically, we allow the cross-sectional observation pattern $P(W_{it}=1|S_i)$ to depend on observed cross-sectional features $S = [S_i] \in \+R^{N \times K}$. These covariates $S_i$ are assumed to be time-invariant. They can be discrete, for example, an indicator variable for gender in the evaluation of a drug treatment or continuous, for example, standardized past test scores in the evaluation of an educational policy change. The cross-sectional features $S_i$ can actually be the estimated latent loadings $\Lambda_i$ themselves. We discuss how the observation probability $P(W_{it}=1|S_i)=p_t(S_i)$ can be estimated with parametric or non-parametric estimators. While this modified estimator requires some changes to Assumption (ref), it provides the same level of generality for the missing pattern, as discussed in Section (ref).
There are two steps to estimate the latent factor model from the partially observed panel data: First, we need to estimate the covariance matrix of the data, and second we estimate the latent factors and loadings based on the eigenvectors of the estimated covariance matrix. The conventional latent factor estimator without missing values applies principal component analysis to the sample covariance matrix. A natural way to deal with the missing values is to set these entries to zero. However, the conventional PCA estimator will then be biased. Our estimator correctly re-weights the entries in the covariance matrix before applying PCA.
We first impute the missing entries by 0 and denote the imputed matrix as $\tilde Y$:\footnote{In matrix notation, we have $\tilde Y = Y \odot W$, where $\odot$ denotes the Hadamard product.} \[\tilde Y_{it} = Y_{it} W_{it}, \quad \text{for } i = 1, 2, \cdots, N \text{ and } t = 1, 2, \cdots, T. \] When some entries in $Y$ are missing, the conventional sample covariance estimator $\frac{1}{T} \tilde Y \tilde Y^\top$ is biased because the actual realizations of the missing values are not equal to zero. We propose a natural estimator of the covariance matrix, where for each entry we only use the time periods when both units are observed. This is equivalent to estimating the sample covariance matrix with $\tilde Y$, but reweighting the entries. Table (ref) is a simple example to illustrate the covariance matrix estimation if the entries are partly missing in the second half of the data. More generally, our sample covariance matrix estimator equals
When the data is fully observed, we can apply PCA to $\frac{1}{NT} Y Y^\top $ to estimate the loadings. Up to rescaling, the eigenvectors of the largest eigenvalues estimate the loadings. Then, we regress $Y$ on the estimated loadings to obtain an estimate of the factors.\footnote{Alternatively, we can apply PCA to $\frac{1}{NT} Y Y^\top $ to estimate the loadings and then regress $Y^{\top}$ on the estimated loadings to estimate the factors. bai2002determining and bai2003inferential show that the estimated factors and loadings from this approach are consistent and asymptotically normally distributed for a fully observed panel.} Similarly, for the partially observed data, we apply PCA to $\frac{1}{N} \tilde\Sigma$ to estimate the loadings.\footnote{We assume that the true number of factors is $r$ and has been consistently estimated as in bai2003inferential.} Under the standard identification assumption $\tilde \Lambda^\top \tilde \Lambda/N = I_r$, we estimate the loadings $\tilde \Lambda$ as $\sqrt{N}$ times the eigenvectors of the $r$ largest eigenvalues of the sample covariance matrix, that is
where $\tilde D$ is a diagonal matrix. Then, for every time period $t$, we regress the observed $Y_{t}$ on $ \tilde \Lambda$ to estimate the factors:
Interestingly, this very simple estimator automatically corrects for the impact of general observation patterns. If we have additional information that allows us to model the observation pattern as $P(W_{it}=1|S_i)$, we propose an alternative weighted regression:
This conditional estimator uses the weights $\frac{1}{P(W_{it} = 1|S_i)}$ in the cross-sectional regression. The estimator for $\tilde F^S_t$ is motivated by the inverse propensity score estimator, which is widely used in causal inference.\footnote{Since horvitz1952generalization, weighting observations by inverse probability has been frequently used to account for missing data in mean estimation david1983nonrandom,little1986survey,little1988missing, regression analysis robins1994estimation,robins1995semiparametric, and causal inference hirano2003efficient.} The rationale is that the re-weighted observations correspond to a model where the data is cross-sectionally missing at random. More specifically, after re-weighting the observed data, the loadings should follow the same distribution as in the complete panel without missing observations. This could be relevant if units, that are exposed to specific factors, are more likely to miss. In the special case for cross-sectional missing at random, i.e., $P(W_{it}=1)=p_t$, the two estimators coincide.\footnote{The cross-sectional reweighting is only affected by cross-sectional differences in the probability $p_{it}$ at time $t$. Therefore we model only the dependency on a set of cross-sectional covariates $S_i$, but allow this dependency to be time-varying. We could without loss of generality consider the dependency on cross-sectional and time-varying covariates $S_{it}$. However, this would not change the cross-sectional reweighting at a specific time $t$, but would come at the cost of a more complicated notation.} We will first study the simple all-purpose estimator $\tilde F_t$ and extend it to the propensity weighted estimator $\tilde F^S_t$ in Section (ref). We show that both estimators are consistent and asymptotically normal. In most cases, $\tilde F_t$ is more efficient than the propensity score estimator, but $\tilde F^S_t$ can have desirable robustness properties under misspecification. The last step is to estimate the common component $C_{it} = \Lambda_i^\top F_t$. We use the plug-in estimator, $\tilde C_{it} = \tilde \Lambda_i^\top \tilde F_t$ respectively $\tilde C^S_{it} = \tilde \Lambda_i^\top \tilde F^S_t$. If $Y_{it}$ is not observed, we impute the missing values with $\tilde C_{it}$ or $\tilde C^S_{it}$.
We illustrate in a simple example how missing observations change the conventional PCA estimator with fully observed data. Assume that we have only one factor, and the factor, loading and residual component are i.i.d. normally distributed with $F_t \overset{i.i.d.}{\sim} \mathcal{N}(0,\sigma_F^2)$, $\Lambda_i \overset{i.i.d.}{\sim} \mathcal{N}(0,1)$ and $ e_{it} \overset{i.i.d.}{\sim} \mathcal{N}(0,\sigma_e^2)$. We assume that the observations for units $1,...,N_0$ and for the times $T_0+1,...,T$ are missing according to the simultaneous adoption pattern of Table (ref). We separate the vector of factor realizations into its first $\mathbf{F}_1=
^{\top} $ and second part $\mathbf{F}_2=
^{\top} $ and similarly for the loadings $\mathbf{\Lambda}_{1}=
^{\top}$ and $\mathbf{\Lambda}_{2}=
^{\top}$. Note that in this simple example $\tilde F$ and $\tilde F^S$ coincide.
We start with the simplest case without error terms $e_t$ to illustrate the logic of reweighting entries. In this case the conventional covariance matrix equals
Obviously, the eigenvector of this matrix is a biased estimate of the loadings. In contrast, the eigenvector of the correctly weighted sample covariance matrix consistently estimates the loadings:
The same logic carries over to the estimator of the factors. Assume that we know the population loadings, which we use here instead of the estimated loadings in the regression to estimate the factors:
which is a biased estimator for the second time period. The regression in Equation (ref) corresponds to a weighted least square regression which provides the correct estimator:
The proper reweighting in the loading and factor estimation leads to an additional correction term in the asymptotic variance of the estimator. As an illustration of this additional challenge, we add the i.i.d. error term $e_{it}$ to our example. In our simplified setup our consistent estimator for the loadings $\tilde \Lambda$ has the following expansion for $i=1,...,N_0$:\footnote{The results are similar for $i > N_0$ with the expansion $\sqrt{T} \left( \tilde \Lambda_i - \Lambda_i \right) =\sqrt{T} \left( \frac{\tilde F^{\top} \tilde F}{T}\right)^{-1} \frac{1}{N} \Big[ \mathbf{\Lambda}_1^{\top}\mathbf{\Lambda}_1 \frac{1}{T_0} \sum_{t=1}^{T_0} F_t e_{it}$ $+ \mathbf{\Lambda}_2^{\top}\mathbf{\Lambda}_2 \frac{1}{T} \sum_{t=1}^{T} F_t e_{it} \Big]+ \sqrt{T} \left( \frac{\tilde F^{\top} \tilde F}{T} \right)^{-1} \frac{1}{N} \mathbf{\Lambda}_1^{\top}\mathbf{\Lambda}_1 \left( \frac{\mathbf{F}_1^{\top} \mathbf{F}_1}{T_0} - \frac{F^{\top}F}{T} \right) \Lambda_i + o_P(1)$ and asymptotic distribution $\sqrt{T} \left( \tilde \Lambda_i - \Lambda_i \right) \overset{d}{\rightarrow} \mathcal{N} \left( 0, \text{plim}\left( \Big( \frac{T-T_0}{T_0} \frac{N_0^2}{N^2} + 1 \Big) \frac{\sigma_e^2}{\sigma_F^2} + 2 \frac{N_0^2}{N^2} \frac{T-T_0}{T_0} \Lambda_i^2 \right) \right)$. }
which results in the asymptotic normal distribution
The second term in the asymptotic expansion is due to averaging over different number of units for different elements of the loadings. This additional variance correction term vanishes for $T_0 \rightarrow T$. Similar terms appear in the distribution of the estimators of the factors and common components. We show under general conditions how these correction terms arise in the asymptotic distribution and how to take them into account for the inferential theory.
We assume an approximate factor structure at the same level of generality as in bai2003inferential. The factors and loadings have non-trivial time-series and cross-sectional dependency. We allow the errors to be weakly correlated in the time-series and cross-sectional dimension. The asymptotic distributions are based on general martingale central limit theorems. The general Assumptions (ref) and (ref) are collected in the Appendix. In the main text, we present a simplified factor model with the stronger Assumptions (ref) and (ref), which substantially simplifies the notation but conveys the main conceptual insights of the general model. It allows us to highlight the effect of missing observations.
The consistency results are based on Assumption (ref) that assumes that all observations are i.i.d. The key elements are that the factors and loadings are systematic in the sense that they lead to exploding eigenvalues, while the error terms are non-systematic with bounded eigenvalues in the covariance matrix of $Y$. These are standard factor model assumptions. The asymptotic distribution results require additional restrictions on the missing patterns, as stated in Assumption (ref).
Assumption (ref) has two key elements. First, the full rank assumption of $\Sigma_{\Lambda, t}$ captures that the factor loadings are systematic for the observed entries. Second, the number of observed units at every time period $t$ is proportional to $N$ and different units share a number of observed entries that is proportional to $T$. The impact of the missing pattern on the asymptotic variances of the estimators is captured by the three key parameters $\omega, \omega_j$ and $\omega_{jj}$. Note that by construction these constants satisfy $\omega_{jj}, \omega_j, \omega \geq 1$. If the observations are missing at random with probability $p$, then $\omega_{jj} = \frac{1}{p}$, $\omega_j = 1$ and $\omega= 1$.
As stated in Proposition (ref) in the Appendix, the simplified model is just a special case of the general approximate factor model specified by Assumptions (ref) and (ref). The simplified Assumption (ref) implies the general Assumption (ref), while Assumption (ref) combined with the other simplified assumptions implies the general Assumption (ref).
We assume that the number of factors $r$ is consistently estimated. For example the criteria developed in bai2002determining can be extended to our case of missing values based on the various bounds and expansions that we derive in this paper. A promising alternative would be to extend the cross-validation estimator of jin2020factor or an eigenvalue ratio argument as in ahn2013 to general missing patterns. Given a consistent estimator for the number of factors, we can treat $r$ as known.
We first show the consistency of our estimators. Our analysis starts with plugging $\tilde Y = (\Lambda^\top F + e) \odot W$ into Equation ((ref)) which yields the following decomposition:
Similar to bai2002determining this decomposition relates the estimated loadings to the population loadings, $\tilde \Lambda_j = H_j \Lambda_j + (a) + (b) + (c)$, up to a rotation matrix $H_j = \frac{1}{NT} \tilde{D}^{-1} \sum_{i=1}^N \tilde{\Lambda}_i \Lambda_i^\top F^\top \mathrm{diag}(W_i \odot W_j) F/q_{ij}$. The key difference to factor analysis with fully observed data is that this rotation matrix can be different for different units $j$. However, the estimation of the factors is based on a projection on the loading space and hence implicitly requires the same rotation matrix for all loadings.
We consider for all units a unified rotation matrix defined as $H = \frac{1}{NT} \tilde{D}^{-1} \tilde{\Lambda}^\top \Lambda F^\top F$ which is essentially the same conventional rotation matrix as in bai2002determining. This yields the decomposition \[\tilde{\Lambda}_j - H \Lambda_j = \tilde{\Lambda}_j - H_j \Lambda_j + (H_j - H) \Lambda_j = (a) + (b) + (c) + (H_j - H) \Lambda_j. \] We show that the cross-section averages of the square of $(a)$, $(b)$ and $(c)$ converge to 0 at the rate $O_P \left( \min\left( \frac{1}{N}, \frac{1}{T} \right)\right)$. The key difference compared with the fully observed factor analysis is the last term. If $\frac{1}{T} F^\top F \xrightarrow{P} \Sigma_F$ and $\frac{1}{|\mathcal{Q}_{ij}|} \sum_{t \in \mathcal{Q}_{ij}} F_t F_t^\top \xrightarrow{P} \Sigma_F$, we can show that $H_j - H = O_P \left( \min \left( \frac{1}{\sqrt{N}}, \frac{1}{\sqrt{T}} \right) \right)$. This rate is sufficiently fast to obtain consistency, but will contribute to the asymptotic normal distribution. Note that the correction term $H_j - H$ is a fundamental problem for any estimator that makes use of all observations.\footnote{The estimator in bai2019matrix can avoid this term by neglecting partially observed entries, which means that in general, they are using less information. The estimator of bai2019matrix is optimized for the block structure of a simultaneous adoption pattern. It runs two PCA estimates for the block with full cross-sectional observations and the block with full time-series observations. Hence, they can infer the “local” rotation matrices for each block and rotate the estimates to avoid the correction term $H_j - H$. cahan2021factor leverage the block structure, and in the first step run PCA to estimate factors on the block with full-series observations, so that the correction term $H_j - H$ is avoided. If we seek to use full observations in the first step, which is what we propose in this paper, then the correction term $H_j - H$ cannot be avoided. }
The next theorem shows the consistency of the estimated loadings.
Theorem (ref) states that the complete loading matrix can be consistently estimated up to an appropriate rotation as $N,T \rightarrow \infty$ even if we only observe an incomplete panel matrix. The convergence rate is the same rate as for the fully observed panel in bai2002determining. Theorem (ref) is based on the assumption that the observed entries are representative of the missing entries and hence provide a consistent estimation. Theorem (ref) is a critical intermediate step to show the asymptotic normality of the estimated factor model in the next section.
The factors, loadings, and common components are asymptotically normally distributed. Indeed, Theorem (ref) states that the asymptotic distributions have two parts: First, we recover the asymptotic variance that is identical to the conventional PCA in bai2003inferential under the same rate conditions. These are the expression when we set the additional correction terms $\Gamma^{\textnormal{miss}}_{\Lambda,j}$ and $\Gamma^{\textnormal{miss}}_{F,t} $ to zero. However, in the presence of missing values, these correction terms are necessary to capture the additional uncertainty. Theorem (ref) also includes the asymptotic expansions that lead to the normal distributions. As stated in the previous section, the difference between the unit-specific rotation $H_j$ and the “global” rotation matrix $H$ contributes to the distribution and leads to the variance correction terms $\Gamma^{\textnormal{miss}}_{\Lambda,j}$ and $\Gamma^{\textnormal{miss}}_{F,t}$. As expected, this variance correction is increasing in the number of missing observations. We want to emphasize again that this type of variance correction is a conceptual issue that cannot be avoided when making use of all observed entries.\footnote{In the asymptotic distribution, we apply the rotation matrices to the estimated loadings and factors instead of their population values as in bai2003inferential. Obviously, these two representations are equivalent and can be easily transformed into each other. Our choice of representation was made for exposition purposes only.}
Importantly, the estimator for the factors has a different convergence rate compared to the conventional estimator on fully observed data. The asymptotic distribution of the factors is determined by two terms with different convergence rates, $\frac{\delta_{NT}}{N} \Gamma^{\textnormal{obs}}_{F,t} +\frac{\delta_{NT}}{T} \Gamma^{\textnormal{miss}}_{F,t}$. With a fully observed panel $ \Gamma^{\textnormal{miss}}_{F,t}$ would disappear, and the factors would converge at a rate of $\sqrt{N}$. However, with observations that are not missing at random, the difference between $H_j$ and $H$, that appears in the loading expansion $ H^{-1} \tilde \Lambda_j - \Lambda_j $ and has a convergence rate of $\sqrt{T}$, also contributes to the asymptotic distribution of estimated factors, which results in the overall rate $\sqrt{\delta_{NT}}$.
The asymptotic distribution of common components depends on the estimation error of the estimated loadings and factors. In the asymptotic distribution of the estimated loadings and factors, the conventional part with asymptotic variances $\Gamma^{\textnormal{obs}}_{\Lambda,j}$ and $\Gamma^{\textnormal{obs}}_{F,t}$ is asymptotically independent as argued in bai2003inferential. However, the second part with the asymptotic variances $\Gamma^{\textnormal{miss}}_{\Lambda,j}$ and $\Gamma^{\textnormal{miss}}_{F,t}$ that captures the difference between $H_j$ and $H$ is in general correlated, and hence their covariance $\Gamma^{\textnormal{miss, cov}}_{\Lambda, F, j, t}$ contributes to the asymptotic variance of common components as stated in Equation (ref).
The correction terms in the asymptotic variances are determined by the functions $h_j(.)$, $g_t(.)$ and $g^{\mathrm{cov}}_{j,t}(.,.)$. These are quadratic functions in the elements of $\Lambda_i$ and $F_t$, which depend on the moments of the factor model. As the loadings and factors are random, it implies that the asymptotic variances themselves are random. This complicates the analysis, but our assumptions ensure that the normalized estimates converge to a standard normal distribution.
The distribution results of Theorem (ref) simplify under Assumptions (ref) and (ref), and we can provide explicit expressions for the asymptotic variances. If we assume in addition that the proportions of observed time-series ($q_{ij}$ and $q_{ij,kl}$) are independent of the second moment of the loadings $\Lambda_i \Lambda_i^\top$, we can further separate the effect of missing patterns from the properties of the factor model.
The simplified model provides a clear interpretation of the effect of missing data. Importantly, the parameters $\omega, \omega_{j}$ and $\omega_{jj}$, that depend only on the missing pattern, but not on the factor model, determine the weights of correction terms. The asymptotic covariance of the loadings is a weighted combination of the variance of an OLS regression of the population factors $F$ on $Y_j$ and the correction term. The weight $\omega_{jj} \geq 1$ depends on the number of the observed entries and the similarities in observation patterns for different units. Without missing data, it equals $\omega_{jj}=1$ and the correction term disappears. If the data is observed uniformly at random with probability $p$, the weight equals $\omega_{jj}=1/p$ which is increasing in the proportion of missing observations.
Similarly, the asymptotic variance of the factors has two components: the variance of an OLS regression of the population loadings on $Y_t$ using only observed entries, and the correction term. The weight $\omega \geq 1$ increases the scale of the correction term. When all entries are observed, or all entries are observed cross-sectionally at random (with either the same or different probabilities), then $\omega = 1$, the correction term vanishes, and the asymptotic variance only depends on $\Sigma^{\textnormal{obs}}_{F,t}$. If the missing pattern does not depend on the loadings, then $\Sigma_{\Lambda,t}=p_t \Sigma_{\Lambda}$ and $\Sigma^{\textnormal{obs}}_{F,t}$ simplifies to $\frac{1}{p_t} \Sigma_{\Lambda}^{-1} \sigma_e^2$ which is the variance of an OLS regression of the population loadings on $Y_t$ scaled by the inverse proportion of observed entries at time $t$.
The distribution of the common component depends on all three parameters $\omega, \omega_{j}$ and $\omega_{jj}$. If all entries are observed at random, then $\omega_j = 1$ and the contribution of the loading and factor distribution to the common component are separated similar to the conventional PCA setup in bai2003inferential. In this case, only the two terms $\omega_{jj} F_t^\top \Sigma^{\textnormal{obs}}_\Lambda F_t$ and $\Lambda_j^\top \Sigma^{\textnormal{obs}}_{F,t} \Lambda_j$ remain in the asymptotic variance.
We provide the assumptions and general distribution theory for the propensity weighted estimator for the factors $\tilde{F}^S_t$ defined in Equation (ref). This conditional estimator uses the weights $\frac{1}{P(W_{it} = 1|S_i)}$ in the cross-sectional regression to obtain the factors. We allow the observation probability to depend on observed cross-sectional features $S = [S_i] \in \+R^{N \times K}$ that explain why certain units are more likely to be observed than other units. This conditional setup requires some modifications of the previous assumptions. In addition to Assumption (ref) we require the following assumption:
We assume $S$ contains all the information in $\Lambda$ that is predictive for the observation pattern. In other words, $W$ is independent of $\Lambda$ conditional on $S$, as stated in Assumption (ref).1. This is closely related to the {unconfoundedness} assumption in causal inference. It also assumes that the conditional probability $P(W_{it}=1|S_i) $ is bounded away from 0, which implies that the number of observed cross-sectional and time-series entries is proportional to $N$ and $T$, respectively. This corresponds to the {overlap assumption} in causal inference.\footnote{See rosenbaum1983central for the connection to unconfoundedness and the overlap assumption. We assume $P(W_{it} = 1|S)$ is bounded away from 0, such that $\frac{1}{P(W_{it}=1|S)}$ does not diverge, which is equivalent to the overlap assumption in causal inference.} Note that it is straightforward to include the covariates of “neighbor units” in $S_i$ to allow for network effects.
We replace Assumptions (ref) and (ref) by their conditional counterpart Assumptions (ref) and (ref) which have a similar level of generality. These are required for the asymptotic normality of $\tilde{F}^S_t$ and $\tilde C^S_{it}$. As before, we collect the Assumptions (ref) and (ref) for a general approximate factor model in the Appendix and present the assumptions for a simplified factor model in the main text, which are sufficient to convey all conceptual insights.
\setcounter{assumpC}{2}
Under Assumption (ref), $\frac{1}{N} \sum_{i = 1}^N \frac{W_{it}}{P(W_{it} = 1|S_i)} \tilde \Lambda_i \tilde \Lambda_i^\top $ converges in probability to an identity matrix which is the same limit as the loading estimates in conventional PCA without missing data. The assumption that $S$ is independent of $F$ and $e$ is conceptually similar to the assumption that $\Lambda$ is independent of $F$ and $e$, where the latter is standard in the literature on large dimensional factor modeling. The additional moment conditions in Assumption (ref) are required for the asymptotic distribution, where $\Sigma_{\Lambda,S,t}$ appears in the asymptotic covariances of $\tilde{F}^S_t$ and $\tilde C^S_{it}$.
Proposition (ref) in the Appendix shows that the simplified model is just a special case of the general approximate factor model specified by Assumptions (ref) and (ref). The simplified Assumption (ref) combined with Assumptions (ref), (ref) and (ref) imply the general conditional Assumption (ref), while Assumption (ref) combined with the other simplified Assumptions (ref), (ref), (ref), (ref).2 and (ref) imply the general Assumption (ref).
The propensity weighted estimator only differs in the distribution of the factors and common components. Both $\tilde F^S_t$ and $\tilde C^S_{it}$ follow a normal distribution, but in most cases have a larger asymptotic variance than the estimators $\tilde F_t$ and $\tilde C_{it}$. The loadings are not affected by the propensity score weighting.
The distribution results have the same general structure as in Theorem (ref). However, there are two key differences. First, the outer matrices in the variance of $\tilde F^S_t$ are $\Sigma_{\Lambda}^{-1}$ while they depend on the observational pattern in $\Sigma_{\Lambda,t}^{-1}$ in Equation (ref). Second, the middle terms $\Gamma^{\textnormal{obs},S}_{F,t}$ and $\Gamma^{\textnormal{miss},S}_{F,t}$ may depend on $P(W_{it}=1|S_i)$. The same structure carries over to the common component. In the case of generalized least squares regressions, it is straightforward to compare the asymptotic covariances for different weights and to determine an efficient estimator. With missing observations, the problem becomes more challenging as the asymptotic covariances depend on two matrices for the factor estimates and three terms for the common components. For the general models in Theorems (ref) and (ref) we cannot state which estimator is more efficient without imposing additional structure. However, for the simplified model, we can rank the efficiency of the two estimators.
An interesting observation is that $\Sigma^{\textnormal{miss},S}_{F,t}$ and $\Sigma^{\textnormal{miss, $S$, cov}}_{\Lambda,F,j,t}$ depend neither on the observation pattern nor on $S$. This is because $\frac{1}{P(W_{it}=1|S_i)}$ removes the asymptotic dependency between $W_{it}$ and $\Lambda_i$. Hence, this part of the asymptotic distribution has a complete separation between the missing observation pattern captured by the weights $\omega, \omega_{j}$ and $\omega_{jj}$ and distribution terms that depend only on the factor model. However, $\Sigma^{\textnormal{obs},S}_{F,t}$ depends on $P(W_{it}=1|S_i)$ as this component comes from a probability weighted least square regression of the population loadings on the observed entries in $Y$, which is different from the corresponding OLS regression in Corollary (ref).1.
The key observation is that $\tilde F_t^S$ and as a consequence also $\tilde C_{it}^S$ seem to be in many cases less efficient than $\tilde F_t$ and $\tilde C_{it}$, which means that the asymptotic variances of the all-purpose estimator are less than or equal to those of the propensity weighted estimator. In the case of only one factor it holds that $ \Sigma^{\textnormal{miss},S}_{F,t} = \Sigma^{\textnormal{miss}}_{F,t} $. Not surprisingly, it holds that $\Sigma^{\textnormal{obs},S}_{F,t} \geq \Sigma^{\textnormal{obs}}_{F,t}$ as in the case of i.i.d. errors, an OLS regression is the most efficient linear estimator. This result can also be derived from $\Sigma_{\Lambda,S,t}/\Sigma_{\Lambda}^2 \geq 1/\Sigma_{\Lambda,t}$, which follows from the Cauchy-Schwartz inequality. In the case of multiple factors, we take advantage of the concavity of the average weighted by $1/P(W_{it}=1|S_i)$ to prove the efficiency relationship. In simulations we confirm that when the loadings depend on $S$, it is possible that $ \Sigma^{\textnormal{miss},S}_{F,t} < \Sigma^{\textnormal{miss}}_{F,t} $, which can result in minor efficiency gains for $\tilde F_t^S$. For a general residual covariance matrix, the efficiency results are more complex. pelgerxiongsparse2020 show that for fully observed data under certain assumptions, the optimal weight in the factor regression is the inverse residual covariance matrix or equivalently PCA, applied to a covariance matrix re-weighted by the square-root of the inverse residual covariance matrix, is the most efficient estimator. Hence, if the propensity weight is close to the inverse residual covariance matrix, it lowers the first term $\Sigma^{\textnormal{obs},S}_{F,t}$. However, the effect on the second term is more complex, and hence there are in general cases where $\tilde F_t^S$ can be more efficient than $\tilde F_t$. In simulations, we show that for a correctly specified model, the estimates of $\tilde F_t$ and $\tilde F_t^S$ are close, but $\tilde F_t$ is generally more precise. However, the “doubly-robust” estimator $\tilde F_t^S$ seems to be less affected by various forms of misspecification, e.g., omitted factors, weak factors, or a nonlinear factor model. Hence, $\tilde F_t^S$ might be appealing because it is more robust but not based on efficiency arguments.
The propensity weighted regressions can be robust to the selection bias from omitting factors. In the causal inference literature regressions weighted by propensity scores have been used in the estimation of causal effects to reduce the bias that arises from omitting regressors or misspecifying the outcome model. However, as the propensity weighted regressions have a larger variance, regressions without the propensity weights seem to be preferred for correctly specified models.\footnote{robins1994estimation,robins1995semiparametric among others discuss the reduction of bias from omitting regressors. robins2000inference,kang2007,robins2007comment show the large variance of propensity weighted regressions, and freedman2008weighting suggests unweighted regressions for correctly specified models} In this section, we illustrate that this logic carries over to our latent factor model setup.
Our setup differs from classical causal inference as we estimate the covariates as latent factors from the data. However, we can have a situation similar to omitted variables if we estimate too few latent factors, the factors are weak, or the population model is nonlinear. In the previous section we have shown that the propensity weighted estimator is in general less efficient than our all-purpose estimator when we use the correct number of factors $r$.\footnote{Given our distribution theory, it is relatively straightforward to extend the consistent estimator of bai2002determining for the number of factors to the more general case with missing data. However, most existing estimators for the number of latent factors explicitly or implicitly depend on choice parameters, which implies that in practice it is possible to use too few factors pelger2019large,lettau2018factors.} However, when we use too few latent factors, our weighted estimator can have a smaller selection bias, and hence be preferable. We illustrate the general logic of this result with an example and confirm it with extensive simulations in Section (ref).
We assume a two-factor model $Y_{it} = \Lambda_{i1} F_{t1} +\Lambda_{i2} F_{t2} + e_{it}$, where
The key assumption is that the observation pattern depends on the loadings. In order to have a transparent example, we assume the following pattern. The first $T_0$ time periods are fully observed. After time $T_0$, whether a unit is observed or not depends on an indicator variable $S_i $, defined as $S_i = \boldsymbol{1}_{\Lambda_{i1} > c_1, \Lambda_{i2} > c_2}$ for some $c_1 > 0$ and $c_2 > 0$. Suppose $P(W_{it} = 1|S_i = 1) = p$ and $P(W_{it} = 1|S_i = 0) = 1 - p$ for some $p > 0.5$. In other words, units with large loadings are more likely to be observed. We would get similar results if large loadings are more likely to be missing. Without loss of generality, we can shuffle the units and obtain the observation pattern in Table (ref). In the following, we estimate the factor model from the shuffled outcome matrix, where the outcomes for the first $N_0$ units after time $T_0$ are missing, and the outcomes for the last $N-N_0$ units are fully observed.
Assume that we omit one factor and estimate only a one-factor model with both, the simple and propensity weighted, estimators. In this case, our factor model is misspecified. Without loss of generality, we can set $\mu_\Lambda^2 + \sigma_{\Lambda}^2 = 1$. The estimated loading vector is consistent and the same for both estimators: \[\tilde{\*\Lambda}_1 = \*\Lambda_{1} + o_p(1).\]
We compare for both approaches the estimates of the first factor and the common components from time $T_0+1$ to $T$ (the two approaches coincide from time 1 to $T_0$ as all units are fully observed). For the simple regression and for $T_0 < t < T$, we can show that
where $\gamma = \lim_{N_0, N \rightarrow \infty} \frac{\sum_{i = N_0+1}^N \Lambda_{i1} \Lambda_{i2} }{\sum_{i = N_0+1}^N \Lambda_{i1}^2 }$. The key element is that $\gamma \neq 0$ because of the dependence of the missing pattern on $\Lambda$. In our example, both $\Lambda_{i1} $ and $ \Lambda_{i2}$ tend to have large values on the observed units, and hence are not asymptotically orthogonal on the subset of observed data. A non-zero $\gamma$ creates a selection bias similar to the conventional omitted variable bias.\footnote{simon1954spurious refers to this type of selection bias as a spurious correlation.}
Since the common component is $C_{it} = \Lambda_{i1} F_{t1} +\Lambda_{i2} F_{t2}$, the estimated common components have an error of $\tilde{C}_{it} - C_{it} =( \gamma \Lambda_{i1} +\Lambda_{i2} ) F_{t2} + o_p(1)$. Hence, the mean squared error from $T_0+1$ to $T$ equals
as $\frac{1}{N} \sum_{i=1}^N \Lambda_{i1} \Lambda_{i2} \rightarrow 0$. In contrast, the propensity weighted regression for $T_0 < t < T$ equals
Thus, the propensity weighted regression corrects for the selection bias and is the same as if $\gamma$ would have been zero. The estimated common components have an error of $\tilde{C}^S_{it} - C_{it} = \Lambda_{i2} F_{t2} + o_p(1)$, which implies that the mean squared error from time $T_0+1$ to $T$ is
Importantly, $\tilde{C}^S_{it} $ has a smaller mean squared error than $\tilde{C}_{it} $. In summary, when the estimated factor model has omitted factor(s) and is misspecified, the propensity weighted estimator could adjust for the selection bias and reduce the estimation error.
We provide feasible estimators for $P(W_{it}=1|S_i)$ which we need in Equation (ref) to estimate the factors, and we show that the asymptotic distribution of factors is not affected by using the estimated weights instead of their population counterpart. While in (stratified) randomized experiments, researchers decide and therefore know the treatment assignment probability given covariates, $P(W_{it}=1|S_i)$, the probability distribution of the missing pattern in observational studies generally needs to be estimated, which can affect the distribution theory for the latent factor model. Here we provide conditions under which the previously derived results continue to hold with a feasible estimator. To simplify notation denote by $p_{it}=P(W_{it}=1|S_i)$ the propensity score and its estimate by $\hat p_{it}=\wh P(W_{it}=1|S_i)$. The feasible estimator for the factors $\hat F_t^S$ replaces $p_{it}$ by $\hat p_{it}$ in Equation (ref), which yields the following decomposition:
Under weak assumptions on $\hat p_{it}$, that are satisfied for feasible estimators of the most empirically relevant observation patterns, the additional term $\hat F^S_t-\tilde F_t^S$ can be neglected in the asymptotic distribution.
We discuss feasible estimators for the most important cases of missing patterns which are summarized in Table (ref). Obviously, we only need to consider the case where $p_{it}$ varies for different cross-sectional units as otherwise the estimator simplifies to our estimator in Equation (ref). For simplicity these examples assume that $S_i$ are $i.i.d.$ and sub-Gaussian but can be generalized to weak dependency patterns. The simplest case is missing at random only in the time-series dimension, that is $P(W_{it}=1|S_i) = p(S_i)$ for some parametric or non-parametric function $p(.)$. A relevant example is the estimation of $p(S_i)$ with a logit model on the full panel $W$ which has the convergence rate $\hat p(S_i)=p(S_i)+ O_P \left( \frac{1}{\sqrt{NT}} \right)$ and a uniform bound of order $\frac{\log(NT)}{\sqrt{NT}}$. Hence, Theorem (ref).2(b) applies. If $p(S_i)$ is estimated non-parametrically with a kernel with bandwidth $h$, the convergence rate is typically $\sqrt{NT h}$ with a uniform bound of order $\frac{\log(NTh)}{\sqrt{NTh}}$, which does not change the distribution results if $Th$ is sufficiently large. In the more complex model $P(W_{it}=1|S_i)=p_t(S_i)$ the observations probability depends on the cross-section and time-series information. A relevant example for a parametric model is a logit model estimated on $W_t$ for each $t$ separately with a convergence rate of $\sqrt{N}$. Under weak assumptions on $S_i$, the uniform convergence bound in Theorem (ref).2(b) holds.
An important special case are discrete values for $S$, that is, the covariates $S$ take only finitely many values. An example for a binary variable $S$ would be gender, when male or female individuals have different probabilities to be treated. If the probabilities for the different discrete outcomes of $S$ are bounded away from zero, then the estimator $P(W_{it}=1|S_i=s)=p_t(s)$ simplifies to $p_t$, but just averaged over the cross-section units for which $S_i=s$. In more detail, consider the estimator $\hat p_t(s) = \frac{|\mathcal{O}_{s,t}|}{N_s}$ where $N_{s} = \sum_{i = 1}^N \mathbf{1}(S = s)$ and $\mathcal{O}_{s,t} = \{i: W_{it} = 1 \text{ and } S = s\}$. Then, $\sqrt{N_s} \left( \hat p_t(s) - p_t(s) \right) \xrightarrow{d} \mathcal{N} \left (0, p_t(s) (1- p_t(s)) \right)$. If $N_s$ is sufficiently large, for example proportional to $N$, then the feasible estimator does not change the distribution results in Theorem (ref). These estimators directly carry over to staggered treatment adoption. The staggered design can also be modeled with a parametric hazard model $P(W_{it}=1|S_i)=p(t,S_i)$, which under appropriate assumptions converges at the rate $\sqrt{N}$ as well. In summary, for all these cases the distribution results are not affected by using a feasible estimator for the propensity score.
As previously mentioned, we allow $S_i=\Lambda_i$. This is appealing as $\Lambda$ is by construction capturing the unit-specific features and hence should account for the differences in cross-sectional observation patterns. As the estimator $\hat \Lambda$ does not depend on the probability weights, it can be used in the estimation of $P(W_{it}|\Lambda_i)$. Theorem (ref) states that the estimation error of $\hat \Lambda_i$ is of the order $O_P \left( \frac{1}{\sqrt{N}} \right)$. While the consistency results for the factors and common components continue to hold, we need some additional weak assumptions on the tail behavior of the loadings and error terms to satisfy the uniform condition in Theorem (ref).2(b).
The key application of our the asymptotic distribution theory is to test causal effects. The fundamental problem in causal inference is that we observe an outcome either for the control or the treated data, but not for both at the same time. The unknown counterfactual of what the treated observations could have been without treatment can be naturally modeled as a data imputation problem. In this section, we consider the case where once a unit adopts the treatment, it stays treated afterward, for example, the simultaneous and staggered treatment designs illustrated in Figure (ref). Given the general missing patterns that we allow for, the generalization to more complex adoption patterns is straightforward. We denote by $T_{0,i}$ and $T_{1,i}$ the number of control and treated time periods for unit $i$ where their sum adds up to $T_{0,i}+T_{1,i}=T$. The superscripts $(0)$ and $(1)$ denote the observations for control and treated observations.
The individual treatment effect measures the difference between the treated and control outcomes:
where by construction for a specific time $t$ and unit $i$ we only observe either $Y_{it}^{(1)}$ or $Y_{it}^{(0)}$, but not both. Average treatment effects can be estimated by an average over time or the cross-section of the individual treatment effects. We assume that the data has a factor structure which results in a model of the form
where $D_{it}=1$ is a treatment indicator. Note that this model is very general and captures many relevant models as special cases. The factor structure includes interactive fixed effects as in Bai (2009). Simple time- and cross-sectional fixed effects are a special case for constant loadings respectively factors. The factors can be either observed covariates or latent factors. One of the main challenges in estimating a treatment effect is to control for all relevant covariates. Failure in doing so results in an omitted variable bias in the treatment effect estimation as discussed among others in gobillon2016regional. The strength of our latent factor model is that we can avoid this problem by automatically including all relevant covariates in a data driven way. Note that our latent factor model can also account for some uncertainty in the functional form of the dependency on the factors. For example, if $Y_{it}$ is a polynomial function of a factor, this could be captured by including additional latent factors as described for example in pelger2018state. A generalization of Equation (ref) adds additional observed covariates $X_{it} \in \mathbbm R^{K}$ to $Y_{it}$ which yields $Y_{it} = \tau_{it} D_{it} + \Lambda_i^{\top} F_t + X_{it}^{\top} b + e_{it}$. If these observed covariates follow a factor structure $X_{it} = {\Lambda^{X}_i}^{\top} F^X_t + e^X_{it}$, it puts us back into the framework of Equation (ref). Otherwise it is straightforward to include general observable covariates $X_{it}$ by studying the residual $Y_{it} - X_{it} \hat b$, where $\hat b$ is estimated by a regression on the control group.\footnote{Using the residuals $Y_{it} - X_{it} \hat b$ for the factor analysis and treatment effect analysis with our method generally adds another covariance term to the asymptotic covariance matrix. This term comes from the regression to obtain $\hat b$ and is straightforward to include. Here we focus on the conceptually more challenging problem of dealing with the unobserved factors.}
We only observe $Y_{it}^{(1)}$ for the treated group and could obtain the counterfactual outcome $Y_{it}^{(0)}$ from the imputed value $\hat Y_{it}^{(0)} =\hat C_{it}^{(0)}$, where $\hat C_{it}^{(0)}$ is the common component estimated only from the untreated control data. This is the same setup as in bai2019matrix. Given our asymptotic distribution theory for the common component, we can provide the asymptotic distribution of the individual and average treatment effects analogously to bai2019matrix. A shortcoming of estimating the individual treatment effect by $Y_{it}^{(1)} - \hat{Y}_{it}^{(0)}$ is that the observed treated observations $Y_{it}^{(1)}$ contain an idiosyncratic error $e_{it}$. Hence, it is not possible to test for individual treatment effects without imposing very strong additional assumptions on the error. For sufficiently large $T-T_0$, this error component can be averaged out in the average treatment effect.
We impose slightly stronger assumptions on the structure of the treatment effect which will allow us to derive substantially stronger results. Assume that the treatment effect has also a factor structure, that is $\tau_{it}=\left( {\Lambda_i^{\tau}} \right)^{\top} F_t^{\tau}$. In this case we can represent the problem as
where the factor structure subsumes the treatment effect. Hence, the individual treatment effect is equivalent to the difference in the common components between the treated and control:
Fundamentally, we are testing if the treatment changes the underlying factor structure. Hence, we can test if the treatment changes interactive fixed effects. This is a very general setup that allows for time and cross-sectional heterogeneity in the treatment effect, while the treatment itself can depend on the latent cross-sectional covariates modeled by the loadings.
In the following we consider three different treatment effects:
Here, $C^{(0)}_{i,(T_{0,i}+1):T} =
^\top \in \+R^{T_{1,i}}$ denotes the observations for $t >T_{0,i}$. The weighted average treatment effect $\tau_{\beta,i}$ generalizes the average treatment effect $\tau_i$, which is a special case for $Z = \vec{1}$. Both tests for the individual treatment effect and the weighted average treatment effect cannot be obtained with conventional estimators, but are important to answer economic questions. For example, in our companion paper pelgerxiongpublication2020, we test if pricing anomalies of investment strategies as measured by their pricing errors persist after these strategies have been published in academic journals. In this problem the treatment is the publication of an investment strategy in a journal and the treatment effect is measured by a change in regression coefficients. More specifically, the pricing error corresponds to the intercept in a regression of the excess returns of the strategies on a set of benchmark risk factors. A simple average treatment effect would not be sufficient to study this question.
For each of the three treatment effects we derive the asymptotic distribution under the null-hypothesis of no effect, which allows us to run one-sided or two-sided hypothesis tests. For example, the two-sided hypothesis test for the weighted average treatment effect takes the form
This is the hypothesis we test in our simulation and the empirical companion paper. The problem formulated in Equation (ref) can be solved by applying our latent factor model estimation twice: First, we estimate $C_{it}^{(1)}$ from the treated data with the control observations as missing values. Second, we estimate $C_{it}^{(0)} $ from the control data, while the treated observations are viewed as missing. The inferential theory follows readily from Theorems (ref) and (ref). The asymptotic variance for the individual treatment effect $\tau_{it}$ is the sum of the asymptotic variances of $\hat C_{it}^{(0)} $ and $\hat C_{it}^{(1)} $ and a covariance term based on the correction terms for the control and treated. While the calculations are tedious, they are a direct consequence of the distribution results that we have derived. The average treatment effects follow then from the results of the individual treatment effects. In this section, we want to focus on a special case, which we consider the most relevant from a practical perspective.
In most causal inference applications, such as the empirical study in our companion paper and abadie2010synthetic,abadie2015comparative, the majority of observations are control observations. Hence, it might be infeasible to estimate a latent factor model only from the treated data as required in Equation (ref). For example in the simultaneous treatment case in Table (ref), we can estimate a latent factor for the control, but not for the treated. Hence, we impose the additional assumption that the control and treated panel share the same underlying factors, while the loadings can be different, that is,
This implies that the treatment can only affect the loadings. This is still a very general setup as the loadings and factors are latent. For example, a model based on Equation (ref) with one factor that changes on the treated data, can be captured in Equation (ref) by a two-factor model where the corresponding loadings change on the treated data.
First, we estimate the factor model from the incomplete control panel $Y^{(0)}$ and obtain $\tilde C^{(0)}=( \tilde \Lambda^{(0)})^{\top} \tilde F$. Second, we use an ordinary least squares regression to estimate the loadings for the treated $\tilde \Lambda_i^{(1)}$,
which yields an estimate for the common components for the treated panel $\tilde{C}^{(1)}_{it} = (\tilde{\Lambda}_i^{(1)})^\top \tilde{F}_t.$\footnote{If units switch between treatment and control, we can modify Equation ((ref)) to $\tilde \Lambda_i^{{(1)}} = \left( \sum_{t \in \mathcal{S}_i} \tilde F_t \tilde F_t^\top \right)^{-1} \sum_{t \in \mathcal{S}_i} \tilde F_t Y^{{(1)}}_{it}$, where $\mathcal{S}_i$ it the set of indices for the treated observations. }
The following theorem shows the asymptotic distributions for $\tilde{C}^{(1)}_{it}$, the individual treatment effect, and the weighted average treatment effect. The asymptotic distributions allow us to construct test statistics for various treatment effects.
The results of Theorem (ref) are a consequence of Theorems (ref) and (ref). The challenge arises from correctly capturing the asymptotic covariance between the estimated treated and control common components. This additional covariance term is due to the correction terms from the missing observations. In Theorem (ref), we impose the additional Assumption (ref) for the general estimator and Assumption (ref) for the probability-weighted estimator. Both simply state that the conventional central limit theorems based on the weak dependencies in the errors apply to the subset of treated time periods. These conditions are automatically satisfied in our simplified model and thus can be neglected, as stated in Proposition (ref) in the Appendix.
Theorems (ref), (ref) and (ref) are formulated with respect to the asymptotic covariances based on the population model. In order to use them in practice we need feasible estimators for the covariance terms. We propose to use the plug-in estimators $\tilde F_t$, $\tilde{\Lambda}_i$ and $\tilde e_{it} = Y_{it} - \tilde{\Lambda}_i^\top \tilde F_t$ for $(H^{-1})^\top F_t $, $H \Lambda_i$ and $e_{it}$. All moments are based on these three objects. For example $\hat \Sigma_F:= \frac{1}{T} \tilde F^{\top} \tilde F$ consistently estimates $(H^{-1})^\top \Sigma_F (H^{-1})$. The rotation matrix $H$ can be ignored in the estimated covariances of the common components and the treatment effects as it cancels out. It is only the distribution of the loadings and factors that are estimated up to a rotation matrix. The challenge is to deal with the time and cross-sectional dependency in the residuals. We impose the additional assumption that the time-series and cross-section covariance matrices of the errors $e_{it}$ are sparse in the sense that only a finite number of row elements are non-zero and we know the indices of the non-zero elements. More specifically we define
and assume that $|\mathcal{E}_t |= O(N)$ and $| \mathcal{E} | = O(NT)$. The estimator for $H \Gamma^{\textnormal{obs}}_{\Lambda,j} H^\top$ and $H \Gamma^{\textnormal{obs}}_{F,t} H^\top$ depend on the dependency structure in the residuals and we propose the plug-in estimator based on only the non-zero moments of the residuals:
The estimators are analogous for the probability-weighted estimator. A special case is the estimation approach in bai2003inferential that assumes independence of the residuals over time and the cross-section and hence only uses the diagonal entries of the residual covariance and autocovariance matrix. Instead of assuming knowledge of the non-zero entries, it is possible to generalize the estimator similar to fan2013large and estimate the non-zero entries with a thresholding estimation approach. We propose a HAC estimator for $\Gamma^{\textnormal{miss}}_{\Lambda,j}$, $\Gamma^{\textnormal{miss}}_{F,t}$ and $\Gamma^{\textnormal{miss, cov}}_{\Lambda, F, j, t}$ to account for the time-series dependency in the factors similar to bai2003inferential.
Hence, the treatment effects normalized by their estimated standard deviations follow asymptotically a standard normal distribution, and we obtain feasible test statistics for the various treatment effects.
Our results can be generalized to the case where the number of observed entries is not proportional to $N$ or $T$ but grows at a strictly smaller rate. The general arguments of the proofs stay the same but we need to carefully account for the convergence rates of each term based on the set $\mathcal{Q}_{ij}$. The mean squared consistency of the estimated loadings in Theorem (ref) generalizes to \[ \frac{1}{N} \sum_{j = 1}^{N} \left\lVert\tilde{\Lambda_j} - H \Lambda_j\right\rVert^2 = O_p \Bigg( \max \bigg( \frac{1}{N}, \,\, \frac{1}{N^2} \sum_{i = 1}^{N} \sum_{j = 1}^N \frac{1}{|\mathcal{Q}_{ij} | } \bigg)\Bigg). \] Moreover, we can show the asymptotic normality of the estimated loadings $\tilde \Lambda$, factors $\tilde F$ from the equally weighted regression (ref), and common components $\tilde C$ under similar assumptions as those in Theorem (ref). The estimated loadings $\tilde{\Lambda}_j$ are asymptotically normal with convergence rate \[ H^{-1} \tilde{\Lambda}_j - \Lambda_j = O_p \Bigg( \bigg[\max \bigg( \frac{1}{N} \sum_{i = 1}^{N} \frac{1}{|\mathcal{Q}_{ij}|}, \,\, \frac{1}{N^2} \sum_{i = 1}^{N} \sum_{l=1}^{N} \frac{|\mathcal{Q}_{ij} \cap \mathcal{Q}_{lj} |}{|\mathcal{Q}_{ij}| |\mathcal{Q}_{lj}| } \bigg) \bigg]^{1/2} \Bigg), \] where the second term is closely related to $\omega_{jj}$ defined in Assumption (ref). The estimated factors $\tilde{F}_t$ are asymptotically normal with convergence rate \[H^\top \tilde{F}_t - F_t = O_p \Bigg( \bigg[\max \bigg( \frac{1}{\sum_{i = 1}^{N} W_{it}}, \,\, \frac{1}{N^4} \sum_{i = 1}^{N} \sum_{j=1}^{N} \sum_{k = 1}^{N} \sum_{l = 1}^{N} \frac{|\mathcal{Q}_{ij} \cap \mathcal{Q}_{kl} |}{|\mathcal{Q}_{ij}| |\mathcal{Q}_{kl}| } \bigg) \bigg]^{1/2} \Bigg), \] where the second term is closely related to $\omega$ defined in Assumption (ref). Similarly, by combining the rates of estimated factors and loadings, the estimated common components $\tilde{C}_{it}$ have an asymptotic normal distribution with rate
The last term is closely related to $\omega_j$ defined in Assumption (ref). The expression for the asymptotic covariances of the estimators become more complex. The proofs for the consistency and asymptotic normality for the general case, when observed entries are not proportional to $N$ and $T$, are very similar to the proofs of Theorems (ref) and (ref), but just require carefully keeping track of the convergence rates of each term.\footnote{The proofs are available upon request.}
We illustrate the more general convergence rates in the simultaneous treatment observation pattern in Table (ref), where we can provide explicit expressions for the different rates. The mean square consistency result of the loadings simplifies to \[ \frac{1}{N} \sum_{j = 1}^{N} \left\lVert\tilde{\Lambda}_j - H \Lambda_j\right\rVert^2 = O_p \Bigg( \max \bigg( \frac{1}{N}, \,\, \frac{N_0}{N T_0}, \,\, \frac{1}{T} \bigg)\Bigg). \] We obtain two different convergence rates for the estimated loadings:
Similarly, the estimated factors have two different convergence rates depending on which time block we consider:
This results in four different convergence rates for each block for the estimated common components:
In this section, we demonstrate the finite sample properties of our asymptotic results for both the observed entries and the missing entries. We confirm the theoretical distribution results for the estimated factor, loadings, common components, and treatment effects. We generate the data from a one-factor model $X_{it} = \Lambda_i F_t + e_{it}$, where $F_t \stackrel{\mathrm{iid}}{\sim} \mathcal{N}(0,1)$, $\Lambda_i \stackrel{\mathrm{iid}}{\sim} \mathcal{N}(0,1)$ and $e_{it} \stackrel{\mathrm{iid}}{\sim} \mathcal{N}(0,1)$. The observation pattern depends on unit-specific characteristics $S_i = \mathbbm{1}(\Lambda_i \geq 0 )$, which are a function of the factor loadings. We study two observation patterns which are illustrated in Figures (ref)(a) and (ref)(b):
To conserve space, we report here the distribution results for the regression based estimator based on Equation (ref), but the results extend to the propensity-weighted estimator. Figure (ref) shows the histograms of standardized factors, loadings, and common components for randomly selected observed entries and missing entries based on Theorem (ref). The histograms match the standard normal density function very well and support the validity of our asymptotic results in finite samples.
Figure (ref) confirms that our treatment test in Theorem (ref) has the correct size. The control data follows our benchmark one-factor model. We assume a constant treatment effect, i.e., $\Lambda_i^{(1)} = \Lambda_i^{(0)} + \tau$, where $\tau$ is set to 0 or 0.25. Figure (ref) shows the histograms of standardized common components for treated and control, the individual treatment effect, and an equally weighted treatment effect for randomly selected units and times. As expected, the histograms support the validity of our asymptotic results in finite samples.
Table (ref) demonstrates the statistical power of our tests for individual and average treatment effects, where the null hypotheses are $\mathcal{H}_0: \beta^{(1)}_i - \beta^{(0)}_i = 0$ with equal weights for all time periods, i.e., $\beta^{(1)}_i=\tau^{(1)}_i$ and $\beta^{(0)}_i=\tau^{(0)}_i$. The power increases with the data dimensionality ($N$ and $T$) and the scale of treatment effect that is determined by the mean of the factor $\mu_F$ and the difference between the control and treated loadings $\Lambda_i^{(1)} - \Lambda_i^{(0)}$. The null hypothesis implies $\Lambda_i^{(1)} - \Lambda_i^{(0)}=0$, which we use in the estimation of the asymptotic variance. This slightly improves the power, but the results in the Internet Appendix show that we also have good power properties without imposing the null hypothesis in the estimation of the asymptotic covariances. Moreover, the statistical power increases with the proportion of observed entries, as shown in the comparison between Tables (ref) and (ref) in the Internet Appendix.
In this section, we show that our benchmark regression-based estimator (denoted as XP) and propensity-weighted estimator (denoted as $\text{XP}_{\text{PROP}}$) perform well under a variety of missing patterns. As reference we also include the estimators of jin2020factor (denoted as JMS) and bai2019matrix (denoted as BN). Each of the two estimators is designed for a specific observation pattern and hence provides a natural reference level for that specific pattern. jin2020factor assume that observations are missing at random, while bai2019matrix is tailored to an observation pattern with a block structure after proper reshuffling. These are the four estimation approaches that provide an inferential theory for imputed common components in an approximate factor model and were available at the time of submission of this paper.
We generate the data from a two-factor model $X_{it} = \Lambda_i^{\top} F_t + e_{it}$, where $F_t \stackrel{\mathrm{iid}}{\sim} \mathcal{N}(0,I_2)$, $\Lambda_i \stackrel{\mathrm{iid}}{\sim} \mathcal{N}(0,I_2)$ and $e_{it} \stackrel{\mathrm{iid}}{\sim} \mathcal{N}(0,1)$. We consider six different observation patterns. The first three cases are (1) missing uniformly at random, (2) simultaneous treatment adoption, and (3) staggered treatment adoption. Then, we allow the observation pattern for these three cases to depend on a unit-specific characteristic defined as $S_i = \mathbbm{1}(\Lambda_{i,2} \geq 0 )$. Hence, case four to six are (4) missing at random conditional on $S_i$, (5) simultaneous treatment adoption conditional on $S_i$, (6) staggered treatment adoption conditional on $S_i$. Table (ref) contains figures showing the observation patterns and their detailed descriptions. Note that these are all practically relevant patterns, in particular the staggered treatment adoption that appears in our empirical companion paper and is prevalent in empirical applications.
Table (ref) compares the performance of estimating the common components. We report the normalized mean squared error (MSE) of the four methods for observed, missing and all units defined as follows: \[\text{relative MSE}_{\mathcal{S}} = \frac{\sum_{(i,t) \in \mathcal{S}} \big(\tilde{C}_{it} - C_{it}\big)^2 }{\sum_{(i,t) \in \mathcal{S}} C_{it}^2}, \] where $\mathcal{S}$ is either the set observed, missing or all observations.
First, and most importantly, our benchmark estimator shows excellent performance for all observation patterns. Our estimator has the smallest or at least a very similar small MSE compared to the other methods, as indicated by the bold numbers. Hence, we view our approach as a simple and reliable all-purpose estimator. Our propensity-weighted estimator is very close to the benchmark estimator but performs slightly worse. This is in line with our theoretical result that propensity weighting is generally less efficient.
In the case of missing at random conditional or unconditional on $S_i$, our methods have the smallest MSE. jin2020factor also have a small MSE as long as the observation pattern does not depend on $S_i$ as their method is designed for missing uniformly at random. Missing at random violates the assumptions of bai2019matrix and their estimator not applicable as there not sufficiently large blocks of fully observed entries.
In the case of simultaneous treatment adoption, bai2019matrix has the smallest MSE as their method is tailored to this case. Interestingly, our method as an all-purpose estimator is very close to bai2019matrix. When the observation pattern depends on $S$, it can shrink the size of fully observed blocks, which increases the importance of using all observed entries resulting in the smallest MSE for our method. In the case of simultaneous treatment adoption, the assumptions in jin2020factor are violated, which is reflected in the larger MSE.
Our methods have the smallest MSE for the case of staggered treatment adoption that is prevalent in empirical applications athey2018design. This holds whether the observation pattern depends or does not depend on $S$. In contrast to bai2019matrix, we use all observed entries in the estimation, which provides a more efficient estimator. Note that in this simulation example the fully observed blocks are very small, and hence, similar to the missing at random case, the assumptions in bai2019matrix might not be satisfied. As the assumptions of jin2020factor are violated, their imputation results in larger errors.
The Internet Appendix shows that the findings are robust to the size of the panel and the parameters of the observation patterns. We also compare the MSE of the various methods after iterations in Tables (ref)-(ref) in the Internet Appendix. In more detail, we first impute the missing values with different methods. In the second step, we apply PCA to the full panel with imputed values to estimate the factor model and update the imputed values with the estimated common components. The observed entries stay the same. This process is repeated for multiple iterations. Note that this iterated estimation approach is actually a different estimation approach by itself. The four methods provide different starting values for the same iterative estimation approach that is based on a fixed-point argument. Importantly, there is no inferential theory for iterative estimators under general patterns.\footnote{While jin2020factor consider iterations, their asymptotic results only hold for missing at random. bai2019matrix provide distribution results for a different iteration that is not making use of all observations and therefore only has a minor effect.} Hence, if the goal is to estimate treatment effects, these iterative estimators cannot be used. Since our methods start with a value that has a smaller MSE, our methods, in general, converge faster (often already after three iterations) and also have a small MSE for a fixed number of iterations. Our results are robust to the choice of $N$ and $T$ and we present the corresponding results for $N = 100$ and $T = 150$ in Tables (ref)-(ref) in the Internet Appendix. In summary, if the goal is to only minimize the imputation error without an inferential theory, the iterative estimation generally improves the results, but the relative performance of the different estimation approaches without iteration carries over to the iteration setup.
In this section, we show that the propensity-weighted estimator can have desirable robustness properties under misspecification. Our results are motivated by insights from causal inference that propose doubly robust estimation procedures for missing values, as discussed, for example, in kang2007. In causal inference, we can either model the relationship between the covariates and the outcome or model the probabilities of missingness to estimate causal effects. Doubly robust procedures combine both by using a propensity weight in regressions to mitigate the selection bias. Their potential advantage is that they can provide reliable estimates in the case of omitted variables. Our setup differs from classical causal inference as we estimate the covariates as latent factors from the data. However, we can have a situation similar to omitted variables if we estimate too few latent factors, the factors are weak, or the population model is nonlinear.
We compare our benchmark estimator (XP) and propensity-weighted estimator ($\text{XP}_{\text{PROP}}$) under two types of model misspecification. In Table (ref), we consider the case of omitted factors. The population model is generated by a two-factor model, but we only estimate one latent factor. In this case, the propensity-weighted estimator can perform better than the benchmark estimator. However, when the model is correctly specified, and we estimate two factors, the benchmark estimator dominates. When the second factor is weak in the sense that its variance and corresponding eigenvalue are very small, the situation is similar to an omitted factor. In this case, it is possible that the propensity-weighted estimator performs better even if we estimate the correct number of latent factors. Note that weak factors are also a form of misspecification, as discussed in onatski2012. In this simulation, observations are more likely to miss if they are exposed to the omitted or weak second factor. Hence, the robustness gains of the propensity-weighted estimator arise for the missing data and the treatment effects.
The case of omitted latent factors shares similarity with the case of a misspecified functional form. In Table (ref) in the Online Appendix we generate the data from a non-linear one-factor model. Under certain assumptions it is possible to approximate a non-linear transformation as a linear function of appropriate basis functions. Such an approximation can be formulated as a linear latent multi-factor model, where the additional factors are non-linear transformations of the underlying one-factor model. Hence, some form of functional model misspecification can be corrected by using more latent factors. In our example, the non-linearity is very well approximated by three latent factors for the simple and propensity-weighted estimator. However, if we use only one or two latent factors, the propensity-weighted estimator is more robust to the misspecification. While we do not provide a formal non-parametric theory, our simulation suggests that a non-linear misspecification can share similar features with the case of omitted factors. Therefore, if a researcher suspects some form of model misspecification, the propensity-weighted estimator can be a useful alternative.\footnote{For a non-linear factor model, feng2020causal proposes a local PCA method that uses a linear model approximation in a local neighborhood. Our argument is based on a global approximation of the non-linear functional relationship, where the additional latent factors serve as additional basis functions.}
This paper develops the inferential theory for latent factor models estimated from large dimensional panel data with missing observations. Our paper stands out by the generality of the missing patterns that we allow for. We propose two estimators for the latent factor model: a simple all-purpose estimator and an extension to a probability-weighted estimator. Our all-purpose estimator is easy to use while it performs well under a variety of missing patterns. The propensity weighted estimator is an alternative that is less efficient for correctly specified models but can be more robust to certain forms of misspecification. The key application of our asymptotic distribution theory is to test causal treatment effects. We provide a test for the point-wise treatment effect that can be heterogeneous and time-dependent under general adoption patterns where the units can be affected by unobserved factors.