The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
76,852 characters
Inference for Low-rank Models without Estimating the Rank
\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1}
\if00
{
\title{\bf Inference for Low-rank Models without Estimating the Rank}
\author[1]{Jungjun Choi}
\author[2]{Hyukjun Kwon}
\author[3]{Yuan Liao}
\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}
\maketitle
} \fi
\if10
{
\bigskip
\bigskip
\bigskip
\begin{center}
\spacingset{1.3}
{\LARGE\bf Inference for Low-rank Models without Estimating the Rank}
\end{center}
\medskip
} \fi
\bigskip
\begin{abstract}
This paper studies the inference about linear functionals of high-dimensional low-rank matrices. While most existing inference methods would require consistent estimation of the true rank, our procedure is robust to rank misspecification, making it a promising approach in applications where rank estimation can be unreliable. We estimate the low-rank spaces using pre-specified weighting matrices, known as diversified projections. A novel statistical insight is that, unlike the usual statistical wisdom that overfitting mainly introduces additional variances, the over-estimated low-rank space also gives rise to a non-negligible bias due to an implicit ridge-type regularization. We develop a new inference procedure and show that the central limit theorem holds as long as the pre-specified rank is no smaller than the true rank. In one of our applications, we study multiple testing with incomplete data in the presence of confounding factors and show that our method remains valid as long as the number of controlled confounding factors is at least as large as the true number, even when no confounding factors are present.
\end{abstract}
\noindent
{\it Keywords:} High-dimensional models; low-rank matrix estimation; rank misspecification; weak factors
\vfill
\newpage
\spacingset{1.9}
\section{Introduction}\label{sec:intro}
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:
\begin{align}
{\bf Y=X \circ M^{\star}+E} \label{eq:model}
\end{align}
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 \cite{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 \citep[][]{chen:2019inference,xia2021statistical}:
\begin{enumerate}
\item[Step 0:] Fix the \textit{true} rank of ${\bf M^{\star}}$ (or consistently estimate the true rank).
\item[Step 1:] Obtain an initial estimator ${\bf \widetilde{M}^{\mathrm{init}}}$.
\item[Step 2:] Subtract a bias correction term: ${\bf \widetilde{M}}^{\mathrm{naive}}= {\bf\widetilde{M}}^{\mathrm{init}}- \mathcal{B}_1.$
\item[Step 3:] Project ${\bf\widetilde{M}}^{\mathrm{naive}}$ onto low-rank spaces that are consistently estimated by using PCA with the ``correct'' knowledge of the true rank.
\end{enumerate}
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 \citep[e.g.,][]{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 \textit{diversified projection} (DP) which was recently proposed by \cite{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 \textit{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 \cite{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 \textit{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, \cite{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 \cite{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]$ \citep[see][]{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{sec:choiceofW}. One of the appealing recommended choices is to use transformations of initial observations, which \text{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 \cite{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. \cite{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.
\section{Model and Estimation} \label{sec:modelandestimation}
In the linear model \eqref{eq:model}, we assume that the matrix $\mathbf{M}^{\star}$ has a low-rank structure:
\begin{align}
\mathbf{Y}=\mathbf{X} \circ {\bf M^{\star}}+\mathbf{E} = \mathbf{X} \circ ({\bm\beta} \mathbf{F}')+ \mathbf{E}\label{eq:model2}
\end{align}
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 \eqref{eq:model2} has numerous applications including varying coefficient models, heterogeneous treatment effects \citep{athey2021matrix}, and matrix completion problems.
\subsection{Diversified projection}
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:
\begin{align*}
\mathbf{Y}={\mathbf{M}^{\star}}+\mathbf{E}= {\bm\beta} \mathbf{F}' + \mathbf{E}.
\end{align*}
To estimate $\mathbf{M}^{\star}$, one approach is to take advantage of
the low-rank structure, by applying projections as follows:
\begin{align*}
\widehat{\mathbf{M}}=\mathbf{P}_{\widetilde{{\bm\beta}}} \mathbf{Y} \mathbf{P}_{\widetilde{\mathbf{F}}}.
\end{align*}
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, \cite{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:
\begin{align*}
\widetilde{\bm\beta}= \frac{1}{T} \mathbf{Y} \mathbf{W}_{\mathbf{F}} \quad \text{and} \quad \widetilde{\mathbf{F}}= \frac{1}{N}\mathbf{Y}' \mathbf{W}_{\bm\beta}
\end{align*}
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{sec:choiceofW} will give specific recommendations for choosing the weighting matrices in applications. One of the appealing recommendations is to use transformations of initial observations, which \text{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, \cite{fan2016projected} and \cite{kelly2020instrumented} use firms' characteristics to study excess returns. In the matrix completion literature, the side information is widely used \citep[e.g.,][]{jain2013provable,xu2013speedup, chiang2015matrix,wang2018high}.
\subsection{Formal procedure}\label{sec:formalprocedure}
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})$.
\begin{breakablealgorithm}
{\raggedright\textbf{\fname@algorithm~\thealgorithm} #\small Formal estimation procedure\par}
\ifx\relax#\relax\relax
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\small Formal estimation procedure}
\else
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\relax}
\fi
\kern2pt\hrule\kern2pt
\label{alg:estimation}
\begin{algorithmic}
\noindent \textbf{Step 1:} \textit{Initialization}. Obtain an initial estimator $\widetilde{\mathbf{M}}^{\mathrm{init}}$ as to be described in Section \ref{sec:initial}, and pre-determine the weighting matrices $\mathbf{W}_{\bm\beta}$ and $\mathbf{W}_\mathbf{F}$ with rank $R$. \\
{\textbf{Step 2:} \textit{First bias correction.} Let $\mathcal{B}_1=\widehat{\bf\Pi}^{-1} \mathbf{X} \circ (\mathbf{X}\circ \widetilde{\mathbf{M}}^{\mathrm{init}}-\mathbf{Y})$. Compute \vspace{-0.7cm}$$\widetilde{\mathbf{M}}^{\mathrm{naive}}=\widetilde{\mathbf{M}}^{\mathrm{init}} - \mathcal{B}_1.$$ } \\\vspace{-0.7cm}
{\textbf{Step 3:} \textit{Diversified projections.} Let $\widetilde{\bm\beta}= T^{-1} \widetilde{\mathbf{M}}^{\mathrm{naive}} \mathbf{W}_\mathbf{F}$ and $\widetilde{\mathbf{F}}= N^{-1} \widetilde{\mathbf{M}}^{\mathrm{naive} \prime} \mathbf{W}_{\bm\beta}.$}
\vspace{-0.7cm}
$$
\widetilde{\mathbf{M}}^{\mathrm{proj}}=\mathbf{P}_{\widetilde{\bm\beta}} \widetilde{\mathbf{M}}^{\mathrm{naive}}\mathbf{P}_{\widetilde{\mathbf{F}}}.
$$ \\
\vspace{-0.7cm}
\textbf{Step 4:} \textit{Second bias correction.} Compute $\widehat{\mathbf{M}}= \widetilde{\mathbf{M}}^{\mathrm{proj}}-\mathcal{B}_2$
where \vspace{-0.7cm}
\begin{align*}
\mathcal{B}_2= \widetilde{\sigma}^2\left[\frac{T}{N} \mathbf{P}_{\widetilde{\bm\beta}}\widehat{\bf\Pi}^{-1} \mathbf{W}_{\bm\beta} ( \widetilde{\mathbf{F}}'\widetilde{\mathbf{F}})^{-1}\widetilde{\mathbf{F}}' + \frac{N}{T} \widetilde{\bm\beta}(\widetilde{\bm\beta}'\widetilde{\bm\beta})^{-1} \mathbf{W}_\mathbf{F}' \widehat{\bf\Psi} \mathbf{P}_{\widetilde \mathbf{F}}\right]
\end{align*} \\
with a variance estimator $\widetilde{\sigma}^2$.\footnote{For the variance estimation, we can use $\widetilde{\mathbf{M}}^{\mathrm{init}}.$ Specifically, when $\mathbf{X}$ is a general regressor matrix, we define $\widetilde{\sigma}^2=(NT)^{-1}\sum_{j=1}^N \sum_{s=1}^T (Y_{js}-X_{js}\widetilde{M}^{\mathrm{init}}_{js})^2$. {When $\mathbf{X}$ is binary, we may use $\widetilde{\sigma}^2=(\sum_{j=1}^N \sum_{s=1}^T X_{js})^{-1} \sum_{j=1}^N \sum_{s=1}^T X_{js}(Y_{js}-\widetilde{M}^{\mathrm{init}}_{js})^2$ instead.}}
\end{algorithmic}
\end{breakablealgorithm}
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 \cite{chen:2019inference} and \cite{xia2021statistical}, who reach the following decomposition after the debias:
\begin{align*}
\widetilde{\mathbf{M}}^{\mathrm{naive}} ={\mathbf{M}^{\star}} +\mathcal{Z}={\mathbf{M}^{\star}} + \underbrace{\widehat{\bf\Pi}^{-1}\mathbf{X} \circ \mathbf{E}}_{\text{dominant noise term}}+ \underbrace{ (\bf{1}_N\bf{1}_T'-\widehat{\bf\Pi}^{-1}\mathbf{X} \circ \mathbf{X}) \circ (\widetilde{\mathbf{M}}^{\mathrm{init}} - {\mathbf{M}^{\star}})}_{\text{higher-order initial estimation error}}.
\end{align*}
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 \textit{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.
\subsection{A new source of bias and Tikhonov-type functions}\label{sec:tikhonov}
After Step 2, we end up with $ \widetilde{\mathbf{M}}^{\mathrm{naive}}$. Write it as
\begin{align*}
\widetilde{\mathbf{M}}^{\mathrm{naive}}= {\mathbf{M}^{\star}} + \mathcal{Z}
\end{align*}
where $\bf\mathcal{Z}$ is the estimation error. Therefore, we have
\begin{align*}
\mathbf{P}_{\widetilde{\bm\beta}} \widetilde{\mathbf{M}}^{\mathrm{naive}} \mathbf{P}_{\widetilde{\mathbf{F}}}-\mathbf{M}^{\star} = (\mathbf{P}_{\widetilde{\bm\beta}} {\mathbf{M}^{\star}} \mathbf{P}_{\widetilde{\mathbf{F}}}-\mathbf{M}^{\star})+\mathbf{P}_{\widetilde{\bm\beta}} \mathcal{Z} \mathbf{P}_{\widetilde{\mathbf{F}}}.
\end{align*}
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
\begin{align*}
\frac{1}{\sqrt{T}} \left\Vert\widetilde{\mathbf{F}}-\mathbf{F}\mathbf{H}'\right\Vert = o_P(1).
\end{align*}
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:
\begin{align*}
\bm{S}_\mathbf{F}\coloneqq \frac{1}{T}\mathbf{H}\mathbf{F}'\mathbf{F}\mathbf{H}', \quad R \times R.
\end{align*}
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 \textit{Tikhonov-regularization function}:
\begin{align*}
f(x) \coloneqq (\mathbf{H}\frac{1}{T}\mathbf{F}' \mathbf{F} \mathbf{H}' +x \mathbf{I}_R )^{-1}.
\end{align*}
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}$ \textit{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 \textit{rescaled Tikhonov-regularization function}:
\begin{align*}
\widetilde{f}(x) \coloneqq \mathbf{H}'(\mathbf{H}\frac{1}{T}\mathbf{F}' \mathbf{F} \mathbf{H}' +x \mathbf{I}_R )^{-1}\mathbf{H}.
\end{align*}
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,
\begin{align*}
\left\|\mathbf{H}'\bm{S}^{-1}_{\widetilde{\mathbf{F}}} \mathbf{H}- \left(\frac{1}{T}\mathbf{F}'\mathbf{F}\right)^{-1}\right\|
\approx \| \widetilde f(x_{NT}) -\widetilde f(0)\|
=O_P( x^2_{NT})
\end{align*}
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:
\begin{align*}
\widehat{\mathbf{M}}:=\widetilde{\mathbf{M}}^{\mathrm{proj}} - \mathcal{B}_2 = \mathbf{P}_{\widetilde{\bm\beta}} {\mathbf{M}^{\star}} \mathbf{P}_{\widetilde{\mathbf{F}}} + \underbrace{\mathbf{P}_{\widetilde{\bm\beta}} \mathcal{Z} \mathbf{P}_{\widetilde{\mathbf{F}}}- \mathcal{B}_2}_{\text{entry-wise negligible}}.
\end{align*}
This is the motivation for introducing $\mathcal{B}_2$ in Step 4.
\section{Asymptotic Results}\label{sec:asympresults}
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}$.
\begin{assumption}[DGP for $\mathbf{E}$]\label{asp:dgpnoise}
\begin{enumerate}
\item[(i)] Conditioning on $({\bm\beta}, \mathbf{F}, \mathbf{W}_{\bm\beta}, \mathbf{W}_\mathbf{F})$, $E_{it}$ is an i.i.d. (across $i$ and $t$) sub-Gaussian random variable with zero mean, a finite variance, $\sigma^2$, and sub-Gaussian norm at most $C \sigma$ for some $C>0.$ Also, $\mathbf{E}$ is independent of $\mathbf{X}$.
\item [(ii)] When $\mathbf{X}$ is a binary matrix taking values in $\{0,1\}$, then condition (i) may be replaced by: There is a noise matrix $\mathbf{E}^{\star}$ that satisfies (i), and $\mathbf{E}$ can be written as $\mathbf{E}=\mathbf{X}\circ \mathbf{E}^{\star}$.
\end{enumerate}
\end{assumption}
We note that Assumption \ref{asp:dgpnoise} (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|$.
\begin{assumption}[DGP for $\mathbf{X}$]\label{asp:dgpX}
\begin{enumerate}
\item[(i)] Conditioning on $({\bm\beta}, \mathbf{F}, \mathbf{W}_{\bm\beta}, \mathbf{W}_\mathbf{F})$, $X_{it}$ is i.i.d. across $t$ for each $i.$ In addition, $X_{it}$ are independent across clusters, and they are allowed to be dependent within clusters. Overall,
$\max_{t \leq T} \max_{j \leq N}\sum_{i=1}^N \left| \mathrm{Cov}(X^2_{it}, X^2_{jt}|{\bm\beta}, \mathbf{F}, \mathbf{W}_{\bm\beta}, \mathbf{W}_\mathbf{F}) \right| \leq C$ for some $C>0.$
\item[(ii)] Let $p_i \coloneqq \mathbb{E}[X_{it}^2|{\bm\beta}, \mathbf{F}, \mathbf{W}_{\bm\beta}, \mathbf{W}_\mathbf{F}]$ and $p_{\min}\coloneqq \min_{i \leq N}p_i$. We assume that $p_{\min}$ is bounded away from zero.
\item[(iii)] $\max_{i\leq N}\max_{t \leq T} X^2_{it} < C$ for some $C>0$ almost surely.
\end{enumerate}
\end{assumption}
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.
\begin{assumption}[Structure of $\mathbf{M}^{\star}$]\label{asp:incoherence}
\begin{enumerate}
\item[(i)] Low-rank factor structure: We assume $\mathbf{M}^{\star}={\bm\beta} \mathbf{F}'$ with an $N \times r$ loading matrix $\bm\beta$ and a $T \times r$ factor matrix $\mathbf{F}$. We assume that $r$ is bounded.
\item[(ii)] Incoherence: The matrix $\mathbf{M}^{\star}$ is incoherent in that:
\begin{align*}
\left\Vert\bm\beta\right\Vert_{2,\infty} \lesssim \frac{\sigma_{\max}({\bm\beta})}{\sqrt{N}} \quad \text{and} \quad \left\Vert\mathbf{F}\right\Vert_{2,\infty} \lesssim \frac{\sigma_{\max}({\mathbf{F}})}{\sqrt{T}}.
\end{align*}
\item[(iii)] $T/N \rightarrow C$ for some $C \in (0,\infty)$.
\end{enumerate}
\end{assumption}
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."
\begin{assumption}[Diversified weights and weak factors]\label{asp:diversifiedweightsandweakfactor}
We define $\mathbf{H}_{\bm\beta}\coloneqq N^{-1} \mathbf{W}_{\bm\beta}'{\bm\beta}$ and $\mathbf{H}_\mathbf{F}\coloneqq T^{-1}\mathbf{W}_\mathbf{F}' \mathbf{F}$.
\begin{enumerate}
\item[(i)] $\mathbf{W}_{\bm\beta}$ is an $N \times R$ matrix and $\mathbf{W}_\mathbf{F}$ is a $T \times R$ matrix, where $R$ is bounded and $R \geq r$. $\mathbf{W}_{\bm\beta}$ and $\mathbf{W}_\mathbf{F}$ are independent of $\mathbf{E}.$\footnote{When Assumption \ref{asp:dgpnoise} (ii) holds with $\mathbf{E}=\mathbf{X} \circ \mathbf{E}^{\star}$, we assume that $\mathbf{W}_{\bm\beta}$ and $\mathbf{W}_\mathbf{F}$ are independent of $\mathbf{E}^{\star}.$} Also, almost surely, $\max\{\left\Vert\mathbf{W}_{\bm\beta}\right\Vert_{2, \infty}, \left\Vert\mathbf{W}_\mathbf{F}\right\Vert_{2, \infty}\}<C$, $\min\{\sigma_{R}(N^{-1} \mathbf{W}_{\bm\beta}'\mathbf{W}_{\bm\beta}), \sigma_{R}(T^{-1} \mathbf{W}_\mathbf{F}'\mathbf{W}_\mathbf{F}) \}>c$ for some $c, C>0$, and, the ranks of $\mathbf{H}_{\bm\beta}$ and $\mathbf{H}_{\mathbf{F}}$ are $r$.
\item[(ii)] $\log N \ll \sigma_{\min}({\bm\beta}) \asymp \sigma_{\max}({\bm\beta}) \lesssim \sqrt{N}$ and $\sigma_{\min}(\mathbf{H}_{\bm\beta}) \asymp \sigma_{\max}(\mathbf{H}_{\bm\beta}) \asymp \sigma_{\max}({\bm\beta})/ \sqrt{N} $.
\item[(iv)]$\sigma_{\min}(\mathbf{H}_\mathbf{F}) \asymp \sigma_{\max}(\mathbf{H}_\mathbf{F}) \asymp C$ for some constant $C>0$.
\end{enumerate}
\end{assumption}
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:
\begin{align*}
{\bm\beta} = a_N {\bm\beta}_0 \quad \text{and} \quad \sigma_{\min}\left(\frac{1}{N}\mathbf{W}_{\bm\beta}' {\bm\beta}_0 \right) \asymp \sigma_{\max}\left(\frac{1}{N}\mathbf{W}_{\bm\beta}' {\bm\beta}_0 \right)\asymp C
\end{align*}
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{asp:diversifiedweightsandweakfactor} (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 \cite{onatski2012asymptotics} in the pure factor model. In contrast, the usual requirement for inference in the PCA-based approaches \citep[see][]{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.
\begin{assumption}[Block shape]\label{asp:blocksize}
We assume i) $\max\{|\mathcal{I}|^2 \vartheta^{3} \log^{3}N,|\mathcal{I}| \vartheta^{6} \log^{6}N \} =o(N)$, ii) $\max\{|\mathcal{T}|^2 \vartheta^{3} \log^{3}T,|\mathcal{T}| \vartheta^{6} \log^{6}T \} = o(N)$, and iii) $\sqrt{\max\{|\mathcal{I}|, |\mathcal{T}|\}}\log N \ll \sigma_{\min} ({\bm\beta})$.
\end{assumption}
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{sec:initial}.
Define the covariance estimator:
\begin{align*}
\widehat{\mathcal{V}}_{\mathcal{G}_{\mathrm{bl}}}
\coloneqq \frac{\widetilde{\sigma}^2}{|\mathcal{T}|^2 N^2}\sum_{t \in \mathcal{T}} \sum_{j=1}^N (\widehat{\mathbf{M}}_{\mathcal{I}, \cdot} \widetilde{\mathbf{F}}(\widetilde{\mathbf{F}}'\widetilde{\mathbf{F}})^{-1} \mathbf{W}_{{\bm\beta},j } \widehat{p}_{j}^{-1}X_{jt})^2 + \frac{\widetilde{\sigma}^2}{|\mathcal{I}|^2 T^2}\sum_{i \in \mathcal{I}} \sum_{s=1}^T (\widehat{\mathbf{M}}_{\cdot, \mathcal{T}}'\widetilde{\bm\beta}(\widetilde{\bm\beta}'\widetilde{\bm\beta})^{-1}\mathbf{W}_{\mathbf{F},s } \widehat{p}_i^{-1} X_{is} )^2
\end{align*}
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.}
\begin{theorem}[Feasible CLT for block average]\label{thm:feasibleclt}
Suppose $R \geq r$ and Assumption \ref{asp:dgpnoise}-\ref{asp:blocksize} hold. In addition, the initial estimator $\widetilde{\mathbf{M}}^{\mathrm{init}}$ is as constructed in Section \ref{sec:initial}. Assume that $\| |\mathcal{I}|^{-1}\sum_{i \in \mathcal{I}} {\bm\beta}_{i }\|> C \sigma_{\min}({\bm\beta})/\sqrt{N}$ and $\||\mathcal{T}|^{-1} \sum_{t \in \mathcal{T}} \mathbf{F}_{t }\|>C$ for some $C>0.$ Then, we have
\vspace{-0.7cm}
\begin{align*}
\widehat{\mathcal{V}}_{\mathcal{G}_{\mathrm{bl}}}^{-\frac{1}{2}} \frac{1}{|\mathcal{G}_{\mathrm{bl}}|}\sum_{(i,t) \in \mathcal{G}_{\mathrm{bl}}}(\widehat{M}_{it}-M^{\star}_{it}) \overset{D}\longrightarrow \mathcal{N} (0,1).
\end{align*}
\end{theorem}
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 \cite{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.
\section{Construction of the Initial Estimator}\label{sec:initial}
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.
\noindent\textbf{Full sample.} We define
\begin{align}
\widetilde{\mathbf{M}}^{\mathrm{full}} \coloneqq \operatorname*{arg\,min}_{\mathbf{M} \in \mathbb{R}^{N \times T}} \frac{1}{2} \sum_{j=1}^N \sum_{s=1}^T \widehat{p}_j^{-1}(Y_{js}-X_{js}M_{js})^2 + \lambda \left\Vert\mathbf{M}\right\Vert_* \label{eq:full}
\end{align}
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 \cite{ma2019missing} and \cite{choi2023inference}.
\noindent\textbf{Restricted 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{asp:blocksize}.\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 \textit{outside of the group}: $i\notin \mathcal I$ and $t\notin\mathcal T$,
\begin{align*}
\widetilde{\mathbf{M}}^{\mathrm{rest}}:= \operatorname*{arg\,min}_{\mathbf{M} } \frac{1}{2} \sum_{i\notin\mathcal{I}} \sum_{t\notin\mathcal{T}} \widehat{p}_i^{-1}(Y_{it}-X_{it}M_{it})^2 + \lambda^{\mathrm{rest}} \left\Vert\mathbf{M}\right\Vert_*.
\end{align*}
\noindent\textbf{Merging two estimates.} Then, we define the initial estimator $\widetilde{\mathbf{M}}^{\mathrm{init}}$ as
\begin{align*}
\widetilde{M}^{\mathrm{init}}_{it}=
\begin{cases}
\widetilde{M}^{\mathrm{rest}}_{it} & \text{if $i\notin \mathcal{I}$ and $t\notin\mathcal{T}$,} \cr
\widetilde{M}^{\mathrm{full}}_{it} & \text{otherwise.}
\end{cases}
\end{align*}
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.
\section{Statistical Applications}\label{sec:application}
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.
\subsection{Heterogeneous Treatment Effects Estimation}\label{sec:treatmenteffect}
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., \cite{rubin:1974, imbens:2015}, we assume that, for each $(i,t),$ there exist two \textit{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:
\begin{align*}
Z^{(\iota)}_{it}= M^{(\iota)}_{it}+ E^{\star}_{it} = {\bm\beta}_{i}' \mathbf{F}_{t}^{(\iota)} + E^{\star}_{it},
\end{align*}
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:
\begin{align}
\mathbf{Y}^{(\iota)}= \mathbf{X}^{(\iota)} \circ \mathbf{Z}^{ (\iota)} = \mathbf{X}^{(\iota)} \circ \mathbf{M}^{(\iota)} +\mathbf{X}^{(\iota)} \circ \mathbf{E}^{\star} = \mathbf{X}^{(\iota)} \circ ({\bm\beta} \mathbf{F}^{(\iota) \prime}) + \underbrace{\mathbf{X}^{(\iota)} \circ \mathbf{E}^{\star}}_{\coloneqq \mathbf{E}^{(\iota)}}, \label{eq:treatment}
\end{align}
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
\begin{align*}
\frac{1}{|\mathcal{G}|} \sum_{(i,t) \in \mathcal{G}} \Gamma_{it} = \frac{1}{|\mathcal{G}|} \sum_{(i,t) \in \mathcal{G}} (M^{(1)}_{it}-M^{(0)}_{it}).
\end{align*}
Then, a natural choice for the estimator for the group treatment effect would be
\begin{align*}
\frac{1}{|\mathcal{G}|} \sum_{(i,t) \in \mathcal{G}} \widehat{\Gamma}_{it} = \frac{1}{|\mathcal{G}|} \sum_{(i,t) \in \mathcal{G}} (\widehat{M}_{it}^{(1)}-\widehat{M}_{it}^{(0)}).
\end{align*}
Assuming that the assumptions in Section \ref{sec:asympresults} 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{sec:asympresults}, 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,$
\begin{align*}
\widehat{\mathcal{V}}_{\mathcal{G}_{\mathrm{bl}}}^{(\iota)}
&\coloneqq \frac{\widetilde{\sigma}^2}{|\mathcal{T}|^2 N^2}\sum_{t \in \mathcal{T}} \sum_{j=1}^N (\widehat{\mathbf{M}}^{(\iota)}_{\mathcal{I}, \cdot} \widetilde{\mathbf{F}}^{(\iota)}(\widetilde{\mathbf{F}}^{(\iota) \prime}\widetilde{\mathbf{F}}^{(\iota)})^{-1} \mathbf{W}_{{\bm\beta},j } (\widehat{p}^{(\iota)}_{j})^{-1}X^{(\iota)}_{jt})^2 \\
&\quad + \frac{\widetilde{\sigma}^2}{|\mathcal{I}|^2 T^2}\sum_{i \in \mathcal{I}} \sum_{s=1}^T (\widehat{\mathbf{M}}^{(\iota) \prime}_{\cdot, \mathcal{T}}\widetilde{\bm\beta}^{(\iota) }(\widetilde{\bm\beta}^{(\iota) \prime}\widetilde{\bm\beta}^{(\iota)})^{-1}\mathbf{W}_{\mathbf{F},s }^{(\iota)} (\widehat{p}^{(\iota)}_i)^{-1} X^{(\iota)}_{is} )^2
\end{align*}
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}$.
\begin{theorem}[Feasible CLT for heterogeneous treatment effect]\label{thm:feasibleclt-treat}
Suppose $R \geq r$ and the assumption in Theorem \ref{thm:feasibleclt} hold for both $(0)$ and $(1).$ Also, the initial estimators $\widetilde{\mathbf{M}}^{\mathrm{init}, (\iota)}$, $\iota=0,1,$ are as constructed in Section \ref{sec:initial}. Then,
\begin{align*}
(\widehat{\mathcal{V}}_{\mathcal{G}_{\mathrm{bl}}}^{(0)}+\widehat{\mathcal{V}}^{(1)}_{\mathcal{G}_{\mathrm{bl}}})^{-\frac{1}{2}} \frac{1}{|\mathcal{G}_{\mathrm{bl}}|}\sum_{(i,t) \in \mathcal{G}_{\mathrm{bl}}}(\widehat{\Gamma}_{it}-\Gamma_{it}) \overset{D}\longrightarrow \mathcal{N} (0,1).
\end{align*}
\end{theorem}
\subsection{ Multiple testing with incomplete data}
In large-scale multiple testing, it is crucial to address the well-known confounding factors and the resulting strong correlations \citep[e.g.][]{leek2008general, friguet2009factor, wang2017confounder,fan2019farmtest}. We consider the model, for $t=1,...,T$:
\begin{align*}
\mathbf{Z}_t= {\bm\mu} + \mathbf{U}_t, \quad \text{where} \quad \mathbf{U}_t={\bm\beta \mathbf{F}_t} + \mathbf{E}^{\star}_t.
\end{align*}
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:
\begin{align*}
H_0^i: \mu_i = 0, \quad i=1, \ldots, N.
\end{align*}
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 \citep[e.g.,][]{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{alg:estimation} on $\mathbf{Y}^{\mathrm{d}}$ to obtain $\widehat{\mathbf{M}}.$\footnote{The initialization in Section \ref{sec:initial}, though necessary for the CLT for $\mathbf{M}^{\star}$, is not required for testing $\bm\mu$. Instead, we simply implement the full-sample estimation \eqref{eq:full} on $\mathbf{Y}^{\mathrm{d}}$ just once and use it as $\widetilde{\mathbf{M}}^{\mathrm{init}}$.} Our proposed estimator for ${\bm\mu}$ is
\begin{align*}
\widehat{\mu}_i&= \bar{Y}_i-\frac{1}{|\mathcal{T}_i|} \sum_{t \in \mathcal{T}_i} \widehat{M}_{it} \quad \text{for $i =1, \ldots, N$}.
\end{align*}
\begin{theorem}\label{thm:multipletesting}
Suppose Assumption \ref{asp:dgpnoise} (ii), \ref{asp:dgpX} hold. In addition, suppose that $\min_i |\mathcal{T}_i| > \epsilon T$ for some $\epsilon>0$ almost surely. Assume the following:
\begin{enumerate}
\item[(i)] When $r>0$, in addition to Assumption \ref{asp:incoherence}, \ref{asp:diversifiedweightsandweakfactor}, we assume $\vartheta^6 \log^7 N \ll N$, and $\log^{\frac{3}{2}}N \ll \sigma_{\min}({\bm\beta})$ hold. Also, $\mathbb E\mathbf{F}_t=\bf{0} $ and $\{\mathbf{F}_t\}_{t \leq T}$ is i.i.d.
\item[(ii)] When $r=0$, $\mathbf{W}_{\bm\beta}$ and $\mathbf{W}_\mathbf{F}$ are independent of $\mathbf{E}^{\star},$ and $\max\{\left\Vert\mathbf{W}_{\bm\beta}\right\Vert_{2, \infty}, \left\Vert\mathbf{W}_\mathbf{F}\right\Vert_{2, \infty}\}<C$, $\min\{\sigma_{R}(N^{-1} \mathbf{W}_{\bm\beta}'\mathbf{W}_{\bm\beta}), \sigma_{R}(T^{-1} \mathbf{W}_\mathbf{F}'\mathbf{W}_\mathbf{F}) \}>c$ for some $c, C>0$ almost surely. Also, $T/N \rightarrow C$ for some $C \in (0, \infty)$.
\end{enumerate}
Then, uniformly for $i=1, \ldots, N$ and all bounded $R\geq r \geq 0$, we have
\begin{align}
\widehat{\mu}_i-\mu_i=\frac{1}{|\mathcal{T}_i|} \sum_{t \in \mathcal{T}_i} E_{it}+o_P\left(\sqrt{\frac{1}{|\mathcal{T}_i|\log N}}\right). \label{eq:FDRresult}
\end{align}
\end{theorem}
The above expansion \eqref{eq:FDRresult} 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 \cite{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{eq:FDRresult}), \cite{liu2014phase} showed that the false discovery rate of the B-H procedure can be controlled below $\tau$ asymptotically.
\section{Choices of Diversified Weights}\label{sec:choiceofW}
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 \textit{informative} with respect to the underlying matrix parameter.
\subsection{Observed characteristics} \label{subsec:obschr}
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 \citep{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, \cite{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 \citep[e.g.,][]{connor2012efficient, fan2016projected}. Also, observed macroeconomic factors or Fama-French factors \citep{fama1993common} can be used as $\mathbf{f}_t$.
\subsection{Exploiting extra samples}\label{subsec:extramsample}
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$.
\subsection{Initial transformation}\label{sec:initialtransformation}
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:
\begin{enumerate}
\item[(i)] For $t_0=1$, $\{E_{i,t_0}: \forall i \}$ is independent of $\{E_{i,t}: \forall i\}$ for all $t\geq 2$; and
\item[(ii)] There is a known individual $i$, say $i=1$, so that the noise $\{E_{1,t}: \forall t\}$ is independent of $\{E_{i,t}: \forall t\}$ for all other $i\in \{1, \ldots, N\}/\{1\}$.
\end{enumerate}
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.
\section{Simulation results}\label{sec:simulation}
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{sec:choiceofW}, 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\}$.
\begin{table}[h]
\begin{center}
\begin{tabular}{c|c|l|lll}
\toprule
\multicolumn{1}{l|}{Sample size} & \multicolumn{1}{l|}{Factor strength} & Diversified weights & Block & CS & Serial \\ \hline \hline
\multirow{6}{*}{$N=T=200$} & \multirow{3}{*}{$\alpha=1$} & Observed & 0.932 & 0.940 & 0.950 \\
& & Extra & 0.913 & 0.940 & 0.946 \\
& & Initial & 0.874 & 0.945 & 0.951 \\ \cline{2-6}
& \multirow{3}{*}{$\alpha=0.5$} & Observed & 0.906 & 0.940 & 0.947 \\
& & Extra & 0.904 & 0.945 & 0.947 \\
& & Initial & 0.915 & 0.949 & 0.950 \\ \hline
\multirow{6}{*}{$N=T=400$} & \multirow{3}{*}{$\alpha=1$} & Observed & 0.942 & 0.954 & 0.952 \\
& & Extra & 0.921 & 0.943 & 0.948 \\
& & Initial & 0.884 & 0.945 & 0.944 \\ \cline{2-6}
& \multirow{3}{*}{$\alpha=0.5$} & Observed & 0.894 & 0.937 & 0.935 \\
& & Extra & 0.891 & 0.949 & 0.950 \\
& & Initial & 0.901 & 0.945 & 0.945 \\ \bottomrule
\end{tabular}
\end{center}
{\raggedright\textbf{\fname@algorithm~\thealgorithm} #\small Coverage probabilities of our debiased estimators in the matrix completion model. The target probability is 0.95. ``Observed,'' ``Extra,'' and ``Initial'' refer to our debiased estimators that use diversified weights of observed characteristics, the extra sample averages, and the initial transformation, respectively. Three averages are estimated: (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\}$. In all specifications, we set $R=4$ while the true rank is $r=2$.\par}
\ifx\relax#\relax\relax
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\small Coverage probabilities of our debiased estimators in the matrix completion model. The target probability is 0.95. ``Observed,'' ``Extra,'' and ``Initial'' refer to our debiased estimators that use diversified weights of observed characteristics, the extra sample averages, and the initial transformation, respectively. Three averages are estimated: (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\}$. In all specifications, we set $R=4$ while the true rank is $r=2$.}
\else
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\relax}
\fi
\kern2pt\hrule\kern2pt
\label{table:MCtable}
\end{table}
Table \ref{table:MCtable} 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 \cite{chernozhukov2018inference,chernozhukov2023inference, choi2023inference}. The choice of $\lambda$ in these studies involves a choice of a small constant $c>0.$ \cite{chernozhukov2018inference,chernozhukov2023inference} set $c=1/10$ while \cite{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.
\begin{figure}[h]
\begin{center}
\includegraphics[width=0.9\textwidth]{MC200200alpha05_revision_third.jpg} \vspace{-0.5cm}
\end{center}
{\raggedright\textbf{\fname@algorithm~\thealgorithm} #\small Histograms of the standardized estimates in the noisy matrix completion model when $N=T=200$ and $\alpha=0.5$. ``Nucl'' in the first row presents the nuclear norm penalized estimators. The solid line represents the standard normal distribution. For our debiased estimators, we set $R=4$ while the true rank is $r=2$.\par}
\ifx\relax#\relax\relax
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\small Histograms of the standardized estimates in the noisy matrix completion model when $N=T=200$ and $\alpha=0.5$. ``Nucl'' in the first row presents the nuclear norm penalized estimators. The solid line represents the standard normal distribution. For our debiased estimators, we set $R=4$ while the true rank is $r=2$.}
\else
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\relax}
\fi
\kern2pt\hrule\kern2pt
\label{fig:MC200200alpha05}
\centering
\end{figure}
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{fig:MC200200alpha05}. 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.
\section{Empirical Study: Impact of the U.S. Presidential Election on the Federal Grants Allocation}\label{sec:empirical}
We apply the inferential theory presented in Section \ref{sec:treatmenteffect} 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 \citep[e.g.,][]{larcinese2006allocating,berry2010president}.
By employing the treatment effect estimation application presented in Section \ref{sec:treatmenteffect}, 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{sec:treatmenteffect}, 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 \eqref{eq:treatment}.
Before proceeding to tests, we present the singular values of the nuclear norm penalized estimators for the treated and control samples (Figure \ref{fig:singularvalues}). 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.
\begin{figure}[H]
\begin{center}
\includegraphics[width= 0.8\textwidth, height=5cm]{singularvalues.jpg} \vspace{-0.5cm}
\end{center}
{\raggedright\textbf{\fname@algorithm~\thealgorithm} #\small The singular values of the (full-sample) nuclear norm penalized estimators for the treated and control sample, i.e., $\widetilde{\mathbf{M}}^{\mathrm{full}, (1)}$ and $\widetilde{\mathbf{M}}^{\mathrm{full}, (0)}$, respectively, in descending order. The left corresponds to the treated sample, while the right corresponds to the control sample.\par}
\ifx\relax#\relax\relax
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\small The singular values of the (full-sample) nuclear norm penalized estimators for the treated and control sample, i.e., $\widetilde{\mathbf{M}}^{\mathrm{full}, (1)}$ and $\widetilde{\mathbf{M}}^{\mathrm{full}, (0)}$, respectively, in descending order. The left corresponds to the treated sample, while the right corresponds to the control sample.}
\else
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\relax}
\fi
\kern2pt\hrule\kern2pt
\label{fig:singularvalues}
\end{figure}
We construct the diversified weighting matrices following the observed characteristics approach explained in Section \ref{subsec:obschr}. 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 \textit{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{fig:individualstateeffects}. 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{fig:map} 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?
\begin{figure}[H]
\begin{center}
\includegraphics[width= 1.0\textwidth, height=8.0cm]{individualstateeffects_revision.jpg} \vspace{-0.5cm}
\end{center}
{\raggedright\textbf{\fname@algorithm~\thealgorithm} #Individual state effects and corresponding t-statistics\par}
\ifx\relax#\relax\relax
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#Individual state effects and corresponding t-statistics}
\else
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\relax}
\fi
\kern2pt\hrule\kern2pt
\label{fig:individualstateeffects}
\end{figure}
We observe that many of the states colored in dark green in Figure \ref{fig:map} 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{tab:numberofswing}.
Based on Table \ref{tab:numberofswing}, we estimate the overall time average treatment effects for the state groups. The first plot in Figure \ref{fig:loyaleffects} 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.}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.8\textwidth, height=7.7cm]{MapChart_Map_revision.png} \vspace{-0.5cm}
\end{center}
{\raggedright\textbf{\fname@algorithm~\thealgorithm} #\small The states in dark green indicate the rejection of the null hypothesis at a significance level of 1\%, whereas the states in light green correspond to a significance level of 5\%. In the gray-colored states, the null hypotheses are not rejected. This figure is created with \textsc{MapChart}.\par}
\ifx\relax#\relax\relax
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\small The states in dark green indicate the rejection of the null hypothesis at a significance level of 1\%, whereas the states in light green correspond to a significance level of 5\%. In the gray-colored states, the null hypotheses are not rejected. This figure is created with \textsc{MapChart}.}
\else
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\relax}
\fi
\kern2pt\hrule\kern2pt
\label{fig:map}
\end{figure}
\begin{table}[H]
\begin{center}
\begin{tabular}{lll}
\toprule
Group & \# of swings & States \\ \hline \hline
Swing states & 8$\sim$ & FL, GA, LA, OH \\ \hline
Weak swing states & 6$\sim$7 & AR, IA, KY, MS, MO, PA, TN, WV, WI \\ \hline
Neutral states & 5 & AL, CO, DE, HI, MD, MI, NV, NH, NM, NY, NC, RI \\ \hline
Weak loyal states & 3$\sim$4 & \begin{tabular}[c]{@{}l@{}}AZ, CA, CT, IL, IN, ME, MA, MN, MT, NJ, OR, \\ SC, TX, VT, VA, WA\end{tabular} \\ \hline
Loyal states & 0$\sim$2 & AK, DC, ID, KS, NE, ND, OK, SD, UT, WY \\ \bottomrule
\end{tabular}
\vspace{-0.5cm}
\end{center}
{\raggedright\textbf{\fname@algorithm~\thealgorithm} #\small The counts of swings\par}
\ifx\relax#\relax\relax
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\small The counts of swings}
\else
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\relax}
\fi
\kern2pt\hrule\kern2pt
\label{tab:numberofswing}
\end{table}
\begin{figure}[H]
\begin{center}
\includegraphics[width=1.0\textwidth, height=4.0cm]{loyaleffects_revision.jpg} \vspace{-0.5cm}
\end{center}
{\raggedright\textbf{\fname@algorithm~\thealgorithm} #\small Loyalty effects and corresponding t-statistics\par}
\ifx\relax#\relax\relax
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\small Loyalty effects and corresponding t-statistics}
\else
\addcontentsline{loa}{algorithm}{\protect\numberline{\thealgorithm}#\relax}
\fi
\kern2pt\hrule\kern2pt
\label{fig:loyaleffects}
\end{figure}
\onehalfspacing
\small
\bibliographystyle{apalike}
\bibliography{reference}