text:
\begin{document}
\begin{frontmatter}

\title{Regression adjustment in completely randomized experiments with a diverging number of covariates}
\runtitle{ANCOVA with a diverging number of covariates}

\begin{aug}
\author{\fnms{Lihua} \snm{Lei}\ead[label=e1]{lihua.lei@berkeley.edu}}
\and
\author{\fnms{Peng} \snm{Ding}\ead[label=e1]{pengdingpku@berkeley.edu}}
\affiliation{University of California, Berkeley}

\runauthor{L. Lei and P. Ding}
\end{aug}

\begin{abstract}
Extending R. A. Fisher and D. A. Freedman's results on the analysis of covariance, \cite{lin13} proposed an ordinary least squares adjusted estimator of the average treatment effect in completely randomized experiments.
We further study its statistical properties under the potential outcomes model in the asymptotic regimes allowing for a diverging number of covariates. We show that \cite{lin13}'s estimator is consistent when $\kappa \log p \rightarrow 0$ and asymptotically normal when $\kappa p \rightarrow 0$ under mild moment conditions, where $\kappa$ is the maximum leverage score of the covariate matrix. In the favorable case where leverage scores are all close together, his estimator is consistent when $p = o(n / \log n)$ and is asymptotically normal when $p = o(n^{1/2})$. In addition, we propose a bias-corrected estimator that is consistent when $\kappa \log p\rightarrow 0$ and is asymptotically normal, with the same variance in the fixed-$p$ regime, when $\kappa^{2} p \log p \rightarrow 0$. In the favorable case, the latter condition reduces to $p = o(n^{2/3} / (\log n)^{1/3})$.
Similar to \cite{lin13}, our results hold for non-random potential outcomes and covariates without any model specification. Our analysis requires novel analytic tools for sampling without replacement, which complement and potentially enrich the theory in other areas such as survey sampling, matrix sketching, and transductive learning.
\end{abstract}

\begin{keyword}[class=MSC]
\kwd[Primary ]{62J05}
\kwd{}
\kwd[secondary ]{62E20}
\end{keyword}

\begin{keyword}
\kwd{causal inference}
\kwd{average treatment effect}
\kwd{regression adjusted estimator}
\kwd{high-dimensional covariates}
\kwd{model misspecification}
\kwd{potential outcomes}
\end{keyword}

\end{frontmatter}

\section{Introduction}\label{sec:intro}

\subsection{Potential outcomes and Neyman's randomization model}
We use potential outcomes to define causal effects \citep{neyman23}. Let $Y_{i}(1)$ and $ Y_{i}(0)$ be the potential outcomes if unit $i\in\{1,\ldots, n \}$ receives the treatment and control, respectively.
\citet{neyman23} treated all the potential outcomes as fixed quantities, and defined the average treatment effect (ATE) as
$
\tau \equiv  n^{-1}\sum_{i=1}^{n}\tau_{i},
$
where $ \tau_{i} = Y_{i}(1) - Y_{i}(0)$ is the individual treatment effect for unit $i.$ In a completely randomized experiment, the experimenter randomly assigns $n_1$ units to the treatment group and $n_0$ units to the control group, with $n=n_1+n_0$. Let $T_{i}$ denote the assignment of the $i$-th unit where $T_{i} = 1$ corresponds to the treatment and $T_i=0$ corresponds to the control. For unit $i$, only $Y_{i}^{\mathrm{obs}}  = Y_{i}(T_{i})$ is observed while the other potential outcome $Y_{i}(1 - T_{i})$ is missing.
Although $(Y_{i}(1), Y_{i}(0))_{i=1}^n$ are fixed, the $Y_{i}^{\mathrm{obs}}$'s are random due to the randomization of the $T_{i}$'s.

\citet[][Chapter 9]{scheffe1960analysis} called the above formulation the \emph{randomization model}, under which all potential outcomes are fixed and the randomness comes solely from the treatment indicators. This finite-population perspective has a long history for analyzing randomized experiments \citep[e.g.][]{neyman23, neyman35, kempthorne52, imbens15, mukerjee2018using, fogarty2018regression, middleton2018unified}.  In contrast, the super-population perspective \citep[e.g.][]{tsiatis2008covariate, berk2013covariance, pitkin2017asymptotically} assumes that the potential outcomes and other individual characteristics are independent and identically distributed (i.i.d.) draws from some distribution. Two perspectives are both popular in the literature, but they are different in the source of randomness: the finite-population perspective quantifies the uncertainty of the sampling procedure in a single ``universe'' of units; in contrast, the super-population perspective also considers the uncertainty across multiple, possibly infinite, ``universes'' of units.

We use the conventional notation $O(\cdot),o(\cdot)$, $O_\P(\cdot)$ and $o_\P(\cdot)$. Let $\one$ denote the vector with all entries $1$, $\id$ denote an identity matrix, and $\projone= \id  -  (\one\tran \one )^{-1} \one \one\tran$ denote the projection matrix orthogonal to $\one$, with appropriate dimensions depending on the context.   Let $\|\cdot\|_{q}$ be the vector $q$-norm, i.e. $\|\alpha\|_{q} = (\sum_{i=1}^{n}|\alpha_{i}|^{q})^{ 1 / q }$ and $\|\alpha\|_{\infty} = \max_{1\leq i \leq n}|\alpha_{i}|$. Let $\|\cdot\|_{\op}$ denote operator norm and $\| \cdot \|_F$ denote the Frobenius norm of matrices. Let $N(0, 1)$ denote the standard normal distribution, and $t(\nu)$ denote standard $t$ distribution with degrees of freedom $\nu$ with $t(1)$ being the standard Cauchy distribution. Let $\stackrel{\text{d}}{\rightarrow}$ and $\stackrel{\P}{\rightarrow}$ denote convergences in distribution and in probability.

\subsection{Regression-adjusted average treatment effect estimates}

Let $\T_{t} = \{i: T_{i} = t\}$ be the indices and $n_{t} = |\T_{t}|$ be the fixed sample size for treatment arm $t\in \{0, 1\}$. We consider a completely randomized experiment in which $\mathcal{T}_1$ is a random size-$n_1$ subset of $\{1,\ldots, n\}$ uniformly over all $\binom{n}{n_1}$ subsets.  The simple difference-in-means estimator
\[
\tauunadj = \frac{1}{n_{1}}\sum_{i\in \T_{1}}Y_{i}^{\mathrm{obs}} - \frac{1}{n_{0}}\sum_{i\in \T_{0}}Y_{i}^{\mathrm{obs}}
= \frac{1}{n_{1}}\sum_{i\in \T_{1}}Y_{i}(1) - \frac{1}{n_{0}}\sum_{i\in \T_{0}}Y_{i}(0)
\]
is unbiased with variance $S_{1}^{2}/n_{1} + S_{0}^{2} / n_{0} - S_{\tau}^{2}/n$ \citep{neyman23}, where $S_{1}^{2}, S_{0}^{2}$ and $ S_{\tau}^{2}$ are the finite-population variances of the $Y_{i}(1)$'s, $Y_{i}(0)$'s and $\tau_{i}$'s, respectively.

The experimenter usually collects pre-treatment covariates. If the covariates are predictive of the potential outcomes, it is intuitive to incorporate them in the analysis to improve the estimation efficiency. Suppose unit $i$ has a $p$-dimensional vector of pre-treatment covariates $x_{i}\in \R^{p}$. Early works on the analysis of covariance assumed constant treatment effects \citep{fisher37, kempthorne52, hinkelmann07}, under which a commonly-used treatment effect estimate is the coefficient of the treatment indicator of the ordinary least squares (OLS) fit of the $Y_{i}^{\mathrm{obs}}$'s on $T_i$'s and $x_i$'s. \citet{freedman08a} criticized this standard approach, showing that (a) it can be even less efficient than $\tauunadj $ in the presence of treatment effect heterogeneity, and (b) the estimated standard error based on the OLS can be inconsistent for the true standard error under the randomization model.

\citet{lin13} proposes a simple solution. Without loss of generality, we center the covariates at $n^{-1}\sum_{i=1}^{n}x_{i} = 0$ because otherwise we can replace $x_{i}$ by $x_{i} -  n^{-1}\sum_{i=1}^{n}x_{i} $. His estimator for the ATE is the coefficient of the treatment indicator in the OLS fit of the $Y_{i}^{\mathrm{obs}}$'s on $T_i$'s, $x_i$'s and the {\it interaction terms} $T_ix_i$'s. He further shows that the Eicker--Huber--White standard error \citep[e.g.][]{mackinnon13} is consistent for the true standard error. \citet{lin13}'s results hold under the finite-population randomization model, without assuming that the linear model is correct.

We use an alternative formulation of regression adjustment and consider the following family of covariate-adjusted ATE estimator:
\begin{equation}\label{eq:family_regadj}
\hat{\tau}(\beta_{1}, \beta_{0}) = \frac{1}{n_{1}}\sum_{i\in \T_{1}}(Y_{i}^{\mathrm{obs}} - x_{i}\tran \beta_{1}) - \frac{1}{n_{0}}\sum_{i\in \T_{0}}(Y_{i}^{\mathrm{obs}} - x_{i}\tran \beta_{0}).
\end{equation}
Because
$  \E \lb   n_{t}^{-1}\sum_{i\in \T_{t}} x_{i}\tran \beta_t \rb  = 0$,
the estimator in \eqref{eq:family_regadj} is unbiased for any fixed coefficient vectors $\beta_t \in \R^{p}$ $(t=0,1)$. It is the difference-in-means estimator with potential outcomes replaced by $(Y_{i}(1) - x_{i}\tran \beta_{1}, Y_{i}(0) - x_{i}\tran \beta_{0})_{i=1}^n $.

Let $Y(t) = (Y_{1}(t), \ldots, Y_{n}(t))\tran \in \R^{n}$ denote the vector of potential outcomes under treatment $t$ $(t=0,1)$, $X = (x_{1} , \ldots, x_{n} )\tran $ denote the matrix of covariates. Without loss of generality, we assume
\begin{equation}\label{eq:centered}
\one\tran X = 0 \text{ and } \rank(X) = p,
\end{equation}
i.e., the covariate matrix has centered columns and full column rank. Otherwise, we transform $X$ to $\projone X$ and remove the redundant columns to ensure the full column rank condition.
Let $\beta_{t}$ be the population OLS coefficient of regressing $Y(t)$ on $X$ with an intercept:
\begin{eqnarray}
(\mu_{t}, \beta_{t}) &=& \argmin_{\mu\in \R, \beta\in \R^{p}}\|Y(t) - \mu\one - X\beta\|_{2}^{2}  \label{eq:betat} \\
&=& \lb \frac{1}{n}\sum_{i=1}^{n}Y_{i}(t) ,     (X\tran X)^{-1}X\tran Y(t)  \rb  , \label{eq:mut-ols}
\end{eqnarray}
where \eqref{eq:mut-ols} holds because $X$ is orthogonal to $\one$. \citet[][Example 9]{li2017general} show that the OLS coefficients $(\beta_{1}, \beta_{0})$ in \eqref{eq:betat} minimize the variance of the estimator defined in \eqref{eq:family_regadj}.

The classical analysis of covariance chooses $\beta_{1} = \beta_{0} = \hat{\beta}$, the coefficient of the covariates in the OLS fit of the $Y_{i}^{\mathrm{obs}}$'s on $T_i$'s and $x_i$'s with an intercept. This strategy implicitly assumes away treatment effect heterogeneity, and can lead to inferior properties when $\beta_1 \neq \beta_0$ \citep{freedman08a}. \cite{lin13} chooses $\beta_{1} = \hat{\beta}_1$ and $\beta_0 = \hat{\beta}_0$, the coefficients of the covariates in the OLS fit of $Y_{i}^{\mathrm{obs}}$'s on $x_i$'s with an intercept, in the treatment and control groups, respectively. Numerically, this is identical to the estimator obtained from the regression with interactions discussed before.

\subsection{Our contributions}\label{subsec:contribution}

In practice, it is common to have many covariates. Therefore, it is important to approximate the sampling distribution with $p$ growing with the sample size $n$ at certain rate. Under the finite-population randomization model, \citet{bloniarz16} discussed a high dimensional regime with possibly larger $p$ than $n$ but assumed that the potential outcomes could be well approximated by a sparse linear combination of the covariates, under the ultra sparse regime (termed, for example, by \citet{cai17}) where the number of non-zero coefficients is many fewer than $n^{1/2} / \log p$. Under a super-population framework, \citet{wager16} discussed covariate adjustment using the OLS and some other machine learning techniques.

We study \citet{lin13}'s estimator under the finite-population perspective in the regime where $p < n$ but $p$ grows with $n$ at certain rate.

We focus on this estimator because (a) it is widely used in practice because of its simplicity, and (b) it does not require any tuning parameter unlike other high dimensional or machine learning methods. As in the classic linear regression, asymptotic properties depend crucially on the maximum leverage score $\kappa = \max_{1\leq i\leq n} H_{ii}$, where the $i$-the leverage score $H_{ii}$ is $i$-th diagonal entry of the hat matrix $H = X(X\tran X)^{-1}X\tran $ \citep{huber73}.

Under the regime $\kappa \log p\rightarrow 0$, we prove the consistency of \citet{lin13}'s estimator under mild moment conditions on the population OLS residuals. In the favorable case where all leverage scores are close to their average $p / n$, the consistency holds if $p = o(n / \log n)$.

In addition, we prove that \citet{lin13}'s estimator is asymptotically normal under $\kappa p\rightarrow 0$ and extra mild conditions, with the same variance formula as the fixed-$p$ regime. Furthermore, we proposed a debiased estimator, which is asymptotically normal under an even weaker assumption $\kappa^{2} p \log p\rightarrow 0$, with the same variance as before. In the favorable case where all leverage scores are close to their average $p / n$, \citet{lin13}'s estimator is asymptotically normal when $p = o(n^{1/2})$, but the debiased estimator is asymptotically normal when $p = o(n^{2/3} / (\log n)^{1/3})$. \citet{lin13}'s estimator may also be asymptotically normal in the latter regime, but it requires an extra condition (See Theorem \ref{thm:asym_normality_htau}). In our simulation, the debiased estimator indeed yields better finite-sample inferences.

For statistical inference, we propose several asymptotically conservative variance estimators, which yield valid asymptotic Wald-type confidence intervals for the ATE. We prove the results under the same regime $\kappa \log p \rightarrow 0$ with the same conditions as required for the asymptotic normality.

Importantly, our theory does not require any modeling assumptions on the {\it fixed} potential outcomes and the covariates. It is nonparametric.

We prove novel vector and matrix concentration inequalities for sampling without replacement. These tools are particularly useful for finite population causal inference, and can also complement and potentially enrich the theory in other areas such as survey sampling \citep[e.g.,][]{cochran07}, matrix sketching \citep[e.g.,][]{woodruff14} and transductive learning \citep[e.g.,][]{el09}.

\section{Regression Adjustment}\label{sec:estimators}

\subsection{Point Estimators}\label{sec::pointestimators}

We reformulate \citet{lin13}'s estimator. The ATE is the difference between the two intercepts of the population OLS coefficients in \eqref{eq:mut-ols}:
$$
\tau =  \frac{1}{n} \sum_{i=1}^n Y_i(1) - \frac{1}{n} \sum_{i=1}^n Y_i(0)  =   \mu_{1} - \mu_{0}.
$$
Therefore, we focus on estimating $\mu_1$ and $\mu_0.$
Let $X_{t}\in \R^{n_{t}\times p}$ denote the sub-matrix formed by the rows of $X$, and $Y_{t}^{\mathrm{obs}}\in \R^{n_{t}}$ the subvector of $Y^{\mathrm{obs}}=(Y_1^{\mathrm{obs}},\ldots, Y_n^{\mathrm{obs}} )\tran $, with indices in $\T_{t}$ $(t=0,1)$.
The regression-adjusted estimator follows two steps.
First, for $t\in \{0, 1\}$, we regress $Y_{t}^{\mathrm{obs}}$ on $X_{t}$ with an intercept, and obtain the fitted intercept $\hat{\mu}_{t}\in \R$ and coefficient of the covariate $\hat{\beta}_{t}\in \R^{p}$. Second, we estimate $\tau$ by
\begin{equation}\label{eq:estimator}
\htau = \hat{\mu}_{1} - \hat{\mu}_{0}.
\end{equation}

In general, $\htau $ is biased in finite samples. Correcting the bias gives stronger theoretical guarantees as our later asymptotic analysis suggests. Here we propose a bias-corrected estimator. Define the potential residuals based on the population OLS as
\begin{equation}
  \label{eq:et}
  e(t) = Y(t) - \mu_{t} - X\beta_{t},\quad (t=0,1) .
\end{equation}
The property of the OLS guarantees that $e(t)$ is orthogonal to $\one$ and $X$:
\begin{equation}
  \label{eq:error_key_property}
  \one\tran e(t) = 0,\quad
  X\tran e(t) = 0 ,\quad (t=0,1).
\end{equation}
Let $\hat{e}\in \R^{n}$ be the vector residuals from the sample OLS:
\begin{equation}
  \label{eq:hatei}
  \hat{e}_{i} = \left\{
    \begin{array}{cc}
      Y_{i}^{\mathrm{obs}} - \hat{\mu}_{1} - x_{i}\tran \hat{\beta}_{1},  & (i\in \T_{1}) ,\\
      Y_{i}^{\mathrm{obs}} - \hat{\mu}_{0} - x_{i}\tran \hat{\beta}_{0} , & (i\in \T_{0}) .
    \end{array}
 \right.
\end{equation}
For any vector $\alpha\in \R^{n}$, let $\alpha_{t}$ denote the subvector of $\alpha$ with indices in $\T_{t}$ (e.g. $Y_{t}(1), e_{t}(1), \hat{e}_{t}$, etc.)

Let $H = X(X\tran X)^{-1}X\tran$ be the hat matrix of $X$, and $H_t= X_{t}(X_{t}\tran X_{t})^{-1}X_{t}\tran $ be the hat matrix of $X_{t}$. Let $H_{ii}$ be the $i$-th diagonal element of $H$, also termed as the \emph{leverage score}. Define
\begin{equation}\label{eq:hatDeltat}
\Delta_{t} = \frac{1}{n}\sum_{i=1}^{n}  e_{i}(t)  H_{ii}, \quad \hat{\Delta}_{t} = \frac{1}{n_{t}}\sum_{i\in \T_{t}} \hat{e}_{i}  H_{ii} .
\end{equation}
We introduce the following debiased estimator:
\begin{equation}\label{eq:estimator_debiased}
\hdtau = \htau - \lb \frac{n_{1}}{n_{0}}\hat{\Delta}_{0} - \frac{n_{0}}{n_{1}}\hat{\Delta}_{1}\rb.
\end{equation}
The bias correction terms in \eqref{eq:estimator_debiased} come from higher order asymptotic expansions.
When $p = 1$, \eqref{eq:estimator_debiased} reduces the bias formula in \citet[][Section 6 point (iv)]{lin13}. Thus \eqref{eq:estimator_debiased} is an extension to the multivariate case.

\subsection{Variance estimators}
For fixed $p$, \cite{lin13} proved that $n^{1/2}(\htau - \tau)$ is asymptotically normal with variance
\begin{eqnarray}
\sigma_{n}^{2}  &=& \frac{1}{n_{1}}\sum_{i=1}^{n}e_i^2(1) + \frac{1}{n_{0}}\sum_{i=1}^{n}e_i^2(0) - \frac{1}{n}\sum_{i=1}^{n}(e_{i}(1)- e_{i}(0))^{2}\label{eq:sigman} \\
&=&\sum_{i=1}^{n}\lb \sqrt{\frac{n_{0}}{n_{1}n}}e_{i}(1)+ \sqrt{\frac{n_{1}}{n_{0}n}}e_{i}(0)\rb^{2} .  \label{eq:sigman2}
\end{eqnarray}
The second form \eqref{eq:sigman2} follows from some simple algebra and shows that $\sigma_{n}^{2}$ is always non-negative. The first form \eqref{eq:sigman} motivates conservative variance estimators. The third term in \eqref{eq:sigman} has no consistent estimator without further assumptions on $e(1)$ and $e(0)$. Ignoring it and estimating the first two terms in \eqref{eq:sigman} by their sample analogues, we have the following variance estimator:
\begin{equation}\label{eq:sigman2hat_simple}
\hat{\sigma}^{2} = \frac{n}{n_{1}(n_{1} - 1)}\sum_{i\in \T_{1}} \hat{e}_i^2 + \frac{n}{n_{0}(n_{0} - 1)}\sum_{i\in \T_{0}} \hat{e}_i^2.
\end{equation}

Although \eqref{eq:sigman2hat_simple} appears to be conservative due to the neglect of the third term in \eqref{eq:sigman2}, we find in numerical experiments that it typically underestimates $\sigma_{n}^{2}$ in the cases beyond our theoretic limit with many covariates or many influential observations. The classic linear regression literature suggests rescaling the residual as
\begin{equation}\label{eq:HC}
\td{e}_{i} = \left\{
    \begin{array}{ll}
   \hat{e}_{i} & \text{(HC0)}\\
     \sqrt{\frac{n - 1}{n - p}}\hat{e}_{i} & \text{(HC1 correction)}\\
     \frac{\hat{e}_{i}}{\sqrt{1 - H_{t,ii}}} & \text{(HC2 correction)}\\
     \frac{\hat{e}_{i}}{1 - H_{t,ii}} & \text{(HC3 correction)}\\
    \end{array}
\right., \quad (i\in \T_t)
\end{equation}
where $H_{t,ii}$ is the diagonal element of $H_t $ corresponding to unit $i$. HC0 corresponds to the estimator \eqref{eq:sigman2hat_simple} without corrections. Previous literature has shown that the above corrections, especially HC3, are effective in improving the finite sample performance of variance estimator in linear regression under independent super-population sampling \citep[e.g.][]{mackinnon13, cattaneo2018inference}. More interestingly, it is also beneficial to borrow these HC$j$'s to the context of a completely randomized experiment. This motivates the following variance estimators
\begin{equation}
  \hat{\sigma}_{\HC j}^{2} = \frac{n}{n_{1}(n_{1} - 1)}\sum_{i\in \T_{1}}   \td{e}_{i, j}^{2} + \frac{n}{n_{0}(n_{0} - 1)}\sum_{i\in \T_{0}} \td{e}_{i, j}^{2}
\end{equation}
where $\td{e}_{i, j}$ is the residual in \eqref{eq:HC} with $j$ corresponding to HC$j$ for $j =0,  1, 2, 3$.

Based on normal approximations, we can construct Wald-type confidence intervals for the ATE based on point estimators $\htau$ and $\hdtau$ with estimated standard errors $\hat{\sigma}_{\HC j}$.

\section{Main Results}\label{sec:main}

\subsection{Regularity conditions}\label{subsec::conditions}

We embed the finite population quantities $\{ (x_i, Y_i(1), Y_i(0))\}_{i=1}^n$ into a sequence, and impose regularity conditions on this sequence. The first condition is on the sample sizes.

\begin{assumption}\label{assumption::A1}
$n / n_{1}  = O(1)$ and $n / n_{0} = O(1)$.
\end{assumption}

Assumption \ref{assumption::A1} holds automatically if treatment and control groups have fixed proportions (e.g., $n_1/n=n_0/n=1/2$ for balanced experiments). It is not essential and can be removed at the cost of complicating the statements.

The second condition is on $\kappa = \max_{1\leq i \leq n}H_{ii}$, the maximum leverage score, which also plays a crucial role in classic linear models \citep[e.g.][]{huber73, mammen89, donoho01}.

\begin{assumption}\label{assumption::A2}
$\kappa \log p = o(1)$.
\end{assumption}

The maximum leverage score satisfies
\begin{equation}
  \label{eq:Hiibound}
  p / n = \tr(H)/ n\le \kappa \le \|H\|_{\op} = 1  \Longrightarrow \kappa \in [p / n, 1] .
\end{equation}

Assumption \ref{assumption::A2} permits influential observations as long as $\kappa = o(1 / \log p)$.
In the favorable case where $\kappa = O(p / n)$, it reduces to $p\log p / n\rightarrow 0$, which permits $p$ to grow almost linearly with $n$. Moreover, it implies
\begin{equation}
  \label{eq:pn}
  \frac{p}{n} \le \kappa = o\lb\frac{1}{\log p}\rb = o(1) \Longrightarrow p = o(n).
\end{equation}

Assumptions \ref{assumption::A1} and \ref{assumption::A2} are useful for establishing consistency. The following two extra conditions are useful for variance estimation and asymptotic normality. The third condition is on the correlation between the potential residuals from the population OLS in \eqref{eq:et}.

\begin{assumption}\label{assumption::A3}
There exist a constant $\eta > 0$ independent of $n$ such that
\[
\rho_{e} \triangleq   \frac{   e(1)\tran e(0)    }{  \| e(1) \|_2 \| e(0) \|_2 }    > -1 + \eta .

\]
\end{assumption}

Assumption \ref{assumption::A3} is mild because it is unlikely to have the perfect negative sample correlation between the treatment and control residual potential outcomes in practice.

The fourth condition is on the following two measures of the potential residuals based on the population OLS in \eqref{eq:et}.
\[
\cure_{2} = n^{-1} \max \left\{\|e(0)\|_{2}^{2} , \|e(1)\|_{2}^{2} \right\}, \quad \cure_{\infty} = \max \left\{\|e(0)\|_{\infty}, \|e(1)\|_{\infty}\right\}.
\]

\begin{assumption}\label{assumption::A4}
$\cure_{\infty}^{2} / (n\cure_{2}) = o(1).$
\end{assumption}

Assumption \ref{assumption::A4} is a Lindeberg--Feller-type condition requiring that no single residual dominates others. A similar form appeared in \citet{hajek1960limiting}'s finite population central limit theorem. Previous works require more stringent assumptions on the fourth moment \citep{lin13, bloniarz16}.

\subsection{Discussion of regularity conditions}

Although the above assumptions are about fixed quantities in the finite population, it is helpful to consider the case where the quantities are realizations of random variables. This approach connects the assumptions to more comprehensible conditions on the data generating process. See \cite{portnoy84, portnoy85, lei16} for examples in other contexts. We emphasize that we do not need the assumptions in this subsection for our main theory but use them to aid interpretation. The readers who believe our assumptions to be mild can skip this subsection at first read.

For Assumption \ref{assumption::A2}, we consider the case where $(x_{i})_{i=1}^{n}$ are realizations of i.i.d. random vectors. \cite{anatolyev17} show that under mild conditions each leverage score concentrates around $p/n$. Here we  further consider the magnitude of the maximum leverage score $\kappa$.

\begin{proposition}\label{prop:kappa_iid}
  Let $Z_{i}$ be i.i.d. random vectors in $\R^{p}$ with arbitrary mean. Assume that $Z_{i}$ has independent entries with
$\max_{1\leq j \leq p}\E |Z_{ij} - \E Z_{ij}|^{\delta} \le M = O(1)$
for some $\delta > 2$.
  Define $Z = (Z_{1}\tran , \ldots, Z_{n}\tran )\tran \in \R^{n\times p}$ and $X =\projone  Z$ so that $X$ has centered columns. If $p = O(n^{\gamma})$ for some $\gamma < 1$, then over the randomness of $Z$,
$$
\max_{1\leq i \leq n}\bigg|H_{ii} - \frac{p}{n}\bigg| = O_\P\lb \frac{p^{2 / \min\{\delta, 4\}}}{n^{(\delta - 2) / \delta}} + \frac{p^{3/2}}{n^{3/2}}\rb,
\quad
\kappa = O_\P\lb\frac{p}{n} + \frac{p^{2 / \min\{\delta, 4\}}}{n^{(\delta - 2) / \delta}}\rb.
$$
\end{proposition}

When $\delta > 4$, Proposition \ref{prop:kappa_iid} implies that $\kappa = O_\P(p / n + n^{-(\delta - 4) / 2\delta}(p / n)^{1/2})$. In this case, Assumption \ref{assumption::A2} holds with high probability if $p = O(n^{\gamma})$ for any $\gamma < 1$. In particular, the fixed-$p$ regime corresponds to $\gamma =0.$

The hat matrix of $X$ is invariant to any nonsingular linear transformation of the columns. Consequently, $X$ and $XA$ have the same leverage scores for any invertible $A\in \R^{n\times n}$. Thus we can extend Proposition \ref{prop:kappa_iid} to random matrices with correlated columns in the form of $\projone ZA$. In particular, when $Z_{i} \stackrel{\text{i.i.d.}}{\sim} N(\mu, \id)$ and $A = \Sigma^{1/2}$, $Z_{i}\tran A\stackrel{\text{i.i.d.}}{\sim} N(\Sigma^{1/2}\mu, \Sigma)$. The previous argument implies that Proposition \ref{prop:kappa_iid} holds for $X = \projone ZA$. We will revisit Proposition \ref{prop:kappa_iid} when imposing further conditions on the $H_{ii}$'s and $\kappa$.

For Assumption \ref{assumption::A4}, we consider the case where the $Y_{i}(t)$'s are realizations of i.i.d. random variables, and make a connection with the usual moment conditions. This helps to understand the growth rates of $\cure_{2}$ and $\cure_{\infty}$.

\begin{proposition}\label{prop:eq}
  Let $Y(t)\in \R^{n}$ be a non-constant random vector with i.i.d. entries, and $X$ be any fixed matrix with centered columns. If for some $\delta > 0$,
$
\E |Y_{i}(t) - \E Y_{i}(t)|^{\delta} < \infty
$
for $t=0,1$,
then
\[\cure_{2} = \left\{
    \begin{array}{ll}
      O_\P(1) & (\delta \ge 2)\\
      o_\P(n^{2 / \delta - 1}) & (\delta < 2)
    \end{array}
\right., \quad \cure_{\infty} = O_{\P}(n^{1/\delta}).\]
Furthermore, $\cure_{2}^{-1} = O_\P(1)$ if $Y_{i}(1)$ or $Y_{i}(0)$ is not a constant.
\end{proposition}

When $\delta > 2$, Proposition \ref{prop:eq} implies $\cure_{\infty}^{2} / (n\cure_{2}) = O_\P(n^{2/\delta - 1}) = o_\P(1)$, and thus Assumption \ref{assumption::A4} holds with high probability.
We will revisit Proposition \ref{prop:eq} for the consistency of $\htau$ and $\hdtau$.

\subsection{Asymptotic Expansions}
We derive an asymptotic expansion of $\htau$.

\begin{theorem}\label{thm:expansion}
  Under Assumptions \ref{assumption::A1} and \ref{assumption::A2},
  \begin{align}
   \htau - \tau & = \lb\frac{\oneone \tran e_{1}(1) }{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}}\rb + \lb\frac{n_{1}}{n_{0}}\Delta_{0} - \frac{n_{0}}{n_{1}}\Delta_{1}\rb \nonumber\\
&\quad + O_\P\lb \sqrt{\frac{\cure_{2} \kappa^{2} p \log p }{n}}   + \sqrt{\frac{\cure_{2}\kappa}{n}}\rb.\label{eq:expansion}
  \end{align}
\end{theorem}

The first term in \eqref{eq:expansion} is the difference-in-means estimator of the residual potential outcomes based on the population OLS. The second term is non-standard and behaves as a ``bias,'' which motivates the debiased estimator $\hdtau$ by subtracting its empirical analogue from $  \htau .$

We need to analyze $\Delta_{t}$ and $\hat{\Delta}_{t} - \Delta_{t}$ to simplify Theorem \ref{thm:expansion} and to derive an asymptotic expansion of $\hdtau$. Define
\begin{equation}\label{eq:Delta}
  \Delta = \max\{|\Delta_{1}|, |\Delta_{0}|\}.
\end{equation}
The Cauchy--Schwarz inequality implies
\begin{equation}\label{eq:Deltat}
|\Delta| = \max_{t=0,1}|\Delta_{t}|\le \sqrt{\frac{1}{n}\sum_{i=1}^{n}H_{ii}} \times \sqrt{\max_{t=0,1}\frac{1}{n}\sum_{i=1}^{n}e_{i}^2(t)  H_{ii}} \le \sqrt{\frac{\cure_{2} \kappa p}{n}} .
\end{equation}
This helps us to obtain the following expansions.

\begin{corollary}\label{cor:expansion}
  Under Assumptions \ref{assumption::A1} and \ref{assumption::A2},
  \begin{eqnarray}
&    \htau - \tau &= \frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}} + O_\P\lb \Delta + \sqrt{\frac{\cure_{2}\kappa^{2} p\log p}{n}} + \sqrt{\frac{\cure_{2}\kappa}{n}}\rb \label{eq:expansion_htau}  \\
& &= \frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}} + O_\P\lb \sqrt{\frac{\cure_{2} \kappa p}{n}}\rb ,  \label{eq:expansion_htau2}
\\
 &   \hdtau - \tau &= \frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}} + O_\P\lb \sqrt{\frac{\cure_{2} \kappa^{2} p \log p}{n}} + \sqrt{\frac{\cure_{2}\kappa}{n}}\rb.  \label{eq:expansion_hdtau}
  \end{eqnarray}
\end{corollary}

Expansion \eqref{eq:expansion_htau} follows from \eqref{eq:expansion} and Assumption \ref{assumption::A1}, and \eqref{eq:expansion_htau2} holds because the upper bound in \eqref{eq:Deltat} dominates the third term of \eqref{eq:expansion}. Expansion \eqref{eq:expansion_hdtau} shows that our de-biasing strategy works because $|\hat{\Delta}_{t} - \Delta_{t}|$ is of higher order compared to the third term of \eqref{eq:expansion_hdtau}. These asymptotic expansions in Corollary \ref{cor:expansion} are crucial for our later analysis.

\subsection{Consistency}
Because the first term in \eqref{eq:expansion} is the difference-in-means of the potential residuals, \citet{neyman23} implies that it has mean $0$ and variance $\sigma_n^2/n$. We then use Chebyshev's inequality to obtain
\begin{equation}\label{eq:first_term}
\frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}}
=O_\P\lb   \sqrt{  \frac{ \sigma_n^2}{ n}  } \rb
= O_\P\lb\sqrt{\frac{\cure_{2}}{n}}\rb.
\end{equation}
Coupled with \eqref{eq:first_term} and $\kappa \le 1$, Corollary \ref{cor:expansion} implies that
  \begin{eqnarray*}
   \htau - \tau &=& O_\P\lb\sqrt{\frac{\cure_{2}(\kappa p + 1)}{n}}\rb , \\
     \hdtau - \tau &=&  O_\P\lb\sqrt{\frac{\cure_{2}(\kappa^{2} p \log p + 1)}{n}}\rb.
  \end{eqnarray*}

These expansions immediately imply the following consistency result. We essentially require the right-hand sides of the above two identities go to $0$.

\begin{theorem}\label{thm:consistency}
   Under Assumptions \ref{assumption::A1} and \ref{assumption::A2},
$\htau$ is consistent if
 $
 \cure_{2} = o\lb n / (\kappa p + 1 ) \rb,
 $
and $\hdtau$ is consistent if
$
\cure_{2} = o\lb  n / ( \kappa^{2} p \log p + 1) \rb.
$
\end{theorem}

In the classical fixed-$p$ regime, Theorem \ref{thm:consistency} implies that both $\htau$ and $\hdtau$ are consistent when $\cure_{2} = o(n)$ because $\kappa \le 1$. From Proposition \ref{prop:eq}, the condition $\cure_{2} = o(n)$ corresponds to the existence of finite first moment under a super-population i.i.d sampling. In the favorable case where $\kappa = O(p / n)$, the same condition $\cure_{2} = o(n)$ is sufficient for the consistency of $\htau$ if $p = O(n^{1/2})$ and for the consistency of $\hdtau$ if $p = O(n^{2/3} / (\log n)^{1/3})$. Thus, both estimators are robust to the heavy-tailedness of the potential residuals.

Moreover, when the residuals are not extremely heavy-tailed such that $\cure_{2} = o(n / p)$, Theorem \ref{thm:consistency} implies that both estimators are always consistent, without any further assumption on $\kappa$ (except Assumption \ref{assumption::A2}). The consistency can hold without a uniformly bounded second moment of the potential residuals.

\subsection{Asymptotic normality}
The first term of \eqref{eq:expansion} is the difference-in-means estimator with potential residuals. We can use the classical finite population central limit theorem to show that it is asymptotically normal with mean $0$ and variance $\sigma_n^2/n$.
Therefore, the asymptotic normalities of $\htau$ and $\hdtau$ hold if the the remainders of \eqref{eq:expansion_htau} and \eqref{eq:expansion_hdtau} are asymptotically vanishing after multiplied by $n^{1/2} / \sigma_{n}$. We first consider $\htau$.
\begin{theorem}\label{thm:asym_normality_htau}
  Under Assumptions \ref{assumption::A1}--\ref{assumption::A4}, $ n^{1/2}(\htau - \tau) / \sigma_{n} \stackrel{\textup{d}}{\rightarrow}N(0, 1)$ if
$\kappa^{2} p \log p = o(1)$ and $ n\Delta^{2} = o(\cure_{2}).$
\end{theorem}

Replacing $\Delta$ in Theorem \eqref{thm:asym_normality_htau} by the upper bound $| \Delta | \leq \sqrt{  \cure_{2} p\kappa / n   }$ in \eqref{eq:Deltat}, we obtain the following looser but cleaner result.

\begin{corollary}\label{cor:asym_normality_htau}
  Under Assumptions \ref{assumption::A1}--\ref{assumption::A4}, $ n^{1/2}(\htau - \tau) / \sigma_{n} \stackrel{\textup{d}}{\rightarrow}N(0, 1)$ if $\kappa p = o(1)$.
\end{corollary}

In the favorable case where $\kappa = O(p / n)$, the condition $\kappa p = o(1)$ reduces to $p^{2} / n\rightarrow 0$, i.e., $p = o(n^{1/2}).$
In this case, Corollary \ref{cor:asym_normality_htau} extends \cite{lin13}'s result to $p = o(n^{1/2})$.

The above result can be sharpened if the leverage scores are well-behaved.
In fact, because $e(t)$ has mean zero, we can rewrite $\Delta_{t}$ as
$$
\Delta_{t} =  n^{-1}\sum_{i=1}^{n}  e_{i}(t)   \lb H_{ii} -  p / n \rb .
$$
The Cauchy--Schwarz inequality implies
\[\Delta = \max_{t=0,1}|\Delta_{t}|\le \max_{1\leq i \leq n}\bigg|H_{ii} - \frac{p}{n}\bigg| \times
\max_{t=0,1} \frac{1}{n}\sum_{i=1}^{n}|e_{i}(t)|
\le \max_{1\leq i \leq n}\bigg|H_{ii} - \frac{p}{n}\bigg|\sqrt{\cure_{2}}.\]
Therefore, the condition $\Delta = o(\sqrt{\cure_{2} / n})$ in Theorem \ref{thm:asym_normality_htau} holds whenever
\begin{equation}\label{eq:Hii-p/n}
\max_{1\leq i \leq n}\bigg|H_{ii} - \frac{p}{n}\bigg| = o\lb n^{-1/2}\rb.
\end{equation}
That is, under \eqref{eq:Hii-p/n}, the asymptotic normality of $\htau$ holds when the other condition in Theorem \ref{thm:asym_normality_htau}  holds, i.e.,
$
\kappa^{2} p \log p\rightarrow 0.
$
In the favorable case where $\kappa = O(p / n)$, the condition reduces to
$p^{3}\log p / n^{2}\rightarrow 0$, which further implies $ p = o(n^{2/3} /(\log n)^{1/3}).$
This relaxes the constraint on the dimension to $n^{2/3}$ up to a log-factor. Under $ p = o(n^{2/3} /(\log n)^{1/3})$, we can use Proposition \ref{prop:kappa_iid} to verify that \eqref{eq:Hii-p/n} holds with high probability if entries of $X$ are independent and have finite $12$-th moments.

Although we relaxes the constraint on the dimension, it is not ideal to impose an extra condition on the leverage scores. In contrast, the debiased estimator is asymptotically normal without any further condition.

\begin{theorem}\label{thm:asym_normality_hdtau}
  Under Assumptions \ref{assumption::A1}--\ref{assumption::A4}, $ n^{1/2}(\hdtau - \tau) / \sigma_{n} \stackrel{\textup{d}}{\rightarrow}N(0, 1)$ if $\kappa^{2} p \log p = o(1)$.
\end{theorem}

Therefore, the debiased estimator has better theoretical guarantees. Under $ p = o(n^{2/3} /(\log n)^{1/3})$, we can use
Proposition \ref{prop:kappa_iid} to verify that the condition $\kappa^{2} p \log p = o(1)$ holds if entries of $X$ are independent and have finite $(6+\eps)$-th moments.

\subsection{Variance estimation}\label{subsec:variance}
The variance estimators $\hat{\sigma}_{\HC j}^{2}$'s are all asymptotically equivalent because the correction terms in \eqref{eq:HC} are negligible under our asymptotic regime. We can prove that the $\hat{\sigma}_{\HC j}^{2}$'s for all $j$ are asymptotically conservative estimators of $\sigma_{n}^{2}$.
\begin{theorem}\label{thm:variance_est}
Under Assumptions \ref{assumption::A1}--\ref{assumption::A4}, $\hat{\sigma}_{\HC j}^{2} / \sigma_{n}^{2}\ge 1 + o_\P(1)$ for all $j\in \{0, 1, 2, 3\}$.
\end{theorem}

Therefore, the Wald-type confidence intervals for the ATE are all asymptotically conservative.

\subsection{Related works}

Theoretical analyses under the finite-population randomization model are challenging due to the lack of probability tools. The closest work to ours is \cite{bloniarz16}, which allows $p$ to grow with $n$ and potentially exceed $n$. However, they assume that the potential outcomes have sparse linear representations based on the covariates, and require $s = o (n^{1/2} / \log p) $ where $s$ is a measure of sparsity. Under additional regularities conditions, they show that $\hat{\tau}(\hat{\beta}_1^{\text{lasso}}, \hat{\beta}_0^{\text{lasso}})$ is consistent and asymptotically normal with $(\hat{\beta}_1^{\text{lasso}}, \hat{\beta}_0^{\text{lasso}})$ being the LASSO coefficients of the covariates. Although the LASSO-adjusted estimator can handle ultra-high dimensional case where $p >\!\!> n$, it has three limitations. First, the requirement $s <\!\!< n^{1/2} / \log p$ is stringent. For instance, the PAC-man dataset considered by \cite{bloniarz16} has $n = 1013 $ and $p = 1172$, so the condition reads $s <\!\!< 4.5$, which implicitly imposes a strong sparse modelling assumption.

Second, the penalty level of the LASSO depends on unobserved quantities. Although they use the cross-validation to select the penalty level, the theoretical properties of this procedure is still unclear. Third, their ``restrictive eigenvalue condition'' imposes certain non-singularity on the submatrices of the covariate matrix. However, (submatrices of) the covariate matrix can be ill-conditioned especially when interaction terms are included in practice. In addition, this condition is computationally challenging to check.

Admittedly, our results cannot deal with the case of $p > n$. Nevertheless, we argue that $p<n$ is an important regime in many applications.

\section{Numerical Experiments}\label{sec:experiment}

We perform extensive numerical experiments to confirm and complement our theory. We examine the performance of the estimators $\htau$ and $\hdtau$ as well as the variance estimators $\hat{\sigma}^{2}_{\textup{HC}j}$ for $j = 0, 1, 2, 3$.

\subsection{Data Generating Process}\label{sec::dgp}

We examine the moderate sample performance of the estimators. We set $n = 2000, n_{1} = n\pi_{1}$ for $\pi_{1}\in \{0.2, 0.5\}$ and generate a matrix $\mathcal{X}\in \R^{n\times n}$ with i.i.d. entries from $t(2)$. We keep the matrix fixed. For each exponent $\gamma \in \{0, 0.05, \ldots, 0.75 \}$, we let $p = \lceil n^{\gamma}\rceil$ and take the first $p$ columns of $\mathcal{X}$ as the covariate matrix. In Supplementary Material III, we also simulate $X$ with $N(0,1)$ and $t(1)$ entries and take $X$ from two real datasets.  We select $t(2)$ distribution for presentation because it is neither too idealized as $N(0, 1)$ (where $\kappa\sim p / n$), nor too irregular as $t(1)$. It is helpful to illustrate and complement our theory.

With $X$, we construct the potential outcomes as
\begin{equation}
  \label{eq:simul_potential_outcomes}
  Y(1) = X\beta_{1}^{*} + \sigma_{1}^{*}\eps(1), \quad Y(0) = X\beta_{0}^{*} + \sigma_{0}^{*}\eps(0),
\end{equation}
with $\beta_{1}^{*} =  \beta_{0}^{*} = 0 \in \R^{p}$, $ \sigma_{1}^{*} = \sigma_{0}^{*} = 1$, and $\eps(1), \eps(0)\in \R^{n}$. Note that for given $\eps(1), \eps(0)$ and $X$, both the ATE estimate ($\htau$ or $\hdtau$) and the variance estimate are invariant to the choices of $\beta_{1}^{*}$ and $\beta_{0}^{*}$. Similarly, we generate $(\eps(1), \eps(0))$ as realizations of random vectors with i.i.d. entries from $N(0,1)$, or $t(2)$, or $t(1)$.

Given $X\in \R^{n\times p}$ and potential outcomes $Y(1), Y(0)\in \R^{p}$, we generate 5000 binary vectors $T\in \R^{n}$ with $n_1$ units assigned to treatment. For each assignment vector, we observe half of the potential outcomes.

\subsection{Repeated Sampling Evaluations}\label{subsec:simulation}

Based on the observe data, we obtain two estimates $\htau$ and $\hdtau$, as well as five variance estimates $\hat{\sigma}^{2}_{\textup{HC}j}\,\, (j = 0,1,2,3)$ and $\sigma_{n}^{2}$. Technically, $\sigma_{n}^{2}$ is not an estimate because it is the theoretical asymptotic variance. Below $\hat{\tau}$ can be either $\htau$ or $\hdtau$, and $\hat{\sigma}^2$ can be any of the five estimates.

Let $\hat{\tau}_{1}, \ldots, \hat{\tau}_{5000}$ denote the estimates in 5000 replicates, and $\tau$ denote the true ATE. The empirical relative absolute bias is
$
|  5000^{-1}\sum_{k=1}^{5000}\hat{\tau}_{k} - \tau | / \sigma_{n}.
$

Similarly, let $\hat{\sigma}_{1}^{2}, \ldots, \hat{\sigma}_{5000}^{2}$ denote the variance estimates obtained in 5000 replicates, and $\hat{\sigma}_{*}^{2}$ denote the empirical variance of $(\hat{\tau}_{1}, \ldots, \hat{\tau}_{5000})$. We compute the standard deviation inflation ratio
$
\mathrm{SDR}(\hat{\sigma})  =  5000^{-1}\sum_{k=1}^{5000}  \hat{\sigma}_{k} / \hat{\sigma}_{*} .
$
Note that $\hat{\sigma}_{*}^{2}$ is an unbiased estimate of true sampling variance of $\hat{\tau}$, which can be different from the theoretical asymptotic variance $\sigma_{n}^{2}$.

For each estimate and variance estimate, we compute the $t$-statistic $n^{1/2}(\hat{\tau} - \tau) / \hat{\sigma}$ and the z-score  $n^{1/2}(\hat{\tau} - \tau) / \hat{\sigma}_{*}$. For each $t$-statistic and the z-score, we estimate the empirical $95\%$ coverage by the proportion within $[-1.96, 1.96]$, the $95\%$ quantile range of $N(0,1)$.

In summary, we compute three measures defined above: relative bias, standard deviation inflation ratios, and $95\%$ coverage.  We repeat $50$ times using different random seeds and record the medians of each measure. Fig. \ref{fig::simulation01} summarizes the results. We emphasize that for each experiment, both $X$ and potential outcomes are fixed and the randomness only comes from treatment assignments.

\subsection{Results}\label{subsec:simulation_results}

\begin{figure}
	\centering
	\begin{subfigure}{0.7\textwidth}
		\includegraphics[width=\textwidth]{synthetic_simul_t2_pi2_bias.pdf}
		\caption{ Relative bias of $\hdtau$ and $\htau$.}\label{fig:t2_bias}
	\end{subfigure}
	\begin{subfigure}{0.7\textwidth}
		\includegraphics[width=\textwidth]{synthetic_simul_t2_pi2_sdinflate.pdf}
		\caption{Ratio of standard deviation between five standard deviation estimates, $\sigma_{n}, \hat{\sigma}_{\textup{HC}0}, \hat{\sigma}_{\textup{HC}1}, \hat{\sigma}_{\textup{HC}2}, \hat{\sigma}_{\textup{HC}3}$, and the true standard deviation of $\htau$.}\label{fig:t2_sdinflate}
	\end{subfigure}
		\begin{subfigure}{0.7\textwidth}
		\includegraphics[width=\textwidth]{synthetic_simul_t2_pi2_coverage.pdf}
		\caption{Empirical $95\%$ coverage of $t$-statistics derived from two estimators and four variance estimators (
                  ``theoretical'' for $\sigma_{n}^{2}$, ``HC2'' for $\hat{\sigma}_{\text{HC}2}^{2}$ and ``HC3''  for $\hat{\sigma}_{\text{HC}3}^{2}$)}\label{fig:t2_coverage}
	\end{subfigure}
	\caption{Simulation with $\pi_{1} = 0.2$. $X$ is a realization of a random matrix with i.i.d. $t(2)$ entries, and $e(t)$ is a realization of a random vector with i.i.d. entries from a distribution corresponding to each column. }\label{fig::simulation01}
\end{figure}

\begin{figure}
	\centering
	\begin{subfigure}{0.7\textwidth}
		\includegraphics[width=\textwidth]{synthetic_simul_t2_pi5_bias.pdf}
		\caption{ Relative bias of $\hdtau$ and $\htau$.}\label{fig:t2_bias-2}
	\end{subfigure}
	\begin{subfigure}{0.7\textwidth}
		\includegraphics[width=\textwidth]{synthetic_simul_t2_pi5_sdinflate.pdf}
		\caption{Ratio of standard deviation between five standard deviation estimates, $\sigma_{n}, \hat{\sigma}_{\textup{HC}0}, \hat{\sigma}_{\textup{HC}1}, \hat{\sigma}_{\textup{HC}2}, \hat{\sigma}_{\textup{HC}3}$, and the true standard deviation of $\htau$.}\label{fig:t2_sdinflate-2}
	\end{subfigure}
		\begin{subfigure}{0.7\textwidth}
		\includegraphics[width=\textwidth]{synthetic_simul_t2_pi5_coverage.pdf}
		\caption{Empirical $95\%$ coverage of $t$-statistics derived from two estimators and four variance estimators (
                  ``theoretical'' for $\sigma_{n}^{2}$, ``HC2'' for $\hat{\sigma}_{\text{HC}2}^{2}$ and ``HC3''  for $\hat{\sigma}_{\text{HC}3}^{2}$)}\label{fig:t2_coverage-2}
	\end{subfigure}
	\caption{Simulation. $X$ is a realization of a random matrix with i.i.d. $t(2)$ entries, and $e(t)$ is a realization of a random vector with i.i.d. entries from a distribution corresponding to each column. }\label{fig::simulation01-2}
\end{figure}

From Figures \ref{fig:t2_bias} and \ref{fig:t2_bias-2}, $\hdtau$ does reduce the bias regardless of the distribution of potential outcomes, especially for moderately large $p$. It is noteworthy that the relative bias is too small ($\le 15\%$) to affect coverage.

For standard deviation inflation ratios, we find that the true sampling variances of $\htau$ and $\hdtau$ are almost identical and thus we set the sampling variance of $\htau$ as the baseline variance $\hat{\sigma}_{*}^{2}$. Figures \ref{fig:t2_sdinflate} and \ref{fig:t2_sdinflate-2} shows an interesting phenomenon that the theoretical asymptotic variance $\sigma_{n}^{2}$ tends to underestimate the true sampling variance for large $p$. Corollary \ref{cor:expansion} partially suggests this. The theoretical asymptotic variance is simply the variance of the first term while the finite sample variance also involves the second term and, more importantly, the error term, which can be large in the presence of high dimensional or influential observations. All variance estimators overestimate $\sigma_{n}^{2}$ because they all ignore the third term of $\sigma_{n}^{2}$. However, all estimators, except the HC3 estimator, tend to underestimate the true sampling variance for large $p$. In contrast, the HC3 estimator does not suffer from anti-conservatism.

Figures \ref{fig:t2_sdinflate} and \ref{fig:t2_sdinflate-2} shows that HC0 and HC1 variance estimates lie between the theoretical asymptotic variance and the HC2 variance estimate. For better visualization, we only plot the $95\%$ coverage of $t$-statistics computed from $\sigma_{n}^{2}, \hat{\sigma}_{\textup{HC}2}^{2}$ and $\hat{\sigma}_{\textup{HC}3}^{2}$  in Figures \ref{fig:t2_coverage} and \ref{fig:t2_coverage-2}.
We draw the following conclusions from Figures \ref{fig:t2_coverage} and \ref{fig:t2_coverage-2}. First, as we pointed out previously, the coverage of two ATE estimates are almost identical because the relative bias is small in these scenarios.
Second, as Figures \ref{fig:t2_sdinflate} and \ref{fig:t2_sdinflate-2} suggest, the $t$-statistic with HC3 variance estimate has the best coverage, and it protects the coverage against the increasing dimension. In contrast, the theoretical asymptotic variance and HC$j$ $(j=0,1,2)$ variance estimates yield significantly lower coverage for large $p$. Therefore, we advocate using $\hat{\sigma}_{\textup{HC}3}^{2}$ for variance estimation.

\subsection{Effectiveness of debiasing}\label{subsec:simulation_debiase}
In the aforementioned settings, the debiased estimator yields almost identical inference as the undebiased estimator. This is not surprising because in the above scenarios the potential outcomes are generated from linear models and thus \cite{lin13}'s estimator has bias close to zero.
However, in practice, the potential outcomes might not have prefect linear relationships with the covariates. To illustrate the potential benefits of debiasing, we consider the ``most-biased'' situation which maximizes the ``bias term'', measured as the second term in the expansion \eqref{eq:expansion}. Specifically, we consider the case where $\eps(0) = \eps$ and $\eps(1) = 2\eps$ for some vector $\eps$ that satisfies \eqref{eq:error_key_property} with sample variance $1$. To maximize the bias term, we take $\eps$ as the solution of
\begin{eqnarray}
\label{eq:worst_pout}
  &&\max_{\eps\in \R^{n}}\bigg|\frac{n_{1}}{n_{0}}\Delta_{0} - \frac{n_{0}}{n_{1}}\Delta_{1}\bigg| = \lb\frac{2n_{0}}{n_{1}} - \frac{n_{1}}{n_{0}}\rb\bigg|\sum_{i=1}^{n}H_{ii}\eps_{i}\bigg|,  \\
  &&\mathrm{s.t.}\,\, \|\eps\|_{2}^{2} / n = 1, X\tran \eps = \one\tran \eps = 0. \nonumber
\end{eqnarray}
We give more details of constructing $\eps$ in Section \ref{app:experiments} of Supplementary Material III. From \eqref{eq:worst_pout}, the bias is amplified when the group sizes are unbalanced. Note that this setting essentially assume a non-linear relationship between the potential outcomes and the covariates.

\begin{figure}[h]
	\centering
	\begin{subfigure}{\textwidth}
          \centering
		\includegraphics[width=0.45\textwidth]{synthetic_simul_debiase_t2_pi2_bias.pdf}

		\includegraphics[width=0.45\textwidth]{synthetic_simul_debiase_t2_pi5_bias.pdf}
		\caption{ Relative bias of $\hdtau$ and $\htau$.}\label{fig:t2_bias_worst}
	\end{subfigure}
		\begin{subfigure}{\textwidth}
                  \centering
		\includegraphics[width=0.45\textwidth]{synthetic_simul_debiase_t2_pi2_coverage.pdf}
		\includegraphics[width=0.45\textwidth]{synthetic_simul_debiase_t2_pi5_coverage.pdf}
		\caption{Empirical $95\%$ coverage of $t$-statistics derived from two estimators and two variance estimators (``HC2'' for $\hat{\sigma}_{\text{HC}2}^{2}$ and ``HC3''  for $\hat{\sigma}_{\text{HC}3}^{2}$)}\label{fig:t2_coverage_worst}
	\end{subfigure}
	\caption{Simulation. $X$ is a realization of a random matrix with i.i.d. $t(2)$ entries and $e(t)$ is defined in \eqref{eq:worst_pout}: (Left) $\pi_{1} = 0.2$; (Right) $\pi_{1} = 0.5$. }\label{fig::simulation02}
\end{figure}

We perform simulation detailed in Section \ref{subsec:simulation} based on potential outcomes in \eqref{eq:worst_pout} and report relative bias and coverage to demonstrate the effectiveness of debiasing. To save space, we only report the coverage for $\hat{\sigma}_{\textup{HC}2}^{2}$ and $\hat{\sigma}_{\textup{HC}3}^{2}$. Fig. \ref{fig::simulation02} summarizes the results.

Unlike the previous settings, the relative bias in this setting is large enough to affect the coverage. From Fig. \ref{fig:t2_bias_worst}, as expected, we see that the relative bias is larger when the group sizes are more unbalanced. The debiased estimator reduces a fair proportion of bias in both cases and improves coverage especially when the dimension is high. We provide experimental results in more settings in Supplementary Material III, which confirm the effectiveness of debiasing.

\subsection{Trimming the Design Matrix}\label{subsec:simulation_regularization}

Our theory suggests that $\kappa$ of the design matrix affects the statistical properties of $\htau$ and $\hdtau$. When there are many influential observations in the data, it is beneficial to reduce $\kappa$ before regression adjustment. Because our theory holds even for mis-specified linear models, any preprocessing of $X$ does not affect the consistency and asymptotic normality if the preprocessing does not depend on $T$ or $Y^{\mathrm{obs}}$. This is a feature of our theory. In contrast, trimming is not applicable to the theory under a super-population perspective assuming a correctly specified regression model.

In Section \ref{sec:experiment}, the entries of $X$ are realizations of heavy-tailed random variables, and $\kappa $ increases even with an infrequent extreme covariate value. For the 50 design matrices used in Section \ref{sec:experiment} with $p = \lceil n^{2/3}\rceil$ and $n = 2000$, the average of $\kappa$ is $0.9558$ with standard error $0.0384$. Now we consider a simple form of trimming which thresholds each column at its 2.5\% and 97.5\% quantiles. Then the average of $\kappa$ reduces dramatically to $0.0704$ with standard error $0.0212$. Fig. \ref{fig::simulation03} shows the coverage of the $t$-statistics derived from $\hdtau$ with and without the trimming. It is clear that the coverage gets drastically improved after trimming.

Since the main goal of this paper is not on statistical methodology, we only propose the above heuristic approach to illustrate the idea of trimming, motivated by our asymptotic theory. The general methodology is an interesting future research topic.

\begin{figure}
	\centering
		\includegraphics[width=0.45\textwidth]{synthetic_simul_thresh_t2_pi2_coverage.pdf}
		\includegraphics[width=0.45\textwidth]{synthetic_simul_thresh_t2_pi5_coverage.pdf}
	\caption{Simulation. Empirical $95\%$ coverage of $t$-statistics derived from the debiased estimator with and without trimming the covariate matrix: (Left) $\pi_{1} = 0.2$; (Right) $\pi_{1} = 0.5$. $X$ is a realization of a random matrix with i.i.d. $t(2)$ entries and $e(t)$ is defined in \eqref{eq:worst_pout}.}\label{fig::simulation03}
\end{figure}

\section{Conclusions and Practical Suggestions}\label{sec:conclusion}

\citet{fisher37} advocated using the analysis of covariance under treatment-unit additivity. \citet{freedman08a} highlighted its dangers under treatment effect heterogeneity. \citet{lin13} proposed a simple OLS estimator with treatment-covariate interactions accounting for potential heterogeneity. We establish the consistency and the asymptotic normality of \citet{lin13}'s estimator allowing for a growing dimension of the covariates. We further propose a debiased estimator which permits valid inference in broader asymptotic regimes.

In summary, we find that the classical inferential procedure tends to be invalid when the design matrix has many covariates or many influential observations. In these scenarios, the bias blows up and the variance estimation becomes anti-conservative. We suggest using the debiased estimator \eqref{eq:estimator_debiased} and the HC3 variance estimator for inference. In addition, we suggest trimming the design matrix to reduce the extreme leverage scores before regression adjustment.

\section{Technical Lemmas}\label{subsec:technicallemmas}

\subsection{Some general results for sampling without replacement}\label{subsec:lemmas1}

Completely randomized experiments have deep connections with sampling without replacement because the treatment and control groups are simple random samples from a finite population of $n$ units. Below we use $\T$ to denote a random size-$m$ subset of $  \{1, \ldots, n\}$ over all $\binom{n}{m}$ subsets, and $\S^{p-1} = \{  (\omega_1, \ldots, \omega_p)\tran :  \omega_1^2+\cdots + \omega_p^2 =1 \}$ to denote the $(p-1)$-dimensional unit sphere.

The first lemma gives the mean and variance of the sample total from sampling without replacement. See \citet[][Theorem 2.2]{cochran07} for a proof.

\begin{lemma}\label{lem:var_sum}
Let $(w_1, \ldots, w_n)$ be fixed scalars with mean $\bar{w} = n^{-1}\sum_{i=1}^{n}w_i$. Then $\sum_{i\in \T} w_i$ has mean $m\bar{w}$ and variance
\[\Var\lb\sum_{i\in \T} w_i  \rb = \frac{m(n-m)}{n(n-1)}\sum_{i=1}^{n}(w_i - \bar{w})^{2}.\]
\end{lemma}

The second lemma gives the Berry--Esseen-type bound for the finite population central limit theorem. See \citet{bikelis69} and \citet{hoglund78} for proofs.

\begin{lemma}\label{prop:finite_CLT}
  Let $(w_{1}, \ldots, w_{n})$ be fixed scalars with
$\bar{w} = n^{-1}\sum_{i=1}^{n}w_{i}$ and $\quad S_w^{2} = \sum_{i=1}^{n}(w_{i} - \bar{w})^{2}.$
Let $m = n f$ for some $f \in (0, 1)$. Then
\begin{align*}
d_{\textup{K}}\lb \frac{\sum_{i\in \T  }(w_{i} - \bar{w})}{ S_w \sqrt{f (1 - f)}}, N(0, 1)\rb
&\le \frac{C}{\sqrt{f(1-f)}}\frac{\sum_{i=1}^{n}(w_{i} - \bar{w})^{2}}{S_w^{3}}\\
& \le \frac{C}{\sqrt{f(1-f)}}\frac{\max_{1\leq i\leq n}|w_{i} - \bar{w}|}{S_w},
\end{align*}
where $d_\textup{K}$ denotes the Kolmogorov distance between two distributions, and $C$ is a universal constant.
\end{lemma}

The following two lemmas give novel vector and matrix concentration inequalities for sampling without replacement.

\begin{lemma}\label{thm:vector_concentration}
Let $(u_{1}, \ldots, u_{n})$ be a finite population of $p$-dimensional vectors with
$\sum_{i=1}^{n}u_{i} = 0.$
Then for any $\delta\in (0, 1)$, with probability $1 - \delta$
\[\left\|\sum_{i\in \T}u_{i}\right\|_{2}\le \|U\|_{F}\sqrt{\frac{m(n - m)}{n(n - 1)}} + \|U\|_{\op}\sqrt{8\log \frac{1}{\delta}}\]
where $u_{i}\tran$ is the $i$-th row of the matrix $U\in\R^{n\times p}$.
\end{lemma}

\begin{lemma}\label{thm:matrix_concentration}
Let $(V_{1}, \ldots, V_{n})$ be a finite population of $(p\times p)$-dimensional Hermittian matrices with
$\sum_{i=1}^{n}V_{i} = 0.$
Let $C(p) = 4(1 + \lceil 2\log p\rceil)$, and
\[\nu^{2} = \left\|\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\right\|_{\op}, \quad \nu_{-}^{2} = \sup_{\omega\in \S^{p-1}}\frac{1}{n}\sum_{i=1}^{n}(\omega\tran V_{i}\omega)^{2}, \quad \nu_{+} = \max_{1\leq i \leq n}\|V_{i}\|_{\op}. \]
Then for any $\delta \in (0, 1)$, with probability $1 - \delta$,
\[\left\|\sum_{i\in \T}V_{i}\right\|_{\op}\le \sqrt{nC(p)}\nu + C(p)\nu_{+} + \sqrt{8n\log \frac{1}{\delta}}\,\nu_{-}.\]
\end{lemma}

The following lemma gives the mean and variance of the summation over randomly selected rows and columns from a deterministic matrix $Q \in \R^{n\times n}.$

\begin{lemma}\label{thm:mean}
Let $Q\in \R^{n\times n}$ be a deterministic matrix, and $Q_{\T} \equiv \sum_{i,j\in \T} Q_{ij}$. Assume $n\ge 4$. Then
 \[\E Q_{\T} = \frac{m(n - m)}{n(n - 1)}\tr(Q) + \frac{m(m - 1)}{n(n - 1)} \one\tran Q\one.\]
 If $Q$ further satisfies $
\one\tran Q = Q\one = 0,
$
then
\[\Var(Q_{\T}) \le \frac{m(n - m)}{n(n - 1)}\|Q\|_{F}^{2}.\]
\end{lemma}

Lemmas \ref{thm:vector_concentration}--\ref{thm:mean} are critical for our proofs. The proofs are relegated to Supplementary Material I. They are novel tools to the best of our knowledge and potentially useful in other contexts such as survey sampling, matrix sketching, and transductive learning.

\subsection{Some results particularly useful for our setting} \label{subsec:lemmas2}

We first give an implication of Assumption \ref{assumption::A3}, a lower bound on $\sigma_{n}^{2}$ under Assumption \ref{assumption::A1}.
\begin{lemma}\label{lem:A3}
  Under Assumptions \ref{assumption::A1} and \ref{assumption::A3},
$
\sigma_{n}^{2} \ge \eta\min\left\{ n_{1} / n_{0},  n_{0} / n_{1} \right\}\cure_{2}.
$
\end{lemma}

Recall $H_t= X_{t}(X_{t}\tran X_{t})^{-1}X_{t}\tran $ and define $\Sigma_t = n_t^{-1} X_{t}\tran X_{t}$ $(t=0,1)$.
The following explicit formula is the starting point of our proof.
\begin{lemma}\label{lem:tauhat}
We have
\begin{equation}
  \label{eq:ATE}
  \htau - \tau = \frac{\oneone \tran e_{1}(1) / n_{1} - \oneone \tran H_1 e_{1}(1) / n_{1}}{1 - \oneone \tran H_1 \oneone  / n_{1}} - \frac{\onezero \tran e_{0}(0) / n_{0} - \onezero \tran H_0e_{0}(0) / n_{0}}{1 - \onezero \tran H_0\onezero  / n_{0}}.
\end{equation}
\end{lemma}

The quantities $\mu_{t}$, $e(t)$, and our estimators $(\htau, \hdtau)$ are all invariant if $X$ is transformed to $X Z$ for any full rank matrix $Z \in \R^{p\times p}$, provided that \eqref{eq:centered} holds. Thus, without loss of generality, we assume
\begin{equation}
  \label{eq:orthonormal}
n^{-1} X\tran X = \id .
\end{equation}
Otherwise, suppose $X$ has the singular value decomposition $U\Sigma V\tran $ with $U\in \R^{n\times p}, \Sigma, V\in \R^{p\times p}$, then we can replace $X$ by $n^{1/2}U = X(n^{1/2}V \Sigma^{-1})$ to ensure \eqref{eq:orthonormal}. We can verify that the key properties in \eqref{eq:error_key_property} still hold. Assuming \eqref{eq:orthonormal}, we can rewrite the hat matrix and the leverage scores as
\begin{equation}\label{eq:hatmatrix}
H = n^{-1} XX\tran , \quad H_{ii} =  n^{-1}   \|x_{i}\|_{2}^{2} ,\quad H_{ij} = n^{-1} x_i \tran x_j .
\end{equation}
Note that the invariance property under the standardization \eqref{eq:orthonormal} is a feature of the OLS-based regression adjustment. It does not hold for many other estimators \citep[e.g.,][]{bloniarz16, wager16}.

We will repeatedly use the following results to obtain the stochastic orders of the terms in \eqref{eq:ATE}. They are consequences of Lemmas \ref{thm:vector_concentration} and \ref{thm:matrix_concentration}.

\begin{lemma}\label{lem:vector_concentration}
Under Assumption \ref{assumption::A1}, for $t = 0, 1$,
\begin{gather*}
\frac{\onet\tran e_t(t)  }{n_{t}} = O_\P\lb \sqrt{\frac{\cure_{2}}{n}}\rb, \quad \left\|\frac{X_{t}\tran \onet }{n_{t}}\right\|_{2} = O_\P\lb \sqrt{\frac{p}{n}}\rb ,  \\
\left\|\frac{X_{t}\tran e_{t}(t)}{n_{t}}\right\|_{2} = O_\P\lb \sqrt{\cure_{2}\kappa}\rb .
\end{gather*}
\end{lemma}

\begin{lemma}\label{lem:matrix_concentration}
Under Assumptions \ref{assumption::A1}, \ref{assumption::A2} and \eqref{eq:orthonormal}, for $t = 0, 1$,
\begin{gather*}
\left\|\Sigma_t - \id \right\|_{\op} = O_\P\lb \sqrt{\kappa\log p}\rb, \quad \left\|\Sigma_t^{-1}\right\|_{\op} = O_\P(1), \\
\left\|\Sigma_t^{-1} - \id \right\|_{\op} = O_\P\lb\sqrt{\kappa\log p}\rb.
\end{gather*}
\end{lemma}

The following lemma states some key properties of an intermediate quantity, which will facilitate our proofs.
\begin{lemma}\label{lemma:Mt}
Define $Q(t)   = H\diag\lb e(t)\rb = \lb H_{ij}e_{j}(t)    \rb_{i, j = 1}^{n}$. It satisfies
\begin{eqnarray*}
\one\tran Q(t) =0, \quad Q(t) \one =0,\quad
\one\tran Q(t) \one = 0,  \\
 \textup{tr}(Q(t) ) =n  \Delta_t,\quad    \|Q(t) \|_{F}^{2} = \sum_{i=1}^{n}e_{i}^2(t)   H_{ii} \le n \cure_{2}\kappa .
\end{eqnarray*}
\end{lemma}

\section{Proofs of the main results}
\subsection{Proof of the asymptotic expansions}\label{subsec:proof_expansion}

\begin{proof}[Proof of Theorem \ref{thm:expansion}]
We need to analyze the terms in \eqref{eq:ATE}. First, by Lemmas \ref{lem:vector_concentration} and \ref{lem:matrix_concentration},
\begin{align*}
  \frac{\onet \tran H_t \onet }{n_{t}} &= \frac{\onet \tran X_{t}}{n_{t}}\Sigma_t^{-1}\frac{X_{t}\tran \onet }{n_{t}}
  \leq    \left\| \Sigma_t^{-1}  \right\|_{\op}  \left \|  \frac{X_{t}\tran \onet }{n_{t}}  \right \|_2^2
  = O_\P\lb \frac{p}{n}\rb.
\end{align*}
Using \eqref{eq:pn} that $p = o(n)$, we obtain that
\begin{equation}\label{eq:expansion1}
\frac{1}{1 - \onet \tran H_t \onet  / n_{t}} = 1 + O_\P\lb\frac{p}{n}\rb.
\end{equation}
Second,
\begin{align}
  &\frac{\onet \tran H_t e_{t}(t)}{n_{t}} = \frac{\onet \tran X_{t}}{n_{t}}\Sigma_t^{-1}\frac{X_{t}\tran e_{t}(t)}{n_{t}}\nonumber\\
= & \frac{\onet \tran X_{t}}{n_{t}}\frac{X_{t}\tran e_{t}(t)}{n_{t}}
+ \frac{\onet \tran X_{t}}{n_{t}}\lb \Sigma_t^{-1} - \id \rb\frac{X_{t}\tran e_{t}(t)}{n_{t}} \equiv   R_{t1} + R_{t2} . \label{eq:R1+R2}
\end{align}
Note that here we do not use the naive bound for $\onet \tran H_t e_{t}(t) / n_t$ as for $\onet \tran H_t \onet / n_t$ in \eqref{eq:expansion1} because this gives weaker results. Instead, we bound $R_{t1}$ and $ R_{t2}$ separately. Lemmas \ref{lem:vector_concentration} and \ref{lem:matrix_concentration} imply
\begin{align}\label{eq:R2}
  R_{t2}  \le \left\|  \Sigma_t^{-1} - \id  \right\|_\op
  \left\| \frac{X_{t}\tran \onet }{n_{t}} \right\|_2
  \left\| \frac{X_{t}\tran e_{t}(t)}{n_{t}} \right\|_2
  = O_\P\lb \sqrt{\frac{\cure_{2}\kappa^{2} p \log p}{n}}\rb.
\end{align}
We apply Chebyshev's inequality to obtain that
\begin{equation}
  \label{eq:Markov}
  R_{t1} = \E R_{t1} + O_\P\lb\sqrt{\Var(R_{t1})}\rb.
\end{equation}
Therefore, to bound $R_{t1}$, we need to calculate its first two moments.
Recalling \eqref{eq:hatmatrix} and the definition of $Q(t)$ in Lemma \ref{lemma:Mt}, we have
\begin{align}
  R_{t1} &= \frac{1}{n_{t}^{2}}\lb\sum_{i\in \T_{t}}x_{i}\tran \rb\lb \sum_{i\in \T_{t}} x_{i}e_{i}(t)\rb
  = \frac{1}{n_{t}^{2}} \sum_{i\in \T_{t}} \sum_{j\in \T_{t}} x_{i}\tran x_{j}e_{j}(t)  \nonumber\\
& = \frac{1}{n_{t}^{2}} \sum_{i\in \T_{t}} \sum_{j\in \T_{t}} nH_{ij}e_{j}(t)   = \frac{n}{n_{t}^{2}}\sum_{i\in \T_{t}} \sum_{j\in \T_{t}} Q_{ij}(t) .\label{eq:R1}
\end{align}
Lemmas \ref{thm:mean} and \ref{lemma:Mt} imply the expectation of $R_{t1}$:
\begin{align}
  \E R_{t1} &= \frac{n}{n_{t}^{2}}\lb\frac{n_{1}{n_{0}}}{n(n - 1)} \tr\lb Q(t) \rb + \frac{n_{t}(n_{t} - 1)}{n (n - 1)}\one\tran  Q(t) \one \rb\nonumber\\
& = \frac{n n_{1}n_{0}}{n_{t}^{2}(n - 1)}\Delta_{t} = \frac{n_{1}n_{0}}{n_{t}^{2}}\Delta_{t} + O\lb\frac{|\Delta_{t}|}{n}\rb.\label{eq:ER1}
\end{align}
We then bound the variance of $ R_{t1} $:
\begin{eqnarray}
  \Var(R_{t1}) &=& \frac{n^{2}}{n_{t}^{4}}\Var\lb\sum_{i,j\in \T_{t}}Q_{ij}(t) \rb
\leq    \frac{n^{2}}{n_{t}^{4}}   \frac{n_1 n_0}{  n(n-1) } \| Q(t) \|_F^2  \label{eq::varianceofQ} \\
&\leq &  \frac{n^{2}}{n_{t}^{4}}   \frac{n_1 n_0}{  n(n-1) } n\mathcal{E}_2 \kappa
=  O\lb\frac{\cure_{2}\kappa}{n}\rb, \label{eq:VarR1}
\end{eqnarray}
 where \eqref{eq::varianceofQ} follows from Lemma \ref{thm:mean}, \eqref{eq:VarR1} follows from Lemma \ref{lemma:Mt} and Assumption \ref{assumption::A1}.
Putting \eqref{eq:R1+R2}--\eqref{eq:ER1} and \eqref{eq:VarR1} together, we obtain that
\begin{align}
  \frac{\onet \tran H_t e_{t}(t)}{n_{t}} &= \frac{n_{1}n_{0}}{n_{t}^{2}}\Delta_{t} + O_\P\lb \sqrt{\frac{\cure_{2}\kappa^{2} p \log p}{n}} + \frac{|\Delta_{t}|}{n} + \sqrt{\frac{\cure_{2}\kappa}{n}}\rb\label{eq:annoyingterm1}
\end{align}
By \eqref{eq:Deltat} and \eqref{eq:pn} that $p = o(n)$, \eqref{eq:annoyingterm1} further simplifies to
\begin{equation}
  \label{eq:annoyingterm2}
  \frac{\onet \tran H_t e_{t}(t)}{n_{t}} = \frac{n_{1}n_{0}}{n_{t}^{2}}\Delta_{t} + O_\P\lb \sqrt{\frac{\cure_{2}\kappa^{2} p\log p}{n}} + \sqrt{\frac{\cure_{2}\kappa}{n}}\rb.
\end{equation}
Using Lemma \ref{lem:vector_concentration}, \eqref{eq:annoyingterm2}, and the fact that $\kappa \le 1$, we have
\begin{align}
  &\frac{\onet \tran e_{t}(t)}{n_{t}} - \frac{\onet \tran H_t  e_{t}(t)}{n_{t}}= O_\P\lb \sqrt{\frac{\cure_{2}}{n}} + \Delta + \sqrt{\frac{\cure_{2}\kappa^{2} p\log p}{n}}\rb.\label{eq:annoyingterm3}
\end{align}

Finally, putting \eqref{eq:expansion1}, \eqref{eq:annoyingterm2} and \eqref{eq:annoyingterm3} together into \eqref{eq:ATE}, we obtain that
\begin{eqnarray}
  \htau - \tau  &=& \lb \frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\oneone \tran H_1 e_{1}(1)}{n_{1}}\rb \lb 1 + O_\P\lb\frac{p}{n}\rb\rb
  \nonumber \\
& &\qquad  - \lb \frac{\onezero \tran e_{0}(0)}{n_{0}} - \frac{\onezero \tran H_0e_{0}(0)}{n_{0}}\rb \lb 1 + O_\P\lb\frac{p}{n}\rb\rb
\nonumber\\
&  =& \frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}} + \frac{\onezero \tran H_0e_{0}(0)}{n_{0}} - \frac{\oneone \tran H_1 e_{1}(1)}{n_{1}}
\nonumber\\
&& \,\, + O_\P\lb \sqrt{\frac{p^{2}\cure_{2}}{n^{3}}} + \frac{p\Delta}{n} + \sqrt{\frac{\cure_{2}\kappa^{2} p^{3}\log p}{n^{3}}}\rb
\nonumber\\
& =& \frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}} + \frac{n_{1}}{n_{0}}\Delta_{0} - \frac{n_{0}}{n_{1}}\Delta_{1}
\nonumber\\
&& \,\, + O_\P\lb \sqrt{\frac{p^{2}\cure_{2}}{n^{3}}} + \frac{p\Delta}{n} + \sqrt{\frac{\cure_{2}\kappa^{2} p\log p}{n}}   + \sqrt{\frac{\cure_{2}\kappa}{n}}\rb. \label{eq::Opterm}
\end{eqnarray}
where \eqref{eq::Opterm} uses \eqref{eq:pn} that $p = o(n)$. The fourth terms dominates the first term in \eqref{eq::Opterm} because $p = o(n)$ and $\kappa\ge p / n$. The third term dominates the second term in \eqref{eq::Opterm}  because, by \eqref{eq:Deltat},
\[\frac{p\Delta}{n} \le \kappa \Delta \le \sqrt{\kappa}\Delta = O\lb \sqrt{\frac{\cure_{2}\kappa^{2} p}{n}}\rb.\]
Deleting the first two terms in \eqref{eq::Opterm}, we complete the proof.
\end{proof}

\begin{proof}[Proof of Corollary \ref{cor:expansion}]
Assumption \ref{assumption::A1} implies
$
  \frac{n_{1}}{n_{0}}\Delta_{0} - \frac{n_{0}}{n_{1}}\Delta_{1} = O\lb \Delta\rb,
$
which, coupled with Theorem \ref{thm:expansion}, implies
\eqref{eq:expansion_htau}.

The key is to prove the result for the debiased estimator. By definition,
\begin{align*}
   \hdtau - \tau & = \frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}} + \frac{n_{1}}{n_{0}}(\Delta_{0} - \hat{\Delta}_{0}) - \frac{n_{0}}{n_{1}}(\Delta_{1} - \hat{\Delta}_{1})\\
&\quad + O_\P\lb \sqrt{\frac{\cure_{2}\kappa^{2} p\log p}{n}}  + \sqrt{\frac{\cure_{2}\kappa}{n}}\rb,
\end{align*}
and therefore, the key is to bound $|\Delta_{t} - \hat{\Delta}_{t}|$.

We introduce an intermediate quantity
$
\td{\Delta}_{t} =  n_{t}^{-1}  \sum_{i\in \T_{t}}H_{ii}e_{i}(t).
$
It has mean
$
\E \td{\Delta}_{t} = \Delta_{t}
$
and variance
\begin{equation}
  \label{eq:VartdDeltat}
\Var(\td{\Delta}_{t}) \le  \frac{1}{n_{t}^{2}}\frac{n_{1}n_{0}}{n(n - 1)}\sum_{i=1}^{n}H_{ii}^{2}e_{i}^2(t)  \le \frac{\cure_{2}\kappa^{2}}{n_{t}^{2}} = O\lb\frac{\cure_{2}\kappa^{2}}{n}\rb,
\end{equation}
from Lemma \ref{lem:var_sum} and Assumption \ref{assumption::A1}. Equipped with the first two moments, we use Chebyshev's inequality to obtain
\begin{equation}
  \label{eq:tdDeltat-Deltat}
  |\td{\Delta}_{t} - \Delta_{t}| = O_\P\lb \sqrt{\frac{\cure_{2}\kappa^{2}}{n}}\rb.
\end{equation}
Next we bound $|\hat{\Delta}_{t} - \td{\Delta}_{t}|$. The Cauchy--Schwarz inequality implies
\begin{equation}
  \label{eq:hatDeltat-tdDeltat1}
  |\hat{\Delta}_{t} - \td{\Delta}_{t}| \le \frac{1}{n_{t}}\sum_{i\in \T_{t}} H_{ii}|\hat{e}_{i} - e_{i}(t)|\le \sqrt{\frac{1}{n_{t}}\sum_{i\in \T_{t}}H_{ii}^{2}}\sqrt{\frac{1}{n_{t}}\sum_{i\in \T_{t}}(\hat{e}_{i} - e_{i}(t))^{2}}.
\end{equation}
First,
\begin{equation}
  \label{eq:hatDeltat-tdDeltat3}
  \frac{1}{n_{t}}\sum_{i\in \T_{t}}H_{ii}^{2}\le \frac{n\kappa}{n_{t}}\lb\frac{1}{n}\sum_{i=1}^{n}H_{ii}\rb = O\lb\frac{\kappa p}{n}\rb.
\end{equation}
Second, using the fact
$
\hat{e}_{t} = (\id - H_{t})e_{t}(t),
$
we have
\begin{align}
  &\frac{1}{n_{t}}\sum_{i\in \T_{t}}(\hat{e}_{i} - e_{i}(t))^{2} = \frac{1}{n_{t}}\|\hat{e}_{t} - e_{t}(t)\|_{2}^{2} =
  \frac{1}{n_t}e_{t}\tran (t)  H_t e_{t}(t)\nonumber\\
 = & \lb\frac{X_{t}\tran e_{t}(t)}{n_{t}}\rb\tran \Sigma_t^{-1}\frac{X_{t}\tran e_{t}(t)}{n_{t}}\nonumber\\
\le & \left\|\Sigma_t\right\|_{\op}^{-1}\left\|\frac{X_{t}\tran e_{t}(t)}{n_{t}}\right\|_{2}^{2}= O_\P(\cure_{2}\kappa),\label{eq:hatDeltat-tdDeltat2}
\end{align}
where the last line follows from Lemma \ref{lem:vector_concentration}.
Putting \eqref{eq:hatDeltat-tdDeltat3} and  \eqref{eq:hatDeltat-tdDeltat2} into \eqref{eq:hatDeltat-tdDeltat1}, we obtain
\begin{equation}
  \label{eq:hatDeltat-tdDeltat}
  |\hat{\Delta}_{t} - \td{\Delta}_{t}| = O_\P\lb \sqrt{\frac{\cure_{2}\kappa^{2} p}{n}}\rb.
\end{equation}
Combining \eqref{eq:tdDeltat-Deltat} and \eqref{eq:hatDeltat-tdDeltat} together, we have
$
|\hat{\Delta}_{t} - \Delta_{t}| = O_\P\lb \sqrt{ \cure_{2}\kappa^{2} p / n} \rb.
$
We complete the proof by invoking Theorem \ref{thm:expansion}.
\end{proof}

\subsection{Proof of asymptotic normality}

\begin{proof}[Proofs of Theorems \ref{thm:asym_normality_htau} and \ref{thm:asym_normality_hdtau}]
We first prove the asymptotic normality of the first term in the expansions:
\begin{eqnarray}
\label{eq::cltfirstterm}
\frac{n^{1/2}}{\sigma_{n}}\lb \frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}}\rb\stackrel{\text{d}}{\rightarrow}N(0, 1).
\end{eqnarray}
Recalling $0 = \one\tran e(0) =  \one \tran e_1(0)  + \one \tran  e_0(0) $, we obtain that
\begin{eqnarray}
&&  \frac{n^{1/2}}{n_{1}}\oneone \tran e_{1}(1) - \frac{n^{1/2}}{n_{0}}\onezero \tran e_{0}(0)
=  \frac{n^{1/2}}{n_{1}}\oneone \tran e_{1}(1) + \frac{n^{1/2}}{n_{0}}\oneone \tran e_{1}(0) \nonumber  \\
&=& \sum_{i\in \T_{t}}\lb \frac{n^{1/2}}{n_{1}} e_{i}(1)+ \frac{n^{1/2}}{n_{0}}e_{i}(0)\rb
\equiv   \sum_{i\in \T_{t}} w_i, \label{eq::srsrepresentation}
\end{eqnarray}
where $w_{i} = \frac{n^{1/2}}{n_{1}} e_{i}(1)+ \frac{ n^{1/2} }{n_{0}}e_{i}(0)$. Based on \eqref{eq:sigman2}, we can verify that
\begin{align*}
&S_w^{2} \equiv  \sum_{i=1}^{n}(w_{i} - \bar{w})^{2} = \sum_{i=1}^{n}w_{i}^{2}
=  n \sum_{i=1}^{n}  \lb \frac{e_{i}(1)}{n_{1}} + \frac{ e_{i}(0) }{n_{0}} \rb^2
= \frac{n^{2}}{n_{1}n_{0}}\sigma_{n}^{2}.
\end{align*}
Applying Lemma \ref{prop:finite_CLT} to the representation \eqref{eq::srsrepresentation}, we have
\[d_{\text{K}}\lb \frac{n^{1/2}}{\sigma_{n}}\lb \frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}}\rb, N(0, 1)\rb = O\lb \frac{\max_{1\leq i \leq n}|w_{i}|}{S_w}\rb.\]
Lemma \ref{lem:A3} and Assumption \ref{assumption::A4} imply
\[S_w^{-1} = O\lb \cure_{2}^{-1 / 2}\rb, \quad \max_{1\leq i \leq n}|w_{i}| = O\lb\frac{\cure_{\infty}}{n^{1/2}}\rb = o\lb\cure_{2}^{1 / 2}\rb.\]
Therefore, \eqref{eq::cltfirstterm} holds because convergence in Kolmogorov distance implies weak convergence.

We then prove the asymptotic normalities of the two estimators.
Corollary \ref{cor:expansion} and Lemma \ref{lem:A3} imply
  \begin{align*}
  &\frac{n^{1/2}(\htau - \tau)}{\sigma_{n}}\\
= &\frac{n^{1/2}}{\sigma_{n}}\lb\frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}}\rb + O_\P\lb \frac{\sqrt{\cure_{2}\kappa^{2} p\log p}}{\sigma_{n}} + \frac{n^{1/2}\Delta}{\sigma_{n}} + \frac{\sqrt{\cure_{2}\kappa}}{\sigma_{n}}\rb\\
= & \frac{n^{1/2}}{\sigma_{n}}\lb\frac{\oneone \tran e_{1}(1)}{n_{1}} - \frac{\onezero \tran e_{0}(0)}{n_{0}}\rb + O_\P\lb \sqrt{\kappa^{2} p\log p} + \sqrt{\frac{n}{\cure_{2}}}\Delta + \sqrt{\kappa}\rb.
  \end{align*}
  We complete the proof  by noting that $\kappa = o(1)$ in \eqref{eq:pn} under Assumption \ref{assumption::A2}. The same proof carries over to $\hdtau$.
\end{proof}

\subsection{Proof of asymptotic conservatism of variance estimators}\label{subsec:proof_var}

\begin{proof}[Proof of Theorem \ref{thm:variance_est}]
First, we prove the result for $j = 0$. Recalling $\hat{e}_{t} = (\id  - H_t )e_{t}(t)$, we have
\begin{align}
&\frac{1}{n_{t}}\sum_{i\in \T_{t}}\hat{e}_{i}^{2} = \frac{1}{n_{t}} e_{t}(t) \tran   (\id  - H_t )e_{t}(t)\nonumber\\
= & \frac{1}{n_{t}}\sum_{i\in \T_{t}}e_{i}^2(t)
 -
 \lb\frac{X_{t}\tran e_{t}(t)}{n_{t}}\rb\tran
 \Sigma_t^{-1}
 \frac{X_{t}\tran e_{t}(t)}{n_{t}}
 \triangleq S_{t1}   - S_{t2}  .\label{eq:S1-S2}
\end{align}
Lemma \ref{lem:vector_concentration} and the fact $\kappa = o(1)$ in \eqref{eq:pn} together imply a bound for $S_{t2}$:
\begin{equation}
  \label{eq:S2t}
  S_{t2} \leq  \left\|   \Sigma_t^{-1} \right\|_\op  \left\|  \frac{X_{t}\tran e_{t}(t)}{n_{t}} \right\|_2^2
  = O_\P\lb \cure_{2}\kappa\rb = o_\P\lb \cure_{2} \rb.
\end{equation}
The first term, $S_{t1}$, has mean $\E S_{t1}   = n^{-1}\sum_{i=1}^{n}e_{i}^2(t)$ and variance
\begin{eqnarray}
  \Var(S_{t1}  ) & \le & \frac{1}{n_{t}^{2}}\frac{n_{1}n_{0}}{n(n  - 1)}\sum_{i=1}^{n}e_{i}^4(t) \label{eq::variance-St1} \\
  & \le& \frac{n}{n_{t}^{2}}\cure_{\infty}^{2} \cure_{2} = O\lb\frac{\cure_{\infty}^{2}\cure_{2}}{n}\rb \label{eq:VarS1t} \\
  &=& o_\P(\cure_{2}^{2}), \label{eq::finalboundvarSt1}
\end{eqnarray}
where \eqref{eq::variance-St1} follows from Lemma \ref{lem:var_sum}, \eqref{eq:VarS1t} follows from the definitions of $\cure_{2}$ and $\cure_{\infty}$ and Assumption \ref{assumption::A1}, and \eqref{eq::finalboundvarSt1} follows from Assumption \ref{assumption::A4} that $\cure_{\infty}^{2} = o(n\cure_{2}  )$. Therefore, applying Chebyshev's inequality, we obtain
\begin{equation}\label{eq:SMarkov}
S_{t1}   = \E S_{t1}   + O_\P\lb\sqrt{\Var(S_{t1}  )}\rb
= \frac{1}{n}\sum_{i=1}^{n}e_{i}^2(t) + o_\P(\cure_{2}) .
\end{equation}
Combining the bounds for $S_{t1} $ in \eqref{eq:SMarkov} and $S_{t2} $ in \eqref{eq:S2t}, we have
\begin{equation}\label{eq:avgeit}
\frac{1}{n_{t}}\sum_{i\in \T_{t}}\hat{e}_{i}^{2} = \frac{1}{n}\sum_{i=1}^{n}e_{i}^2(t)   + o_\P\lb \cure_{2}\rb.
\end{equation}
Using the formula of $\hat{\sigma}^{2}$ in \eqref{eq:sigman2hat_simple} and Assumption \ref{assumption::A1}, we have
\begin{align*}
  \hat{\sigma}_{\HC 0}^{2}
&= \frac{n}{n_1 - 1} \lb     \frac{1}{n}\sum_{i=1}^{n}e_{i}^2(1)   + o_\P\lb \cure_{2}\rb \rb
+ \frac{n}{n_0 - 1} \lb \frac{1}{n}\sum_{i=1}^{n}e_{i}^2(0)   + o_\P\lb \cure_{2}\rb \rb \\
  & = \frac{1}{n_{1}}\sum_{i=1}^{n} e_i^2(1) + \frac{1}{n_{0}}\sum_{i=1}^{n}e_i^2(0) + o_\P(\cure_{2}).
\end{align*}
Using the formula of $\sigma_n^2$ in \eqref{eq:sigman}, we have
$$
 \hat{\sigma}_{\HC 0}^{2}  \ge \sigma_{n}^{2} + \frac{1}{n}\sum_{i=1}^{n}(e_{i}(1)- e_{i}(0))^{2} + o_\P(\cure_{2}) \ge \sigma_{n}^{2} + o_\P(\cure_{2}),
$$
which, coupled with Lemma \ref{lem:A3}, implies that $ \hat{\sigma}_{\HC 0}^{2} /  \sigma_{n}^{2} \geq 1+o_\P(1) .$

Next we prove that the $\hat{\sigma}^{2}_{\HC j}$'s are asymptotically equivalent. It suffices to show that
\begin{equation}\label{eq:ehat_equivalence}
\min_{j=1,2,3}\min_{1\leq i \leq n }|\td{e}_{i, j}| / |\hat{e}_{i}| = 1 + o_\P(1).
\end{equation}
The proof for $j = 1$ follows from $p/n = o(1)$ in \eqref{eq:pn}. To prove \eqref{eq:ehat_equivalence} for $j=2,3$, we need to prove that
$
\max_{t=0,1} \max_{i\in \T_{t}}H_{t,ii} = o_\P(1).
$
This follows from Lemma \ref{lem:matrix_concentration} and the fact that $\kappa = o(1)$ in \eqref{eq:pn}:
\[
\max_{i\in \T_{t}}H_{t,ii} = \max_{i\in \T_{t}}  n_{t}^{-1}x_{i}\tran \Sigma_t^{-1}x_{i}
= O_\P\left(  n_t^{-1} \max_{1\leq i \leq n}\|x_{i}\|_{2}^{2} \right) = O_\P\lb\kappa\rb.
\]
\end{proof}

\section*{Acknowledgment}
We thank Luke Miratrix at Harvard for helpful suggestions. Two reviewers gave us critical comments which improved the quality of our paper.

\section*{Supplementary Materials}
Supplementary Material I gives the proofs of the lemmas in Section \ref{subsec:technicallemmas}. Supplementary Material II gives the proofs of the propositions in Section \ref{subsec::conditions}. Supplementary Material III gives extensive numerical experiments. We post the programs to replicate all the experimental results at \url{https://github.com/lihualei71/RegAdjNeymanRubin/}.

\bibliography{high_dim_Neyman_Rubin}
\bibliographystyle{plainnat}

\newpage

\appendix

\renewcommand {\theequation} {S\arabic{equation}}
\renewcommand {\thefigure} {S\arabic{figure}}

\begin{center}
  \begin{LARGE}
    Supplementary Material I: Technical Tools
  \end{LARGE}
\end{center}

\section{Concentration Inequalities for Sampling Without Replacement}\label{sec:concentration}

\subsection{Some existing tools}
The proofs rely on concentration inequalities for sampling without replacement. \citet[][Theorem 4]{hoeffding63} proved the following result that sampling without replacement is more concentrated in convex ordering than i.i.d. sampling.

\begin{proposition}\label{prop:hoeffding}
Let $C = (c_{1}, \ldots, c_{n})$ be a finite population of fixed elements. Let $Z_{1}, \ldots, Z_{m}$ be a random sample with replacement from $C$ and $W_{1}, \ldots, W_{m}$ be a random sample without replacement from $C$. If the function $f(x)$ is continuous and convex, then
\[\E f\lb\sum_{i=1}^{m}Z_{i}\rb\ge \E f\lb\sum_{i=1}^{m}W_{i}\rb.\]
\end{proposition}

From Proposition \ref{prop:hoeffding}, most concentration inequalities for independent sampling carry over to sampling without replacement. Later a line of works, in different contexts, showed an even more surprising phenomenon that sampling without replacement can have strictly better concentration than independent sampling \citep[e.g.,][]{serfling74, diaconis87, lee98, bobkov04, cortes09, el09, bardenet15, tolstikhin17}. In particular, \citet[][Theorem 9]{tolstikhin17} proved a useful concentration inequality for the empirical processes for sampling without replacement.

\begin{proposition}\label{prop:emp_process}
Let $C = (c_{1}, \ldots, c_{n})$ be a finite population of fixed elements, and $W_{1}, \ldots, W_{m}$ be a random sample without replacement from $C$. Let $\F$ be a class of functions on $C$, and
\[S(\F) = \sup_{f\in \F} \sum_{i=1}^{m}f(W_{i}),\quad
\nu(\F)^{2} = \sup_{f\in \F}\Var(f(W_{1})). \]
Then
\[\P(S(\F) - \E[S(\F)]\ge t)\le \exps{-\frac{(n + 2)t^{2}}{8n^{2}\nu(\F)^{2}}} . \]
\end{proposition}

Proposition \ref{prop:emp_process} gives a sub-gaussian tail of $S(\F)$ with the sub-gaussian parameter depending solely on the variance. In contrast, the concentration inequalities in the standard empirical process theory for independent sampling usually requires the functions in $\F$ to be uniformly bounded and the tail is either sub-gaussian with the sub-gaussian parameter being the uniform bound on $\F$ or sub-exponential with Bernstein-style behaviors; see \cite{boucheron13} for instance. Therefore, Proposition \ref{prop:emp_process} provides a more precise statement that sampling without replacement is more concentrated than independent sampling for a large class of statistics.

We need the following result from \citet[][Theorem 5.1.(2)]{tropp16} to prove the matrix concentration inequality.
\begin{proposition}\label{prop:norm_mean}
Let $\td{V}_{1}, \ldots, \td{V}_{m}$ be independent $p\times p$ random matrices with $\E \td{V}_{i} = 0$ for all $i$. Let $C(p) = 4(1 + \lceil 2\log p\rceil)$. Then
\[\lb\E \left\|\sum_{i=1}^{n}\td{V}_{i}\right\|_{\op}^{2}\rb^{\frac{1}{2}}\le \sqrt{C(p)} \left\|\sum_{i=1}^{n}\E \td{V}_{i}^{2}\right\|_{\op}^{\frac{1}{2}} + C(p)\lb \E \max_{1\leq i \leq n}\|\td{V}_{i}\|_{\op}^{2}\rb^{\frac{1}{2}}.\]
\end{proposition}

We will also use the facts that for any $u\in \R^{p}$ and Hermitian $V\in \R^{p\times p}$,
\[\|u\|_{2} = \sup_{\omega\in \S^{p-1}}u\tran \omega, \quad \|V\|_{\op} = \sup_{\omega\in \S^{p-1}}\omega\tran V\omega.\]

\subsection{Proofs of Lemmas \ref{thm:vector_concentration} and \ref{thm:matrix_concentration}}

\begin{proof}[Proof of Lemma \ref{thm:vector_concentration}]
Let
\[C = (u_{1}, \ldots, u_{n}), \quad \text{and}\quad \F = \{f_{\omega}(u) = u\tran \omega: \omega\in \S^{p-1}\}.\]
Let $u$ be a vector that is randomly sampled from $C$. Then
\begin{align*}
  \nu^{2}(\F) &= \sup_{\omega\in \S^{p-1}}\Var(u\tran \omega) \le \sup_{\omega\in \S^{p-1}}\E(u\tran \omega)^{2} \\
& = \sup_{\omega\in \S^{p-1}}  \frac{1}{n}\sum_{i=1}^{n}(u_{i}\tran \omega)^{2}
= \sup_{\omega\in \S^{p-1}}  \omega\tran \lb\frac{1}{n}\sum_{i=1}^{n}u_{i}u_{i}\tran \rb\omega\\
& = \sup_{\omega\in \S^{p-1}} \omega\tran \lb\frac{U\tran U}{n}\rb\omega \le \frac{\|U\|_{\op}^{2}}{n}.
\end{align*}
By Proposition \ref{prop:emp_process},
\begin{align*}
  \P\lb \left\|\sum_{i\in \T}u_{i}\right\|_{2} \ge \E\left\|\sum_{i\in \T}u_{i}\right\|_{2} + t \rb\le \exps{-\frac{(n + 2)t^{2}}{8n \|U\|_{\op}^{2}}}\le \exps{-\frac{t^{2}}{8\|U\|_{\op}^{2}}},
\end{align*}
or, equivalently, with probability $1 - \delta$,
\begin{equation}
  \label{eq:vector1}
  \left\|\sum_{i\in \T}u_{i}\right\|_{2} \le \E\left\|\sum_{i\in \T}u_{i}\right\|_{2} + \|U\|_{\op}\, \sqrt{8\log \frac{1}{\delta}}.
\end{equation}
By the Cauchy--Schwarz inequality,
\[\lb\E\left\|\sum_{i\in \T}u_{i}\right\|_{2}\rb^{2}\le \E\left\|\sum_{i\in \T}u_{i}\right\|_{2}^{2} = \sum_{j=1}^{p}\E\lb\sum_{i\in \T}u_{ij}\rb^{2}.\]
Lemma \ref{lem:var_sum} implies
\[\E\lb\sum_{i\in \T}u_{ij}\rb^{2} = \frac{m(n - m)}{n(n - 1)}\sum_{i=1}^{n}u_{ij}^{2}.\]
As a result,
\begin{equation}\label{eq:vector2}
\lb\E\left\|\sum_{i\in \T}u_{i}\right\|_{2}\rb^{2}\le \frac{m(n - m)}{n(n - 1)}\sum_{i=1}^{n}\|u_{i}\|_{2}^{2} = \|U\|_{F}^{2}\,\frac{m(n - m)}{n(n - 1)}.
\end{equation}
We complete the proof by using \eqref{eq:vector1} and \eqref{eq:vector2}.
\end{proof}

\begin{proof}[Proof of Lemma \ref{thm:matrix_concentration}]
Let
\[C = (V_{1}, \ldots, V_{n}), \quad \text{and}\quad \F = \{f_{\omega}(V) = \omega\tran V\omega: \omega\in \S^{p-1}\}.\]
Let $V$ be a vector that is randomly sampled from $C$. Then
$$
  \nu^{2}(\F) = \sup_{\omega\in \S^{p-1}}\Var(\omega\tran V\omega) \le \sup_{\omega\in \S^{p-1}}\E(\omega\tran V\omega)^{2}
 =  \sup_{\omega\in \S^{p-1}} \frac{1}{n}\sum_{i=1}^{n}(\omega\tran V_{i}\omega)^{2} = \nu^{2}_{-}.
$$
By Proposition \ref{prop:emp_process},
\begin{align*}
  \P\lb \left\|\sum_{i\in \T}V_{i}\right\|_{\op} \ge \E\left\|\sum_{i\in \T}V_{i}\right\|_{\op} + t \rb\le \exps{-\frac{(n + 2)t^{2}}{8n^{2}\nu^{2}_{-}}}\le \exps{-\frac{t^{2}}{8n \nu^{2}_{-}}},
\end{align*}
or, equivalently, with probability $1 - \delta$,
\begin{equation}
  \label{eq:matrix1}
  \left\|\sum_{i\in \T}V_{i}\right\|_{\op} \le \E\left\|\sum_{i\in \T}V_{i}\right\|_{\op} + \sqrt{8n\log \frac{1}{\delta}}\,\nu_{-}.
\end{equation}

We then bound $\E\left\|\sum_{i\in \T}V_{i}\right\|_{\op} .$
Let $\td{V}_{1}, \ldots, \td{V}_{m}$ be an i.i.d. random sample with replacement from $C$. We have
\begin{align}
  &\E\left\|\sum_{i\in \T}V_{i}\right\|_{\op} \le \E\left\|\sum_{i=1}^{m}\td{V}_{i}\right\|_{\op}
  \le \ \lb  \E\left\|\sum_{i=1}^{m}\td{V}_{i}\right\|_{\op}^{2}\rb^{\frac{1}{2}} \le\sqrt{nC(p)}\nu + C(p)\nu_{+},\label{eq:matrix2}
\end{align}
where the first inequality follows from Proposition \ref{prop:hoeffding} due to the convexity of $\|\cdot\|_{\op}$, the second inequality follows from the Cauchy--Schwarz inequality, and the third inequality follows from Proposition \ref{prop:norm_mean}.

Combining \eqref{eq:matrix1} and \eqref{eq:matrix2}, we complete the proof.
\end{proof}

\section{Mean and Variance of the Sum of Random Rows and Columns of a Matrix}\label{app:mean_variance}

We give a full proof of Lemma \ref{thm:mean}. When $m = 0$ or $m = n$, $Q_{\T}$ is deterministic with zero variance and the inequality holds automatically. Thus we assume $1\le m\le n-1$.

Let $\sum_{[i_1, \ldots, i_k]}$ denote the sum over all $(i_1, \ldots, i_k)$ with mutually distinct elements in $\{1,\ldots, n\}$. For instance, $\sum_{[i, j]}$ denotes the sum over all pairs $(i, j)$ with $i\not = j$.
We first state a basic result for sampling without replacement.

\begin{lemma}\label{lem:trivial_comb}
  Let $i_{1}, \ldots, i_{k}$ be distinct indices in $\{1,\ldots, n  \}$ and $\T$ be a uniformly random subset of $\{ 1,\ldots, n\}$ with size $m$. Then
\[\P\lb i_{1}, \ldots, i_{k}\in \T\rb = \frac{m\cdots (m - k + 1)}{n\cdots (n - k + 1)}.\]
\end{lemma}

  By definition,
\begin{equation}\label{eq:QT}
Q_{\T} = \sum_{i=1}^{n}Q_{ii}I(i\in \T) + \sum_{[i, j]} Q_{ij}I(i,j\in \T).
\end{equation}
The mean of $Q_\T$ follows directly from Lemma \ref{lem:trivial_comb}:
\begin{align*}
\E Q_{\T} &= \sum_{i=1}^{n}Q_{ii}\cdot \frac{m}{n} + \sum_{[i, j]} Q_{ij}\cdot \frac{m(m-1)}{n(n-1)} \\
& = \frac{m(n - m)}{n(n - 1)}\tr(Q) + \frac{m(m - 1)}{n(n - 1)}(\one\tran Q\one).
\end{align*}

The rest of this section proves the result of the variance.
Let
\begin{align*}
&  c_1 = \frac{m(n - m)}{n(n - 1)}, \quad  c_2 = \Var\lb I(1, 2\in \T)\rb = c_1 \frac{(m-1)(n + m - 1)}{n(n-1)},\\
&  c_3 = \Cov\lb I(1, 2\in \T), I(1, 3\in \T)\rb = c_1\frac{(m - 1)(mn - 2m - 2n + 2)}{n(n-1)(n-2)},\\
&  c_4 = \Cov\lb I(1, 2\in \T), I(3, 4\in \T)\rb = c_1\frac{(m - 1)(-4mn + 6n + 6m - 6)}{n(n-1)(n - 2)(n - 3)},\\
& c_5 = \Cov\lb I(1 \in \T), I(1, 2\in \T)\rb = c_1\frac{m - 1}{n},\\
& c_6 = \Cov\lb I(1 \in \T), I(2, 3\in \T)\rb = -c_1\frac{2(m - 1)}{n(n - 2)}.
\end{align*}
Using \eqref{eq:QT}, we have
\begin{align}\label{eq:var_decompose}
\Var(Q_{\T})& =  \underbrace{  \Var\lb \sum_{i=1}^{n}Q_{ii}I(i\in \T)\rb }_{V_\textup{I}}
+ \underbrace{  \Var\lb \sum_{[i, j]} Q_{ij}I(i, j\in \T)\rb}_{V_\textup{II}}\\
& \quad + 2  \underbrace{  \Cov\lb \sum_{i=1}^{n}Q_{ii}I(i\in \T), \sum_{[i, j]} Q_{ij}I(i, j\in \T)   \rb  }_{V_\textup{III}}.  \nonumber
\end{align}

The next subsection deals with the three terms in \eqref{eq:var_decompose}, separately.

\subsection{Simplifying  \eqref{eq:var_decompose}}

\paragraph{Term $V_\textup{I}$}
Lemma \ref{lem:var_sum} implies
\begin{align}
V_\textup{I}=  \Var\lb \sum_{i=1}^{n}Q_{ii}I(i\in \T)\rb& = \frac{m(n - m)}{n(n-1)}\sum_{i=1}^{n}\lb Q_{ii} - \frac{1}{n}\sum_{i=1}^{n}Q_{ii}\rb^{2} \nonumber\\
& = c_1\sum_{i=1}^{n}Q_{ii}^{2} - \frac{c_1}{n}(\tr(Q))^{2}. \label{eq:var_decompose1}
\end{align}

\paragraph{Term $V_\textup{II}$}
We expand $V_\textup{II}$ as
\begin{align*}
  &V_\textup{II} = \Var\lb \sum_{[i, j]}Q_{ij}I(i, j\in \T)\rb = \Cov\lb \sum_{[i, j]}Q_{ij}I(i, j\in \T), \sum_{[i', j']}Q_{i'j'}I(i', j'\in \T)\rb \\
&  =\sum_{[i, j]} \lb Q_{ij}^{2} + Q_{ij}Q_{ji}\rb\Var(I(i, j\in \T)) + \sum_{[i, j, k, \ell]}Q_{ij}Q_{k\ell} \Cov\lb I(i, j\in \T), I(k, \ell\in \T)\rb  \\
& \quad + \sum_{[i, j, k]}\lb Q_{ij}Q_{ik} + Q_{ij}Q_{ki}\rb \Cov\lb I(i, j\in \T), I(i, k\in \T)\rb  \\
& \quad + \sum_{[i, j, k]}\lb Q_{ij}Q_{jk} + Q_{ij}Q_{kj}\rb \Cov\lb I(i, j\in \T), I(j, k\in \T)\rb  \\
 &   =  c_2 \sum_{[i, j]}\lb Q_{ij}^{2} + Q_{ij}Q_{ji}\rb + c_4 \sum_{[i, j, k, \ell]}Q_{ij}Q_{k\ell} \\
& \quad  + c_3\sum_{[i, j, k]}\lb Q_{ij}Q_{ik} + Q_{ij}Q_{ki} + Q_{ij}Q_{jk} + Q_{ij}Q_{kj}\rb .
\end{align*}
We then reduce the summation over $[i,j,k,l]$ to summations over fewer indices.
First,
\begin{align*}
  \lb \sum_{[i, j]} Q_{ij}\rb^{2} &= \sum_{[i, j]}\lb Q_{ij}^{2} + Q_{ij}Q_{ji}\rb + \sum_{[i, j, k, \ell]}Q_{ij}Q_{k\ell}\\
& \quad + \sum_{[i, j, k]}\lb Q_{ij}Q_{ik} + Q_{ij}Q_{ki} + Q_{ij}Q_{jk} + Q_{ij}Q_{kj}\rb.
\end{align*}
Second, $\one\tran Q\one = 0$ implies
$\sum_{[i, j]} Q_{ij} = -\sum_{i=1}^{n}Q_{ii} = -\tr(Q),$
which further implies
\begin{align*}
\sum_{[i, j, k, \ell]}Q_{ij}Q_{k\ell} & = \lb\tr(Q)\rb^{2} - \sum_{[i, j]}\lb Q_{ij}^{2} + Q_{ij}Q_{ji}\rb \\
& \quad  - \sum_{[i, j, k]}\lb Q_{ij}Q_{ik} + Q_{ij}Q_{ki} + Q_{ij}Q_{jk} + Q_{ij}Q_{kj}\rb.
\end{align*}
The above two facts simplify $V_{\textup{II}}$ to
\begin{align}
V_\textup{II} = & \,\, c_{4}(\tr(Q))^{2} + (c_2 - c_4)\sum_{[i, j]}\lb Q_{ij}^{2} + Q_{ij}Q_{ji}\rb \nonumber\\
 + &(c_3 - c_4)\sum_{[i, j, k]}\lb Q_{ij}Q_{ik} + Q_{ij}Q_{ki} + Q_{ij}Q_{jk} + Q_{ij}Q_{kj}\rb\label{eq:var_decompose2_3}.
\end{align}
We then reduce the summation over $[i,j,k]$ to summations over fewer indices.
Note that $\one\tran Q = Q\one = 0$ implies
$  \sum_{j=1}^{n}Q_{ij} = \sum_{i=1}^{n}Q_{ij} = 0,$
which further implies
\begin{align*}
  \sum_{[i, j, k]} Q_{ij}Q_{ik} &= \sum_{[i, j]}Q_{ij}\sum_{k\not= i,j}Q_{ik} = -\sum_{[i, j]}Q_{ij}(Q_{ii} + Q_{ij})\\
& = -\sum_{i=1}^nQ_{ii}\sum_{j\not = i}Q_{ij} - \sum_{[i, j]}Q_{ij}^{2} = \sum_{i=1}^{n}Q_{ii}^{2} - \sum_{[i, j]}Q_{ij}^{2}.
\end{align*}
Similarly,
\[\sum_{[i, j, k]} Q_{ij}Q_{kj} = \sum_{i=1}^{n}Q_{ii}^{2} - \sum_{[i, j]}Q_{ij}^{2},\]
\[\sum_{[i, j, k]} Q_{ij}Q_{ki} = \sum_{[i, j, k]} Q_{ij}Q_{jk} = \sum_{i=1}^{n}Q_{ii}^{2} - \sum_{[i, j]}Q_{ij}Q_{ji}.\]
Using the above three identities to simplify the third term in \eqref{eq:var_decompose2_3}, we obtain
\begin{align}
V_\textup{II}  =  c_{4}(\tr(Q))^{2} + 4(c_3 - c_4) \sum_{i=1}^n Q_{ii}^{2}
 + (c_2 - 2c_3 + c_4)\sum_{[i, j]}\lb Q_{ij}^{2} + Q_{ij}Q_{ji}\rb . \label{eq:var_decompose2_4}
\end{align}

\paragraph{Term $V_\textup{III}$}
The covariance term is
\begin{eqnarray*}
V_\textup{III} &=   &\Cov\lb \sum_{i=1}^{n}Q_{ii}I(i\in \T), \sum_{[i, j]} Q_{ij}I(i, j\in \T)\rb \\
 & =& \sum_{[i, j]}Q_{ii}(Q_{ij} + Q_{ji})\Cov\lb I(i\in \T), I(i, j\in \T)\rb  \\
& & \qquad + \sum_{[i, j, k]}Q_{ii}Q_{jk}\Cov\lb I(i\in \T), I(j, k\in \T)\rb \\
&  =& c_5 \sum_{[i, j]}Q_{ii}(Q_{ij} + Q_{ji}) + c_6 \sum_{[i, j, k]}Q_{ii}Q_{jk}.
\end{eqnarray*}
Similar to previous arguments,
\begin{eqnarray*}
  \sum_{[i, j]}Q_{ii}(Q_{ij} + Q_{ji}) & = & \sum_{i=1}^nQ_{ii}\sum_{j\not = i}(Q_{ij} + Q_{ji}) = -2\sum_{i=1}^nQ_{ii}^{2} ,\\
  \sum_{[i, j, k]}Q_{ii}Q_{jk} & =& \sum_{[i, j]}Q_{ii}\sum_{k\not = i,j}Q_{jk} = -\sum_{[i, j]}Q_{ii}(Q_{jj} + Q_{ji})\\
& =& -\sum_{i=1}^nQ_{ii} \sum_{j\not = i}(Q_{jj} + Q_{ji}) = -\sum_{i=1}^nQ_{ii} \lb \tr(Q) - Q_{ii} - Q_{ii}\rb\\
& =& -(\tr(Q))^{2} + 2\sum_{i=1}^nQ_{ii}^{2}.
\end{eqnarray*}
Using the above two identities, we can simplify $V_\textup{III}$ to
\begin{align}
 V_\textup{III}   = -c_6 (\tr(Q))^{2} - 2(c_5 - c_6)\sum_{i=1}^nQ_{ii}^{2}.  \label{eq:var_decompose3_2}
\end{align}

Putting \eqref{eq:var_decompose1}, \eqref{eq:var_decompose2_4} and \eqref{eq:var_decompose3_2} together, we obtain that
\begin{align}
\Var(Q_\T) & =   \underbrace{  (c_1 + 4c_3 - 4c_4 - 4c_5 + 4c_6)}_{C_\textup{I}}  \sum_{i=1}^{n}Q_{ii}^{2}
+  \underbrace{   \lb c_4 - \frac{c_1}{n} - 2c_6\rb  }_{C_\textup{II}} (\tr(Q))^{2}\nonumber\\
& \quad + \underbrace{   (c_2 - 2c_3 + c_4)}_{C_\textup{III}}  \sum_{[i, j]}(Q_{ij}^{2} + Q_{ij}Q_{ji}).\label{eq:exact_var}
\end{align}
We simplify \eqref{eq:exact_var} in the next subsection by deriving bounds for the coefficients.

\subsection{Bounding the coefficients $C_\textup{I}, C_\textup{II}$ and $C_\textup{III}$ in \eqref{eq:exact_var}}
\paragraph{Bounding $C_\textup{I}$}
We have
\begin{align}
&C_\textup{I} = c_1 + 4c_3 - 4c_4 - 4c_5 + 4c_6\nonumber\\
 = & c_1 + 4c_1\frac{m - 1}{n}\lb \frac{mn - 2m - 2n + 2}{(n - 1)(n - 2)} + \frac{4mn - 6m - 6n + 6}{(n - 1)(n - 2)(n - 3)} - 1 - \frac{2}{n - 2}\rb.\nonumber
\end{align}
Through tedious calculation, we obtain that
\begin{align*}
  &\frac{mn - 2m - 2n + 2}{(n - 1)(n - 2)} + \frac{4mn - 6m - 6n + 6}{(n - 1)(n - 2)(n - 3)} - 1 - \frac{2}{n - 2} = -\frac{(n - m - 1)n}{(n - 2)(n - 3)}.
\end{align*}
Thus,
$
  C_\textup{I} = c_1\lb 1 - \frac{4(m - 1)(n - m - 1)}{(n - 2)(n - 3)}\rb.
$

\paragraph{Bounding $C_\textup{II}$}
We have
\begin{align*}
C_\textup{II} = & \,\, c_4 - \frac{c_1}{n} - 2c_6 = -\frac{c_1}{n} + c_1\frac{m - 1}{n(n - 2)}\lb \frac{-4mn + 6m + 6n - 6}{(n - 1)(n - 3)} + 4\rb\\
= & -\frac{c_1}{n} + c_1\frac{(m - 1)(4n^2 - 4mn + 6m - 10n + 6)}{n(n - 1)(n - 2)(n - 3)}\\
= & -\frac{c_1}{n}\lb 1 - \frac{(m - 1)(n - m - 1)(4n - 6)}{(n - 1)(n - 2)(n - 3)}\rb\\
\le & \,\, c_1 \frac{(m - 1)(n - m - 1)(4n - 6)}{n(n - 1)(n - 2)(n - 3)} \le  \,\, \frac{c_1}{n} \frac{4(m - 1)(n - m - 1)}{n(n - 2)(n - 3)}.
\end{align*}

\paragraph{Bounding $C_\textup{III}$}
We consider four cases.

\noindent $\bullet$ If $m = 1$, then $c_2 = c_3 = c_4 = 0$ and
$C_{\textup{III}}\le \frac{c_{1}}{2}.$

\noindent $\bullet$ If $m = 2$, then
\begin{align*}
  C_{\textup{III}} &= c_1\lb\frac{n + 1}{n(n - 1)} - \frac{-4}{n(n - 1)(n - 2)} - \frac{2}{n(n - 1)(n - 2)}\rb\\
& = c_1\lb \frac{n + 1}{n(n - 1)} + \frac{2}{n(n - 1)(n - 2)}\rb \le \frac{c_1}{2}.
\end{align*}

\noindent $\bullet$ If $m = 3$, then
\begin{align*}
  C_{\textup{III}} &= c_1\lb\frac{2(n + 2)}{n(n - 1)} - \frac{4(n - 4)}{n(n - 1)(n - 2)} - \frac{12n - 24}{n(n - 1)(n - 2)(n - 3)}\rb\\
& = c_1\lb\frac{2(n + 2)}{n(n - 1)} - \frac{4(n - 4)}{n(n - 1)(n - 2)} - \frac{12}{n(n - 1)(n - 3)}\rb.
\end{align*}
If $n\ge 7$,
\[C_{\textup{III}}\le c_1\frac{2(n + 2)}{n (n - 1)}\le \frac{c_1}{2}.\]
For $n = 4, 5, 6$, we can also verify that
$C_{\textup{III}}\le\frac{c_1}{2}$.

\noindent $\bullet$ If $m \ge 4$, then
\[4mn - 6m - 6n + 6 = (2m - 6)(n - 3) + 2(mn - 6)\ge (2m + 2)(n - 3).\]
and thus
\[c_{4}\le c_{1}\frac{(2m + 2)(m - 1)}{n(n - 1)(n - 2)}.\]
Then we have
\begin{align}
  C_{\textup{III}} & \le c_{1}\frac{m - 1}{n(n - 1)}\lb n + m - 1 - \frac{2(mn - 2m - 2n + 2)}{n - 2} - \lihua{\frac{2m + 2}{n - 2}}\rb\nonumber\\
& = c_{1}\frac{m - 1}{n(n - 1)}\lb n + m - 1 - \lihua{\frac{2mn - 4n - 2m + 6}{n - 2}}\rb\nonumber\\
& = c_{1}\frac{m - 1}{n(n - 1)}\lb n - m + 3 - \lihua{\frac{2m - 2}{n - 2}}\rb\nonumber\\
& \le c_{1}\lb \frac{(m - 1)(n - m + 3)}{n(n - 1)} - \frac{2(m - 1)^{2}}{n(n - 1)(n - 2)}\rb \nonumber\\
& \le c_{1}\lb \frac{(n + 2)^{2}}{4n(n - 1)} - \frac{2(m - 1)^{2}}{n(n - 1)(n - 2)}\rb \nonumber\\
& \le c_{1}\lb \frac{(n + 2)^{2}}{4n(n - 1)} - \frac{18}{n(n - 1)(n - 2)}\rb.

\end{align}
If $n\ge 7$,
\[C_{\textup{III}}\le c_{1}\frac{(n + 2)^{2}}{4n(n - 1)}\le \frac{81c_{1}}{168}\le \frac{c_1}{2}.\]
For $n = 4, 5, 6$, we can also verify that
$C_{\textup{III}}\le \frac{c_1}{2}.$

Therefore, we always have
$
  C_{\textup{III}}\le \frac{c_{1}}{2}.
$

\medskip
Using the above bounds for $(C_{\textup{I}}, C_{\textup{II}}, C_{\textup{III}})$ in \eqref{eq:exact_var}, we obtain that
\begin{align*}
\Var\lb Q_\T \rb&\le c_1\lb 1 - \frac{4(m - 1)(n - m - 1)}{(n - 2)(n - 3)}\rb\sum_{i=1}^{n}Q_{ii}^{2}\\
& + c_1\frac{4(m - 1)(n - m - 1)}{(n - 2)(n - 3)} \frac{(\tr(Q))^{2}}{n} + \frac{c_1}{2}\sum_{[i, j]}(Q_{ij}^{2} + Q_{ij}Q_{ji}).
\end{align*}
Because $(\tr(Q))^{2} \le n\sum_{i=1}^{n}Q_{ii}^{2}$ and $Q_{ij}Q_{ji}\le  (Q_{ij}^{2} + Q_{ji}^{2} ) / 2$, we conclude that
$\Var\lb Q_\T\rb \le c_1 \|Q\|_{F}^{2}.$

\section{Proofs of the Lemmas in Section 6.2}\label{subapp:lemmas}

\begin{proof}[Proof of Lemma \ref{lem:A3}]
  Using the definitions of $ \sigma_{n}^{2}$ and $\rho_{e}$, we have
  \begin{align*}
    \sigma_{n}^{2} &= \lb\frac{1}{n_{1}} - \frac{1}{n}\rb\sum_{i=1}^{n}e_i^2(1) + \lb\frac{1}{n_{0}} - \frac{1}{n}\rb\sum_{i=1}^{n}e_i^2(0) + \frac{2}{n}\sum_{i=1}^{n}e_{i}(1)   e_{i}(0)\\
& = \frac{n_{0}}{n_{1}n}\sum_{i=1}^{n}e_i^2(1) + \frac{n_{1}}{n_{0}n}\sum_{i=1}^{n}e_i^2(0) + \frac{2\rho_{e}}{n}\sqrt{\sum_{i=1}^{n}e_i^2(1)}\sqrt{\sum_{i=1}^{n}e_i^2(0)}.
  \end{align*}
If $\rho_{e} \ge 0$, then
\[\sigma_{n}^{2}\ge \frac{n_{0}}{n_{1}n}\sum_{i=1}^{n}e_i^2(1) + \frac{n_{1}}{n_{0}n}\sum_{i=1}^{n}e_i^2(0)\ge \min\left\{\frac{n_{1}}{n_{0}}, \frac{n_{0}}{n_{1}}\right\}\cure_{2}.\]
If $\rho_{e} < 0$, then using the fact
\[
\lb   \sqrt{\frac{n_0}{n_1}} a - \sqrt{\frac{n_1}{n_0}} b   \rb^2    \geq 0 \Longleftrightarrow
2ab \le \frac{n_{0}}{n_{1}}a^{2} + \frac{n_{1}}{n_{0}}b^{2},
\]
we obtain that
\begin{align*}
  \sigma_{n}^{2} &\ge (1 + \rho_{e})\lb\frac{n_{0}}{n_{1}n}\sum_{i=1}^{n}e_i^2(1) + \frac{n_{1}}{n_{0}n}\sum_{i=1}^{n}e_i^2(0)\rb \ge \eta\min\left\{\frac{n_{1}}{n_{0}}, \frac{n_{0}}{n_{1}}\right\}\cure_{2}.
\end{align*}
Putting the pieces together, we complete the proof.
\end{proof}

\begin{proof}[Proof of Lemma \ref{lem:tauhat}]
Recall that $\hat{\mu}_{t}$ is the intercept from the OLS fit of $Y_{t}^{\mathrm{obs}}$ on $\onet $ and $X_{t}.$ From the Frisch--Waugh Theorem, it is identical to the coefficient of the OLS fit of the residual $ (\id - H_t ) Y_{t}^{\mathrm{obs}}$ on the residual $ (\id - H_t ) \onet $, after projecting onto $X_t$:
\[\hat{\mu}_{t}
= \frac{  \onet \tran  (\id - H_t ) \tran     (\id - H_t ) Y_{t}^{\mathrm{obs}}    }{  \|  (\id - H_t ) \onet  \|_2^2  }
= \frac{\onet \tran  (\id - H_t ) Y_{t}^{\mathrm{obs}} }{\onet \tran  (\id - H_t ) \onet }.\]
Using the definition \eqref{eq:et} and the fact that $ (\id - H_t ) X_{t} = 0$, we have
\begin{eqnarray*}
(\id - H_t ) Y_{t}^{\mathrm{obs}} =  (\id - H_t ) (\mu_{t}\onet  + X_{t}\beta_{t} + e_{t}(t)) = \mu_{t} (\id - H_t ) \onet  +  (\id - H_t ) e_{t}(t),\\
\Longrightarrow
 \hat{\mu}_{t} = \mu_{t} + \frac{\onet \tran  (\id - H_t )   e_t(t)  }{\onet \tran  (\id - H_t ) \onet } = \mu_{t} + \frac{\onet \tran e_{t}(t) / n_{t} - \onet \tran H_t e_{t}(t) / n_{t}}{1 - \onet \tran H_t \onet  / n_{t}}.
\end{eqnarray*}
Recalling that $\tau = \mu_1 - \mu_0$, we complete the proof.
\end{proof}

\begin{proof}[Proof of Lemma \ref{lem:vector_concentration}]
Because $\|  U \|_\op \leq \|  U \|_F$, Lemma \ref{thm:vector_concentration} implies that with probability $1-\delta$,
\[\left\|\sum_{i\in \T}u_{i}\right\|_{2}  \Big /  \|U\|_{F} \le \sqrt{\frac{m(n - m)}{n(n - 1)}} +  \sqrt{8\log \frac{1}{\delta}},
\]
which further implies
$
\left\|\sum_{i\in \T}u_{i}\right\|_{2} = O_\P \lb  \|U\|_{F}  \rb.
$
This immediately implies the three results in Lemma \ref{lem:vector_concentration} by choosing appropriate $U.$

Let $u_{i} = e_{i}(t)$ with $\sum_{i=1}^{n} u_{i} = 0$, $U = (u_{1}, \ldots,u_{n})\tran \in \R^{n\times 1}$, and
$ \|U\|_{F}^{2} = \sum_{i=1}^{n}u_{i}^{2} =\sum_{i=1}^{n} e_{i}^2(t) .$
Therefore,
\[ \onet\tran e_{t}(t)   = \left\|   \sum_{i\in \mathcal{T}_t} u_i \right\|_{2} =O_\P( \|  U \|_F)
= O_\P\lb    \sqrt{\sum_{i=1}^{n} e_{i}^2(t)} \rb = O_\P\lb \sqrt{n\cure_{2}}\rb.\]

Let $u_{i} = x_{i}$ with $\sum_{i=1}^{n} u_{i} = 0$, $U=X$, and $\| U \|_F =  \|X\|_{F}  = \sqrt{\tr(X\tran X)} = \tr(n \id ) = np$. Therefore,
$$
\|X_{t}\tran \onet \|_{2} = \left \|   \sum_{i\in \mathcal{T}_t}  u_i  \right  \|_2   = O_\P\lb \|U\|_{F}\rb = O_\P\lb\sqrt{np}\rb.
$$

Let $u_{i} = x_{i}e_{i}(t)$ with $\sum_{i=1}^{n} u_{i} = 0$ due to \eqref{eq:error_key_property}. Therefore,
\[\|X_{t}\tran e_{t}(t)\|_{2}
= \left \|   \sum_{i\in \mathcal{T}_t}  u_i  \right  \|_2
= O_\P\lb \sqrt{\sum_{i=1}^{n}\|x_{i}\|^{2}e_{i}^2(t)  }\rb.\]
Recalling \eqref{eq:hatmatrix} that
$  \|x_{i}\|^{2}_2 = n H_{ii}\le n\kappa,$
we have
$\|X_{t}\tran e_{t}(t)\|_{2} = O_\P\lb n\sqrt{\cure_{2}\kappa}\rb.$

\end{proof}

We need the following proposition to prove Lemma \ref{lem:matrix_concentration}.
\begin{proposition}\label{prop:matrix_inv_diff}
$A$ and $ B$ are two symmetric matrices. $A$ is positive definite， and $A + B$ is invertible. Then
\[\|(A + B)^{-1} - A^{-1}\|_{\op}\le \frac{\|A^{-1}\|_{\op}^{2}\cdot \|B\|_{\op}}{1 - \min\{1, \|A^{-1}\|_{\op}\cdot \|B\|_{\op}\}}.\]
\end{proposition}

\begin{proof}[Proof of Proposition \ref{prop:matrix_inv_diff}]
Let $M = A^{-\frac{1}{2}}BA^{-\frac{1}{2}}$ and $\Lambda(M)$ be its spectrum. By definition, $\|M\|_{\op}\le \|A^{-1}\|_{\op}\cdot \|B\|_{\op}.$ If $\|A^{-1}\|_{\op}\cdot \|B\|_{\op} \ge 1$, the inequality is trivial because the right-hand side of it is $\infty$. Without loss of generality, we assume $\|A^{-1}\|_{\op}\cdot \|B\|_{\op} < 1$, which implies $\|M\|_{\op} < 1.$

Proposition \ref{prop:matrix_inv_diff} follows by combining
  \begin{align*}
\|(A + B)^{-1} - A^{-1}\|_{\op} & = \|A^{-\frac{1}{2}}((\id  + M)^{-1} - \id )A^{-\frac{1}{2}}\|_{\op}\\
& \le \|A^{-1}\|_{\op}\cdot \|\id  - (\id  + M)^{-1}\|_{\op}
  \end{align*}
and
\[\|\id  - (\id  + M)^{-1}\|_{\op}\le \sup_{\lambda\in \Lambda(M)}\bigg|\frac{\lambda}{1 + \lambda}\bigg| =  \frac{\|M\|_{\op}}{1 - \|M\|_{\op}} \leq \frac{\|A^{-1}\|_{\op} \cdot \|B\|_{\op}}{1 - \|A^{-1}\|_{\op}\cdot \|B\|_{\op}}   .\]

\end{proof}

\begin{proof}[Proof of Lemma \ref{lem:matrix_concentration}]
  Let $V_{i} = x_{i}x_{i}\tran  - \id $, then $\sum_{i=1}^{n}V_{i} = 0$. In the following, we will repeatedly use the basic facts: $n^{-1} X\tran X = \id$, $\| x_i\|_2^2 = nH_{ii}$, and $\sum_{i=1}^{n}x_{i}x_{i}\tran = XX\tran = nH$.
  Recalling the definitions of $\nu, \nu_{+}$ and $\nu_{-}$ in Lemma \ref{thm:matrix_concentration}, we have
\begin{align*}
   \nu^{2} &= \left\|\frac{1}{n}\sum_{i=1}^{n}V_{i}^{2}\right\|_{\op} = \left\|\frac{1}{n}\sum_{i=1}^{n}\lb \|x_{i}\|_{2}^{2} x_{i}x_{i}\tran  - 2x_{i}x_{i}\tran  + \id \rb \right\|_{\op} \\
& = \left\|\lb\frac{1}{n}\sum_{i=1}^{n}\|x_{i}\|_{2}^{2} x_{i}x_{i}\tran \rb - \id \right\|_{\op} = \left\|\lb \sum_{i=1}^{n}H_{ii} x_{i}x_{i}\tran \rb - \id \right\|_{\op}\\
& \le \left\|\sum_{i=1}^{n}H_{ii} x_{i}x_{i}\tran \right\|_{\op} + 1 \le \kappa \left\|\sum_{i=1}^{n}x_{i}x_{i}\tran \right\|_{\op} + 1
= n \kappa  \| H\|_\op + 1 = n \kappa + 1,\\
  \nu_{+} &= \max_{1\leq i \leq n}\|x_{i}x_{i}\tran  - \id \|_{\op} \le \max_{1\leq i \leq n}\|x_{i}\|_{2}^{2} + 1
  =  n\max_{1\leq i \leq n} H_{ii}  +1  = n\kappa + 1,\\
  \nu_{-}^{2} & = \sup_{\omega\in \S^{p-1}}\frac{1}{n}\sum_{i=1}^{n}(\omega\tran V_{i}\omega)^{2} = \sup_{\omega\in \S^{p-1}}\frac{1}{n}\sum_{i=1}^{n}((x_{i}\tran \omega)^{2} - 1)^{2}\\
& = \sup_{\omega\in \S^{p-1}}\frac{1}{n}\sum_{i=1}^{n}\left[(x_{i}\tran \omega)^{4} - 2(x_{i}\tran \omega)^{2} + 1\right]\\
& = \sup_{\omega\in \S^{p-1}}\frac{1}{n}\sum_{i=1}^{n}(x_{i}\tran \omega)^{4} - 2\omega\tran \lb\frac{X\tran X}{n}\rb\omega + 1\\
& = \sup_{\omega\in \S^{p-1}}\frac{1}{n}\sum_{i=1}^{n}(x_{i}\tran \omega)^{4} - 1 \le \sup_{\omega\in \S^{p-1}}\frac{1}{n}\sum_{i=1}^{n}(x_{i}\tran \omega)^{4}\\
& \le \sup_{\omega\in \S^{p-1}}\frac{1}{n}\sum_{i=1}^{n}\|x_{i}\|_{2}^{2}(x_{i}\tran \omega)^{2} = \left\|\sum_{i=1}^{n}H_{ii}x_{i}x_{i}\tran \right\|_{\op} \le n\kappa .
\end{align*}
By Lemma \ref{thm:matrix_concentration},
\begin{align*}
  \left\|\Sigma_t - \id \right \|_{\op}
  = \frac{1}{n_{t}}\left\|\sum_{i \in \mathcal{T}_t} V_{i}\right\|_{\op}
  =  &O_\P\lb \frac{1}{n_{t}}\left[n\sqrt{C(p)\kappa} + nC(p)\kappa + n\sqrt{\kappa}\right]\rb\\
  = & O_\P\lb \sqrt{\kappa\log p} + \kappa\log p\rb.
\end{align*}
By Assumption \ref{assumption::A2}, $\kappa \log p = o(1)$, and therefore the first result holds:
\begin{equation}\label{eq:inv1}
\left\|\Sigma_t - \id \right \|_{\op} = O_\P\lb \sqrt{\kappa\log p} \rb = o_\P(1).
\end{equation}

 Thus with probability $1 - o(1)$,
\begin{equation}\label{eq:original_concentration}
\left\|\Sigma_t - \id \right \|_{\op} \le \frac{1}{2}\Longrightarrow \left\|\Sigma_t\right\|_{\op}\ge \frac{1}{2},
\end{equation}
where we use the convexity of $\|\cdot\|_{\op}$. Note that for any Hermitian matrix $A$, $\|A^{-1}\|_{\op} = \lambda_{\min}(A)^{-1}$ where $\lambda_{\min}$ denotes the minimum eigenvalue. Thus with probability $1 - o(1)$,
\begin{equation}\label{eq:le2}
\left\|\Sigma_t^{-1}\right\|_{\op}\le 2.
\end{equation}
Therefore, the second result holds:
$
\left\|\Sigma_t^{-1}\right\|_{\op} = O_\P(1).
$

To prove the third result, we apply Proposition \ref{prop:matrix_inv_diff} with $A = \id $ and $B = \Sigma_t  -  \id $. By \eqref{eq:original_concentration} and \eqref{eq:le2}, with probability $1 - o(1)$, $A + B$ is invertible and
$\|B\|_{\op} \le 1/2.$
Together with \eqref{eq:inv1}, we have
\[\left\|\Sigma_t^{-1} - \id \right \|_{\op} = O_\P\lb\frac{\|B\|_{\op}}{1 - \|B\|_{\op}}\rb = O_\P(\|B\|_{\op}) = O_\P(\sqrt{\kappa \log p}).\]
\end{proof}

\begin{proof}[Proof of Lemma \ref{lemma:Mt}]
First, \eqref{eq:error_key_property} implies
\[\one\tran Q(t)   = \one\tran H \diag(e(t)) = \one\tran X(X\tran X)^{-1}X\tran \diag(e(t)) = 0,\]
\[Q(t) \one = H \diag(e(t))\one = He(t) = X(X\tran X)^{-1}X\tran e(t) = 0,\]
which further imply $\one\tran Q(t) \one = 0.$
Second, \eqref{eq:hatDeltat} implies
$
\text{tr}(Q(t) ) =n  \Delta_t .
$
Third,
$$
  \|Q(t) \|_{F}^{2} = \sum_{i=1}^{n} \sum_{j=1}^{n} H_{ij}^{2}e_{j}^2(t)    = \sum_{j=1}^{n}e_{j}^2(t)    \lb\sum_{i=1}^{n}H_{ij}^{2}\rb.
$$
Because $H$ is idempotent,
$ H\tran H = H  \Longrightarrow \sum_{i=1}^{n}H_{ij}^{2} = H_{jj}$ for all $j$. Thus,
$
  \|Q(t) \|_{F}^{2} = \sum_{j=1}^{n}e_{j}^2(t)   H_{jj} \le n \cure_{2}\kappa.
$
\end{proof}

\newpage

\begin{center}
  \begin{LARGE}
    Supplementary Material II: Propositions 3.1 and 3.2
  \end{LARGE}
\end{center}

The proofs rely on the following results.

\begin{proposition}\label{prop:yaskov}[modified version of Corollary 3.1 of \cite{yaskov14}]
  Let $Z_{i}$ be i.i.d. random vectors in $\R^{p}$ with mean $0$ and covariance $\id $. Suppose
\[L(\delta)\triangleq   \sup_{ \nu \in \mathcal{S}^{p-1} }
\E |\nu\tran Z_{i}|^{\delta} < \infty\]
for some $\delta > 2$. For any constant $C > 0$, with probability $1 - e^{-Cp}$,
\[\lambda_{\min}\lb \frac{Z\tran Z}{n}\rb\ge 1 - 5\lb\frac{pC}{n}\rb^{\frac{\delta}{\delta + 2}}L(\delta)^{\frac{2}{\delta + 2}}\lb 1+ \frac{1}{C}\rb.\]
\end{proposition}

\begin{proof}[Proof of Proposition \ref{prop:yaskov}]
Write $y = p / n$ and $L = L(\delta)$. The proof of Corollary 3.1 of \citet[][page 6]{yaskov14} showed that for any $a > 0$,
\[\P\lb \lambda_{\min}\lb \frac{Z\tran Z}{n}\rb < 1 - 4La^{- \delta / 2} - 5ay\rb\le \exps{-La^{-1 - \delta / 2}n}.\]
Let $a = \lb Cy / L\rb^{-2/ (\delta + 2)}$. Then the right-hand side is $1 - e^{-Cp}$. Thus with probability $1 - e^{-Cp}$,
\begin{align*}
 \lambda_{\min}\lb \frac{Z\tran Z}{n}\rb
 &\ge 1 - y^{\frac{\delta}{\delta + 2}}L^{\frac{2}{\delta + 2}}\lb 5 C^{-\frac{2}{\delta + 2}} + 4C^{\frac{\delta}{\delta + 2}}\rb
\\
& \ge 1 - 5\lb Cy\rb^{\frac{\delta}{\delta + 2}}L^{\frac{2}{\delta + 2}}\lb 1+ \frac{1}{C}\rb.
\end{align*}
\end{proof}

\begin{proposition}[Theorem 1 of \cite{tikhomirov17}]\label{prop:operator_norm}
  Let $Z_{i}$ be i.i.d. random vectors in $\R^{p}$ with mean $0$ and covariance $\id$. Suppose
\[L(\delta)\triangleq \sup_{ \nu \in \mathcal{S}^{p-1} }
\E |\nu\tran Z_{i}|^{\delta} < \infty\]
for some $\delta > 2$. Then with probability at least $1 - 1/n$,
\[\nu(\delta)^{-1}\left\|\frac{Z\tran Z}{n} - \id \right \|_{\op} \le \frac{\max_{1\leq i \leq n}\|Z_{i}\|_{2}^{2}}{n} + L(\delta)^{\frac{2}{\delta}}\left\{\lb \frac{p}{n}\rb^{\frac{\delta - 2}{\delta}}\log^{4}\lb \frac{n}{p}\rb + \lb\frac{p}{n}\rb^{\frac{\min\{\delta - 2, 2\}}{\min\{\delta, 4\}}}\right\},\]
for some constant $\nu(\delta)$ depending only on $\delta$.
\end{proposition}

\begin{proposition}[Theorem 2 of \cite{vonbahr65}]\label{prop:vonbahresseen}
  Let $Z_{i}$ be independent mean-zero random variables. Then for any $r \in [1, 2)$,
\[\E \bigg|\sum_{i=1}^{n}Z_{i}\bigg|^{r}\le 2\sum_{i=1}^{n}\E |Z_{i}|^{r}.\]
\end{proposition}

\section{Proof of Proposition ~$\lowercase{\ref{prop:kappa_iid}}$}

\subsection{A lemma}

First we prove a more general result.
\begin{lemma}\label{prop:kappa}
Let $Z_{i}$ be i.i.d. random vectors in $\R^{p}$ with mean $\mu \in \R^p $ and covariance matrix $\Sigma\in \R^{p\times p}$. Let $\td{Z}_{i} = \Sigma^{-1/2}(Z_{i} - \mu)$, and assume
\[\sup_{\nu\in \S^{p-1}}\E |\nu\tran \td{Z}_{i}|^{\delta} = O(1), \quad \text{and}\quad
\max_{1\leq i \leq n}\big|\|\td{Z}_{i}\|_{2}^{2} - \E \|\td{Z}_{i}\|_{2}^{2}\big| = O_\P(\omega(n, p)),\]
for some $\delta> 2$ and some function $\omega(n, p)$ increasing in $n$ and $p$. Further let $Z = (Z_{1}\tran , \ldots, Z_{n}\tran )\tran $ and $X = \projone  Z$ so that $X$ has centered columns. If $p = O(n^{\gamma})$ for some $\gamma < 1$, then over the randomness of $Z$,
\[\kappa = \frac{p}{n} + O_\P\lb \frac{\omega(n, p)}{n} + \lb\frac{p}{n}\rb^{\frac{2\delta - 2}{\delta}}\log^{4}\lb \frac{n}{p}\rb + \lb\frac{p}{n}\rb^{\frac{\min\{2\delta - 2, 6\}}{\min\{\delta, 4\}}}\rb.\]
\end{lemma}

\begin{proof}[Proof of Lemma \ref{prop:kappa}]
Let $\td{Z} = (\td{Z}_{1}\tran , \ldots, \td{Z}_{n})\tran $ and
$\td{X} = \projone  \td{Z}.$
Then
$\td{X} = \projone  \lb Z - \one  \mu\tran \rb\Sigma^{-\frac{1}{2}} = \projone  Z\Sigma^{-\frac{1}{2}},$
and thus
\begin{align*}
  \td{X}(\td{X}\tran \td{X})^{-1}\td{X}\tran  &= \projone  Z\lb Z\tran \projone  Z\rb^{-1}Z\tran \projone
 = X(X\tran X)^{-1}X\tran .
\end{align*}
Therefore, we can assume $\mu = 0$ and $ \Sigma = \id $ without loss of generality, in which case $Z_{i} = \td{Z}_{i}$ has mean 0 and covariance matrix $\id$.

By definition,
$H_{ii} = x_{i}\tran (X\tran X)^{-1}x_{i}$, and therefore
\begin{align}
  H_{ii} & = \frac{1}{n}x_{i}\tran \lb \lb\frac{X\tran X}{n}\rb^{-1} - \id \rb  x_{i} + \frac{\|x_{i}\|_{2}^{2}}{n}\nonumber\\
& \le \frac{\|x_{i}\|_{2}^{2}}{n}\lb 1 + \left\|\lb\frac{X\tran X}{n}\rb^{-1} - \id \right \|_{\op}\rb.\label{eq:kappaHii1}
\end{align}
To bound $\kappa$, we need to bound two key terms below.

\paragraph{Bounding $\left\|\lb  n^{-1} X\tran X \rb^{-1} - \id \right \|_{\op}$} Let $\bar{Z} = n^{-1}  \sum_{i=1}^{n}Z_{i}$. Note that
\[\E \|\bar{Z}\|_{2}^{2} = \frac{1}{n^{2}}\sum_{i=1}^{n}\E \|Z_{i}\|_{2}^{2} = \frac{1}{n}\E \|Z_{1}\|_{2}^{2} = \frac{p}{n}.\]
By Markov's inequality,
\begin{equation}\label{eq:barZ}
  \|\bar{Z}\|_{2}^{2} = O_{\P}\lb\frac{p}{n}\rb,
\end{equation}
and more precisely,
\begin{equation}\label{eq:E1}
  \P\lb \|\bar{Z}\|_{2}^{2} \le \sqrt{\frac{p}{n}}\rb = 1 - \P\lb \|\bar{Z}\|_{2}^{2} > \sqrt{\frac{p}{n}}\rb \ge 1 - \sqrt{\frac{p}{n}}.
\end{equation}
Let $\mathcal{A}_{1}$ denote the above event that $\|\bar{Z}\|_{2}^{2} \le \sqrt{ p / n}$. Then
$
\P(\mathcal{A}_{1}) \ge 1 - \sqrt{ p / n }.
$

By Proposition \ref{prop:operator_norm},
\[\left\|\frac{Z\tran Z}{n} - \id \right \|_{\op}= O_\P\lb \frac{\max_{1\leq i \leq n}\|Z_{i}\|_{2}^{2}}{n} + \lb\frac{p}{n}\rb^{\frac{\delta - 2}{\delta}}\log^{4}\lb \frac{n}{p}\rb + \lb\frac{p}{n}\rb^{\frac{\min\{\delta - 2, 2\}}{\min\{\delta, 4\}}}\rb  .\]
By the condition of Lemma \ref{prop:kappa},
\begin{equation}\label{eq:max2norm}
\frac{\max_{1\leq i \leq n}\|Z_{i}\|_{2}^{2}}{n} = \frac{p}{n}  +    \max_{1\leq i \leq n}\frac{\|Z_{i}\|_{2}^{2}    - \E \|Z_{i}\|_{2}^{2}}{n}  = \frac{p}{n} + O_\P\lb \frac{\omega(n, p)}{n}\rb.
\end{equation}
Combining the above three equations, we have
\begin{align}
  & \left\|\frac{X\tran X}{n} - \id \right \|_{\op} = \left\|\frac{Z\tran Z}{n} - \id - \bar{Z}\bar{Z}\tran \right\|_{\op}\le \left\|\frac{Z\tran Z}{n} - \id \right \|_{\op} + \|\bar{Z}\|_{2}^{2}\nonumber\\
= & O_\P\lb \frac{p}{n} + \frac{\omega(n, p)}{n} + \lb\frac{p}{n}\rb^{\frac{\delta - 2}{\delta}}\log^{4}\lb \frac{n}{p}\rb + \lb\frac{p}{n}\rb^{\frac{\min\{\delta - 2, 2\}}{\min\{\delta, 4\}}}\rb\nonumber\\
= & O_\P\lb\frac{\omega(n, p)}{n} + \lb\frac{p}{n}\rb^{\frac{\delta - 2}{\delta}}\log^{4}\lb \frac{n}{p}\rb + \lb\frac{p}{n}\rb^{\frac{\min\{\delta - 2, 2\}}{\min\{\delta, 4\}}}\rb,\label{eq:XTX-I}
\end{align}
where the last line uses the fact that the third term dominates the first term due to $p / n \rightarrow 0$. On the other hand, by Proposition \ref{prop:yaskov} with $C = \sqrt{n / p}$, with probability $1 - e^{-\sqrt{np}}$,
\begin{align}
  \lambda_{\min}\lb\frac{Z\tran Z}{n}\rb&\ge 1 - 5\lb\sqrt{\frac{p}{n}}\rb^{\frac{\delta}{\delta + 2}}L(\delta)^{\frac{2}{\delta + 2}}\lb 1 + \sqrt{\frac{p}{n}}\rb \nonumber\\
&\ge 1 - 10\lb\frac{p}{n}\rb^{\frac{\delta}{2(\delta + 2)}}L(\delta)^{\frac{2}{\delta + 2}}. \label{eq::A2}
\end{align}
Let $\mathcal{A}_{2}$ denote the event in \eqref{eq::A2}. Then
$
\P(\mathcal{A}_{2})\ge 1 - e^{-\sqrt{np}}.
$

Note that for any Hermitian matrices $A$ and $B$, the convexity of $\|\cdot\|_{\op}$ implies that
\[|\lambda_{\min}(A) - \lambda_{\min}(B)| = |\lambda_{\max}(-A) - \lambda_{\max}(-B)| \le \|-A - (-B)\|_{\op} = \|A - B\|_{\op}.\]
 We have
\[\lambda_{\min}\lb\frac{X\tran X}{n}\rb\ge \lambda_{\min}\lb\frac{Z\tran Z}{n}\rb - \|\bar{Z}\bar{Z}\tran \|_{\op} = \lambda_{\min}\lb\frac{Z\tran Z}{n}\rb - \|\bar{Z}\|_{2}^{2}.\]
Let $\mathcal{A} = \mathcal{A}_{1}\cup \mathcal{A}_{2}$. Then on $\mathcal{A}$,
\begin{equation*}
  \lambda_{\min}\lb\frac{X\tran X}{n}\rb \ge 1 - 10\lb\frac{p}{n}\rb^{\frac{\delta}{2(\delta + 2)}}L(\delta)^{\frac{2}{\delta + 2}} - \sqrt{\frac{p}{n}}.
\end{equation*}
Since $p / n\rightarrow 0$, for sufficiently large $n$,
\[\lambda_{\min}\lb\frac{X\tran X}{n}\rb \ge \frac{1}{2}\]
with probability
\[\P(\mathcal{A})\ge \P(\mathcal{A}_{1}) + \P(\mathcal{A}_{2}) - 1\ge 1 - e^{-\sqrt{np}L(\delta)} - \sqrt{\frac{p}{n}} = 1 - o(1).\]

Finally, using Proposition \ref{prop:matrix_inv_diff} with $A=\id$ and $B = n^{-1} X\tran X - \id$, by Slusky's lemma, we have that
\begin{equation}
  \label{eq:XTXinv-I}
  \left\|\lb\frac{X\tran X}{n}\rb^{-1} - \id \right \|_{\op} = O_\P\lb\frac{\omega(n, p)}{n} + \lb\frac{p}{n}\rb^{\frac{\delta - 2}{\delta}}\log^{4}\lb \frac{n}{p}\rb + \lb\frac{p}{n}\rb^{\frac{\min\{\delta - 2, 2\}}{\min\{\delta, 4\}}}\rb.
\end{equation}
Because  $p = O(n^{\gamma})$ for some $\gamma < 1$,
\begin{equation}
  \label{eq:XTXinv-I2}
  \left\|\lb\frac{X\tran X}{n}\rb^{-1} - \id \right \|_{\op} = O_\P\lb\frac{\omega(n, p)}{n}\rb + o_\P(1).
\end{equation}

\paragraph{Bounding $  \max_{1\leq i\leq n} \|x_{i}\|_{2}^{2} $}
Because  $x_{i} = Z_{i} - \bar{Z}$, the Cauchy--Schwarz inequality implies
\[\|x_{i}\|_{2}^{2} = \|Z_{i}\|_{2}^{2} - 2Z_{i}\tran \bar{Z} + \|\bar{Z}\|_{2}^{2}\le \|Z_{i}\|_{2}^{2} + 2\|Z_{i}\|_{2}\|\bar{Z}\|_{2} + \|\bar{Z}\|_{2}^{2}.\]
By \eqref{eq:max2norm} and \eqref{eq:barZ},
\begin{align}
  \frac{\max_{1\leq i \leq n}\|x_{i}\|_{2}^{2}}{n}& = \lihua{\frac{\E \|Z_{i}\|_{2}^{2}}{n} + \frac{\max_{i}\|Z_{i}\|_{2}^{2} - \E \|Z_{i}\|_{2}^{2}+ 2 \|Z_{i}\|_{2}\|\bar{Z}\|_{2} + \|\bar{Z}\|_{2}^{2}}{n}} \nonumber\\
& = \frac{p}{n} + O_\P\lb \frac{\omega(n, p)}{n} + \sqrt{\frac{(p + \omega(n, p))p}{n^{3}}} + \frac{p}{n^{2}}\rb\nonumber\\
& = \frac{p}{n} + O_\P\lb \frac{\omega(n, p)}{n} + \sqrt{\frac{\omega(n, p)p}{n^{3}}} + \frac{p}{n^{3/2}} \rb\nonumber.
\end{align}
Because  $\omega(n, p)$ is increasing and $p / n\rightarrow 0$, we have
\[\sqrt{\frac{\omega(n, p)p}{n^{3}}} = O\lb\frac{\omega(n, p)}{n}\lb\frac{p}{n}\rb^{1/2}\rb = o\lb\frac{\omega(n, p)}{n}\rb.\]
Thus, we obtain that
\begin{equation}
  \label{eq:max2normx}
  \frac{\max_{1\leq i \leq n}\|x_{i}\|_{2}^{2}}{n} = \frac{p}{n} + O_\P\lb \frac{\omega(n, p)}{n}  + \frac{p}{n^{3/2}} \rb.
\end{equation}

Putting \eqref{eq:kappaHii1}, \eqref{eq:XTXinv-I2} and \eqref{eq:max2normx} together and using some tedious cancellations, we have
\begin{align}
&\kappa  = \frac{p}{n} + O_\P\lb \frac{\omega(n, p)}{n}  + \frac{p}{n^{3/2}}\rb \nonumber\\
& \quad + O_\P\lb\frac{\omega^{2}(n, p)}{n^{2}} + \lb\frac{p}{n}\rb^{1 + \frac{\delta - 2}{\delta}}\log^{4}\lb \frac{n}{p}\rb + \lb\frac{p}{n}\rb^{1 + \frac{\min\{\delta - 2, 2\}}{\min\{\delta, 4\}}}\rb.\label{eq:Hii1}
\end{align}
Because
\[
\lb\frac{p}{n}\rb^{1 + \frac{\min\{\delta - 2, 2\}}{\min\{\delta, 4\}}}\ge \lb\frac{p}{n}\rb^{3/2}\ge \frac{p}{n^{3/2}},
\]
\eqref{eq:Hii1} further simplifies to
\begin{align*}
  \kappa &= \frac{p}{n} + O_\P\lb \frac{\omega(n, p)}{n} + \frac{\omega^{2}(n, p)}{n^{2}} + \lb\frac{p}{n}\rb^{\frac{2\delta - 2}{\delta}}\log^{4}\lb \frac{n}{p}\rb + \lb\frac{p}{n}\rb^{\frac{\min\{2\delta - 2, 6\}}{\min\{\delta, 4\}}}\rb\\
& = \frac{p}{n} + O_\P\lb \frac{\omega(n, p)}{n}\max\left\{ \frac{\omega(n, p)}{n}, 1\right\} + \lb\frac{p}{n}\rb^{\frac{2\delta - 2}{\delta}}\log^{4}\lb \frac{n}{p}\rb + \lb\frac{p}{n}\rb^{\frac{\min\{2\delta - 2, 6\}}{\min\{\delta, 4\}}}\rb .
\end{align*}
We complete the proof using $\kappa\le 1$.
\end{proof}

\subsection{Use Lemma \ref{prop:kappa} to prove Proposition \ref{prop:kappa_iid}}

We have argued in the proof of Proposition \ref{prop:kappa} that we can assume $\mu = 0$ without loss of generality. Because the hat matrix is invariant to rescaling, we further assume $\E Z_{ij}^{2} = 1$ without loss of generality. Based on Proposition \ref{prop:kappa}, it suffices to verify
\begin{eqnarray}
  \label{eq:kappa_goal1}
\sup_{\nu\in \S^{p-1}}\E |\nu\tran Z_{i}|^{\delta} &=& O(1),
\\
  \label{eq:kappa_goal2}
\max_{1\leq i \leq n} \big| \|Z_{i}\|_{2}^{2} - \E \|Z_{i}\|_{2}^{2}\big| &=& O_\P\lb n^{\frac{2}{\delta}}p^{\frac{2}{\min\{\delta, 4\}}}\rb.
\end{eqnarray}
If \eqref{eq:kappa_goal1} and \eqref{eq:kappa_goal2} hold, by Proposition \ref{prop:kappa}, we have that
\[\kappa = \frac{p}{n} + O_\P\lb \frac{p^{2/\min\{\delta, 4\}}}{n^{(\delta - 2) / \delta}} + \lb\frac{p}{n}\rb^{\frac{2\delta - 2}{\delta}}\log^{4}\lb \frac{n}{p}\rb + \lb\frac{p}{n}\rb^{\frac{\min\{2\delta - 2, 6\}}{\min\{\delta, 4\}}}\rb.\]
Then we can prove Proposition \ref{prop:kappa_iid} for two cases.

\paragraph{Case 1}
If $\delta > 4$, then $\frac{2\delta - 2}{\delta} < \frac{3}{2} = \frac{\min\{2\delta - 2, 6\}}{\min\{\delta, 4\}}$. Thus the third term dominates the second term in the above $O_\P(\cdot)$, implying
\begin{align*}
  \kappa = \frac{p}{n} + O_\P\lb \frac{p^{1/2}}{n^{(\delta - 2) / \delta}} + \lb\frac{p}{n}\rb^{\frac{3}{2}}\rb.
\end{align*}

\paragraph{Case 2}
If $\delta \le 4$, then
\[\kappa = \frac{p}{n} + O_\P\lb\frac{p^{2 / \delta}}{n^{(\delta - 2) / \delta}} + \lb\frac{p}{n}\rb^{\frac{2\delta - 2}{\delta}} \log^{4}\lb\frac{n}{p}\rb\rb.\]
Because
\[\lb\frac{p}{n}\rb^{\frac{2\delta - 2}{\delta}} = \frac{p^{2 / \delta}}{n^{(\delta - 2) / \delta}}\frac{p^{(2\delta - 4) / \delta}}{n}\le \frac{p^{2 / \delta}}{n^{(\delta - 2) / \delta}}\frac{p}{n} , \]
the first term dominates in the above $O_\P(\cdot)$, implying
\[
\kappa = \frac{p}{n} + O_\P\lb\frac{p^{2 / \delta}}{n^{(\delta - 2) / \delta}}\rb
=  \frac{p}{n} + O_\P\lb\frac{p^{2 / \delta}}{n^{(\delta - 2) / \delta}} + \lb\frac{p}{n}\rb^{3/2}\rb .
\]
The last identity holds because $p^{3/2}/n^{3/2}$ is of smaller order and thus we can add it back.

We will prove \eqref{eq:kappa_goal1} and \eqref{eq:kappa_goal2} below.

\subsubsection{Proving \eqref{eq:kappa_goal1}}
By \citet{rosenthal70}'s inequality,
\begin{align*}
  \E |\nu\tran Z_{i}|^{\delta} &= \E \bigg|\sum_{j=1}^{p}\nu_{j}Z_{ij}\bigg|^{\delta} \le C \lb \sum_{j=1}^{p}|\nu_{j}|^{\delta} \E |Z_{ij}|^{\delta} + \lb \sum_{j=1}^{p}\nu_{j}^{2} \E Z_{ij}^{2}\rb^{\delta / 2}\rb
\end{align*}
where $C$ is a constant depending only on $\delta$. Because  $\|\nu\|_{2} = 1$, we have $\max_{1\leq j \leq p}|\nu_{j}|\le 1$ and thus
\[\sum_{j=1}^{p}|\nu_{j}|^{\delta} \E |Z_{ij}|^{\delta}\le M \sum_{j=1}^{p}|\nu_{j}|^{\delta} \le M \sum_{j=1}^{p}|\nu_{j}|^{\delta} =  M.\]
H\"{o}lder's inequality implies
$\E Z_{ij}^{2} \le \lb \E |Z_{ij}|^{\delta}\rb^{2/\delta}\le M^{2 / \delta},$
which further implies
\[\lb \sum_{j=1}^{p}\nu_{j}^{2} \E Z_{ij}^{2}\rb^{\delta / 2}\le \lb M^{2 / \delta}\rb^{\delta / 2} = M.\]
Because the above two bounds hold regardless of $\nu$, we conclude that
\[ \sup_{\nu\in \S^{p-1}} \E |\nu\tran Z_{i}|^{\delta}\le 2CM = O(1).\]

\subsubsection{Proving \eqref{eq:kappa_goal2}}
Let
$W_{ij} = Z_{ij}^{2} - \E Z_{ij}^{2}.
$
Using Jensen's inequality that $\E |(X + Y) / 2|^{r}\le (\E |X|^{r} + \E |Y|^{r}) / 2$ for any random variables $X, Y$ and $r > 1$, we obtain that
\begin{align*}
  \E |W_{ij}|^{\delta / 2}\le 2^{\delta / 2 - 1}\lb\E |Z_{ij}|^{\delta} + (\E Z_{ij}^{2})^{\delta / 2}\rb
 \le 2^{\delta / 2}\E |Z_{ij}|^{\delta}\le 2^{\delta / 2}M \triangleq \td{M}.
\end{align*}
We consider two cases.

\paragraph{Case 1: $\delta \ge 4$}
By H\"{o}lder's inequality,
$\E W_{ij}^{2}\le \td{M}^{4 / \delta}.$
By \citet{rosenthal70}'s inequality,
\begin{align*}
  \E |\|Z_{i}\|_{2}^{2} - \E \|Z_{i}\|_{2}^{2}|^{\delta / 2} &= \E \bigg|\sum_{j=1}^{p}W_{ij}\bigg|^{\delta / 2}
\le C\lb \sum_{j=1}^{p}\E |W_{ij}|^{\delta / 2} + \lb \sum_{j=1}^{p}\E W_{ij}^{2}\rb^{\delta / 4}\rb\\
&\le  C\lb p\td{M} + p^{\delta / 4}\td{M}\rb \le C\td{M}p^{\delta / 4},
\end{align*}
which implies
$
\E \big| \|Z_{i}\|_{2}^{2} - \E \|Z_{i}\|_{2}^{2}\big|^{\delta / 2} = O\lb p^{\delta / 4}\rb.
$
As a result,
\begin{align*}
  &\E   \left\{   \max_{1\leq i \leq n}\big| \|Z_{i}\|_{2}^{2} - \E \|Z_{i}\|_{2}^{2}\big|^{\delta / 2} \right\}
   \le \sum_{i=1}^{n}\E \big| \|Z_{i}\|_{2}^{2} - \E \|Z_{i}\|_{2}^{2}\big|^{\delta / 2} = O\lb n p^{\delta / 4}\rb.
\end{align*}
By Markov's inequality,
$
\max_{1\leq i \leq n}\big| \|Z_{i}\|_{2}^{2} - \E \|Z_{i}\|_{2}^{2}\big| = O_\P\lb n^{2/\delta}p^{1/2}\rb.
$

\paragraph{Case 2: $\delta < 4$}
By Proposition \ref{prop:vonbahresseen}, with $\delta / 2\in (1, 2)$,
\[\E |\|Z_{i}\|_{2}^{2} - \E \|Z_{i}\|_{2}^{2}|^{\delta / 2} = \E \bigg|\sum_{j=1}^{p}W_{ij}\bigg|^{\delta / 2}\le 2\sum_{j=1}^{p}\E |W_{ij}|^{\delta / 2}\le 2p\td{M}.\]
Similar to Case 1,
$
\max_{1\le i\le n}\big| \|Z_{i}\|_{2}^{2} - \E \|Z_{i}\|_{2}^{2}\big| = O_\P\lb n^{2 / \delta}p^{2 / \delta}\rb.
$

\section{Proof of Proposition ~$\lowercase{\ref{prop:eq}}$}\label{subapp:cure}

Let $\bar{Y}(t) = n^{-1} \sum_{i=1}^{n}Y_{i}(t)$. Note that $H\one = X(X\tran X)^{-1}X\tran \one = 0$. By definition,
$
e(t) = (\id  - H)\{  Y(t) - \bar{Y}(t)\one \}  = (\id  - H)\{  Y(t) - \E Y_{i}(t)\one \} .
$
Throughout the rest of the proof, we assume that $\E Y_{i}(t) = 0$ without loss of generality, and define
$M(\delta)\triangleq \max_{t=0,1}\E |Y_{i}(t)|^{\delta}.$

\subsection{Bounding $\cure_2$}
Let $Z_{i} = Y_{i}(t)^{2}$. Then the moment condition reads
$\E |Z_{i}|^{\delta / 2} < \infty.$
The Kolmogorov--Marcinkiewicz--Zygmund strong law of large number \citep[][Theorem 4.23]{kallenberg06} implies
\begin{eqnarray*}
\frac{1}{n}\sum_{i=1}^{n}Z_{i}\stackrel{\text{a.s.}}{\rightarrow} \E Z_{1} = O_\P(1), \quad \text{ if } \delta \ge 2, \\
\frac{1}{n^{2/\delta}}\sum_{i=1}^{n}Z_{i} = o(1) \Longrightarrow \frac{1}{n}\sum_{i=1}^{n}Z_{i} = o_\P(n^{2/\delta - 1}),\quad
\text{ if }  \delta < 2.
\end{eqnarray*}
On the other hand,
\[\frac{1}{n}\|e(t)\|_{2}^{2} = \frac{1}{n}Y(t)\tran (\id  - H)Y(t)\le \frac{1}{n}\|Y(t)\|_{2}^{2} = \frac{1}{n}\sum_{i=1}^{n}Z_{i},\]
which further implies the bound for $\cure_{2}$.

\subsection{Bounding $\cure_2^{-1}$}
Without loss of generality, we assume that $Y_{i}(1)$ is not a constant with probability $1$. First we show that
\[
\frac{Y(1)\tran H Y(1)}{Y(1)\tran Y(1)} = o_{\P}(1).
\]
For any permutation $\pi$ on $\{1, \ldots, n\}$, let $H(\pi)$ denote the matrix with
\[H(\pi)_{ij} = H_{\pi(i), \pi(j)}.\]
Because the $Y_{i}(1)$'s are i.i.d., for any $\pi$,
\[(Y_{1}(1), \ldots, Y_{n}(1))\stackrel{\text{d}}{=}(Y_{\pi^{-1}(1)}(1), \ldots, Y_{\pi^{-1}(n)}(1)),\]
and thus
\begin{align*}
 & \frac{Y(1)\tran H(\pi) Y(1)}{Y(1)\tran Y(1)} = \frac{\sum_{i = 1}^{n}\sum_{j=1}^nH_{\pi(i), \pi(j)}Y_{i}(1) Y_{j}(1)}{\sum_{i=1}^{n}Y_{i}(1)^{2}}\\
 = & \frac{\sum_{i = 1}^{n}\sum_{j=1}^nH_{i, j}Y_{\pi^{-1}(i)}(1) Y_{\pi^{-1}(j)}(1)}{\sum_{i=1}^{n}Y_{\pi^{-1}(i)}(1)^{2}}\stackrel{\text{d}}{=} \frac{Y(1)\tran H Y(1)}{Y(1)\tran Y(1)}.
\end{align*}
Furthermore,
\[\frac{Y(1)\tran H Y(1)}{Y(1)\tran Y(1)}\le 1\]
and thus it has finite expectation. This implies that
\[
\E \frac{Y(1)\tran H Y(1)}{Y(1)\tran Y(1)} = \frac{1}{n!}\sum_{\pi} \frac{Y(1)\tran H(\pi) Y(1)}{Y(1)\tran Y(1)} = \frac{1}{n!}\frac{Y(1)\tran H^{*} Y(1)}{Y(1)\tran Y(1)},
\]
where $H^{*} = \sum_{\pi}H(\pi) / n!$ with the summation over all possible permutations. We can show that
\[H_{ii}^{*} = \frac{1}{n}\sum_{i=1}^{n}H_{ii} = \frac{p}{n}, \quad H_{ij}^{*} = \frac{1}{n(n - 1)}\sum_{i\not = j}H_{ij} = -\frac{1}{n(n - 1)}\sum_{i=1}^{n}H_{ii} = - \frac{p}{n (n - 1)},\]
where the last equality uses the fact that $\sum_{i = 1}^{n}\sum_{j=1}^nH_{ij} = 0$. Therefore,
\begin{align*}
  &\E \frac{Y(1)\tran H Y(1)}{Y(1)\tran Y(1)} = \E \frac{Y(1)\tran H^{*} Y(1)}{Y(1)\tran Y(1)}\\
= &\E\frac{\frac{p}{n}Y(1)\tran Y(1) - \frac{p}{n(n - 1)}\sum_{i\not = j}Y_{i}(1)Y_{j}(1)}{Y(1)\tran Y(1)}\\
= & \frac{p}{n-1} - \frac{p}{n(n - 1)}\E\frac{(\sum_{i=1}^{n}Y_{i}(1))^{2}}{Y(1)\tran Y(1)}\le \frac{p}{n - 1}.
\end{align*}
By Markov's inequality, with probability $1 - \frac{2p}{n - 1} = 1 - o(1)$,
\[\frac{Y(1)\tran H Y(1)}{Y(1)\tran Y(1)} \le \frac{1}{2}.\]
Let $\mathcal{A}$ denote this event. Then
\[ \P(\mathcal{A}^{c}) = o(1),\]
and on $\mathcal{A}$,
\[\frac{1}{n}\|e(1)\|_{2}^{2} = \frac{1}{n}Y(1)\tran (I - H)Y(1)\ge \frac{1}{2n}\|Y(1)\|_{2}^{2}.\]

~\\
\noindent  On the other hand, fix $k > 0$, and let
$\td{Z}_{i} = Y_{i}(1)I(|Y_{i}(1)|\le k).$
For sufficiently large $k$,
$\E \td{Z}_{i} > 0.$
By the law of large numbers,
$
 n^{-1}\sum_{i=1}^{n}\td{Z}_{i} = \E \td{Z}_{i} \times  (1 + o_\P(1)).
$
Thus on $\mathcal{A}$,
\[
\cure_{2}\ge \frac{1}{2n}\sum_{i=1}^{n}Y_{i}(1)^{2} \ge \frac{1}{2n}\sum_{i=1}^{n}\td{Z}_{i} = \E \td{Z}_{i} \times (1 + o_\P(1))

\]
Since $\P(\mathcal{A}^{c}) = o(1)$, we conclude that
$
\cure_{2}^{-1} = O_\P(1).
$

\subsection{Bounding $\cure_{\infty}$}
We apply the triangle inequality to obtain
\[\|e(t)\|_{\infty}\le \|Y(t)\|_{\infty} + \|HY(t)\|_{\infty}.\]
We bound the first term using a standard technique and Markov's inequality:
\begin{equation}\label{eq:Ytinfty1}
\E \|Y(t)\|_{\infty}^{\delta}\le \sum_{i=1}^{n}\E |Y_{i}(t)|^{\delta} = n M(\delta)
\Longrightarrow
\|Y(t)\|_{\infty} = O_\P\lb n^{1/\delta}\rb.
\end{equation}
Next we bound the second term $\|HY(t)\|_{\infty}$. Define $\td{Y}(t) = HY(t)$ with
\[\td{Y}_{i}(t) = \sum_{j=1}^{n}H_{ij}Y_{j}(t), \quad (i=1,\ldots, n).\]
Fix $\eps > 0$ and define
\[D = \lb\frac{M(\delta)}{\eps}\rb^{1/\delta}.\]
 We decompose $\td{Y}_{i}(t)$ into two parts:
  \begin{align*}
  \td{Y}_{i}(t) &= \sum_{j=1}^{n}H_{ij}Y_{j}(t)I(|Y_{j}(t)|\le Dn^{1/\delta}) + \sum_{j=1}^{n}H_{ij}Y_{j}(t)I(|Y_{j}(t)|> Dn^{1/\delta})\\
&\triangleq R_{1, i}(t) + R_{2, i}(t).
  \end{align*}
The second term $R_{2, i}(t)$ satisfies
\begin{align}
  \P\lb \exists i, R_{2, i}(t) \not = 0\rb &\le \P\lb \exists j, \,\, |Y_{j}(t)| > Dn^{1/\delta}\rb\le \sum_{j=1}^{n} \P\lb |Y_{j}(t)| > Dn^{1/\delta}\rb\nonumber\\
& \le \sum_{j=1}^{n}\frac{1}{D^{\delta} n}\E |Y_{j}(t)|^{\delta} \le \frac{M(\delta)}{D^{\delta}} = \eps.\label{eq:R2i}
\end{align}
To deal with the first term $R_{1, i}(t)$, we define
\[w_j  (t) = Y_j (t)I(|Y_j (t)|\le Dn^{1/\delta}) - \E \left\{  Y_j (t)I(|Y_j (t)|\le Dn^{1/\delta}) \right\}  ,\]
with $\E w_j (t) = 0$. Because
\[\one\tran H = 0 \Longrightarrow
\sum_{j=1}^{n}H_{ij} = 0\Longrightarrow
\sum_{j=1}^{n}H_{ij} \E \left\{  Y_j (t)I(|Y_j (t)|\le Dn^{1/\delta}) \right\}  = 0.\]
we can rewrite $R_{1,i}(t)$ as
\[R_{1,i}(t) = \sum_{j=1}^{n}H_{ij}w_{j}(t).\]
The rest of the proof proceeds based on two cases.

\paragraph{Case 1: $\delta < 2$}
First, the $w_{j}(t)$'s are i.i.d. with second moment bounded by
  \begin{align*}
    \E w_j (t)^{2} &\le \E \left\{   Y_j ^{2}(t)I(|Y_j (t)|\le Dn^{1/\delta}) \right\} \\
&\le (Dn^{1/\delta})^{2 - \delta}\E |Y_j (t)|^{\delta}\\
& \le n^{(2 - \delta) / \delta} D^{2-\delta} M(\delta)
  = n^{(2 - \delta) / \delta}\eps^{-(2 - \delta) / \delta} M(\delta)^{2/\delta}.
  \end{align*}
Second, using the fact that $\sum_{j=1}^{n}H_{ij}^{2} = H_{ii}$, we obtain
\[\E R_{1, i}(t)^{2} = \sum_{j=1}^{n}H_{ij}^{2}\E w_{j}(t)^{2} = \lihua{\E w_{1}(t)^{2}\lb\sum_{j=1}^{n}H_{ij}^{2}\rb = }H_{ii}\E w_{1}(t)^{2}.\]
Let $R_{1}(t)$ denote the vector $(R_{1, i}(t))_{i=1}^{n}$. Then
\begin{align*}
  \E \|R_{1}(t)\|_{\infty}^{2} &\le \sum_{i=1}^{n} \E R_{1, i}(t)^{2}  = \lb\sum_{i=1}^{n}H_{ii}\rb \E w_{1}(t)^{2}
 \le  pn^{(2 - \delta) / \delta}\eps^{-(2 - \delta) / \delta} M(\delta)^{2/\delta} .
\end{align*}
By Markov's inequality, with probability $1 - \eps$,
\begin{align}\label{eq::r1case1}
  \|R_{1}(t)\|_{\infty} \le \lb\frac{\E \|  R_{1}(t)  \|_{\infty}^{2}}{\eps}\rb^{1/2}
 = p^{1/2}n^{(2 - \delta) / 2\delta}\eps^{-(4-\delta) / 2\delta} M(\delta)^{1/\delta}.
\end{align}
Combining \eqref{eq:R2i} and \eqref{eq::r1case1}, we obtain that with probability $1 - 2\eps$,
\[\|HY(t)\|_{\infty}\le p^{1/2}n^{(2 - \delta) / 2\delta}\eps^{-(4-\delta) / 2\delta} M(\delta)^{1/\delta}.\]
Because  this holds for arbitrary $\eps$, we conclude that if $\delta < 2$,
\[\|HY(t)\|_{\infty} = O_\P\lb p^{1/2}n^{1 / \delta - 1 / 2}\rb = o_{\P}(n^{1/\delta}) . \]

\paragraph{Case 2: $\delta \ge 2$}
Using the convexity of the mapping $|\cdot|^{\delta}$, we have
    \begin{align*}
      \E \bigg|\frac{w_{j}(t)}{2}\bigg|^{\delta} & \le \frac{\E \left\{  |Y_{j}(t)|^{\delta} I(|Y_{j}(t)|\le Dn^{1/\delta}) \right\}
      + |\E \left\{  Y_{j}(t)I(|Y_{j}(t)|\le Dn^{1/\delta}) \right\}|^{\delta}}{2}.
    \end{align*}
Applying Jensen's inequality on the second term, we have
$$
  \E |w_{j}(t)|^{\delta}  \le 2^{\delta}\E \left\{   |Y_{j}(t)|^{\delta} I(|Y_{j}(t)|\le Dn^{1/\delta}) \right\}
\le 2^{\delta}\E |Y_{j}(t)|^{\delta} \le 2^{\delta}M(\delta).
$$
By \citet{rosenthal70}'s inequality, there exists a constant $C$ depending only on $\delta$, such that
\begin{align*}
& \E |R_{1, i}(t)|^{\delta} \le C\lb\sum_{j=1}^{n} \E|H_{ij}w_j (t)|^{\delta} + \lb\sum_{j=1}^{n}\E |H_{ij}w_j (t)|^{2}\rb^{\delta / 2}\rb\\
& \le C\lb \lihua{2^{\delta}}M(\delta) \sum_{j=1}^{n} |H_{ij}|^{\delta} + \lb \lihua{2^{2}}M(2)\sum_{j=1}^{n}H_{ij}^{2}\rb^{\delta / 2}\rb\\
& \le  C\lihua{2^{\delta}}\lb M(\delta)H_{ii}^{\delta / 2 - 1}\sum_{j=1}^{n}H_{ij}^{2} + M(2)^{\delta / 2}H_{ii}^{\delta / 2}\rb\\
&  = C2^{\delta}(M(\delta) + M(2)^{\delta / 2})H_{ii}^{\delta/2} \leq   C2^{\delta}(M(\delta) + M(2)^{\delta / 2})H_{ii}.
\end{align*}
where the last two lines use $\sum_{j=1}^{n}H_{ij}^{2} = H_{ii}$, $H_{ij}^{2}\le H_{ii}$, and $ H_{ii}^{\delta/2} \leq H_{ii}$ due to $H_{ii}\leq 1$ and $\delta/2>1$. As a result,
\begin{align*}
  \E \|R_{1}(t)\|_{\infty}^{\delta} & \le \sum_{i=1}^{n}\E |R_{1, i}(t)|^{\delta}  \le C\lihua{2^{\delta}}(M(\delta) + M(2)^{\delta / 2})\sum_{i=1}^{n}H_{ii}\\
& = C\lihua{2^{\delta}}(M(\delta) + M(2)^{\delta / 2})p.
\end{align*}
Markov's inequality implies that with probability $1-\eps$,
\begin{align}\label{eq::r1case2}
  \|R_{1}(t)\|_{\infty} \le \lb\frac{\E \|      R_1(t)    \|_{\infty}^{\delta}}{\eps}\rb^{1/\delta}
 = p^{1/\delta} \lb C2^{\delta}(M(\delta) + M(2)^{\delta / 2})\rb^{1/\delta}.
\end{align}
Combining \eqref{eq:R2i} and \eqref{eq::r1case2}, we obtain that with probability $1 - 2\eps$,
\[\|HY(t)\|_{\infty}\le p^{1/\delta} \lb C2^{\delta}(M(\delta) + M(2)^{\delta / 2})\rb^{1/\delta} . \]
Because  this holds for arbitrary $\eps$, we conclude that if $\delta \ge 2$,
\[\|HY(t)\|_{\infty} = O_\P\lb p^{1/\delta}\rb = o_{\P}(n^{1/\delta}).\]

\newpage
\begin{center}
  \begin{LARGE}
    Supplementary Material III: Experiments
  \end{LARGE}
\end{center}

\section{Additional Experiments}\label{app:experiments}

Using the following proposition, we know that the solution of $\eps$ in Section \ref{sec::dgp} is the rescaled OLS
residual vector obtained by regressing the leverage scores $(H_{ii})_{i=1}^n$ on $X$ with an intercept.

\begin{proposition}\label{prop:lpqp}
Let $a\in \R^{n}$ be any vector, and $A\in \R^{n\times m}$ be any matrix with $H_{A} = A(A\tran A)^{-1} A\tran $ being its projection matrix. Define $e = (\id - H_{A})a$. Then $x^{*} = n^{1/2}e / \|e\|_{2} $ is the optimal solution of
\[
\max_{x\in \R^{n}}|a\tran x|\quad \mathrm{s.t.}\,\, \|x\|_{2}^{2} / n = 1, A\tran x = 0  .
\]
\end{proposition}

\begin{proof}[Proof of Proposition \ref{prop:lpqp}]
  The constraint $A\tran x = 0$ implies $H_{A}x = 0$. Thus,
$
|a\tran x| =   |a\tran x - a\tran H_{A}x|  = |  a\tran (\id - H_A) x |
= |e\tran x|.
$
The Cauchy--Schwarz inequality implies
$ |e\tran x|\le \|e\|_{2}\|x\|_{2} = n^{1/2}\|e\|_{2}$, with the maximum objective value achieved by $x = n^{1/2} e / \|e\|_{2}$.
\end{proof}

We present more simulation results in the rest of this section.

\subsection{Other experimental results on synthetic datasets}\label{subapp:experiments2}

Section \ref{sec:experiment} shows the results for $X$ contains i.i.d. $t(2)$ entries. Here we plot the results for $X$ containing i.i.d. entries from $N(0, 1)$ and $t(1)$, analogous to the results in Sections \ref{subsec:simulation_results}--\ref{subsec:simulation_regularization}.

The case with $N(0, 1)$ entries exhibits almost the same qualitative pattern; see Fig. \ref{fig::simulation_normal_01}  and Fig. \ref{fig::simulation_normal_02}. However, for the case with $t(1)$ entries, the bias reduction is less effective and none of the variance estimates, including HC3 estimate, is able to protect against undercoverage when $p > n^{1/2}$; see Fig. \ref{fig::simulation_t1_01} and Fig. \ref{fig::simulation_t1_02}.

\begin{figure}[htp]
	\centering
	\begin{subfigure}{0.99\textwidth}
          \centering
		\includegraphics[width=0.48\textwidth]{synthetic_simul_normal_pi2_bias.pdf}
		\includegraphics[width=0.48\textwidth]{synthetic_simul_normal_pi5_bias.pdf}
		\caption{ Relative bias of $\hdtau$ and $\htau$.}
	\end{subfigure}
	\begin{subfigure}{0.99\textwidth}
          \centering
		\includegraphics[width=0.48\textwidth]{synthetic_simul_normal_pi2_sdinflate.pdf}
		\includegraphics[width=0.48\textwidth]{synthetic_simul_normal_pi5_sdinflate.pdf}
		\caption{Ratio of standard deviation between five standard deviation estimates, $\sigma_{n}, \hat{\sigma}_{\textup{HC}0}, \hat{\sigma}_{\textup{HC}1}, \hat{\sigma}_{\textup{HC}2}, \hat{\sigma}_{\textup{HC}3}$, and the true standard deviation of $\htau$.}
	\end{subfigure}
		\begin{subfigure}{0.99\textwidth}
                  \centering
		\includegraphics[width=0.48\textwidth]{synthetic_simul_normal_pi2_coverage.pdf}
		\includegraphics[width=0.48\textwidth]{synthetic_simul_normal_pi5_coverage.pdf}
		\caption{Empirical $95\%$ coverage of $t$-statistics derived from two estimators and four variance estimators (
                  ``theoretical'' for $\sigma_{n}^{2}$, ``HC2'' for $\hat{\sigma}_{\text{HC}2}^{2}$ and ``HC3''  for $\hat{\sigma}_{\text{HC}3}^{2}$)}
	\end{subfigure}
	\caption{Simulation. $X$ is a realization of a random matrix with i.i.d. $N(0, 1)$ entries and $e(t)$ is a realization of a random vector with i.i.d. entries: (Left) $\pi_{1} = 0.2$; (Right) $\pi_{1} = 0.5$. Each column corresponds to a distribution of $e(t)$. }\label{fig::simulation_normal_01}
\end{figure}

\begin{figure}
	\centering
	\begin{subfigure}{0.95\textwidth}
          \centering
		\includegraphics[width=0.48\textwidth]{synthetic_simul_debiase_normal_pi2_bias.pdf}
		\includegraphics[width=0.48\textwidth]{synthetic_simul_debiase_normal_pi5_bias.pdf}
		\caption{ Relative bias of $\hdtau$ and $\htau$.}\label{fig:normal_bias_worst}
	\end{subfigure}
		\begin{subfigure}{0.95\textwidth}
                  \centering
		\includegraphics[width=0.48\textwidth]{synthetic_simul_debiase_normal_pi2_coverage.pdf}
		\includegraphics[width=0.48\textwidth]{synthetic_simul_debiase_normal_pi5_coverage.pdf}
		\caption{Empirical $95\%$ coverage of $t$-statistics derived from two estimators and two variance estimators (``HC2'' for $\hat{\sigma}_{\text{HC}2}^{2}$ and ``HC3''  for $\hat{\sigma}_{\text{HC}3}^{2}$)}
	\end{subfigure}
	\caption{Simulation. $X$ is a realization of a random matrix with i.i.d. $N(0, 1)$ entries and $e(t)$ is defined in \eqref{eq:worst_pout}: (Left) $\pi_{1} = 0.2$; (Right) $\pi_{1} = 0.5$. }\label{fig::simulation_normal_02}
\end{figure}

\begin{figure}[htp]
	\centering
	\begin{subfigure}{0.99\textwidth}
          \centering
		\includegraphics[width=0.48\textwidth]{synthetic_simul_t1_pi2_bias.pdf}
		\includegraphics[width=0.48\textwidth]{synthetic_simul_t1_pi5_bias.pdf}
		\caption{ Relative bias of $\hdtau$ and $\htau$.}
	\end{subfigure}
	\begin{subfigure}{0.99\textwidth}
          \centering
		\includegraphics[width=0.48\textwidth]{synthetic_simul_t1_pi2_sdinflate.pdf}
		\includegraphics[width=0.48\textwidth]{synthetic_simul_t1_pi5_sdinflate.pdf}
		\caption{Ratio of standard deviation between five standard deviation estimates, $\sigma_{n}, \hat{\sigma}_{\textup{HC}0}, \hat{\sigma}_{\textup{HC}1}, \hat{\sigma}_{\textup{HC}2}, \hat{\sigma}_{\textup{HC}3}$, and the true standard deviation of $\htau$.}
	\end{subfigure}
		\begin{subfigure}{0.99\textwidth}
                  \centering
		\includegraphics[width=0.48\textwidth]{synthetic_simul_t1_pi2_coverage.pdf}
		\includegraphics[width=0.48\textwidth]{synthetic_simul_t1_pi5_coverage.pdf}
		\caption{Empirical $95\%$ coverage of $t$-statistics derived from two estimators and four variance estimators (
                  ``theoretical'' for $\sigma_{n}^{2}$, ``HC2'' for $\hat{\sigma}_{\text{HC}2}^{2}$ and ``HC3''  for $\hat{\sigma}_{\text{HC}3}^{2}$)}
	\end{subfigure}
	\caption{Simulation. $X$ is a realization of a random matrix with i.i.d. $t(1)$ entries and $e(t)$ is a realization of a random vector with i.i.d. entries: (Left) $\pi_{1} = 0.2$; (Right) $\pi_{1} = 0.5$. Each column corresponds to a distribution of $e(t)$. }\label{fig::simulation_t1_01}
\end{figure}

\begin{figure}
	\centering
	\begin{subfigure}{0.85\textwidth}
          \centering
		\includegraphics[width=0.48\textwidth]{synthetic_simul_debiase_t1_pi2_bias.pdf}
		\includegraphics[width=0.48\textwidth]{synthetic_simul_debiase_t1_pi5_bias.pdf}
		\caption{ Relative bias of $\hdtau$ and $\htau$.}
	\end{subfigure}
		\begin{subfigure}{0.85\textwidth}
                  \centering
		\includegraphics[width=0.48\textwidth]{synthetic_simul_debiase_t1_pi2_coverage.pdf}
		\includegraphics[width=0.48\textwidth]{synthetic_simul_debiase_t1_pi5_coverage.pdf}
		\caption{Empirical $95\%$ coverage of $t$-statistics derived from two estimators and two variance estimators (``HC2'' for $\hat{\sigma}_{\text{HC}2}^{2}$ and ``HC3''  for $\hat{\sigma}_{\text{HC}3}^{2}$)}
	\end{subfigure}
	\caption{Simulation. $X$ is a realization of a random matrix with i.i.d. $t(1)$ entries and $e(t)$ is defined in \eqref{eq:worst_pout}: (Left) $\pi_{1} = 0.2$; (Right) $\pi_{1} = 0.5$. }\label{fig::simulation_t1_02}
\end{figure}

\newpage
\subsection{Experimental results on real datasets}\label{subapp:experiments3}

\subsubsection{The LaLonde data}
We use the dataset from a randomized experiment on evaluating the impact of National Supported Work Demonstration, a labor training program, on postintervention income levels \citep{lalonde86, dehejia99}. It is available at \url{http://users.nber.org/~rdehejia/data/nswdata2.html}, and has $n = 445$ units with $n_{1} = 185$ units assigned in the program. It has 10 basic covariates: age, education, Black (1 if black, 0 otherwise), Hispanic (1 if Hispanic, 0 otherwise), married (1 if married, 0 otherwise), nodegree (1 if no degree, 0 otherwise), RE74/RE75 (earnings in 1974/1975), u74/u75 (1 if RE74/RE75 = 0, 0 otherwise). We form a $445\times 49$ $X$ by including all covariates and two-way interaction terms, and removing the ones perfectly collinear with others. We generate potential outcomes which mimics the truth. Specifically, we first regress the observed outcomes on the covariates in each group separately to obtain the coefficient vectors $\hat{\beta}_{1}, \hat{\beta}_{0}\in \R^{49}$ and the estimates $\hat{\sigma}_{1}, \hat{\sigma}_{0}$ of error standard deviation.

For each $p \in \{1, 2, \ldots, 49\}$, we randomly extract $p$ columns to form a $445 \times p$ submatrix. Then we generate potential outcomes from \eqref{eq:simul_potential_outcomes} by setting $\hat{\beta}_{1}, \hat{\beta}_{0}$ to be the subvector of $\beta_{1}, \beta_{0}$ corresponding to the positions of selected columns and setting $\sigma_{1} = \hat{\sigma}_{1} / 2$ and $\sigma_{0} = \hat{\sigma}_{0} / 2$.
Then we perform all steps as for the synthetic datasets before. For each $p$ we repeat the above procedure using $50$ random seeds and report the median of all measures. Fig. \ref{fig::real_lalonde_01} and Fig. \ref{fig::real_lalonde_02} show the results.

Compared to the synthetic dataset in Section \ref{sec:experiment}, this dataset is more adversarial to our theory in that even the HC3 variance estimate suffers from undercoverage for large $p$. It turns out that $\kappa = 0.887$ in this dataset while $\kappa = 0.184$ for random matrices with i.i.d. $N(0,1)$ entries.

\subsubsection{The STAR data}
The second dataset is from the Student Achievement and Retention (STAR) Project, a randomized evaluation of academic services and incentives on college freshmen. It has 974 units with 118 units assigned to the treatment group. \cite{angrist09} give more details. We include gender, age, high school GPA, mother language, indicator on whether living at home, frequency on putting off studying for tests, education, mother education, father education, intention to graduate in four years and indicator whether being at the preferred school. We also include the interaction terms between age, gender, high school GPA and all other variables. This ends up with $53$ variables. Fig. \ref{fig::real_star_01} and  Fig. \ref{fig::real_star_02} show the results.

\begin{figure}[htp]
	\centering
	\begin{subfigure}{0.99\textwidth}
          \centering
		\includegraphics[width=0.85\textwidth]{real_simul_lalonde_bias.pdf}
		\caption{ Relative bias of $\hdtau$ and $\htau$.}
	\end{subfigure}
	\begin{subfigure}{0.99\textwidth}
          \centering
		\includegraphics[width=0.85\textwidth]{real_simul_lalonde_sdinflate.pdf}
		\caption{Ratio of standard deviation between five standard deviation estimates, $\sigma_{n}, \hat{\sigma}_{\textup{HC}0}, \hat{\sigma}_{\textup{HC}1}, \hat{\sigma}_{\textup{HC}2}, \hat{\sigma}_{\textup{HC}3}$, and the true standard deviation of $\htau$.}
	\end{subfigure}
		\begin{subfigure}{0.99\textwidth}
                  \centering
                  \includegraphics[width=0.85\textwidth]{real_simul_lalonde_coverage.pdf}
		\caption{Empirical $95\%$ coverage of $t$-statistics derived from two estimators and four variance estimators (
                  ``theoretical'' for $\sigma_{n}^{2}$, ``HC2'' for $\hat{\sigma}_{\text{HC}2}^{2}$ and ``HC3''  for $\hat{\sigma}_{\text{HC}3}^{2}$)}
	\end{subfigure}
	\caption{Simulation on Lalonde dataset. $e(t)$ is a realization of a random vector with i.i.d. entries. Each column corresponds to a distribution of $e(t)$. }\label{fig::real_lalonde_01}
\end{figure}

\begin{figure}[htp]
	\centering
	\begin{subfigure}{0.7\textwidth}
          \centering
		\includegraphics[width= \textwidth]{real_simul_debiase_lalonde_bias.pdf}
		\caption{ Relative bias of $\hdtau$ and $\htau$.}
	\end{subfigure}
		\begin{subfigure}[t]{0.7\textwidth}
                  \centering
		\includegraphics[width= \textwidth]{real_simul_debiase_lalonde_coverage.pdf}
		\caption{Empirical $95\%$ coverage of $t$-statistics derived from two estimators and two variance estimators (``HC2'' for $\hat{\sigma}_{\text{HC}2}^{2}$ and ``HC3''  for $\hat{\sigma}_{\text{HC}3}^{2}$)}
	\end{subfigure}

	\caption{Simulation on Lalonde dataset. $e(t)$ is defined in \eqref{eq:worst_pout}.}\label{fig::real_lalonde_02}
\end{figure}

\begin{figure}[htp]
	\centering
	\begin{subfigure}{0.99\textwidth}
          \centering
		\includegraphics[width=0.85\textwidth]{real_simul_star_bias.pdf}
		\caption{ Relative bias of $\hdtau$ and $\htau$.}
	\end{subfigure}
	\begin{subfigure}{0.99\textwidth}
          \centering
		\includegraphics[width=0.85\textwidth]{real_simul_star_sdinflate.pdf}
		\caption{Ratio of standard deviation between five standard deviation estimates, $\sigma_{n}, \hat{\sigma}_{\textup{HC}0}, \hat{\sigma}_{\textup{HC}1}, \hat{\sigma}_{\textup{HC}2}, \hat{\sigma}_{\textup{HC}3}$, and the true standard deviation of $\htau$.}
	\end{subfigure}
		\begin{subfigure}{0.99\textwidth}
                  \centering
                  \includegraphics[width=0.85\textwidth]{real_simul_star_coverage.pdf}
		\caption{Empirical $95\%$ coverage of $t$-statistics derived from two estimators and four variance estimators (
                  ``theoretical'' for $\sigma_{n}^{2}$, ``HC2'' for $\hat{\sigma}_{\text{HC}2}^{2}$ and ``HC3''  for $\hat{\sigma}_{\text{HC}3}^{2}$)}
	\end{subfigure}
	\caption{Simulation on STAR dataset. $e(t)$ is a realization of a random vector with i.i.d. entries. Each column corresponds to a distribution of $e(t)$. }\label{fig::real_star_01}
\end{figure}

\begin{figure}[htp]
	\centering
	\begin{subfigure}{0.7\textwidth}
          \centering
		\includegraphics[width= \textwidth]{real_simul_debiase_star_bias.pdf}
		\caption{ Relative bias of $\hdtau$ and $\htau$.}
	\end{subfigure}
		\begin{subfigure}[t]{0.7\textwidth}
                  \centering
		\includegraphics[width= \textwidth]{real_simul_debiase_star_coverage.pdf}
		\caption{Empirical $95\%$ coverage of $t$-statistics derived from two estimators and two variance estimators (``HC2'' for $\hat{\sigma}_{\text{HC}2}^{2}$ and ``HC3''  for $\hat{\sigma}_{\text{HC}3}^{2}$)}
	\end{subfigure}

	\caption{Simulation on STAR dataset. $e(t)$ is defined in \eqref{eq:worst_pout}.}\label{fig::real_star_02}
\end{figure}

\end{document}

config: arxiv