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.
136,135 characters · 33 sections · 36 citation commands
Longitudinal or panel data consists of a collection of observations of units (e.g., individuals, companies, countries) that are collected over time. In many applications, a subset of units are exposed to a “treatment” (e.g., drugs, regulations, or governmental interventions) beginning at some time. Given data of this type, it is frequently of interest to estimate counterfactual quantities, such as what would have been their response if they had not been treated. Such estimates underpin inferential methods for the treatment effect, corresponding to the difference between the treated and untreated responses. From the perspective of the untreated observations, performing treatment can be modeled as inducing missing data: we no longer have observations of a given unit's untreated response at all times after treatment. This problem and its variants, known as (causal) inference with panel data, find wide applications in economics, social sciences, and biomedical research (e.g., hu2008ownership,lewis2008economics,huang2008causal,imbens2015causal).
In this paper, we assume that the treatment assignments follow the staggered adoption design athey2022design,shaikh2021randomization, meaning that units may begin treatment at possibly different times, but that once initiated, the treatment is irreversible. We can collect data for the untreated unit/period pairs into a matrix, whose rows index units and columns index periods. Performing treatment can be viewed as inducing a form of missingness in this matrix; we refer the reader to (ref) for an illustration of the induced missing pattern.
Broadly speaking, in the literature on causal panel data, there are at least two main approaches for imputing missing entries in panel data: those based on unconfoundedness (e.g., rosenbaum1983central,imbens2015causal), and those based on synthetic controls (e.g., abadie2003economic,abadie2015comparative,abadie2021using). In order to understand how these methods work, suppose that we are interested in estimating the missing outcome for unit $i$ at period $t$. Approaches based on unconfoundedness seek to identify a subset of untreated units whose outcomes (before unit $i$ was treated) are similar to those of unit $i$, and use their observed outcomes at period $t$ to estimate unit $i$'s missing outcome. On the other hand, approaches based on synthetic controls estimate unit $i$'s missing outcome at time $t$ by using a weighted average of all untreated units' observed outcomes at that time, where the weights are typically determined by solving a re problem using a set of predictors that are unaffected by the treatment.
\newlength{\mywidth} {0.55cm}
Athey et al. athey2021matrix observed that both the unconfoundedness and synthetic control methods can be related to the literature on low-rank matrix completion (e.g., ExactMC09,Negahban2012restricted). Based on this connection, they proposed to impute the missing potential outcomes via the standard convex relaxation for low-rank matrix completion---namely, minimizing a least-squares objective with a nuclear norm regularization. When the potential outcome matrix is approximately low-rank, this approach can be equipped with attractive guarantees. Following this general avenue, the past few years have witnessed substantial progress in the development and analysis of matrix completion algorithms tailored for panel data with missingness (e.g., bai2021matrix,cahan2023factor,agarwal2023causal,choi2023matrix). Matrix completion methods based on nuclear norm relaxation are known to be optimal when entries are missing completely at random (e.g., wainwright2019high), but there are currently no such optimality guarantees for missingness induced by staggered adoption.
With this context, the main contribution of this paper is to propose a simple procedure, substantially less computationally intensive than nuclear norm relaxation, for inferring individual treatment effects in panel data. Most importantly, we are able to show that our procedure---despite its simplicity---is unimprovable in a sharp instance-optimal sense. More precisely:
Our procedure and theory is built upon an inferential toolbox developed for the SVD algorithm in the context of low-rank matrix denoising, which might be of independent interest. See (ref) for details.
\paragraph{A preview:} (ref) provides a high-level preview of some of the consequences of our results, and the sharpness of our theoretical predictions. In particular, let $\Mstar \in \real^{N \times T}$ be the underlying matrix of outcomes for the untreated group.
Our methodological innovation is to propose an estimate $\Mhat$ that is attractively simple: its computation involves only SVDs and other basic matrix operations. For each unobserved entry $(i, t)$, we define the rescaled error term
where $\ensuremath{\widehat{\ensuremath{\gamma}}}_{i,t}$ is a data-dependent scaling factor. On the achievable side, we provide a non-asymptotic decomposition of the rescaled error $\ensuremath{R}_{i,t}$, a particular consequence of which is that it converges in distribution to a standard $N(0,1)$ variate; see (ref) for statements of these results. Panel (a) in (ref) provides an empirical demonstration of this predicted entrywise asymptotic normality. Moreover, we prove that these estimation-theoretic and inferential guarantees are optimal in a strong sense: our variance estimates $\ensuremath{\widehat{\ensuremath{\gamma}}}_{i,t}$ converge to a population quantity $\ensuremath{\ensuremath{\gamma}^*}_{i,t}$ that is defined via a non-asymptotic Bayesian Cram\'{e}r--Rao lower bound; see (ref) and the discussion following (ref) for these claims. Panel (b) provides empirical confirmation of the instance optimality of our method; we construct an ensemble of problems whose difficulty is indexed by an “incoherence” parameter. Panel (b) compares Monte Carlo estimates of the mean-squared error obtained by our procedure with the lower bound predicted by (ref); note the excellent agreement. We refer the reader to (ref) for full details on the simulation set-ups that underlie these empirical studies.
So as to put our contributions in context, let us discuss some related lines of work. Low-rank matrix completion has been studied in great detail when entries are assumed to be missing completely at random; the behavior of the standard convex relation based on the nuclear norm minimization is now very well-understood. Initial investigations ExactMC09,recht2010guaranteed,CanTao10 primarily focused on the noiseless scenario and examined the minimal sample size required for exact recovery. In the context of noisy matrix completion (where the observed entries are corrupted by random noise), the last decade has seen the establishment of optimal estimation guarantees for convex relaxation Negahban2012restricted,klopp2014noisy,chen2019noisy. A more recent line of work chen2019inference,xia2021statistical has focused on debiasing techniques for convex relaxation, showing how it is possible to construct confidence intervals for each entry of the unknown matrix. We refer to the reader to the survey chi2018nonconvex for an overview of other non-convex approaches to matrix completion (e.g., burer2003nonlinear,KesMonSew2010,srebro2004learning).
Athey et al. athey2021matrix pioneered the use of convex matrix completion for imputing potential outcomes in panel data. They provided Frobenius norm bounds on the error in estimating the full matrix; such a bound can be viewed as providing estimation guarantees for averaged treatment effects. More recently, Choi and Yuan choi2023matrix proposed a more refined approach that first divides the missing entries into smaller groups, and then applies the convex relaxation to each group. They derived entrywise error bounds that are sharper than the Frobenius norm error bound from the paper athey2021matrix, and also established asymptotic normality for some statistics of interest. Agarwal et al. agarwal2023causal proposed an algorithm based on synthetic nearest-neighbors, for which they established error bounds and asymptotic normality. There is also a line of work bai2021matrix,cahan2023factor on models satisfying certain factor constraints, to which we compare in more detail in (ref). Other work on related but distinct estimation problems for causal panel data include the papers farias2021learning,lei2023estimating,choi2023inference,agarwal2020synthetic.
Our algorithm involves spectral techniques, and it is natural that our analysis of it has connections to past work on spectral methods abbe2017entrywise,cai2019subspace,yan2021inference,zhou2023deflated. Notably, some past work el2015impact,sur2017likelihood,ma2017implicit,chen2020bridging has made effective use of “leave-one-out” methods; see the monograph chen2020spectral for an overview. Among the technical contributions of this paper is a natural generalization of this idea, which we refer to as “leave-one-block-out”. As will be clarified by our analysis, this extension is essential in providing a sharp characterization of the subspace perturbation error that arises in the panel data model.
\paragraph{Notation:} For a positive integer $n$, we adopt the shorthand $[n] \coloneqq \{ 1 , \ldots , n \}$, For a matrix $\bm{A} \in \mathbb{R}^{n_1 \times n_2}$ and subsets $\mathcal{I} = \{ i_1, i_2, \ldots, i_I\} \subseteq [n_1]$ and $\mathcal{J} = \{ j_1, j_2, \ldots, j_J \} \subseteq [n_2]$, define $\bm{A}_{\mathcal{I} , \mathcal{J}} \in \mathbb{R} ^ { I \times J}$ to be the submatrix of $\bm{A}$ comprising the intersection of its $I$ rows indexed by $\mathcal{I}$ and $J$ columns indexed by $\mathcal{J}$, specifically $ (\bm{A}_{\mathcal{I}, \mathcal{J}})_{k, l}=A_{i_k, j_l} $ for any $k \in [ I ]$ and $l \in [ J ]$. In addition, we abbreviate $\bm{A}_{\mathcal{I}, \cdot} \coloneqq \bm{A}_{ \mathcal{I} , [ n_2 ] }$ and $\bm{A}_{ \cdot, \mathcal{J} } \coloneqq \bm{A}_{ [n_1] , \mathcal{J}}$. If a subset $\mathcal{I} = \{ i \}$ is a singleton, we may directly use the index $i$ to replace $\mathcal{I}$ for simplicity. We use $\Vert \bm{A} \Vert$, $\Vert \bm{A} \Vert_{\mathrm{F}}$ and $\Vert \bm{A} \Vert_{\infty}$ (respectively) to denote the spectral norm, Frobenius norm and entrywise $\ell_{\infty}$ norms of a matrix $\bm{A}$. For any symmetric matrices $\bm{A}$ and $\bm{B} \in \mathbb{R}^{n\times n}$, the relation $\bm{A} \succeq\bm{B}$ (resp. $\bm{A}\preceq \bm{B}$) means that $\bm{A}-\bm{B}$ (resp. $\bm{B}-\bm{A}$) is positive semidefinite. For any invertible matrix $\bm{H} \in \mathbb{R}^{r\times r}$ with SVD $\bm{U} \bm{\Sigma} \bm{V}^\top$, define its sign matrix $\mathsf{sgn} ( \bm{H} ) \coloneqq \bm{U} \bm{V}^\top$. We also define the function $\mathsf{svds}(\bm{A},r)$ to output the truncated rank-$r$ SVD $(\bm{U}, \bm{\Sigma}, \bm{V})$ of $\bm{A}$, where the columns of $\bm{U} \in \mathbb{R}^{n_1 \times r}$ and $\bm{V} \in \mathbb{R}^{n_2 \times r}$ are the top-$r$ left/right singular vectors, and $\bm{\Sigma} \in \mathbb{R}^{r \times r}$ is a diagonal matrix containing the top-$r$ singular values.
In this section, we provide a more precise description of the observation model, and its reformulation in terms of matrix estimation with missing entries.
We consider a panel data setting in which there are $\Nunit$ units over $\Time$ periods. For each unit $i \in [N]$ and at each time period $t \in [T]$, we observe a response $Y_{i,t}$. The interpretation of this observation depends on whether or not unit $i$ has undergone treatment by time $t$. We assume that treatment follows the staggered adoption design, meaning that each unit $i \in [\Nunit]$ may differ in the time that it is first exposed to treatment, and the treatment is irreversible (e.g., athey2022design,athey2021matrix,shaikh2021randomization).
More precisely, for each unit $i$, we let the integer $t_i$ denote the time at which unit $i$ was first exposed to the treatment; for completeness, we set $t_i = \infty$ if unit $i$ never undergoes the treatment. Using the notation of potential outcomes, let $Y_{i,t}(0)$ be the mean outcome of unit $i$ at time $t$ if the unit is never exposed to the treatment. For a given unit $i$ and for all times $t = 1, 2, \ldots, t_i-1$, we assume that the observed response $Y_{i,t}$ is a noisy version of this mean outcome---that is
where $\{ E_{i,t} \}$ are independent noise variables. On the other hand, for any time index $t \geq t_i$, we do not model any joint structure between the observation $Y_{i,t}$ and $Y_{i,t}(0)$, because we expect that the former might also depend on the adoption time $t_i$.
For each unit $i$ and each time $t \geq t_i$, our goal is to estimate the individual treatment effect (ITE)
corresponding to the difference between the treated response $Y_{i,t}$ (that we observe) and and the (unobserved) untreated mean $Y_{i,t}(0)$. Without further assumptions, none of our observations give any information about $Y_{i,t}(0)$ for all $t \geq t_i$; thus, the ITEs are not identifiable in general. Additional assumptions are required to ensure identifiability, and following a line of previous work, out approach is to assume that the untreated mean outcomes have low-rank structure.
In order to state this assumption more precisely, we introduce an $\Nunit \times \Time$ matrix that is populated by the untreated mean outcomes---that is, matrix $\bm{M}^\star \in \RR^{N \times T}$ with entries $M_{i,t}^\star = Y_{i,t}(0)$. Given the previously described observation model, we observe
Consequently, the problem of estimating each individual treatment effect $\tau_{i,t}$ can reduced to estimating the missing entries $\{M_{i,t}^\star \mid i\in[N], t \geq t_i\}$. As in a line of past work athey2021matrix,choi2023matrix, we adopt the assumption that $\bm{M}^\star$ is a low-rank matrix (i.e., has rank much smaller than $\min \{N, T \}$), under which it becomes possible to estimate $\bm{M}^\star$ even in the missing data setting given here.
We now give a precise formulation in terms of low-rank matrix completion.
We assume that the matrix $\MatTre{}$ has rank $r \ll \min \{\Nunit,\Time \}$. Define the subset $\Omega \subseteq [N] \times [T]$ of indices where $(i, t) \in \Omega$ if and only if unit $i$ has not yet been treated at time $t$. In terms of this notation, we can write our observation model in the compact form
where $\mathcal{P}_{\Omega}(\bm{X})$ denotes the Euclidean projection of a matrix $\bm{X}$ onto the subspace of matrices supported on $\Omega$, and $\bm{E}$ is a random noise matrix with i.i.d. $\mathcal{N} \big(0, \ensuremath{\omega}^{2}/(NT) \big)$ entries. Our scaling of the noise variance ensures that $\bm{E}$ has Frobenius norm close to $\ensuremath{\omega}$ with high probability, and facilitates interpretation of results in the sequel.
Under the staggered adoption design, we can sort the $N$ units according to the time $t_{i}$ at which they were first exposed to the treatment. With this sorting, we have $t_{1} \geq t_{2} \geq \cdots \geq t_{N}$, and we are guaranteed the existence of an integer $k \geq 1$ and partitions $N_{1}, \ldots N_{k}$ and $T_{1}, \ldots, T_{k}$ with $\sum_{i = 1}^{k}N_{i} = N$ and $\sum_{j = 1}^{k}T_{j} = T$ such that the Boolean matrix
associated with the subset $\Omega$ has the block structure as in the left panel of (ref). We also partition $\MatTre{}$ according to the pattern in $\bm{\Omega}$ to get the submatrices $\bm{M}_{i, j}^{\star} \in \mathbb{R}^{N_{i} \times T_{j}}$, as shown in the right panel of (ref). For each $i\in[k]$ and $1 \leq j \leq k-i$ define $\bm{M}_{i, j} = \bm{M}_{i, j}^{\star} + \bm{E}_{i, j}$, where $\bm{E}_{i,j}$ is the corresponding submatrix of $\bm{E}$.
Our goal is to design an algorithm that provides point estimates, along with confidence intervals, for each unseen entry $M_{i,t}^\star$ as well as corresponding ITE $\tau_{i,t} = Y_{i,t} - M_{i,t}^\star$, where $(i,t) \notin \Omega$. Since $Y_{i,t}$ is known, the estimation and inference for these two sets of quantities are equivalent, hence we will only present algorithms and results for $M_{i,t}^\star$ in the following sections.
Our estimator for the general case can be understood most easily by understanding how it applies to a special four-block structure. More precisely, suppose that a subset of $N_{1}$ units never receive the treatment, while the other $N_{2} = N - N_{1}$ units are exposed to the irreversible treatment at time $T_{1} + 1$. Then our general observation model can be written in the four-block form
Here the subscript $a$ denotes matrices of size $N_{1}$ by $T_{1}$; the subscript $b$ denotes matrices of size $N_1$ by $T - T_1$; and so on for the subscripts $c$ and $d$. As shown in the sequel, it suffices to understand how to perform estimation and inference in this scenario because our algorithm for the general staggered design first reduces the problem to a simpler problem with four-block structure. In view of the discussion at the end of the previous section, our goal is to estimate and construct confidence intervals for the entries of $\bm{M}_{d}^{\star}$, since they correspond to the counterfactual outcomes associated with the control group.
In this section, we first design an algorithm that is applicable to the simpler four-block model (ref). This algorithm is relatively easy to describe, and is the basic building block for our algorithm that applies to the general staggered design.
We first define some necessary notation for the four-block model. We begin by writing our observation matrix $\bm{M}$ in the form
The unknown matrix $\bm{M}^{\star}$ has a (compact) singular value decomposition, which we write in the form $\bm{M}^\star = \bm{U}^{\star} \bm{\Sigma}^{\star} \bm{V}^{\star\top}$, where $\bm{U}^{\star} \in \mathbb{R}^{N\times r}$ and $\bm{V}^{\star} \in \mathbb{R}^{T\times r}$ have orthonormal columns, and $\bm{\Sigma}^{\star} = \mathsf{diag}(\sigma_{1}^{\star}, \ldots, \sigma_{r}^{\star})$ is a diagonal matrix consisting of the ordered singular values $\sigma_{1}^{\star} \geq \sigma_{2}^{\star} \geq \cdots \geq \sigma_{r}^{\star}$. We partition the matrices of singular vectors $\bm{U}^{\star} \in \real^{N \times r}$ and $\bm{V}^{\star} \in \real^{T \times r}$ into two blocks as
It is also convenient to define
To motivate our algorithm, consider the ideal setting where there is no noise (i.e. $ \bm{E} = \bm{0}$). Given that $\bm{M}_{a}^{\star}$ is rank-$r$, we can first compute the matrix $\bm{U}^{\star}$ (up to a rotation) by using the SVD of the left submatrix $\bm{M}_{\mathsf{left}}^{\star}$, and then recover $\bm{M}_{d}^{\star}$ by
The idea behind the above formula is that, each column of $ \bm{M}^{\star}$ can be written uniquely as a linear combination of the columns of $ \bm{U}^{\star}$. When $ \bm{U}_{1}^{\star}$ has full rank, we can find the coefficients by regressing each column of $ \bm{M}_{b}^{\star}$ over the columns of $ \bm{U}_{1}^{\star}$. These coefficients in turn allow reconstruction of $\bm{M}_{d}^{\star}$ using the columns of $ \bm{U}_{2}^{\star}$.
In the presence of random noise, the na\"{i}ve procedure needs to be adjusted. Doing so leads to the following three-step procedure:
This procedure is summarized in (ref). As we will see in the following sections, under mild assumptions the estimate $\widehat{ \bm{M}}_{d}$ achieves full statistical efficiency for estimating $ \bm{M}_{d}^{\star}$ (even including the preconstant).
We move on to discuss the general staggered adoption design. At a high level, our procedure is based on reducing to the four-block case. Concretely, for any block $(i_{0}, j_{0})$ in (ref), we define an associated four-block estimation problem. The solution to this problem---which can be obtained using the algorithm from the preceding section---allows us to the block $\bm{M}_{i_{0}, j_{0}}^{\star}$. We obtain the four-block problem by removing a subset of data so as to simplify the observation structure; importantly, our theory to be described in the next section shows that the left-out data is ultimately not of significant statistical utility.
Let us now describe how to estimate a block $ \bm{M}_{i_{0}, j_{0}}^{\star}$, for some index $(i_{0}, j_{0})$ such that $j_{0} > k + 1-i_{0}$, so that it is not observed. Defining the indices $k_{1} =k + 1 - j_{0}$ and $k_{2} = k + 1 - i_{0}$, we construct a four-block data matrix $\bm{M}^{(i_{0}, j_{0})}$ as follows
Next we call (ref) with input $\bm{M}^{(i_{0}, j_{0})}$ to estimate its bottom-right block. Although the output provides estimation jointly for all blocks $ \bm{M}_{i, j}^{\star}$ where $k_{1} < i \leq i_{0}$ and $k_{2} < j \leq j_{0}$, we only use it to estimate the block $ \bm{M}_{i_{0}, j_{0}}^{\star}$. (For other blocks, this estimate may not be statistically efficient.) (ref) and its caption provide two examples to help illustrate how we construct the four-block data matrix for estimating each unobserved block. The complete estimation procedure is summarized in (ref).
As noted above, we show in the sequel that this algorithm achieves full statistical efficiency (including the constant pre-factors). Consequently, the data discarded while estimating each block are indeed unimportant for estimating that block, due to the lack of increase in statistical error.
We are now positioned to present our main theoretical guarantees for the algorithms described in (ref). We begin with theory for the four-block design in (ref), before extending it to the general case in (ref).
As noted previously, low-rank matrix recovery is an under-determined problem, and certain regularity conditions are required to provide theoretical guarantees. We begin by stating the conditions used in our analysis; see the discussion in (ref) for some interpretation and intuition.
\paragraph{Sub-block conditioning:} Our first set of conditions involve the top sub-blocks $\ensuremath{\bm{U}^\star}_1 \in\real^{N_1 \times r}$ and $\ensuremath{\bm{V}^\star}_1 \in \real^{T_1 \times r}$ of $\ensuremath{\bm{U}^\star}$ and $\ensuremath{\bm{V}^\star}$, respectively. In particular, we assume that there are constants $0 < \ensuremath{c_\ell} \leq \ensuremath{c_u} < \infty$ such that
\paragraph{Noise level:} Throughout, so as to avoid degenerate corner cases, we assume that $\min \{\Nunit_{1}, T_{1}\} \geq c \log (N+ T)$. In addition, some of our results require an upper bound on the noise level, as measured by the ratio $\ensuremath{\omega}/\sigma^*_r$. In particular, we define the noise functional
For a given target error level $\delta > 0$, we require that $\min \{\Nunit_1, \Time_1\}$ is sufficiently large so as to ensure that
for sufficiently small constant $c_{\mathsf{noise}} > 0$.
\paragraph{Local incoherence:} For each $i \in [N]$ and $t \in [T]$, we define the local incoherence parameters
We refer the reader to (ref) for some analysis not requiring this last condition.
With these conditions in hand, we are now ready to state an achievable result for our estimator. It provides a non-asymptotic bound, for any unobserved index $(i,t)$, on the rescaled error of our estimator in terms of a standard Gaussian variable $G_{i,t} \sim \ensuremath{\mathcal{N}}(0,1)$.
Let us explore some consequences of (ref), beginning with its uses in constructing confidence intervals. Observe that the bound (ref) guarantees that the rescaled error is close to a standard Gaussian variable. The scaling factor $\ensuremath{\ensuremath{\gamma}^*}_{i,t}$, corresponding to the error variance, depends on unknown population quantities. However, it can be estimated accurately by standard plug-in approach, leading to data-driven construction of entrywise confidence intervals. In particular, we define
Here all of the matrices $\bm{U} = \bm{U}_{\mathsf{left}}$, $\bm{V} = \bm{V}_{\mathsf{left}}$ and $\widehat{\bm{M}}_{\mathsf{upper}} = \bm{U}_{\mathsf{upper}} \bm{\Sigma}_{\mathsf{upper}} \bm{V}_{\mathsf{upper}}^\top$ can be computed from (ref). Letting $\Phi$ denote the CDF of standard Gaussian distribution, we then define the interval
See (ref) for the proof. \\
Second, while (ref) is stated in a way that facilitates its inferential use, it also allows us to derive entrywise error bounds for $\widehat{\bm{M}}_d$. If the conditions of (ref) hold with some $\delta$ decreasing to zero as $\min \{N_1, T_1 \}$ grows, then it can be verified that
In fact, even without the incoherence conditions (ref) and with a relaxation of the noise condition (ref), our analysis still provides bounds on the estimation error; see (ref) for details.
Is it possible to improve the variance term (ref) in our achievable result? In this section, we provide a negative answer to this question by showing that it matches with a Bayesian Cram\'er--Rao lower bound.
We decompose the unknown matrix as $\bm{M}^{\star} = \bm{X}^{\star} \bm{Y}^{ \star \top }$, where $\bm{X}^{\star} \coloneqq \bm{U}^{\star}( \bm{\Sigma}^{\star} )^{1/2}$ and $\bm{Y}^{\star} \coloneqq \bm{V}^{\star}( \bm{\Sigma}^{\star} )^{1/2}$ are $N \times r$ and $T \times r$ matrices, respectively. Now suppose that there is a genie that reveals all rows of $\bm{X}^{\star}$ except for the $i$-th one $\xstar \coloneqq \bm{X}_{i, \cdot}^{\star}$, and also all rows of $\bm{Y}^{\star}$ except for its $t$-th one $\ystar \coloneqq \bm{Y}_{t, \cdot}^{\star}$. Given this information, the statistician is left with the following problem:
Thus, with the aid from the genie, the causal panel data model reduces to a linear regression model with dimension $2 r$.
We use this genie-aided problem to compute lower bounds for the estimation task---both the Cram\'er-Rao lower bound as well as a local minimax risk. The local minimax risk at scale $\varepsilon > 0$ is given by
Here the infimum is taken over all estimators of the scalar $\inprod{\bm{x}}{\bm{y}}$ based on samples from the model (ref), except that the true parameters are $\bm{x}$ and $\bm{y}$ instead of $\bm{x}^\star$ and $\bm{y}^\star$.
See (ref) for the proof. \\
Under the noise assumption (ref), we have $\ensuremath{\ensuremath{\rho}_{N,T}}^2 \leq \frac{c^2_{\mathsf{noise}} \delta^2}{r + \log(N + T)}$, whence $\tfrac{R_{\mathsf{loc}} (\varepsilon_{N,T})}{\gamma_{i,t}^\star} \geq 1 - O(\delta^2)$ as $\delta$ goes to zero. Thus, we see that our our estimator (cf. (ref)) matches the local minimax risk aided by the genie, including a sharp pre-factor. Thus, our procedure is optimal in a strong instance-dependent sense.
Let us interpret the structural conditions---namely, the sub-block and local incoherence conditions---that underlie (ref).
\paragraph{Sub-block condition:} Since $\bm{M}^{\star}_a = \ensuremath{\bm{U}^\star}_1 \ensuremath{\bm{\Sigma}^\star} \ensuremath{\bm{V}^\star}_1$, the existence of the constant $\ensuremath{c_\ell} > 0$ ensures that $\MatTre{a}$ has full rank. This condition is necessary, because otherwise it would be impossible to recover the missing block $\MatTre{d}$ even in the noiseless setting.\footnote{For example, if the row rank of $\MatTre{a}$ is less than $r$, then we may find a row of $\MatTre{b}$ that is not in the subspace spanned by the rows of $\MatTre{a}$. In this case, the corresponding row of $\MatTre{d}$ is unidentifiable.}. Otherwise, we observe that $\ensuremath{\bm{U}^\star}_1$ is an $N_1 \times T_1$ sub-matrix of the orthonormal matrix $\ensuremath{\bm{U}^\star} \in \real^{N \times T}$. Since $\ensuremath{\bm{U}^\star}$ has orthonormal columns, we expect that the typical size of each entry is $1/\sqrt{N}$. Thus, the typical squared norm of each column of the sub-matrix $\ensuremath{\bm{U}^\star}_1$ is of the order $N_1/N$, which accounts for the normalization in the condition (ref) on $\ensuremath{\bm{U}^\star}_1$. Similar comments apply to the condition on $\ensuremath{\bm{V}^\star}_1$.
\paragraph{Local incoherence conditions:} Since the Frobenius norm of $\bm{U}^{\star}$ and $\bm{V}^{\star}$ are both $\sqrt{r}$, the “typical size” of their rows (in $\ell_{2}$ norm) are $\sqrt{r/N}$ and $\sqrt{r/T}$ respectively. We can see that $\mu_{i}$ (resp. $\nu_{i}$) evaluates the ratio between $\Vert \bm{U}_{i, \cdot}^{\star} \Vert_{2}$ (resp. $\Vert \bm{V}_{t, \cdot}^{\star} \Vert_{2}$) and its typical size. This definition connects to the (global) incoherence condition that is widely assumed in the low-rank matrix completion literature ExactMC09, chen2015incoherence, chi2018nonconvex, where the (global) incoherence parameters $\mu \coloneqq\max_{1 \leq i \leq N} \mu_{i}$ and $\nu \coloneqq \max_{1\leq t\leq T}\nu_{t}$ defined therein cannot be too large in order for the recovery to be possible. In contrast, our theory does not assume any global incoherence condition to hold. If our goal is to estimate and infer the counterfactual outcome or treatment effect of unit $i$ at time $t$, we just need the local incoherence parameters $\mu_{i}$ and $\nu_{t}$ to be reasonably well-behaved.
Our constraints involving the local incoherence conditions (ref) are reasonably mild. When $\min\{\Nunit_{1}, T_{1}\}\gg r$ and the noise condition ((ref)) is satisfied, we can allow both $\Vert\bm{U}_{i, \cdot}^{\star}\Vert_{2}$ and $\Vert \bm{V} _ {t , \cdot}^{\star} \Vert_{2}$ to be much larger than their typical size. The other condition (ref) requires one of $\Vert\bm{U}_{i, \cdot}^{\star}\Vert_{2}$ and $\Vert\bm{V}_{t, \cdot}^{\star}\Vert_{2}$ to be not vanishingly small compared to their typical size. In fact, even when ((ref)) does not hold, we can still characterize the entrywise distribution of $\widehat{\bm{M}}_{d}-\bm{M}_{d}^{\star}$, although it is not approximately Gaussian; we refer interested readers to (ref) for this general result.
In this section, we extend our theoretical analysis to the general (non-four-block) setting to which (ref) applies. Without loss of generality, we focus on estimating and inferring entries of an unobserved submatrix $\bm{M}_{i_{0}, j_{0}}^{\star}$. In (ref), we described how to construct a data matrix $\bm{M}^{(i_{0}, j_{0})}$ for estimating $\bm{M}_{i_{0}, j_{0}}^{\star}$ that admits the four-block design; see equation (ref) for details. (ref) provides a simple example so as to help visualize the set-up and to facilitate later discussion.
At the left of (ref), we depict a panel data problem with units separated into $k = 6$ subsets, as indicated by the large blocks. The $i$-th subset has $N_i$ units that all received the treatment at time $1+\sum_{j=1}^{i} T_j$; recall that a time larger than $T$ means that the unit never received the treatment. All the colored blocks are observed, while the uncolored blocks are missing. Suppose we want to estimate the sub-matrix $\bm{M}_{i_0,j_0}^\star$ indexed by the pair $(i_0, j_0)=(5,4)$, as indicated by a the question mark in the figure. The data matrix $\bm{M}_{(i_{0}, j_{0})}^\star$ consists of the $i_0=5$ by $j_0=4$ block matrix on the top left corner of $\bm{M}$. The gray blocks, although observed, are discarded to enforce the four-block structure. We use three colors (red, blue and yellow) to show the sub-matrices that play the role of $\bm{M}_a^\star$, $\bm{M}_b^\star$ and $\bm{M}_c^\star$ in the four-block model. The red component consists of the $k_1=3$ by $k_2=2$ block matrix on the top left corner; in general, these two integers satisfy the relations $k_1 = k + 1 - j_0$ and $k_2 = k + 1 - i_0$.
Now we define the dimension parameters $\bar{N}_{1}$, $\bar{N}_{2}$, $\bar{T}_{1}$ and $\bar{T}_{2}$ with respect to the four-block structure, given by
We also define $\bar{N} \coloneqq \bar{N}_{1}+\bar{N}_{2}$ and $\bar{T} \coloneqq \bar{T}_{1}+\bar{T}_{2}$. Let $\bm{U}_{1}^{\star} \in \mathbb{R}^{\bar{N}_{1}\times r}$ and $\bm{U}_{2}^{\star} \in \mathbb{R}^{\bar{N}_{2}\times r}$ be submatrices of $\bm{U}^{\star} \in \mathbb{R}^{N\times r}$, consisting of its first $\bar{N}_{1}$ rows and the next $\bar{N}_{2}$ rows respectively. Similarly, let $\bm{V}_{1}^{\star} \in \mathbb{R}^{\bar{T}_{1}\times r}$ and $\bm{V}_{2}^{\star} \in \mathbb{R}^{\bar{T}_{2}\times r}$ be submatrices of $\bm{V}^{\star} \in \mathbb{R}^{T\times r}$, consisting of its first $\bar{T}_{1}$ rows and the next $\bar{T}_{2}$ rows respectively. These dimension parameters are also illustrated in (ref).
In analogy with the four-block design, we first introduce the conditions required for establishing theoretical guarantees in the general staggered adoption design.
\paragraph{Sub-block conditioning:} We assume that there are constants $0 < \ensuremath{c_\ell} \leq \ensuremath{c_u} < \infty$ such that
Recall the sub-block condition (ref) for the four-block design and the discussion in (ref), here we are assuming in addition that $\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star} + \bm{U}_{2}^{\star\top} \bm{U}_{2}^{\star}$ and $\bm{V}_{1}^{\star\top} \bm{V}_{1}^{\star} + \bm{V}_{2}^{\star\top} \bm{V}_{2}^{\star}$ do not deviate too much from their typical size. This additional requirement connects the spectrum of the ground truth of $\bm{M}^{(i_{0}, j_{0})}$ (which is a submatrix of $\bm{M}^{\star}$) with that of the full ground truth $\bm{M}^{\star}$, which allows us to present cleaner results that depend on the spectrum of $\bm{M}^{\star}$ instead of its submatrices.
\paragraph{Noise level:} The noise level in this general setting takes the following form:
For any prescribed target error level $\delta > 0$, we assume that it is upper bounded by
for some sufficiently small constant $c_{\mathsf{noise}} > 0$.
\paragraph{Local incoherence:} For each $i \in [N]$ and $t \in [T]$, the local incoherence parameters $\mu_i$ and $\nu_t$ are defined as in equation (ref). For any given target error level $\delta>0$, we assume that there is some sufficiently small constant $\ensuremath{c_{\mathsf{inc}}}>0$ such that,
With this set-up, we let us describe how our theory applies to this more general setting, in particular by showing proximity of the rescaled error to a standard normal $G_{i,t} \sim \ensuremath{\mathcal{N}}(0,1)$.
See (ref) for the proof. \\
Observe that this result is the natural generalization of our previous achievable result ((ref)) for the simpler four-block setting. This correspondence makes sense, since our algorithm for the general case involves “reducing” the problem to instances of the simpler four-block case. Of course, one might wonder whether or not this reduction---which involves discarding some data that might be relevant---leads to a sub-optimal guarantee.
Interestingly, our lower bound---establishing the optimality of the variance $\gamma^\star_{i,t}$, as stated in (ref)---also extends to this multi-block setting. (We do not state this lower bound formally here, but note that our proof is actually given for the multi-block setting; see (ref) for details.) As a consequence, under the conditions of (ref), there is no appreciable loss of accuracy in solving the problem by successive reduction to the four-block setting.
Having given precise statements of our results, let us now provide a detailed comparison with related results established in past work athey2021matrix,bai2021matrix,agarwal2023causal,choi2023matrix. To simplify discussion, we ignore any log factors in making these comparisons, and use $\lesssim$ to reflect inequalities that ignore constants and such log factors. For making comparisons, we record here the estimation error bound that can can be deduced from (ref): for any entry $(i,t)$, the error is upper bounded as
with high probability.
In this section, we report the results of some numerical experiments to corroborate the computational simplicity of our algorithm, and to study the sharpness of our theoretical predictions. For certain random ensembles, we draw matrices $\bm{U} \in \mathbb{R}^{N\times r}$ from the Haar distribution over the Grassmann manifold; for simplicity, we say that any such $\bm{U}$ is a random orthonormal matrix.
All experiments in this section were based on matrices with dimensions $N = T = 500$, rank $r = 3$, and using the staggered adoption design with $k = 3$ groups with dimension parameters $N_{1} = T_{1} = 200$, $N_2 = T_2 = 200$ and $N_3 = T_3 = 100$. We used the fixed noise level $\omega = 1$.
Our goal in this set of experiments was to examine the necessity and role of incoherence conditions (ref). Focusing on estimation of the entry $(i, t) =(500, 500)$ for concreteness, we generated random problems of the following type. For an incoherence parameter $\mu \in (0, \sqrt{N/r})$, we generated an orthonormal matrix as
where $\bm{V}^{-} \in \mathbb{R}^{(T-r) \times r}$ and $\bm{O}_2 \in \mathbb{R}^{r\times r}$ are random orthonormal matrices. We consider two cases:
For each case, we performed a total of $100$ Monte Carlo trials, each time computing our estimate and forming a Monte Carlo approximation of the mean-squared error $\mathbb{E}[(\widehat{M}_{i,t}-M_{i,t}^\star)^2]$ associated with estimating entry $(i,t)$ of $M_{i,t}^\star$. Panels (a) and (b), respectively, of (ref) give plots of these empirical MSEs versus the incoherence parameter for Cases 1 and 2, respectively. For comparison, we also plot the theoretically predicted variance $\gamma^\star_{i,t}$ from equation (ref) for these two cases. For Case 1 in panel (a), we see excellent agreement between the empirical behavior and the theoretical prediction for the full range of incoherence parameters. For Case 2 in panel (b), the agreement is excellent for moderate and large values of the incoherence parameters $\mu_i$ and $\nu_t$, and we observe some discrepancies for very small values of incoherence (below $10^{-2}$). This difference arises because in this extreme regime, certain “higher-order” terms---not part of the theoretical prediction---become dominant. Overall, we see the entrywise estimation error matches the theoretical scaling (hence is statistically optimal) even when $\mu_i$ and/or $\nu_t$ are large, or at least one of $\mu_i$ and $\nu_t$ is not vanishingly small. These results support the practical validity of our algorithm and theory over a wide range of problem instances.
In this section, we devised some experiments that expose the effect of varying the matrix rank on our procedures. In order to do so, we set $N = T = 500$, the noise level $\omega=1$, and consider the same staggered adoption design with $k = 3$ groups with $N_{1} = T_{1} = 200$, $N_2 = T_2 = 200$ and $N_3 = T_3 = 100$.
As before, we consider estimating the missing entry $(i,t)=(500,500)$. Under the sub-block condition (ref), the theoretical scaling $\gamma_{i,t}^\star$ defined from equation (ref) is of the order
Although the rank $r$ does not appear explicitly, it arises implicitly via the $\ell_2$-norm of two vectors $\bm{U}_{i,\cdot}^\star$ and $\bm{V}_{t,\cdot}^\star \in \mathbb{R}^r$. Recall the definition of $\bm{U}^\star (\cdot) \in \mathbb{R}^{N\times r}$ and $\bm{V}^\star (\cdot) \in \mathbb{R}^{T \times r}$ from equation (ref), we consider the following two cases:
(ref) compares the mean-squared estimation $\mathbb{E}[(\widehat{M}_{i,t}-M_{i,t}^\star)^2]$ and the theoretical scaling $\gamma^\star_{i,t}$ for these two cases. Panel (a) shows results for Case 1, where we expect the scaling with rank to be roughly linear; here for the given dimensions $N = T = 500$, we see good agreement between empirical and theory up to around rank $r \approx 25$. On the other hand, panel (b) shows results for the more challenging Case 2, where the theoretical effect of the rank is more severe (growing as $r^{3/2}$). In this case, we see good agreement between empirical and theoretical prediction up until around rank $r \approx 20$. Given that the matrix dimension $N=T=500$ and $N_1 = T_1 = 200$, we can see that our algorithm and theory remain valid for over a reasonably large range of ranks $r$.
Finally, we did some simple experiments to explore the computational efficiency of our methods. Both of our procedures (cf. Algorithms (ref) and (ref)) involve only the computation of partial SVDs, and solving least-squares problems. Viewing the rank $r$ as an order one quantity, under the four-block setting, the computational complexity of Algorithm (ref) is dominated by the cost of computing the top-$r$ SVD of a $N$ by $T_1$ matrix and a $N_1$ by $T$ matrix; the cost of doing so scales $O(N T_1+N_1 T)$, apart from log factors. For the general staggered adoption setting with $k$ groups, the computational complexity is
In rough terms, this upper bound scales as $O(k^2 NT)$, so that when the number of groups $k$ is small, the computational cost is of the same order as cost of reading the data (i.e., the $N T$ entries of the matrix).
To give some sense of the typical runtime of our algorithm, we conducted some experiments on square matrices ($N = T$) with the side length $N$ varying over the set $\{2^9,2^{10},\ldots,2^{15}\}$, and with fixed rank $r = 3$. In all cases, we ran trials on a staggered adoption design with $k = 3$ groups with $N_{1} = T_{1} = \lfloor 0.4N \rfloor$ and $N_2 = T_2 = \lfloor 0.4N \rfloor$. Our code is Python-based, making use of the standard numpy and scipy routines for linear algebra.\footnote{One could likely speed up our code substantially by incorporating various speed-ups, e.g., from randomized numerical linear algebra; the current results show runtimes without any such code optimization.} Based on our discussion in the previous paragraph, we expect that the runtime of Algorithm (ref) scales quadratically as $O(N^2)$; this is confirmed by the simulation results shown in (ref). We note that this quadratic scaling is much faster than SDP-based algorithms like convex relaxation considered in previous literature (e.g., athey2021matrix), which become computationally burdensome even for moderate matrix dimensions (e.g., when $N > 3000$).
In this section, we provide the proofs of certain results, including our achievable guarantee in the four-block case ((ref)) in (ref); our data-dependent procedure for computing confidence intervals ((ref)) in (ref); and our lower bounds ((ref)) in (ref).
In all cases, we defer the proofs of various intermediate results, of a more technical nature, to the appendices. In addition, we defer the proof of (ref) to (ref).
Throughout this and other proofs, we use the shorthand notation $\spicy \coloneqq \log(N + T)$.
Our proof is based on a key decomposition, which we state as a separate lemma in its own right. It involves the matrix
as well as the error terms
See (ref) for the proof. \\
Armed with this result, we are now equipped to prove (ref). We begin with the additive decomposition $Z_{i,t} = \alpha_{i,t} + \beta_{i,t} + \lambda_{i,t}$, where
Observe that $\alpha_{i,t}$ and $\beta_{i,t}$ are independent zero-mean Gaussian random variables with variances
where (i) and (ii) both follow from (ref). Then we can check that
Here step (i) follows from (ref), step (ii) utilizes the bound (ref), whereas step (iii) follows from the bounds (ref) and (ref). Conditioned on the random matrix $\bm{E}_{b}$, we have
Therefore, for any fixed $(i,t)$, with probability at least $1 - O((N + T)^{-10})$, we have
Here the constant $5$ in step (i) comes from the fact that $\mathbb{P}(X\geq 5\sqrt{\spicy} ) = O((N+T)^{-10})$ for $X\sim\mathcal{N}(0,1)$, step (ii) uses the concentration bound \[ \big\Vert \bm{\zeta} ( \bm{E}_{b})_{\cdot,t} \big\Vert_{2} \leq C_g \sigma \big\Vert \bm{\zeta} \big\Vert \sqrt{ r + \spicy }, \quad \text{where} \quad \bm{\zeta} = \bm{V}_{1}^{\star}( \bm{V}_{1}^{\star\top} \bm{V}_{1}^{\star})^{-1}( \bm{\Sigma}^{\star})^{-1}( \bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star})^{-1} \bm{U}_{1}^{\star\top}, \] which is a consequence of (ref) in (ref), whereas step (iii) follows from (ref). Taking the bound on $|\lambda_{i,t}|$ and (ref) together yields
Let $G_{i,t} = \alpha_{i,t} + \beta_{i,t}$ and $\gamma_{i,t}^\star = \mathsf{var} (\alpha_{i,t}+\beta_{i,t})$, we can check that
where step (i) follows from (ref) and step (ii) follows from equations (ref) and (ref).
Our proof is based on decomposing the estimation error as
where
With these definitions in hand, we now state two auxiliary lemmas that bound these two terms.
See (ref) for the proof.
See (ref) for the proof.
Combining equation (ref) with these two lemmas yields $ \bm{M}_{d} - \bm{M}_{d}^{\star} = \bm{Z} + \bm{\Delta}$, where the random matrix $\bm{Z}$ was previously defined (ref), and $\bm{\Delta} \mydefn \bm{\xi}_{\bm{A}} + \bm{\xi}_{\bm{B}}$. The residual matrix $\bm{\Delta}$ satisfies, with probability exceeding $1-O((N+T)^{-9})$, that $ | \Delta_{ i , t} | \leq C_4 \sum_{k=1}^4 b_k (i,t) $ holds for each $(i,t)$ with some universal constant $C_4 = C_A + C_B$.
Naturally, a central step in the proof is to bound the error incurred by using $\widehat{\gamma}_{i,t}$ as an estimate of $\gamma_{i,t}^\star$. Let us summarize this auxiliary claim here: with probability at least $1 - O((N+T)^{-10})$, there exists some constant $\ensuremath{C_{\gamma}}>0$ such that
Here $\delta \in (0,1)$ is the tolerance parameter in the statement of (ref). We return to prove this result in (ref).
Taking the auxiliary claim (ref) as given, let us prove the coverage guarantee stated in (ref). Define the random variable
where $G_{i,t}$ is the standard Gaussian variable defined in (ref).
By applying (ref), we find that with probability at least $1-O((N+T)^{-10})$,
Next we observe that
and
where claim (i) follows from (ref), and step (ii) follows from the auxiliary claim (ref) and its consequence $\widehat{\gamma}_{i,t} \geq \gamma_{i,t}^\star / 2$, provided that $\delta$ is sufficiently small and $\min\{N_1,T\}$ is sufficiently large.
Combining these facts with our earlier bound (ref) and defining $C_v = 12 \ensuremath{C_{\gamma}}$, we find that
where step (iii) holds provided that $\min\{ N_1, T \} \geq \delta^{-2} r \spicy$.
We now use this high probability bound on $\ensuremath{V_{i,t}}$ to establish our coverage guarantee. For any $x \in \mathbb{R}$, from the definition of $\ensuremath{V_{i,t}}$, we have
where step (c) follows from equality $(\star)$ in equation (ref); and step (d) follows from the Lipschitz property of $\Phi$.
Finally, by taking $x=\pm \Phi^{-1} (1-\alpha/2)$ and combining the upper and lower bounds (ref) and (ref), we find that
as claimed in (ref).
To simplify presentation, we introduce the shorthand $\sigma = \omega/\sqrt{NT}$, along with the associated estimate $\widehat{\sigma}^2 = \widehat{\omega}^2 / NT$. Using this notation, we write
Next we bound each of the terms $\beta_1$ and $\beta_2$ in turn.
\paragraph{Bounding $\beta_1$:} In order to bound $\beta_1$, we need to control the estimation error associated with $\widehat{\sigma}^2$. The following result provides such a bound:
See (ref) for the proof. \\
We now apply the relation (ref) to bound $\beta_1$, thereby obtaining
Collecting together the two bounds (ref) and (ref) on $\beta_1$ and $\beta_2$, respectively, yields the upper bound
Let $\ensuremath{C_{\gamma}} = \sqrt{2} \ensuremath{C_{\sigma}} +1$ to achieve the claimed bound (ref).
We first decompose $\widehat{\sigma}^{2}$ into a sum of three terms as
As our analysis will show, the first quantity concentrates around $\sigma^2$. Thus, it is natural to bound the estimation error associated with $\widehat{\sigma}^2$ in the following way:
By applying Bernstein's inequality vershynin2016high, we can conclude that for some constant $\ensuremath{C_{\mathsf{b}}}>0$, as long as $N_1 T \geq \spicy$, we find that
provided that $N_1 T \geq \spicy$.
In order to bound the terms $\alpha_2$ and $\alpha_3$, we make use of the following lemma.
See (ref) for the proof.\\
Combining (ref) with the bounds (ref) yields
Putting together the upper bounds for $\alpha_1$, $\alpha_2$ and $\alpha_3$ yields the claimed bound (ref).
In this section, we prove the Cram\'{e}r--Rao and local minimax bounds stated in the main text. While our formal statement ((ref)) only covers the four-block case, here we actually provide a proof for the multi-block case, so that we match our general set of achievable results ((ref)).
Let us recall the set-up for our genie-aided problem associated with estimating entry $\bm{M}_{i,t}^{\star}$ of the matrix. The only unknown quantities are the vectors $\xstar \coloneqq \bm{X}_{i,\cdot}^{\star}$ and $\ystar \coloneqq \bm{Y}_{t,\cdot}^{\star}$, and our goal is to estimate the inner product $\bm{M}_{i,t}^{\star} = \inprod{\xstar}{\ystar}$ based on linear measurements of the form
\paragraph{Cram\'{e}r--Rao lower bound:} We first compute the Cram\'{e}r--Rao lower bound (CRLB) for the problem. Using $\bftheta = (\xstar, \ystar)$ to denote the unknown parameters, our goal is to estimate the functional $\Psi(\bftheta) \coloneqq \inprod{\xstar}{\ystar}$. The CRLB is given by $\nabla \Psi(\bftheta)^T \FishMat^{-1}(\bftheta) \nabla \Psi(\bftheta)$, where $\FishMat(\bftheta)$ is the Fisher information matrix, and the gradient takes the form $\nabla \Psi(\bftheta) =
^T$. Given that our measurements are linear with Gaussian observation noise with variance $\ensuremath{\omega}^2/(N T)$, the Fisher information matrix takes the form
Putting together the pieces, we find that the CRLB is given by
Here step (i) uses the fact that $\bm{X}^\star = \bm{U}^\star (\bm{\Sigma}^\star)^{1/2}$ and $\bm{Y}^\star = \bm{V}^\star (\bm{\Sigma}^\star)^{1/2}$, along with some algebra.
\paragraph{Computation of the local minimax bound:} We now turn to computation of the local minimax lower bound, where we make use of a Bayesian form of the CRLB. We begin by observing that for any $\varepsilon > 0$, the local minimax risk for estimating $\bm{M}^\star_{i,t}$ can be lower bounded as
where $\pi$ is any prior supported on the Cartesian product space. $\ensuremath{\mathcal{B}}_\infty (\bm{x}^\star , \varepsilon ) \times \ensuremath{\mathcal{B}}_\infty (\bm{y}^\star , \varepsilon )$.
Following a standard avenue for computing Bayesian CRLBs polyanskiy2023information, we first define the squared cosine density function $g(u) = \cos^2 (\pi u / 2)$ for $u \in [-1 , 1]$. Using this building block, we then define the product prior distribution $\pi(\theta_1, \ldots, \theta_{2r}) = \prod_{j=1}^{2 r} \pi_j(\theta_j)$ with marginals
With this choice of prior, we can compute the lower bound
where step (i) makes use of some standard computations with squared cosine priors (e.g., \S 29.2 in the book polyanskiy2023information); and step (ii) follows from the spectral decomposition of $\bm{M}^\star$, along with some algebra.
Thus far, we proven that $R^\star(\varepsilon) \geq \widetilde{\gamma}_{i,t} (\varepsilon) $. It remains to show that with the choice $\varepsilon = \varepsilon_{N,T} \equiv \sqrt{\sigma_r^\star} \max \{1/\sqrt{N}, 1/\sqrt{T} \}$ we have
from which it will follow that $R^\star(\varepsilon) \geq \widetilde{\gamma}_{i,t} (\varepsilon) = \gamma^\star_{i,t} - \Delta^\star_{i,t} \geq \big(1 - \frac{\ensuremath{c_u} \pi^2}{c_\ell} \ensuremath{\ensuremath{\rho}_{N,T}}^2 \big) \gamma^\star_{i,t}$, as claimed.
From the definitions of $\widetilde{\gamma}_{i,t}(\varepsilon)$ and $\gamma_{i,t}^\star$, we have
Now observe that $\Vert (\bm{\Sigma}^\star)^{-1} \Vert \leq \frac{1}{\sigma_r^\star}$, and that the sub-block condition condition (ref) ensures that
Combining these facts, we find that
Combining inequalities (ref) and (ref) yields the bound
where step (i) follows from our choice of $\varepsilon$, and the definition of $\ensuremath{\ensuremath{\rho}_{N,T}}$. This establishes the auxiliary claim (ref), and thus the claimed lower bound.
In this paper, we proposed and analyzed a simple algorithm for estimation and inference for panel data with missingness induced by staggered adoption design. It is appealing from a computational point of view, since it is non-iterative, requiring only elementary matrix operations and singular value decomposition. At the same time, we demonstrated that its statistical properties are also attractive, in that it achieves estimation accuracy, in an elementwise sense, that matches non-asymptotic lower bounds applicable to any estimator. Moreover, we showed how this theory enables data-driven construction of confidence intervals for the missing potential outcomes. We note that our development relies on a more general inferential toolbox for the SVD algorithm in the matrix denoising model. It introduces the idea of “leave-one-block-out”, an analysis technique that we suspect will be useful for other problems.
Moving forward, let us outline some open questions, along with broader directions for future study. Our current theory relies on a noise condition (ref) and an incoherence condition (ref); it remains unclear whether or not these conditions are improvable, for instance in terms of dependence on the rank $r$. Moreover, it would be interesting to extend our guarantees to accommodate matrices that are only approximately low-rank, along with more general noise models. Finally, while the the staggered adoption design arises frequently, it would be interesting to understand optimal procedures for more general mechanisms of missing data.
Y. Yan was supported by the Norbert Wiener Postdoctoral Fellowship from MIT. M. J. Wainwright was partially supported by ONR grant N00014-21-1-2842 and NSF grant DMS-2311072. We thank Eric Xia for helpful discussions.
\printbibliography