Matrix Completion, Counterfactuals, and Factor Analysis of Missing Data

Introduction

Missing observations are prevalent in empirical work, and it is not surprising that solutions have been proposed by researchers in many disciplines. The classic econometric solution is some variant of the EM algorithm. In the case of factor analysis with missing data, the EM approach is to predict the missing values using initial estimates of the factors obtained from a balanced panel and iterate. While convergence of the algorithm can be established, the asymptotic properties of the converged estimates are not well understood.

Progress can be made if the panel of incompletely observed data $X$ is large in both dimensions and have a strong factor structure. This means in particular that $X$ has a common component $C$ of reduced rank $r$, and whose population covariance matrix has $r$ eigenvalues that increase with the size of the panel. We show in this paper that in spite of missing values in $X$, every entry of $C$ can be consistently estimated using a tall-wide (tw) algorithm that involves two applications of principal components. We provide an asymptotic characterization of the estimation error for each $C_{it}$ and show that there will be four convergence rates depending on observability of $X_{it}$. The approach can be used to construct missing values of potential outcomes satisfying a factor structure.

Our factor imputation approach contrasts with those used in matrix completions. In that literature, the probabilistic structure of the data is not specified; nuclear norm regularization via singular-value thresholding is crucial, and successful matrix recovery typically requires that the low rank component is incoherent, and that the data are missing uniformly at random. See, for example, Cai et al. (2008). Instead, we impose moment conditions to ensure that the assumed factor structure is strong and identifiable. We also require that $N_{o}\rightarrow\infty$ and $T_{o}\rightarrow\infty$, where $N_{o}$ is the number of units with data observed over the entire span, and $T_{o}$ is the length of the sample that data are available for all $N$ units. Exploiting the commonality between the observed and missing data allows recovery of the entire matrix $C$ using principal components with as few as $N_{o}T+T_{o}N$ observations. Provided that there are enough observed data to estimate the factors and the loadings consistently, a large fraction $1-\frac{N_{o}}{N}-\frac{T_{o}}{T}$ of the data can potentially be missing. And while regularization can yield fewer factors, it is not needed for consistent estimation of the missing values. In independent work, Jin et al. (2021) and Xiong and Pelger (2019) suggest alternative factor-based approaches. What makes our theory distinct is that we analyze factors estimated directly from the data without preliminary adjustment or iteration. The estimates that emerge from the tw algorithm are already consistent and asymptotically normal, though one re-estimation using imputed data will provide a faster convergence rate for one sub-block.

An immediate application of the tw approach is program evaluation in which the object of interest is the effect of treatment in the potential outcome model. The method of synthetic control pioneered in Abadie and Gardeazabal (2003) estimates potential outcomes from a weighted average of the control units, assuming that the difference between the treated and the control group is constant in the absence of treatment. Common variations between the treated and the control groups are crucial in estimation of counterfactual outcomes, and the factor model is a natural framework for modeling them. Treating potential outcomes as missing values, we develop a factor-based estimator of the individual as well as the average treatment effect. A distribution theory that permits tests of hypothesis is developed assuming that the size of the control group is large.

The rest of the paper is structured as follows. After presentation of the preliminaries in Section 2, Section 3 presents the least squares version of Algorithm tw and studies the asymptotic properties of the factor estimates that the algorithm delivers. Section 4 studies a factor-based estimation of treatment effect. Section 5 concludes.

Preliminaries

We use $i=1,\ldots N$ to index cross-section units and $t=1,\ldots T$ to index time series observations. Let $X_{i}=(X_{i1},\ldots X_{iT})^{\prime}$ be a $T\times 1$ vector of random variables and $X=(X_{1},X_{2},\ldots,X_{N})$ be a $T\times N$ matrix. In practice, $X_{i}$ is transformed to be stationary, demeaned, and often standardized. The normalized data $Z=\frac{X}{\sqrt{NT}}$ has singular value decomposition (svd)

In the above, $D_{r}$ is a diagonal matrix of $r$ singular values, $U_{r},V_{r}$ are the corresponding left and right singular vectors respectively. Without loss of generality, the singular values in the diagonal entries of $D_{r}$ are ordered such that $d_{1}\geq d_{2}\ldots\geq d_{r}$. Note that while the $r$ largest singular values of $X$ diverge and the remaining $N-r$ ones are bounded, the $r$ largest singular values of $Z$ are bounded and the remaining ones tend to zero because the singular values of $Z$ are those of $X$ divided by $\sqrt{NT}$. The Eckart and Young (1936) theorem posits that the best rank $k$ approximation of $Z$ is $U_{k}D_{k}V_{k}{{}^{\prime}}$. The svd is also the goto algorithm for solving matrix factorization problems that seek to represent a matrix $Z$ as a product of two low rank matrices. These results can be obtained without an assumed data generating process for $Z$.

We are interested in the principal components of $X$ viewed from the perspective of a factor model. Let $F$ be a $T\times r$ matrix of common factors, $\Lambda$ be a $N\times r$ matrix of factor loadings, and $e$ be a $T\times N$ matrix of idiosyncratic errors $e$. The data $X$ are assumed to have a factor structure

where $(F^{0},\Lambda^{0})$ are the true values of $(F,\Lambda)$. The common component $C^{0}=F^{0}\Lambda^{0{{}^{\prime}}}$ has reduced rank $r$ because $F^{0}$ and $\Lambda^{0}$ both have rank $r$. The econometrics literature on large dimensional factor analysis has largely adopted as framework the approximate factor model due to Chamberlain and Rothschild (1983) according to which the largest $r$ population eigenvalues of $X$ will increase with $N$ and $T$ while the remaining ones are bounded. The following assumptions formalize these ideas in a statistical setting:

There exists a constant $M<\infty$ not depending on $N,T$ such that

(Factors and Loadings): (i) $E\|F_{t}^{0}\|^{4}\leq M$, $\|\Lambda_{i}^{0}\|\leq M$, (ii) $\frac{F^{0{{}^{\prime}}}F^{0}}{T}\smash{\mathop{\longrightarrow}\limits^{p}}\Sigma_{F}>0$, and $\frac{\Lambda^{0{{}^{\prime}}}\Lambda^{0}}{N}\smash{\mathop{\longrightarrow}\limits^{p}}\Sigma_{\Lambda}>0$, and (iii) the eigenvalues of $\Sigma_{F}\Sigma_{\Lambda}$ are distinct.

$E(\frac{1}{N}\sum_{i=1}^{N}e_{it}e_{is})=\gamma_{N}(s,t)$, $\sum_{t=1}^{T}|\gamma_{N}(s,t)|\leq M$, $\forall s$;

$E(e_{it}e_{jt})=\tau_{ij,t}$, $|\tau_{ij,t}|\leq|\tau_{ij}|$ for some $\tau_{ij}$ $\forall t$, and $\sum_{j=1}^{N}|\tau_{ij}|\leq M$, $\forall i$;

$E(e_{it}e_{js})=\tau_{ij,st}$ and $\frac{1}{NT}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t=1}^{T}\sum_{s=1}^{T}|\tau_{ij,ts}|<M$;

$E|N^{-1/2}\sum_{i=1}^{N}[e_{is}e_{it}-E(e_{is}e_{it})]^{4}\leq M$ for every $(t,s)$;

$E(\frac{1}{N}\sum_{i=1}^{N}\|\frac{1}{\sqrt{T}}\sum_{t=1}^{T}F_{t}^{0}e_{it}\|^{2})\leq M$.

(Central Limit Theorems): for each $i$ and $t$, $\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\Lambda^{0}_{i}e_{it}\smash{\mathop{\longrightarrow}\limits^{d}}N(0,\Gamma_{t})$ as $N\rightarrow\infty$, and $\frac{1}{\sqrt{T}}\sum_{t=1}^{T}F^{0}_{t}e_{it}\smash{\mathop{\longrightarrow}\limits^{d}}N(0,\Phi_{i})$ as $T\rightarrow\infty$.

Assumption A ensures that the factor structure is strong and can be separated from the idiosyncratic errors. The requirements that $E\|F_{t}^{0}\|^{4}\leq M$ and $\|\Lambda_{i}^{0}\|\leq M$ ensure that the eigenvectors of the low rank (common) component are sufficiently spread out and play the role analogous to incoherence conditions. Positive definiteness of $\Sigma_{F}$ and $\Sigma_{\Lambda}$ ensure that each factor has a non-trivial contribution to the low rank component. It is possible for A(a.i) to hold but not A(a.ii) and vice versa, and identification will fail. The assumption of distinct eigenvalues is used to separately identify the factors and factor loadings but is not needed to identify the common components. Part (b) of Assumption A allows the errors to be weakly correlated both in the time and cross-section dimensions. For example, b(ii) requires the sum of autocovariances be bounded, and b(iii) is an analogous assumption for cross-sectional weak dependence. Assumption A(vi) assumes weak dependence between the factors and the errors. Part (c) is needed for asymptotic distribution of the factor estimates.

We observe $X$ but not $F^{0}$ or $\Lambda^{0}$, and $F$ and $\Lambda$ are not separately identifiable. We use as in Stock and Watson (2002); Bai and Ng (2002); Bai (2003) the normalizations $\frac{F{{}^{\prime}}F}{T}=I_{r}$ and $\Lambda^{\prime}\Lambda$ diagonal. The method of asymptotic principal components (APC) then constructs the factor estimates as

For each $t\in[1,T]$ and for each $i\in[1,N]$, $(\widetilde{F}_{t}$, $\widetilde{\Lambda}_{i})$ consistently estimate $(F_{t}^{0}$, $\Lambda^{0}_{i}$) up to rotation matrices $H$ and $G$ where

When Assumption A is satisfied, the factors and loadings are consistently estimable up to a rotation, and the low rank component is recoverable. Furthermore, the number of factors can be consistently estimated using, for example, the criteria developed in Bai and Ng (2002, 2019). Hence the number of factors $r$ can be treated as known.

Missing Data

Missing data is a problem that researchers frequently encounter. As Zhu et al. (2019) points out, we can expect more occurrence of incomplete observations in the era of big data. Data can be missing for a variety of reasons: non-response in surveys, lack of economic activity, and staggered releases by statistical agencies to name a few. One can always work with a balanced panel but this effectively throws away information in many series and cannot be efficient. This has led to development of simple methods that replace the missing values with zero or the mean as well as sophisticated methods that fully specify the data generating process and the missing data mechanism. For example, Rubin (1987) suggests a Bayesian approach that fills in missing values by repeatedly sampling from the predictive distribution of the missing values. Kamakura and Wedel (2000) suggests a simulation based approach that is aimed at handling different types of missing data in survey responses within a likelihood setup. See Horton and Kieinman (2007) for a survey of the literature. Rubin (1976) obtains two sufficient conditions for unbiased estimation. First, missingness cannot depend on the missing values after conditioning on the observed data (a condition known as missing at random), and second, the parameters of the model must not depend on the missingness mechanism. While missing at random can be a reasonable characterization in observational studies and surveys, there are situations when the assumption is not appropriate. For example, high income survey respondents may be more likely to ignore questions with tax consequences.

The EM algorithm of Dempster et al. (1977) imputes missing values by alternating between an E-step that computes the expected log-likelihood using the most recent parameter estimates, and an M-step that maximizes the expected log-likelihood. In cases when the expected log-likelihood is difficult to compute, the Expectation-Conditional Maximization algorithm of Meng and Rubin (1993) can be considered. Schneider (2001) considers a ridge-regression based regularized EM algorithm for imputing missing values in climate data. But as Honaker and King (2010) noted, methods that work well in a cross-section setting tend not to work well in panel data that exhibit dependence across units and over time. Imputing the missing values using the Kalman filter such as discussed in Shumway and Stoffer (1982) remains a popular time series approach, but it is fully parametric. Dropping a series altogether because of partially missing data could lead to a huge loss of information. In the FRED-MD database of over 130 series for example, as many as 30 series can be discarded even though some are missing for only a handful of months (about 2% of the observations in the panel) over the sample 1960:01 to 2018:12.

For estimation of strict factor models with missing data, more options are available. Banbura and Modugno (2014), Jungbacker and Koopman (2011), Jungbacker et al. (2011) consider likelihood estimation which is conceptually appealing but non-linear filters are needed to compute the likelihood as $F_{t}$ and $\Lambda_{i}$ are both random. For approximate factor models, Stock and Watson (1998) suggests to fill missing values in $X$ with the most recent estimate of the common component. Though the precise implementation may differ, using estimates from the balanced panel as initial values is the most common, see Stock and Watson (2016) and Giannone et al. (2008).

While a variety of methods have been proposed for factor analysis with missing data, there are surprisingly few theoretical results until recently. Jin et al. (2021) puts zeros to observations assumed to be missing at random and rescales the asymptotic principal components by the probability of missing data. It is shown that these estimates are consistent but not asymptotically normal in general, though iteration can restore normality. Xiong and Pelger (2019) estimates the factors from a weighted covariance matrix. The estimates are robust to the unknown missing pattern at the expense of larger variances. Our objective is similar to that of Jin et al. (2021) and Xiong and Pelger (2019) and also use a factor model for imputation, but we re-organize the data instead of re-weigh them, and we do not make assumptions about the missing data mechanism.

Suppose we rearrange the data such that the observed ones are ordered first. Rearrangement is not necessary in practice but it makes the idea easier to grasp. Consider the following example:

The transformation from $Z_{0}$ to $Z_{1}$ shuffles the columns so that those that are observed at all times are ordered first. The transformation from $Z_{1}$ to $Z_{2}$ shuffles the rows so that time periods with complete data for all units are ordered first.

We will use ‘o’ to denote the size of the observed and ‘m’ for the size of the missing samples, respectively. The northwest block, labeled bal, is a subpanel of complete data of dimension $T_{o}\times N_{o}$. The wide block extends the bal block in the cross-section dimension to include data of all $N$ units with data for $T_{o}$ periods. The tall block extends the bal block in the time dimension to include all $N_{o}$ units with complete time series observations. The southeast block collects the missing data into a $T_{m}\times N_{m}$ matrix where $T_{m}=T-T_{o}$ and $N_{m}=N-N_{o}$. This block, labeled miss, is the sub-block bordered by the rows $T_{o}+1:T$ and columns $N_{o}+1:N$. As drawn, miss is a “largest possible” block of missing data since some points in it are actually observed.

To give some economic content, we can think of the tall block as data for developed countries, the wide block for newly developed countries which have complete data over a shorter span, while the miss block consists of data for the less developed countries for which missing data are more prevalent. For financial data, acquisitions and mergers can yield a block structure. In macroeconomic settings, it is not uncommon for statistical agencies to stagger the release of data groups, but bunch the release of series within a group. In other cases, data may not collected in early years and terminated in later years due to attrition. As will be seen below, missing data also play a role in estimation of treatment effects, and as discussed in Athey et al. (2018), treatment may be given at the end of the sample, or it may be staggered or bunched by design of the experiment. In some of these cases, the assumption of missing at random may be inappropriate.

The estimation issue is that principal components cannot be directly applied when there are missing values. If we initialize using estimates from the balanced block, the factor estimate will always be spanned by factors that are originally in the balanced block. If this block is small in size, the information loss can be significant. Reorganizing data into four blocks reveals that we can make better use of the data to estimate the factors.

Our estimator is based on the idea that $(\widetilde{F}_{\text{tall}},\widetilde{\Lambda}_{\text{tall}})$ can be obtained from the tall block, while $(\widetilde{F}_{\text{wide}},\widetilde{\Lambda}_{\text{wide}})$ can be obtained from the wide block by APC, hence the acronym tall-wide, or tw for short. Results from our previous work can be used to show that

where $H_{\text{tall}}$ and $H_{\text{wide}}$ are unknown rotation matrices. The $o_{p}(1)$ terms are uniform in $i$ and $t$. To make further progress, some assumptions are needed. Our working assumptions are firstly, that there is an identifiable rank $r$ component in $X$ and in each of its four blocks so that there is commonality between blocks, and secondly that the observed blocks are sufficiently large so that consistent estimation by APC is possible.

(a) The conditions in Assumption A hold for the full $T\times N$ matrix $X$ if it were observed, as well as the four sub-blocks: bal, tall, wide, and miss. (b) The order conditions $TN_{o}>r(T+N_{o})$ and $T_{o}N>r(T_{o}+N)$ are satisfied for any $N,T,N_{o},T_{o}$. Furthermore, $\frac{\sqrt{N}}{\min\{N_{o},T_{o}\}}\rightarrow 0$ and $\frac{\sqrt{T}}{\min\{N_{o},T_{o}\}}\rightarrow 0$ as $N_{o},N\rightarrow\infty$ and $T_{o},T\rightarrow\infty$.

Assumption C:

(strong sub-block factor and factor loadings)

$\frac{\Lambda_{o}^{0\prime}\Lambda_{o}^{0}}{N_{o}}\smash{\mathop{\longrightarrow}\limits^{p}}\Sigma_{\Lambda,o}>0,\frac{\Lambda_{m}^{0\prime}\Lambda_{m}^{0}}{N_{m}}\smash{\mathop{\longrightarrow}\limits^{p}}\Sigma_{\Lambda,m}>0,\frac{1}{\sqrt{N_{o}}}\sum_{i=1}^{N_{o}}\Lambda_{i}^{0}e_{it}\smash{\mathop{\longrightarrow}\limits^{d}}N(0,\Gamma_{ot}),$

$\frac{F_{o}^{0\prime}F_{o}^{0}}{T_{o}}\smash{\mathop{\longrightarrow}\limits^{p}}\Sigma_{F,o}>0,\frac{F_{m}^{0\prime}F_{m}^{0}}{T_{m}}\smash{\mathop{\longrightarrow}\limits^{p}}\Sigma_{F,m}>0,\frac{1}{\sqrt{T_{o}}}\sum_{s=1}^{T_{o}}F_{s}^{0}e_{is}\smash{\mathop{\longrightarrow}\limits^{d}}N(0,\Phi_{oi}).$

Assumption D:

(block stationarity) Let $\Sigma_{\Lambda}$, $\Sigma_{F}$, $\Gamma_{t}$ and $\Phi_{i}$ be defined as in Assumption A.

$\Sigma_{\Lambda,o}=\Sigma_{\Lambda,m}=\Sigma_{\Lambda}$, and $\Gamma_{ot}=\Gamma_{t}$,

$\Sigma_{F,o}=\Sigma_{F,m}=\Sigma_{F}$, and $\Phi_{oi}=\Phi_{i}$.

Assumption B allows $p_{N}=N_{o}/N\rightarrow 0$ and $p_{T}=T_{o}/T\rightarrow 0$ as $N_{o},T_{o}\rightarrow\infty$ but requires an order condition to hold for the tall block, and one for the wide block. Assumption C requires the subsample moment matrices for the factor and factor loadings to be positive definite so that the factors from the subsamples are identifiable. Assumption D restricts these subsample limits to be the same as the limits for the whole sample. Results without Assumption D will also be stated.

Though there is no explicit restriction on the missing pattern, it is implicitly imposed through the factor model. If the miss block has nothing in common with the tall block, the missing values cannot be predicted. For example, if $N_{o}$ units are driven by a factor $F_{1}$ and the remaining $N-N_{o}$ units are driven by a different factor $F_{2}$, our method will not help with the imputation. Our method is also not suited for situations when $N_{o}$ or $T_{o}$ is small. We certainly cannot handle the case of $T_{o}=0$ or $N_{o}=0$. For example, if each $(i,t)$ entry is randomly observed with probability $p\leq\frac{1}{2}$ so that when $N\rightarrow\infty$, $(1-p^{N})^{N}\rightarrow 1$, $N_{o}$ will tend to zero and tall block will not be available for estimation.

Suppose that there are $N_{o}$ units with complete data in all $T$ rows, and $T_{o}$ periods when data are available for all $N$ units. Let $(\widetilde{F}_{\text{tall}},\widetilde{\Lambda}_{\text{tall}})$ and $(\widetilde{F}_{\text{wide}},\widetilde{\Lambda}_{\text{wide}})$ be obtained by applying the APC estimator to the tall and wide blocks respectively. Let $G_{c}=H_{c}^{-1}$ for $c=\text{tall},\text{wide}$. Suppose that $N_{o}\rightarrow\infty$ as $N\rightarrow\infty$ and $T_{o}\rightarrow\infty$ as $T\rightarrow\infty$. Under Assumptions A-D,

Lemma 2 is a direct implication of Lemma 1. The asymptotic variances are the same as defined in Lemma 1; the convergence rates differ because estimation is no longer based on the full sample.

Under our maintained assumptions, it immediately follows from Lemma 2 that

Obtaining an estimate of $C_{it}$ for $(i,t)$ in miss requires a bit more work because $\widetilde{F}_{\text{tall}}$ and $\widetilde{\Lambda}_{\text{wide}}$ are estimated from different blocks of data. However, if Assumptions A-D hold, the data in the block bal will contain information shared by both the tall and wide and can be used to re-rotate the data. Specifically, for any $i\in\textsc{bal}$, it holds that $\Lambda_{i}^{0}=H_{\text{wide}}\widetilde{\Lambda}_{\text{wide},i}+o_{p}(1)=H_{\text{tall}}\widetilde{\Lambda}_{\text{tall},i}+o_{p}(1)$, and

Define a new $r\times r$ non-singular rotation matrix

This matrix can be estimated by regressing the $N_{o}\times r$ matrix $\widetilde{\Lambda}_{\text{tall}}$ on the $N_{o}\times r$ matrix $\widetilde{\Lambda}_{\text{wide}}$, which is the $N_{o}\times r$ sub-matrix of $\widetilde{\Lambda}_{\text{wide}}$ associated with the balanced block. As shown in the appendix, $H_{\text{tall}}\widetilde{H}_{\text{miss}}H_{\text{wide}}^{-1}=I_{r}+O_{p}(1/N_{o}+1/T_{o})$. This suggests the following;

Algorithm TW

Let $\Omega$ be the $T\times N$ matrix that is one in positions when the data are observed, i.e. $\Omega_{it}=1$ if $X_{it}$ is observed and zero otherwise.

From the tall block of $X$, obtain $(\widetilde{F}_{\text{tall}},\widetilde{\Lambda}_{\text{tall}})$ by APC where $\widetilde{F}_{\text{tall}}$ is $T\times r$.

From the wide block of $X$, obtain $(\widetilde{F}_{\text{wide}},\widetilde{\Lambda}_{\text{wide}})$ by APC where $\widetilde{\Lambda}_{\text{wide}}$ is $N\times r$.

Let $\widetilde{C}_{\text{miss}}=\widetilde{F}_{\text{tall}}\widetilde{H}_{\text{miss}}\widetilde{\Lambda}_{\text{wide}}^{\prime}$ where $\widetilde{H}_{\text{miss}}$ is obtained by regressing $\widetilde{\Lambda}_{\text{tall}}$ on a submatrix of $\widetilde{\Lambda}_{\text{wide}}$.

Output $\widetilde{X}_{it}=X_{it}$ if $\Omega_{it}=1$ and $X_{it}=\widetilde{C}_{it}$ if $\Omega_{it}=0.$

Steps (1) and (2) imply that a complete set of estimates of the low rank component can be obtained from $TN_{o}+T_{o}N>T_{o}\times N_{o}$ data points. When the number of factors in tall and wide do not coincide, we let $r=\max(r_{\text{tall}},r_{\text{wide}})$ in Step (3).

Properties of C~~𝐶\widetilde{C} and Re-estimation

This section studies the properties of $\widetilde{C}_{it}$ which will be used to replace the missing $X_{it}$. The estimation error can be decomposed into four terms:

Proposition 1 shows that the estimates of the entire $C_{it}$ matrix are consistent and asymptotically normal without explicit restrictions on the nature of missingness. However, the convergence rate of $\widetilde{C}_{it}$ depends on whether $X_{it}$ is observed. Note that the bal block can use estimates from tall or from wide block, so convergence rate for this block is the faster of the two rates, ie.

In contrast, the convergence rate of $C_{it}$ in the miss block is always the slowest possible. These convergence rates and asymptotic distributions are obtained without iteration. It is possible for other estimators to achieve a better convergence rate. But to our knowledge, few (if any) consistent estimator exists that does not require iteration. Algorithm tw is best suited for cases when a large tall blow is availble, and the number of missing values is similar across units. In the event that a few units have many more missing values, $T_{o}$ will be the smallest sample size possible and could result in information loss. In Cahan et al. (2021), we propose an alternative algorithm that estimates the loadings and the rotation matrix jointly by a series of projections.

While Algorithm tw produces factor estimates that are mutually orthogonal within the four blocks of data, they are not mutually orthogonal over the entire data matrix. Furthermore, the estimates in tall do not use all information available, and similarly for the estimates in wide. Re-estimation using $\widetilde{X}$ provides an opportunity to use the imputed entries not previously available. However, embedded in $\widetilde{X}$ are imputation errors which will propagate to other blocks in re-estimation because the APC is a weighted average of $\widetilde{X}$. A formal analysis is needed to determine whether re-estimation using $\widetilde{X}$ can be justified.

Since $\widetilde{X}$ was constructed using estimates of $F$ constructed from the tall block, and of $\Lambda$ constructed from the wide block, we partition the matrices as follows:

where $T=T_{o}+T_{m}$ and $N=N_{o}+N_{m}$. We show in the Appendix that if $X_{it}$ is missing,

where $r_{it}=O_{p}(\delta^{-2}_{T_{o},N_{o}})$ uniformly in $(i,t)$ and

The dependence of these elements on $(T_{o},N_{o})$ is suppressed to simplify notation. Now

Imputation injects three errors into $\widetilde{X}_{it}$ when $X_{it}$ is not observed:- a quantity $r_{it}$ that is negligible, an error from estimating $F_{t}$, and one from estimating $\Lambda_{i}$. As a consequence, $u_{it}+v_{it}+r_{it}$ will differ from the true error $e_{it}$.

For $\widetilde{X}=\widetilde{U}\widetilde{D}\widetilde{V}{{}^{\prime}}$, let $(\widetilde{F}^{+},\widetilde{\Lambda}^{+})=(\sqrt{T}\widetilde{U}_{r},\sqrt{N}\widetilde{V}_{r}\widetilde{D}_{r})$ be the APC estimates based on $\widetilde{X}$ with the normalization $\frac{\widetilde{F}^{+\prime}\widetilde{F}^{+}}{T}=I_{r}$. Let $H^{+}=(\Lambda^{0\prime}\Lambda^{0}/N)(F^{0\prime}\widetilde{F}^{+}/T)\widetilde{D}_{r}^{-2}$. Then, under Assumptions A and B,

$\frac{1}{T}\sum_{t=1}^{T}\|\widetilde{F}^{+}_{t}-H^{+\prime}F^{0}_{t}\|^{2}=O_{p}(\delta^{-2}_{N_{o},T_{o}}).$

Bai and Ng (2002) shows that in the complete data case, $\frac{1}{T}\sum_{t=1}^{T}\|\widetilde{F}_{t}-H^{\prime}F_{t}^{0}\|^{2}=O_{p}(\delta_{NT}^{-2})$. Lemma 3 says that when the factors are estimated from $\widetilde{X}$, the convergence rate depends on the size of the balanced panel $T_{o}$ and $N_{o}$, which is evidently slower than when all data are observed.

To obtain a distribution theory for the factor estimates, we also need the representation for $\widetilde{F}^{+}$ and $\widetilde{\Lambda}^{+}$. We show in the Appendix that

where $\hat{\xi}_{NT,t}=O_{p}(\delta^{-2}_{N_{o},T_{o}})$ uniformly in $t$. The first representation is for those estimates of $F_{t}$ when $t\leq T_{o}$. Except for the $\hat{\xi}_{NT,t}$ term that is asymptotically negligible, the representation is the same as the case when all data are observed. More interesting is the $t>T_{o}$ case when $\widetilde{X}_{it}$ has imputation error. In appendix, we show

where the $r\times r$ matrix ${\mathbf{B}_{\Lambda}}$ is defined as

Thus the convergence rate for $\widetilde{F}^{+}_{t}-H^{\prime}F_{t}^{0}$ is $\sqrt{N_{o}}$. Similar derivations show that

where $G^{+}=(H^{+})^{-1}$, $\hat{\eta}_{NT,i}=O_{p}(\delta_{N_{o},T_{o}}^{-2})$ uniformly in $i$ and

Thus the convergence rate for $\widetilde{\Lambda}^{+}_{i}-G^{+}\Lambda_{i}^{0}$ is $\sqrt{T}$ for $i\leq N_{o}$ and the rate is $\sqrt{T_{o}}$ for $i>N_{o}$. Note that under Assumption D,

Given the asymptotic representations, it is relatively easy to derive the limiting distributions.

Under Assumptions $A-D$, the following holds as $N\rightarrow\infty$ and $T\rightarrow\infty$, with $G^{+}=(H^{+})^{-1}$.

The main implication of the proposition is that if $X_{it}$ is in the balanced block, then $\widetilde{F}^{+}_{t}$ is $\sqrt{N}$ consistent while $\widetilde{\Lambda}^{+}_{i}$ is $\sqrt{T}$ consistent. These are faster rates than those stated in Proposition 1 without re-estimation. Note that there is only a single rotation matrix for the factor estimates (instead of one for tall and one for wide) which are mutually orthogonal. This is a consequence of the fact that the factors are now estimated from the entire $\widetilde{X}$ matrix instead of the sub-blocks.

The four convergence rates in Proposition 2 are derived under the assumption that no data from the miss block are available. Parts (a) and (c) already have the same rate of convergence as in complete data and are unaffected by this assumption. It can be shown that the convergence rates for parts (b) and (d) will be faster when some observations in the miss block are available. The assumption that no observations in the miss block is consistent with the usual counterfactual matrix $Y(0)$ in the program evaluation analysis to follow.

Let $\widetilde{C}^{+}_{it}=\widetilde{F}^{+}_{t}{{}^{\prime}}\widetilde{\Lambda}^{+}_{i}$ be the common component estimated from $\widetilde{X}$ in which missing values of $X$ are replaced by the tw estimates of the common component $C$. Under Assumptions $A-D$, it holds that as $N\rightarrow\infty$ and $T\rightarrow\infty$,

The highlight of Proposition 3 is that three of the convergence rates are the same as in Proposition 1, but if $X_{it}$ is in bal, the convergence rate is now $\min(N,T)$, the same as when $X$ were completely observed. This improvement is due to the simple fact that $\widetilde{C}_{it}$ is based on data in the tall and wide blocks only, while $\widetilde{C}^{+}_{it}$ also also exploits information in the miss block. Bai and Ng (2019) considers a robust principal components (rpc) estimator $(\hat{F},\hat{\Lambda})=(\widetilde{F}(D^{\gamma}_{r})^{1/2},\widetilde{\Lambda}(D^{\gamma}_{r})^{-1/2})$ where $D_{ii}^{\gamma}=(D_{i}i-\gamma)_{+}$, $\gamma>0$ is a regualrization parameter. If we replace the apc part of Algorithm tw by rpc, $\bar{C}_{it}$ will also have the four convergence rates, unaffected by regularization.

From the proof of Proposition 3 in Appendix, the asymptotic representation for $\widetilde{C}_{it}^{+}-C_{it}^{0}$ implies the following error average rate in Frobenius norm (denoted $\|\cdot\|$) for the four blocks:

For the block defined by $i\leq N_{o},T\leq T_{o}$: $\frac{\|\widetilde{C}_{1}^{+}-C_{1}^{0}\|}{\sqrt{N_{o}T_{o}}}=O_{p}(\frac{1}{\sqrt{N}})+O_{p}(\frac{1}{\sqrt{T}})+O_{p}(\delta_{N_{o},T_{o}}^{-2})$.

For the block defined by $i\leq N_{o},t>T_{o}$: $\frac{\|\widetilde{C}_{2}^{+}-C_{2}^{0}\|}{\sqrt{N_{o}T_{m}}}=O_{p}(\frac{1}{\sqrt{N_{o}}})+O_{p}(\frac{1}{\sqrt{T}})+O_{p}(\delta_{N_{o},T_{o}}^{-2})$.

For block defined by $i>N_{o},T\leq T_{o}$: $\frac{\|\widetilde{C}_{3}^{+}-C_{3}^{0}\|}{\sqrt{N_{m}T_{o}}}=O_{p}(\frac{1}{\sqrt{N}})+O_{p}(\frac{1}{\sqrt{T_{o}}})+O_{p}(\delta_{N_{o},T_{o}}^{-2})$.

For the block defined by $i>N_{o},t>T_{o}$: $\frac{\|\widetilde{C}_{4}^{+}-C_{4}^{0}\|}{\sqrt{N_{m}T_{m}}}=O_{p}(\frac{1}{\sqrt{N_{o}}})+O_{p}(\frac{1}{\sqrt{T_{o}}})+O_{p}(\delta_{N_{o},T_{o}}^{-2})$.

This in turn implies an average squared error for the entire common components matrix

where the weights are the proportions of block size:

The sum of the first four terms in the average squared error is $O_{p}(\frac{1}{N})+O_{p}(\frac{1}{T})+(1-p_{N})(1-p_{T})O_{p}(\frac{1}{N_{o}}+\frac{1}{T_{o}})$ where $p_{N}=N_{o}/N$ and $p_{T}=T_{o}/T$. We have the following.

The first term in the square bracket is present even in the complete data case. The second term is due entirely to missing data, and the magnitude depends on the fraction of missing data but does not depend on the missing data mechanism.

Remark 1: In the machine learning literature, matrix completion problems are typically solved by nuclear norm regularization, and algorithmic errors bounds are given for fixed $N$ and $T$. Athey et al. (2018) finds that for $\sigma$-sub-Gaussian data, the worse case bound for the average error $\widetilde{C}$ (analogous to Corollary 1) depends on the regularization parameter and the unspecified distribution that generates $\Omega$. Our analysis complements the algorithmic error with a theory for asymptotic inference and shows that regularization is not needed for consistent estimation of $C$. In fact, we are able to characterize the sampling error of each $\widetilde{C}_{it}$, not just the average over over $i$ and $t$, while allowing $NT>>N_{o}T_{o}$. This is made possible by more fully exploiting the factor structure.

Compared to Jin et al. (2021) and Xiong and Pelger (2019), the missing pattern plays a less important role in our analysis because we look at the problem from the perspective of re-arranged data. Our results also differ in that our first step estimate is already consistent and asymptotically normal. Re-estimation in our setup accelerates the convergence rate of a sub-block rather than restores asymptotic normality. Furthermore, we do not need to assume that the number of missing data points as a fraction of $T\cdot N$ is bounded away from zero. If the balanced block is of dimension $T_{o}\times N_{o}$, such an assumption would have implicitly required that $T_{o}$ and $T$ are of the comparable order, and likewise for $N_{o}$ and $N$.

Lemma 2: the convergence rates are the same, but

Proposition 1: the convergence rates are the same, but for part i, replace $(\Sigma_{\Lambda},\Gamma_{t})$ by $(\Sigma_{\Lambda,o},\Gamma_{ot})$ in $V_{it}$; for part ii, replace $(\Sigma_{F},\Phi_{i}$) by $(\Sigma_{F,o},\Phi_{oi})$ in $W_{it}$; and for part iv, replace $(\Sigma_{\Lambda},\Sigma_{F},\Gamma_{t},\Phi_{i}$) by $(\Sigma_{\Lambda,o},\Sigma_{F,o},\Gamma_{ot},\Phi_{oi}$).

Proposition 2: parts $a$ and $c$ remain the same, the other two parts will read

where $\Gamma_{ot}$ and $\Phi_{oi}$ are given in Assumption C, and

are the limits of $\mathbf{B}_{\Lambda}$ and $\mathbf{B}_{F}$, respectively, $p_{1}$ is the limit of $N_{o}/N$ and $p_{2}$ is the limit of $T_{o}/T$.

2 Finite Sample Properties

Simulations are used to compare the performance of tw with and without updating. For comparison, we also consider an iterative EM algorithm considered in Stock and Watson (2016), which will be denoted em. Using estimates from the balanced panel as initial values, the algorithm repeatedly regresses $X$ on $F$ and then $X$ on $\Lambda$ by ols till convergence. Note, however, that the converged factor estimates produced by em may not be mutually orthogonal.

Data are generated from $F\sim N(0,D_{r})$ and $\Lambda\sim N(0,D_{r})$ with $r=2$, the diagonal entries in $D_{r}$ are equally spaced between 1 and $1/r$, and $e_{it}\sim N(0,1)$. We report results for $N=T=200$ only as those for $(N,T)=(200,400)$ and $(N,T)=(400,200)$ are similar. For each replication, $\|\widetilde{C}-C^{0}\|$ is computed for the four blocks. Also reported are results (labeled complete) for the infeasible case when all data are observable.

Our theory is silent about how to compute principal components. In the case of complete data, a common practice is to first standardize the data. But with missing data, it is unclear whether this is still desirable. Hence, we consider three versions of the estimator: one applied to the standardized data $X$, one to the demeanend data, one to the raw data. These are labeled TW(0,1,2) in the tables reported. The mean and standard deviation used in the centering and normalization are computed using the observations available for each series.

Table 1 compares the mean error over 5000 replications, normalized by the size of the corresponding block. Not surprisingly, the error in estimating the low rank component is inflated by missing data. As seen from Table 1, re-estimation always reduces the error, and imputation of the raw data (ie. method (2)) always has smaller errors than re-estimation using standardized data (ie. method (0)). One possible explanation is that when there are few observations from which to estimate the sample means and standard deviations, noise could be injected into the factor estimates. The results for em also favor imputation of the raw data.

The Frobenius normed error strongly favors $\widetilde{C}(\widetilde{X})$, but this is based on averaging the error over all $T\times N$ estimates of $C$. Table 2 reports the root-mean-square-error for four chosen $(i,t)$ pairs, one in each of the four blocks. Evidently, the estimation error is largest if $X_{it}$ is in the miss block and smallest when $X_{it}$ is in the bal block. As in Table 1, the error is smallest when the factors are re-estimated using $\widetilde{X}$. Both results are consistent with the theory.

In results not reported, the squared correlation between $\widetilde{C}_{i}$ and $C_{t}$ averaged over $i$ is over 0.93 when all data are observed. For the four data points considered in Table 1, the squared correlations are 0.89, 0.85, 0.85, 0.81 using TW, and 0.92, 0.93, 0.90, and 0.86 upon re-estimation. Regardless of re-estimation, $\widetilde{C}_{it}$ is well approximated by the normal distribution.

Factor Based Estimation of the Treatment Effects

If $X$ is a panel of data on GDP growth, a macroeconomist may be interested in a counterfactual prediction for some $i\in[1,N]$ at some time $t^{*}>T$. The theory in Bai and Ng (2006) can then be used to construct prediction intervals. We now show that Algorithm tw can also be used to estimate microeconomic type counterfactuals.

Program evaluation is widely used in economic analysis. Let $\mathcal{T}$ denote the treated group and $\mathcal{C}$ be the control group that is never exposed to treatment. To conform with the notation in this literature, we now index the treatment group by 1 and the control group by 0. The group size are $N_{1}$ and $N_{0}$ respectively with $N=N_{1}+N_{0}$. Unit $i$ receives treatment in period $T_{0,i}+1$ and thus $T_{0,i}$ is number of pretreatment periods for unit $i$. In this paper, we assume that $T_{0i}=T_{0}$ for all $i$ and thus $T_{1}=T-T_{0}$ is also constant across $i$.

Let $Y_{it}(1)$ be the potential outcome if individual $i$ receives the treatment, and $Y_{it}(0)$ be the potential outcome of individual $i$ without treatment in period $t$. The individual treatment effect is

The sample average treatment effect on the treated at $t>T_{0}$ (often referred to as sample $\text{ATT}_{t}$) is:

A variety of methods have been proposed to estimate these quantities. As discussed in Athey et al. (2018), the unconfoundedness regressions literature tends to use a single-treated period to impute the missing potential outcomes in the last period from the control units with similar lagged outcomes. The synthetic control literature pioneered in Abadie and Gardeazabal (2003) uses a weighted average of the control units ( $\sum_{i\in\mathcal{C}}w_{i}Y_{it}$) as estimate of the counterfactual outcome of the treated unit in the absence of treatment. The method can be specialized to produce difference-in-difference estimates under some conditions, see also Doudchenko and Imbens (2016). A ‘parallel-trend’ condition is needed to ensure that the sample path of the weighted average is parallel to the path of the treated in the absence of treatment.

Hsiao et al. (2012) assumes that potential outcome has a factor structure and considers estimation when the sample size is too small for estimation of the common factors. Their two step least squares-based estimator can be understood as replacing the latent factors by the outcome of the control units. To motivate the use of synthetic control methods in comparative case study research, Abadie et al. (2010) assumes that the outcome variable is driven by common factors. Increasingly, the synthetic control approach is studied from a factor model perspective. Gobillon and Magnac (2016) shows that if the true model is a linear factor model, synthetic controls are equivalent to factor models if the factor loadings and exogenous covariates for the treated are in the support of these variables for the control group. These conditions are analogous to requiring that the observed and missing blocks are driven by some factors in common.

In terms of factor-based estimation of treatment effects, Xu (2017) directly estimates the factors by principal components when $N$ and $T$ are large but the theoretical analysis is incomplete. Li (2018) suggests a procedure to determine $r$ and provides some asymptotic results for the treatment effect of a single treated unit in the absence of exogenous covariates. Amjad et al. (2018) analyses the mean-squared error of a robust synthetic control procedure. Xiong and Pelger (2019) allows the probability of missing data to depend on observed variables, but requires the fraction of observed data to bounded away from zero.

We will provide a distribution theory for estimates of the individual and average treatment effect without assuming a missing data mechanism. Indeed, when potential outcome is assumed to have a factor structure, estimation of treatment effect is very much related to factor analysis with missing data in which the outcomes for the treated group had they not been treated are the missing values to be recovered. Let $x_{it}$ be a $K\times 1$ vector of observed covariates. Let $D_{it}$ be the treatment indicator for individual $i$ is treated in period $t$. Then

where $F_{t}$ is $r\times 1$ vector of latent common factors. This is precisely the interactive fixed effect model developed in Bai (2009) where $C_{it}=\Lambda_{i}^{\prime}F_{t}$ is interactive fixed effect. The covariates $x_{it}$ also helps control for known sources of missingness and allows for more general missing patterns. We only observe $Y_{it}(1)$ on the treated and thus need to impute the corresponding counterfactual outcome in the absence of treatment. In the terminology of the previous section, we now have:

(IFE): Interactive fixed effect estimation of $\beta$ using observations in the control group. Let $R$ be $T\times N$ matrix of residuals where $R_{it}=y_{it}-x_{it}^{\prime}\hat{\beta}$

(tw): Estimate $F$ from $R_{\text{tall}}$ and $\Lambda$ from the $R_{\text{wide}}$; compute $\widetilde{C}=\widetilde{F}_{\text{tall}}\widetilde{H}_{miss}\widetilde{\Lambda}_{\text{wide}}{{}^{\prime}}$ and $\widetilde{C}^{+}=\widetilde{F}^{+}\widetilde{\Lambda}^{+}{{}^{\prime}}.$

replace the $(i,t)$th entry in $Y_{\text{miss}}(0)$ by $\hat{Y}_{it}(0)=x_{it}^{\prime}\hat{\beta}+\widetilde{C}_{it}$, OR

replace the $(i,t)$th entry in $Y_{\text{miss}}(0)$ by $\hat{Y}_{it}(0)=x_{it}^{\prime}\hat{\beta}+\widetilde{C}_{it}^{+}.$

Compute the average treatment effect $\widehat{\theta}_{t}=\frac{1}{N_{1}}\sum_{i\in\mathcal{T}}\hat{\theta}_{it}$.

The analysis to follow assumes re-estimation of the factors and so takes $\hat{C}_{it}=\widetilde{C}_{it}^{+}$ from (3b).

1 The Average Treatment Effect

Under the assumed factor structure, we see that for $(i,t)\in$ miss

There are three errors in the counterfactual $\hat{Y}_{it}(0)$, one from estimation of $\beta$, one from estimation of interactive fixed effects, and an idiosyncratic noise $e_{it}$. Since $Y_{it}(1)-\hat{Y}_{it}(0)=\theta_{it}+x_{it}(\beta-\hat{\beta})+C_{it}-\hat{C}_{it}+e_{it}$, it follows that

As $\frac{1}{\sqrt{N_{1}}}\sum_{i\in\mathcal{T}}e_{it}=O_{p}(1)$, the convergence rate for $\widehat{\theta}_{t}-\theta_{t}$ is at most $\sqrt{N_{1}}$. Now $\beta$ is homogeneous across $i$ and $t$ by assumption, and from Bai (2009), $\hat{\beta}-\beta=O_{p}(\frac{1}{\sqrt{T_{0}N_{0}}})$. The first term is thus $O_{p}(1/\sqrt{T_{0}N_{0}})$ and is dominated. By (A.20) in Appendix (or equation (A.4) if Step (3a) is used)

where $\bar{\Lambda}_{\mathcal{T}}=\frac{1}{N_{1}}\sum_{i\in\mathcal{T}}\Lambda_{i}$ is the average of factor loadings in the treatment group. If $N_{1}$ is large, the first term on the right is $O_{p}(1/\sqrt{T_{0}N_{1}})$ which is also dominated. This leads to the asymptotic representation for the case when $N_{1}$ large:

The following shows that $\widehat{\theta}_{t}$ is asymptotically normal.

Suppose Assumptions A-D and those in Bai (2009) hold. Then as $N_{0},T_{0},N_{1}\rightarrow\infty$,

The estimation of $\theta_{t}$ by tw presented above can be generalized to allow $T_{0}$ to vary with $i$, as in Xu (2017). This approach was first considered in the unpublished dissertation of Cahan (2013), and which we analyze further in Cahan et al. (2021).

2 Treatment Effect on a Single Unit

Consider now the estimation of treatment effect on a single unit $j$ for some $j>N_{0}$. Then

As before, we can ignore the error $x_{jt}^{\prime}(\beta-\hat{\beta})$. Imputation error from factor estimation (that is, the second and third terms on the right) are $O(1/\sqrt{T_{0}})$ and $O(1/\sqrt{N_{0}})$, respectively. But unlike $\widehat{\theta}_{t}$, no averaging is taken over $i=N_{0}+1,\ldots,N$ and, as a consequence, $e_{jt}$ now dominates the composite estimation error. The distribution of $\widehat{\theta}_{jt}-\theta_{jt}$ thus depends on the distribution of $e_{jt}$. If one is willing to assume $e_{jt}$ is identically distributed across $i$ and $t$, its distribution can be estimated using the residuals $\hat{e}_{it}=Y_{it}(0)-\hat{Y}_{it}(0).$ The estimated individual treatment effect has variance

It is also of interest to consider the average treatment effect over the treatment period for a single unit, defined as $\theta_{j}=\frac{1}{T_{1}}\sum_{s>T_{0}}\theta_{js}$. Let $\hat{\theta}_{j}=\frac{1}{T_{1}}\sum_{s>T_{0}}\hat{\theta}_{js}$ and $\bar{F}=\frac{1}{T_{1}}\sum_{s>T_{0}}F_{s}$. It can be shown that a result similar to Proposition 4 holds:

Table 3 reports the bias, root-mean-squared-error, and the probability that the true treatment effect is within two standard errors of the estimated effect (labeled covr). The first panel evaluates $\hat{\theta}_{j,t}$ at $(j,t)=(N_{0}+1,T_{0}+5)$, the second panel reports $\widehat{\theta}_{t}$ at $t=T_{0}+5$. The estimation error is larger for $\hat{\theta}_{j,t}$ than $\hat{\theta}_{T_{0}}$ since averaging reduces noise. The error decreases with $N_{0}$ for given $N_{1}$, but is quite insensitive to changes in $N_{1}$ for given $N_{0}$ because the convergence rate is $\min(N_{0},N_{1})$. Except when $N_{0}$ and $T_{0}$ are both small, the coverage is quite close to 95%.

Conclusion

Missing data is prevalent in empirical work. There is a presumption that iteration is needed to impute missing values, and successful matrix recovery requires solving a regularized problem under a missing at random assumption. This paper shows that if we are willing to impose a strong factor structure, then the entire low rank component of the data can be consistently estimated by our proposed tw procedure without iteration or regularization. The methodology can be used within the potential outcomes framework to estimate the effect of treatment on the treated, and a distribution theory is provided.

References

Appendix A

This appendix provides proofs to the results in the main text along with more general results that are of independent interest. Throughout, $\|A\|$ denotes the Frobenius norm of matrix $A$.

Recall the notation $T=T_{o}+T_{m};N=N_{o}+N_{m}$ and

where $F_{o}^{0}$ is $T_{o}\times r$, and $F_{m}^{0}$ is $T_{m}\times r$; $\Lambda_{o}^{0}$ is $N_{0}\times r$ and $\Lambda_{m}^{0}$ is $N_{m}\times r$. The partitions of $\widetilde{F}_{tall}$ and $\widetilde{\Lambda}_{wide}$ are the same. We write $(i,t)\in\Omega$ if $\Omega_{it}=1$, and $(i,t)\in\Omega_{\bot}$ if $\Omega_{it}=0$.

Similar to the rotation matrix $H$ given in Section 2, let $H_{tall}^{\prime}=D_{tall}^{-2}(\widetilde{F}_{tall}^{\prime}F^{0}/T)(\Lambda_{o}^{\prime}\Lambda_{o}/N_{o})$, then the tall estimator satisfies (see Bai (2003), Theorem 1)

where $\xi_{NT,t}=O_{p}(1/N_{o}+1/T)$ uniformly in $t$. Similarly, there is a rotation matrix $H_{wide}$ such that for each $i$, the wide estimator satisfies (see Bai (2003), Theorem 2)

where $\eta_{NT,i}=O_{p}(1/T_{o}+1/N)$ uniformly in $i$. Let $H_{miss}=H_{tall}^{-1}H_{wide}$, and define for $(i,t)\in\Omega_{\bot}$

where $\widetilde{H}_{miss}$ is an estimator for $H_{miss}$ obtained by regressing $\widetilde{\Lambda}_{tall,i}$ on $\widetilde{\Lambda}_{wide,i}$ for $i=1,2,...,N_{o}$.

Under the assumptions of Proposition 1: (i)

where $r_{NT,it}=O_{p}(1/T_{o}+1/N_{o})$ uniformly in $(i,t)\in\Omega_{\bot}$.

Consider part (i). Rewrite the representation in (A.2) as

This follows from [see Bai (2003), p 166) and Bai and Ng (2019)],

Similarly, the Tall estimator has the asymptotic representation

This implies $\widetilde{\Lambda}_{tall,i}=H_{miss}\widetilde{\Lambda}_{wide,i}+\varphi_{NT,i},\quad i=1,2...,N_{o}$, where

Note $\frac{1}{N_{o}}\sum_{i=1}^{N_{o}}\widetilde{\Lambda}_{wide,i}\widetilde{\Lambda}_{wide,i}^{\prime}$ converges in probability to a positive definite matrix. Consider the numerator,

Replace $\widetilde{\Lambda}_{wide,i}^{\prime}$ by $\Lambda_{i}^{0\prime}H_{wide}^{\prime}$ (ignore higher orders), we see the two terms on the right hand side are each $O(1/\sqrt{N_{o}T_{o}})$, thus dominated by $O_{p}(1/\delta_{N_{o}T_{o}}^{2})$. This proves (A.3).

We next proof part (ii). We can rewrite the representations in (A.1) by

This follows by multiplying (A.1) by $H_{tall}^{-1\prime}$ and using

where $\widetilde{F}_{tall}=(\widetilde{F}_{tall,1},...,\widetilde{F}_{tall,T})^{\prime}$. The second equality uses the definition of $H_{tall}^{\prime}$. Rewrite (A.5) as

Multiply (A.6) by $\Lambda_{i}^{0\prime}$ and multiply (A.7) by $F^{0\prime}$, we have

where $u_{it}$ and $v_{it}$ are defined earlier. Thus

Proof of Lemma 2.

Proof of Proposition 1.

This proposition is also a direct consequence of Lemma 1. Applying Lemma 1 to the tall and wide blocks separately, we obtain the results in (i), (ii), and (iii). For the missing block, that is, for $(i,t)\in\Omega_{\bot}$, the limiting distribution for the estimated common component follows from the asymptotic representation in (A.4). By Assumption C,

Analysis based on imputed data matrix

Let $\widetilde{X}_{it}=X_{it},(i,t)\in\Omega$ and $\widetilde{X}_{it}=\widetilde{C}_{it},(i,t)\in\Omega_{\bot}$ so the missing values are replaced by the estimated common components $\widetilde{C}_{it}$. We have $X_{it}=\Lambda_{i}^{0\prime}F_{t}^{0}+e_{it},\quad(i,t)\in\Omega$ and $\widetilde{X}_{it}=\Lambda_{i}^{0\prime}F_{t}^{0}+u_{it}+v_{it}+r_{NT,it},\quad(i,t)\in\Omega_{\bot}.$ Consider estimating the factor and factor loadings using the $T\times N$ matrix $\widetilde{X}=(\widetilde{X}_{it})$. Let $\widetilde{F}^{+}$ be the first $r$ eigenvectors corresponding to the first $r$ largest eigenvalues (arranged in decreasing order) of the matrix $\widetilde{X}\widetilde{X}^{\prime}/(NT)$ with the normalization $\widetilde{F}^{+\prime}\widetilde{F}^{+}/T=I_{r}$, that is,