Large Dimensional Latent Factor Modeling with Missing Observations and Applications to Causal Inference

Introduction

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, Bai and Ng (2002) show that the factor model can be estimated with PCA applied to the covariance matrix of the data. Bai (2003) and Fan, Liao, and Mincheva (2013) 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 Bai (2009), sparse and interpretable latent factors in Pelger and Xiong (2021a), time-varying loadings in Fan, Liao, and Wang (2016) and Pelger and Xiong (2021b), high-frequency estimation in Pelger (2019) and including additional moments to estimate weak factors as in Lettau and Pelger (2020a). 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 Jin, Miao, and Su (2021), Bai and Ng (2021), and Cahan, Bai, and Ng (2021). 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. Jin, Miao, and Su (2021) 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 (Rubin, 1976; Dempster, Laird, and Rubin, 1977).Stock and Watson (2002b); Bańbura and Modugno (2014); Negahban and Wainwright (2012) propose to use EM algorithms to estimate the factor model from panel data with missing observations. Giannone, Reichlin, and Small (2008); Doz, Giannone, and Reichlin (2011); Jungbacker, Koopman, and Van der Wel (2011); Stock and Watson (2016) propose to use the state-space framework and Kalman Filtering to estimate the factor model with missing observations. Gagliardini, Ossola, and Scaillet (2019) propose a simple diagnostic criterion for an approximate factor structure in large (unbalanced) panel data sets. Bai and Ng (2021) 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. Cahan, Bai, and Ng (2021) 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 Jin, Miao, and Su (2021) for data missing at random, and as Bai and Ng (2021) 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 (Mazumder, Hastie, and Tibshirani, 2010; Negahban and Wainwright, 2011, 2012). 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, Chen, Fan, Ma, and Yan (2019) 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 Chen, Fan, Ma, and Yan (2019). 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. Abadie, Diamond, and Hainmueller (2010, 2015), Doudchenko and Imbens (2016), Xu (2017), Li and Bell (2017) and Li (2019) propose to match each treated unit by weighted averages of all control units using the pretreatment observations. Li and Bell (2017), Li (2019) and Masini and Medeiros (2018) 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. Athey, Bayati, Doudchenko, Imbens, and Khosravi (2021) 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 Kang and Schafer (2007) 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 Hirano, Imbens, and Ridder (2003).

The rest of the paper is organized as follows. Section 2 introduces the model and provides the simple all-purpose estimator for factors, loadings, and common components. Section 3 states the necessary assumptions for the asymptotic distribution results that are presented in Section 4. Sections 5 and 6 extend the results to the propensity weighted estimator. Section 7 shows how to apply our model to test treatment effects. We discuss the feasible estimation in Section 8 and how to relax the rate conditions in Section 9. The extensive simulation in Section 10 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.

Model and Estimation

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.

2 Missing Observations

We allow for very general patterns in the missing observations. Figure 1 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.See Candès and Recht (2009); Zhou, Wilkinson, Schreiber, and Pan (2008) for the Netflix challenge and Athey, Bayati, Doudchenko, Imbens, and Khosravi (2021); Athey and Imbens (2021) 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 S1 states the conditions on the observation pattern.

For a given observation matrix $W$, $\frac{|\mathcal{Q}_{ij}|}{T}\geq\underline{q}>0$ and there exist constants $q_{ij}$ and $q_{ij,kl}$ for all $i,j,k,l$ such that $q_{ij}=\lim_{T\rightarrow\infty}\frac{|\mathcal{Q}_{ij}|}{T}$ and $q_{ij,kl}=\lim_{T\rightarrow\infty}\frac{|\mathcal{Q}_{ij}\cap\mathcal{Q}_{kl}|}{T}$.

Assumption S1 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 Rosenbaum and Rubin (1983). Assumption S1 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:

Missing at random: $P(W_{it}=1)=p$ for all $i$ and $t$. In this case all units and times are equally likely to be observed.

Cross-section missing at random: $P(W_{it}=1)=p_{t}$. For each $t$ each cross-sectional unit is equally likely to miss.

Time-series missing at random: $P(W_{it}=1)=p_{i}$. For each $i$ each time observation is equally likely to miss.

Cross-section and time-series dependency: $P(W_{it}=1)=p_{it}$, which allows for different probabilities for each unit and time.

Staggered treatment adoption: If $W_{it}=0$ then $W_{it^{\prime}}=0$ for all $t^{\prime}\geq t$. This is a special case of 4. with $P(W_{it}=1)=p_{it}$. For the special case that the probability does not depend on $i$, the staggered design is a special case of cross-section missing at random $P(W_{it}=1)=p_{t}$.

Mixed frequency observations: Each cross-section unit has a fixed known observation pattern over time. This can be modeled as one random draw for each cross-section unit to assign it to a specific pattern. A feasible model approach uses $P(W_{it}=1)=p_{t}$ as this is another special case of cross-section missing at random.

3 Estimator

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.

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:

4 Illustration

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 (3) corresponds to a weighted least square regression which provides the correct estimator:

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.

Assumptions

We assume an approximate factor structure at the same level of generality as in Bai (2003). 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 G2 and G3 are collected in the Appendix. In the main text, we present a simplified factor model with the stronger Assumptions S2 and S3, 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 S2 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 S3.

There exists a positive constant $M<\infty$ such that:

Independence: $F$, $\Lambda$ and $e$ are independent.

Eigenvalues: The eigenvalues of $\Sigma_{\Lambda}\Sigma_{F}$ are distinct.

Systematic loadings: $\frac{1}{N}\sum_{i=1}^{N}\Lambda_{i}\Lambda_{i}^{\top}W_{it}\overset{P}{\rightarrow}\Sigma_{\Lambda,t}$ for some positive definite matrix $\Sigma_{\Lambda,t}$ for any $t$.

Dependency in missing pattern: $\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{l=1}^{N}\frac{q_{ij,lj}}{q_{ij}q_{lj}}\overset{P}{\rightarrow}\omega_{jj}$, $\lim_{N\rightarrow\infty}\frac{1}{N^{3}}\sum_{i=1}^{N}\sum_{l=1}^{N}\sum_{k=1}^{N}\frac{q_{li,kj}}{q_{li}q_{kj}}\overset{P}{\rightarrow}\omega_{j}$ and $\lim_{N\rightarrow\infty}\frac{1}{N^{4}}\sum_{i=1}^{N}\sum_{l=1}^{N}\sum_{j=1}^{N}\sum_{k=1}^{N}\frac{q_{li,kj}}{q_{li}q_{kj}}\overset{P}{\rightarrow}\omega$ for all $j$ and some constants $\omega_{jj},\omega_{j},\omega$.

Assumption S3 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 3 in the Appendix, the simplified model is just a special case of the general approximate factor model specified by Assumptions G2 and G3. The simplified Assumption S2 implies the general Assumption G2, while Assumption S3 combined with the other simplified assumptions implies the general Assumption G3.

We assume that the number of factors $r$ is consistently estimated. For example the criteria developed in Bai and Ng (2002) 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 Jin, Miao, and Su (2021) or an eigenvalue ratio argument as in Ahn and Horenstein (2013) to general missing patterns. Given a consistent estimator for the number of factors, we can treat $r$ as known.

Asymptotic Results

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.The estimator in Bai and Ng (2021) can avoid this term by neglecting partially observed entries, which means that in general, they are using less information. The estimator of Bai and Ng (2021) 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$. Cahan, Bai, and Ng (2021) 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.

Define $\delta_{NT}=\min(N,T)$. Under Assumptions S1 and G2 it holds that

Theorem 1 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 Bai and Ng (2002). Theorem 1 is based on the assumption that the observed entries are representative of the missing entries and hence provide a consistent estimation. Theorem 1 is a critical intermediate step to show the asymptotic normality of the estimated factor model in the next section.

2 Asymptotic Normality

The factors, loadings, and common components are asymptotically normally distributed. Indeed, Theorem 2 states that the asymptotic distributions have two parts: First, we recover the asymptotic variance that is identical to the conventional PCA in Bai (2003) 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 2 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.In the asymptotic distribution, we apply the rotation matrices to the estimated loadings and factors instead of their population values as in Bai (2003). Obviously, these two representations are equivalent and can be easily transformed into each other. Our choice of representation was made for exposition purposes only.

Under Assumptions S1, G2 and G3 and for $N,T\rightarrow\infty$ we have for each $j$ and $t$:

For $\sqrt{T}/N\rightarrow 0$ the asymptotic distribution of the loadings is

with $\Gamma^{\textnormal{miss}}_{\Lambda,j}=h_{j}(\Lambda_{j})$. $\Gamma^{\textnormal{obs}}_{\Lambda,j}$ and the function $h_{j}(\cdot)$ are defined in Assumptions G3.3 and G3.5.

For $\sqrt{N}/T\rightarrow 0$ and $\sqrt{T}/N\rightarrow 0$, the asymptotic distribution of the factors is

with $\Gamma^{\textnormal{miss}}_{F,t}=g_{t}(F_{t})$. $\Gamma^{\textnormal{obs}}_{F,t}$ and the function $g_{t}(\cdot)$ are defined in Assumptions G3.4 and G3.5.

The asymptotic distribution of the common component is

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 Bai (2003). 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 (8).

The distribution results of Theorem 2 simplify under Assumptions S2 and S3, 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.

Suppose Assumptions S1, S2 and S3 hold and $N,T\rightarrow\infty$. Then Theorem 2 holds. If in addition, $q_{ij}$ and $q_{ij,kl}$ are independent of $\Lambda_{m}\Lambda_{m}^{\top}$ for all $i,j,k,l,m$, then the asymptotic variances simplify as follows with the weights $\omega,\omega_{j}$ and $\omega_{jj}$ defined in Assumption S3:

The asymptotic variance of the loadings in formula (6) simplifies to

The asymptotic variance of the factors in formula (7) simplifies to

The asymptotic variance of the common component in formula (8) simplifies to

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 Bai (2003). 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.

If all entries are observed at random with equal probability, we can use the approach of Jin, Miao, and Su (2021) to estimate the factor model and impute the missing entries. We compare the efficiency of our approach with the one of Jin, Miao, and Su (2021). For a direct comparison, we follow the order of estimation in Jin, Miao, and Su (2021) and switch the role of factors and loadings in our all-purpose estimator: We first estimate the factors from the time-series sample covariance matrix, and then estimate the loadings from a time-series regression of the observed outcomes on the estimated factors.

Suppose Assumptions S1, S2 and S3 hold and that every entry is randomly missing with observed probability $p$. We switch the role of factors and loadings in the all-purpose estimator. As $N,T\rightarrow\infty$, it holds that:

Propensity Weighted Estimator

$W$ is independent of $\Lambda$ conditional on $S$.

For any $i$ and $j$ satisfying $i\neq j$, and for any $t$ and $s$, $W_{it}$ is independent of $W_{js}$ conditional on $S_{i}$ and $S_{j}$ where $t$ and $s$ can be the same. The probability of $W_{it}=1$ depends on $S_{i}$ and satisfies $P(W_{it}=1|S_{i})\geq\underline{p}>0$.

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 C1.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.See (Rosenbaum and Rubin, 1983) 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.

For any $i$, $\Lambda_{i}$ is independent of $S_{j}$ conditional on $S_{i}$ for $j\neq i$. Moreover, for any $i$ and $j$ satisfying $i\neq j$, $\Lambda_{i}$ is independent of $\Lambda_{j}$ conditional on $S_{i}$ and $S_{j}$.

Proposition 3 in the Appendix shows that the simplified model is just a special case of the general approximate factor model specified by Assumptions GC2 and GC3. The simplified Assumption S2 combined with Assumptions S1, C1 and S2 imply the general conditional Assumption G2, while Assumption S3 combined with the other simplified Assumptions S1, C1, S2, S3.2 and C2 imply the general Assumption GC3.

Under Assumptions S1, C1, G2, GC2 and GC3 and for $N,T\rightarrow\infty$ we have for each $j$ and $t$:

The asymptotic distribution of the loadings is the same as in Theorem 2.

For ${\sqrt{N}}/{T}\rightarrow 0$ and $\sqrt{T}/N\rightarrow 0$, the asymptotic distribution of the factors is

with $\Gamma^{\textnormal{miss},S}_{F,t}=g^{S}_{t}(F_{t})$. $\Gamma^{\textnormal{obs},S}_{F,t}$ and $g^{S}_{t}(\cdot)$ are defined in Assumptions GC3.4 and GC3.5.

The asymptotic distribution of the common components is

Suppose Assumptions S1, C1, S2, S3.2, C2, and C3 hold and $N,T\rightarrow\infty$. Then Theorem 3 holds. If in addition, $q_{ij}$ and $q_{ij,kl}$ are independent of $\Lambda_{m}\Lambda_{m}^{\top}$ for all $i,j,k,l,m$, then the asymptotic variances simplify as follows with the weights $\omega,\omega_{j}$ and $\omega_{jj}$ defined in Assumption S3:

The asymptotic variance of the factors in formula (9) simplifies to

The asymptotic variance of the common component in formula 10 simplifies to

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 1.1.

2 Robustness to Model Misspecification

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.Robins, Rotnitzky, and Zhao (1994); Robins and Rotnitzky (1995) among others discuss the reduction of bias from omitting regressors. Robins and Wang (2000); Kang and Schafer (2007); Robins, Sued, Lei-Gomez, and Rotnitzky (2007) show the large variance of propensity weighted regressions, and Freedman and Berk (2008) 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$.Given our distribution theory, it is relatively straightforward to extend the consistent estimator of Bai and Ng (2002) 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 (Pelger, 2019; Lettau and Pelger, 2020b). 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 10.3.

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}=\bm{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 1(a). 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:

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.Simon (1954) refers to this type of selection bias as a spurious correlation.

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

Feasible Estimator of the Probability Weighting

We provide feasible estimators for $P(W_{it}=1|S_{i})$ which we need in Equation (4) 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}=\widehat{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 (4), which yields the following decomposition:

We replace the propensity score in $p_{it}$ in Equation (4) by its estimate $\hat{p}_{it}$.

The estimates of the loadings do not depend on the propensity score. Hence, Theorem 1 and the asymptotic distribution of the loadings in Theorem 3 continue to hold independently of $\hat{p}_{it}$.

The following holds for the distribution of the factors and common components.

If $\max_{i}|\hat{p}_{it}-p_{it}|=o_{P}(1)$, then the factors and common components are estimated consistently pointwise under the assumptions of Theorem 3.

If $\max_{i}|\hat{p}_{it}-p_{it}|=o_{P}\left(\frac{1}{N^{1/4}}\right)$, then Theorem 3 continues to hold as it is.

We discuss feasible estimators for the most important cases of missing patterns which are summarized in Table 2. 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 (2). 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 4.2(b) applies. If $p(S_{i})$ is estimated non-parametrically with a kernel with bandwidth $h$, the convergence rate is typically $\sqrt{NTh}$ 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 4.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 3. 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 3 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 4.2(b).

Tests of Treatment Effects

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

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 Bai and Ng (2021). Given our asymptotic distribution theory for the common component, we can provide the asymptotic distribution of the individual and average treatment effects analogously to Bai and Ng (2021). 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:

Individual treatment effect: $\tau_{it}=C_{it}^{(1)}-C_{it}^{(0)}$

Average treatment effect over time: $\tau_{i}=\frac{1}{T_{1,i}}\sum_{t=T_{0,i}+1}^{T}\tau_{it}$

Weighted average treatment effect: $\tau_{\beta,i}=\beta_{i}^{(1)}-\beta_{i}^{(0)}$ where $\beta_{i}$ are the regression coefficients on some covariates $Z$:

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 (12) 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 2 and 3. 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 Abadie, Diamond, and Hainmueller (2010, 2015), 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 (12). For example in the simultaneous treatment case in Table 1, 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 (12) with one factor that changes on the treated data, can be captured in Equation (14) by a two-factor model where the corresponding loadings change on the treated data.

Suppose Assumptions S1, G2, G3 and G4 hold and the control and treated panel share the same factors. For $\delta_{NT_{i}}=\min(N,T_{1,i})$, as $\delta_{NT_{i}}\rightarrow\infty$ the following holds:

The asymptotic distribution for the common component is

with $\Gamma^{\textnormal{obs}}_{F,t}$ and $\Gamma^{\textnormal{miss}}_{F,t}$ given in Theorem 2, $\Gamma^{\textnormal{obs},(1)}_{\Lambda,i}=\Sigma_{F,e_{i}}$, \Gamma^{\textnormal{miss},(1)}_{\Lambda,i}=\Sigma_{\Lambda}^{-1}\Big{[}\frac{1}{T_{1,i}^{2}}\sum_{u,s=T_{0,i}+1}^{T}g_{u,s}(\Sigma_{\Lambda,u}^{-1}\Lambda_{i}^{(1)},\Sigma_{\Lambda,s}^{-1}\Lambda_{i}^{(1)})\Big{]}\Sigma_{\Lambda}^{-1}, \Gamma^{\textnormal{miss, cov},(0),(1)}_{\Lambda,F,i,t}=\Sigma_{\Lambda}^{-1}\Big{[}\frac{1}{T_{1,i}}\sum_{u=T_{0,i}+1}^{T}g_{u,s}(\Sigma_{\Lambda,u}^{-1}\Lambda_{i}^{(1)},\Sigma_{\Lambda}^{-1}\Sigma_{F}^{-1}F_{t})\Big{]}, and the function $g_{u,s}(\cdot,\cdot)$ is defined in Assumption G4.

The asymptotic distribution for the individual treatment effect is

The results of Theorem 5 are a consequence of Theorems 2 and 3. 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 5, we impose the additional Assumption G4 for the general estimator and Assumption GC4 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 3 in the Appendix.

Feasible Estimation and Testing

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 Bai (2003) 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 Fan, Liao, and Mincheva (2013) 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 Bai (2003).

Suppose that the assumptions of Theorems 2, 3 or 5 hold. In addition, we assume that the time-series and cross-section covariance matrices of the errors $e_{it}$ are sparse in the sense that $|\mathcal{E}_{t}|=O(N)$ and $|\mathcal{E}|=O(NT)$ and we know the non-zero elements. Then, the plug-in estimators of the asymptotic covariances in Theorems 2, 3 and 5 are consistent and the asymptotic statements in the respective theorems continue to hold with the estimated covariance matrices.

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.

Generalization of the Missing Patterns

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 1 generalizes to

The last term is closely related to $\omega_{j}$ defined in Assumption S3. 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 1 and 2, but just require carefully keeping track of the convergence rates of each term.The proofs are available upon request.

We illustrate the more general convergence rates in the simultaneous treatment observation pattern in Table 1(a), where we can provide explicit expressions for the different rates. The mean square consistency result of the loadings simplifies to

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:

Simulation

Missing at random: Entries are observed independently with probability 0.75 if $S_{i}=1$, and 0.5 if $S_{i}=0$.

Simultaneous treatment adoption: Once a unit adopts treatment, it stays treated afterward. For the units with $S_{i}=1$, $25\%$ randomly selected units adopt the treatment from time $0.75\cdot T$ and the remaining $75\%$ units stay in the control group until the end. For the units with $S_{i}=0$, $62.5\%$ randomly selected units adopt the treatment from time $0.375\cdot T$ and the remaining $37.5\%$ units stay in the control group until the end. We model the treated data as missing.

To conserve space, we report here the distribution results for the regression based estimator based on Equation (3), but the results extend to the propensity-weighted estimator. Figure 2(d) shows the histograms of standardized factors, loadings, and common components for randomly selected observed entries and missing entries based on Theorem 2. The histograms match the standard normal density function very well and support the validity of our asymptotic results in finite samples.

Figure 3 confirms that our treatment test in Theorem 5 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 3 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 3 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 3 and 14 in the Internet Appendix.

2 Robustness to Missing Patterns

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 Jin, Miao, and Su (2021) (denoted as JMS) and Bai and Ng (2021) (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. Jin, Miao, and Su (2021) assume that observations are missing at random, while Bai and Ng (2021) 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.

Table 4 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:

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. Jin, Miao, and Su (2021) 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 Bai and Ng (2021) and their estimator not applicable as there not sufficiently large blocks of fully observed entries.

In the case of simultaneous treatment adoption, Bai and Ng (2021) has the smallest MSE as their method is tailored to this case. Interestingly, our method as an all-purpose estimator is very close to Bai and Ng (2021). 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 Jin, Miao, and Su (2021) 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 (Athey and Imbens, 2021). This holds whether the observation pattern depends or does not depend on $S$. In contrast to Bai and Ng (2021), 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 Bai and Ng (2021) might not be satisfied. As the assumptions of Jin, Miao, and Su (2021) 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 6-8 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.While Jin, Miao, and Su (2021) consider iterations, their asymptotic results only hold for missing at random. Bai and Ng (2021) 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 9-12 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.

3 Misspecification and Robustness of Propensity-Weighted Estimator

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 Kang and Schafer (2007). 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 5, 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 Onatski (2012). 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 19 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.For a non-linear factor model, Feng (2020) 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.

Conclusion

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.

Appendix

Let $M<\infty$ denote a generic constant. Let $\left\lVert v\right\rVert$ denote the vector norm and $\left\lVert A\right\rVert=trace(A^{\top}A)^{1/2}$ the Frobenius norm of matrix $A$.

General Assumptions

Time and cross-section dependence and heteroskedasticity of errors: There exists a positive constant $M<\infty$, such that for all $N$ and $T$:

Weak dependence between factor and idiosyncratic errors: for every $(i,j)$,

Eigenvalues: The eigenvalues of $\Sigma_{\Lambda}\Sigma_{F}$ are distinct.

$\frac{\sqrt{T}}{N}\sum_{i=1}^{N}\Lambda_{i}\Lambda_{i}^{\top}\frac{1}{|\mathcal{Q}_{ij}|}\sum_{t\in\mathcal{Q}_{ij}}F_{t}e_{it}\xrightarrow{d}\mathcal{N}(0,\Gamma^{\textnormal{obs}}_{\Lambda,j})$ for every $j$.

$\frac{1}{\sqrt{N}}\sum_{i=1}^{N}W_{it}\Lambda_{i}e_{it}\xrightarrow{d}\mathcal{N}(0,\Gamma^{\textnormal{obs}}_{F,t})$ for every $t$.

We define the filtration $\mathcal{G}^{t}=\sigma(\cup_{s=1}^{T}\mathcal{G}^{t}_{Ts})$ with $\mathcal{G}^{t}_{Ts}=\sigma(\{W_{ij},j\leq s,\text{all }i\},\Lambda,v_{t})$ generated by $\{W_{ij},j\leq s,\text{all }i\}$, $\Lambda$ and $v_{t}$, which is given by $v_{t}=\Sigma_{\Lambda}^{-1}\Sigma_{F}^{-1}F_{t}$. For every $i$ and $t$, and $u_{i}=\Lambda_{i}$, it holds

where X_{i}=\frac{1}{N}\sum_{l=1}^{N}\Lambda_{l}\Lambda_{l}^{\top}\Big{(}\frac{1}{|\mathcal{Q}_{li}|}\sum_{s\in\mathcal{Q}_{li}}F_{s}F_{s}^{\top}-\frac{1}{T}\sum_{s=1}^{T}F_{s}F_{s}^{\top}\Big{)} and $\mathbf{X}_{t}=\frac{1}{N}\sum_{i=1}^{N}W_{it}X_{i}\Lambda_{i}\Lambda_{i}^{\top}$.

$\frac{1}{\sqrt{T_{1,i}}}\sum_{T_{0,i}+1}^{T}F_{t}e_{it}\xrightarrow{d}\mathcal{N}(0,\Sigma_{F,e_{i}})$.

Assumption G3.5 holds for $v_{t}$ equal to $\Sigma_{\Lambda,t}^{-1}\Lambda^{(1)}_{i}$ and $\Sigma_{\Lambda,t}^{-1}(\Lambda^{(1)}_{i}-\Lambda^{(0)}_{i})$ under the filtration $\mathcal{G}=\sigma(\cup_{s=1}^{T}\mathcal{G}_{Ts})$ with $\mathcal{G}_{Ts}=\sigma(\{W_{ij},j\leq s,\text{all }i\},\Lambda)$ generated by $\{W_{ij},j\leq s,\text{all }i\}$ and $\Lambda$.

$\frac{\sqrt{T}}{N}\sum_{i=1}^{N}\Lambda_{i}\Lambda_{i}^{\top}\frac{1}{|\mathcal{Q}_{ij}|}\sum_{t\in\mathcal{Q}_{ij}}F_{t}e_{it}\xrightarrow{d}\mathcal{N}(0,\Gamma^{\textnormal{obs}}_{\Lambda,j})$ for every $j$.

We define the filtration $\mathcal{G}^{t}=\sigma(\cup_{s=1}^{T}\mathcal{G}^{t}_{Ts})$ with $\mathcal{G}^{t}_{Ts}=\sigma(\{W_{ij},j\leq s,\text{all }i\},\Lambda,v_{t})$ generated by $\{W_{ij},j\leq s,\text{all }i\}$, $\Lambda$ and $v_{t}$, which is given by $v_{t}=\Sigma_{\Lambda}^{-1}\Sigma_{F}^{-1}F_{t}$. For every $i$ and $t$, and $u_{i}=\Lambda_{i}$, it holds

where X_{i}=\frac{1}{N}\sum_{l=1}^{N}\Lambda_{l}\Lambda_{l}^{\top}\Big{(}\frac{1}{|\mathcal{Q}_{li}|}\sum_{s\in\mathcal{Q}_{li}}F_{s}F_{s}^{\top}-\frac{1}{T}\sum_{s=1}^{T}F_{s}F_{s}^{\top}\Big{)} and $\mathbf{X}_{t}=\frac{1}{N}\sum_{i=1}^{N}\frac{W_{it}}{P(W_{it}=1|S_{i})}X_{i}\Lambda_{i}\Lambda_{i}^{\top}$.

$\frac{1}{\sqrt{T_{1,i}}}\sum_{T_{0,i}+1}^{T}F_{t}e_{it}\xrightarrow{d}\mathcal{N}(0,\Sigma_{F,e_{i}})$.

Assumption GC3.5 holds for $v_{t}$ equal to $\Sigma_{\Lambda,t}^{-1}\Lambda^{(1)}_{i}$ and $\Sigma_{\Lambda,t}^{-1}(\Lambda^{(1)}_{i}-\Lambda^{(0)}_{i})$ under the filtration $\mathcal{G}=\sigma(\cup_{s=1}^{T}\mathcal{G}_{Ts})$ with $\mathcal{G}_{Ts}=\sigma(\{W_{ij},j\leq s,\text{all }i\},\Lambda)$ generated by $\{W_{ij},j\leq s,\text{all }i\}$ and $\Lambda$.

Assumption G2 describes an approximate factor structure and is at a similar level of generality as Bai (2003): (1) Assumption G2.1 ensures that each factor has a nontrivial contribution to the variation in $X$. (2) We assume loadings are random but independent of factors and errors in Assumption G2.2. We could study a factor model conditioned on some particular realization of the loadings, and the analysis would essentially be equivalent to that under the assumption that loadings are nonrandom. (3) Assumption G2.3 allows errors to be time-series and cross-sectionally weakly correlated. (4) Assumption G2.4 allows factors and idiosyncratic errors to be weakly correlated. (5) Assumption G2.5 guarantees that each loading and factor can be uniquely identified up to some rotation matrix. Additionally, we assume that these aspects also hold if we look at a subset of all time periods (the subset is denoted as $\mathcal{Q}_{ij}$ in Assumption G2). Together with Assumption C1.2, our covariance matrix estimator (1) using incomplete observations has similar properties as the conventional covariance matrix estimator $\frac{1}{T}XX^{\top}$ using full observations. For example, both $\frac{1}{|\mathcal{Q}_{ij}|}\sum_{t\in\mathcal{Q}_{ij}}X_{it}X_{jt}$ and $\frac{1}{T}\sum_{t=1}^{T}X_{it}X_{jt}$ are consistent estimators for $\Sigma_{ij}$. Moreover, the top $r$ eigenvalues estimated from both matrices are consistent as shown in Lemma 4 in the Internet Appendix, which is the foundation for developing the inferential theory of the factor model estimated from Equation (1).

is not sufficient as $X_{i}$ and $\mathbf{X}_{t}$ are multiplied with the random variables $u_{i}$ and $v_{t}$ in $\sqrt{T}\begin{bmatrix}(X_{i}u_{i})^{\top}&(\mathbf{X}_{t}^{S}v_{t})^{\top}\end{bmatrix}$. The asymptotic variances of these products are quadratic functions in the elements of those random variables given by $h_{i}(u_{i})$ and $g_{t}(v_{t})$ and take the form of $h_{i}(u_{i})=(u_{i}^{\top}\otimes I_{r})\Phi_{i}(u_{i}\otimes I_{r})$ and $g_{t}(v_{t})=(v_{t}^{\top}\otimes I_{r})\mathbf{\Phi}_{t}(v_{t}\otimes I_{r})$ respectively.

Assumption G3.5 requires a central limit theorem for stable convergence in law which is stronger than the conventional central limit theorem for convergence in distribution. The reason is that the asymptotic variance in Assumption G3.5 depends on both $\Lambda_{i}$ and $F_{t}$, which are random variables. Hence, we deal with a mixed normal limit and stable convergence in law ensures that the normal distribution of the central limit theorem will be independent of $\Lambda_{i}$ and $F_{t}$. Because of the stable convergence in law result, the estimated factors and common components normalized by their random standard deviation will converge to a standard normal distribution. In more detail, Assumption G3.5 implies that $X_{i}$ and $\mathbf{X}_{t}$ jointly converge $\mathcal{G}^{t}$-stably for $(N,T)\rightarrow\infty$ to a mixed normal distribution, whose asymptotic variance is random but measurable with respect to the sigma-field $\mathcal{G}^{t}$. Assumption G3.5 is used in Theorem 2 to show the asymptotic distribution of the variance correction term whose asymptotic variance is random. Our simplified factor model specified by Assumption S2 is sufficient to guarantee a central limit theorem for stable convergence in law. Proposition 3 shows that the simplified model implies Assumption G3.5.

Assumptions GC2 and GC3 are the corresponding assumptions for the propensity-weighted estimator with a similar level of generality. The additional Assumptions G4 and GC4 are only needed for the treatment effect tests. The simplified assumptions imply the general assumptions as stated in Proposition 3.

The simplified model is a special case of the general model:

Assumptions G2 and G3 are satisfied in the simplified model:

Assumptions S1 and S2 imply Assumption G2.

Assumptions S1, S2 and S3 imply Assumption G3.

Assumptions GC2 and GC3 are satisfied in the simplified conditional model:

Assumptions S1, C1, S2, S3.2, C2 and C3 imply Assumption GC3.

Assumptions G4 and GC4 are satisfied in the simplified model. Specifically,

Assumptions S1, S2 and S3 imply Assumption G4.

Assumptions S1, C1, S2, S3.2, C2 and C3 imply Assumption GC4.

References

Simulation Results

2 Asymptotic Distribution

3 Statistical Power of Treatment Effect Tests

4 Estimation under Misspecification

Proofs

In this note, we assume every entry is randomly missing with observed probability $q$. For a direct comparison, we follow the order of estimation in Jin, Miao, and Su (2021) and switch the role of factors and loadings in our all-purpose estimator: We first estimate the factors from the time-series sample covariance matrix, and then estimate the loadings from a time-series regression of the observed outcomes on the estimated factors.

Let the time-series sample covariance matrix be $\Sigma^{\dagger}$, where

Our approach to estimate the factors differs from Jin, Miao, and Su (2021) in that we adjust each entry in the sample covariance by the number of units that are observed in both time periods (i.e., the denominator $\mathcal{O}_{st}$ in Equation 18), while Jin, Miao, and Su (2021) adjust for the overall observed proportion (i.e., the denominator in Equation 18 is replaced by $N\hat{q}$ where $\hat{q}=\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}W_{it}$).

The corresponding term in the initial estimator in Jin, Miao, and Su (2021) is \big{(}\frac{1}{T}\sum_{t=1}^{T}F_{t}F_{t}^{\top}\big{)}\cdot\big{(}\frac{1}{\sqrt{N}q}\sum_{i=1}^{N}W_{is}\Lambda_{i}e_{is}\big{)}.

When every entry is missing at random, we have $\frac{\sqrt{N}}{|\mathcal{O}_{st}|}=\frac{1}{\sqrt{N}\cdot(|\mathcal{O}_{st}|/N)}\rightarrow\frac{1}{N\sqrt{q}}$. Jin, Miao, and Su (2021) differ from our term in that Jin, Miao, and Su (2021) disentangle the average over the time dimension from the average over the unit dimension. Under the simplified factor model, we can disentangle two averages and show that our term is asymptotically equivalent to the corresponding term in Jin, Miao, and Su (2021).

Under the simplified factor model, and using a similar proof as the third step in the proof of Proposition 3.2 on p23 in the online appendix, we can show that this term is asymptotically normal with the asymptotic variance

where $\omega_{tt}=\frac{1}{q}$. Therefore, the asymptotic variance coincides with the asymptotic variance of the initial estimator of factors in Jin, Miao, and Su (2021).

The corresponding term in the initial estimator in Jin, Miao, and Su (2021) is $\left(\frac{1}{T}\sum_{s=1}^{T}F_{s}F_{s}^{\top}\right)\cdot\left(\frac{1}{\sqrt{N}q}\sum_{i=1}^{N}W_{it}\Lambda_{i}\Lambda_{i}^{\top}-\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\Lambda_{i}\Lambda_{i}^{\top}\right)\cdot F_{t}$.

Under the simplified factor model, we can disentangle the average over the time dimension from the average over the unit dimension, and show that our term is asymptotically equivalent to the corresponding term in Jin, Miao, and Su (2021). Using a similar proof as the fourth step in the proof of Proposition 3.1(b), we can show that this term has a stable limiting distribution with the asymptotic variance

and $\omega_{tt}-1=\frac{1-q}{q}$ when every entry is missing at random. Therefore, the asymptotic variance coincides with the asymptotic variance of the initial estimator in Jin, Miao, and Su (2021).

Since both terms have the same asymptotic variance as those in the initial estimator in Jin, Miao, and Su (2021), our estimated factors are asymptotically the same as the initial estimates of factors in Jin, Miao, and Su (2021). Since the iterated estimator is more efficient than the initial estimator in Jin, Miao, and Su (2021), our estimated factors are asymptotically less efficient than the iterated estimates of factors in Jin, Miao, and Su (2021).

1.2 Proof of Proposition 1.2

where $\mathbf{X}_{i}=\frac{1}{T^{2}}\sum_{s=1}^{T}\sum_{t=1}^{T}F_{s}F_{s}^{\top}\left(\frac{1}{|\mathcal{O}_{st}|}\sum_{i=1}^{N}W_{it}\Lambda_{i}\Lambda_{i}^{\top}-\frac{1}{N}\sum_{i=1}^{N}\Lambda_{i}\Lambda_{i}^{\top}\right)W_{it}F_{t}F_{t}^{\top}$.

where \frac{1}{T}\sum_{t=1}^{T}W_{it}F_{t}F_{t}^{\top}=q\Sigma_{F}+O_{P}\big{(}\frac{1}{\sqrt{T}}\big{)} holds under the assumption that the observation pattern is exogenous and does not depend on the value of $F_{t}$. Therefore this term is asymptotically the same as the corresponding term in the initial estimator in Jin, Miao, and Su (2021).

The corresponding term in the initial estimator in Jin, Miao, and Su (2021) is $\frac{H^{\top}}{\sqrt{T}q}\cdot\sum_{t=1}^{T}W_{it}F_{t}F_{t}^{\top}\Lambda_{i}(W_{it}-1)$.

2 Proof of Proposition 3: Simplified Model

In this section, we prove Proposition 3.1, 3.2 and 3.3 step by step.

In this proof, we suppose Assumption S1 holds without further statement. We show Assumption G2.1 holds under Assumption S2.1, G2.2 holds under S2.2, G2.3 holds under S2.3, and G2.4 holds under S2.4.

Step 1.1: Show that $\frac{1}{T}\sum_{t=1}^{T}F_{t}F_{t}^{\top}\xrightarrow{P}\Sigma_{F}$ holds under Assumption S2.1.

Let $F_{t,j}$ be the $j$-th entry of $F_{t}$ and $\Sigma_{F,jk}$ be the $(j,k)$-th entry of $\Sigma_{F}$.

2.2 Proof of Proposition 3.1(b)

In this proof suppose Assumptions S1 and S2 hold without further statement. We show that each part in Assumption G3 holds under Assumption S3. For notation simplicity, denote $q_{ij}=q_{ij,ij}$.

Since $F$, $\Lambda$ and $e$ are independent, and $W$ is independent of $F$ and $e$, then $\phi_{it}$ is independent of $F$ and $e$ for any $i$ and $t$ and we have

Step 3.2: Show that $\frac{\sqrt{T}}{N}\sum_{i=1}^{N}\Lambda_{i}\Lambda_{i}^{\top}\frac{1}{|\mathcal{Q}_{ij}|}\sum_{t\in\mathcal{Q}_{ij}}F_{t}e_{jt}$ is asymptotically normal under Assumption S3.1. Assumption S3.1 and Assumption S2 together with the CLT imply

The CLT implies that $\begin{bmatrix}\frac{\sqrt{T}}{|\mathcal{Q}_{1j}|}\sum_{t\in\mathcal{Q}_{1j}}F_{t}^{\top}e_{jt}&\frac{\sqrt{T}}{|\mathcal{Q}_{2j}|}\sum_{t\in\mathcal{Q}_{2j}}F_{t}e_{jt}&\cdots&\frac{\sqrt{T}}{|\mathcal{Q}_{Nj}|}\sum_{t\in\mathcal{Q}_{Nj}}F_{t}e_{jt}\end{bmatrix}^{\top}$ is jointly asymptotic normal, and for $i\neq l$,

which follows from the results in Step 3.1.

$\Lambda_{i}e_{it}$ is independent across $i$. The CLT and the independence of $\Lambda$ and $e$ yield

as the last term is 0. Hence, by Chebyshev’s inequality

where the last equality follows from Assumption S3.2.

The vectorized form of $\Lambda_{i}\Lambda_{i}^{\top}\left(\frac{1}{|\mathcal{Q}_{ij}|}\sum_{t\in\mathcal{Q}_{ij}}F_{t}F_{t}^{\top}-\frac{1}{T}\sum_{t=1}^{T}F_{t}F_{t}^{\top}\right)$ is $(I_{r}\otimes\Lambda_{i}\Lambda_{i}^{\top})v^{(i,j)}$. Then Assumption G3.5 holds with

We can show the stable convergence in law similar to Jin, Miao, and Su (2021). Note that

where $\hat{q}_{ij}=\frac{|\mathcal{Q}_{ij}|}{N}$. Define the sigma-field $\mathcal{G}_{Tt}=\sigma(\{W_{is},s\leq t,\text{all }i\},\Lambda)$ that is generated from $\{W_{is},s\leq t,\text{all }i\}$ and $\Lambda$. Let $\mathcal{G}=\sigma(\cup_{t=1}^{T}\mathcal{G}_{Tt})$. Let $\omega\in R^{r}$ be a nonrandom vector with $\left\lVert\omega\right\rVert=1$. Let

where the last equality follows from Step 5.3.

We can rewrite $\frac{\sqrt{T}}{N^{2}}\sum_{l=1}^{N}\sum_{i=1}^{N}\Lambda_{l}\Lambda_{l}^{\top}\left(\frac{1}{|\mathcal{Q}_{li}|}\sum_{t\in\mathcal{Q}_{li}}F_{t}F_{t}^{\top}-\frac{1}{T}\sum_{t=1}^{T}F_{t}F_{t}^{\top}\right)W_{it}\Lambda_{i}\Lambda_{i}^{\top}$ as

where $\hat{q}_{il}=\frac{|\mathcal{Q}_{il}|}{N}$. Define the sigma-field $\mathcal{G}^{t}_{Ts}=\sigma(\{W_{iu},u\leq s,\text{all }i\},\Lambda,F_{t})$ that is generated from $\{W_{iu},u\leq s,\text{all }i\}$, $\Lambda$ and $F_{t}$. Let $\mathcal{G}^{t}=\sigma(\cup_{s=1}^{T}\mathcal{G}^{t}_{Ts})$. Let $\omega\in R^{r}$ be a nonrandom vector with $\left\lVert\omega\right\rVert=1$. Let

where $g_{t}(v_{t})=\mathbf{M}_{F,t}\mathbf{\Phi}_{t}\mathbf{M}_{F,t}^{\top}$. Step 5.5: Show that $\frac{\sqrt{T}}{N}\sum_{i=1}^{N}\Lambda_{i}\Lambda_{i}^{\top}\left(\frac{1}{|\mathcal{Q}_{ij}|}\sum_{t\in\mathcal{Q}_{ij}}F_{t}F_{t}^{\top}-\frac{1}{T}\sum_{t=1}^{T}F_{t}F_{t}^{\top}\right)$ and $\frac{\sqrt{T}}{N^{2}}\sum_{l=1}^{N}\sum_{i=1}^{N}\Lambda_{l}\Lambda_{l}^{\top}\left(\frac{1}{|\mathcal{Q}_{li}|}\sum_{t\in\mathcal{Q}_{li}}F_{t}F_{t}^{\top}-\frac{1}{T}\sum_{t=1}^{T}F_{t}F_{t}^{\top}\right)W_{it}\Lambda_{i}\Lambda_{i}^{\top}$ are jointly asymptotically normal and their asymptotic covariance converges. The randomness of these two terms come both from $v^{(l,i)}$, which is asymptotically normal. These two terms are weighted average of $v^{(l,i)}$ and therefore they are jointly asymptotically normal as well. Next, we show that their aymptotic covariance converges and we provide the limit.