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,225 characters · 17 sections · 56 citation commands
Inference for Low-rank Completion without Sample Splitting with Application to Treatment Effect Estimation
\def\spacingset#1{ {#1}} \spacingset{1}
\onehalfspacing \if11 { \affil[1]{{ Columbia University, 1255 Amsterdam Avenue New York, NY 10027, United States}} \affil[2,3]{{ Rutgers University, 75 Hamilton St, New Brunswick, NJ 08901, United States}} } \fi
\spacingset{2} \if01 {
} \fi
\spacingset{1.3}
{\it Keywords:} Matrix completion; Nuclear norm penalization; Two-step least squares estimation; Approximate factor model; Causal inference
{\it JEL Classification:} C12, C14, C33, C38, C55
\spacingset{1.5}
{8pt} {8pt} \setlength\intextsep{8pt} {4pt}
The task of imputing the missing entries of a partially observed matrix, often dubbed as matrix completion, is widely applicable in various areas. In addition to the well-known application to recommendation systems (e.g., the Netflix problem), this problem is applied in a diverse array of science and engineering such as collaborative filtering, system identification, social network recovery, and causal inference.
In this paper, we focus on the following approximate low-rank model with a factor structure:
where $Y$ is an $N \times T$ data matrix which is subject to missing, $M$ is a latent matrix of interest, and $\mathcal{E}$ represents a noise contamination. Importantly, $M$ is assumed to be an approximate low-rank matrix having an approximate factor structure $M \approx \beta F'$, where $\beta$ is factor loadings and $F$ is latent factors. In addition, we allow some entries of $Y$ to be unobserved by defining an indicator $\omega_{it}$, which equals one if the $(i,t)$ element of $Y $ is observed, and zero otherwise. In this practical setting, we provide the inferential theory for each entry of $M$, regardless of whether its corresponding entry in $Y$ is observed or not.
One of the widely used methods for the low-rank matrix completion is the nuclear norm penalization and it has been intensively studied in the last decade. candes2009exact, candes2010matrix, koltchinskii:2011, negahban2012restricted, and chen2020noisy provide statistical rates of convergence for the nuclear norm penalized estimator and a branch of studies including beck:2009, cai:2010, mazumder:2010, ma:2011, and parikh2014proximal provide algorithms to compute the nuclear norm penalized estimator. However, research on inference is still limited. This is because the shrinkage bias caused by the penalization, as well as the lack of the closed-form expression of the estimator, hinders the distributional characterization of the estimator.
We contribute to the literature by providing an inferential theory of the low-rank estimation without sample splitting. Our estimation procedure consists of the following main steps:
Note that steps 2-3 are only conducted once without further iterations.
An important contribution is that we do not rely on the sample splitting to make inference, but simply use the full (observed) sample in every step of our procedure. There are at least three advantages to avoid sample splitting. First, the resulting estimator using sample splitting is unstable and random even conditioning on the data. Second, sample splitting requires relatively large $T$ in practice, because it practically works with only $T/2$ observations. This is demanding in applied micro applications when $T$ is just a few decades. In the simulation study, we show that the performance of the estimator using sample splitting is worse than that of the estimator without sample splitting when $T$ is relatively small. Lastly, sample splitting increases computational costs in multiple tests because for each target time `$t$', we need to use different sample splitting.
Technically, we apply a new approach to showing the negligibility of the potential bias terms, by making use of a hypothetically defined auxiliary leave-one-out (ALOO) estimator. We emphasize the word “auxiliary" because it is only introduced in the technical argument, but not implemented in the estimation. So it is a hypothetical estimator, which is to be shown that it is
i) asymptotically equivalent to the initial estimator for $\beta$ in Step 1 and
ii) independent of the sample used in the least squares estimation, namely, the sample in period $t$. Using the ALOO estimator, we can separate out the part in the initial estimator for $\beta$, which is correlated with the sample in period $t$. Once we separate out the correlated part, we can enjoy a similar effect to the sample splitting. And we show the separated correlated part is sufficiently small. Importantly, the leave-one-out estimator only appears in the proof as an auxiliary point of the initial estimator for $\beta$, so we do not need to compute it in the estimation procedure, which allows us to remove the sample splitting step without implementing any additional steps.
Empirically, we apply the proposed procedure to making inference for the impact of the presidential vote on allocating the U.S. federal budget to the states. We find the states that supported the incumbent president in past presidential elections tend to receive more federal funds and this tendency is stronger for the loyal states than the swing states. In addition, this tendency is stronger after the 1980s.
Very recently, some studies proposed the ways of achieving unbiased estimation for the inference of the nuclear norm penalized estimator. chernozhukov2019inference, chernozhukov2021inference propose a two-step least square procedure with sample splitting, which estimates the factors and loadings successively using the least square estimations. As we discussed earlier, sampling splitting comes with several undesirable costs.
The idea of the ALOO estimator has been employed in other recent works such as ma2019implicit,chen:2019inference,chen2020nonconvex,chen2020noisy,yan2021inference as well. Among them, in particular, chen:2019inference pioneered using this idea to convex relaxation of low-rank inference. This paper has some important contributions compared to chen:2019inference.
Other related works on inference include xia2021statistical, xiong2020large, and jin2021factor. We compare these methods with ours in simulations.
Lastly, a comparison with other literature that takes advantage of a low-rank model to estimate the treatment effect would be helpful. The close connection between low-rank completion and treatment effect estimation was first made formal by athey2021matrix who showed that the nuclear norm regularization can be useful for causal panel data by presenting the convergence rate of the estimator. Another line of research proposes inferential theories under weaker assumptions on the treatment assignment with other restrictions. farias2021learning allow the assignment of the treatment that can depend on historical observations while focusing on the estimation of the average treatment effect. agarwal2021causal and bai2021matrix consider the case where the assignment is not random but has a certain block structure that often occurs in causal panel data.\footnote{In agarwal2021causal, a certain submatrix for estimation has a block structure.} In addition, arkhangelsky2021synthetic propose an estimator that is more robust than the conventional difference-in-differences and synthetic control methods by using a low-rank fixed effect model with the homogeneous treatment effect assumption.
This paper is organized as follows. Section (ref) provides the model and the estimation procedure as well as our strategy for achieving the unbiased estimation. Section (ref) gives the asymptotic results of our estimator. Section (ref) provides the inferential theory for the average treatment effect estimator as an application. Section (ref) presents an empirical study about the impact of the president on allocating the U.S. federal budget to the states to illustrate the use of our inferential theory. Section (ref) includes the simulation studies. Section (ref) concludes.
There are a few words on our notation. For any matrix $A$, we use $\left\VertA\right\Vert_F$, $\left\VertA\right\Vert$, and $\left\VertA\right\Vert_*$ to denote the Frobenius norm, operator norm, and nuclear norm respectively. $\left\VertA\right\Vert_{2,\infty}$ denotes the largest $l_2$ norm of all rows of a matrix $A$. $\mathrm{vec}(A)$ is the vector constructed by stacking the columns of the matrix $A$ in order. Also, $\psi_r(A)$ is $r$th largest singular value of $A$. $\psi_{\max}(A)$ and $\psi_{\min}(A)$ are the largest and the smallest nonzero singular value of A. For any vector $B$, $\mathrm{diag}(B)$ is the diagonal matrix whose diagonal entries are $B$. $a\asymp b$ means $a/b $ and $b/a$ are $O_P(1)$.
We consider the following nonparametric panel model subject to missing data problem: \[ y_{it} = h_t\left( \zeta_{i} \right) + \varepsilon_{it}, \] where $y_{it}$ is the scalar outcome for a unit $i$ in a period $t$, $h_t(\cdot)$ is a time-varying nonparametric function, $ \zeta_{i}$ is a unit-specific latent state variable, $\varepsilon_{it}$ is the noise, and $\omega_{it} =1{\{y_{it} \text{ is observed}\}}$. Here, $\{h_t( \cdot ), \zeta_{i}, \varepsilon_{it}\}$ are unobservable. In the model, the (latent) unit states $ \zeta_{i}$ have a time-varying effect on the outcome variable through $h_t(\cdot)$. This model can be written in (ref) using the sieve representation. Suppose the function $h_t(\cdot)$ has the following sieve approximation: \[ h_t( \zeta_{i}) = \sum_{r=1}^{K} \kappa_{t,r} \phi_{r}( \zeta_{i}) + M^{R}_{it} = \beta_{i}^{\prime}F_t + M^{R}_{it} = M^{\star}_{it} + M^{R}_{it} , \] where $\beta_{i} = (\phi_{1}( \zeta_{i}),\dots, \phi_{K}( \zeta_{i}))^{\prime}$ and $F_t = (\kappa_{t,1},\dots,\kappa_{t,K})^{\prime}$. Here, $M_{it}^{R}$ is the sieve approximation error and, for all $1\leq r \leq K$, $\phi_{r}( \zeta_{i})$ is the sieve transformation of $\zeta_i$ using the basis function $\phi_{r}( \cdot)$ and $\kappa_{t,r}$ is the sieve coefficient. Then, $$ M=[M_{it}]_{N\times T},\quad M_{it}= h_t(\zeta_{i}) $$ can be successfully represented as the approximate factor structure.
In matrix form, we can represent the model as
where we denote $Y = [y_{it}]_{N\times T}$, $M = [M_{it}]_{N\times T}$, $M^{\star} = [M^{\star}_{it}]_{N\times T}$, $M^{R} = [M^{R}_{it}]_{N\times T}$, $\beta = [\beta_{1}, \ldots , \beta_{N} ]'$, $F = [F_1, \ldots , F_T ]'$, and $\mathcal{E} = [\varepsilon_{it}]_{N\times T}$. Note that $Y$ is an incomplete matrix that has missing components.
Let $\mathcal{M} \coloneqq (\beta, F,M^{R})$ be the triplet of random matrices that compose $M$. In the paper, we allow the heterogeneous observation probability, i.e., $P(\omega_{it}=1)=p_i$ and denote $\Pi = \mathrm{diag}(p_1,\dots,p_N)$. Here, we shall assume the sieve dimension $K$ is pre-specified by researchers and propose some data-driven ways of choosing $K$ in Section (ref) of Appendix.
To accommodate the heterogeneous observation probability, this paper uses the inverse probability weighting scheme, referred to as inverse propensity scoring (IPS) or inverse probability weighting in causal inference literature (e.g., imbens:2015, little2019statistical, schnabel2016recommendations), in addition to the nuclear norm penalization:
where $\widehat{\Pi} = \mathrm{diag}(\widehat{p}_1, \dots, \widehat{p}_N)$, and $\widehat{p}_i = \frac{1}{T}\sum_{t=1}^{T}\omega_{it}$ for each $i\leq N$, $\Omega = [\omega_{it}]_{ N\times T}$ and $\circ$ denotes the Hadamard product. As noted in ma2019missing, this inverse probability weighting debiases the objective function itself. If there is heterogeneity in the observation probability, $\|\Pi^{-\frac{1}{2}}\Omega \circ\left( A - Y \right)\|_F^2$ is an unbiased estimate of $\left\VertA-Y\right\Vert_F^2$, which we would use if there is no missing entry, in the sense that $\mathbb{E}_{\Omega}[\|\Pi^{-\frac{1}{2}}\Omega \circ\left( A - Y \right)\|_F^2]=\left\VertA-Y\right\Vert_F^2$, while $\left\Vert\Omega \circ(A-Y)\right\Vert_F^2$ is biased.
Although the inverse probability weighting enhances the estimation quality, the weighting alone cannot guarantee the asymptotic normality of the estimator because of the shrinkage bias. To achieve the unbiased estimation having the asymptotic normality, we run the two-step least squares procedure. As noted previously, our estimation does not have the sample splitting steps. Our estimation algorithm is as follows:
After deriving the initial estimator of loadings from the nuclear norm penalized estimator $\widetilde{M}$, we estimate latent factors and loadings using the two-step least squares procedure. The final estimator of $M$ is then the product of the estimates for latent factors and loadings.
It is well-known that the nuclear-norm penalized estimator $\widetilde{M}$, like other penalized estimators, is subject to shrinkage bias which complicates statistical inference. To resolve this problem, we use the two-step least squares procedure, i.e., Steps 3 and 4 in Algorithm (ref). In showing the asymptotic normality of the resulting estimator $\widehat{M}$, a key challenge is to show the following term is asymptotically negligible: \[ R_t = \frac{1}{\sqrt{N}} \sum_{j=1}^N \omega_{jt}\varepsilon_{jt}(\widetilde{\beta}_j - H_1'\beta_j) \] where $H_1$ is some rotation matrix.\footnote{Another term $\frac{1}{\sqrt{N}} \sum_{j=1}^N (\omega_{jt}-p_j)\beta_j F_t'H_1^{\prime -1}(\widetilde{\beta}_j - H_1'\beta_j)$ is also to be shown negligible, but the argument is similar to that of $R_t$.} This term represents the effect of the bias of the nuclear-norm penalization since $\widetilde{\beta}_j$ is derived from the nuclear-norm penalized estimator. chernozhukov2019inference,chernozhukov2021inference resort to sample splitting to show the asymptotic negligibility of $R_t.$
Motivated by chen2020noisy, we show the asymptotic negligibility of $R_t$ without sample splitting by using two hypothetical estimators which are asymptotically equivalent to the nuclear norm penalized estimator $\widetilde{\beta}$. Namely, we consider a hypothetical non-convex iteration procedure for the low-rank regularization, where singular vectors are iteratively solved as the solution and show that this procedure can be formulated as the following two problems:
Here, $\| \cdot \|_{F,(-t)}$ denotes the Frobenius norm which is computed ignoring $t$-th column and $\| \cdot \|_{F,t}$ is the Frobenius norm of only $t$-th column. Note that the only difference between (ref) and (ref) is that the $t$-th column of the goodness of fit part in (ref) is replaced by its conditional expectation in (ref). So, $\{\omega_{jt}, \varepsilon_{jt}\}_{j \leq N}$ is excluded from the problem (ref).
We emphasize that (i) both problems defined above are non-convex; (ii) both problems are “auxiliary", meaning that they are introduced only for proofs, not actually implemented. (iii) Optimizing $L^{(-t)}(B, F) $ is an auxiliary leave-one-out (ALOO) problem, leading to the ALOO estimator $\breve{\beta}^{(-t)} $ to be discusssed below.
Because of the non-convexity, both hypothetical problems should be computed iteratively until the gradients of the non-convex loss functions become “sufficiently small.” However, the gradients do not monotonically decrease as iteration proceeds since the problem is non-convex. So, one cannot let it iterate until convergence is reached, but has to stop at the point where the gradient is small enough. The choice of this “stoping point” is crucial in the analysis of the residual terms. chen:2019inference define the stopping point using the full sample problem (ref), which potentially causes dependence problem of the leave-one-out estimators. We propose two approaches of addressing this issue.
In either way, we can successfully show the negligibility of $R_t$ uniformly in $t$ without resorting to sample splitting. We highlight that the first approach is more natural in the sense that it automatically ensures the independence that we need for term $b$. Our first approach, while technically more involved, is potentially more applicable to general machine learning inferences that rely on auxiliary leave-one-out estimators, because of the natural independence. In contrast, it is unclear whether the second approach is still applicable in other cases.
It is natural to ask why would not we define the ALOO estimator more naturally as the original estimator $\widetilde\beta$, but simply dropping the $t$ th column from the data matrix in the optimization? One of the key differences between $L^{(-t)}(B, F)$ in ((ref)) and the “more natural dropping-$t$" loss, is that the $t$ th column in the least squares part of $L^{(-t)}(B, F)$ is not simply dropped, but is replaced by its expectation: $$ \mathbb E \|\Pi^{-\frac{1}{2}}\Omega \circ\left( B F' - Y \right)\|_{F,t}^2 = \| BF'-M^{\star}\|_{F,t}^2 + C $$ where the constant $C$ does not depend on $(B,F)$. The reason for defining the ALOO loss function in this way is to gain “hypothetical efficiency", so that the ALOO estimator would be closer to the full-sample estimator.
It is easier to understand the issue using a simple example. Consider estimating the mean $\mathbb EY_t$ using iid data $Y_t$. The full-sample estimator $\widehat\mu$ is the solution to $$ \widehat\mu=\arg\min_{\mu}L(\mu),\quad \text{where } L(\mu)=\sum_{s=1}^T(Y_s- \mu)^2. $$ Now consider the ALOO version of this problem. Our definition of $L^{(-t)}(\mu)$ is not dropping $Y_t$, but replacing $(Y_t-\mu)^2$ with its expectation: $$ \breve{\mu}^{(-t)}=\arg\min_{\mu}L^{(-t)}(\mu),\quad \text{where } L^{(-t)}(\mu)=\sum_{s\neq t}(Y_s- \mu)^2+ \mathbb E (Y_t-\mu)^2. $$ The solution is then $\breve{\mu}^{(-t)}= \frac{1}{T}(\sum_{s\neq t}Y_s+\mathbb EY_t)$. Then straightforward calculations can verify that $\breve{\mu}^{(-t)}$ (although infeasible) is more efficient and “closer" to the full-sample average $\widehat\mu$ than the naive dropping-$t$ estimator $\bar Y_{-t}:=\frac{1}{T-1}\sum_{s\neq t} Y_s$. For instance, $$ \frac{\mathrm{Var} (\breve{\mu}^{(-t)} )}{\mathrm{Var}(\bar Y_{-t})} =\left(\frac{T-1}{T}\right)^2<1,\quad \frac{\mathbb E(\breve{\mu}^{(-t)}- \widehat\mu)^2}{\mathbb E(\bar Y_{-t}-\widehat\mu )^2} =\frac{T-1}{T}<1. $$
The definitions of $L^{(-t)}(B, F)$ and $L^{(-t)}(\mu)$ also fulfill the intuition of the EM algorithm, which imputes the missing data in the loss function by its conditional expectations before optimizations, rather than simply dropping the missing values.
From Algorithm (ref), we see that there is no explicit debias step. In fact, in terms of estimating the singular vector space, the singular vector estimator from the least square estimation following the nuclear norm penalization, $\widehat{F}_t$, is unbiased (up to a rotation).
To see this, note that the estimation of $F_t$ has the following maximization problem: \[\widehat{F}_t \coloneqq \operatorname*{arg\,max}_{f \in \mathbb{R}^K} Q_t(f, \widetilde{\beta})\] where $Q_t(f,B) = -N^{-1}\sum_{j=1}^N \omega_{jt}(y_{jt}-f' b_j)^2$, $B=(b_1,\dots,b_N)^{\prime}$ and $b_j$ are $K$ dimensional vectors. In this step, $\beta$ is the nuisance parameter and $F_t$ is the parameter of interest. By Taylor expansion, we have, for some invertible matrix $A$,
The first term is the score which leads to the asymptotic normality and the second term represents the effect of the $\beta$ estimation which is subject to the shrinkage bias. The second term, while is the “usual bias" in a generic machine learning inference problem, can be shown to take the form: $$ d=\sqrt{N} \varphi H_1^{-1} F_t + o_P(1) $$ for some $\varphi=o_P(1)$. It has a useful feature of being on the space of $F_t$. Making use of this fact, (ref) can be re-written as follows: \[ \sqrt{N} (\widehat{F}_t-H_2 F_t) = - \underbrace{\sqrt{N} A^{-1} \frac{\partial Q_t(H_1^{-1}F_t, \beta H_1)}{\partial f}}_{\text{asymptotically normal}} + o_P(1) \] by defining $H_2 \coloneqq (I_K + \varphi) H_1^{-1}$. Note that the non-negligible bias term in $d$ is absorbed by the rotation matrix $H_2$, and thus $\widehat{F}_t$ can unbiasedly estimate $F_t$ up to this new rotation. Then, in Step 4 of Algorithm (ref), $\widehat{\beta}$, the least square estimator using $\widehat{F}$ as a regressor, can unbiasedly estimate $\beta_i$ up to the rotation since $\widehat{F}_t$ has only a higher order bias now. As a result, the product of them estimates $M_{it}$ unbiasedly:
which allows us to conduct inference successfully. This is how the two-step least squares procedure works.
This section presents the inferential theory. We provide the asymptotic normality of the estimator of the group average of $M_{it}$. Our assumptions allow the rank $K$ to grow, but slowly. Remind the following notation: \[ h_t( \zeta_{i}) = \sum_{r=1}^{K} \kappa_{t,r} \phi_{r}( \zeta_{i}) + M^{R}_{it} = \beta_{i}^{\prime}F_t + M^{R}_{it}, \] where $\beta_{i} = (\phi_{1}( \zeta_{i}),\dots, \phi_{K}( \zeta_{i}))^{\prime}$ and $F_t = (\kappa_{t,1},\dots,\kappa_{t,K})^{\prime}$. Let $S_{\beta} = N^{-1}\sum_{i=1}^N \beta_i \beta_i'$, $S_{F} = T^{-1}\sum_{s=1}^T F_s F_s'$, and $Q = S_{\beta}^{1/2}S_{F}^{1/2}$.
First, we present some assumptions for the sieve representation. Assumption (ref) (ii) is well satisfied with a large $\nu$ if the functions $\{h_t\left( \cdot \right)\} $ are sufficiently smooth. For example, consider $h_t$ belonging to a H$\ddot{o}$lder class: for some $a,b,C>0$, $\left\lbrace h : \|D^{b}h(x_1) - D^{b}h(x_2) \| \leq C \| x_1 - x_2\|^{a}\right\rbrace.$ In addition, suppose that we take a usual basis like polynomials, trigonometric polynomials, and B-splines. Then, $\max_{i,t}|M^{R}_{it}| \leq C K^{-\nu},$ and $\nu = 2(a + b)/\text{dim}(\zeta_i).$ So, Assumption (ref) (ii) is satisfied with very large $\nu$ if $\{h_t\left( \cdot \right)\} $ are smooth. In addition, the first part of Assumption (ref) (iii) can be satisfied if the basis is a bounded basis like trigonometric basis or $\zeta_i$ has a compact support. Assumption (ref) (iv) and (v) are not restrictive, and have been verified by chernozhukov2021inference.
We assume the heterogeneous observation probability across $i$. It generalizes the homogeneous observation probability assumption which is a typical assumption in the matrix completion literature. The sub-gaussian assumption in Assumption (ref) (iii) helps us to bound $\left\Vert\Omega \circ \mathcal{E} \right\Vert$ and $\left\Vert\Omega- \Pi \bold{1}_N\bold{1}_T^{\prime}\right\Vert$.
While the serial independence of the missing data indicators $\omega_{it}$ is assumed, we allow they are cross-sectional dependence among $i$. In doing so, we assume a cluster structure in $\{1, \dots, N\}$, i.e., there is a family of nonempty disjoint clusters, $\mathcal{C}_1,\dots, \mathcal{C}_{\rho}$ such that $\cup_{g=1}^{\rho} \mathcal{C}_g = \{1, \dots, N\}$. So we divide units $\{1,...,N\}$ into $\rho$ disjoint clusters. In addition, denote the size of the largest cluster by $\vartheta$. That is, $\vartheta=\max_{g} |\mathcal{C}_g|_o$. We highlight that $\vartheta$ is allowed to increase as $N$ and $T$ increase.
Due to the cluster structure in Assumption (ref) (i), we can construct a “leave-cluster-out” estimator $\breve{\beta}^{\{-i\}}$ which is independent of the sample of unit $i$. Similarly to the idea of (ref) and (ref), we can rule out the samples of the cluster that includes unit $i$. The difference from (ref) is that we identify all the units which are in the same cluster as unit $i$ and replace their rows of the goodness of fit part by their conditional expectations.\footnote{For the formal definitions of the estimators, please refer to Section (ref) of Appendix and Remark (ref) in the section.} Together with the leave-one-out estimator $\breve{\beta}^{(-t)}$, the leave-cluster-out estimator $\breve{\beta}^{\{-i\}}$ plays a pivotal role in showing the solution of (ref) is close to that of (ref).
The parameter for the cluster size $\vartheta$ is bounded by Assumption (ref). For instance, in the case where $N \asymp T$ and $\{h_{t}(\cdot) \}_{t\leq T}$ are smooth enough, if we estimate the cross-sectional average of a certain period, the assumption requires $\vartheta \approx o(\sqrt{N / \log N} )$ since $K$ is allowed to grow very slowly when $\{h_{t}(\cdot) \}_{t\leq T}$ are smooth.
We are interested in making inference about group-averaged effects. Let $\mathcal G $ be a particular group; the object of interest is $$ \frac{1}{|\mathcal{G}|_o}\sum_{(i,t)\in\mathcal{G}} {M}_{it}=\frac{1}{|\mathcal{G}|_o}\sum_{(i,t)\in\mathcal{G}} h_t(\zeta_i). $$ Here the group of interest as $\mathcal{G} = \mathcal{I} \times \mathcal{T}$ where $\mathcal{I} \subseteq \{1, \ldots, N\}$ and $\mathcal{T} \subseteq \{1, \ldots, T\}$. We impose the following assumption on the rates of parameters. Define a sequence $\psi_{NT}$ as $\psi_{NT} \asymp \sqrt{K^{-(2\delta +1)} \sum_{i=1}^N \sum_{t=1}^Th_t^2(\zeta_i)}$. It is a lower bound of $\psi_{\min}(\beta F')$ and works as the parameter for signal. Recall that $K$ denotes the sieve dimension.
Assumption (ref) (ii) is used to bound the sieve approximation error. For this condition to be satisfied, the smoothness of $\{h_{t}(\cdot) \}_{t\leq T}$ is crucial. If $\{h_{t}(\cdot) \}_{t\leq T}$ are smooth enough, $\nu = 2(a + b)/\text{dim}(\zeta_i)$ can be arbitrarily large. Hence, Assumption (ref) (ii) can be easily satisfied with a slowly increasing $K$ as long as $\{h_{t}(\cdot) \}_{t\leq T}$ is smooth.
Assumptions (ref) (i) is the conditions about sample complexity and signal-to-noise ratio. As long as $K,\eta,\gamma$ are bounded or increase sufficiently slowly, it would be satisfied. Note that, in the cases like the cross-sectional average of a certain period t or the time average of a certain unit i, $ \min\{|\mathcal{I}|_o^{\frac{1}{2}},|\mathcal{T}|_o^{\frac{1}{2}}\} =1$. In many interesting cases, $\min\{|\mathcal{I}|_o^{\frac{1}{2}},|\mathcal{T}|_o^{\frac{1}{2}}\}$ is usually not that large. However, due to Assumption (ref) (i), we cannot derive the inferential theory in the case where both $|\mathcal{I}|_o$ and $|\mathcal{T}|_o$ are large like $|\mathcal{I}|_o=N$ and $|\mathcal{T}|_o=T$. In this case, the asymptotically normal distribution part cannot dominate other residual parts, since the convergence rate of the asymptotically normal distribution part is roughly $\frac{1}{\sqrt{N|\mathcal{T}|_o}} + \frac{1}{\sqrt{T|\mathcal{I}|_o}}$, while that of the residual term is similar to or greater than $\frac{1}{\sqrt{NT}}$ regardless of the group size. For inference, at least one part of the asymptotically normal term should dominate other residual terms. On the other hand, in terms of the convergence rate, the large sizes of $|\mathcal{I}|_o$ and $|\mathcal{T}|_o$ are beneficial, as noted in Section (ref) in Appendix. In addition, for comparison with the conditions of other low-rank literature, it would be helpful to refer to Assumption (ref) in Appendix where we consider the general low-rank model.
Under the above assumptions, Theorem (ref) shows that the estimator for the group average of $M_{it}$ has the asymptotic normality: $$\mathcal{V}_{\mathcal{G}}^{-\frac{1}{2}}\left( \frac{1}{|\mathcal{G}|_o}\sum_{(i,t)\in\mathcal{G}} \widehat{M}_{it} - \frac{1}{|\mathcal{G}|_o}\sum_{(i,t)\in\mathcal{G}} M_{it} \right) \overset{D}\longrightarrow \mathcal{N}(0,1),$$ where the asymptotic variance $\mathcal{V}_{\mathcal{G}}$ is given in the statement of Theorem (ref), and needs to be estimated. In this result, $\mathcal{G}$ can consist of either multiple columns with multiple rows or solely a certain $(i,t)$, implying that we can conduct inference for one specific element of the matrix.
To make the estimation of $\mathcal{V}_{\mathcal G}$ feasible, we consider the case of $\mathbb{E}[\varepsilon_{it}^2|\mathcal{M}] = \sigma^2_{i}$. Let $U_{i}'$ is the $i$-th row of the left singular vector of $\beta F'$ and $V_{t}'$ is the $t$-th row of the right singular vector of $\beta F'$. The following theorem gives the feasible asymptotic normality.
We now establish the semiparametric efficiency of our estimator, following a similar approach as in jankova2018semiparametric. In order to make calculation tractable, we suppose that $\omega_{it}\sim \mathrm{Bernoulli}(p)$ and $\varepsilon_{it}\sim \mathcal{N}(0,\sigma^2)$ are independent across $(i,t)$. We will focus on the case of block group, where both $|\mathcal{I}|_o$ and $|\mathcal{T}|_o$ are finite or growing slowly, satisfying $N|\mathcal{T}|_o \ll T^2 |\mathcal{I}|_o^2$ and $T|\mathcal{I}|_o \ll N^2 |\mathcal{T}|_o^2$. The other cases like cross-sectional and serial groups (e.g., $|\mathcal{I}|_o =N$ and $|\mathcal{T}|_o$ is finite or slowly growing, or vice versa) can also be attained, which are very similar to Theorem 4.2 in chernozhukov2021inference. Hence, we omit them. The novelty of our efficiency theorem is that it is for estimating the general block group.
As specified in Theorem (ref), the asymptotic variance in this case is
The following theorem shows that $s^2_*(M,p, \sigma)$ is the asymptotic Cram\'er-Rao bound for asymptotically unbiased estimators.
In this section, we propose the inference procedure for treatment effects by utilizing the asymptotic results in Section (ref). Following the causal potential outcome setting (e.g., rubin:1974, imbens:2015), we assume that for each of $N$ units and $T$ time periods, there exists a pair of potential outcomes, $y_{it}^{(0)}$ and $y_{it}^{(1)}$ where $y_{it}^{(0)}$ denotes the potential outcome of the untreated situation and $y_{it}^{(1)}$ denotes the potential outcome of the treated situation. Importantly, among potential outcomes $y_{it}^{(0)}$ and $y_{it}^{(1)}$, we can observe only one realized outcome $y_{it}^{(\Upsilon_{it})}$ where $\Upsilon_{it} = 1\{\text{unit $i$ is treated at period $t$}\}$. Hence, we have two incomplete potential outcome matrices, $Y^{(0)}$ and $Y^{(1)}$, having missing components, and the problem of estimating the treatment effects can be cast as a matrix completion problem because of the missing components in the two matrices.
Specifically, we consider the nonparametric model such that for each $\iota \in \{0,1\}$, \[y_{it}^{(\iota )} = M_{it}^{(\iota )} + \varepsilon_{it} = h_{t} ^{(\iota )} (\zeta_{i}) + \varepsilon_{it}, \] where $\varepsilon_{it}$ is the noise, $\zeta_{i}$ is a vector of unit specific latent state variables. We regard $h_t^{(\iota)}(\cdot)$ as a deterministic function while $\zeta_{i}$ is a random vector. In the model, the treatment effect comes from the difference between the time-varying treatment function $h_t^{(1)}(\cdot)$ and the control function $h_t^{(0)}(\cdot)$. Let $\omega_{it}^{(\iota )} =1{\{y_{it}^{(\iota )} \text{ is observed}\}}$. Then, $\omega_{it}^{(1)} = \Upsilon_{it}$ and $\omega_{it}^{(0)} = 1 - \Upsilon_{it}$ because we observe $y_{it}^{(1)}$ when there is a treatment on $(i,t)$ and observe $y_{it}^{(0)}$ when there is no treatment on $(i,t)$.
We suppose the following seive representation for $h_t^{(\iota)}$ : \[ h_t^{(\iota)} ( \zeta_{i})= \sum_{r=1}^{K} \kappa_{t,r}^{(\iota)} \phi_{r}( \zeta_{i}) + M^{R(\iota)}_{it},\ \ \ \ \ \ \iota \in \{0,1\} \] where $\kappa^{(\iota)}_{t,r}$ is the sieve coefficient, $\phi_{r}( \zeta_{i})$ is the sieve transformation of $\zeta_i$ using the basis function $\phi_{r}( \cdot)$ and $M_{it}^{R(\iota)}$ is the sieve approximation error. Then, by representing $\sum_{r=1}^{K} \kappa_{t,r}^{(\iota)} \phi_{r}( \zeta_{i}) $ as $\beta_{i}^{\prime}F_t^{(\iota)}$ where $\beta_{i} = [\phi_{1}( \zeta_{i}),\dots, \phi_{K}( \zeta_{i})]^{\prime}$ and $F_t^{(\iota)} = [ \kappa^{(\iota)}_{t,1},\dots,\kappa^{(\iota)}_{t,K}]^{\prime}$, $h^{(\iota)}_t(\zeta_{i})$ can be successfully represented as the approximate factor structure.
We make inference about the average treatment effect for a particular group of interest $(i,t)\in\mathcal{G}$: $$\frac{1}{|\mathcal{G}|_o}\sum_{(i,t)\in\mathcal{G}} \Gamma_{it},\quad\text{where } \Gamma_{it} = M_{it}^{(1)} - M_{it}^{(0)}.$$ The individual treatment effect $\Gamma_{it}$ is estimated by $\widehat{\Gamma}_{it} = \widehat{M}_{it}^{(1)} -\widehat{M}_{it}^{(0)}$ where $\widehat{M}_{it}^{(0)}$ and $\widehat{M}_{it}^{(1)}$ are estimators of $M_{it}^{(0)}$ and $M_{it}^{(1)}$, respectively. Hence, by implementing the estimation steps in Algorithm (ref) for each $\iota \in \{0,1\}$, we can derive the estimators for the group average of $M_{it}^{(0)}$ and $M_{it}^{(1)}$, and construct the average treatment effect estimator.
The notations are essentially the same as those in Section (ref), and we just put the superscript $(\iota)$ to all notations to distinguish the pair of potential realizations.
To illustrate the use of our inferential theory, we present an empirical study about the impact of the president on allocating the U.S. federal budget to the states. The allocation of the federal budget in the U.S. is the outcome of a complicated process involving diverse institutional participants. However, the president plays a particularly important role among the participants. Ex-ante, the president is responsible for composing a proposal, which is to be submitted to Congress and initiates the actual authorization and appropriations processes. Ex post, once the budget has been approved, the president has a veto power that can be overridden only by a qualified majority equal to two-thirds of Congress. In addition, the president exploits extra additional controls over agency administrators who distribute federal funds.
There is a vast theoretical and empirical literature about the impact of the president on allocating the federal budget to the states (e.g., cox1986electoral, anderson1991congressional, mccarty2000presidential, larcinese2006allocating, berry2010president). In particular, cox1986electoral provide a theoretical model which supports the idea that more funds are allocated where the president has larger support because of the ideological relationship between voters and the president, and larcinese2006allocating have found that states which supported the incumbent president in past presidential elections tend to receive more funds empirically. We contribute by showing the impact using our inferential theory for the heterogeneous treatment effect with a wider set of data.
Here, the hypothesis we want to test is whether federal funds are disproportionately targeted to states where the incumbent president is supported in the past presidential election. We use data on federal outlays for the 50 U.S. states with the District of Columbia from 1953 to 2018. The data are obtained from websites of the U.S. Census Bureau, NASBO (National Association of State Budget Officers), and SSA (Social Security Administration).
Following Section (ref), we set the treatment indicator as $\Upsilon_{it} = 1$ if the state $i$ supported the president of year $t$ in the presidential election, and $\Upsilon_{it} = 0$ otherwise. If the candidate whom the state $i$ supported in the previous presidential election is the same as the president at year $t$, we consider it as “treated" and otherwise, we consider it as “untreated". While applying our inferential procedure, we adopt the assumption that the treatment (whether state $i$ supported the resident in the election) is exogenously assigned, which is probably not practical, but we take our stand on this assumption in this study, and do not claim a causal interpretation of the treatment effect.
In addition, for the outcome variable $y_{it}$, we use the following ratio: $y_{it} = (\tilde{y}_{it}/\sum_{i} \tilde{y}_{it}) \times 100$ where $\tilde{y}_{it}$ is the per-capita federal grant in state $i$ at year $t$. Note that the outcome variable, $y_{it}$, is a proportion so that $\sum_i y_{it}=100$ for all $t,$ which is to treat each period equally.
Our inferential theory allows novel approaches to study the following effects:
First, Figure (ref) presents the State Effects and the corresponding t-statistics. The results suggest significantly positive treatment effects in most states. To investigate the reason of differences, we categorize states according to the number of times a state swung the party it supports in the presidential elections as in Table (ref). Together with Figure (ref), it shows that most states with large t-statistics are in “Loyal states” while other states are generally in “Swing state” or “Weak swing state”. It suggests that the treatment effect is closely related to the loyalty of states to parties.
In addition, the results for the Region Effects in Figure (ref) show that, at the 1% significant level, New England, Mid Atlantic, Plains, Rocky Mountain, and Far West have the positive treatment effects while Great Lakes, South East, and South West do not. Note that Many states in Great Lakes, South East, and South West are in “Swing states” or “Weak swing states.” As we can see in Figure (ref), “Swing states” do not have statistically significant positive treatment effects while “Loyal states” have significant positive treatment effects. This result is in line with the empirical study of larcinese2006allocating finding that states with loyal supports tend to receive more funds, while swing states are not rewarded. In addition, it is aligned with the assertion of cox1986electoral that the targeting of loyal voters can be seen as a safer investment as compared to aiming for swing voters and risk-averse political actors may allocate more funds to loyal states.
Figure (ref) shows the President Effects and the Party Effects. Despite some exceptions, there are no statistically significant positive treatment effects before Carter, while there are significant positive treatment effects after Reagan. Figure (ref) shows that before 1980, there is no significant positive treatment effect in most states, while there are significant positive treatment effects in most states after 1981. Hence, there is a substantial difference between `before 1980' and `after 1981' and the tendency that incumbent presidents reward states that showed their support in the presidential elections became significant after Reagan, that is, after the 1980s. It suggests that after the 1980s, the presidents show more influence on the allocation of federal funds to reward their supporters. Evidence is that starting from the 1980s, all presidents have put forward proposals for the introduction of presidential line-item veto and tried to increase the power of the president to control federal spending.
Finally, when testing for the treatment effects of multiple states, the tests may subject to the issue of multiple testing problems, with undesirable false discovery rates (FDR). We also address this issue by adopting the procedure of benjamini1995controlling to control the FDR at 5%. We find that the list of states with significant treatment effects is unchanged.
This section provides the finite sample performances of the estimators. We first study the performances of the estimators of $M_{it}$ and $|\mathcal{G}|_o^{-1}\sum_{(i,t)\in\mathcal{G}} M_{it}$, and then study performances of the average treatment effect estimators. To save space, some results are relegated to Appendix.
First of all, in order to check the estimation quality of our estimator, we compare the Frobenius norms of the estimation errors for several existing estimators of $M$. Our two-step least squares is labelled as “TLS". We also consider the debiased nuclear norm penalized estimators from xia2021statistical, “(Hetero) XY,” and chen:2019inference, “(Hetero) CFMY.” “(Hetero)” represents that they are modified to allow the heterogeneous observation probabilities. The comparison also includes the inverse probability weight based estimator, “IPW,” from xiong2020large, and the EM algorithm based estimator, “EM,” from jin2021factor. The plain nuclear norm penalized estimator, “Plain Nuclear,” and the TLS estimator using sample splitting, “TLS with SS,” are also considered. For the data-generating designs, we consider the following three models:
Here, $U_{t,r}$ is generated from $\mathcal{N}(2,1)$ and $\zeta_{i}$ is generated from $\text{Uniform}[0,1]$. In addition, $\varepsilon_{it}$ is generated from the standard normal distribution independently across $i$ and $t$. The observation pattern follows a heterogeneous missing-at-random mechanism where $\omega_{it}\sim \text{Bernoulli}(p_i)$ and $p_i$ is generated from Uniform $[0.3,0.7]$.
Table (ref) reports $\|\widehat{M} - M \|_F/\sqrt{NT}$ averaged over 100 replications. We highlight that the TLS shows the best performance in almost all scenarios. Only the EM is comparable to ours, but it computes much slower since it requires multi-step iterations. In contrast, our proposed method does not iterate. Also, our method always outperforms the TLS with SS. The (Hetero) XY and (Hetero) CFMY are slightly worse than ours in this experiment. Lastly, both the IPW and the Plain Nuclear show the worst performances uniformly. The IPW, being non-statistically efficient, is only slightly better than the Plain Nuclear.
Additionally, to show the relative advantage of TLS over TLS with sample splitting, Table (ref) reports $(\widehat{M}_{it} - M_{it})^2$ in the case where $T$ is small. Here, we choose $(i,t)$ randomly and fix it during replications. As we can check in the table, when $T$ is relatively small, the performance of TLS with sample splitting is much worse than that of TLS without sample splitting. Especially, in the factor model, the difference in performance is quite large.
Second, we study the finite sample distributions for standardized estimates defined as $(\widehat{M}_{it} - M_{it})/se( \widehat{M}_{it})$. For comparison, we report the results of the Plain Nuclear and the TLS with SS, in addition to the TLS. For the Plain Nuclear, we use the sample standard deviation obtained from the simulations for $se( \widetilde{M}_{it} )$ because the theoretical variance of it is unknown. For the TLS with SS, we construct the standard error following chernozhukov2019inference. Here, we consider the nonparametric models in (ref). Hereinafter, the number of replications is 1,000, and the sample size is $N=T=200.$
Figure (ref) plots the scaled histograms of the standardized estimates with the standard normal density. As we expected in theory, it shows that the standardized TLS and the standardized TLS with SS fit the standard normal distribution well, while the standardized Plain Nuclear is biased. Without sample splitting, the TLS itself provides a good approximation to the standard normal distribution so that it can be used for the inference successfully. The coverage probabilities of confidence interval in Appendix also show similar results.
Next, we study the finite sample performance of the average treatment effect estimator. Following Section (ref), for each $\iota \in \{0,1\}$, we generate the data from $y_{it}^{( \iota)} = h_t^{(\iota)}(\zeta_{i}) + \varepsilon_{it}$, where $ h_t^{(0)}(\zeta) = \sum_{r=1}^{\infty} |U_{t,r}| r^{-a}\sin(r\zeta)$, $h_t^{(1)}(\zeta) = \sum_{r=1}^{\infty} (|U_{t,r}|+2)r^{-a}\sin(r\zeta)$. The power parameter $a>1$ controls the decay speed of the sieve coefficients. The forms of the above functions and the treatment effect $\Gamma_{it}= h_t^{(1)}(\zeta_i) - h_t^{(0)}(\zeta_i)$ are in Figure (ref).
Here, $\varepsilon_{it}$ and $U_{t,r}$ are independently generated from the standard normal distribution and $\zeta_{i}$ is independently generated from $\text{Uniform}[0,1]$. The treatment pattern follows $\Upsilon_{it}\sim \text{Bernoulli}( p_i^{(1)}) $ and $p_i^{(1)} \sim \text{Uniform}[0.3,0.7]$.
Figure (ref) presents the scaled histograms of the standardized estimates of the average treatment effect estimators for the groups $\mathcal{G}_1 = \{(i,t)\}$, $\mathcal{G}_2 = \{(j,t):1\leq j \leq N \}$, and $\mathcal{G}_3 = \{(i,s): 1\leq s \leq T\}$. Here, the standard estimates are given as $$\frac{\frac{1}{|\mathcal{G}|_o}\sum_{(i,t)\in\mathcal{G}} \widehat{\Gamma}_{it} - \frac{1}{|\mathcal{G}|_o}\sum_{(i,t)\in\mathcal{G}} \Gamma_{it} }{se\left( \frac{1}{|\mathcal{G}|_o}\sum_{(i,t)\in\mathcal{G}} \widehat{\Gamma}_{it} \right) }.$$ As expected in theory, the standardized estimates of the average treatment effect estimators of all groups approximately show the standard normal distribution. In addition, the coverage probabilities of the confidence interval in Appendix also show similar results.
This paper studies the inferential theory for low-rank matrices and provides an inference method for the average treatment effect as an application. Without the aid of sample splitting, our estimation procedure successfully resolves the problem of the shrinkage bias, and the resulting estimator attains the asymptotic normality. Unlike chernozhukov2019inference, chernozhukov2021inference which exploit sample splitting, our estimation step is simple, and we can avoid some undesirable properties of sample splitting. In addition, this paper allows the heterogeneous observation probability and uses inverse probability weighting to control the effect of the heterogeneous observation probability.
For the sake of brevity, some of the technical proofs are relegated to the Supplement.