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.
76,852 characters · 16 sections · 35 citation commands
Inference for Low-rank Models without Estimating the Rank
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 { \affil[1]{Department of Statistics, Columbia University} \affil[2]{Department of Operations Research and Financial Engineering, Princeton University} \affil[3]{Department of Economics, Rutgers University} } \fi
\if10 {
} \fi
{\it Keywords:} High-dimensional models; low-rank matrix estimation; rank misspecification; weak factors
\spacingset{1.9}
The objective of this paper is to make inference about a low rank matrix without consistently estimating its rank. We consider the following linear model:
where $\bf{Y}$, $\mathbf{X}$, ${\mathbf{M}^{\star}}$, and $\mathbf{E}$ are $N \times T$ matrices, with $N$ and $T$ approaching infinity, and $\circ$ represents the matrix entry-wise product. The outcome matrix $\mathbf{Y}$ and the regressor matrix $\mathbf{X}$ are observed, and $\mathbf{E}$ is the matrix of noise. The target of interest is the matrix ${\mathbf{M}^{\star}}$, whose rank (denoted by $r$) is unknown and low compared to its dimensions. Our model incorporates various models such as noisy matrix completion, heterogeneous treatment effects estimation, and varying coefficients models.
The inferential theory for low rank matrices has been developed in the recent literature, as in chernozhukov2018inference, chernozhukov2023inference, chen:2019inference, xia2021statistical,yan2024inference. Most of these methods utilize the estimated eigenvectors for rank reductions, so we call them principal components analysis (PCA)-based methods. A standard PCA-based inference procedure can be outlined as follows chen:2019inference,xia2021statistical:
Here $\mathcal{B}_1$ captures a shrinkage bias one often encounters from the initial estimator. Clearly, one needs to start by taking the true rank (or its consistent estimator), which however requires that the signal-to-noise ratio be sufficiently high. This assumption often breaks down in finite sample applications. In fact, the rank estimation is threatened by the presence of weak factors, raising severe concerns in applications. For instance, in financial applications it is often the case that the first eigenvalue is much larger than the remaining eigenvalues, so consistent rank estimators often identify only one factor, which however violates empirical practices in asset pricing. As another example in the forecast practice, empirical eigenvalues do not decay as fast as the theory requires, so it is often difficult to determine the cut-off value to separate “spiked eigenvalues" from the remaining ones.
A simple solution is to over-estimate the rank. To avoid the risk of under-estimating the rank, one can select a sufficiently large rank, and use it throughout the inference procedure. However, when the PCA-based debiasing methods are employed with the over-estimated rank, they fail to exclude the eigenvectors that are generated from and thus strongly correlated with the noise. In addition, the over-estimated eigenvectors correspond to non-spiked eigenvalues, which are well known to be inconsistent in high dimensional settings johnstone2009consistency.
This paper makes a novel contribution to the low-rank inference literature by proposing a procedure robust to the rank over-estimation. In order to circumvent the aforementioned issues of the PCA-based method, we adopt a non-PCA based method, known as diversified projection (DP) which was recently proposed by fan2022learning in the pure factor model. Similar to the PCA, the DP is a dimension reduction method that projects the original high-dimensional space to a low-rank space. But it is much simpler than the PCA: it does not require calculating the eigenvalue/eigenvectors, so demands much weaker conditions on the eigen-structures of the low-rank matrix. We employ the DP to estimate a “larger” low-rank space, which is inconsistent to the true singular vector spaces. Nevertheless, the simple structure of the DP enables us to characterize the over-estimated rank spaces relatively easily. Importantly, we do not require correctly specifying (or consistently estimating) the true rank of the underlying matrix.
Although our adoption of the DP approach is inspired by fan2022learning, the primary difference lies in our objective. Their underlying model is a pure factor model with no missing data, and the focus is exclusively on the extracted factors. In contrast, motivated by the study of treatment effects, our main goal is to make inferences about the low-rank matrix itself. When the objective shifts from focusing on factors to directly examining the low-rank matrix, a novel statistical insight emerges: a new source of bias from over-estimating the low-rank space when the rank exceeds its true value. Contrary to the conventional view that overfitting primarily introduces additional variance, we show that in the context of low-rank inference, an over-estimated low-rank space can be characterized by a Tikhonov-type function, leading to an implicit ridge-type regularization bias. Our statistical interpretation is that the over-estimated low-rank space is highly correlated with the model's noise, and such correlation depends on a non-stochastic second moment of the noise, giving rise to the new bias. This issue would not arise if the objective were solely focused on the extracted factors for other types of inference.\footnote{For instance, fan2022learning applied estimated factors to “factor-augmented regression" problems (e.g., forecasts and post-selection inference). In such cases, additional steps often involve regressing external time series on the extracted factors, which can “average out" the effect of over-estimating the low-rank space.}
Our condition is more robust to the strength of the singular values than the PCA-based approach: we allow weak singular values of ${\mathbf{M}^{\star}}$, which just need to be stronger than $\sqrt{T\log^2 N}$. This condition is only slightly stronger than the setting of “weak factors" in onatski2012asymptotics in the pure factor model, whose condition was $\sqrt{T }$. On the contrary, the usual PCA-based approaches would require the singular values be $\sqrt{T N^c}$ for some $c\in(1/2,1]$ bai2023approximate.
A practical reason for allowing the true rank to be over-estimated is to account for confounding factors (CF) in variable selection and multiple testing. Standard statistical inference can be significantly affected by strong correlations among variables due to hidden confounding factors, so needs to be adjusted. The key practical question is determining how many CF to account for, including the special case that there are none---though this is unknown in practice, so as a precaution, statisticians often account for some CF regardless. In one of our statistical applications, we address this issue within the context of multiple testing. We show that our method remains valid across all bounded numbers of CF, including the special case where none are present.
We recommend a few specific choices for the weighting matrices in Section (ref). One of the appealing recommended choices is to use transformations of initial observations, which does not require extra data. Moreover, we also note that the use of extra information is not uncommon in a wide range of literature.
There is a large literature on estimation of high-dimensional low-rank matrices, such as candes2009exact,keshavan2010matrix,negahban2011estimation,klopp2014noisy,cai2018rate,cape2019two,abbe2020entrywise,koltchinskii2020efficient,cai2021subspace,zhu2022high among many others. In this literature, profound theories on optimal rates of convergence have been developed, whereas as we commented earlier, the distributional theory has been studied in the more recent literature. barigozzi2020consistent proposed a PCA-based method to over-estimate the rank and achieved insightful rate results.
We adopt the following notations. Let $\left\Vert\cdot\right\Vert$ and $\left\Vert\cdot\right\Vert_*$ denote the matrix operator norm and nuclear norm, respectively. Also, we use $\left\Vert\cdot\right\Vert_{2,\infty}$ to denote the largest $l_2$ norm of all rows of a matrix. We write $\sigma_{\max}(\cdot)$ and $\sigma_{\min}(\cdot)$ to represent the largest and smallest singular values of a matrix, respectively, and $\sigma_{j}(\cdot)$ to represent the $j$th largest singular value of a matrix. For a matrix $\mathbf{A}$, define $\mathrm{span}(\mathbf{A})$ as the linear space spanned by the columns of matrix $\mathbf{A}$. When $\mathbf{A}'\mathbf{A}$ is invertible, define $\mathbf{P}_\mathbf{A} = \mathbf{A} (\mathbf{A}'\mathbf{A})^{-1}\mathbf{A}'.$ For a vector $\mathbf{v}$, $\mathrm{diag}(\mathbf{v})$ represents the diagonal matrix whose diagonal entries are $\mathbf{v}$ in order. For two sequences $a_{NT}$ and $b_{NT}$, we denote $a_{NT} \ll b_{NT}$ (or $b_{NT} \gg a_{NT}$) if $a_{NT} =o(b_{NT})$, $a_{NT} \lesssim b_{NT}$ (or $b_{NT} \gtrsim a_{NT}$) if $a_{NT} =O(b_{NT})$, and $a_{NT} \asymp b_{NT}$ if $a_{NT} \lesssim b_{NT}$ and $a_{NT} \gtrsim b_{NT}$ (almost surely if random). Finally, due to page limits, all proofs and some simulation studies are provided in the appendix.
In the linear model (ref), we assume that the matrix $\mathbf{M}^{\star}$ has a low-rank structure:
where $\bm\beta$ is an $N\times r$ matrix of rescaled left singular vectors, and $\mathbf{F}$ is a $T\times r$ matrix of the rescaled right singular vectors; rescaled by the singular values. Model (ref) has numerous applications including varying coefficient models, heterogeneous treatment effects athey2021matrix, and matrix completion problems.
The diversified projection (DP) is a dimension reduction technique that projects a high-dimensional object onto a low-rank space. To illustrate the idea, consider the high-dimensional factor model:
To estimate $\mathbf{M}^{\star}$, one approach is to take advantage of the low-rank structure, by applying projections as follows:
Here, $\mathbf{P}_{\widetilde{{\bm\beta}}}$ and $\mathbf{P}_{\widetilde{\mathbf{F}}}$ are the projection matrices that estimate the true projections $\mathbf{P}_{\bm\beta}$ and $\mathbf{P}_{\mathbf{F}}$, respectively. So, $ \widehat{\mathbf{M}}$ reduces to the “intrinsic dimension" of the parameters by projecting the data matrix onto the low-dimensional subspaces, which are respectively spanned by $\widetilde{\bm\beta}$ and $\widetilde \mathbf{F}$. Usually, this is accomplished by employing PCA, where columns of $\widetilde{\bm\beta}$ and $\widetilde \mathbf{F}$ respectively denote the top left and right singular vectors of $\mathbf{Y}$. However, one major limitation of PCA-based low-rank projection is that it requires $\mathrm{rank}(\widetilde{\bm\beta})$ and $\mathrm{rank}(\widetilde{\mathbf{F}})$ be either equal to the true rank or a consistent estimator for it, which is often a strong assumption. As we commented in the introduction, the consistent estimation of the true rank is threatened by the strength of factors, raising a severe concern in practical applications.
In the context of pure factor model, fan2022learning proposed diversified projection (DP), as an alternative low-rank projection. This method begins by specifying two weighting matrices, an $N \times R$ matrix $\mathbf{W}_{\bm\beta}$ and a $T \times R$ matrix $\mathbf{W}_\mathbf{F}$, and defining:
where the weighting matrices consist of “diversified elements”, but not necessarily eigenvectors. These weighting matrices should satisfy:
(i) they are uncorrelated with the noise $\mathbf{E}$.
(ii) they are correlated with the actual $\bm\beta$ and $\mathbf{F}$.
Instead of the correct estimation of the number of factors, it is only required that the rank of the weighting matrix, $R$, be no smaller than the true number of factors, $r$.
Section (ref) will give specific recommendations for choosing the weighting matrices in applications. One of the appealing recommendations is to use transformations of initial observations, which does not require extra data. Moreover, we also note that the use of extra information is not uncommon in a wide range of literature. To name a few, fan2016projected and kelly2020instrumented use firms' characteristics to study excess returns. In the matrix completion literature, the side information is widely used jain2013provable,xu2013speedup, chiang2015matrix,wang2018high.
We utilize the diversified projections in the general form of low-rank inference problems. To begin with, we pre-determine weighting matrices, $\mathbf{W}_{\bm\beta}$ and $\mathbf{W}_{\mathbf{F}}$, whose rank $R$ is pre-determined, but not necessarily equal to the true rank $r$. The theory holds as long as $R\geq r$. To account for heterogeneity, let $\widehat{\bf\Pi}=\mathrm{diag}(\widehat{p}_1, \ldots, \widehat{p}_N)$ where $\widehat{p}_i=T^{-1} \sum_{t=1}^T X_{it}^2$, and $\widehat{\bf\Psi}=N^{-1}\mathrm{diag}(\sum_{j=1}^N X^2_{j1}\widehat{p}_j^{-2}, \ldots, \sum_{j=1}^N X^2_{jT}\widehat{p}_j^{-2})$.
Our Steps 1-3 are in spirit similar to that of the existing procedure outlined in the introduction. In particular, the bias $\mathcal{B}_1$ in Step 2 is similar to the one derived by chen:2019inference and xia2021statistical, who reach the following decomposition after the debias:
But we also have three key differences: First, we do not need to specify the true rank or its consistent estimator. Secondly, the low-rank projections in Step 3, $\mathbf{P}_{\widetilde{\bm\beta}}$ and $\mathbf{P}_{\widetilde{\mathbf{F}}}$, are inconsistent. Due to the rank over-estimation, they are “larger” than the true low-rank projections, $\mathbf{P}_{\bm\beta}$ and $\mathbf{P}_\mathbf{F}$ (we will make this insight precise later). Finally, the inconsistency of the estimated projections gives rise to a new bias correction in Step 4, which is the main novel statistical insight of this paper. In the next subsection we explain the source of this new bias in details.
After Step 2, we end up with $ \widetilde{\mathbf{M}}^{\mathrm{naive}}$. Write it as
where $\bf\mathcal{Z}$ is the estimation error. Therefore, we have
The first term on the right hand side yields the asymptotic normality of the estimator. For now we will focus on the second term $\mathbf{P}_{\widetilde{\bm\beta}} \mathcal{Z} \mathbf{P}_{\widetilde{\mathbf{F}}}$. Recall that $R$ denotes the rank of $\mathbf{P}_{\widetilde{\bm\beta}}$ and $\mathbf{P}_{\widetilde{\mathbf{F}}}$. In the usual case $R=r$, this term is asymptotically negligible. But when $R>r,$ $\mathrm{span}(\widetilde{\bm\beta})$ and $\mathrm{span}(\widetilde{\mathbf{F}})$ would also encompass additional noise components that are orthogonal to $\mathrm{span}(\bm\beta)$ and $\mathrm{span}(\mathbf{F})$, but are strongly correlated with the noise $\mathcal{Z} $. This strong correlation renders the second term asymptotically non-stochastic, introducing a bias through an implicit regularization, as we explain in detail below.
To understand the intuition, note that the diversified projection yields an $R \times r$ rotation matrix $\mathbf{H}$ such that
In the usual setting when $R=r$, the rotation matrix is well invertible (all its singular values are bounded away from zero), but this is no longer the case when $R>r$. In this case, the following rank-$r$ matrix becomes degenerate:
On the other hand, define $$\bm{S}_{\widetilde{\mathbf{F}}}=\frac{1}{T} \widetilde{\mathbf{F}}' \widetilde{\mathbf{F}}, \quad R \times R.$$ We shall show that when $R>r$, $\bm{S}_{\widetilde{\mathbf{F}}}$ is still asymptotically invertible, but its eigenvalues may decay very fast. Nevertheless, the projection matrix $\mathbf{P}_{\widetilde \mathbf{F}}=\widetilde \mathbf{F}(\widetilde \mathbf{F}'\widetilde \mathbf{F})^{-1}\widetilde \mathbf{F}'$ is still well defined with probability approaching one.
The asymptotic property of the projection matrix critically depends on the following ridge-type projection function, also known as Tikhonov-regularization function:
In fact, there is a rate $x_{NT} \rightarrow 0$ such that, $\bm{S}_{\widetilde{\mathbf{F}}}^{-1} \approx f(x_{NT})$. The key challenge, however, is that $f(x)$ is discontinuous at $x=0$ and $\lim_{x \rightarrow 0}f(x)$ does not exist. Therefore when $R>r$, the inverse matrix $\bm{S}_{\widetilde{\mathbf{F}}}^{-1}$ does not converge in probability to the generalized inverse $\bm{S}_\mathbf{F}^+$, its population counterpart.
The discontinuity challenge can be avoided in our context by considering the following rescaled Tikhonov-regularization function:
Unlike $f(x)$, the rescaled Tikhonov function is continuous in neighborhoods of zero and has $(\frac{1}{T}\mathbf{F}'\mathbf{F})^{-1}$ as its limit when $x \rightarrow 0$. Fortunately, in the low-rank inference problem, it is sufficient to study the behavior of $\widetilde{f}(x) $ instead of $f(x)$ because the projection matrix $\mathbf{P}_{\widetilde{\mathbf{F}}}$ asymptotically depends on $\bm{S}_{\widetilde{\mathbf{F}}}$ through $\mathbf{H}'\bm{S}_{\widetilde{\mathbf{F}}}^{-1}\mathbf{H}$. Let $\bm{S}_\mathbf{F}^{+}$ denote the generalized inverse of $\bm{S}_\mathbf{F}$. We shall show that while $\|\bm{S}^{-1}_{\widetilde{\mathbf{F}}}-\bm{S}^+_{\mathbf{F}}\| \neq o_P(1),$ when rescaled by $\mathbf{H}$, we have,
for some sequence $x_{NT}\rightarrow 0$.
Above all, when $R>r$, the correlation between the estimated low-rank projection matrices and $ \mathcal{Z} $ can be characterized by a Tikhonov-type regularization function, akin to ridge regression, which acts as an implicit regularization. Therefore, a key new statistical insight of this paper is that, unlike the conventional understanding where overfitting primarily results in additional variance, in the context of low-rank inference, over-estimating the rank introduces a novel source of asymptotic bias through $\mathbf{P}_{\widetilde{\bm\beta}} \mathcal{Z} \mathbf{P}_{\widetilde{\mathbf{F}}}$. We identify and address this bias and reach the final estimator:
This is the motivation for introducing $\mathcal{B}_2$ in Step 4.
The objective is to establish the asymptotic normality of our estimator for the group average, given by $|\mathcal{G}|^{-1} \sum_{(i,t) \in \mathcal{G}} \widehat{M}_{it}$, where $\mathcal{G} \subset \{1, \ldots, N\} \times \{1, \ldots, T\}$ represents the group we are interested in.
The following assumption formalizes the data generating process (DGP) of the noise $\mathbf{E}$.
We note that Assumption (ref) (ii) accommodates the noisy matrix completion problem where $X_{it}$ indicates if the entry $(i,t)$ is observed.
The next assumption specifies the DGP of the regressor matrix $\mathbf{X}$. We allow heterogeneity in $\mathbf{X}$ across units. In the matrix completion application, it can accommodate the heterogeneous missing probabilities across $i$. Also, we allow cross-sectional weak dependence in $\mathbf{X}$ through the cluster structure, where the size of the largest cluster is allowed to grow. Let $\mathcal{C}_1,\ldots, \mathcal{C}_{\rho} \subset \{1, \ldots, N\}$ be non-empty and disjoint clusters such that $\cup_{g=1}^\rho \mathcal{C}_g=\{1, \ldots, N\}$. We denote $\vartheta \coloneqq \max_{g=1, \ldots, \rho}|\mathcal{C}_g|$.
Next, the following assumption specifies the class of high-dimensional matrices we are interested in. We require $\mathbf{M}^{\star}$ to be of low-rank and incoherent.
The following assumption specifies conditions on the diversified weights and the strength of factors. Regarding the factor strength, without loss of generality, we fix $\sigma_{\min}(\mathbf{F}) \asymp \sigma_{\max}(\mathbf{F}) \asymp \sqrt{T}$ and accommodate the weak factor by allowing the aggregated factor loadings to be “small."
Our condition on the factor strength is relatively weak. For example, suppose that there exist a bounded sequence $a_N=O(1)$, an $N \times r$ matrix ${\bm\beta}_0$, and diversified weights $\mathbf{W}_{\bm\beta}$ such that:
for some constant $C>0$. Here, we assume $\sigma_{\min}(N^{-1}{\bm\beta}_0'{\bm\beta}_0) \asymp \sigma_{\max}(N^{-1}{\bm\beta}_0'{\bm\beta}_0)\asymp C$ for some constant $C>0$, so that ${\bm\beta}_0$ can be regarded as the “standardized direction” of ${\bm\beta}$. Therefore, the strength of ${\bm\beta}$ is governed by the sequence $a_N$. Then, Assumption (ref) (ii) and (iii) can be simplified to: $ \log N \ll \sqrt{N} a_N. $ This implies the condition on the factor strength as: $$\sigma_{\min}\left( {\bm\beta}'{\bm\beta}\right) \asymp Na_N^2 \gg \log^2 N ,$$ which is only slightly stronger than the definition of “weak factors” in onatski2012asymptotics in the pure factor model. In contrast, the usual requirement for inference in the PCA-based approaches bai2023approximate corresponds to $\sigma_{\min}\left( {\bm\beta}'{\bm\beta}\right)\gtrsim N^{c}$ for some $c \in (1/2,1].$
We now introduce notations for groups and assumptions for them. This paper explores three distinct types of group averages of $\mathbf{M}^{\star}$ for inference: i) block averages, ii) serial averages, and iii) cross-sectional averages. For the block, denoted as $\mathcal{G}_{\mathrm{bl}}$, we define $\mathcal{G}_{\mathrm{bl}} = \mathcal{I} \times \mathcal{T}$ where $\mathcal{I} \subset \{1, \ldots, N\}$ and $\mathcal{T} \subset \{1, \ldots, T\}$. Similarly, for cross-sectional groups, denoted as $\mathcal{G}_{\mathrm{cs}}$, let $\mathcal{G}_{\mathrm{cs}} = \{1, \ldots, N\} \times \mathcal{T}$ where $\mathcal{T} \subset \{1, \ldots, T\}$, and for serial groups, denoted as $\mathcal{G}_{\mathrm{serial}}$, let $\mathcal{G}_{\mathrm{serial}} = \mathcal{I} \times \{1, \ldots, T\}$ where $\mathcal{I} \subset \{1, \ldots, N\}$. For brevity, we present our main results for the block average case only. The results for other cases, which are similar, are provided in the appendix.
Finally, it is critical for the initial estimator $\widetilde{\mathbf{M}}^{\mathrm{init}}$ to have desirable properties, such as reasonable convergence rates and weak correlation with the noise. We provide the construction of such $\widetilde{\mathbf{M}}^{\mathrm{init}}$, that is based on the nuclear norm penalized estimation, in Section (ref).
Define the covariance estimator:
where $\widehat{\mathbf{M}}_{\mathcal{I}, \cdot}\coloneqq |\mathcal{I}|^{-1} \sum_{i \in \mathcal{I}}\widehat{\mathbf{M}}_{i, \cdot}$, and $\widehat{\mathbf{M}}_{\cdot, \mathcal{T}}\coloneqq|\mathcal{T}|^{-1} \sum_{t \in \mathcal{T}}\widehat{\mathbf{M}}_{\cdot,t}$.\footnote{For a matrix $\mathbf{A}$, let $\mathbf{A}_{j,\cdot}$ and $\mathbf{A}_{\cdot, j}$ denote $j$th row and $j$th column of $\mathbf{A}$, respectively.}
The asymptotic variance depends on the choice of the DP weighting matrices. In the special case that $R$ is a consistent estimator of the true rank, these weights can be chosen optimally as the usual eigenvectors. Then the asymptotic variance attains the efficiency bound as achieved by chernozhukov2023inference,choi2023inference. But more general weighting matrices will lead to efficiency loss, which is the cost of being flexible of choosing these weighting matrices to be robust to over-estimating the rank.
We now formally characterize the initial estimator $\widetilde{\mathbf{M}}^{\mathrm{init}}$ in Step 1. The central limit theorem arises from the entries of $\mathbf{P}_{\widetilde{\bm\beta}} {\mathbf{M}^{\star}} \mathbf{P}_{\widetilde{\mathbf{F}}}-{\mathbf{M}^{\star}}$, which can be shown as weighted row and column sums of the noise matrix $\mathbf{E}$. Not so surprisingly, however, the weights in these sums are correlated with $(\widetilde{\bm\beta}'\widetilde{\bm\beta})^{-1}$ and $(\widetilde{\mathbf{F}}'\widetilde{\mathbf{F}})^{-1}$ which, by construction, depend on the initial estimator. To establish the CLT, a technical challenge arises from the correlation between $\widetilde{\mathbf{M}}^{\mathrm{init}}$ and the rows and the columns of $\mathbf{E}$. Hence the initial estimator should be constructed in a way such that the correlation can be well controlled.
We adopt a standard approach that artificially creates independence that is inspired by the sample splitting idea. Specifically, $\widetilde{\mathbf{M}}^{\mathrm{init}}$ is constructed through two steps, the first step with the full sample and the second step with restricted samples.
\noindentFull sample. We define
where $\lambda>0$ is a tuning parameter. The objective function incorporates inverse weights, $\widehat{p}^{-1}_j=(T^{-1} \sum_{t=1}^T X_{jt}^2)^{-1}$ for $j=1, \ldots, N$, to accommodate the heterogeneity in $\mathbf{X}$. Similar weighting techniques have been employed in previous works such as ma2019missing and choi2023inference.
\noindentRestricted sample. Next, to remove the effect of correlations between $\widetilde{\mathbf{M}}^{\mathrm{full}}$ and entries in $\mathbf{E}$, we then replace the majority of entries in $\widetilde{\mathbf{M}}^{\mathrm{full}}$ with the estimates that are independent of the noises appearing in the weighted sum. This section will focus on the block group $\mathcal{G}_{\mathrm{bl}}$, where we are interested in making inference for the group average over $i\in \mathcal I$ and $t\in\mathcal T$: $$ \frac{1}{|\mathcal{T}||\mathcal{I}|}\sum_{i\in\mathcal{I}}\sum_{t\in\mathcal{T}} M^{\star}_{it}, $$ under Assumption (ref).\footnote{The construction of $\widetilde{\mathbf{M}}^{\mathrm{init}}$ slightly differs for other group types. Refer to the appendix for other types.} Compute the nuclear norm penalized estimation using the sample outside of the group: $i\notin \mathcal I$ and $t\notin\mathcal T$,
\noindentMerging two estimates. Then, we define the initial estimator $\widetilde{\mathbf{M}}^{\mathrm{init}}$ as
In the block group case, the weighted sum of noises, that leads to asymptotic normality, will consist of rows $i\in\mathcal{I}$ and columns $t\in\mathcal{T}$ of $\mathbf{E}$. By construction, the majority of entries ($i\notin \mathcal I$ and $t\notin\mathcal T$) in $\widetilde{\mathbf{M}}^{\mathrm{init}}$ are independent of these noises, as intended.
We present two specific statistical applications: treatment effect estimation and multiple testing. In both applications the true rank is typically unknown and it is critical to develop inferential methods that are robust to over-specifying the rank.
This section elaborates on the application of our theory to heterogeneous treatment effect estimation and presents formal asymptotic results. Following the causal inference literature, e.g., rubin:1974, imbens:2015, we assume that, for each $(i,t),$ there exist two potential outcomes, $Z^{(0)}_{it}$ and $Z^{(1)}_{it}$, where $Z^{(0)}_{it}$ represents the outcome that would be observed if $(i,t)$ is controlled and $Z^{(1)}_{it}$ represents the outcome that would be observed when $(i,t)$ is treated. We observe only one of the potential outcomes for each $(i,t)$, which basically defines two incomplete matrices.
We assume that the two potential outcomes have the following structure:
for each $\iota \in \{0,1\}$. By defining $D_{it}=\mathbf{1}\{\text{$(i,t)$ is treated}\}$, $\mathbf{X}^{(0)}=[1-D_{it}]_{i \leq N, t \leq T}$, and $\mathbf{X}^{(1)}=[D_{it}]_{i \leq N, t \leq T}$, we can represent the two sets of observed data in the following way:
for each $\iota \in \{0,1\}$. Therefore, by applying the matrix completion method to each of $\mathbf{Y}^{(0)}$ and $\mathbf{Y}^{(1)}$, we will obtain $\widehat{\mathbf{M}}^{(0)}$ and $\widehat{\mathbf{M}}^{(1)}$ that estimate $\mathbf{M}^{(0)}$ and $\mathbf{M}^{(1)}$ respectively.
Our goal is to perform inference about group average treatment effects for group $\mathcal{G}$. We denote the treatment effect for each $(i,t)$ as $\Gamma_{it}= M^{(1)}_{it}-M^{(0)}_{it}$. Then, the group average treatment effect is defined as
Then, a natural choice for the estimator for the group treatment effect would be
Assuming that the assumptions in Section (ref) hold for each superscript $(0)$ and $(1)$, we can establish the asymptotic normality for $|\mathcal{G}|^{-1} \sum_{(i,t) \in \mathcal{G}} \widehat{\Gamma}_{it}- |\mathcal{G}|^{-1} \sum_{(i,t) \in \mathcal{G}} \Gamma_{it}$. As in Section (ref), we provide the result for the block average (the heterogeneous treatment effect). The results for other group types are provided in the appendix.
We define, for each $\iota=0,1,$
with $\widehat{\mathbf{M}}_{\mathcal{I}, \cdot}^{(\iota)}\coloneqq |\mathcal{I}|^{-1} \sum_{i \in \mathcal{I}}\widehat{\mathbf{M}}_{i, \cdot}^{(\iota)}$, and $\widehat{\mathbf{M}}^{(\iota)}_{\cdot, \mathcal{T}}\coloneqq|\mathcal{T}|^{-1} \sum_{t \in \mathcal{T}}\widehat{\mathbf{M}}^{(\iota)}_{\cdot,t}$.
In large-scale multiple testing, it is crucial to address the well-known confounding factors and the resulting strong correlations leek2008general, friguet2009factor, wang2017confounder,fan2019farmtest. We consider the model, for $t=1,...,T$:
Here ${\bm\mu}= (\mu_1, \ldots, \mu_N)'$ is the mean vector and $\mathbb E \mathbf{U}_t=0$. The noise $ \mathbf{U}_t$ consists of the independent part $\mathbf{E}^{\star}_t$ and the confounding factor part ${\bm\beta \mathbf{F}_t} $.
We consider a practical situation where data is not fully observable. Let $\mathbf{X}_t$ be an $N$-dimensional vector, whose element $X_{it}=\mathbf{1}\{Z_{it} \text{ is observed}\}$. Then, the observed data is $\mathbf{Y}_t:=\mathbf{X}_t \circ \mathbf{Z}_t$, satisfying: $$ \mathbf{Y}_t= \mathbf{X}_t \circ ({\bm\mu} + \mathbf{U}_t). $$ The objective is to test $N$ hypotheses:
This model differs from the usual multiple testing model in two ways: (i) the noise in $\mathbf{U}_t$ are strongly dependent due to the presence of confounding factors $\bm\beta \mathbf{F}_t$; and (ii) the “data" is observed subjected to missing values, indicated by the binary vector $\mathbf{X}_t$.
While the importance of addressing confounding correlations in multiple testing has been widely recognized wang2017confounder, a critical yet unresolved question remains: How much confounding correlation should be accounted for? This question raises two concerns, both of which are highly relevant to practical applications.
First, researchers have relied on consistent estimation of the rank in $\mathbf{U}_t$, though ensuring its accuracy has always been a concern. Secondly, methods that explicitly allow $\bm\beta \mathbf{F}_t$ inherently assume the presence of at least one confounding factor. But what happens if we account for confounding factors when, in reality, there are none? This corresponds to a very special case of over-estimating the rank, where $r=0$ but $R>0$. The fact that whether $r$ equals zero is unknown in practice, As a precaution, statisticians often account for $R>0$ “factors" regardless. In this subsection, we prove that the results are uniformly valid for all cases where $R \geq r$, even when $r=0$. This finding is empirically significant: one should always account for potential confounding correlations, as doing so does no harm even when none exist, at least asymptotically. \footnote{An alternative practice is to pretest whether $\bm\beta \mathbf{F}_t$ exists. But the power of such tests are inherently affected by the strength of the factors, which may not be detectable if factors are weak in finite sample. }
We define $\bar{Y}_i=|\mathcal{T}_i|^{-1} \sum_{t \in \mathcal{T}_i} Y_{it}$ where $\mathcal{T}_i$ is the set of observed indices in the $i$th row. We also define the demeaned data $Y^{\mathrm{d}}_{it}= Y_{it}- \bar{Y}_i$ if $X_{it}=1$, and $Y^{\mathrm{d}}_{it}=0$ otherwise. Let $(\mathbf{Y}^{\mathrm{d}},\mathbf{X}, \mathbf{E}^{\star}) $ be the matrices of $(Y_{it}^{\mathrm{d}}, X_{it}, E_{it}^{\star})$. Then $ \mathbf{Y}^{\mathrm{d}}\approx \mathbf{X}\circ\mathbf{M}^{\star} +\mathbf{E} $, where $\mathbf{M}^{\star} $ denotes the matrix of $\bm\beta \mathbf{F}_t$ and $\mathbf{E} \coloneqq \mathbf{X}\circ \mathbf{E}^{\star}$. We implement Algorithm (ref) on $\mathbf{Y}^{\mathrm{d}}$ to obtain $\widehat{\mathbf{M}}.$\footnote{The initialization in Section (ref), though necessary for the CLT for $\mathbf{M}^{\star}$, is not required for testing $\bm\mu$. Instead, we simply implement the full-sample estimation (ref) on $\mathbf{Y}^{\mathrm{d}}$ just once and use it as $\widetilde{\mathbf{M}}^{\mathrm{init}}$.} Our proposed estimator for ${\bm\mu}$ is
The above expansion (ref) shows that the leading terms in the expansion of $\widehat{\mu}_i-\mu_i$ are cross-sectionally weakly correlated, e.g., the confounding correlations have been successfully removed, regardless of whether the confounding factors are present. From here, one can apply the standard multiple testing approach, such as benjamini1995controlling (B-H procedure). The standard B-H procedure implements the test as follows: let $p_{(1)}\leq...\leq p_{(N)}$ denote the sorted p-values for each test. Then $H_0^i$ is rejected if $p_i\leq p_{(k)}$ where $k=\max\{i\leq N: p_{(i)}\leq \tau i/N\}$. Building on ((ref)), liu2014phase showed that the false discovery rate of the B-H procedure can be controlled below $\tau$ asymptotically.
In this section, several choices of the diversified weighting matrices are proposed. We emphasize that the requirement for the constructed weights is quite mild: they do not need to consistently estimate the true parameters. Instead, it suffices that they are informative with respect to the underlying matrix parameter.
Suppose $\beta_{ik}= g_k(\mathbf{b}_i,\eta_{ik})$ where $g_k(\cdot)$ are unknown functions, $\mathbf{b}_i$ is a vector of observable characteristics, and $\eta_{ik}$ is noise. We define $\mathbf{W}_{\bm\beta}$ as the transformations of $\mathbf{b}_i$, i.e,. $W_{{\bm\beta}, ik}=\phi_k(\mathbf{b}_i)$ with transformation functions $\phi_k(\cdot)$. Similarly, we assume that there are observable characteristics $\mathbf{f}_t$ for factors. Namely, $F_{tk}=h_k(\mathbf{f}_t, \xi_{tk})$, where $\xi_{tk}$ is noise. We then define $W_{\mathbf{F},tk}=\varphi_k(\mathbf{f}_t)$ with a set of transformation functions $\varphi_k(\cdot).$
It is not uncommon to observe individual-specific characteristics, which bring additional information regarding the singular vectors. For example, the “Netflix Challenge” called for a matrix completion problem bennett2007netflix, where $\mathbf{b}_i$ denotes customers' demographic characteristics such as age and sex, which might be related to their preferences. In addition, films are classified according to their genres: action, romance, sci-fi, and drama, which are denoted by $\mathbf{f}_t$. For example, harper2015movielens provide user ratings of movies along with characteristics of users and movies.\footnote{Additionally, Yale University's library has documented over 40 film genres, styles, categories, and series in its Film Studies Research Guide. For further reference, see \url{https://guides.library.yale.edu/c.php?g=295800&p=1975072}.} In financial applications such as asset pricing, $\mathbf{b}_i$ may consist of firm characteristics as numerous studies connor2012efficient, fan2016projected. Also, observed macroeconomic factors or Fama-French factors fama1993common can be used as $\mathbf{f}_t$.
Suppose we have access to additional data $(Y_{it},X_{it})$ for $i\in \mathcal I_1$ and $t\in\mathcal T_1$, where $\mathcal I_1$ and $\mathcal T_1$ are the index sets for extra data.
Consider the subsample on $t\in\mathcal T_1$. We can write the model as $ \mathbf{Y}_{t} = \mathbf{X}_{t} \circ ({\bm\beta} \mathbf{F}_t) + \text{noise}, $ for $t\in\mathcal T_1$, where $\mathbf{Y}_t$ and $\mathbf{X}_t$ are the vectors of the same subjects but observed at the extra time $t\in \mathcal T_1$. This provides extra information about the factor loading ${\bm\beta}$, and the size of $\mathcal T_1$ is sufficient as long as it is larger than $R$. The diversified weight $\mathbf{W}_{\bm\beta}$ can be constructed based on: $$ \bar{b}_i= \sum_{t\in\mathcal T_1} X_{it}Y_{it}\slash \sum_{t\in\mathcal T_1} X_{it}^2. $$ We then construct $W_{{\bm\beta},ik}= \phi_{k}(\bar{b}_i)$ for $i \leq N$ and $k \leq R$ using the nonlinear transformations of $\bar{b}_i$ to span $R$-dimensional space. Similarly, we can construct $\mathbf{W}_\mathbf{F}$ using the additional sample $i\in\mathcal I_1$.
If neither characteristics nor the extra sample are available, we can still implement DP via transformations of the initial observations. Appealingly, it does not require extra data.
We suppose there are “initial observations" to satisfy the following conditions:
These conditions state that the idiosyncratic noise of the initial period and of some individual are independent of the rest. Then we can use transformations of first observation and the observation of $i=1$ as the diversified weights: let $$ \mathbf{b}= (Y_{1,t_0},...,Y_{N,t_0})',\quad N\times 1; \quad \mathbf{Y}_1= (Y_{1,1},...,Y_{1,T})',\quad T\times 1, \quad \text{and} $$ $$ \mathbf{W}_{\bm\beta}= (\phi_1(\mathbf{b}),...,\phi_R(\mathbf{b}));\quad \mathbf{W}_\mathbf{F}= (\psi_1(\mathbf{Y}_1),...,\psi_R(\mathbf{Y}_1)) $$ where $\{\phi_k:k\leq R\}$ and $\{\psi_k:k\leq R\}$ are the sets of transformation functions. Then apply it to data except for $t=t_0$ and $i=1$. The cost would be the loss of $N+T$ observations in the panel, which is mild.
We evaluate the finite sample performance of our estimator in matrix completion design. In Section B of the supplement, we also examine the performance in other two designs: varying coefficient model and heterogeneous treatment effect.
The DGP for the matrix completion design is as follows: $$\mathbf{Y}= \mathbf{X} \circ ({\bm\beta} \mathbf{F}' + \mathbf{E}^{\star} ) = \mathbf{X} \circ ({\bm\beta} \mathbf{F}')+\underbrace{\mathbf{X}\circ\mathbf{E}^{\star}}_{\coloneqq \mathbf{E}}$$ where $\beta_{ik} =(2*\cos^k(b_i)+0.5*u_{ik})* N^{-(1-\alpha)/2}$, and $F_{tk}=2*\cos^k(f_t) +0.5*\nu_{tk}$ with $b_i$, $f_t$, $u_{ik}$, and $\nu_{tk}$ drawn from $\mathcal{N}(0,1)$ independently. The noise $E^{\star}_{it}$ is generated from $\mathcal{N}(0,1)$ independently. We note that the constant $\alpha \in (0,1]$ determines the strength of factors. A larger (smaller) $\alpha$ implies stronger (weaker) factors. $\mathbf{X}$ is the binary matrix where $X_{it} \sim \mathrm{Bernoulli}(p_i)$ with $p_i$ drawn from $\mathrm{Unif}[0.5, 0.8].$ Throughout simulations, we fix the true rank at $r=2$ and conduct $1,000$ replications.
We construct the diversified weights, $\mathbf{W}_{\bm\beta}$ and $\mathbf{W}_\mathbf{F}$ with $R=4$, following Section (ref), which are (1) Observed characteristics, (2) Extra sample averages, and (3) Initial transformation. See Section B in the supplement for detailed implementation of each choice, the choice of the tuning parameter, and experiments with other values of $R$. We are interested in estimating three types of averages of the low-rank matrix: (I) “Block," where $\mathcal{G}=\mathcal{I} \times \mathcal{T}$ and $|\mathcal{I}|=|\mathcal{T}|=5$; (II) “CS," where $\mathcal{G}=\{1, \ldots, N\} \times \{T\}$, and (III) “Serial," where $\mathcal{G}=\{N\} \times \{1, \ldots, T\}$.
Table (ref) displays coverage probabilities of the confidence intervals, with the target probability 0.95. The weak factor setting ($\alpha=0.5$) is particularly relevant to our theory as the direct rank estimation fails in this DGP.\footnote{The rank of the nuclear norm penalized estimator depends on the tuning parameter $\lambda.$ For the choice of $\lambda$, we follow the methods from chernozhukov2018inference,chernozhukov2023inference, choi2023inference. The choice of $\lambda$ in these studies involves a choice of a small constant $c>0.$ chernozhukov2018inference,chernozhukov2023inference set $c=1/10$ while choi2023inference set $c=1/7.$ With both values, the rank of the nuclear norm penalized estimator is almost always one, whereas the true rank $r$ is 2.} The coverage probabilities are reasonably good especially when estimating the “CS" and “Serial" averages, although we do observe size distortions when estimating the “Block" average.
To compare among choices of the weighting matrices, we fix $N=T=200$ and $\alpha=0.5$ (weak factors) and present the histograms of the standardized estimates, with the standard normal density, in Figure (ref). For comparison, we include the nuclear norm penalized estimator (“Nucl") whose standard error is simulation-based. The standard error is estimated for all other methods. It is evident that “Nucl" suffers from bias, and our estimators approximate the standard normal distribution well for all group types and diversified weights.
In the supplement we also compare the MSE of these methods. Not surprisingly, “Observed" and “Extra" yield smaller MSE as they use extra information to specify the DP weights. But appealingly, the initial transformation, which does not require extra information, performs reasonably well.
We apply the inferential theory presented in Section (ref) to investigate the influence of the U.S. presidential election on the federal grant allocation to states. Although the specific allocation of federal funds is carried out by Congress, the U.S. president also wields notable influence in the grant allocation process. The role of the president in the allocation involves proposing annual federal budget proposals to Congress, signing or vetoing bills, and supervising the executive agencies. In the allocation process, the presidents may seek to allocate more federal funds to the states that supported them in the election, for political incentives. This practice is commonly referred to as “pork-barrel politics.” For decades, this practice has been investigated in numerous studies both theoretically and empirically larcinese2006allocating,berry2010president. By employing the treatment effect estimation application presented in Section (ref), we aim to test whether the pork-barrel politics exist or not, in the history of the U.S. presidential elections.
We use the data of the U.S. federal grants, which cover fiscal years from 1953 to 2021 and include 50 U.S. states in addition to the District of Columbia.\footnote{In these data sets, the years 1960, 1972, and 1977 to 1980 are missing. As a result, our analysis does not include the Carter administration.} These data are publicly accessible on the websites of the U.S. Census Bureau, the National Associate of State Budget Officers (NASBO), and the Social Security Administration (SSA).
Following the notation in Section (ref), we define state $i$ is “treated" if it supported the incumbent president of year $t$ in the preceding election. We note that the treatment assignments in this case could potentially be endogenous. Nevertheless, we presume that they are randomly assigned treatments and proceed to apply our approach. We calculate the per-capital federal grant, denoted as $Z_{it}$, for each state-year pair $(i,t)$. In order to detrend the data, we define $Y^{(\iota)}_{it}=X^{(\iota)}_{it} \times (Z_{it}/\sum_{i=1}^N Z_{it}) \times 100$ for each $\iota=0,1,$ and assume that $\mathbf{Y}^{(\iota)}$ follows the model (ref).
Before proceeding to tests, we present the singular values of the nuclear norm penalized estimators for the treated and control samples (Figure (ref)). The singular value plots do not provide a very definitive understanding of the true rank. While several data-driven methods choose $r=2$ as presented in Section A.3, we find that it does not lead to optimal out-of-sample performance in Section A.2. This motivates our method which does not rely on estimating the rank.
We construct the diversified weighting matrices following the observed characteristics approach explained in Section (ref). To be specific, the columns of $\mathbf{W}_\mathbf{F}$ and $\mathbf{W}_{\bm\beta}$ consist of polynomial transformations (up to the second power) of the annual data of U.S. GDP growth rates and unemployment rates, and the state-by-state data on the averages of the population and the annual per-capita personal income from 1953 to 2021.
To begin with, we estimate individual state effects, i.e., the overall time average treatment effects for each state. The individual state effects and their t-statistics are presented in Figure (ref). For each state, the null hypothesis is that there is no pork-barrel politics, i.e., the average treatment effect is non-positive. We reject the null hypothesis in 30 states at the 5% significance level and in 20 states at the 1% significance level. Figure (ref) depicts the distribution of the states where the null hypotheses are rejected. These test results suggest that the pork-barrel politics exist in a larger number of states. We now turn to a very natural question: Why are some states enjoying the “pork,” while others are not?
We observe that many of the states colored in dark green in Figure (ref) are known as “loyal” states in that they have supported one party over decades. For example, DC has exclusively supported the Democratic Party since 1964, while AK, ID, OK, SD, UT, and WY have consistently supported the Republican Party since 1968. Inspired by this observation, we classify all 51 states based on their “loyalty” to a particular party. Specifically, we count the number of times that a state switches the party it supports, referred to as a “swing,” since the 1952 U.S. presidential election, in Table (ref).
Based on Table (ref), we estimate the overall time average treatment effects for the state groups. The first plot in Figure (ref) indicates a positive relation between the loyalty and treatment effects: stronger loyalty leads to more substantial treatment effects. This “rewarding-loyalty” pattern is even clearer with t-statistics in the second plot. Also, these t-statistics show that the pork-barrel politics exist in almost all groups. Only the t-statistic of the swing states is slightly lower than 2.33, and all other t-statistics are much larger. \footnote{In the appendix, we provide additional test results for other group averages: group averages for each Party governance and each presidential administration.}
\onehalfspacing