Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
83,914 characters · 28 sections · 48 citation commands
Matrix Completion Methods for Causal Panel Data Models
{\it Keywords:} Causality, Synthetic Controls, Unconfoundedness, Interactive Fixed Effects, Low-Rank Matrix Estimation
In this paper we develop new methods for estimating average causal effects in settings with panel or longitudinal data, where some units are exposed to a binary treatment during some periods. To estimate the average causal effect of the treatment on the treated units in this setting, we impute the missing potential control outcomes.
The statistics and econometrics causal inference literatures have taken two general approaches to this problem. The literature on unconfoundedness (rosenbaum1983central, imbens2015causal) can be interpreted as imputing missing potential control outcomes for treated units using observed control outcomes for control units with similar values for observed outcomes in previous periods. In contrast, the recent synthetic control literature (abadie2003, abadie2010synthetic, abadie2014, doudchenko, ben2018augmented, li2019statistical, ferman2019synthetic, arkhangelsky2019synthetic, chernozhukov2017exact, see abadie2019using for a review) imputes missing control outcomes for treated units using weighted average outcomes for control units with the weights chosen so that the weighted lagged control outcomes match the lagged outcomes for treated units. Although at first sight similar, the two approaches are conceptually quite different in terms of the correlation patterns in the data they exploit to impute the missing potential outcomes. The unconfoundedness approach assumes that patterns over time are stable across units, and the synthetic control approach assumes that patterns across units are stable over time. In empirical work the two sets of methods have primarily been applied in settings with different structures on the missing data or assignment mechanism. In the case of the unconfoundedness literature the typical setting is one with the treated units all treated in the same periods, typically only the last period, and with a substantial number of control and treated units. The synthetic control literature has primarily focused on the setting with one or a small number of treated units observed prior to the treatment over a substantial number of periods. We argue that once regularization methods are used, the two approaches, unconfoundedness and synthetic controls, are applicable in the same settings, leaving the researcher with a real choice in terms of methods. In addition this insight allows for a more systematic comparison of their perfomance than has been appreciated in the literature.
In this study we draw on the econometric literature on factor models and interactive fixed effects, and the computer science and statistics literatures on matrix completion, to take an approach to imputing the missing potential outcomes that is different from the unconfoundedness and synthetic control approaches. In fact, we show that it can be viewed as nesting both. In the literature on factor models and interactive effects (bai2002determining, bai2003inferential) researchers model the observed outcome as the sum of a linear function of covariates and an unobserved component that is a low rank matrix plus noise. Estimates are typically based on minimizing the sum of squared errors given the rank of the matrix of unobserved components, sometimes with the rank estimated. xu2017generalized extends these ideas to causal settings where a subset of units is treated from a common period onward, so that complete data methods for estimating the factors and factor loadings can be exploited. The matrix completion literature (candes2009exact, candes2010matrix, mazumder2010spectral) focuses on imputing missing elements in a matrix assuming that: $(i)$ the complete matrix is the sum of a low rank matrix plus noise and $(ii)$, the missingness is completely at random (except gamarnik2016note that study a stylized rank one case). The rank of the matrix is implicitly determined by the regularization through the addition of a penalty term to the objective function. Especially with complex missing data patterns using the nuclear norm as the regularizer is attractive for computational reasons.
In the current paper we make three contributions. First, we present formal results for settings where the missing data patterns are not completely at random and have a structure that allows for correlation over time, generalizing the results from the matrix completion literature. In particular we allow for the possibility of staggered adoption ({\it e.g.,} athey2018design, shaikh2019randomization), where units are treated from some initial adoption date onwards, but the adoption dates vary between units. We also modify the estimators from the matrix completion and factor model literatures to allow for unregularized unit and time fixed effects. Although these can be incorporated in the low rank matrix, in practice the performance of the estimator with the unregularized two-way fixed effects is substantially better. Compared to the factor model literature in econometrics the proposed estimator focuses on nuclear norm regularization to avoid the computational difficulties that would arise for complex missing data patterns with the fixed-rank methods in bai2002determining and xu2017generalized, similar to the way LASSO (or $\ell_1$ regularization, tibshirani1996regression) is computationally attractive relative to subset selection (or $\ell_0$ regularization) in linear regression models. The second contribution is to show that the synthetic control and unconfoundedness approaches, as well as our proposed method, can all be viewed as matrix completion methods based on matrix factorization, all with the same objective function based on the Fr\"obenius norm for the difference between the latent matrix and the observed matrix. Given this common objective function, the unconfoundedness and synthetic control approaches impose different sets of restrictions on the factors in the matrix factorization. In contrast, the proposed method does not impose any restrictions but uses regularization to characterize the estimator. In our third contribution we apply our methods to two real data sets where we observe the complete matrix. We artificially designate outcomes for some units and time periods to be missing, and then compare the performance of different imputation estimators. We find that the nuclear norm matrix completion estimator does well in a range of cases, including when $T$ is small relative to $N$, when $T$ is large relative to $N$, and when $T$ and $N$ are comparable. In contrast, the unconfoundedness and synthetic control approaches break down in some of these settings in the expected pattern (the unconfoundedness approach does not work very well if $T\gg N$, and the synthetic control approach does not work very well if $N\gg T$).
We discuss some extensions in the final part of the paper. In particular we consider extensions to settings where the probability of assignment to the treatment may vary systematically with observed characteristics. In the program evaluation literature such settings have often been addressed using inverse propensity score weighting (rubin2006matched, hirano2003efficient), which can be applied here as well.
We start by stating the causal problem. Consider a setting with $N$ units observed over $T$ periods. In each period each unit is characterized by two potential outcomes, $Y_{it}(0)$ and $Y_{it}(1)$. In period $t$ unit $i$ is exposed or not to a binary treatment, with $W_{it}=1$ indicating that the unit is exposed to the treatment and $W_{it}=0$ otherwise. We observe for each unit and period the pair $(W_{it},Y_{it})$ where the realized outcome is $Y_{it}=Y_{it}(W_{it})$. In addition to observing the matrix ${\bf Y}$ of realized outcomes and the matrix of treatment assignments ${\bf W}$, we may also observe covariate matrices ${\bf X}\in{\mathbb{R}}^{N\times P}$ and ${\bf Z}\in{\mathbb{R}}^{T\times Q}$ where columns of ${\bf X}$ are unit-specific covariates, and columns of ${\bf Z}$ are time-specific covariates. We may also observe unit/time specific covariates $V_{it}\in{\mathbb{R}}^J$. Implicit in our set up is that we rule out dynamic effects and make the stable-unit-treatment-value assumption (rubin2006matched, imbens2015causal): the potential outcomes are indexed only by the contemporaneous treatment for that unit and not by past treatments or treatments for other units. Cases where such assumptions are restrictive include those analyzed in the dynamic treatment regime literature (chamberlain1993feedback, hernan2010causal). In the case where units are only exposed to the treatment in the last period this issue is not material. Also, in the case with staggered adoption violations of the no-dynamics assumption simply changes the interpretation of the estimand, but does not in general invalidate a causal interpretation.
Here we focus on estimating the average effect for the treated, $\tau=\sum_{(i,t):W_{it}=1} [Y_{it}(1)-Y_{it}(0)]/\sum_{i,t} W_{it}$, although other averages such as the overall average causal effect, $\sum_{i,t} [Y_{it}(1)-Y_{it}(0)]/(NT)$, could be of interest too. In order to estimate such average treatment effects, one approach is to impute the missing potential outcomes. Because we focus on estimating the average effect for the treated, all the relevant values for $Y_{it}(1)$ are observed, and thus we only need to impute the missing entries in the ${\bf Y}(0)$ matrix for treated units with $W_{it}=1$. For the moment we focus on the problem of imputing the missing entries in ${\bf Y}(0)$ given the observed values of ${\bf Y}(0)$ and the observed matrix ${\bf W}$. To ease the notation and facilitate the connection to the matrix completion literature we drop from here on the $(0)$ part of ${\bf Y}(0)$ and simply focus on imputing the missing values of a partially observed matrix ${\bf Y}$ (with the understanding that this would be the matrix of control outcomes ${\bf Y}(0)$), with ${\bf W}$ the matrix of missing data (treatment assignment) indicators. One may also wish to use the observed values of ${\bf Y}(1)$ for imputing the missing values for ${\bf Y}(0)$, but we do not do so here. In setting with few values of ${\bf Y}(1)$ observed it is unlikely that the information in these values is important. (In particular in the case we focus on for part of this study, with only a single treated unit/period pair there would be no information in this value.). Extension to the cases that leverage also data from ${\bf Y}(1)$ require assumptions on the treatment effect and are briefly discussed in \S (ref).
For any positive integer $n$, we use notation $[n]$ to refer to the set $\{1,\ldots,n\}$ and use ${\mathbf{1}}_n$ to denote the $n$ by $1$ vector of all ones. We define ${\cal M}$ to be the set of pairs of indices $(i,t)$, $i\in[N]$, $t\in[T]$, corresponding to the missing entries with $W_{it}=1$ and ${\cal O}$ to be the set os pairs of indices corresponding to the observed entries in ${\bf Y}$ with $W_{it}=0$. Putting aside the covariates for the time being, the data can be thought of as consisting of two $N\times T$ matrices, one incomplete and one complete, {
} where \[W_{it}=\left\{
\right.\] is an indicator for the event that the corresponding component of ${\bf Y}$, that is, $Y_{it}$, is missing. The main part of the paper is about the statistical problem of imputing the missing values in ${\bf Y}$. Once these are imputed we can then estimate the average causal effect of interst, $\tau$.
In this section, we discuss a number of particular configurations of the matrices ${\bf Y}$ and ${\bf W}$ that are the focus of distinct parts of the general literature. This discussion serves to put in context the problem, and to motivate previously developed methods from the literature on causal inference under unconfoundedness, the synthetic control literature, and the interactive fixed effect literature, and subsequently to develop formal connections between all three and the matrix completion literature. Note that the matrix completion literature has focused primarily on the case where ${\bf W}$ is completely random, as in Equation ((ref)), and where both dimensions of ${\bf Y}$ and ${\bf W}$ are large. First, we consider patterns of missing data, that is, distributions for ${\bf W}$ that differ from completely random. Second, we consider different shapes of the matrix ${\bf Y}$ where the relative size of the dimensions $N$ and $T$ may be very different and one or both may be modest in magnitude. Third, we consider a number of specific analyses in the econometrics literature that focus on particular combinations of missing data patterns and shapes of the matrices.
In the statistics and computer science literatures on matrix completion the focus is typically on settings with randomly missing values, allowing for general patterns on the matrix of missing data indicators candes2010thepower,recht2011simpler. In contrast in causal social science applications the missingness arises from treatment assignments and the choices that lead to these assignments. As a result are often specific structures on the missing data that depart substantially from complete randomness.
A leading example is a block structure, with a subset of the units adopting an irreversible treatment at a particular point in time $T_0+1$. In the example below the $\checkmark$ marks indicate observed values and the $?$ indicate missing values. { \[ {\bf Y}_{N\times T}=\left(
\right)\,. \] } There are two special cases of the block structure. Much of the literature on estimating average treatment effects under unconfoundedness ({\it e.g.,} imbens2015causal) focuses on the case where $T_0=T-1$, so that the only treated units are in the last period. We will refer to this as the single-treated-period block structure. In contrast, the synthetic control literature ({\it e.g.}, abadie2010synthetic, abadie2019using) focuses primarily on the case of with a single treated unit which are treated for a number of periods from period $T_0+1$ onwards, the single-treated-unit block structure: { \[ {\bf Y}=\left(
\right)\hskip0.5cm{\rm and}\ \ {\bf Y}=\left(
\right)\,.\] } A special case that fits in both these settings is that with a single missing unit/time pair:
{ \[ {\bf Y}=\left(
\right).\] } This specific setting is useful to contrast methods developed for the single-treated period (unconfoundedness) case with those developed for the single-treated unit (synthetic control) case because both sets of methods are potentially applicable.
Another setting that has received attention is the staggered adoption design (athey2018design, shaikh2019randomization). Here units may differ in the time they are first exposed to the treatment, but the treatment is irreversible. This naturally arises in settings where the treatment is some new technology that units can choose to adopt (e.g., athey2002impact). Here: { \[ {\bf Y}_{N\times T}=\left(
\right)\,. \] }
A second classification of the problem concerns the shape of the matrix ${\bf Y}$. Relative to the number of time periods, we may have many units, few units, or a comparable number. These data configurations may make particular analyses more attractive partly by removing the need for regularization. For example, ${\bf Y}$ may be a thin matrix, with $N\gg T$, or a fat matrix, with $N\ll T$, or an approximately square matrix, with $N\approx T$: { \[ {\bf Y}=\left(
\right)\hskip0.3cm ({\bf thin}) \hskip1cm {\bf Y}=\left(
\right)\hskip0.2cm ({\bf fat}),\] } or { \[{\bf Y}=\left(
\right)\hskip0.2cm ({\bf approximately\ square}).\] }
Two special combinations of missing data patterns and matrix shape deserve particular attention because they are the focus of large, mostly separate, literatures.
The unconfoundedness literature (rosenbaum1983central, rubin2006matched, imbenswooldridge, abadie2018econometric) focuses primarily on the single-treated-period block structure with a thin matrix ($N\gg T$), a substantial number of treated and control units, and imputes the missing potential outcomes in the last period using control units with similar lagged outcomes: { \[ {\bf Y}=\left(
\right),\] } A simple version of the unconfoundedness approach is to regress the last period outcome on the lagged outcomes and use the estimated regression to predict the missing potential outcomes. That is, for the units with $(i,T)\in{\cal M}$, the predicted outcome is
We refer to this as a {\bf horizontal} regression, where the rows of the ${\bf Y}$ matrix form the units of observation. A more flexible, nonparametric, version of this estimator would correspond to matching where we find for each treated unit $i$ a corresponding control unit $j$ with $Y_{jt}$ approximately equal to $Y_{it}$ for all pre-treatment periods $t=1,\ldots,T-1$.
The synthetic control literature (abadie2010synthetic) focuses primarily on the single-treated-unit block structure with a relatively fat $(T\gg N$) or approximately square matrix $(T\approx N$), and a substantial number of pre-treatment periods: { \[ {\bf Y}=\left(
\right).\] } doudchenko and ferman2019synthetic show how the Abadie-Diamond-Hainmueller synthetic control method can be interpreted as regressing the outcomes for the treated unit prior to the treatment on the outcomes for the control units in the same periods. That is, for the treated unit in period $t$, for $t=T_0,\ldots,T$, the predicted outcome is
We refer to this as a {\bf vertical} regression, where the columns of the ${\bf Y}$ matrix form the units of observation. As shown in doudchenko, this is generalization of the original abadie2010synthetic synthetic control estimator, relaxing two restriction: $(i)$ that the coefficients are nonnegative and $(ii)$ that the intercept in this regression is zero. Note that these restrictions may well be substantively plausible and they can greatly improve precision.
Although this does not appear to have been pointed out previously, a matching version of this estimator would correspond to finding, for each period $t$ where unit $N$ is treated, a corresponding period $s\in\{1,\ldots,T_0\}$ such that $Y_{is}$ is approximately equal to $Y_{Ns}$ for all control units $i=1,\ldots,N-1$. This matching version of the synthetic control estimator may serve to clarify the link between the treatment effect literature under unconfoundedness and the synthetic control literature.
Suppose that the only missing entry is in the last period for unit $N$, period $T$. In that case if we estimate the horizontal regression in ((ref)), it is still the case that imputed $\hat Y_{NT}$ is linear in the observed $Y_{1T},\ldots,Y_{N-1,T}$, just with different weights than those obtained from the vertical regression. Similarly, if we estimate the vertical regression in ((ref)), it is still the case that $\hat Y_{NT}$ is linear in $Y_{N1},\ldots,Y_{N,T-1}$, just with different weights from the horizontal regression in ((ref)). Note also that the restrictions that the coefficients are nonnegative and sum to one are common in the synthetic control literature, but could also be imposed in the unconfoundedness literature, although they do not appear to have been used there.
Juxtaposing the unconfoundedness and synthetic control approaches as we have done here raises the question how they are related, and whether there is an approach that avoids the choice between focusing on the cross-section and time-series correlation patterns. We further elaborate on the connection between the horizontal and vertical regression in \S (ref) after introducing a third approach.
The horizontal regression focuses on a pattern in the time path of the outcome $Y_{it}$, specifically the relation between $Y_{iT}$ and the lagged $Y_{it}$ for $t=1,\ldots,T-1$, for the units for whom these values are observed, and assumes that this pattern is the same for units with missing outcomes. The vertical regression focuses on a pattern between units at times when we observe all outcomes, and assumes this pattern continues to hold for periods when some outcomes are missing. However, by focusing on only one of these patterns, cross-section or time series, these approaches ignore alternative patterns that may help in imputing the missing values. An alternative is to consider approaches that allow for the exploitation of both stable patterns over time, and stable patterns accross units. Such methods have a long history in the panel data literature, including the literature on two-way fixed effects, and more generally, factor and interactive fixed effect models (e.g., chamberlain1984panel, angristpischke, arellano2001panel, liang1986longitudinal, bai2003inferential, bai2009panel, bai2002determining, pesaran2006estimation, moon2015linear, moon2017dynamic). In the absence of covariates (although in this literature the coefficients on these covariates are typically the primary focus of the analyses), the common two-way fixed effect model is
More general factor models can be written as
where ${\bf U}$ is $N\times R$ and ${\bf V}$ is $T\times R$. Most of the early literature, anderson1958introduction and goldberger1972structural), focused on the thin matrix case, with $N\gg T$, where asymptotic approximations are based on letting the number of units increase with the number of time periods fixed. In the modern part of this literature bai2003inferential, bai2009panel, pesaran2006estimation, moon2015linear, moon2017dynamic, bai2017principal researchers allow for more complex asymptotics with both $N$ and $T$ increasing, at rates that allow for consistent estimation of the factors ${\bf V}$ and loadings ${\bf U}$ after imposing normalizations. In this literature it is typically assumed that the number of factors $R$ is fixed, although it is not necessarily known to the researcher. Methods for estimating the rank $R$ are discussed in bai2002determining and moon2015linear.
xu2017generalized adapts this interactive fixed effect approach to the matrix completion problem in the special case with blocked assignment, with additional applications in gobillon2013regional, kim2014divorce and hsiao2012panel. The block structure greatly simplifies the computation of the fixed rank estimators. However, this approach is not efficient, nor computationally attractive, in settings with more complex missing data patterns.
A closely related literature has emerged in machine learning and statistics on matrix completion srebro2005generalization,candes2009exact,candes2010thepower,keshavan2010matrixFew,keshavan2010matrixNoisy,gross2011recovering,recht2011simpler, rohde2011estimation, negahban2011estimation, negahban2012restricted, koltchinskii2011nuclear,klopp2014noisy, wang2018deconfounded. In this literature the starting point is an incompletely observed matrix ${\bf Y}$, and researchers have proposed low-rank matrix models as the basis for matrix completion, similar to ((ref)). The focus is not on estimating ${\bf U}$ and ${\bf V}$ consistently, but on imputing the missing elements of ${\bf Y}$. Instead of fixing the rank $R$ of the underlying matrix, a family of these estimators rely on regularization, and in particular nuclear norm regularization.
In the absence of covariates we model the $N\times T$ matrix of complete outcomes data matrix ${\bf Y}$ as
The $\varepsilon_{it}$ can be thought of as measurement error.
The goal is to estimate the matrix ${\bf L}^*$. We note that here the fixed effects are absorbed in ${\bf L}^*$ since they are two rank 1 matrices and their addition does not affect our low-rank assumption on ${\bf L}^*$.
To facilitate the characterization of the estimator, define for any matrix ${\bf A}$, and given a set of pairs of indices ${\cal O}$, the two matrices ${\bf P}_{\cal O}({\bf A})$ and ${\bf P}_{\cal O}^\perp({\bf A})$ with typical elements: \[ {\bf P}_{\cal O}({\bf A})_{it}= \left\{
\right.\hskip0.5cm {\rm and}\ \ {\bf P}_{\cal O}^\perp({\bf A})_{it}= \left\{
\right. \] A critical role is played by various matrix norms, summarized in Table (ref). Some of these depend on the singular values, where, given the full Singular Value Decomposition (SVD) $ {\bf L}_{N\times T}={\bf S}_{N\times N} {\bf \Sigma}_{N\times T} {\bf R}_{T\times T}^\top,$ the singular values $\sigma_i({\bf L})$ are the ordered diagonal elements of ${\bf \Sigma}$.
Now consider the problem of estimating ${\bf L}^*$. Directly minimizing the sum of squared differences,
does not lead to a useful estimator: if $(i,t)\in{\cal M}$ the objective function does not depend on $L_{it}$, and for pairs $(i,t)\in{\cal O}$ the estimator would simply be $Y_{it}$. To address this we regularize the problem by adding a penalty term $\lambda\|{\bf L}\|$ to the objective function in ((ref)), for some choice of the norm $\|\cdot\|$ and a scalar $\lambda$. However, since we don not wish to regularize the fixed effects (that are included in ${\bf L}^*$), we estimate them explicitly by introducing variables $\Gamma\in{\mathbb{R}}^{N\times 1}$ and $\Delta\in{\mathbb{R}}^{T\times 1}$, and the variable ${\bf L}$ will be used for estimating the remaining low-rank components of ${\bf L}^*$. This is conceptually similar to not regularizing the intercept term in LASSO estimator, to reduce the bias created by the regularization term hastie2009elements.
\paragraph{The estimator:} The general form of our proposed estimator for ${\bf L}^*$ is $\hat{\bf L}+\hat\Gamma{\mathbf{1}}_T^\top+{\mathbf{1}}_N\hat\Delta^\top$ where
Compared to the setting studied by candes2009exact, candes2010matrix, mazumder2010spectral, we include the fixed effects $\Gamma$ and $\Delta$. Although formally the fixed effects can be subsumed in the matrix ${\bf L}$ ($\Gamma{\mathbf{1}}_T^\top$ and ${\mathbf{1}}_N \Delta^\top$ are both rank one matrices), in practice, including these fixed effects separately and not regularizing them greatly improves the quality of the imputations. This is partly because compared to the settings studied in the matrix completion literature the fraction of observed values is relatively high, and so these fixed effects can be estimated accurately. The penalty factor $\lambda$ will be chosen through cross-validation that will be described at the end of this section. We will call this the Matrix-Completion with Nuclear Norm Minimization (MC-NNM) estimator.
Other commonly used Schatten norms would not work as well for this specific problem. For example, the Fr\"obenius norm on the penalty term would not have been suitable for estimating ${\bf L}^*$ in the case with missing entries because the solution for $L_{it}$ for $(i,t)\in{\cal M}$ is always zero (which follows directly from the representation of $\|{\bf L}\|_F=\sum_{(i,t)\in[N]\times[T]}L_{it}^2$). The rank norm is not computationally feasible for large $N$ and $T$ if the cardinality and complexity of the set ${\cal M}$ are substantial. Formally, the optimization problem withs the rank norm is NP-hard. In contrast, a major advantage of using the nuclear norm is that the resulting estimator can be computed using fast convex optimization programs, e.g. the {\sc soft-impute} algorithm by mazumder2010spectral that will be described next.
\paragraph{Calculating the Estimator:} For simplicity, let us first assume that there are no fixed effects (so that we do not need to estimate $\Gamma$ and $\Delta$). The algorithm for calculating our estimator goes as follows. Given the SVD for ${\bf A}$, ${\bf A}={\bf S}{\bf \Sigma}{\bf R}^\top$, with singular values $\sigma_1({\bf A}),\ldots,\sigma_{\min(N,T)}({\bf A})$, define the matrix shrinkage operator
where $\tilde{\bf \Sigma}$ is equal to ${\bf \Sigma}$ with the $i$-th singular value $\sigma_i({\bf A})$ replaced by $\max(\sigma_i({\bf A})-\lambda,0)$. Now start with the initial choice ${\bf L}_1(\lambda,{\cal O})={\bf P}_{\cal O}({\bf Y})$. Then for $k=1,2,\ldots,$ define,
until the sequence $\left\{{\bf L}_k(\lambda,{\cal O})\right\}_{k\ge 1}$ converges. The limiting matrix $\hat{\bf L}(\lambda,{\cal O})=\lim_{k\rightarrow\infty}{\bf L}_k(\lambda,{\cal O}) $ is our estimator given the regularization parameter $\lambda$. For the case that we are estimating fixed effects, after each iteration of obtaining ${\bf L}_{k+1}$, we can estimate $\Gamma_{k+1}$ and $\Delta_{k+1}$ by using the first order conditions since they only appear in the squared error term. We would also replace the ${\bf P}_{\cal O}({\bf Y})$ term in (ref) by ${\bf P}_{\cal O}({\bf Y}-\Gamma_k{\mathbf{1}}_T^\top-{\mathbf{1}}_N\Delta_k^\top)$.
\paragraph{Cross-validation:} The optimal value of $\lambda$ is selected through cross-validation. We choose $K$ (e.g., $K=5)$ random subsets ${\cal O}_k\subset{\cal O}$ with cardinality $\lfloor|{\cal O}|^2/NT\rfloor$ to ensure that the fraction of observed data in the cross-validation data sets, $|{\cal O}_k/|{\cal O}|$, is equal to that in the original sample, $|{\cal O}|/(NT)$. We then select a sequence of candidate regularization parameters $\lambda_1>\cdots>\lambda_L=0$, with a large enough $\lambda_1$, and for each subset ${\cal O}_k$ calculate $ \hat{\bf L}(\lambda_1,{\cal O}_k),\ldots,\hat{\bf L}(\lambda_L,{\cal O}_k) $ and evaluate the average squared error on ${\cal O}\setminus{\cal O}_k$. The value of $\lambda$ that minimizes the average squared error (among the $K$ produced estimators corresponding to that $\lambda$) is the one chosen. It is worth noting that one can expedite the computation by using $\hat{\bf L}(\lambda_i,{\cal O}_k)$ as a warm-start initialization for calculating $\hat{\bf L}(\lambda_{i+1},{\cal O}_k)$ for each $i$ and $k$.
\paragraph{Confidence intervals:} Studying asymptotic distribution of ${\bf L}^*-\hat{\bf L}$ in order to construct confidence intervals is beyond the scope of this paper and is an interesting future research question. However, one can use re-sampling methods to view statistical fluctuations of the imputed matrix. For example, one can again choose $K$ random subsets ${\cal O}_k\subset{\cal O}$ and construct a cross-validated estimator $\hat{\bf L}^{(k)}$ for each set ${\cal O}_k$. Then, for each entry $(i,t)$ use statistical fluctuations of $\{\hat L_{it}^{(k)}\}_{k\in[K]}$ to construct a confidence interval for $L^*_{it}$, related to the use of permutation metods in the synthetic control literature (abadie2010synthetic).
In the second contribution of this paper we discuss the relation between the matrix completion estimator and the horizontal (unconfoundedness), vertical (synthetic control) and difference-in-differences approaches. To faciliate the discussion, we focus on the case with the set of missing pairs ${\cal M}$ containing a single pair, unit $N$ in period $T$, ${\cal M}=\{(N,T)\}$. In that case the various previously proposed versions of the vertical and horizontal regressions are both directly applicable, although estimating the coefficients may require regularization depending on the relative magnitude of $N$ and $T$.
The observed data are ${\bf Y}$, an $N\times T$ matrix with the $(N,T)$ entry missing. We can partition this matrix as \[ {\bf Y}=\left(
\right),\] where ${\mathbf{Y}_0}$ is a $(N-1)\times (T-1)$ matrix, and $\mathbf{y}_1$ and $\mathbf{y}_2$ are $(N-1)$ and $(T-1)$ component vectors, respectively.
In this case the matrix completion, horizontal regression, vertical regression, synthetic control regression, the elastic net version, and difference-in-differences estimators are very closely related. They can all be charactized as focusing on the exact same objective function, but differing in the regularization and additional restrictions imposed on the parameters of the objective function.
To make this precise, define for a given positive integer $R$, an $N\times R$ matrix ${\bf A}$, an $T\times R$ matrix ${\bf B}$, an $N$-vector $\gamma$ and a $T$-vector $\delta$ the objective function
For any pair of positive integers $K$ and $L$, let $\mathbb{M}^{K,L}$ be the set of all $K\times L$ real-valued matrices. When $R=0$, we take the product ${\bf A}{\bf B}^\top$ to be the $N\times T$ matrix with all elements equal to zero. First note that simply minimizing $Q({\bf Y};R,{\bf A},{\bf B},\gamma,\delta)$ over the rank $R$, the matrices ${\bf A}$, ${\bf B}$ and the vectors $\gamma$ and $\delta$, \[ \min_{R\in\{0,1,\ldots,\min(N,T)\}}~~\min_{{\bf A}\in\mathbb{M}^{N,R},{\bf B}\in\mathbb{M}^{T,R},\gamma\in\mathbb{M}^{N,1},\delta\in\mathbb{M}^{T,1}} Q({\bf Y};R,{\bf A},{\bf B},\gamma,\delta),\] has multiple solutions for the imputations $\hat{\bf Y}_{NT}$ where $\hat{\bf Y}={\bf A}{\bf B}^\top+\gamma{\mathbf{1}}_T^\top+{\mathbf{1}}_N\delta^\top$. By choosing the rank $R$ to the minimum of $N$ and $T$, we can find for any pair $\gamma$ and $\delta$ a solution for ${\bf A}$ and ${\bf B}$ such that $P_{\cal O} \left({\bf Y}-{\bf A}{\bf B}^\top-\gamma{\mathbf{1}}_T^\top-{\mathbf{1}}_N\delta^\top\right)$ has all elements equal to zero, with different values for $\hat{\bf Y}_{NT}$.
The implication is that we need to add some structure to the optimization problem. The next result shows that horizontal regression, vertical regression, the Abadie-Diamond-Hainmueller synthetic control estimator, the difference-in-differences estimator, and the nuclear norm minimization matrix completion can all be expressed as minimizing $Q({\bf Y};R,{\bf A},{\bf B},\gamma,\delta)$ under different restrictions on, or with different approaches to regularization of the unknown parameters $(R,{\bf A},{\bf B},\gamma,\delta)$. The following theorem lays out these differences in hard restrictions and regularization approaches. Here the minimization for $R$ is over the set $\{0,1,2,\ldots,\min(T,N)\}$, and the minimization for ${\bf A}$ and ${\bf B}$ is over the sets $\mathbb{M}^{N,R}$ and $\mathbb{M}^{T,R}$ respectively.
The proof for this result is in \S (ref).
In this section we focus on the case that there are no covariates or fixed effects, and provide theoretical results for the estimation error. Let $L_{\max}$ be a positive constant such that $\|{\bf L}^*\|_{\textrm{max}}\le L_{\max}$ (recall that $\|{\bf L}^*\|_{\textrm{max}}=\max_{i,t}|{\bf L}_{it}^*|$). We also assume that ${\bf L}^*$ is a deterministic matrix. Then consider the following estimator for ${\bf L}^*$.
First, we start by introducting some new notation. Recall that for each positive integer $n$ notation $[n]$ refers to the set of integers $\{1,2,\ldots,n\}$. For any two real numbers $a$ and $b$, we denote their maximum by $a\vee b$. In addition, for any pair of integers $i,n$ with $i\in[n]$ define $e_i(n)$ to be the $n$ dimensional column vector with all of its entries equal to $0$ except the $i^{th}$ entry that is equal to $1$. In other words, $\{e_1(n), e_2(n),\ldots,e_n(n)\}$ forms the standard basis for ${\mathbb{R}}^n$. For any two matrices ${\bf A},{\bf B}$ of the same dimensions define the inner product $\langle {\bf A},{\bf B} \rangle\equiv\mathrm{trace}({\bf A}^\top {\bf B})$. Note that with this definition, $\langle {\bf A},{\bf A}\rangle=\|{\bf A}\|_F^2$.
Next, we describe a random observation process that defines the set ${\cal O}$. Consider $N$ independent random variables $\{t_i\}_{i\in[N]}$ on $[T]$ with distributions $\{\pi^{(i)}\}_{i\in[N]}$. Specifically, for each $(i,t)\in[N]\times[T]$, define $\pi_{t}^{(i)}\equiv\mathbb{P}[t_i=t]$. We also use the short notation $\mathbb{E}_\pi$ when taking expectation with respect to all distributions $\{\pi^{(i)}\}_{i\in[N]}$. Now, ${\cal O}$ can be written as $ {\cal O} = \bigcup_{i=1}^N \Big\{(i,1),(i,2),\ldots,(i,t_i)\Big\}$. The equivalent of the unconfoundedness assumption in the program evaluation literature is that the adoption dates are independent of each other and of the idiosyncratic part of the outcomes, conditional on the systematic part. Formally, we make the following assumption:
Also, for each $(i,t)\in{\cal O}$, we use the notation ${\bf A}_{it}$ to refer to $e_{i}(N)e_{t}(T)^\top$ which is a $N$ by $T$ matrix with all entries equal to zero except the $(i,t)$ entry that is equal to $1$. The data generating model can now be written as \[ Y_{it}=\langle {\bf A}_{it}, {\bf L}^* \rangle + \varepsilon_{it}\,,~~~ \forall ~(i,t)\in{\cal O}\,, \] where noise variables $\varepsilon_{it}$ satisfy Assumptions (ref)-(ref).
Note that the number of control units ($N_{\text{c}}$) is equal to the number of rows that have all entries observed (i.e., $N_{\text{c}}=\sum_{i=1}^N\mathbb{I}_{\{t_i=T\}}$). Therefore, the expected number of control units can be written as $\mathbb{E}_\pi[N_{\text{c}}]=\sum_{i=1}^N\pi_T^{(i)}$. Defining \[ p_{\text{c}} \equiv \min_{1\le i\le N}\pi_T^{(i)}\,, \] we expect to have (on average) at least $N p_{\text{c}}$ control units. The parameter $p_{\text{c}}$ will play an important role in our main theoretical results. To provide some intuition, assume ${\bf L}^*$ is a matrix that is zero everywhere except in its $i^{th}$ row. Such ${\bf L}^*$ is clearly low-rank. But recovering the entry $L^*_{iT}$ is impossible when $i_t<T$ which means $\pi_T^{(i)}$ cannot be too small. Since $i$ is arbitrary, in general, $p_{\text{c}}$ cannot be too small.
The main result of this section is the next theorem (proved in \S (ref)) that provides an upper bound for $\|{\bf L}^*-\hat{\bf L}\|_F/\sqrt{NT}$, the root-mean-squared-error (RMSE) of the estimator $\hat{\bf L}$.
\paragraph{Interpretation of Theorem (ref):} In order to see when the RMSE of $\hat{\bf L}$ converges to zero as $N$ and $T$ grow, we note that the right hand side of (ref) converges to $0$ when ${\bf L}^*$ is low-rank ($R$ is constant), { and $p_{\text{c}}\gg (\sqrt{1/T}\vee \sqrt{1/N})\log^{3/2}(N+T)$. For example, when $T$ is the same order as $N$, a sufficient condition for the latter is that the lower bound for the average number of control units ($Np_{\text{c}}$) grows faster than a constant multiple of $\sqrt{N}\log^{3/2}(N)$.} In \S (ref) we will discuss how the estimator $\hat{\bf L}$ should be modified to obtain a sharper result that would hold for a smaller number of control units.
\paragraph{Comparison with existing theory on matrix-completion:} Our estimator and its theoretical analysis are motivated by and generalize existing research on matrix-completion srebro2005generalization, mazumder2010spectral, candes2009exact,candes2010thepower,keshavan2010matrixFew,keshavan2010matrixNoisy,gross2011recovering,recht2011simpler, rohde2011estimation, negahban2011estimation, negahban2012restricted, koltchinskii2011nuclear,klopp2014noisy. The main difference is in our observation model ${\cal O}$. Existing papers assume that entries $(i,t)\in{\cal O}$ are independent random variables whereas we allow for a time series dependency structure. In particular this includes the staggered adoption setting where if $(i,t)\in{\cal O}$ then $(i,t')\in{\cal O}$ for all $t'< t$. The impact of this additional correlation is that the estimation error deteriorates significantly, compared to the ones in prior literature. For example, as discussed above, in the case of { $N=T$}, in order to have a consistent estimation we need more data. Specifically, a factor { $\sqrt{N}$} (up to logarithmic factors) more entries per column should be observed, than in the matrix completion literature.
The objective of this section is to compare the accuracy of imputation for the matrix completion method with previously used methods. In particular, in a real data matrix ${\bf Y}$ where no unit is treated (no entries in the matrix are missing), we choose a subset of units as hypothetical treated units and aim to predict their values (for time periods following a randomly selected initial time). Then, we report the average root-mean-squared-error (RMSE) of each algorithm on values for the pseudo-treated (time, period) pairs. In these cases there is not necessarily a single right algorithm. Rather, we wish to assess which of the algorithms generally performs well, and which ones are robust to a variety of settings, including different adoption regimes and different configurations of the data.
We compare the following five estimators:
The comparison between MC-NNM and the two versions of the elastic net estimator, HR-EN and VT-EN, is particularly salient. In much of the literature researchers choose ex ante between vertical and horizontal type regressions. The MC-NNM method allows one to sidestep that choice in a data-driven manner.
We use the control units from the California smoking data studied in abadie2010synthetic with $N=38, T=31$. Note that in the original data set there are $39$ units but one of them (state of California) is treated which will be removed in this section since the untreated values for that unit are not available. We then artificially designate some units and time periods to be treated, and compare predicted values for those unit/time-periods to the actual values.
We consider two settings for the treatment adoption:
In each case, the average RMSE, for different ratios $T_0/T$, is reported in Figure (ref). For clarity of the figures, for each $T_0/T$, while all 95% sampling intervals of various methods are calculated using the same ratio $T_0/T$, in the figure they are slightly jittered to the left or right. In the simultaneous adoption case, DID generally does poorly, suggesting that the data are rich enough to support more complex models. For small values of $T_0/T$, SC-ADH and HR-EN perform poorly while VT-EN is superior. As $T_0/T$ grows closer to one, VT-EN, HR-EN, SC-ADH and MC-NNM methods all do well. The staggered adoption results are similar with some notable differences; VT-EN performs poorly (similar to DID) and MC-NNM is the superior approach. The performance improvement of MC-NNM can be attributed to its use of additional observations (pre-treatment values of treatment units).
In the next illustration we use a financial data set -- daily returns for $2453$ stocks over 10 years ($3082$ days). Since we only have access to a single instance of the data, in order to observe statistical fluctuations of the RMSE, for each $N$ and $T$ we create $50$ sub-samples by looking at the first $T$ daily returns of $N$ randomly sampled stocks for a range of pairs of $(N,T)$, always with $N\times T=4900$, ranging from very thin to very fat, $(N,T)=(490,10)$, $\ldots$, $(N,T)=(70,70)$, $\ldots$, $(N,T)=(10,490)$, with in each case the second half the entries missing for a randomly selected half the units (so 25% of the entries missing overall), in a block design. Here we focus on the comparison between the HR-EN, VT-EN, and MC-NNM estimators as the shape of the matrix changes. We report the average RMSE. Figure (ref) shows the results.
In the $T\ll N$ case the VT-EN estimator does poorly, not surprisingly because it attempts to do the vertical regression with too few time periods to estimate that well. When $N\ll T$, the HR-EN estimator does poorly for the same reason: it is trying to do the horizontal regression with too few observations relative to the number of regressors. The most interesting finding is that the proposed MC-NNM method adapts well to both regimes and does as well as the best estimator in both settings, and better than both in the approximately square setting.
The bottom graph in Figure (ref) shows that MC-NNM approximates the data with a matrix of rank 4 to 12, where smaller ranks are used as $N$ grows relative to $T$. This validates the fact that there is a stronger correlation between daily return of different stocks than between returns for different time periods of the same stock.
Here we provide a brief discussion on how our estimator and its analysis should be adapted to more general settings.
In \S (ref) we described the basic model, and discussed the specification and estimation for the case without covariates. In this section we extend that to the case with unit-specific, time-specific, and unit-time specific covariates. For unit $i$ we observe a vector of unit-specific covariates denoted by $X_i$, and ${\bf X}$ denoting the $N\times P$ matrix of covariates with $i$th row equal to $X_i^\top$. Similarly, $Z_t$ denotes the time-specific covariates for period $t$, with ${\bf Z}$ denoting the $T\times Q$ matrix with $t^{\text{th}}$ row equal to $Z_t^\top$. In addition we allow for a unit-time specific $J$ by $1$ vector of covariates $V_{it}$.
The model we consider is
the $\varepsilon_{it}$ is random noise. We are interested in estimating the unknown parameters ${\bf L}^*$, ${\bf H}^*$, $\gamma^*$, $\delta^*$ and $\beta^*$. This model allows for traditional econometric fixed effects for the units (the $\gamma_i^*$) and time effects (the $\delta_t^*$). It also allows for fixed covariate (these have time varying coefficients) and time covariates (with individual coefficients) and time varying individual covariates. Note that although we can subsume the unit and time fixed effects into the matrix ${\bf L}^*$, we do not do so because we regularize the estimates of ${\bf L}^*$, but do not wish to regularize the estimates of the fixed effects.
The model can be rewritten as
Here ${\bf L}^*$ is in ${\mathbb{R}}^{N\times T}$, ${\bf H}^*$ is in ${\mathbb{R}}^{P\times Q}$, $\Gamma^*$ is in ${\mathbb{R}}^{N\times 1}$ and $\Delta^*$ is in ${\mathbb{R}}^{T\times 1}$. An slightly richer version of this model that allows linear terms in covariates can be defined as by
where $\tilde{{\bf X}}=[{\bf X}|{\bf I}_{N\times N}]$, $\tilde{{\bf Z}}=[{\bf Z}|{\bf I}_{T\times T}]$, and \[ \tilde{{\bf H}}^* = \left[
\right] \] where ${\bf H}_{XZ}^*\in{\mathbb{R}}^{P\times Q}$, ${\bf H}_{Z}^*\in{\mathbb{R}}^{N\times Q}$, and ${\bf H}_{X}^*\in{\mathbb{R}}^{P\times T}$. In particular,
From now on, we will use the richer model (ref) but abuse the notation and use notation ${\bf X},{\bf H}^*,{\bf Z}$ instead of $\tilde{{\bf X}},\tilde{{\bf H}}^*,\tilde{{\bf Z}}$. Therefore, the matrix ${\bf H}^*$ will be in ${\mathbb{R}}^{(N+P)\times(T+Q)}$.
We estimate ${\bf H}^*$, ${\bf L}^*$, $\delta^*$, $\gamma^*$, and $\beta^*$ by solving the following convex program,
Here $\|{\bf H}\|_{1,e}=\sum_{i,t} |H_{it}|$ is the element-wise $\ell_1$ norm. We choose $\lambda_L$ and $\lambda_H$ through cross-validation.
Solving this convex program is similar to the covariate-free case. In particular, by using a similar operator to ${\rm shrink}_\lambda$, defined in \S (ref), that performs coordinate descent with respect to ${\bf H}$. Then we can apply this operator after each step of using ${\rm shrink}_\lambda$. Coordinate descent with respect to $\gamma$, $\delta$, and $\beta$ is performed similarly but using a simpler operation since the function is smooth with respect to them.
In previous sections we only focused on imputing ${\bf Y}(0)$ to solve the treatment effect estimation problem. We note that this approach allows for very general assumptions on the treatment effect. For example if treatment effect has no (low-dimensional) patterns, imputing ${\bf Y}(0)$ is the best one can do because ${\bf Y}(1)$ would not have any pattern that can be used for imputation. We also note that in many of the applications there are very few treated unit/periods, so imputing the missing entries in ${\bf Y}(1)$ would be much more challenging in practice.
However, when the treatment effect is constant or has a low-rank pattern we can extend our approach and leverage the additional data from ${\bf Y}(1)$. We describe these next.
One drawback of MC-NNM is that it does not take into account the time series nature of the observations. It is likely that the ${\boldsymbol{\varepsilon}}_{it}$ are correlated over time. We can take this into account by modifying the objective function. Let us consider this in the case without covariates, and, for illustrative purposes, let us use an autoregressive model of order one. Let ${\bf Y}_{i\cdot}$ and ${\bf L}_{i\cdot}$ be the $i^{th}$ row of ${\bf Y}$ and ${\bf L}$ respectively. The original objective function for ${\cal O}=[N]\times[T]$ is \[ \frac{1}{|{\cal O}|}\sum_{i=1}^N\sum_{t=1}^T (Y_{it}-L_{it})^2+\lambda_L \|{\bf L}\|_* =\frac{1}{|{\cal O}|} \sum_{i=1}^N (Y_{i\cdot}-L_{i\cdot})(Y_{i\cdot}-L_{i\cdot})^\top+\lambda_L \|{\bf L}\|_* .\] We can modify this to $\sum_{i=1}^N (Y_{i\cdot}-L_{i\cdot}){\bf \Omega}^{-1}(Y_{i\cdot}-L_{i\cdot})^\top/|{\cal O}|+\lambda_L \|{\bf L}\|_*$, where the choice for the $T\times T$ matrix ${\bf \Omega}$ would reflect the autocorrelation in the ${\boldsymbol{\varepsilon}}_{it}$. For example, with a first order autoregressive process, we would use $\Omega_{ts}=\sigma^2\rho^{|t-s|}$, with $\rho$ an estimate of the autoregressive coefficient. Similarly, for the more general version ${\cal O}\subset[N]\times[T]$, we can use the function \[ \frac{1}{|{\cal O}|} \sum_{(i,t)\in{\cal O}}\sum_{(i,s)\in{\cal O}} (Y_{it}-L_{it})[{\bf \Omega}^{-1}]_{ts}(Y_{is}-L_{is})+\lambda_L \|{\bf L}\|_*\,. \]
Another limitation of MC-NNM is that it puts equal weight on all observed elements of the difference ${\bf Y}-{\bf L}$ (ignoring the covariates). Ultimately we care solely about predictions of the model for the missing elements of ${\bf Y}$, and for that reason it is natural to emphasize the fit of the model for elements of ${\bf Y}$ that are observed, but that are similar to the elements that are missing. In the program evaluation literature this is often achieved by weighting the fit by the propensity score, the probability of outcomes for a unit being missing.
We can do so in the current setting by modelling this probability in terms of the covariates and a latent factor structure. Let the propensity score be $e_{it}=\mathbb{P}(W_{it}=1|X_i,Z_t,V_{it})$, and let ${\bf E}$ be the $N\times T$ matrix with typical element $e_{it}$. Let us again consider the case without covariates. In that case we may wish to model the assignment ${\bf W}$ as \[ {\bf W}_{N\times T}={\bf E}_{N\times T}+\boldsymbol{\eta}_{N\times T}.\] We can estimate this using the same matrix completion methods as before, now without any missing values: \[\hat{{\bf E}}=\arg\min_{{\bf E}}\frac{1}{N T} \sum_{(i,t)} \left(W_{it} - e_{it} \right)^2+\lambda_L \|{\bf E}\|_*\,. \] Given the estimated propensity score we can then weight the objective function for estimating ${\bf L}^*$: \[\hat{{\bf L}}=\arg\min_{{\bf L}} \frac{1}{|{\cal O}|} \sum_{(i,t) \in {\cal O}} \frac{\hat e_{it}}{1-\hat e_{it}}\left(Y_{it} - L_{it} \right)^2+\lambda_L \|{\bf L}\|_*\,. \]
Recall from \S (ref) that the average number of control units is $\sum_{i=1}^N\pi_T^{(i)}$. Therefore, the fraction of control units is $\sum_{i=1}^N\pi_T^{(i)}/N$. However, the estimation error in Theorem (ref) depends on $p_{\text{c}}=\min_{1\le i\le N}\pi_T^{(i)}$ rather than $\sum_{i=1}^N\pi_T^{(i)}/N$. The reason for this, as discussed in \S (ref) is due to special classes of matrices ${\bf L}^*$ where most of the rows are nearly zero (e.g, when only one row is non-zero). In order to relax this constraint we would need to restrict the family of matrices ${\bf L}^*$. An example of such restriction is given by negahban2012restricted where they assume ${\bf L}^*$ is not too spiky. Formally, they assume the ratio $\|{\bf L}^*\|_\textrm{max}/\|{\bf L}^*\|_F$ should be of order $1/\sqrt{NT}$ up to logarithmic terms. To see the intuition for this, in a matrix with all equal entries this ratio is $1/\sqrt{NT}$ whereas in a matrix where only the $(1,1)$ entry is non-zero the ratio is $1$. While both matrices have rank $1$, in the former matrix the value of $\|{\bf L}^*\|_F$ is obtained from most of the entries. In such situations, one can extend our results and obtain an upper bound that depends on $\sum_{i=1}^N\pi_T^{(i)}/N$.
Another possible extension of Theorem (ref) is to the cases where ${\bf L}^*$ may have high rank, but most of its singular values are small. More formally, if $\sigma_1\ge \cdots > \sigma_{\min(N,T)}$ are singular values of ${\bf L}^*$, one can obtain upper bounds that depend on $k$ and $\sum_{r=k+1}^{\min(N,T)}\sigma_{r}$ for any $k\in[\min(N,T)]$. One can then optimize the upper bound by selecting the best $k$. In the low-rank case such optimization leads to selecting $k$ equal to $R$. This type of more general upper bound has been proved in some of prior matrix completion literature, e.g. negahban2012restricted. We expect their analyses would be generalize-able to our setting (when entries of ${\cal O}$ are not independent).
In \S (ref) we assumed that all entries $(i,t)$ of ${\bf Y}$ for $t\le t_i$ are observed. However, it may be possible that some such values are missing due to lack of data collection. This does not mean that any treatment occurred in the pre-treatment period. Rather, such scenario can occur when measuring outcome values is costly. In this case, one can extend Theorem (ref) to the setting with $ {\cal O} = \left[\bigcup_{i=1}^N \Big\{(i,1),(i,2),\ldots,(i,t_i)\Big\}\right]\setminus {\cal O}_{\text{miss}}$, where each $(i,t)\in \cup_{i=1}^N \{(i,1),(i,2),\ldots,(i,t_i)\}$ can be in ${\cal O}_{\text{miss}}$, independently, with probability $p$ for $p$ that is not too large.
We present new results for estimation of causal effects in panel or longitudinal data settings. The proposed estimator, building on the interactive fixed effects and matrix completion literatures has attractive computational properties in settings with large $N$ and $T$, and allows for a relatively large number of factors. We show how this set up relates to the program evaluation and synthetic control literatures. In illustrations we show that the method adapts well to different configurations of the data, and find that generally it outperforms the synthetic control estimators proposed abadie2010synthetic and the elastic net estimators proposed by doudchenko.