# Adjustments with Many Regressors under Covariate-Adaptive
Randomizations††thanks: Ke Miao thanks the National Natural Science Foundation
of China for financial support under grant number 72103046. Yichong Zhang
acknowledges the financial support from the Lee Kong Chian fellowship and the
NSFC under the grant No. 72133002. Any and all errors are our own.

Liang Jiang Fudan University. E-mail address<EMAIL_ADDRESS>Liyao
Li The corresponding author. East China Normal University. E-mail address:
<EMAIL_ADDRESS> Ke Miao Fudan University. E-mail address:
<EMAIL_ADDRESS>Yichong Zhang Singapore Management University. E-mail
address<EMAIL_ADDRESS>
###### Abstract

Our paper discovers a new trade-off of using regression adjustments (RAs) in
causal inference under covariate-adaptive randomizations (CARs). On one hand,
RAs can improve the efficiency of causal estimators by incorporating
information from covariates that are not used in the randomization. On the
other hand, RAs can degrade estimation efficiency due to their estimation
errors, which are not asymptotically negligible when the number of regressors
is of the same order as the sample size. Ignoring the estimation errors of RAs
may result in serious over-rejection of causal inference under the null
hypothesis. To address the issue, we construct a new ATE estimator by
optimally linearly combining the adjusted and unadjusted estimators. We then
develop a unified inference theory for this estimator under CARs. It has two
features: (1) the Wald test based on it achieves the exact asymptotic size
under the null hypothesis, regardless of whether the number of covariates is
fixed or diverges no faster than the sample size; and (2) it guarantees weak
efficiency improvement over both the adjusted and unadjusted estimators.

Keywords: Covariate-adaptive randomization, many regressors, regression
adjustment.

JEL codes: C14, C21, D14, G21

##  1 Introduction

Covariate-adaptive randomizations (CARs) have recently seen growing use in a
wide variety of randomized experiments in economic research.111See, for
example, Chong et al. (2016); Greaney et al. (2016); Jakiela and Ozier (2016);
Burchardi et al. (2019); Anderson and McKenzie (2021), and etc. Under CARs,
units are first stratified using some baseline covariates, and within each
stratum, the treatment status is assigned independently to achieve the balance
between the numbers of treated and control units. Given data from CARs,
researchers often estimate and infer various treatment effect parameters via
regressions with strata dummies and baseline covariates as controls, a
practice known as regression adjustments (RAs). However, Freedman (2008a, b)
showed that OLS regression with covariates, as in the analysis of covariance
(ANCOVA), can actually decrease the precision of the average treatment effect
(ATE) estimator. Lin (2013) further examined Freedman’s critique and found
that for a linear RA to be “no-harm” (i.e., guaranteed efficiency improvement
over the unadjusted estimator), it must include a full set of interactions
between treatment status and covariates. Since then, many “no-harm” RAs have
been developed for various randomization schemes and causal parameters under
both finite- and super-population asymptotics.222See, for example, Shao et al.
(2010); Bloniarz et al. (2016); Fogarty (2018); Liu et al. (2023); Li and Ding
(2020); Ma et al. (2022); Cohen and Fogarty (2023); Lei and Ding (2021); Jiang
et al. (2023, 2023); Reluga et al. (2022); Ye et al. (2023, 2022); Bai et al.
(2023); Cytrynbaum (2023); Chiang et al. (2023), and etc.

“No-harm” RAs are usually coupled with high dimensional regressors in causal
inference with economic experiments, as noted by Duflo (2018).333Examples of
such experiments include Bursztyn et al. (2019); Banerjee et al. (2020); Dhar
et al. (2022), among others. Moreover, when the regressors are sieve bases
with a growing dimension, causal estimators with “no-harm” RAs can potentially
achieve the semiparametric efficiency bound, as shown by Jiang et al. (2023)
and Bai et al. (2023). However, it is found that including many regressors in
RAs can distort the size of causal inference in finite sample. This issue of
dimensionality is further exacerbated under CARs, as the “no-harm” RAs under
CARs are better performed stratum by stratum for efficiency gains ((Ye et al.,
2022)).

This paper studies linear RAs for the estimation and inference of the average
treatment effect (ATE) under CARs in the “many regressors” regime–an
asymptotic regime where the number of regressors can diverge at most as fast
as the sample size. We discover a new efficiency trade-off of using “no-harm”
RAs due to the dimensionality of regressors. On one hand, RAs can improve the
estimation efficiency by incorporating information from covariates that are
not used in the randomization; on the other hand, they can degrade estimation
efficiency due to their estimation errors, which are not asymptotically
negligible when the number of regressors grows at the same rate of the sample
size. This trade-off poses potentially serious challenges for inference.
First, ignoring the cost of RAs in the asymptotic variance estimation may lead
to substantial over-rejection under the null. Second, a consistent variance
estimator that accounts for the cost of RAs can potentially be larger than
that of the simple estimator without RAs, even asymptotically, if the cost
outweighs the benefit. Therefore, using a “no-harm” RA that incurs such costs
can actually be harmful.

Under the “many regressors” asymptotic regime, we derive the joint asymptotic
distribution for both the adjusted and unadjusted estimators. We also propose
a consistent estimator of the asymptotic covariance matrix, which accounts for
both the benefit and cost of RAs. We use this matrix to construct a new ATE
estimator by optimally linearly combining the adjusted and unadjusted
estimators. We show that (1) the Wald test based on this new optimal linear
combination estimator achieves the exact asymptotic size under the null and
(2) the optimal linear combination estimator is weakly more efficient than
both the adjusted and unadjusted estimators. Furthermore, we conduct a local
asymptotic power analysis and establish that the Wald test based on the
optimal linear combination estimator is asymptotically the uniformly most
powerful (UMP) test against two-sided local alternatives over a class of
unbiased tests that only depend on the adjusted and unadjusted estimators.
This implies our new estimator is optimal locally among not only linear but
also nonlinear combination estimators. These results collectively address the
challenge of the efficiency trade-off unearthed in the paper.

We further consider an alternative asymptotic regime where the dimension of
regressors is fixed. Under mild regularity conditions, we show that the same
optimal linear combination estimator with the same covariance matrix estimator
achieves the same properties (1) and (2) as above. This implies that the new
inference method proposed in this paper has uniform size control, and
empirical researchers can use our inference method without predetermining
whether the number of regressors is fixed or diverges to infinity.

Our results build on the work of Cattaneo et al. (2018), who studied the
inference of linear coefficients in OLS regression with many regressors and
independent or cluster independent observations. We extend their analysis in
three new directions. First, we consider observations generated by CARs, which
introduce cross-sectional dependence among treatment assignments and outcomes.
Second, we estimate the ATE in a two-step procedure by linearly combining
intercepts from multiple OLS regressions with many regressors. Because the
first step intercept estimators are asymptotically dependent and the second
step of the procedure introduces additional estimation errors, the asymptotic
distribution of our ATE estimator is not a simple application of the result by
Cattaneo et al. (2018). Instead, we develop a new distributional theory for
our ATE estimator using techniques from Bugni et al. (2018), the Yurinskii’s
coupling, and an anti-concentration inequality from Chernozhukov et al.
(2014). Third, we propose a new optimal linear combination estimator and
discuss its efficiency improvement under the many regressors framework, which
is novel in the literature. Cattaneo et al. (2019) also considered a two-step
procedure, and their knife-edge order of the number of regressors is the
square root of the sample size. Our paper differs from theirs because our
first step intercept estimators enter the second step in a linear fashion,
resulting in a knife-edge rate of the same order of the sample size.

Our covariance matrix estimator is built upon the cross-fit estimators of both
the variance of regression coefficients and the “variance component” with many
regressors, as proposed by Jochmans (2022) and Kline et al. (2020),
respectively. Moreover, to account for the dependence among the first step
intercept estimators, we extend Kline et al. (2020) by proposing a new
estimator for the “covariance component” with many regressors, which involves
coefficients from two distinctive linear regressions.

This paper relates to other recent works on RAs in randomized experiments. Liu
et al. (2023), Ma et al. (2022), and Ye et al. (2022) considered linear RAs
for ATE under CARs when the (effective) number of covariates is either fixed
or diverges at a rate less than the square root of the sample size.444Liu et
al. (2023) consider the ultra-high dimension that the number of regressors can
be higher than the sample size and use Lasso to select the effective
regressors, which is assumed to be less than the square root of the sample
size. In their settings, the estimation errors of RAs are asymptotically
negligible. Under a finite-population framework with the complete random
sampling, Li and Ding (2020) and Chiang et al. (2023) considered bias-
correction for the ATE allowing the number of regressors to grow at a rate
that is slower than the sample size. In contrast, we consider a super-
population framework with CARs, in which the number of regressors can grow at
the same rate as the sample size.

This paper also adds to the literature on RAs with high-dimensional sparse
covariates. For example, Wager et al. (2016), Liu et al. (2023), Jiang et al.
(2023), and Jiang et al. (2023) consider the regularized regression
adjustments for the inference of average, quantile, and local average
treatment effects, respectively, under a sparsity condition. However, as
argued by Li and Müller (2021), the sparsity condition (1) may not hold in
social science applications as it is unclear why the large majority of control
coefficients should be very nearly zero, (2) is not invariant to linear
transformations of the controls, and (3) depends on some tuning parameter
(e.g., the penalty term in Lasso), which complicates small sample
interpretation of the resulting inference. Our approach does not need the
sparsity condition, and thus, avoids these issues.

The rest of this paper is organized as follows. Section 2 introduces the
setup. Section 3 identifies the efficiency trade-off in RAs with many
regressors. In Section 4, we present the optimal linear combination estimator
for ATE and its asymptotic properties. Section 5 discusses the properties of
the optimal linear combination estimator with fixed number of regressors.
Sections 6 and 7 provide simulations and an empirical application,
respectively. Proofs and additional figures are given in the Online Appendix.

Notation. For any positive integer $m$, let $1_{m}$ and $I_{m}$ be the
$m\times 1$ vector of ones and the $m\times m$ identity matrix, respectively.
Let $||\cdot||_{2}$ and $||\cdot||_{op}$ denote the $\ell_{2}$-norm for a
vector and operator norm for a matrix, respectively. For a symmetric and
positive semi-definite matrix $\Upsilon$, $\lambda_{\max}(\Upsilon)$ denotes
the maximum eigenvalue of $\Upsilon$. We define $W^{(n)}$ as the sample of
$W$’s, i.e., $W^{(n)}=\\{W_{i}\\}_{i\in[n]}$, where $[n]=1,\cdots,n$. We write
$U\stackrel{{\scriptstyle d}}{{=}}V$ for two random variables $U$ and $V$ if
they share the same distribution.

##  2 Setup

Potential outcomes for treated and control groups are denoted by $Y(1)$ and
$Y(0)$, respectively. Treatment status is denoted by $A$, with $A=1$
indicating treated and $A=0$ untreated. The stratum indicator is denoted by
$S$, based on which the researcher implements the covariate-adaptive
randomization. The support of $S$ is denoted by $\mathcal{S}$, a finite set.
After randomization, the researcher can observe the data
$\\{Y_{i},S_{i},A_{i},X_{i}\\}_{i\in[n]}$ where
$Y_{i}=Y_{i}(1)A_{i}+Y_{i}(0)(1-A_{i})$ is the observed outcome, and $X_{i}$
contains covariates besides $S_{i}$ in the dataset. The support of $X$ is
denoted as $\text{Supp}(X)$. In this paper, we allow $X_{i}$ and $S_{i}$ to be
dependent. For $i\in[n]$, let $n_{s}=\sum_{i\in[n]}1\\{S_{i}=s\\}$,
$n_{1,s}=\sum_{i\in[n]}A_{i}1\\{S_{i}=s\\}$, and $n_{0,s}=n_{s}-n_{1,s}$. Let
$\aleph_{a,s}=\\{i\in[n]:A_{i}=a,S_{i}=s\\}$ denote the set of individuals in
stratum $s$ with treatment status $a$ and
$\aleph_{s}=\aleph_{1,s}\cup\aleph_{0,s}$. Our parameter of interest is the
average treatment effect defined as $\tau=\mathbb{E}(Y(1)-Y(0))$.

We make the following assumptions on the data generating process (DGP) and the
treatment assignment rule.

######  Assumption 1.

  1. (i)

$\\{Y_{i}(1),Y_{i}(0),S_{i},X_{i}\\}_{i\in[n]}$ are independent and
identically distributed (i.i.d.).

  2. (ii)

$\\{Y_{i}(1),Y_{i}(0),X_{i}\\}_{i\in[n]}\perp\\!\\!\\!\perp\\{A_{i}\\}_{i\in[n]}|\\{S_{i}\\}_{i\in[n]}$.

  3. (iii)

Suppose $p_{s}=\mathbb{P}(S_{i}=s)$ is fixed with respect to (w.r.t.) $n$ and
is positive for every $s\in\mathcal{S}$.

  4. (iv)

Let $\pi_{s}$ denote the target fraction of treatment for stratum $s$. Then,
$c<\min_{s\in\mathcal{S}}\pi_{s}\leq\max_{s\in\mathcal{S}}\pi_{s}<1-c$ for
some constant $c\in(0,0.5)$ and $\frac{D_{n,s}}{n_{s}}=o_{P}(1)$ for
$s\in\mathcal{S}$, where
$D_{n,s}=\sum_{i\in[n]}(A_{i}-\pi_{s})1\\{S_{i}=s\\}$.

######  Remark 1.

Several remarks are in order. First, Assumption 1(i) assumes
$(Y(1),Y(0),S,X)^{(n)}$ are independent but allows for $A^{(n)}$, and thus,
$Y^{(n)}$ to be cross-sectionally dependent, which will be the case under
CARs. The identical distribution assumption can be relaxed by using a set of
more complex notation. Second, Assumption 1(ii) implies that the treatment
assignment $A^{(n)}$ are generated only based on strata indicators. Third,
Assumption 1(iii) imposes that the number of strata is bounded and the strata
sizes are approximately balanced. Fourth, Bugni et al. (2018) show that
Assumption 1(iv) holds under several specific covariate-adaptive treatment
assignment rules such as simple random sampling (SRS), biased-coin design
(BCD), adaptive biased-coin design (WEI) and stratified block randomization
(SBR). For completeness, we provide brief descriptions below. Note that the
requirement $D_{n,s}/n_{s}=o_{P}(1)$ is weaker than the assumption imposed by
Bugni et al. (2018), but it is the same as that imposed by Bugni et al. (2019)
and Zhang and Zheng (2020).

######  Example 1 (SRS).

Let $\\{A_{i}\\}_{i=1}^{n}$ be drawn independently across $i$ and of
$\\{S_{i}\\}_{i=1}^{n}$ as Bernoulli random variables with success rate
$\pi_{s}$, i.e., for $k=1,\ldots,n$,

$\mathbb{P}\left(A_{k}=1\big{|}\\{S_{i}\\}_{i=1}^{n},\\{A_{j}\\}_{j=1}^{k-1}\right)=\mathbb{P}(A_{k}=1|S_{k})=\pi_{S_{k}}.$

######  Example 2 (WEI).

This design was first proposed by Wei (1978). Let
$n_{k-1}(S_{k})=\sum_{i=1}^{k-1}1\\{S_{i}=S_{k}\\}$,
$B_{k-1}(S_{k})=\sum_{i=1}^{k-1}\left(A_{i}-\frac{1}{2}\right)1\\{S_{i}=S_{k}\\}$,
and

$\displaystyle\mathbb{P}\left(A_{k}=1\big{|}\\{S_{i}\\}_{i=1}^{k},\\{A_{i}\\}_{i=1}^{k-1}\right)=f\biggl{(}\frac{2B_{k-1}(S_{k})}{n_{k-1}(S_{k})}\biggr{)},$

where $f(\cdot):[-1,1]\mapsto[0,1]$ is a pre-specified non-increasing function
satisfying $f(-x)=1-f(x)$. Here, $\frac{B_{0}(S_{1})}{n_{0}(S_{1})}$ and
$B_{0}(S_{1})$ are understood to be zero.

######  Example 3 (BCD).

The treatment status is determined sequentially for $1\leq k\leq n$ as

$\displaystyle\mathbb{P}\left(A_{k}=1|\\{S_{i}\\}_{i=1}^{k},\\{A_{i}\\}_{i=1}^{k-1}\right)=\begin{cases}\frac{1}{2}&\text{if
}B_{k-1}(S_{k})=0\\\ \lambda&\text{if }B_{k-1}(S_{k})<0\\\ 1-\lambda&\text{if
}B_{k-1}(S_{k})>0,\end{cases}$

where $B_{k-1}(s)$ is defined as above and $\frac{1}{2}<\lambda\leq 1$.

######  Example 4 (SBR).

For each stratum, $\lfloor\pi_{s}n_{s}\rfloor$ units are assigned to treatment
and the rest are assigned to control.

######  Remark 2.

We note that SRS and SBR allow for the target faction of treatment to be
different for different strata. For WEI and BCD, we have $\pi_{s}=1/2$ for
$s\in\mathcal{S}$. Bugni et al. (2018) further shows BCD and SBR achieve
strong balance in the sense that $D_{n,s}=o_{P}(n^{1/2})$.

##  3 Efficiency Trade-off in RAs with Many Regressors

This section derives the joint asymptotic distribution of the fully saturated
adjusted and unadjusted estimators under the assumption that the number of
regressors grows no faster than the sample size. We then identify the
efficiency trade-off in RAs based on the comparison of two asymptotic
variances. The variance-covariance matrix will also be used to construct the
optimal linear combination estimator in Section 4.

The fully saturated regression adjusted estimator is denoted by
$\hat{\tau}^{adj}$ and defined as

$\displaystyle\hat{\tau}^{adj}=\sum_{s\in\mathcal{S}}\hat{p}_{s}(\hat{\tau}_{1,s}-\hat{\tau}_{0,s}),$

where $\hat{\tau}_{a,s}$ is computed as the intercept in the OLS regression of

$\displaystyle Y_{i}\sim 1+\breve{X}_{i},$ (3.1)

using all observations in $\aleph_{a,s}$, $\hat{p}_{s}=n_{s}/n$, and
$\breve{X}_{i}=X_{i}-\overline{X}_{S_{i}}$ with
$\overline{X}_{s}=\frac{1}{n_{s}}\sum_{i\in\aleph_{s}}X_{i}$.

######  Remark 3.

We consider this adjusted estimator for two reasons. First, it is a natural
extension to the fully saturated estimator proposed by Bugni et al. (2019). Ye
et al. (2022) have also considered this same estimator (their
$\hat{\theta}_{A}$) and showed that it is “no-harm” when the dimension of
$X_{i}$ is fixed. Second, we can write $\hat{\tau}^{adj}$ as an augmented
inverse propensity score weighted (AIPW) estimator:

$\displaystyle\hat{\tau}^{adj}$
$\displaystyle=\frac{1}{n}\sum_{i\in[n]}\frac{A_{i}(Y_{i}-X_{i}^{\top}\hat{\beta}_{1,S_{i}})}{\hat{\pi}_{S_{i}}}-\frac{1}{n}\sum_{i\in[n]}\frac{(1-A_{i})(Y_{i}-X_{i}^{\top}\hat{\beta}_{0,S_{i}})}{1-\hat{\pi}_{S_{i}}}+\frac{1}{n}\sum_{i\in[n]}X_{i}^{\top}(\hat{\beta}_{1,S_{i}}-\hat{\beta}_{0,S_{i}}),$

where $\hat{\pi}_{s}=n_{1,s}/n_{s}$. It is important to note that, if an
individual $i$ belongs to stratum $s$, then $\breve{X}_{i}$ is defined as
$X_{i}$ demeaned by the stratum-level mean
$\frac{1}{n_{s}}\sum_{i\in\aleph_{s}}X_{i}$, not the stratum-treatment level
mean $\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}X_{i}$. This is the key for the
AIPW representation to hold.

The unadjusted estimator is denoted by $\hat{\tau}^{unadj}$ and defined as

$\displaystyle\hat{\tau}^{unadj}=\frac{1}{n}\sum_{i\in[n]}\frac{A_{i}Y_{i}}{\hat{\pi}_{S_{i}}}-\frac{1}{n}\sum_{i\in[n]}\frac{(1-A_{i})Y_{i}}{1-\hat{\pi}_{S_{i}}}=\sum_{s\in\mathcal{S}}\hat{p}_{s}(\hat{\tau}_{1,s}^{unadj}-\hat{\tau}_{0,s}^{unadj}),$

where
$\hat{\tau}_{a,s}^{unadj}=\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}Y_{i}$.

######  Remark 4.

It is natural to consider $\hat{\tau}^{unadj}$ in our context. It is simply
the fully saturated regression adjusted estimator with only the constant as
the regressor, which is first proposed by Bugni et al. (2019). It is also just
the AIPW estimator with an empty set of covariates. Therefore, the unadjusted
estimator $\hat{\tau}^{unadj}$ is the polar opposite of $\hat{\tau}^{adj}$ in
the sense that it does not use the information of $X_{i}$ at all, and thus,
neither benefits nor suffers from RAs.

###  3.1 Asymptotic Properties

Denote the dimension of $X$ as $k$. In the following, we derive the asymptotic
properties of $\hat{\tau}^{adj}$ and $\hat{\tau}^{unadj}$ jointly in the case
that $k=k_{n}$ increases no faster than the sample size $n$.

To clearly state our assumptions, we need to introduce more notation. Let
$\breve{X}_{\aleph_{a,s}}$ be an $n_{a,s}\times k_{n}$ matrix which is
constructed by stacking $\breve{X}_{i}^{\top}$ for $i\in\aleph_{a,s}$. We
define $A_{\aleph_{a,s}}\in\Re^{n_{a,s}}$ similarly. Then, let
$P_{a,s}=\breve{X}_{\aleph_{a,s}}(\breve{X}_{\aleph_{a,s}}^{\top}\breve{X}_{\aleph_{a,s}})^{-1}\breve{X}_{\aleph_{a,s}}^{\top}$
be the projection matrix of $\breve{X}_{\aleph_{a,s}}$, and
$M_{a,s}=I_{n_{a,s}}-P_{a,s}$. Let $M_{a,s,i,j}$ and $P_{a,s,i,j}$ be the
$(i,j)$th entry of $M_{a,s}$ and $P_{a,s}$, respectively. Define

$\displaystyle\varepsilon_{i}(a)=Y_{i}(a)-\mathbb{E}(Y_{i}(a)|X_{i},S_{i}),$
$\displaystyle\gamma_{a,s,n}=(1_{n_{a,s}}^{\top}M_{a,s}1_{n_{a,s}})/n_{a,s},$
$\displaystyle\sigma^{2}_{a,s,n}=\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}\left[\left(\sum_{j\in\aleph_{a,s}}M_{a,s,i,j}\right)^{2}\mathbb{E}(\varepsilon_{i}^{2}(a)|X_{i},S_{i}=s)\right],\quad\text{and}$
$\displaystyle\rho_{a,s,n}=\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}\left[\left(\sum_{j\in\aleph_{a,s}}M_{a,s,i,j}\right)\mathbb{E}(\varepsilon_{i}^{2}(a)|X_{i},S_{i}=s)\right].$

######  Assumption 2.

  1. (i)

$\max_{i\in[n]}\mathbb{E}\left[\varepsilon_{i}^{4}(a)|X_{i},S_{i}\right]=O_{P}(1).$

  2. (ii)

There exists a constant $b>0$ such that
$\min_{a=0,1,s\in\mathcal{S},i\in[n]}\mathbb{E}\left[\varepsilon_{i}^{2}(a)|X_{i},S_{i}=s\right]\geq
b.$

  3. (iii)

$k_{n}/n_{a,s}\rightarrow\kappa_{a,s}\in\left[0,1\right)$ and
$\limsup_{n\rightarrow\infty}\max_{a=0,1,s\in\mathcal{S}}P_{a,s,i,i}\leq
1-\delta$ for some constant $\delta\in(0,1)$.

  4. (iv)

For $a=0,1$ and $s\in\mathcal{S}$, we have

$\displaystyle\mathbb{E}(Y_{i}(a)|X_{i},S_{i}=s)=\alpha_{a,s}+X_{i}^{\top}\beta_{a,s}+e_{i,s}(a),$

such that $\mathbb{E}(e_{i,s}^{2}(a)|S_{i}=s)=o(n^{-1})$.

  5. (v)

$\max_{a=0,1,s\in\mathcal{S}}\max_{i\in\aleph_{a,s}}\left|\sum_{j\in\aleph_{a,s}}M_{a,s,i,j}\right|=o_{P}(n^{1/2})$.

######  Remark 5.

Assumption 2(i) and 2(ii) are mild regularity conditions. Assumption 2(iii)
implies that we allow the number of covariates to diverge at the rate of the
sample size $n$. Because we run the stratum-treatment level regression with an
effective sample size of $n_{a,s}$, we require $\kappa_{a,s}<1$ to avoid
multicollinearity. Li and Ding (2020) and Chiang et al. (2023) consider the
inference of ATE under complete randomization in the finite-population setup
and require $\kappa_{a,s}=0$ (in fact, $k_{n}=o(n^{1/3})$) to eliminate the
asymptotic bias in their analysis. Instead, we consider the super-population
framework and assume the approximately correct specification of the linear
regression (Assumption 2(iv)). Under this assumption, the issue of
dimensionality of regressors is an inflation of asymptotic variance instead of
the bias. The approximately correct specification of the linear regression
(Assumption 2(iv)) is commonly assumed in linear regression analyses with many
regressors, as seen in Cattaneo et al. (2018, Assumption 3) and Jochmans
(2022, Assumption 3). Kline et al. (2020) further require Assumption 2(iv) to
hold with no approximation error, i.e., $e_{i,s}(a)=0$. This condition is
reasonable because we allow the dimension $k_{n}$ to diverge to infinity. For
example, if baseline covariates, denoted as $Z_{i}$, are fixed-dimensional and
continuous, then Assumption 2(iv) holds when $X_{i}$ contains sieve basis
functions of $Z_{i}$, and $\mathbb{E}(Y_{i}(a)|Z_{i}=z,S_{i}=s)$ is
sufficiently smooth in $z$.555In this case, we have
$\mathbb{E}(Y_{i}(a)|Z_{i},S_{i}=s)=\alpha_{a,s}+X_{i}^{\top}\beta_{a,s}+\tilde{e}_{i,s}(a)$
such that $\mathbb{E}(\tilde{e}_{i,s}^{2}(a)|S_{i}=s)=o(n^{-1})$. Then, we
have
$\mathbb{E}(Y_{i}(a)|X_{i},S_{i}=s)=\alpha_{a,s}+X_{i}^{\top}\beta_{a,s}+e_{i,s}(a)$
where $e_{i,s}(a)=\mathbb{E}(\tilde{e}_{i,s}(a)|X_{i},S_{i}=s)$. Then, by
Jansen’s inequality, we have
$\mathbb{E}(e_{i,s}^{2}(a)|S_{i}=s)\leq\mathbb{E}(\tilde{e}_{i,s}^{2}(a)|S_{i}=s)=o(n^{-1})$.
When all baseline covariates are categorical (discrete), Assumption 2(iv)
holds when $X_{i}$ contains the fully saturated dummies for all categories.
This condition is similar in spirit to the approximate sparsity condition that
is widely assumed in the literature of causal inference with ultra-high-
dimensional data. See, for example, Assumption 4.2 in the seminal work by
Belloni et al. (2017). Although this assumption is reasonable when the
dimension of $X_{i}$ ($k_{n}$) diverges to infinity, it is not when the
dimension is fixed. However, in Section 5, we show that, in the case with a
fixed $k$, our estimator and inference methods are still valid even when the
linear regression is misspecified, and the approximation error is not
asymptotically negligible. We will provide more discussion on this point
below. Assumption 2(v) is the same as Cattaneo et al. (2018, Assumption 3) and
Jochmans (2022, Assumption 3),666In fact, $\hat{v}_{i,n}$ in both Cattaneo et
al. (2018, Assumption 3) and Jochmans (2022, Assumption 3) is just
$\sum_{j\in\aleph_{a,s}}M_{a,s,i,j}$. which we recommend readers refer to for
more discussion.

######  Assumption 3.

Suppose $\gamma_{a,s,n}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\gamma_{a,s,\infty}>0$,
$\gamma_{a,s,n}^{-2}\sigma^{2}_{a,s,n}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\omega^{2}_{a,s,\infty}>0$, and
$\gamma_{a,s,n}^{-1}\rho_{a,s,n}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\varpi_{a,s,\infty}$ for some deterministic constants
$(\gamma_{a,s,\infty},\omega_{a,s,\infty},\varpi_{a,s,\infty})$.

######  Remark 6.

If $k_{n}\log k_{n}=o(n)$ (which means $\kappa_{a,s}=0$ for all $(a,s)$), then
under general conditions on the distribution of $X_{i}$, it is possible to
show that Assumption 3 holds with $\gamma_{a,s,\infty}=1$ and
$\omega^{2}_{a,s,\infty}=\varpi_{a,s,\infty}=\mathbb{E}(\varepsilon_{i}^{2}(a)|S_{i}=s)$.
See Belloni et al. (2015) and Cattaneo et al. (2020) for more discussion and
examples. If $\kappa_{a,s}>0$, the eigenvalues of gram matrix
$n^{-1}\breve{X}_{\aleph_{a,s}}^{\top}\breve{X}_{\aleph_{a,s}}$ do not
converge to those of its expectation. Instead, they will converge to some
spectral density which determines the quantities
$(\gamma_{a,s,\infty},\omega_{a,s,\infty},\varpi_{a,s,\infty})$. Specifically,
in Section A in the Online Supplement, we provide detailed calculation of
$\gamma_{a,s,\infty}$ when the $k_{n}\times 1$ vector $X_{i}$ given $S_{i}=s$
is jointly Gaussian with a covariance matrix $\Sigma_{s,n}$. Based on the
Marčenko-Pastur theorem ((Marčenko and Pastur, 1967)), we can show that

$\displaystyle\gamma_{a,s,\infty}=\frac{1}{1+(a(1-\pi_{s})+(1-a)\pi_{s})\zeta_{a,s}}<1,$

where
$\zeta_{a,s}=\int_{\lambda_{-}}^{\lambda_{+}}\frac{\sqrt{(\lambda_{+}-\lambda)(\lambda-\lambda_{-})}}{2\pi\lambda^{2}}d\lambda$
and $\lambda_{\pm}=(1\pm\sqrt{\kappa_{a,s}})^{2}$. If we further assume
homoskedasticity in the sense that
$\mathbb{E}(\varepsilon_{i}^{2}(a)|X_{i},S_{i}=s)=\mathbb{E}(\varepsilon_{i}^{2}(a)|S_{i}=s)$,
then, we have

$\displaystyle\omega^{2}_{a,s,\infty}=\varpi_{a,s,\infty}=\gamma_{a,s,\infty}^{-1}\mathbb{E}(\varepsilon_{i}^{2}(a)|S_{i}=s)>\mathbb{E}(\varepsilon_{i}^{2}(a)|S_{i}=s).$

We believe that Assumption 3 holds for non-Gaussian distributions of $X_{i}$
and heteroskedastic errors. However, establishing general preliminary
conditions for this assumption requires the use of advanced tools from random
matrix theory.777See, for example, Bai (2008) for a survey. Such a discussion
is beyond the scope of this paper. Finally, we emphasize that researchers do
not need to know the exact values of
$(\gamma_{a,s,\infty},\omega^{2}_{a,s,\infty},\varpi_{a,s,\infty})$ to
implement our inference method.

######  Theorem 3.1.

Suppose Assumptions 1–3 hold. Then, we have

$\displaystyle\sqrt{n}\begin{pmatrix}\hat{\tau}^{adj}-\tau\\\
\hat{\tau}^{unadj}-\tau\end{pmatrix}\rightsquigarrow\begin{pmatrix}\mathcal{U}^{adj}\\\
\mathcal{U}^{unadj}\\\ \end{pmatrix}+\begin{pmatrix}\mathcal{V}^{adj}\\\
\mathcal{V}^{unadj}\\\ \end{pmatrix}+\begin{pmatrix}\mathcal{W}\\\
\mathcal{W}\end{pmatrix},$

where

$\displaystyle\mathcal{U}=\begin{pmatrix}\mathcal{U}^{adj}\\\
\mathcal{U}^{unadj}\\\ \end{pmatrix}\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}\left(\begin{pmatrix}0\\\
0\end{pmatrix},\Sigma_{\mathcal{U}}\right),~{}\Sigma_{\mathcal{U}}=\begin{pmatrix}\mathbb{E}\left(\frac{\omega_{1,S_{i},\infty}^{2}}{\pi_{S_{i}}}+\frac{\omega_{0,S_{i},\infty}^{2}}{1-\pi_{S_{i}}}\right)&\mathbb{E}\left(\frac{\varpi_{1,S_{i},\infty}}{\pi_{S_{i}}}+\frac{\varpi_{0,S_{i},\infty}}{1-\pi_{S_{i}}}\right)\\\
\mathbb{E}\left(\frac{\varpi_{1,S_{i},\infty}}{\pi_{S_{i}}}+\frac{\varpi_{0,S_{i},\infty}}{1-\pi_{S_{i}}}\right)&\mathbb{E}\left(\frac{\mathbb{E}(\varepsilon_{i}^{2}(1)|S_{i})}{\pi_{S_{i}}}+\frac{\mathbb{E}(\varepsilon_{i}^{2}(0)|S_{i})}{1-\pi_{S_{i}}}\right)\end{pmatrix},$
$\displaystyle\mathcal{V}=\begin{pmatrix}\mathcal{V}^{adj}\\\
\mathcal{V}^{unadj}\\\ \end{pmatrix}\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}\left(\begin{pmatrix}0\\\
0\end{pmatrix},\Sigma_{\mathcal{V}}\right),~{}\Sigma_{\mathcal{V}}=\begin{pmatrix}\mathbb{E}var(\phi_{i}(1)-\phi_{i}(0)|S_{i})&\mathbb{E}var(\phi_{i}(1)-\phi_{i}(0)|S_{i})\\\
\mathbb{E}var(\phi_{i}(1)-\phi_{i}(0)|S_{i})&\mathbb{E}\left[\frac{var(\phi_{i}(1)|S_{i})}{\pi_{S_{i}}}+\frac{var(\phi_{i}(0)|S_{i})}{(1-\pi_{S_{i}})}\right]\end{pmatrix},$
$\displaystyle\mathcal{W}\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\Sigma_{\mathcal{W}}),~{}\Sigma_{\mathcal{W}}=var(\mathbb{E}(Y_{i}(1)-Y_{i}(0)|S_{i})),$

$\phi_{i}(a)=\mathbb{E}(Y_{i}(a)|X_{i},S_{i})-\mathbb{E}(Y_{i}(a)|S_{i})$ for
$a=0,1$, and $(U,V,W)$ are independent.

######  Remark 7.

Denote
$\Sigma=\Sigma_{\mathcal{U}}+\Sigma_{\mathcal{V}}+\Sigma_{\mathcal{W}}1_{2}1_{2}^{\top}$.
We will focus on the variance components of $\Sigma$ to identify the
efficiency trade-off in RAs herein. The whole matrix will be used to construct
the optimal linear combination estimator in Section 4.

###  3.2 Efficiency Trade-off

Theorem 3.1 decomposes the variance of ATE estimators into
$\Sigma_{\mathcal{U}}$, $\Sigma_{\mathcal{V}}$, and $\Sigma_{\mathcal{W}}$,
which represent the variation of $Y(a)$ given $X_{i}$ and $S_{i}$, the
variation of $X_{i}$ given $S_{i}$, and the variation of $S_{i}$,
respectively. Following the argument in Remark 6, we can show that
$Var(\mathcal{U}^{adj})=Var(\mathcal{U}^{unadj})$ if $k_{n}\log k_{n}=o(n)$.
Moreover, it can be shown that $Var(\mathcal{V}^{adj})\leq
Var(\mathcal{V}^{unadj})$, which indicates the benefit of RA from
incorporating information of $X_{i}$. Therefore, when $k_{n}\log k_{n}=o(n)$,
$\hat{\tau}^{adj}$ is weakly more efficient than $\hat{\tau}^{unadj}$, a
result that the previous literature has established under more stringent
conditions.

However, when $k_{n}$ is of the same order of $n$, $Var(\mathcal{U}^{adj})$
can be larger than $Var(\mathcal{U}^{unadj})$. To illustrate this, let’s
continue with the Gaussian covariates example mentioned in Remark 6. Suppose
for $a=0,1$ and $s\in\mathcal{S}$, $\pi_{s}=1/2$, $\kappa_{a,s}=\kappa$, and
$\varepsilon_{i}(a)$ is homoskedastic in the sense that
$\mathbb{E}(\varepsilon_{i}^{2}(a)|X_{i},S_{i}=s)=\mathbb{E}(\varepsilon_{i}^{2}(a)|S_{i}=s)$,
then we have

$\displaystyle\text{Variance Inflate
Factor}\equiv\frac{\text{Var}(\mathcal{U}^{adj})-\text{Var}(\mathcal{U}^{unadj})}{\text{Var}(\mathcal{U}^{unadj})}=\frac{\zeta}{2},$

where
$\zeta=\int_{\lambda_{-}}^{\lambda_{+}}\frac{\sqrt{(\lambda_{+}-\lambda)(\lambda-\lambda_{-})}}{2\pi\lambda^{2}}d\lambda$
and $\lambda_{\pm}=(1\pm\sqrt{\kappa})^{2}$. Figure 1 plots the values of the
variance inflation factor (VIF) as a function of $\kappa$. We can see that the
VIF is more than $12.5\%$, $25\%$, $50\%$, and $100\%$ if $\kappa$ is greater
than $1/5$, $1/3$, $1/2$, $2/3$, respectively.

Figure 1: The Variance Inflation Factor

Consequently, an efficiency trade-off in RAs between $Var(\mathcal{U}^{adj})$
and $Var(\mathcal{V}^{adj})$ can be identified in Theorem 3.1. This trade-off
can be elucidated by decomposing the difference between the asymptotic
variances of $\hat{\tau}^{adj}$ and $\hat{\tau}^{unadj}$ as

$\displaystyle(var(\mathcal{U}^{adj})-var(\mathcal{U}^{unadj}))+(var(\mathcal{V}^{adj})-var(\mathcal{V}^{unadj})),$

where the second term is always non-positive but the first term can be
positive when $k_{n}$ is of the same order of the sample size. Under some
circumstances, the first term may dominate the second, resulting in the
estimator with the theoretically “no-harm” adjustment being even less
efficient than the unadjusted estimator. For example, see Model 1 of Panel (b)
in Figures 2 and 3. This defeats the purpose of using RAs.

This also implies that the asymptotic variance estimator of $\hat{\tau}^{adj}$
in the literature, e.g., the one proposed by Ye et al. (2022), tends to under-
estimate the true asymptotic variance when there are many regressors, as it
neglects the VIF induced by the dimensionality. Therefore, when there are many
regressors, the causal inference using such a variance estimator may over-
reject under the null hypothesis, as illustrated in Section 6.

##  4 Optimal Linear Combination Estimator

As shown above, ignoring the VIF can result in size distortion in causal
inference. When we account for the VIF, however, the adjusted estimator may
become even less efficient than the unadjusted one due to the estimation
errors incurred in RAs with many regressors. In this section, we propose an
optimal linear combination estimator, denoted by $\hat{\tau}^{*}$, designed to
address both issues. We show that $\hat{\tau}^{*}$ is consistent,
asymptotically normal, and weakly more efficient than both $\hat{\tau}^{adj}$
and $\hat{\tau}^{unadj}$. We also present the local asymptotic power analysis
for the standard Wald test based on $\hat{\tau}^{*}$. Finally, we provide a
consistent estimator $\hat{\Sigma}$ for
$\Sigma=\Sigma_{\mathcal{U}}+\Sigma_{\mathcal{V}}+\Sigma_{\mathcal{W}}1_{2}1_{2}^{\top}$
in Section 4.3 below.

###  4.1 Asymptotic Properties of $\hat{\tau}^{*}$

We construct the optimal linear combination estimator as

$\displaystyle\hat{\tau}^{*}=\hat{w}\hat{\tau}^{adj}+(1-\hat{w})\hat{\tau}^{unadj},$
(4.1)

which is a weighted average of $\hat{\tau}^{adj}$ and $\hat{\tau}^{unadj}$
with the weight $\hat{\omega}$ minimizing the asymptotic variance.

Denote $\hat{\Sigma}=\begin{pmatrix}\hat{\Sigma}_{1,1}&\hat{\Sigma}_{1,2}\\\
\hat{\Sigma}_{1,2}&\hat{\Sigma}_{2,2}\end{pmatrix}$ as a consistent estimator
for $\Sigma=\begin{pmatrix}\Sigma_{1,1}&\Sigma_{1,2}\\\
\Sigma_{1,2}&\Sigma_{2,2}\end{pmatrix}$ and suppose $\hat{\tau}^{adj}$ and
$\hat{\tau}^{unadj}$ are not asymptotically equivalent, i.e.,
$\Sigma_{1,1}+\Sigma_{2,2}-2\Sigma_{1,2}>0$. Then, the optimal weight is

$\displaystyle\hat{w}=\frac{\hat{\Sigma}_{2,2}-\hat{\Sigma}_{1,2}}{\hat{\Sigma}_{1,1}+\hat{\Sigma}_{2,2}-2\hat{\Sigma}_{1,2}},$
(4.2)

and the corresponding estimator of the asymptotic variance for
$\hat{\tau}^{*}$ is
$\hat{\sigma}^{2}_{*}=(\hat{w},1-\hat{w})\hat{\Sigma}(\hat{w},1-\hat{w})^{\top}$.

######  Remark 8.

When $\Sigma_{1,1}-2\Sigma_{1,2}+\Sigma_{2,2}>0$, the adjusted estimator
$\hat{\tau}^{adj}$ and the unadjusted estimator $\hat{\tau}^{unadj}$ are not
asymptotically equivalent. In this case, $(w,1-w)\Sigma(w,1-w)^{\top}$ has a
unique minimizer
$w^{*}=\frac{\Sigma_{2,2}-\Sigma_{1,2}}{\Sigma_{1,1}-2\Sigma_{1,2}+\Sigma_{2,2}}$,
and our $\hat{w}$ is a consistent estimator of $w^{*}$. If the estimation
error of the RA is asymptotically negligible and $X_{i}$ does not contain
useful information of $(Y_{i}(1),Y_{i}(0))$ in the sense that

$\displaystyle\mathbb{E}(Y_{i}(a)|X_{i},S_{i})=\mathbb{E}(Y_{i}(a)|S_{i}),$

then we have $\Sigma_{1,1}-2\Sigma_{1,2}+\Sigma_{2,2}=0$ (which implies
$\Sigma_{1,1}=\Sigma_{1,2}=\Sigma_{2,2}$). In this case, all linear
combinations are asymptotically equivalent. If this case is a realistic
concern, then we can always safeguard against the zero denominator by
introducing a small positive constant $\lambda$ and constructing the final
estimator as
$\tilde{\tau}=\tilde{w}\hat{\tau}^{adj}+(1-\tilde{w})\hat{\tau}^{unadj}$ with
$\tilde{w}=\frac{\hat{\Sigma}_{2,2}-\hat{\Sigma}_{1,2}}{\hat{\Sigma}_{1,1}-2\hat{\Sigma}_{1,2}+\hat{\Sigma}_{2,2}+\lambda}$.
It is possible to show that $\tilde{\tau}$ is always weakly more efficient
than $\hat{\tau}^{unadj}$ as long as $\lambda>0$, regardless of whether
$\Sigma_{1,1}-2\Sigma_{1,2}+\Sigma_{2,2}=0$ or
$\Sigma_{1,1}-2\Sigma_{1,2}+\Sigma_{2,2}>0$. Finally, we note that the
condition $\Sigma_{1,1}-2\Sigma_{1,2}+\Sigma_{2,2}>0$ is testable.

######  Theorem 4.1.

Suppose Assumptions 1–3 hold. If $\Sigma_{1,1}-2\Sigma_{1,2}+\Sigma_{2,2}>0$
and $\hat{\Sigma}$ is an consistent estimator for $\Sigma$, then we have

$\displaystyle\sqrt{n}\left[(\hat{w},1-\hat{w})\hat{\Sigma}(\hat{w},1-\hat{w})^{\top}\right]^{-1/2}(\hat{\tau}^{*}-\tau)\rightsquigarrow\mathcal{N}(0,1),$

where $\hat{\tau}^{*}$ and $\hat{\omega}$ are defined in (4.1) and (4.2),
respectively. In this case, $\hat{\tau}^{*}$ is asymptotically weakly more
efficient than both $\hat{\tau}^{adj}$ and $\hat{\tau}^{unadj}$.

###  4.2 Local Power

Suppose we want to test the null hypothesis $H_{0}:\tau=\tau_{0}$ versus the
two-sided alternative $H_{1}:\tau\neq\tau_{0}$ with level $\alpha\in(0,1)$,
and we are under the local alternative in the sense that
$\tau=\tau_{0}+\Delta/\sqrt{n}$. Then, the hypothesis can be rewritten as
$H_{0}:\Delta=0$ versus $H_{1}:\Delta\neq 0$. In addition, by Theorem 3.1, we
have

$\displaystyle\sqrt{n}\Sigma^{-1/2}\begin{pmatrix}(\hat{\tau}^{adj}-\tau_{0})\\\
(\hat{\tau}^{unadj}-\tau_{0})\\\
\end{pmatrix}\rightsquigarrow\mathcal{N}\left(\begin{pmatrix}a\Delta\\\
b\Delta\end{pmatrix},I_{2}\right),$

where $\begin{pmatrix}a\\\ b\end{pmatrix}=\Sigma^{-1/2}\begin{pmatrix}1\\\
1\end{pmatrix}$. Let us consider a limit experiment in which a researcher
observes

$\displaystyle\begin{pmatrix}N_{1}\\\
N_{2}\end{pmatrix}\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}\left(\begin{pmatrix}a\Delta\\\
b\Delta\end{pmatrix},I_{2}\right),$ (4.3)

knows the values of $(a,b)$, and wants to test $\Delta=0$ versus the two-sided
alternative. Then, by the factorization criterion, it is easy to see that
$N^{*}=\frac{aN_{1}+bN_{2}}{\sqrt{a^{2}+b^{2}}}\sim
N(\sqrt{a^{2}+b^{2}}\Delta,1)$ is a sufficient statistic for $\Delta$. Define
the usual two-sided Wald test based on $\hat{\tau}^{*}$ as

$\displaystyle\mathbb{W}_{n}=1\left\\{n\left[(\hat{w},1-\hat{w})\hat{\Sigma}(\hat{w},1-\hat{w})^{\top}\right]^{-1}(\hat{\tau}^{*}-\tau_{0})^{2}\geq\mathcal{C}_{\alpha}\right\\},$
(4.4)

where $\mathcal{C}_{\alpha}$ is the $(1-\alpha)$ quantile of the chi-squared
distribution with one degree of freedom. Then, the following corollary implies
$\mathbb{W}_{n}\rightsquigarrow 1\\{(N^{*})^{2}\geq\mathcal{C}_{\alpha}\\}$.
In addition, $1\\{(N^{*})^{2}\geq\mathcal{C}_{\alpha}\\}$ is known as the
uniformly most powerful (UMP) unbiased test (see, for example, Lehmann and
Romano (2006, Section 4.2)) in the sense that it is more powerful than any
level $\alpha$ unbiased test that is a (potentially nonlinear) combination of
$(N_{1},N_{2})$. Consequently, it means $\mathbb{W}_{n}$ is asymptotically the
UMP test over a class of level $\alpha$ unbiased tests that only depend on

$\displaystyle\begin{pmatrix}\hat{N}_{1}\\\
\hat{N}_{2}\end{pmatrix}=\hat{\Sigma}^{-1/2}\begin{pmatrix}\sqrt{n}(\hat{\tau}^{adj}-\tau_{0})\\\
\sqrt{n}(\hat{\tau}^{unadj}-\tau_{0})\end{pmatrix}.$ (4.5)

This implies $\hat{\tau}^{*}$ is optimal among not only linear but also
nonlinear combination estimators.

To formalize the statement, we need to introduce the following notation.
Suppose $(N_{1},N_{2})$ follow the joint distribution in (4.3) and let

$\displaystyle\Psi_{\alpha}^{U}=\begin{Bmatrix}\psi(\cdot):\mathbb{E}\psi(N_{1},N_{2})\leq\alpha\quad\text{if
$\Delta=0$},\\\
\mathbb{E}\psi(\mathcal{N}_{1},\mathcal{N}_{2})\geq\alpha\quad\text{if
$\Delta\neq 0$},\\\ \text{the discontinuities of $\psi(\cdot)$ has zero
Lebesgue measure}\end{Bmatrix}$

be a class of unbiased tests with level $\alpha$.

######  Corollary 4.1.

Suppose all conditions in Theorem 4.1 hold, and we are under the local
alternative in the sense that $\tau=\tau_{0}+\Delta/\sqrt{n}$. Then, we have

$\displaystyle\sqrt{n}\left[(\hat{w},1-\hat{w})\hat{\Sigma}(\hat{w},1-\hat{w})^{\top}\right]^{-1/2}(\hat{\tau}^{*}-\tau_{0})\rightsquigarrow
N^{*}\stackrel{{\scriptstyle d}}{{=}}\mathcal{N}(\sqrt{a^{2}+b^{2}}\Delta,1),$

where $\begin{pmatrix}a\\\ b\end{pmatrix}=\Sigma^{-1/2}\begin{pmatrix}1\\\
1\end{pmatrix}$. In addition, suppose $\breve{\psi}_{n}$ is a generic test
such that $\breve{\psi}_{n}=\psi(\hat{N}_{1},\hat{N}_{2})+o_{P}(1)$ for some
$\psi\in\Psi_{\alpha}^{U}$ and the sequence $\\{\breve{\psi}_{n}\\}_{n\geq 1}$
is uniformly integrable, where $(\hat{N}_{1},\hat{N}_{2})$ are defined in
(4.5). Then, we have

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{E}\mathbb{W}_{n}=\sup_{\psi\in\Psi_{\alpha}^{U}}\lim_{n\rightarrow\infty}\mathbb{E}\psi(\hat{N}_{1},\hat{N}_{2})\geq\lim_{n\rightarrow\infty}\mathbb{E}\breve{\psi}_{n}.$

###  4.3 Covariance Matrix Estimator

We now propose a consistent estimator $\hat{\Sigma}$ for
$\Sigma=\Sigma_{\mathcal{U}}+\Sigma_{\mathcal{V}}+\Sigma_{\mathcal{W}}1_{2}1_{2}^{\top}$.
We first note that

$\Sigma_{2,2}=\mathbb{E}\left[\frac{var(Y_{i}(1)|S_{i})}{\pi_{S_{i}}}+\frac{var(Y_{i}(0)|S_{i})}{1-\pi_{S_{i}}}\right]+\Sigma_{\mathcal{W}}$

is the asymptotic variance of the unadjusted estimator and does not depend on
$X_{i}$. Now let

$\displaystyle\hat{\Sigma}_{2,2}=$
$\displaystyle\frac{1}{n}\sum_{s\in\mathcal{S}}\left[\sum_{i\in\aleph_{1,s}}\left(\frac{Y_{i}}{\hat{\pi}_{s}}-\frac{1}{n_{1,s}}\sum_{i\in\aleph_{1,s}}\frac{Y_{i}}{\hat{\pi}_{s}}\right)^{2}+\sum_{i\in\aleph_{0,s}}\left(\frac{Y_{i}}{1-\hat{\pi}_{s}}-\frac{1}{n_{0,s}}\sum_{i\in\aleph_{0,s}}\frac{Y_{i}}{1-\hat{\pi}_{s}}\right)^{2}\right]$
$\displaystyle+\sum_{s\in\mathcal{S}}\hat{p}_{s}\left(\frac{1}{n_{1,s}}\sum_{i\in\aleph_{1,s}}Y_{i}-\frac{1}{n_{0,s}}\sum_{i\in\aleph_{0,s}}Y_{i}-\hat{\tau}^{unadj}\right)^{2}.$

This estimator is consistent as long as Assumption 1 holds. See, for example,
Bugni et al. (2019) for details. Therefore, it suffices to construct
consistent estimators for the $(1,1)$ and $(1,2)$ elements of
$\Sigma_{\mathcal{U}}$ and $\Sigma_{\mathcal{V}}$, and $\Sigma_{\mathcal{W}}$.

To this end, denote $\hat{\beta}_{a,s}$ as the coefficient of $\breve{X}_{i}$
in the linear regression (3.1) and recall the intercept is $\hat{\tau}_{a,s}$.
The residual of the regression and its leave-one-out version are denoted as
$\hat{\varepsilon}_{a,s,i}=Y_{i}-\hat{\tau}_{a,s}-\breve{X}_{i}^{\prime}\hat{\beta}_{a,s}$
and $\acute{\varepsilon}_{a,s,i}=\hat{\varepsilon}_{a,s,i}/M_{a,s,i,i}$ for
$i\in\aleph_{a,s}$, respectively. Then, we construct our estimator
$\hat{\Sigma}_{\mathcal{U}}$ following the cross-fit approach developed by
Jochmans (2022). It is also possible to construct the estimator based on the
sample splitting approach proposed by Cattaneo et al. (2018). The cross-fit
estimator, as shown below, is asymptotically valid as long as
$\kappa_{a,s}<1$, but is not guaranteed to be positive definite in finite
sample. In contrast, the sample splitting estimator is guaranteed to be
positive semidefinite but requires $\kappa_{a,s}<1/2$. When both of them are
valid, they are asymptotically equivalent. Interested readers can refer to
Jochmans (2022) for more discussions. For the rest of the paper, we focus on
cross-fit estimators of the $(1,1)$ and $(1,2)$ elements of
$\hat{\Sigma}_{\mathcal{U}}$, which is defined as

$\displaystyle\hat{\Sigma}_{\mathcal{U}}=\begin{pmatrix}\sum_{s\in\mathcal{S}}\sum_{a=0,1}\frac{n_{s}^{2}}{nn_{a,s}}\hat{\omega}_{a,s}^{2}&\sum_{s\in\mathcal{S}}\sum_{a=0,1}\frac{n_{s}^{2}}{nn_{a,s}}\hat{\varpi}_{a,s}\\\
\sum_{s\in\mathcal{S}}\sum_{a=0,1}\frac{n_{s}^{2}}{nn_{a,s}}\hat{\varpi}_{a,s}&\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\end{pmatrix},$

where

$\displaystyle\hat{\omega}_{a,s}^{2}=\frac{1}{n_{a,s}}\gamma_{a,s,n}^{-2}\sum_{i\in\aleph_{a,s}}\left(\sum_{j\in\aleph_{a,s}}M_{a,s,i,j}\right)^{2}Y_{i}\acute{\varepsilon}_{a,s,i},$
$\displaystyle\hat{\varpi}_{a,s}=\frac{1}{n_{a,s}}\gamma_{a,s,n}^{-1}\sum_{i\in\aleph_{a,s}}\left(\sum_{j\in\aleph_{a,s}}M_{a,s,i,j}\right)Y_{i}\acute{\varepsilon}_{a,s,i}.$

To define $\hat{\Sigma}_{\mathcal{V}}$, we note that
$\mathbb{E}(Y_{i}(a)|X_{i},S_{i})-\mathbb{E}(Y_{i}(a)|S_{i})$ can be
approximated by $\breve{X}_{i}^{\top}\beta_{a,s}$ for $i\in\aleph_{a,s}$,
where $\beta_{a,s}$ is defined in Assumption 2(iv). Kline et al. (2020)
considered the estimation and inference for quadratic functions of the
coefficient in one linear regression with many regressors. However, in our
setting, $\Sigma_{\mathcal{V}}$ depends on multiple linear coefficients
$\beta_{a,s}$ which are estimated from different strata of observations, and
each regression model may only be approximately linear. Therefore, the
estimator $\hat{\Sigma}_{\mathcal{V}}$ includes both the quadratic and cross-
product terms of the estimated coefficients $\hat{\beta}_{a,s}$. Denote
$\Gamma_{a,s}=\sum_{i\in\aleph_{a,s}}\breve{X}_{i}\breve{X}_{i}^{\top}$ and
$\Gamma_{s}=\sum_{i\in\aleph_{s}}\breve{X}_{i}\breve{X}_{i}^{\top}$. Following
the spirit of Kline et al. (2020), our estimator of the $(1,1)$ and $(1,2)$
elements of $\hat{\Sigma}_{\mathcal{V}}$ is defined as

$\displaystyle\hat{\Sigma}_{\mathcal{V}}=\begin{pmatrix}\hat{\Sigma}_{\mathcal{V}}^{adj}&\hat{\Sigma}_{\mathcal{V}}^{adj}\\\
\hat{\Sigma}_{\mathcal{V}}^{adj}&\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\end{pmatrix},$

where

$\displaystyle\hat{\Sigma}_{\mathcal{V}}^{adj}=\sum_{s\in\mathcal{S}}\hat{p}_{s}\left[\sum_{a=0,1}\left(\frac{1}{n_{a,s}}\hat{\beta}_{a,s}^{\top}\Gamma_{a,s}\hat{\beta}_{a,s}-\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}P_{a,s,i,i}Y_{i}\acute{\varepsilon}_{a,s,i}\right)-\frac{2}{n_{s}}\hat{\beta}_{1,s}^{\top}\Gamma_{s}\hat{\beta}_{0,s}\right].$

The estimator of $\Sigma_{\mathcal{W}}$ is standard:

$\displaystyle\hat{\Sigma}_{\mathcal{W}}=\sum_{s\in\mathcal{S}}\hat{p}_{s}\left(\hat{\tau}_{1,s}-\tau_{0,s}-\hat{\tau}^{adj}\right)^{2}.$

######  Assumption 4.

Suppose
$\max_{a=0,1,s\in\mathcal{S}}\left\|\Gamma_{s}^{1/2}\Gamma_{a,s}^{-1}\Gamma_{s}^{1/2}\right\|_{op}=o_{P}(n).$

######  Remark 9.

By the random matrix theory, under general regularity conditions, one can show
that all the eigenvalues of $\Gamma_{s}/n_{s}$ and $\Gamma_{a,s}/n_{a,s}$ are
bounded and bounded away from zero, even if $\kappa_{a,s}>0$. This then
implies
$\max_{a=0,1,s\in\mathcal{S}}\left\|\Gamma_{s}^{1/2}\Gamma_{a,s}^{-1}\Gamma_{s}^{1/2}\right\|_{op}=O_{P}(1)$,
and thus, Assumption 4 holds.

######  Theorem 4.2.

Suppose Assumptions 1–4 hold. Let

$\hat{\Sigma}=\begin{pmatrix}\sum_{s\in\mathcal{S}}\sum_{a=0,1}\frac{n_{s}^{2}}{nn_{a,s}}\hat{\omega}_{a,s}^{2}+\hat{\Sigma}_{\mathcal{V}}^{adj}+\hat{\Sigma}_{\mathcal{W}}&\sum_{s\in\mathcal{S}}\sum_{a=0,1}\frac{n_{s}^{2}}{nn_{a,s}}\hat{\varpi}_{a,s}+\hat{\Sigma}_{\mathcal{V}}^{adj}+\hat{\Sigma}_{\mathcal{W}}\\\
\sum_{s\in\mathcal{S}}\sum_{a=0,1}\frac{n_{s}^{2}}{nn_{a,s}}\hat{\varpi}_{a,s}+\hat{\Sigma}_{\mathcal{V}}^{adj}+\hat{\Sigma}_{\mathcal{W}}&\hat{\Sigma}_{2,2}\end{pmatrix}.$

Then, we have

$\displaystyle\hat{\Sigma}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\Sigma.$

###  4.4 Further Discussions

Efficiency. Our inference procedure is ”optimal” within the class of unbiased
tests constructed by both adjusted and unadjusted estimators
($\hat{\tau}^{adj}$ and $\hat{\tau}^{unadj}$). If researchers have prior
knowledge about a subset of adjustment regressors and use them to form another
regression-adjusted estimator $\hat{\tau}^{prior}$, it is feasible to
optimally combine $(\hat{\tau}^{adj},\hat{\tau}^{unadj},\hat{\tau}^{prior})$
in the same manner. The crucial step involves constructing a consistent
estimator for the three-dimensional covariance matrix in the same approach
outlined in Section 4.3.

Theoretically, the approach can extend to optimally combine any fixed number
of regression-adjusted estimators. However, two problems arise. Firstly,
combining multiple estimators may make the covariance matrix nearly singular
and cause numerical issues. Secondly, using data-driven methods (e.g., Lasso)
to choose significant regressors may cause size distortion in the inference,
especially when regressors have weak but dense effects, as in Models 1-3 in
our simulation. Due to these problems, we focus on an inference procedure that
combines $\hat{\tau}^{adj}$ and $\hat{\tau}^{unadj}$ only.

Other Forms of RAs. Our linear RA method can capture non-linearity by
incorporating higher-order polynomials of baseline covariates. It can be shown
that, with the sieve bases of a fixed number of baseline covariates as
adjustment regressors, the adjusted estimator $\hat{\tau}^{adj}$ achieves the
semiparametric efficiency bound under the condition that $k_{n}$, the
dimension of regressors, is $o(\sqrt{n})$. This implies it is unnecessary to
consider other forms of RAs in the standard setup since there is no efficiency
gain. However, nonparametric sieve RAs often distort the size in finite
samples. In an asymptotic regime where $k_{n}$ is proportional to $n$, this
paper provides a theoretical explanation for this documented size distortion.
This result matters for empirical work too, since whether $k_{n}$ is
$o(\sqrt{n})$ or proportional to $n$ is unclear in practice. Through the
analysis in an alternative asymptotic regime, we provide a new inference
method that has better finite-sample performance.

Alternatively, we can use a single linear regression with all observations
instead of stratum-specific RA as in (3.1), which mitigates concerns related
to having many regressors. However, Liu et al. (2023) demonstrate that this RA
is less efficient than ours in the fixed dimension case. Consequently,
researchers face the challenge of choosing the appropriate asymptotic regime
for their model and, based on this choice, deciding whether to implement the
full-sample or stratum-specific RA. We aim to avoid this choice. In Section 5
below, we demonstrate that our inference method remains valid and is the most
efficient among linear RAs in an alternative asymptotic regime with a fixed
number of regressors. This implies that researchers can apply our method
without the need to choose the asymptotic regime.

Indeed, stratum-specific RAs reduce the effective sample size, exacerbating
the dimensionality issue. When $k_{n}$ exceeds the effective sample size
$n_{a,s}$, researchers can employ Ridge-regularized linear regression for
adjustment. Studying the statistical properties of the corresponding
regression-adjusted ATE estimator is an interesting avenue for future
research.

Other Randomization Schemes. We focus on CARs, which are widely used in
empirical research. See Chong et al. (2016), Jakiela and Ozier (2016),
Burchardi et al. (2019), Anderson and McKenzie (2021), etc. Theoretically, Bai
(2022), Cytrynbaum (2023), and Bai et al. (2023) highlight the efficiency of
“finely stratified” experiments with fixed-dimension baseline covariates. But
“fine stratification” may have increasing covariate differences within strata
as the dimension of the matching variable grows, potentially reducing
efficiency. This creates a similar trade-off in “fine stratification” based on
a subset of covariates in the experiment design stage, a topic deserving
further investigation beyond the scope of this paper.

Unlike Bai (2022) and Tabord-Meehan (2023), who aim to determine the “optimal”
experiment in the design stage based on a pilot experiment, our approach,
aligned with the literature on RAs, takes the randomization scheme (i.e.,
CARs) as given and seeks more efficient estimators with covariate information.
Our proposed method does not require a pilot experiment and is particularly
well-suited for the analysis of experiments conducted in the past or by
others, thus complementing their methods.

##  5 Fixed Number of Regressors

In this section, we consider the CARs specified in Assumption 1 with a fixed
dimension of $X_{i}$, i.e., $k$ is fixed. We then introduce another set of
regularity conditions to replace Assumptions 2–4. Under these new conditions,
the Wald test based on the same estimator $\hat{\tau}^{*}$ \- constructed
using the same $\hat{\tau}^{adj}$, $\hat{\tau}^{unadj}$, and $\hat{\Sigma}$ as
defined above - maintains exact asymptotic size under the null and is weakly
more powerful than the Wald test based on the unadjusted estimator
$\hat{\tau}^{unadj}$.

Particularly, in this fixed dimension regime, there’s no need to assume that
the RAs are approximately correctly specified, as in Assumption 2(iv). So the
results here are not merely special cases of those in Section 3. Additionally,
they are also novel to the literature of RAs with fixed number of regressors
(e.g., Ye et al. (2022)). This is because the cross-fit estimator of the
covariance matrix and the optimal linear combination estimator, motivated by
our high-dimensional analysis, are new. In Section B of the Online Supplement,
we establish the uniform size control of our inference method across
asymptotic regimes with fixed or diverging (up to the sample fraction) numbers
of regressors. This means that researchers can use our method without choosing
the asymptotic regime.

######  Assumption 5.

Denote the dimension of $X_{i}$ as $k$ which is fixed w.r.t. the sample size.
In addition, there exist constants $c,C$ such that

$\displaystyle
0<c\leq\lambda_{\min}\left(\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}\breve{X}_{i}\breve{X}_{i}^{\top}\right)\leq\lambda_{\max}\left(\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}\breve{X}_{i}\breve{X}_{i}^{\top}\right)\leq
C<\infty$

and $\max_{i\in[n]}||\breve{X}_{i}||_{2}=o_{P}(n^{1/2})$.

######  Remark 10.

Assumption 5 is standard when the dimension of $\breve{X}_{i}$ is fixed.

######  Theorem 5.1.

Suppose Assumptions 1 and 5 holds. Then, we have

$\displaystyle\sqrt{n}\begin{pmatrix}\hat{\tau}^{adj}-\tau\\\
\hat{\tau}^{unadj}-\tau\end{pmatrix}\rightsquigarrow\mathcal{N}\left(\begin{pmatrix}0\\\
0\end{pmatrix},\Omega\right),$ $\displaystyle\Omega$
$\displaystyle=\left\\{\mathbb{E}\left[\frac{Var(Y_{i}(1)|X_{i},S_{i})}{\pi_{S_{i}}}+\frac{Var(Y_{i}(0)|X_{i},S_{i})}{1-\pi_{S_{i}}}\right]+\mathbb{E}(m_{1}(X_{i},S_{i})-m_{0}(X_{i},S_{i})-\tau)^{2}\right\\}1_{2}1_{2}^{\top}$
$\displaystyle+\sum_{s\in\mathcal{S}}\frac{p_{s}}{\pi_{s}(1-\pi_{s})}\begin{pmatrix}V_{s}&V_{s}\\\
V_{s}&V_{s}^{\prime}\end{pmatrix},$

and

$\displaystyle\hat{\Sigma}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\Omega,$

where $m_{a}(x,s)=\mathbb{E}(Y_{i}(s)|X_{i}=x,S_{i}=s)$,
$\overline{\beta}_{s}^{*}=(1-\pi_{s})\beta_{1,s}^{*}+\pi_{s}\beta_{0,s}^{*}$,

$\displaystyle\beta_{a,s}^{*}=Var(X_{i}\mid
S_{i}=s)^{-1}Cov(X_{i},Y_{i}(a)|S_{i}=s),$ $\displaystyle
V_{s}=Var((1-\pi_{s})m_{1}(X_{i},s)+\pi_{s}m_{0}(X_{i},s)-X_{i}^{\top}\overline{\beta}_{s}^{*}\mid
S_{i}=s),$ $\displaystyle
V_{s}^{\prime}=Var((1-\pi_{s})m_{1}(X_{i},s)+\pi_{s}m_{0}(X_{i},s)\mid
S_{i}=s),$

If we further assume $\pi_{s}\beta_{0,s}^{*}+(1-\pi_{s})\beta_{1,s}\neq 0$,
then

$\displaystyle\sqrt{n}\left[(\hat{w},1-\hat{w})\hat{\Sigma}(\hat{w},1-\hat{w})^{\top}\right]^{-1/2}(\hat{\tau}^{*}-\tau)\rightsquigarrow\mathcal{N}(0,1).$

In this case, $\hat{\tau}^{*}$ is asymptotically equivalent to
$\hat{\tau}^{adj}$ in sense that
$\hat{\tau}^{*}=\hat{\tau}^{adj}+o_{P}(n^{-1/2})$ and weakly more efficient
than $\hat{\tau}^{unadj}$.

######  Remark 11.

The limit distribution of $\hat{\tau}^{adj}$ has already been derived by Ye et
al. (2022). Here, our main contributions are (1) deriving the joint
distribution of $\hat{\tau}^{adj}$ and $\hat{\tau}^{unadj}$ and (2)
establishing the consistency of our cross-fit covariance matrix estimator
$\hat{\Sigma}$ under the fixed dimension framework. Theorem 5.1 shows that the
Wald test based on our optimal linear combination estimator $\hat{\tau}^{*}$
and the covariance matrix estimator $\hat{\Sigma}$ proposed in Section 4.3
still controls asymptotic size under the null and has the same power as
$\hat{\tau}^{adj}$ under local alternatives with a fixed dimension of $X_{i}$.
In fact, when the dimension of $X_{i}$ is fixed, $\hat{\tau}^{*}$ is
equivalent to $\hat{\tau}^{adj}$, which is the optimally linearly adjusted
estimator in the sense that it achieves the minimum asymptotic variance among
a class of estimators adjusted by linear functions of $X_{i}$. Therefore, it
is weakly more efficient that $\hat{\tau}^{unadj}$. See, for example, Liu et
al. (2023), Ma et al. (2022), and Ye et al. (2022) for more discussion.

##  6 Simulations

We conduct a simulation study to evaluate the finite sample performance of our
proposed inference method. The nominal level is $\alpha=5\%$ throughout the
section.

###  6.1 Data Generating Processes

For $a\in\\{0,1\\}$, we generate potential outcomes according to the equation

$Y_{i}(a)=\mu_{a}+m_{a}(Z_{i})+\sigma_{a}(Z_{i})\varepsilon_{i}(a),$ (6.1)

where $\mu_{a}$, and $m_{a}\left(Z_{i}\right)$ and $\sigma_{a}(Z_{i})$ are
specified as follows. In each of the following specifications,
$\\{Z_{i},\varepsilon_{i}(1),\varepsilon_{i}(0)\\}$ are i.i.d. We have
considered four models following the designs by Cattaneo et al. (2018). For
the first three models, the dimension of $Z_{i}$ is denoted as $d_{n}$, and we
let $X_{i}$ be either the full set or the first $k_{n}$ elements of $Z_{i}$.
For the last model, $Z_{i}$ is of fixed dimensional, and we let $X_{i}$ be
either the full set or the first $k_{n}$ elements of power series of $Z_{i}$.
The power series will be specified below. By varying the choice of $k_{n}$, we
can illustrate the uniformity of our inference method.

Model 1: Linear model with many dummy variables.

The dimension of $Z_{i}$ is set to $d_{n}=0.2n/|\mathcal{S}|$. The first entry
of $Z_{i}$ is uniformly distributed in $[-1,1]$, i.e., $Z_{1i}\sim U[-1,1]$.
The other entries of $Z_{i}$ are dummies:
$[Z_{2i},Z_{3i},\cdots,Z_{d_{n},i}]^{\top}=\textbf{1}(\textbf{v}_{i}\geq\Phi^{-1}(0.8))$
with $\textbf{v}_{i}\sim\mathcal{N}(0,\textbf{I}_{d_{n}-1})$. We set
$m_{a}(\cdot)$ and $\sigma_{a}(\cdot)$ respectively as

$m_{1}(Z_{i})=m_{0}(Z_{i})=Z_{1i}+2\sum_{j=2}^{d_{n}}Z_{ji}/\sqrt{d_{n}-1}$

and

$\sigma_{0}(Z_{i})=\sigma_{1}(Z_{i})=c_{\varepsilon}\left[1+\left(Z_{1i}+\sum_{j=2}^{d_{n}}Z_{ji}/\sqrt{d_{n}-1}\right)^{2}\right]^{1/2}.$

Let $(\varepsilon_{i}(1),\varepsilon_{i}(0))\sim\mathcal{N}(0,\mathbf{I}_{2})$
be independent of $Z_{i}$. We set $c_{\varepsilon}$ as the normalizing
constant such that $\mathbb{E}\sigma_{a}^{2}(Z_{i})=1$.

Model 2: Linear model with many independent continuous variables.

This model is the same as model 1, except that $Z_{ji}$’s are $U[-1,1]$ and

$m_{1}(Z_{i})=m_{0}(Z_{i})=2\sum_{j=2}^{d_{n}}Z_{ji}/\sqrt{d_{n}-1},$

where $Z_{1i}$ does not appear in $m_{1}(Z_{i})$ and $m_{0}(Z_{i})$.

Model 3: Linear model with many correlated continuous variables.

This model is the same as Model 2, except that
$[Z_{2i},\cdots,Z_{d_{n},i}]^{\top}=\Sigma_{\rho}^{1/2}\mathbf{v}_{i}$, where
$\mathbf{v}_{i}\in\mathbb{R}^{d_{n}-1}$ and $\Sigma_{\rho}$ is a Toeplitz
matrix with $[\Sigma_{\rho}]_{i,j}=\rho^{|i-j|}$ and $\rho=0.6$. We set
$\mathbf{v}_{i}$ to contain i.i.d. $U[-1,1]$ entries.

Model 4: Nonlinear model with many power series.

Let $d_{n}=6$, $Z_{ji}\sim U[-1,1]$ for $j=1,\cdots,6$,

$m_{1}(Z_{i})=m_{0}(Z_{i})=2\exp\left(\sqrt{\frac{1}{6}\sum_{j=1}^{6}Z_{ji}^{2}}\right),$

and

$\sigma_{0}(Z_{i})=\sigma_{1}(Z_{i})=c_{\varepsilon}\left[1+\left(Z_{1i}+\sum_{j=2}^{6}Z_{ji}/\sqrt{5}\right)^{2}\right]^{1/2}.$

We use the following polynomials or some of them as regressors:

  * •

First order terms: $Z_{1i},Z_{2i},Z_{3i},Z_{4i},Z_{5i},Z_{6i}$,

  * •

Second order terms:
$Z_{1i}^{2},Z_{2i}^{2},Z_{3i}^{2},Z_{4i}^{2},Z_{5i}^{2},Z_{6i}^{2}$ and first-
order interactions,

  * •

Third-order terms:
$Z_{1i}^{3},Z_{2i}^{3},Z_{3i}^{3},Z_{4i}^{3},Z_{5i}^{3},Z_{6i}^{3}$ and
second-order interactions,

  * •

Fourth-order terms:
$Z_{1i}^{4},Z_{2i}^{4},Z_{3i}^{4},Z_{4i}^{4},Z_{5i}^{4},Z_{6i}^{4}$ and third-
order interactions.

For each model, strata are determined by dividing the support of $Z_{1i}$ into
$|\mathcal{S}|$ intervals of equal length. Specifically, we have
$S_{i}=\sum_{j=1}^{|\mathcal{S}|}1\\{Z_{1i}\leq g_{j}\\},$ where
$g_{j}=2j/|\mathcal{S}|-1$. For all strata, we set $\pi_{s}=1/2$. The
treatment status is determined according to one of the following four CAR
schemes:

  1. 1.

SRS: Treatment assignment is generated as in Example 1,

  2. 2.

BCD: Treatment assignment is generated as in Example 2 with $\lambda=0.75$,

  3. 3.

WEI: Treatment assignment is generated as in Example 3 with $\phi(x)=(1-x)/2$,

  4. 4.

SBR: Treatment assignment is generated as in Example 4.

###  6.2 Simulation Results

Rejection probabilities are computed using $10,000$ replications. We compare
our methods with three existing methods recently introduced by Bugni et al.
(2018), Ye et al. (2022), and Liu et al. (2023). We list these methods below.

$\hat{\tau}^{adj}$:

The inference method based on the fully saturated regression adjusted
estimator introduced in Section 2 and the variance estimator
$\hat{\Sigma}_{1,1}$ as in Theorem 4.2.

$\hat{\tau}^{\ast}$ :

The inference method based on the optimal linear combination estimator
introduced in Section 4 and Theorem 4.1.

BCS:

The inference method introduced by Bugni et al. (2019). The estimator for
$\tau$ is exactly the unadjusted IPW estimator $\hat{\tau}^{unadj}$.

YYS:

The inference method introduced by Ye et al. (2022). The estimator for $\tau$
is exactly the same as $\hat{\tau}^{adj}$. However, its variance estimator
does not take into account the issue of many covariates.

LTM:

The inference method introduced by Liu et al. (2023), which uses Lasso to
select relevant regressors in RA. Specifically, we use the ‘HDM’ package in R
with the default choice of tuning parameter to obtain the corresponding Lasso
estimates.

Tables 1 reports the size results for $n=400$ and $800$, under
$H_{0}:\mu_{1}=\mu_{0}=0$ and $S=2$. In the first three models, we include all
elements of $Z_{i}$ as regressors, i.e., $X_{i}=Z_{i}$. For the last model, we
include all power series of $Z_{i}$ as regressors. Therefore, the number of
covariates used in each regression is approximately $\mathbf{40\%}$ of the
effective sample size (i.e., $\kappa_{a,s}=0.4$). We find that our methods
$\hat{\tau}^{adj}$ and $\hat{\tau}^{\ast}$ have good size control under four
models, improving as $n$ increases. In contrast, YYS has high rejection rates,
over 10%, in all cases, even when $n$ is as large as 800. Our theory in
Section 3 shows this size distortion is caused by ignoring the estimation
errors from many regressors. Furthermore, LTM has rejection probabilities over
7% in many cases. In particular, in Model 3 with correlated covariates, LTM
has rejection rates over 9% in all randomization schemes, even when $n$ is
800. This is because in this model, the sparsity condition does not hold.
Finally, BCS does not use any regressors, and thus, is immune to the size
distortion due to the many regressors.

Table 2 shows the power results under $H_{1}:\mu_{1}-\mu_{0}=0.2$. We find
that, first, $\hat{\tau}^{\ast}$ always has higher power than
$\hat{\tau}^{adj}$ and BCS, as our theory predicts. Second, YYS ignores the
cost of RAs, resulting in a smaller variance estimator, and thus, higher
rejection rates under both the null and alternative. Third, LTM has lower
power than $\hat{\tau}^{\ast}$ in Models 1-3, where the sparsity condition for
the regression coefficients fails. In fact, in this case, Lasso excludes many
informative covariates in the adjustments. In model 4, LTM has better power
than $\hat{\tau}^{\ast}$ since the sparsity condition is satisfied. By
exploiting the sparsity, it is possible to show that the ATE estimator with
the $\ell_{1}$ regularized adjustments achieves the semiparametric efficiency
bound. However, even in this case that favors the Lasso regression
adjustments, the power of our method based on $\hat{\tau}^{*}$ is very close
to LTM.

Table 1: Rejection rate (in percent) under $H_{0}:\mu_{1}-\mu_{0}=0$ |  | n=400, $k_{n}=40$ | n=800, $k_{n}=80$
---|---|---|---
Model |  Method |  SRS |  BCD |  WEI |  SBR | SRS | BCD | WEI | SBR
1 |  $\hat{\tau}^{adj}$ |  5.78 |  5.89 |  5.67 |  5.89 |  5.43 |  5.32 |  5.46 |  5.32
|  $\hat{\tau}^{\ast}$ |  5.72 |  5.96 |  5.81 |  5.96 |  5.37 |  5.41 |  5.28 |  5.41
|  BCS |  5.05 |  5.41 |  5.22 |  5.41 |  5.12 |  5.30 |  5.18 |  5.30
|  YYS |  12.04 |  12.23 |  12.56 |  12.23 |  11.73 |  11.53 |  11.53 |  11.53
|  LTM |  6.64 |  7.08 |  7.16 |  7.08 |  5.95 |  6.33 |  6.07 |  6.33
2 |  $\hat{\tau}^{adj}$ |  5.51 |  5.57 |  5.32 |  5.57 |  5.10 |  5.65 |  5.38 |  5.65
|  $\hat{\tau}^{\ast}$ |  5.42 |  5.53 |  5.18 |  5.53 |  5.06 |  5.63 |  5.51 |  5.63
|  BCS |  5.04 |  4.88 |  5.12 |  4.88 |  4.82 |  5.18 |  5.14 |  5.18
|  YYS |  11.53 |  11.49 |  11.53 |  11.49 |  11.56 |  11.63 |  11.35 |  11.63
|  LTM |  7.36 |  7.33 |  7.29 |  7.33 |  7.25 |  7.43 |  7.02 |  7.43
3 |  $\hat{\tau}^{adj}$ |  5.90 |  6.12 |  5.79 |  5.85 |  5.44 |  5.75 |  5.48 |  5.75
|  $\hat{\tau}^{\ast}$ |  5.76 |  6.23 |  5.83 |  5.98 |  5.39 |  5.75 |  5.60 |  5.75
|  BCS |  4.79 |  4.84 |  5.15 |  5.27 |  5.12 |  5.04 |  5.23 |  5.04
|  YYS |  11.47 |  11.32 |  11.33 |  11.26 |  11.38 |  11.68 |  11.3 |  11.68
|  LTM |  9.82 |  9.47 |  9.22 |  9.53 |  9.69 |  9.69 |  9.62 |  9.69
4 |  $\hat{\tau}^{adj}$ |  6.63 |  7.11 |  7.25 |  7.09 |  6.28 |  6.45 |  6.68 |  6.45
|  $\hat{\tau}^{\ast}$ |  7.21 |  7.57 |  7.57 |  7.64 |  6.29 |  6.67 |  6.52 |  6.67
|  BCS |  4.83 |  5.02 |  5.09 |  5.37 |  4.93 |  5.28 |  4.93 |  5.28
|  YYS |  12.52 |  12.87 |  13.40 |  12.62 |  12.46 |  12.86 |  13.08 |  12.86
|  LTM |  6.07 |  6.24 |  6.14 |  6.57 |  5.76 |  5.87 |  5.92 |  5.87

Figures 2 and 3 show how the rejection rates of $\hat{\tau}^{adj}$,
$\hat{\tau}^{\ast}$, YYS and LTM vary with the number of regressors under the
null and alternative hypotheses. We consider Models 1-4 with 80 and 40
regressors for $n=800$ and $n=400$, respectively. We set the number of
regressors $k_{n}=0,2,\cdots,80$ ($k_{n}=0,1,\cdots,40)$ for $n=800$
($n=400$), corresponding to $\kappa_{a,s}=0,0.02,\cdots,0.4$. The x-axis is
$\kappa_{a,s}$. The y-axis is the rejection probabilities. We display only the
figures under the SBR scheme. Similar patterns are found under the other
randomization schemes, which are shown in the Online Supplement.

Panels (a) of the figures show that our methods (i.e., $\hat{\tau}^{adj}$ and
$\hat{\tau}^{*}$) have uniform size control over $\kappa_{a,s}$. Even when the
dimensions of regressors vary from 0 to 80 (40) for $n=800$ ($n=400$), our
methods have close-to-nominal rejection rates under $H_{0}$, consistent with
Theorem B.1. In contrast, YYS has rejection probabilities that increase
linearly with $\kappa_{a.s}$ in all models. Similarly, LTM displays this size
distortion in Models 2 and 3, in which the sparsity condition fails.

Table 2: Rejection rate (in percent) under $H_{1}:\mu_{1}-\mu_{0}=0.2$ |  | n=400, $k_{n}=40$ | n=800, $k_{n}=80$
---|---|---|---
Model |  Method | SRS | BCD | WEI | SBR | SRS | BCD | WEI | SBR
1 |  $\hat{\tau}^{adj}$ |  41.54 |  41.60 |  42.16 |  41.34 |  68.42 |  68.99 |  68.12 |  68.99
|  $\hat{\tau}^{\ast}$ |  44.66 |  44.93 |  45.37 |  44.44 |  72.31 |  72.44 |  72.15 |  72.44
|  BCS |  32.53 |  33.67 |  33.17 |  32.68 |  57.23 |  57.86 |  57.59 |  57.86
|  YYS |  55.59 |  55.82 |  56.66 |  55.33 |  80.26 |  80.62 |  80.57 |  80.62
|  LTM |  40.68 |  41.23 |  41.07 |  40.20 |  65.89 |  65.79 |  65.51 |  65.79
2 |  $\hat{\tau}^{adj}$ |  41.88 |  41.99 |  41.99 |  42.28 |  68.33 |  68.82 |  68.64 |  68.82
|  $\hat{\tau}^{\ast}$ |  43.70 |  43.80 |  44.05 |  43.87 |  70.56 |  71.22 |  70.69 |  71.22
|  BCS |  25.66 |  25.50 |  25.71 |  25.33 |  44.38 |  45.34 |  44.88 |  45.34
|  YYS |  56.10 |  55.84 |  55.35 |  55.93 |  79.83 |  80.32 |  80.66 |  80.32
|  LTM |  35.74 |  35.58 |  34.97 |  35.44 |  55.46 |  56.22 |  56.00 |  56.22
3 |  $\hat{\tau}^{adj}$ |  42.56 |  42.14 |  42.39 |  42.59 |  68.24 |  68.30 |  68.34 |  68.30
|  $\hat{\tau}^{\ast}$ |  43.13 |  43.05 |  42.99 |  42.77 |  68.60 |  68.98 |  69.04 |  68.98
|  BCS |  12.74 |  12.89 |  13.01 |  12.83 |  20.00 |  20.91 |  20.10 |  20.91
|  YYS |  55.62 |  55.80 |  55.43 |  55.63 |  79.62 |  79.80 |  80.05 |  79.80
|  LTM |  43.79 |  43.41 |  44.14 |  43.87 |  65.27 |  66.16 |  66.14 |  66.16
4 |  $\hat{\tau}^{adj}$ |  41.02 |  42.41 |  41.59 |  41.52 |  65.65 |  65.67 |  65.11 |  65.55
|  $\hat{\tau}^{\ast}$ |  45.59 |  46.78 |  45.88 |  46.14 |  72.20 |  71.57 |  71.58 |  72.21
|  BCS |  33.66 |  34.90 |  34.13 |  34.35 |  59.12 |  59.06 |  59.20 |  58.55
|  YYS |  55.25 |  56.03 |  55.64 |  55.51 |  78.73 |  78.85 |  78.31 |  78.58
|  LTM |  47.38 |  48.09 |  47.42 |  47.47 |  78.35 |  77.38 |  77.86 |  77.59

(a) $H_{0}:\mu_{1}-\mu_{0}=0$

(b) $H_{1}:\mu_{1}-\mu_{0}=0.2$

Figure 2: Rejection probabilities under SBR when $n=400$

(a) $H_{0}:\mu_{1}-\mu_{0}=0$

(b) $H_{1}:\mu_{1}-\mu_{0}=0.2$

Figure 3: Rejection probabilities under SBR when $n=800$

Panels (b) of the figures show the power of the five methods under
$H_{1}:\mu_{1}-\mu_{0}=0.2$. Several observations emerge. First, BCS’s power,
our benchmark, does not vary with $\kappa_{a,s}$ because it uses no
regressors. Second, for all models, $\hat{\tau}^{\ast}$ is always more
powerful than BCS and usually more powerful than $\hat{\tau}^{adj}$, as our
theory predicts. Third, as the number of regressors increases, the
misspecification error decreases, and $\hat{\tau}^{\ast}$ becomes more
powerful. However, $\hat{\tau}^{adj}$’s power drops even below that of the
unadjusted estimator as $k_{n}$ increases for Model 1. This reflects that many
covariates do not always improve estimation accuracy, and the “no-harm”
regression adjustments can in fact harm the estimation precision. Fourth, for
Models 1-3, $\hat{\tau}^{\ast}$ is generally more powerful than LTM as the
dimension of regressors increases. For Model 4, LTM slightly outperforms
$\hat{\tau}^{\ast}$ when $\kappa_{a,s}$ is large for $n=800$ due to the
sparsity of the model.

##  7 Empirical Application

How do moral incentives affect economic decisions? Using a field experiment
with the credit card customers of a large bank in Indonesia, Bursztyn et al.
(2019) studied the role of moral text messages in an important economic
decision: the decision to repay one’s debts. In the main treatment, clients
received text messages emphasizing the ethical consequences of non-repayment.
They also conducted a subsequent experiment with the clients who had received
the moral message in the main experiment but were still late in their
payments. Bursztyn et al. (2019) estimated the average treatment effects
(ATEs) of repeated moral appeals on debt repayment.

We focus on this follow-up experiment, which involved 898 recurrent late
payers who were randomly assigned to either the treatment or control group.
The randomization was conducted under CARs. Specifically, late payers were
first stratified by the time elapsed since they received the first moral
message. Two strata were formed, one with 364 late payers and the other with
534 payers. Within each stratum, half of the payers were treated and the rest
were in the control. It is straightforward to verify that Assumption 1(i),
(ii), and (iii) are satisfied. Assumption 1(iv) is also satisfied because
$\max_{s\in\mathcal{S}}\left|\frac{D_{n,s}}{n_{s}}\right|\approx 0.002$.

We consider three main outcome variables in our analysis: delinquency, which
is a binary indicator of whether customers repaid before the end of the grace
period; log repayment, which is the logarithm of the amount repaid; and
repayment ratio, which is the ratio of the amount repaid to the minimum
repayment required by the bank. We follow Bursztyn et al. (2019) to estimate
the same ATEs on delinquency and report the additional results for log
repayment and repayment ratio.

To estimate the treatment effects on debt repayment, Bursztyn et al. (2019)
used regression models with a large set of covariates, including client and
account characteristics and wave fixed effects. The dataset contains 4
continuous variables and 35 dummy variables. We set $X_{i}$ to include all
these covariates and some interaction terms.888Denote
$X_{1},X_{2},X_{3},X_{4}$ as the four continuous variables and $B_{j}$ as the
cubic spline basis with three knots for continuous variable $j$. The
interaction terms include $\\{X_{1}\times X_{2},X_{1}\times X_{3},X_{1}\times
X_{4},X_{2}\times X_{3},X_{2}\times X_{4},X_{3}\times
X_{4},(X_{2},X_{3},X_{4})\times B_{1},(X_{1},X_{3},X_{4})\times
B_{2},(X_{1},X_{2},X_{4})\times B_{3},(X_{1},X_{2},X_{3})\times B_{4}\\}$ The
dimension of $X_{i}$ is therefore 57, $\kappa_{1,1}$ and $\kappa_{0,1}$ are
about 0.31, and $\kappa_{1,2}$ and $\kappa_{0,2}$ are about 0.21.999The first
stratum has 182 payers in each group. The second stratum has 268 and 266
payers in the treatment and control, respectively.

Table 3 presents the ATEs and the standard errors (in parentheses) estimated
by different methods.101010The description of these methods is similar to that
in Section 6.2. To maintain comparability, we use the same $X_{i}$ for all
three outcome variables and across all the adjusted estimators.  We make two
observations from the results. First, for all three outcome variables, the
variance estimates of the YYS method are smaller than those of
$\hat{\tau}^{adj}$. Their difference is significant in some cases. For
example, for repayment ratio, although the YYS estimates are identical to
$\hat{\tau}^{adj}$, their standard errors are only about 66% of the respective
values of $\hat{\tau}^{adj}$. This indicates that YYS method does not account
for the estimation errors arising from using many regressors.

Second, for all three outcome variables, the standard errors of
$\hat{\tau}^{\ast}$ are smaller than those of both $\hat{\tau}^{adj}$ and BCS.
This is expected because $\hat{\tau}^{\ast}$ is a weighted average of
$\hat{\tau}^{adj}$ and the unadjusted estimator (BCS), where the weight is
chosen to minimize the asymptotic variance. The standard errors of
$\hat{\tau}^{\ast}$ are close to those of LTM. In particular, when the outcome
variable is the repayment ratio, the standard error of $\hat{\tau}^{\ast}$ is
about 5.9% lower than that of LTM.

Table 3: Impacts of Moral Incentives on Credit Card Debt Repayment

Y |  $\hat{\tau}^{adj}$ |  $\hat{\tau}^{\ast}$ |  BCS |  YYS |  LTM
---|---|---|---|---|---
Delinquency |  -0.034 |  -0.039 |  -0.041 |  -0.034 |  -0.040
|  (0.032) |  (0.030) |  (0.031) |  (0.031) |  (0.030)
Log Repayment |  0.441 |  0.528 |  0.578 |  0.441 |  0.541
|  (0.438) |  (0.416) |  (0.423) |  (0.424) |  (0.415)
Repayment Ratio |  0.030 |  0.065 |  0.054 |  0.030 |  0.055
|  (0.326) |  (0.209) |  (0.222) |  (0.216) |  (0.221)

Notes: The table reports the ATE estimates of the effect of moral incentives
on credit card debt repayment. Standard errors are in parentheses.

##  8 Conclusion

In this paper, we show that the estimation error of the “no-harm” regression
adjustments can contaminate the ATE estimator and degrade estimation
efficiency when the number of regressors is of the same order as the sample
size. We then propose a new ATE estimator which is guaranteed to be weakly
more efficient than the unadjusted estimator. We also propose a consistent
estimator of its asymptotic variance, construct the corresponding Wald-
statistic, and show its uniform size control no matter whether the dimension
of the regressors is fixed or grows at a rate no faster than the sample size.

##  Appendix A Verifying Assumption 3 for Gaussian Covariates

For $i\in\aleph_{s}$, we assume that $X_{i}=\Sigma_{s,n}^{1/2}Z_{i}$ where
$Z_{i}=(Z_{i,1},\cdots,Z_{i,k_{n}})$ are $k_{n}$ independent standard normal
random variables and the covariance matrix $\Sigma_{s,n}\in\Re^{k_{n}\times
k_{n}}$ is symmetric and positive definite. The assumption that $Z_{i}$ has
zero mean is without loss of generality because the means will be canceled in
the definition of $\breve{X}_{i}$. We focus on the case with $a=1$. By an
abuse of notation, we assume $\aleph_{1,s}=\\{1,\cdots,n_{1,s}\\}$ and
$\aleph_{0,s}=\\{1+n_{1,s},\cdots,n_{s}\\}$ and denote
$Z_{\aleph_{1,s}}=(Z_{1},\cdots,Z_{n_{1,s}})^{\top}\in\Re^{n_{1,s}\times
k_{n}}$ and $Z_{\aleph_{s}}=(Z_{1},\cdots,Z_{n_{s}})^{\top}\in\Re^{n_{s}\times
k_{n}}$. Then, we have

$\breve{X}_{\aleph_{1,s}}=X_{\aleph_{1,s}}-1_{n_{1,s}}1_{n_{s}}^{\top}X_{\aleph_{s}}/n_{s}\stackrel{{\scriptstyle
d}}{{=}}(Z_{\aleph_{1,s}}-1_{n_{a,s}}1_{n_{s}}^{\top}Z_{\aleph_{s}}/n_{s})\Sigma_{s,n}^{1/2}=\Gamma_{n}Z_{\aleph_{s}}\Sigma_{s,n}^{1/2}$

and

$\displaystyle
P_{1,s}=\breve{X}_{\aleph_{1,s}}\left[\breve{X}_{\aleph_{1,s}}^{\top}\breve{X}_{\aleph_{1,s}}\right]^{-1}\breve{X}_{\aleph_{1,s}}^{\top}\stackrel{{\scriptstyle
d}}{{=}}\Gamma_{n}Z_{\aleph_{s}}\left[Z_{\aleph_{s}}^{\top}\Gamma_{n}^{\top}\Gamma_{n}Z_{\aleph_{s}}\right]^{-1}Z_{\aleph_{s}}^{\top}\Gamma_{n}^{\top},$

where
$\Gamma_{n}=(I_{n_{1,s}},0_{n_{1,s},n_{0,s}})-1_{n_{1,s}}1_{n_{s}}^{\top}/n_{s}$
and $0_{n_{1,s},n_{0,s}}$ is a $n_{1,s}\times n_{0,s}$ matrix of zeros. For
the $l$-th column of $\Gamma_{n}Z_{\aleph_{s}}$ denoted as
$\Gamma_{n}Z_{\aleph_{s},l}$, we have

$\displaystyle\Gamma_{n}Z_{\aleph_{s},l}\stackrel{{\scriptstyle
d}}{{=}}W_{n}^{1/2}Z_{\aleph_{1,s},l},$

where $Z_{\aleph_{1,s},l}$ is an $n_{1,s}\times 1$ standard normal vector and
$W_{n}$ is a $n_{1,s}\times n_{1,s}$ matrix with $1-1/n_{s}$ and $-1/n_{s}$
being the diagonal and off-diagonal elements, respectively. Given that
$Z_{\aleph_{s}}$ are independent across columns, we have

$\displaystyle\Gamma_{n}Z_{\aleph_{s}}\stackrel{{\scriptstyle
d}}{{=}}W_{n}^{1/2}Z_{\aleph_{1,s}},$

which implies

$\displaystyle 1-\gamma_{1,s,n}$ $\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\frac{1_{n_{1,s}}^{\top}\Gamma_{n}Z_{\aleph_{s}}}{n_{1,s}}\left[\frac{Z_{\aleph_{s}}^{\top}\Gamma_{n}^{\top}\Gamma_{n}Z_{\aleph_{s}}}{n_{1,s}}\right]^{-1}\frac{Z_{\aleph_{s}}^{\top}\Gamma_{n}^{\top}1_{n_{1,s}}}{n_{1,s}}$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\frac{1_{n_{1,s}}^{\top}W_{n}^{1/2}Z_{\aleph_{1,s}}}{n_{1,s}}\left[\frac{Z_{\aleph_{1,s}}^{\top}W_{n}Z_{\aleph_{1,s}}}{n_{1,s}}\right]^{-1}\frac{Z_{\aleph_{1,s}}^{\top}W_{n}^{1/2}1_{n_{1,s}}}{n_{1,s}}$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\frac{1_{n_{1,s}}^{\top}U_{n}D_{n}^{1/2}Z_{\aleph_{1,s}}}{n_{1,s}}\left[\frac{Z_{\aleph_{1,s}}^{\top}D_{n}Z_{\aleph_{1,s}}}{n_{1,s}}\right]^{-1}\frac{Z_{\aleph_{1,s}}^{\top}D_{n}^{1/2}U_{n}^{\top}1_{n_{1,s}}}{n_{1,s}}$
$\displaystyle=\frac{(1-n_{1,s}/n_{s})}{n_{1,s}}Z_{n_{1,s}}^{\top}\left[\frac{1}{n_{1,s}}\sum_{i\in[n_{1,s}-1]}Z_{i}Z_{i}^{\top}+\left(\frac{1}{n_{1,s}}-\frac{1}{n_{s}}\right)Z_{n_{1,s}}Z_{n_{1,s}}^{\top}\right]^{-1}Z_{n_{1,s}}$
$\displaystyle=\frac{(1-n_{1,s}/n_{s})}{n_{1,s}}Z_{n_{1,s}}^{\top}\left[\mathcal{A}_{n}^{-1}-\frac{\mathcal{A}_{n}^{-1}(\left(\frac{1}{n_{1,s}}-\frac{1}{n_{s}}\right)Z_{n_{1,s}}Z_{n_{1,s}}^{\top})\mathcal{A}_{n}^{-1}}{1+\left(\frac{1}{n_{1,s}}-\frac{1}{n_{s}}\right)Z_{n_{1,s}}^{\top}\mathcal{A}_{n}^{-1}Z_{n_{1,s}}}\right]Z_{n_{1,s}},$
(A.1)

where $W_{n}=U_{n}D_{n}U_{n}^{\top}$ is the eigenvalue decomposition of
$W_{n}$, $D_{n}=\text{diag}(1,\cdots,1,1-n_{1,s}/n_{s})$ is the eigenvalues,
$U_{n}$ is the corresponding eigenvectors,
$\mathcal{A}_{n}=\frac{1}{n_{1,s}}\sum_{i\in[n_{1,s}-1]}Z_{i}Z_{i}^{\top}$,
the fourth equality is by
$1_{n_{1,s}}^{\top}U_{n}=(0,\cdots,0,\sqrt{n_{1,s}})$, and the last equality
follows the Sherman-Morrison matrix inverse formula. Consider the eigenvalue
decomposition of
$\frac{n_{1,s}}{n_{1,s}-1}\mathcal{A}_{n}=\mathcal{U}_{n}\Lambda_{n}\mathcal{U}_{n}^{\top}$
where $\Lambda_{n}=\text{diag}(\lambda_{1},\cdots,\lambda_{k_{n}})$. We note
that $Z_{n_{1,s}}$ and $\mathcal{A}_{n}$, and thus, $\mathcal{U}_{n}$ and
$\Lambda_{n}$, are independent by construction. Then, for any bounded and
smooth function $f:\Re^{k_{n}}\mapsto\Re$, we have

$\displaystyle\mathbb{E}[f(\mathcal{U}_{n}^{\top}Z_{n_{1,s}}\mid\Lambda_{n}]=\mathbb{E}[\mathbb{E}(f(\mathcal{U}_{n}^{\top}Z_{n_{1,s}}\mid\mathcal{U}_{n},\Lambda_{n})\mid\Lambda_{n}]=\mathbb{E}f(Z_{n_{1,s}})=\mathbb{E}f(\mathcal{U}_{n}^{\top}Z_{n_{1,s}}),$

where we use the fact that the distribution of $Z_{n_{1,s}}$ ($k_{n}\times 1$
standard normal random vector) is invariant to rotations. This implies
$\mathcal{U}_{n}^{\top}Z_{n_{1,s}}$ and $\Lambda_{n}$ are independent. Denote
$\mathcal{U}_{n}^{\top}Z_{n_{1,s}}$ as
$\mathcal{G}=(g_{1},\cdots,g_{K_{n}})^{\top}$ which are independent of
$\Lambda_{n}$ and suppose $k_{n}/n_{1,s}\rightarrow\kappa_{1,s}$. Then, we
have

$\displaystyle\frac{1}{n_{1,s}}Z_{n_{1,s}}^{\top}\left(\frac{n_{1,s}}{n_{1,s}-1}\mathcal{A}_{n}\right)^{-1}Z_{n_{1,s}}$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{n_{1,s}}\sum_{l=1}^{k_{n}}\lambda_{l}^{-1}g_{l}^{2}$
$\displaystyle=\frac{1}{n_{1,s}}\sum_{l=1}^{k_{n}}\lambda_{l}^{-1}+\frac{1}{n_{1,s}}\sum_{l=1}^{k_{n}}\lambda_{l}^{-1}(g_{l}^{2}-1)$
$\displaystyle=\frac{k_{n}}{n_{1,s}}\frac{1}{k_{n}}\sum_{l=1}^{k_{n}}\lambda_{l}^{-1}+o_{P}(1)$
$\displaystyle\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\int_{\lambda_{-}}^{\lambda_{+}}\frac{\sqrt{(\lambda_{+}-\lambda)(\lambda-\lambda_{-})}}{2\pi\lambda^{2}}d\lambda,$
(A.2)

where $\lambda_{\pm}=(1\pm\sqrt{\kappa_{1,s}})^{2}$. We note that the third
equality in (A) holds because $\lambda_{k_{n}}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}1-\sqrt{\kappa_{1,s}}$ which is bounded away from zero,

$\displaystyle\mathbb{E}\left(\frac{1}{n_{1,s}}\sum_{l=1}^{k_{n}}\lambda_{l}^{-1}(g_{l}^{2}-1)\mid\Lambda_{n}\right)=0,$

and

$\displaystyle
Var\left(\frac{1}{n_{1,s}}\sum_{l=1}^{k_{n}}\lambda_{l}^{-1}(g_{l}^{2}-1)\mid\Lambda_{n}\right)\leq\frac{k_{n}}{n_{1,s}^{2}\lambda_{k_{n}}^{2}}\rightarrow
0,$

and the last convergence in probability (A) holds by the Marčenko-Pastur
theorem (see, for example, Bai (2008)). Denote
$\int_{\lambda_{-}}^{\lambda_{+}}\frac{\sqrt{(\lambda_{+}-\lambda)(\lambda-\lambda_{-})}}{2\pi\lambda^{2}}d\lambda$
as $\zeta_{1,s}$, then we have

$\displaystyle\frac{1}{n_{1,s}}Z_{n_{1,s}}^{\top}\mathcal{A}_{n}^{-1}Z_{n_{1,s}}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\int_{\lambda_{-}}^{\lambda_{+}}\frac{\sqrt{(\lambda_{+}-\lambda)(\lambda-\lambda_{-})}}{2\pi\lambda^{2}}d\lambda,$

and thus, by (A),

$\displaystyle\gamma_{1,s,n}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\frac{1}{1+(1-\pi_{s})\zeta_{1,s}}\equiv\gamma_{1,s,\infty}.$

##  Appendix B Uniform Size Control

Recall the Wald test $\mathbb{W}_{n}$ defined in (4.4). By Theorems 4.1 and
5.1, we can show that this Wald test has the uniform asymptotic size control
under the null for a wide range of DGPs that have either diverging or fixed
dimension of covariates.

Let $P_{n}$ be a DGP with sample size $n$ in a class of DGPs
$\mathcal{P}_{n}$. Then, all parameters such as $\Sigma$ from Theorem 3.1 and
$\beta_{a,s}^{*}$ for $a=0,1,s\in\mathcal{S}$ from Theorem 5.1 are implicitly
indexed by $P_{n}$. We consider the sequence of classes of DGPs
($\\{\mathcal{P}_{n}\\}_{n\geq 1}$) that satisfies the following condition:

$\displaystyle\begin{pmatrix}&\text{for any subsequence $\\{n_{l}\\}$ of
$\\{n\\}$, there exists a subsequence $\\{n_{l^{\prime}}\\}$ of
$\\{n_{l}\\}$}\\\ &\text{such that along $\\{P_{n_{l^{\prime}}}\\}$ for
$P_{n_{l^{\prime}}}\in\mathcal{P}_{n_{l^{\prime}}}$},\\\ &\text{either
$k_{n_{l^{\prime}}}\rightarrow\infty$, Assumptions
\ref{ass:assignment1}--\ref{ass:variance} hold, and
$\liminf_{l^{\prime}\rightarrow\infty}\Sigma_{1,1}-2\Sigma_{1,2}+\Sigma_{2,2}>0$}\\\
&\text{or $k_{n_{l^{\prime}}}=k$ is fixed, Assumptions \ref{ass:assignment1}
and \ref{ass:reg_fixed} hold, and
$\liminf_{l^{\prime}\rightarrow\infty}|\pi_{s}\beta_{0,s}^{*}+(1-\pi_{s})\beta_{1,s}|>0$
.}\end{pmatrix}$ (B.1)

######  Theorem B.1.

Suppose $\\{\mathcal{P}_{n}\\}_{n\geq 1}$ satisfies (B.1). Then, under the
null hypothesis that $\tau_{0}=\tau$, we have

$\displaystyle\liminf_{n\rightarrow\infty}\inf_{P_{n}\in\mathcal{P}_{n}}\mathbb{E}_{P_{n}}(\mathbb{W}_{n})=\limsup_{n\rightarrow\infty}\sup_{P_{n}\in\mathcal{P}_{n}}\mathbb{E}_{P_{n}}(\mathbb{W}_{n})=\alpha.$

##  Appendix C Proof of Theorem 3.1

Denote $\hat{\tau}_{s}=\hat{\tau}_{1,s}-\hat{\tau}_{0,s}$,
$\tau_{a,s}=\mathbb{E}(Y(a)|S=s)$ for $a=0,1$, and
$\tau_{s}=\mathbb{E}(Y(1)-Y(0)|S=s)$. Note that

$\displaystyle Y_{i}(1)$
$\displaystyle=\alpha_{1,S_{i}}+X_{i}^{\top}\beta_{1,S_{i}}+e_{i,S_{i}}(1)+\varepsilon_{i}(1)$
$\displaystyle=\alpha_{1,S_{i}}+\overline{X}_{S_{i}}^{\top}\beta_{1,S_{i}}+\breve{X}_{i}^{\top}\beta_{1,S_{i}}+e_{i,S_{i}}(1)+\varepsilon_{i}(1),$

which implies

$\displaystyle\hat{\tau}_{1,s}$
$\displaystyle=(1_{n_{1,s}}^{\top}M_{1,s}1_{n_{1,s}})^{-1}(1_{n_{1,s}}^{\top}M_{1,s}Y_{\aleph_{1,s}}(1))$
$\displaystyle=\alpha_{1,s}+\overline{X}_{s}^{\top}\beta_{1,s}+(1_{n_{1,s}}^{\top}M_{1,s}1_{n_{1,s}})^{-1}(1_{n_{1,s}}^{\top}M_{1,s}(e_{\aleph_{1,s}}(1)+\varepsilon_{\aleph_{1,s}}(1)))$
$\displaystyle=\frac{1}{n_{s}}\sum_{i\in[n]}1\\{S_{i}=s\\}\left(\mathbb{E}(Y_{i}(1)|X_{i},S_{i}=s)-e_{i,s}(1)\right)$
$\displaystyle+(1_{n_{1,s}}^{\top}M_{1,s}1_{n_{1,s}})^{-1}(1_{n_{1,s}}^{\top}M_{1,s}(e_{\aleph_{1,s}}(1)+\varepsilon_{\aleph_{1,s}}(1))).$

By Assumption 2(iv), we have

$\displaystyle\left|\frac{1}{n_{s}}\sum_{i\in[n]}1\\{S_{i}=s\\}e_{i,s}(1)\right|\leq\left(\frac{1}{n_{s}}\sum_{i\in[n]}1\\{S_{i}=s\\}e_{i,s}^{2}(1)\right)^{1/2}=o_{P}(n^{-1/2})$

and

$\displaystyle\left|n_{1,s}^{-1}(1_{n_{1,s}}^{\top}M_{1,s}e_{\aleph_{1,s}}(1))\right|$
$\displaystyle\leq
n_{1,s}^{-1}||M_{1,s}||_{op}||1_{n_{1,s}}||_{2}||e_{\aleph_{1,s}}(1)||_{2}$
$\displaystyle\leq\left(\frac{1}{n_{s}}\sum_{i\in[n]}1\\{S_{i}=s\\}e_{i,s}^{2}(1)\right)^{1/2}=o_{P}(n^{-1/2}).$

This implies the linear expansion

$\displaystyle\hat{\tau}_{1,s}-\tau_{1,s}$
$\displaystyle=\frac{1}{n_{s}}\sum_{i\in[n]}1\\{S_{i}=s\\}\left(\mathbb{E}(Y_{i}(1)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(1)|S_{i}=s)\right)$
$\displaystyle+\gamma_{1,s,n}^{-1}n_{1,s}^{-1}(1_{n_{1,s}}^{\top}M_{1,s}\varepsilon_{\aleph_{1,s}}(1))+o_{P}(n^{-1/2}).$

Similarly, we have

$\displaystyle\hat{\tau}_{0,s}-\tau_{0,s}$
$\displaystyle=\frac{1}{n_{s}}\sum_{i\in[n]}1\\{S_{i}=s\\}\left(\mathbb{E}(Y_{i}(0)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(0)|S_{i}=s)\right)$
$\displaystyle+\gamma_{0,s,n}^{-1}n_{0,s}^{-1}(1_{n_{0,s}}^{\top}M_{0,s}\varepsilon_{\aleph_{0,s}}(0))+o_{P}(n^{-1/2}),$

and thus

$\displaystyle\hat{\tau}_{s}-\tau_{s}$
$\displaystyle=\frac{1}{n_{s}}\sum_{i\in[n]}1\\{S_{i}=s\\}\left[\mathbb{E}(Y_{i}(1)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(0)|X_{i},S_{i}=s)-\tau_{s}\right]$
$\displaystyle+\gamma_{1,s,n}^{-1}n_{1,s}^{-1}(1_{n_{1,s}}^{\top}M_{1,s}\varepsilon_{\aleph_{1,s}}(1))-\gamma_{0,s,n}^{-1}n_{0,s}^{-1}(1_{n_{0,s}}^{\top}M_{0,s}\varepsilon_{\aleph_{0,s}}(0))+o_{P}(n^{-1/2}).$

Therefore, we have

$\displaystyle\hat{\tau}^{adj}-\tau$
$\displaystyle=\sum_{s\in\mathcal{S}}\hat{p}_{s}\left[\tilde{\gamma}_{1,s}^{-1}n_{1,s}^{-1}(1_{n_{1,s}}^{\top}M_{1,s}\varepsilon_{\aleph_{1,s}}(1))-\tilde{\gamma}_{1,s}^{-1}n_{0,s}^{-1}(1_{n_{0,s}}^{\top}M_{0,s}\varepsilon_{\aleph_{0,s}}(0))\right]$
$\displaystyle+\sum_{s\in\mathcal{S}}\frac{\hat{p}_{s}}{n_{s}}\sum_{i\in[n]}1\\{S_{i}=s\\}\left[\mathbb{E}(Y_{i}(1)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(1)|S_{i}=s)\right]$
$\displaystyle-\sum_{s\in\mathcal{S}}\frac{\hat{p}_{s}}{n_{s}}\sum_{i\in[n]}1\\{S_{i}=s\\}\left[\mathbb{E}(Y_{i}(0)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(0)|S_{i}=s)\right]$
$\displaystyle+\frac{1}{n}\sum_{i\in[n]}\left[\mathbb{E}(Y_{i}(1)|S_{i})-\mathbb{E}(Y_{i}(0)|S_{i})-\tau\right]+o_{P}(n^{-1/2}).$

In addition, we have

$\displaystyle\hat{\tau}^{unadj}-\tau$
$\displaystyle=\sum_{s\in\mathcal{S}}\frac{\hat{p}_{s}}{n_{1,s}}\sum_{i\in[n]}A_{i}1\\{S_{i}=s\\}\varepsilon_{i}(1)-\sum_{s\in\mathcal{S}}\frac{\hat{p}_{s}}{n_{0,s}}\sum_{i\in[n]}(1-A_{i})1\\{S_{i}=s\\}\varepsilon_{i}(0)$
$\displaystyle+\sum_{s\in\mathcal{S}}\frac{\hat{p}_{s}}{n_{1,s}}\sum_{i\in[n]}A_{i}1\\{S_{i}=s\\}(\mathbb{E}(Y_{i}(1)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(1)|S_{i}=s))$
$\displaystyle-\sum_{s\in\mathcal{S}}\frac{\hat{p}_{s}}{n_{0,s}}\sum_{i\in[n]}(1-A_{i})1\\{S_{i}=s\\}(\mathbb{E}(Y_{i}(0)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(0)|S_{i}=s))$
$\displaystyle+\frac{1}{n}\sum_{i\in[n]}(\mathbb{E}(Y_{i}(1)|S_{i})-\mathbb{E}(Y_{i}(0)|S_{i})-\tau).$

Therefore, we have

$\displaystyle\sqrt{n}\begin{pmatrix}\hat{\tau}^{adj}-\tau\\\
\hat{\tau}^{unadj}-\tau\end{pmatrix}=\mathcal{U}_{n}+\mathcal{V}_{n}+\begin{pmatrix}\mathcal{W}_{n}\\\
\mathcal{W}_{n}\end{pmatrix}+o_{P}(1),$

where

$\displaystyle\mathcal{U}_{n}=\sum_{s\in\mathcal{S}}\hat{p}_{s}\sqrt{n}\begin{pmatrix}\gamma_{1,s,n}^{-1}n_{1,s}^{-1}(1_{n_{1,s}}^{\top}M_{1,s}\varepsilon_{\aleph_{1,s}}(1))-\gamma_{0,s,n}^{-1}n_{0,s}^{-1}(1_{n_{0,s}}^{\top}M_{0,s}\varepsilon_{\aleph_{0,s}}(0))\\\
n_{1,s}^{-1}(1_{n_{1,s}}^{\top}\varepsilon_{\aleph_{1,s}}(1))-n_{0,s}^{-1}(1_{n_{0,s}}^{\top}\varepsilon_{\aleph_{0,s}}(0))\end{pmatrix},$
(C.1)
$\displaystyle\mathcal{V}_{n}=\sum_{s\in\mathcal{S}}\hat{p}_{s}\sqrt{n}\begin{pmatrix}\frac{1}{n_{s}}\sum_{i\in[n]}1\\{S_{i}=s\\}\left(\phi_{i}(1)-\phi_{i}(0)\right)\\\
\frac{1}{n_{1,s}}\sum_{i\in\aleph_{1,s}}\phi_{i}(1)-\frac{1}{n_{0,s}}\sum_{i\in\aleph_{0,s}}\phi_{i}(0)\end{pmatrix},$
(C.2)
$\displaystyle\mathcal{W}_{n}=\frac{1}{\sqrt{n}}\sum_{i\in[n]}(\mathbb{E}(Y_{i}(1)|S_{i})-\mathbb{E}(Y_{i}(0)|S_{i})-\tau),$
(C.3)

and $\phi_{i}(a)=\mathbb{E}(Y_{i}(a)|X_{i},S_{i})-\mathbb{E}(Y_{i}(a)|S_{i})$.
Then, the desired result holds by Lemma I.1.

##  Appendix D Proof of Theorem 4.1

We note that

$\displaystyle\hat{w}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\frac{\Sigma_{2,2}-\Sigma_{1,2}}{\Sigma_{1,1}+\Sigma_{2,2}-2\Sigma_{1,2}}=\operatorname*{arg\,min}_{w}(w,1-w)\Sigma(w,1-w)^{\top}.$

This implies $\hat{\tau}^{*}$ is weakly more efficient than both
$\hat{\tau}^{adj}$ and $\hat{\tau}^{adj}$ which correspond to $\hat{w}=1$ and
0, respectively. Then, by the continuous mapping theorem, we have

$\displaystyle\sqrt{n}\left[(\hat{w},1-\hat{w})\hat{\Sigma}(\hat{w},1-\hat{w})^{\top}\right]^{-1/2}(\hat{\tau}^{*}-\tau)\rightsquigarrow\mathcal{N}(0,1).$

##  Appendix E Proof of Corollary 4.1

For the weak convergence, we note that

$\displaystyle\sqrt{n}(\hat{\tau}^{*}-\tau_{0})=(\hat{\omega},1-\hat{\omega})\Sigma^{1/2}\begin{pmatrix}\hat{N}_{1}\\\
\hat{N}_{2}\end{pmatrix}\rightsquigarrow(\omega^{*},1-\omega^{*})\Sigma^{1/2}\begin{pmatrix}N_{1}\\\
N_{2}\end{pmatrix}.$

Recall $(a,b)=(1,1)\Sigma^{-1/2}$. Given that both
$N_{1}^{*}=\frac{aN_{1}+bN_{2}}{\sqrt{a^{2}+b^{2}}}$ and the limit
distribution of
$\sqrt{n}\left[(\hat{w},1-\hat{w})\hat{\Sigma}(\hat{w},1-\hat{w})^{\top}\right]^{-1/2}(\hat{\tau}^{*}-\tau_{0})$
have variance 1, it suffices to show that $a/b=\tilde{a}/\tilde{b}$ where
$(\tilde{a},\tilde{b})=(\omega^{*},1-\omega^{*})\Sigma^{1/2}$ and
$\omega^{*}=(\Sigma_{2,2}-\Sigma_{1,2})/(\Sigma_{1,1}+\Sigma_{2,2}-2\Sigma_{1,2})$.
To see this, we note that

$\displaystyle\Sigma^{1/2}=\frac{1}{t}\begin{pmatrix}\Sigma_{1,1}+\eta&\Sigma_{1,2}\\\
\Sigma_{1,2}&\Sigma_{2,2}+\eta\end{pmatrix}\quad\text{and}\quad\Sigma^{-1/2}=\frac{t}{(\Sigma_{1,1}+\eta)(\Sigma_{2,2}+\eta)-\Sigma^{2}_{1,2}}\begin{pmatrix}\Sigma_{2,2}+\eta&-\Sigma_{1,2}\\\
-\Sigma_{1,2}&\Sigma_{1,1}+\eta\end{pmatrix},$

where $\eta=\sqrt{\Sigma_{1,1}\Sigma_{2,2}-\Sigma_{1,2}^{2}}$ and
$t=\sqrt{\Sigma_{1,1}+\Sigma_{2,2}+2\eta}$. Then, it is obvious that

$\displaystyle a$
$\displaystyle=\frac{t(\Sigma_{2,2}-\Sigma_{1,2}+\eta)}{(\Sigma_{1,1}+\eta)(\Sigma_{2,2}+\eta)-\Sigma^{2}_{1,2}},$
$\displaystyle b$
$\displaystyle=\frac{t(\Sigma_{1,1}-\Sigma_{1,2}+\eta)}{(\Sigma_{1,1}+\eta)(\Sigma_{1,1}+\eta)-\Sigma^{2}_{1,2}},$
$\displaystyle\tilde{a}$
$\displaystyle=\frac{(\Sigma_{2,2}-\Sigma_{1,2})(\Sigma_{1,1}+\eta)+(\Sigma_{1,1}-\Sigma_{1,2})\Sigma_{1,2}}{t(\Sigma_{1,1}+\Sigma_{2,2}-2\Sigma_{1,2})}$
$\displaystyle=\frac{(\Sigma_{2,2}\Sigma_{1,1}-\Sigma_{1,2}^{2})+\eta(\Sigma_{2,2}-\Sigma_{1,2})}{t(\Sigma_{1,1}+\Sigma_{2,2}-2\Sigma_{1,2})}$
$\displaystyle=\frac{\eta(\eta+\Sigma_{2,2}-\Sigma_{1,2})}{t(\Sigma_{1,1}+\Sigma_{2,2}-2\Sigma_{1,2})},\quad\text{and}$
$\displaystyle\tilde{b}$
$\displaystyle=\frac{\eta(\eta+\Sigma_{1,1}-\Sigma_{1,2})}{t(\Sigma_{1,1}+\Sigma_{2,2}-2\Sigma_{1,2})},$

which implies $a/b=\tilde{a}/\tilde{b}$.

For the optimality result, note that

$\displaystyle(\mathbb{W}_{n},\psi(\hat{N}_{1},\hat{N}_{2}))\rightsquigarrow(1\\{(N^{*})^{2}\geq\mathcal{C}_{\alpha}\\},\psi(N_{1},N_{2})),$

where $(N_{1},N_{2})^{\top}\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}\left(\begin{pmatrix}a\Delta\\\
b\Delta\end{pmatrix},I_{2}\right)$ and
$N^{*}=\frac{aN_{1}+bN_{2}}{\sqrt{a^{2}+b^{2}}}$. Consider a limit experiment
in which a researcher observes $(N_{1},N_{2})$, knows the values of $(a,b)$,
and wants to test $\Delta=0$ versus the two-sided alternative. Lehmann and
Romano (2006, Section 4.2) show $1\\{(N^{*})^{2}\geq\mathcal{C}_{\alpha}\\}$
is the UMP unbiased test for this limit experiment, implying that

$\displaystyle\mathbb{E}1\\{(N^{*})^{2}\geq\mathcal{C}_{\alpha}\\}\geq\mathbb{E}\psi(N_{1},N_{2})$

for any $\psi\in\Psi_{\alpha}^{U}.$ This implies

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{E}\mathbb{W}_{n}=\mathbb{E}1\\{(N^{*})^{2}\geq\mathcal{C}_{\alpha}\\}=\sup_{\psi\in\Psi_{\alpha}^{U}}\mathbb{E}\psi(N_{1},N_{2})=\sup_{\psi\in\Psi_{\alpha}^{U}}\lim_{n\rightarrow\infty}\mathbb{E}\psi(\hat{N}_{1},\hat{N}_{2}).$

The inequality that

$\displaystyle\sup_{\psi\in\Psi_{\alpha}^{U}}\lim_{n\rightarrow\infty}\mathbb{E}\psi(\hat{N}_{1},\hat{N}_{2})\geq\lim_{n\rightarrow\infty}\mathbb{E}\breve{\psi}_{n}$

holds due to the uniform integrability of $\breve{\psi}_{n}$, the fact that
$\breve{\psi}_{n}=\psi(\hat{N}_{1},\hat{N}_{2})+o_{P}(1)$ for some
$\psi\in\Psi^{U}_{\alpha}$, and the dominated convergence theorem.

##  Appendix F Proof of Theorem 4.2

We derive the limits of $\hat{\Sigma}_{\mathcal{U}}$,
$\hat{\Sigma}_{\mathcal{V}}$, and $\hat{\Sigma}_{\mathcal{W}}$ in the
following three steps.

Step 1: Limit of $\hat{\Sigma}_{\mathcal{U}}$. Following Bugni et al. (2018),
we define $\\{(X_{i}^{s},\varepsilon_{i}^{s}(1),\varepsilon_{i}^{s}(0)):1\leq
i\leq n\\}$ as a sequence of i.i.d. random variables with marginal
distributions equal to the distribution of
$(X_{i},\varepsilon_{i}(1),\varepsilon_{i}(0))|S_{i}=s$ and
$N_{s}=\sum_{i=1}^{n}1\\{S_{i}<s\\}$. We order units by strata and then by
$A_{i}=1$ first and $A_{i}=0$ second within each stratum. This means for the
$s$th stratum, units indexed from $N_{s}+1$ to $N_{s}+n_{1,s}$ are treated and
units indexed from $N_{s}+n_{1,s}+1$ to $N_{s}+n_{s}$ are untreated. Then, we
have $(Y_{i}(a)\mid S_{i}=s)\stackrel{{\scriptstyle d}}{{=}}Y_{i}^{s}(a)$ and
$(\tilde{e}_{i}(a)\mid S_{i}=s)\stackrel{{\scriptstyle d}}{{=}}e_{i,s}(a)$
where

$\displaystyle
Y_{i}^{s}(a)=\alpha_{a,s}+(X_{i}^{s})^{\top}\beta_{a,s}+\tilde{e}_{i}^{s}(a)+\varepsilon_{i}^{s}(a)\quad\text{and}\quad\tilde{e}_{i}^{s}(a)=\mathbb{E}(Y_{i}(a)|X_{i}=x,S_{i}=s)-\alpha_{a,s}-(X_{i}^{s})^{\top}\beta_{a,s}.$

We further define $\tilde{M}_{a,s,i,j}$ as the $(i,j)$th entry of the
$n_{a,s}\times n_{a,s}$ matrix,
$\tilde{M}_{a,s}=I_{n_{a,s}}-\Xi_{a,s}(\Xi_{a,s}^{\top}\Xi_{a,s})^{-1}\Xi_{a,s}^{\top}$,
$\Xi_{1,s}=((\breve{X}_{N_{s}+1}^{s})^{\top},\cdots,(\breve{X}_{N_{s}+n_{1,s}}^{s})^{\top})^{\top}$
is an $n_{1,s}\times k_{n}$ matrix,
$\Xi_{0,s}=((\breve{X}_{N_{s}+1+n_{1,s}}^{s})^{\top},\cdots,(\breve{X}_{N_{s}+n_{s}}^{s})^{\top})^{\top}$
is an $n_{0,s}\times k_{n}$ matrix,

$\displaystyle\breve{X}_{i}^{s}=X_{i}^{s}-\frac{1}{n_{s}}\sum_{j=N_{s}+1}^{N_{s}+n_{s}}X_{j}^{s},\quad\text{if}\quad
N_{s}+1\leq i\leq N_{s}+n_{s},$

$\tilde{\gamma}_{a,s,n}=1_{n_{a,s}}^{\top}\tilde{M}_{a,s}1_{n_{a,s}}$, and
$\phi_{i}^{s}(a)=\mathbb{E}(Y_{i}(a)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(a)|S_{i}=s)$.

Step 1.1: Limit of $\hat{\omega}_{a,s}^{2}$. We consider the case with $a=1$.
The result for $a=0$ can be established in the same manner. Given the above
definitions, we see that, conditionally on $(A^{(n)},S^{(n)})$, and thus,
unconditionally,

$\displaystyle\hat{\omega}_{1,s}^{2}\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{n_{1,s}}\tilde{\gamma}_{1,s,n}^{-2}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\left(\sum_{j=N_{s}+1}^{N_{s}+n_{1,s}}\tilde{M}_{1,s,i,j}\right)^{2}Y_{i}^{s}(1)\tilde{\acute{\varepsilon}}_{1,s,i},$

where
$\tilde{\acute{\varepsilon}}_{1,s,i}=\tilde{\varepsilon}_{1,s,i}/\tilde{M}_{1,s,i,i}$
and $\tilde{\varepsilon}_{1,s,i}$ is the residual from the linear regression
of $Y_{i}^{s}(1)$ on $(1,\breve{X}_{i}^{s})$ with observations $N_{s}+1\leq
i\leq N_{s}+n_{1,s}$. Because $\\{(X_{i}^{s},Y_{i}^{s}(1)):1\leq i\leq n\\}$
are i.i.d. and independent of $(A^{(n)},S^{(n)})$ we can directly apply
Jochmans (2022, Theorem 1) and conclude that

$\displaystyle\hat{\omega}_{1,s}^{2}-\tilde{\gamma}_{1,s,n}^{-2}\sigma_{1,s,n}^{2}=o_{P}(1).$

Then, because $\tilde{\gamma}_{a,s,n}\stackrel{{\scriptstyle
d}}{{=}}\gamma_{a,s,n}$ and Assumption 3 holds, we have

$\displaystyle\hat{\omega}_{1,s}^{2}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\omega_{1,s,\infty}^{2}.$

Step 1.2: Limit of $\hat{\varpi}_{a,s}$. Again, we focus on
$\hat{\varpi}_{1,s}$. Following the argument in Step 1.1, we have

$\displaystyle\hat{\varpi}_{1,s}\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{n_{1,s}}\tilde{\gamma}_{1,s,n}^{-1}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\left(\sum_{j=N_{s}+1}^{N_{s}+n_{1,s}}\tilde{M}_{1,s,i,j}\right)Y_{i}^{s}(1)\tilde{\acute{\varepsilon}}_{1,s,i}.$

Note that

$\displaystyle\tilde{\acute{\varepsilon}}_{1,s,i}$
$\displaystyle=\frac{\sum_{k=N_{s}+1}^{N_{s}+n_{1,s}}\tilde{M}_{1,s,i,k}(Y_{k}^{s}(1)-\tilde{\tau}_{1,s})}{\tilde{M}_{1,s,i,i}}$
$\displaystyle=\frac{\sum_{k=N_{s}+1}^{N_{s}+n_{1,s}}\tilde{M}_{1,s,i,k}(\tilde{e}_{k}^{s}(1)+\alpha_{1,s}+\widetilde{X}_{s}^{\top}\beta_{1,s}-\tilde{\tau}_{1,s}+\varepsilon_{k}^{s}(1))}{\tilde{M}_{1,s,i,i}},$

where $\tilde{\tau}_{1,s}$ is the intercept of the OLS regression of
$Y_{i}^{s}(1)$ on $(1,\breve{X}_{i}^{s})$ using observations $N_{s}+1\leq
i\leq N_{s}+n_{1,s}$ and
$\widetilde{X}_{s}=\frac{1}{n_{s}}\sum_{i=N_{s}+1}^{N_{s}+n_{s}}X_{i}^{s}$. By
construction, we have $\tilde{\tau}_{1,s}\stackrel{{\scriptstyle
d}}{{=}}\hat{\tau}_{1,s}$. Let
$\theta_{i}=\sum_{k=N_{s}+1}^{N_{s}+n_{1,s}}\tilde{M}_{1,s,i,k}$. Therefore,
we have

$\displaystyle\frac{1}{n_{1,s}}\tilde{\gamma}_{1,s,n}^{-1}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\left(\sum_{j=N_{s}+1}^{N_{s}+n_{1,s}}\tilde{M}_{1,s,i,j}\right)Y_{i}^{s}(1)\tilde{\acute{\varepsilon}}_{1,s,i}$
$\displaystyle=\frac{1}{n_{1,s}}\tilde{\gamma}_{1,s,n}^{-1}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\sum_{k=N_{s}+1}^{N_{s}+n_{1,s}}\frac{\theta_{i}Y_{i}^{s}(1)\tilde{M}_{1,s,i,k}\varepsilon_{k}^{s}(1)}{\tilde{M}_{1,s,i,i}}+\frac{1}{n_{1,s}}\tilde{\gamma}_{1,s,n}^{-1}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\sum_{k=N_{s}+1}^{N_{s}+n_{1,s}}\frac{\theta_{i}Y_{i}^{s}(1)\tilde{M}_{1,s,i,k}\tilde{e}_{k}^{s}(1)}{\tilde{M}_{1,s,i,i}}$
$\displaystyle+\frac{1}{n_{1,s}}\tilde{\gamma}_{1,s,n}^{-1}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\frac{\theta_{i}^{2}Y_{i}^{s}(1)}{\tilde{M}_{1,s,i,i}}(\alpha_{1,s}+\widetilde{X}_{s}^{\top}\beta_{1,s}-\tilde{\tau}_{1,s}))$
$\displaystyle\equiv T_{1}+T_{2}+T_{3}.$

For $T_{2}$, we have

$\displaystyle|T_{2}|$
$\displaystyle\leq\tilde{\gamma}_{1,s,n}^{-1}\left[\frac{1}{n_{1,s}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\frac{\theta_{i}^{2}(Y_{i}^{s}(1))^{2}}{\tilde{M}_{1,s,i,i}^{2}}\right]^{1/2}\left[\frac{1}{n_{1,s}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}(\tilde{e}_{i}^{s}(1))^{2}\right]^{1/2}$
$\displaystyle\leq\gamma_{1,s,n}^{-1}\frac{\max_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}|\theta_{i}|}{\min_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}\tilde{M}_{1,s,i,i}}\left[\frac{1}{n_{1,s}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}(Y_{i}^{s}(1))^{2}\right]^{1/2}\left[\frac{1}{n_{1,s}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}(\tilde{e}_{i}^{s}(1))^{2}\right]^{1/2}$
$\displaystyle=o_{P}(1),$

where the first inequality is by $||\tilde{M}_{1,s}||_{op}\leq 1$ and the last
equality holds because by Assumption 2, we have

$\displaystyle\min_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}\tilde{M}_{1,s,i,i}\stackrel{{\scriptstyle
d}}{{=}}\min_{N_{s}+1\leq i\in\aleph_{1,s}}M_{1,s,i,i}\geq\delta>0,$
$\displaystyle\frac{1}{n_{1,s}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}(\tilde{e}_{i}^{s}(1))^{2}\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{n_{1,s}}\sum_{i\in\aleph_{1,s}}e_{i,s}^{2}(1)=o_{P}(n^{-1}),\quad\text{and}$
$\displaystyle\max_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}|\theta_{i}|\stackrel{{\scriptstyle
d}}{{=}}\max_{i\in\aleph_{1,s}}|\sum_{j\in\aleph_{1,s}}M_{1,s,i,j}|=o_{P}(n^{1/2}).$

For $T_{3}$, we have

$\displaystyle|T_{3}|$
$\displaystyle\leq\tilde{\gamma}_{1,s,n}^{-1}\left[\frac{1}{n_{1,s}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\frac{\theta_{i}^{4}}{\tilde{M}_{1,s,i,i}^{2}}\right]^{1/2}\left[\frac{1}{n_{1,s}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}(Y_{i}^{s}(1))^{2}\right]^{1/2}|\alpha_{1,s}+\widetilde{X}_{s}^{\top}\beta_{1,s}-\tilde{\tau}_{1,s}|$
$\displaystyle\leq\gamma_{1,s,n}^{-1}\frac{\max_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}|\theta_{i}|}{\min_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}\tilde{M}_{1,s,i,i}}\left[\frac{1}{n_{1,s}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\theta_{i}^{2}\right]^{1/2}$
$\displaystyle\times\left[\frac{1}{n_{1,s}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}(Y_{i}^{s}(1))^{2}\right]^{1/2}|\alpha_{1,s}+\widetilde{X}_{s}^{\top}\beta_{1,s}-\tilde{\tau}_{1,s})|$
$\displaystyle=o_{P}(1),$

where the last equality holds because

$\displaystyle\frac{1}{n_{1,s}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\theta_{i}^{2}=\frac{1}{n_{1,s}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\sum_{j=N_{s}+1}^{N_{s}+n_{1,s}}\tilde{M}_{1,s,i,j}=\frac{1}{n_{1,s}}1_{n_{1,s}}^{\top}\tilde{M}_{1,s}1_{n_{1,s}}\leq||\tilde{M}_{1,s}||_{op}\leq
1$

and

$\displaystyle|\alpha_{1,s}+\widetilde{X}_{s}^{\top}\beta_{1,s}-\tilde{\tau}_{1,s})|\stackrel{{\scriptstyle
d}}{{=}}|\alpha_{1,s}+\overline{X}_{s}^{\top}\beta_{1,s}-\hat{\tau}_{1,s})|$
$\displaystyle=|(1_{n_{1,s}}^{\top}M_{1,s}1_{n_{1,s}})^{-1}(1_{n_{1,s}}^{\top}M_{1,s}\varepsilon_{\aleph_{1,s}}(1))|+o_{P}(n^{-1/2})$
$\displaystyle=O_{P}(n^{-1/2}).$

Next, we focus on $T_{1}$. Denote
$\mu_{i}^{s}(1)=\mathbb{E}(Y_{i}^{s}(1)|X_{i}^{s})\stackrel{{\scriptstyle
d}}{{=}}\mathbb{E}(Y_{i}(1)|X_{i},S_{i}=s)$. Then, we have

$\displaystyle T_{1}$
$\displaystyle=\frac{1}{n_{1,s}}\tilde{\gamma}_{1,s,n}^{-1}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\sum_{k=N_{s}+1}^{N_{s}+n_{1,s}}\frac{\theta_{i}\mu_{i}^{s}(1)M_{1,s,i,k}\varepsilon_{k}^{s}(1)}{\tilde{M}_{1,s,i,i}}$
$\displaystyle+\frac{1}{n_{1,s}}\tilde{\gamma}_{1,s,n}^{-1}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\sum_{N_{s}+1\leq
k\leq N_{s}+n_{1,s},k\neq
i}\frac{\theta_{i}\varepsilon_{i}^{s}(1)M_{1,s,i,k}\varepsilon_{k}^{s}(1)}{\tilde{M}_{1,s,i,i}}$
$\displaystyle+\frac{1}{n_{1,s}}\tilde{\gamma}_{1,s,n}^{-1}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\theta_{i}[(\varepsilon_{i}^{s}(1))^{2}-\tilde{V}_{1,s,i}^{2}]+\frac{1}{n_{1,s}}\tilde{\gamma}_{1,s,n}^{-1}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\theta_{i}\tilde{V}_{1,s,i}^{2}$
$\displaystyle\equiv T_{1,1}+T_{1,2}+T_{1,3}+T_{1,4},$

where
$\tilde{V}_{1,s,i}^{2}=\mathbb{E}[(\varepsilon_{i}^{s}(1))^{2}|X_{i}^{s}]\stackrel{{\scriptstyle
d}}{{=}}\mathbb{E}(\varepsilon_{i}^{2}(1)|X_{i},S_{i}=s)$. Recall
$\Xi_{1,s}=((\breve{X}_{N_{s}+1}^{s})^{\top},\cdots,(\breve{X}_{N_{s}+n_{1,s}}^{s})^{\top})^{\top}$.
We have $\mathbb{E}(T_{1,1}\mid\Xi_{1,s})=0$ and

$\displaystyle Var(T_{1,1}\mid\Xi_{1,s})$
$\displaystyle=\frac{1}{n_{1,s}^{2}\tilde{\gamma}_{1,s,n}^{2}}\sum_{k=N_{s}+1}^{N_{s}+n_{1,s}}\tilde{V}_{1,s,k}^{2}\left[\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\frac{\theta_{i}\mu_{i}^{s}(1)M_{1,s,i,k}}{\tilde{M}_{1,s,i,i}}\right]^{2}$
$\displaystyle\leq\frac{\max_{i\in[n]}\tilde{V}_{1,s,i}^{2}}{n_{1,s}^{2}\tilde{\gamma}_{1,s,n}^{2}}\sum_{k=N_{s}+1}^{N_{s}+n_{1,s}}\left[\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\frac{\theta_{i}\mu_{i}^{s}(1)M_{1,s,i,k}}{\tilde{M}_{1,s,i,i}}\right]\left[\sum_{j=N_{s}+1}^{N_{s}+n_{1,s}}\frac{\theta_{j}\mu_{j}^{s}(1)M_{1,s,j,k}}{\tilde{M}_{1,s,j,j}}\right]$
$\displaystyle=\frac{\max_{i\in[n]}\tilde{V}_{1,s,i}^{2}}{n_{1,s}^{2}\tilde{\gamma}_{1,s,n}^{2}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\sum_{j=N_{s}+1}^{N_{s}+n_{1,s}}\left[\frac{\theta_{i}\mu_{i}^{s}(1)\theta_{j}\mu_{j}^{s}(1)M_{1,s,i,j}}{\tilde{M}_{1,s,i,i}\tilde{M}_{1,s,j,j}}\right]$
$\displaystyle\leq\frac{\max_{i\in[n]}\tilde{V}_{1,s,i}^{2}}{n_{1,s}^{2}\tilde{\gamma}_{1,s,n}^{2}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\left[\frac{\theta_{i}\mu_{i}^{s}(1)}{\tilde{M}_{1,s,i,i}}\right]^{2}$
$\displaystyle\leq\left[\frac{\max_{i\in[n]}\tilde{V}_{1,s,i}^{2}}{n_{1,s}\tilde{\gamma}_{1,s,n}^{2}\min_{N_{s}+1\leq
i\leq
N_{s}+n_{1,s}}\tilde{M}_{1,s,i,i}^{2}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}(\mu_{i}^{s}(1))^{2}\right]\left[\frac{\max_{N_{s}+1\leq
i\leq N_{s}+n_{1,s}}\theta_{i}^{2}}{n_{1,s}}\right]$ $\displaystyle=o_{P}(1),$

where the second equality holds because the matrix $\tilde{M}_{1,s}$ is
idempotent, the second inequality holds because $||\tilde{M}_{1,s}||_{op}\leq
1$, and the last equality holds because by Assumption 2, we have

$\displaystyle\min_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}\tilde{M}_{1,s,i,i}\stackrel{{\scriptstyle
d}}{{=}}\min_{N_{s}+1\leq
i\in\aleph_{1,s}}M_{1,s,i,i}\geq\delta>0,\quad\max_{i\in[n]}\tilde{V}_{1,s,i}^{2}=O_{P}(1),\quad\text{and}$
$\displaystyle\max_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}|\theta_{i}|\stackrel{{\scriptstyle
d}}{{=}}\max_{i\in\aleph_{1,s}}\left|\sum_{j\in\aleph_{1,s}}M_{1,s,i,j}\right|=o_{P}(n^{1/2}).$

This implies $T_{1,1}=o_{P}(1)$.

For $T_{1,2}$, we have $\mathbb{E}(T_{1,2}\mid\Xi_{1,s})=0$ and

$\displaystyle Var(T_{1,2}\mid\Xi_{1,s})$
$\displaystyle=\frac{1}{n_{1,s}^{2}}\tilde{\gamma}_{1,s,n}^{-2}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\sum_{N_{s}+1\leq
k\leq N_{s}+n_{1,s},k\neq
i}\frac{\theta_{i}^{2}\tilde{M}_{1,s,i,k}^{2}\tilde{V}_{1,s,i}^{2}\tilde{V}_{1,s,k}^{2}}{\tilde{M}_{1,s,i,i}^{2}}$
$\displaystyle+\frac{1}{n_{1,s}^{2}}\tilde{\gamma}_{1,s,n}^{-2}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\sum_{N_{s}+1\leq
k\leq N_{s}+n_{1,s},k\neq
i}\frac{|\theta_{i}\theta_{k}|\tilde{M}_{1,s,i,k}^{2}\tilde{V}_{1,s,i}^{2}\tilde{V}_{1,s,k}^{2}}{\tilde{M}_{1,s,i,i}^{2}}$
$\displaystyle\leq\frac{\max_{i\in[n]}\tilde{V}_{1,s,i}^{4}}{n_{1,s}^{2}\tilde{\gamma}_{1,s,n}^{2}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\frac{\theta_{i}^{2}\sum_{k=N_{s}+1}^{N_{s}+n_{1,s}}\tilde{M}_{1,s,i,k}^{2}}{\tilde{M}_{1,s,i,i}^{2}}$
$\displaystyle+\frac{\max_{i\in[n]}\tilde{V}_{1,s,i}^{4}}{n_{1,s}^{2}\tilde{\gamma}_{1,s,n}^{2}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\sum_{k=N_{s}+1}^{N_{s}+n_{1,s}}\frac{|\theta_{i}\theta_{k}|\tilde{M}_{1,s,i,k}^{2}}{\tilde{M}_{1,s,i,i}^{2}}$
$\displaystyle\leq 2\left[\frac{\max_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}\tilde{V}_{1,s,i}^{4}\max_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}\theta_{i}^{2}}{n_{1,s}\tilde{\gamma}_{1,s,n}^{2}\min_{N_{s}+1\leq
i\leq N_{s}+n_{1,s}}\tilde{M}_{1,s,i,i}}\right]$ $\displaystyle=o_{P}(1),$

where the last equality is by $\max_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}\theta_{i}^{2}=o_{P}(n)$.

For $T_{1,3}$, we have $\mathbb{E}(T_{1,3}\mid\Xi_{1,s})=0$ and

$\displaystyle Var(T_{1,3}\mid\Xi_{1,s})$
$\displaystyle=\frac{\tilde{\gamma}_{1,s,n}^{-2}}{n_{1,s}^{2}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\theta_{i}^{2}\mathbb{E}\left([(\varepsilon_{i}^{s}(1))^{2}-\tilde{V}_{1,s,i}^{2}]^{2}\mid\Xi_{1,s}\right)$
$\displaystyle\leq\frac{\tilde{\gamma}_{1,s,n}^{-2}\max_{N_{s}+1\leq i\leq
N_{s}+n_{1,s}}\mathbb{E}\left([(\varepsilon_{i}^{s}(1))^{2}-\tilde{V}_{1,s,i}^{2}]^{2}\mid\Xi_{1,s}\right)}{n_{1,s}^{2}}\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\theta_{i}^{2}=O_{P}(n^{-1}),$

where we use the fact that
$\sum_{i=N_{s}+1}^{N_{s}+n_{1,s}}\theta_{i}^{2}=1_{n_{1,s}}^{\top}\tilde{M}_{1,s}1_{n_{1,s}}\leq
n_{1,s}$. This implies $T_{1,3}=o_{P}(1)$.

Last, by Assumption 3, we have

$\displaystyle T_{1,4}\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{n_{1,s}}\gamma_{1,s,n}^{-1}\sum_{i\in\aleph_{1,s}}(\sum_{j\in\aleph_{1,s}}M_{1,s,i,j})\mathbb{E}(\varepsilon_{i}^{2}(1)|X_{i},S_{i}=s)\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\varpi_{a,s,\infty}.$

This concludes the proof of Step 1.2.

Step 1.3: Consistency of $\hat{\Sigma}_{\mathcal{U}}$. We note that
$n_{s}^{2}/(nn_{1,s})\stackrel{{\scriptstyle
p}}{{\longrightarrow}}p_{s}/\pi_{s}$ and
$n_{s}^{2}/(nn_{0,s})\stackrel{{\scriptstyle
p}}{{\longrightarrow}}p_{s}/(1-\pi_{s})$. Therefore, we have

$\displaystyle\hat{\Sigma}_{\mathcal{U}}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\sum_{s\in\mathcal{S}}p_{s}\begin{pmatrix}\frac{\omega_{1,s,\infty}^{2}}{\pi_{s}}+\frac{\omega_{0,s,\infty}^{2}}{1-\pi_{s}}&\frac{\omega_{1,s,\infty}^{2}}{\pi_{s}}+\frac{\omega_{0,s,\infty}^{2}}{1-\pi_{s}}\\\
\frac{\omega_{1,s,\infty}^{2}}{\pi_{s}}+\frac{\omega_{0,s,\infty}^{2}}{1-\pi_{s}}&\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\end{pmatrix}.$

Step 2: Limit of $\hat{\Sigma}_{\mathcal{V}}^{adj}$. It suffices to consider
the limits of
$(\hat{\beta}_{a,s}^{\top}\Gamma_{a,s}\hat{\beta}_{a,s}-\sum_{i\in\aleph_{a,s}}P_{a,s,i,i}Y_{i}\acute{\varepsilon}_{a,s,i})/n_{a,s}$
and $\hat{\beta}_{1,s}^{\top}\Gamma_{s}\hat{\beta}_{0,s}/n_{s}$, which are
derived in the following two steps.

Step 2.1: Limit of
$(\hat{\beta}_{a,s}^{\top}\Gamma_{a,s}\hat{\beta}_{a,s}-\sum_{i\in\aleph_{a,s}}P_{a,s,i,i}Y_{i}\acute{\varepsilon}_{a,s,i})/n_{a,s}$.
We note that

$\displaystyle\hat{\beta}_{a,s}$
$\displaystyle=\Gamma_{a,s}^{-1}\sum_{i\in\aleph_{a,s}}\breve{X}_{i}(Y_{i}-\hat{\tau}_{a}(s))$
$\displaystyle=\Gamma_{a,s}^{-1}\sum_{i\in\aleph_{a,s}}\breve{X}_{i}\left[(\alpha_{a,s}+\overline{X}_{s}^{\top}\beta_{a,s}-\hat{\tau}_{a}(s))+\breve{X}_{i}^{\top}\beta_{a,s}+e_{i,s}(a)+\varepsilon_{i}(a)\right]$
$\displaystyle=\beta_{a,s}+\Gamma_{a,s}^{-1}\sum_{i\in\aleph_{a,s}}\breve{X}_{i}\varepsilon_{i}(a)+\Gamma_{a,s}^{-1}\sum_{i\in\aleph_{a,s}}\breve{X}_{i}\left[(\alpha_{a,s}+\overline{X}_{s}^{\top}\beta_{a,s}-\hat{\tau}_{a}(s))+e_{i,s}(a)\right].$
(F.1)

By (F.1), we have

$\displaystyle\hat{\beta}_{a,s}^{\top}\Gamma_{a,s}\hat{\beta}_{a,s}-\sum_{i\in\aleph_{a,s}}P_{a,s,i,i}Y_{i}\acute{\varepsilon}_{a,s,i}$
$\displaystyle=\beta_{a,s}^{\top}\Gamma_{a,s}\beta_{a,s}+\left(\varepsilon_{\aleph_{a,s}}^{\top}(a)P_{a,s}\varepsilon_{\aleph_{a,s}}(a)-\sum_{i\in\aleph_{a,s}}P_{a,s,i,i}Y_{i}\acute{\varepsilon}_{a,s,i}\right)+2\sum_{i\in\aleph_{a,s}}\breve{X}_{i}^{\top}\beta_{a,s}\varepsilon_{i}(a)$
$\displaystyle+2\sum_{i\in\aleph_{a,s}}\breve{X}_{i}^{\top}\beta_{a,s}(\alpha_{a,s}+\overline{X}_{s}^{\top}\beta_{a,s}-\hat{\tau}_{a}(s))+\varepsilon_{\aleph_{a,s}}^{\top}(a)P_{a,s}\left[1_{n_{a,s}}(\alpha_{a,s}+\overline{X}_{s}^{\top}\beta_{a,s}-\hat{\tau}_{a}(s))+e_{\aleph_{a,s}}(a)\right]$
$\displaystyle+\left[1_{n_{a,s}}(\alpha_{a,s}+\overline{X}_{s}^{\top}\beta_{a,s}-\hat{\tau}_{a}(s))+e_{\aleph_{a,s}}(a)\right]^{\top}P_{a,s}\left[1_{n_{a,s}}(\alpha_{a,s}+\overline{X}_{s}^{\top}\beta_{a,s}-\hat{\tau}_{a}(s))+e_{\aleph_{a,s}}(a)\right]$
$\displaystyle\equiv\beta_{a,s}^{\top}\Gamma_{a,s}\beta_{a,s}+\sum_{l\in[5]}I_{l},$
(F.2)

where $e_{\aleph_{a,s}}(a)$ and $\varepsilon_{\aleph_{a,s}}(a)$ are
$n_{a,s}\times 1$ vectors of $\\{e_{i,s}(a)\\}_{i\in\aleph_{a,s}}$ and
$\\{\varepsilon_{i,s}(a)\\}_{i\in\aleph_{a,s}}$, respectively. We note that

$\displaystyle\sum_{i\in\aleph_{a,s}}(\breve{X}_{i}^{\top}\beta_{a,s})^{2}$
$\displaystyle=\sum_{i\in\aleph_{a,s}}\left[\mathbb{E}(Y_{i}(a)|X_{i},S_{i}=s)-e_{i,s}(a)-\frac{1}{n_{s}}\sum_{j\in\aleph_{s}}\left(\mathbb{E}(Y_{j}(a)|X_{j},S_{j}=s)-e_{j,s}(a)\right)\right]^{2}$
$\displaystyle\lesssim\sum_{i\in\aleph_{a,s}}\left[\mathbb{E}(Y_{i}(a)|X_{i},S_{i}=s)-e_{i,s}(a)\right]^{2}+\frac{n_{a,s}}{n_{s}}\sum_{j\in\aleph_{s}}\left[\left(\mathbb{E}(Y_{j}(a)|X_{j},S_{j}=s)-e_{j,s}(a)\right)\right]^{2}$
$\displaystyle=O_{P}(n),$ (F.3)
$\displaystyle\left|\left[\frac{1}{n_{a,s}}\beta_{a,s}^{\top}\Gamma_{a,s}\beta_{a,s}\right]^{1/2}-\left[\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}(\mathbb{E}(Y_{i}(a)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(a)|S_{i}=s))^{2}\right]^{1/2}\right|$
$\displaystyle=\left|\left[\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}(\breve{X}_{i}^{\top}\beta_{a,s})^{2}\right]^{1/2}-\left[\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}(\mathbb{E}(Y_{i}(a)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(a)|S_{i}=s))^{2}\right]^{1/2}\right|$
$\displaystyle\leq\left[\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}\left(e_{i,s}(a)-\frac{1}{n_{s}}\sum_{j\in\aleph_{s}}e_{j,s}(a)\right)^{2}\right]^{1/2}+\left|\frac{1}{n_{s}}\sum_{j\in\aleph_{s}}\mathbb{E}(Y_{j}(a)|X_{j},S_{j}=s)-\mathbb{E}(Y_{i}(a)|S_{i}=s)\right|$
$\displaystyle=o_{P}(1),$

and

$\displaystyle\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}(\mathbb{E}(Y_{i}(a)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(a)|S_{i}=s))^{2}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}var(\phi_{i}(a)|S_{i}=s),$

where
$\phi_{i}(a)=\mathbb{E}(Y_{i}(a)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(a)|S_{i}=s)$.
This implies

$\displaystyle\frac{1}{n_{a,s}}\beta_{a,s}^{\top}\Gamma_{a,s}\beta_{a,s}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}var(\phi_{i}(a)|S_{i}=s).$

Next, we show $I_{1},\cdots,I_{5}$ are $o_{P}(n)$. Lemma I.2 shows
$I_{1}=o_{P}(n)$. For $I_{2}$, we note that, conditional on
$(A^{(n)},S^{(n)},X^{(n)})$, $\\{\varepsilon_{i}(a)\\}_{i\in\aleph_{a,s}}$ is
independent across $i$ and mean zero, which implies

$\displaystyle
Var\left(\sum_{i\in\aleph_{a,s}}\breve{X}_{i}^{\top}\beta_{a,s}\varepsilon_{i}(a)\mid
A^{(n)},S^{(n)},X^{(n)}\right)$
$\displaystyle\leq\left[\max_{i\in\aleph_{a,s}}\mathbb{E}(\varepsilon_{i}^{2}(a)\mid
A^{(n)},S^{(n)},X^{(n)})\right]\left[\sum_{i\in\aleph_{a,s}}(\breve{X}_{i}^{\top}\beta_{a,s})^{2}\right]=O_{P}(n),$

where the equality is by (F). Therefore, we have
$I_{2}=O_{P}(n^{1/2})=o_{P}(n)$. Because
$\alpha_{a,s}+\overline{X}_{s}^{\top}\beta_{a,s}-\hat{\tau}_{a}(s)\stackrel{{\scriptstyle
p}}{{\longrightarrow}}o_{P}(1)$ by the proof of Theorem 3.1,
$\sum_{i\in\aleph_{a,s}}e_{i,s}^{2}(a)=o(n)$, and $\lambda_{\max}(P_{a,s})\leq
1$, we have $I_{l}=o_{P}(n)$ for $l=3,4,5$. This means

$\displaystyle\frac{1}{n_{a,s}}\left(\hat{\beta}_{a,s}^{\top}\Gamma_{a,s}\hat{\beta}_{a,s}-\sum_{i\in\aleph_{a,s}}P_{a,s,i,i}Y_{i}\acute{\varepsilon}_{a,s,i}\right)\stackrel{{\scriptstyle
p}}{{\longrightarrow}}var(\phi_{i}(a)|S_{i}=s).$

Step 2.2: Limit of
$\hat{\beta}_{1,s}^{\top}\Gamma_{s}\hat{\beta}_{0,s}/n_{s}$. By (F.1), we have

$\displaystyle\hat{\beta}_{1,s}^{\top}\Gamma_{s}\hat{\beta}_{0,s}$
$\displaystyle=\beta_{1,s}^{\top}\Gamma_{s}\beta_{0,s}+\beta_{1,s}^{\top}\Gamma_{s}\Gamma_{0,s}^{-1}\sum_{i\in\aleph_{0,s}}\breve{X}_{i}\varepsilon_{i}(0)$
$\displaystyle+\beta_{1,s}^{\top}\Gamma_{s}\Gamma_{0,s}^{-1}\left[\sum_{i\in\aleph_{0,s}}\breve{X}_{i}\left(\alpha_{0,s}+\overline{X}_{s}^{\top}\beta_{0,s}-\hat{\tau}_{0}(s)+e_{i,s}(0)\right)\right]$
$\displaystyle+\left[\sum_{i\in\aleph_{1,s}}\breve{X}_{i}\varepsilon_{i}(1)\right]^{\top}\Gamma_{1,s}^{-1}\Gamma_{s}\beta_{0,s}+\left[\sum_{i\in\aleph_{1,s}}\breve{X}_{i}\varepsilon_{i}(1)\right]^{\top}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}\sum_{i\in\aleph_{0,s}}\breve{X}_{i}\varepsilon_{i}(0)$
$\displaystyle+\left[\sum_{i\in\aleph_{1,s}}\breve{X}_{i}\varepsilon_{i}(1)\right]^{\top}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}\left[\sum_{i\in\aleph_{0,s}}\breve{X}_{i}\left(\alpha_{0,s}+\overline{X}_{s}^{\top}\beta_{0,s}-\hat{\tau}_{0}(s)+e_{i,s}(0)\right)\right]$
$\displaystyle+\left[\sum_{i\in\aleph_{1,s}}\breve{X}_{i}\left(\alpha_{1,s}+\overline{X}_{s}^{\top}\beta_{1,s}-\hat{\tau}_{1}(s)+e_{i,s}(1)\right)\right]^{\top}\Gamma_{1,s}^{-1}\Gamma_{s}\beta_{0,s}$
$\displaystyle+\left[\sum_{i\in\aleph_{1,s}}\breve{X}_{i}\left(\alpha_{1,s}+\overline{X}_{s}^{\top}\beta_{1,s}-\hat{\tau}_{1}(s)+e_{i,s}(1)\right)\right]^{\top}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}\sum_{i\in\aleph_{0,s}}\breve{X}_{i}\varepsilon_{i}(0)$
$\displaystyle+\left[\sum_{i\in\aleph_{1,s}}\breve{X}_{i}\left(\alpha_{1,s}+\overline{X}_{s}^{\top}\beta_{1,s}-\hat{\tau}_{1}(s)+e_{i,s}(1)\right)\right]^{\top}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}$
$\displaystyle\times\left[\sum_{i\in\aleph_{0,s}}\breve{X}_{i}\left(\alpha_{0,s}+\overline{X}_{s}^{\top}\beta_{0,s}-\hat{\tau}_{0}(s)+e_{i,s}(0)\right)\right]$
$\displaystyle\equiv\sum_{l\in[9]}I_{l}.$

Following the previous argument, we can show that

$\displaystyle\frac{I_{1}}{n}$
$\displaystyle=\frac{1}{n}\sum_{i\in\aleph_{s}}(\mathbb{E}(Y_{i}(1)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(1)|S_{i}=s))(\mathbb{E}(Y_{i}(0)|X_{i},S_{i}=s)-\mathbb{E}(Y_{i}(0)|S_{i}=s))+o_{P}(1)$
$\displaystyle\stackrel{{\scriptstyle
p}}{{\longrightarrow}}cov(\phi_{i}(1),\phi_{i}(0)|S_{i}=s).$

For $I_{2}$, we have

$\displaystyle\mathbb{E}(I_{2}|A^{(n)},S^{(n)},X^{(n)})=0$

and

$\displaystyle
var(I_{2}|A^{(n)},S^{(n)},X^{(n)})\lesssim\beta_{1,s}^{\top}\Gamma_{s}\Gamma_{0,s}^{-1}\Gamma_{s}\beta_{1,s}\max_{i\in\aleph_{0,s}}\mathbb{E}(\varepsilon_{i}^{2}(0)|A^{(n)},S^{(n)},X^{(n)})=o_{P}(n^{2}),$

where we use the facts that (F) holds and

$\displaystyle\beta_{1,s}^{\top}\Gamma_{s}\Gamma_{0,s}^{-1}\Gamma_{s}\beta_{1,s}\leq\left\|\Gamma_{s}^{1/2}\Gamma_{0,s}^{-1}\Gamma_{s}^{1/2}\right\|_{op}||\breve{X}_{\aleph_{s}}\beta_{1,s}||_{2}^{2}=o_{P}(n^{2}).$
(F.4)

This implies $I_{2}/n=o_{P}(1)$.

For $I_{3}$, by Cauchy’s inequality, we have

$\displaystyle|I_{3}|\leq\left\|\beta_{1,s}^{\top}\Gamma_{s}\Gamma_{0,s}^{-1}\breve{X}_{\aleph_{0,s}}^{\top}\right\|_{2}\times
O_{P}(1)=o_{P}(n),$

where the first inequality holds by Assumption 2(iv) and the facts that

$\displaystyle\alpha_{0,s}+\overline{X}_{s}^{\top}\beta_{0,s}-\hat{\tau}_{0}(s)=O_{P}(n^{-1/2}),\quad\sum_{i\in\aleph_{0,s}}e_{i,s}^{2}(a)=o_{P}(1)$

and the equality holds by Assumption 4, (F.4), and the fact that

$\displaystyle\left\|\beta_{1,s}^{\top}\Gamma_{s}\Gamma_{0,s}^{-1}\breve{X}_{\aleph_{0,s}}^{\top}\right\|_{2}=\left(\beta_{1,s}^{\top}\Gamma_{s}\Gamma_{0,s}^{-1}\Gamma_{s}\beta_{1,s}\right)^{1/2}=o_{P}(n).$

We can show $I_{4}=o_{P}(n)$ following the same argument used to bound
$I_{2}$.

For $I_{5}$, we note that the two index sets $\aleph_{0,s}$ and $\aleph_{1,s}$
do not overlap. Therefore, we have

$\displaystyle\mathbb{E}(I_{5}|A^{(n)},S^{(n)},X^{(n)})=0$

and

$\displaystyle
var(I_{5}|A^{(n)},S^{(n)},X^{(n)})\leq\sum_{i\in\aleph_{1,s}}\sum_{j\in\aleph_{0,s}}B_{i,j}^{2}\left[\max_{a=0,1}\max_{i\in\aleph_{a,s}}\mathbb{E}(\varepsilon_{i}^{2}(0)|A^{(n)},S^{(n)},X^{(n)})\right]^{2}=o_{P}(n^{2}),$

where we have
$B_{i,j}=\breve{X}_{i}^{\top}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}\breve{X}_{j}$
and

$\displaystyle\sum_{i\in\aleph_{1,s}}\sum_{j\in\aleph_{0,s}}B_{i,j}^{2}$
$\displaystyle=\sum_{i\in\aleph_{1,s}}\text{trace}(\breve{X}_{i}^{\top}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}\Gamma_{s}\Gamma_{1,s}^{-1}\breve{X}_{i})$
$\displaystyle=\text{trace}(\Gamma_{s}^{1/2}\Gamma_{0,s}^{-1}\Gamma_{s}\Gamma_{1,s}^{-1}\Gamma_{s}^{1/2})$
$\displaystyle=\text{trace}(\Gamma_{s}^{1/2}\Gamma_{0,s}^{-1}\Gamma_{s}^{1/2})+\text{trace}(\Gamma_{s}^{1/2}\Gamma_{1,s}^{-1}\Gamma_{s}^{1/2})=o_{P}(n^{2}),$

where we use the fact that $\Gamma_{s}=\Gamma_{0,s}+\Gamma_{1,s}$.

For $I_{6}$, we have

$\displaystyle|I_{6}|$
$\displaystyle\leq\left|\left[\sum_{i\in\aleph_{1,s}}\breve{X}_{i}\varepsilon_{i}(1)\right]^{\top}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}\left[\sum_{i\in\aleph_{0,s}}\breve{X}_{i}\right]\right||\alpha_{0,s}+\overline{X}_{s}^{\top}\beta_{0,s}-\hat{\tau}_{0}(s)|$
$\displaystyle+\left|\left[\sum_{i\in\aleph_{1,s}}\breve{X}_{i}\varepsilon_{i}(1)\right]^{\top}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}\left[\sum_{i\in\aleph_{0,s}}\breve{X}_{i}e_{i,s}(0)\right]\right|$
$\displaystyle\equiv|I_{6,1}||\alpha_{0,s}+\overline{X}_{s}^{\top}\beta_{0,s}-\hat{\tau}_{0}(s)|+|I_{6,2}|.$

We note that

$\displaystyle\mathbb{E}(I_{6,1}|A^{(n)},S^{(n)},X^{(n)})=0$

and

$\displaystyle var(I_{6,1}|A^{(n)},S^{(n)},X^{(n)})$
$\displaystyle\leq\left[\sum_{i\in\aleph_{0,s}}\breve{X}_{i}\right]^{\top}\Gamma_{0,s}^{-1}\Gamma_{s}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}\left[\sum_{i\in\aleph_{0,s}}\breve{X}_{i}\right]\left[\max_{i\in\aleph_{a,s}}\mathbb{E}(\varepsilon_{i}^{2}(1)|A^{(n)},S^{(n)},X^{(n)})\right]$
$\displaystyle=o_{P}(n^{3}).$

where the second last equality holds because

$\displaystyle\left[\sum_{i\in\aleph_{0,s}}\breve{X}_{i}\right]^{\top}\Gamma_{0,s}^{-1}\Gamma_{s}\Gamma_{1,s}^{-1}\Gamma_{s}S_{0}^{-1}\left[\sum_{i\in\aleph_{0,s}}\breve{X}_{i}\right]$
$\displaystyle\leq
n_{0,s}\left\|\breve{X}_{\aleph_{0,s}}\Gamma_{0,s}^{-1}\Gamma_{s}\Gamma_{1,s}^{-1}\Gamma_{s}S_{0}^{-1}\breve{X}_{\aleph_{0,s}}\right\|_{op}$
$\displaystyle\leq
n_{0,s}\left\|\Gamma_{s}^{1/2}\Gamma_{1,s}^{-1}\Gamma_{s}^{1/2}\right\|_{op}\left\|\breve{X}_{\aleph_{0,s}}\Gamma_{0,s}^{-1}\Gamma_{s}^{1/2}\right\|_{op}^{2}$
$\displaystyle\leq
n_{0,s}\left\|\Gamma_{s}^{1/2}\Gamma_{1,s}^{-1}\Gamma_{s}^{1/2}\right\|_{op}\left\|\Gamma_{s}^{1/2}\Gamma_{0,s}^{-1}\Gamma_{s}^{1/2}\right\|_{op}=o_{P}(n^{3}).$

This implies
$|I_{6,1}||\alpha_{0,s}+\overline{X}_{s}^{\top}\beta_{0,s}-\hat{\tau}_{0}(s)|=o_{P}(n)$
because
$|\alpha_{0,s}+\overline{X}_{s}^{\top}\beta_{0,s}-\hat{\tau}_{0}(s)|=O_{P}(n^{-1/2})$.
Similarly, for $I_{6,2}$, we have

$\displaystyle\mathbb{E}(I_{6,2}|A^{(n)},S^{(n)},X^{(n)})=0$

and

$\displaystyle var(I_{6,2}|A^{(n)},S^{(n)},X^{(n)})$
$\displaystyle\leq\left[\sum_{i\in\aleph_{0,s}}\breve{X}_{i}e_{i,s}(0)\right]^{\top}\Gamma_{0,s}^{-1}\Gamma_{s}\Gamma_{1,s}^{-1}\Gamma_{s}S_{0}^{-1}\left[\sum_{i\in\aleph_{0,s}}\breve{X}_{i}e_{i,s}(0)\right]\left[\max_{i\in\aleph_{a,s}}\mathbb{E}(\varepsilon_{i}^{2}(1)|A^{(n)},S^{(n)},X^{(n)})\right]$
$\displaystyle\leq\left\|\breve{X}_{\aleph_{0,s}}\Gamma_{0,s}^{-1}\Gamma_{s}\Gamma_{1,s}^{-1}\Gamma_{s}S_{0}^{-1}\breve{X}_{\aleph_{0,s}}\right\|_{op}\left\|e_{\aleph_{0,s}}(0)\right\|_{2}^{2}\left[\max_{i\in\aleph_{a,s}}\mathbb{E}(\varepsilon_{i}^{2}(1)|A^{(n)},S^{(n)},X^{(n)})\right]=o_{P}(n^{2}),$

where we use the fact that
$\left\|e_{\aleph_{0,s}}(0)\right\|_{2}^{2}=O_{P}(1)$. This implies
$I_{6}=o_{P}(n)$.

We can show $I_{7}=o_{P}(n)$ and $I_{8}=o_{P}(n)$ following the same argument
used to bound $I_{3}$ and $I_{6}$, respectively.

For $I_{9}$, we have

$\displaystyle|I_{9}|$
$\displaystyle\leq\left\|\breve{X}_{\aleph_{1,s}}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}\breve{X}_{\aleph_{0,s}}\right\|_{op}\left\|1_{n_{1,s}}(\alpha_{1,s}+\overline{X}_{s}^{\top}-\hat{\tau}_{1}(s))+e_{\aleph_{1,s}}(1)\right\|_{2}$
$\displaystyle\times\left\|1_{n_{0,s}}(\alpha_{0,s}+\overline{X}_{s}^{\top}-\hat{\tau}_{0}(s))+e_{\aleph_{0,s}}(0)\right\|_{2}=o_{P}(n),$

where we use the fact that

$\displaystyle\left\|\breve{X}_{\aleph_{1,s}}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}\breve{X}_{\aleph_{0,s}}\right\|_{op}$
$\displaystyle=\left\|\breve{X}_{\aleph_{1,s}}\Gamma_{1,s}^{-1}\Gamma_{s}\Gamma_{0,s}^{-1}\Gamma_{s}\Gamma_{1,s}^{-1}\breve{X}_{\aleph_{1,s}}^{\top}\right\|_{op}^{1/2}$
$\displaystyle\leq\left[\left\|\Gamma_{s}^{1/2}\Gamma_{0,s}^{-1}\Gamma_{s}^{1/2}\right\|_{op}\left\|\Gamma_{s}^{1/2}\Gamma_{1,s}^{-1}\breve{X}_{\aleph_{1,s}}\right\|_{op}^{2}\right]^{1/2}$
$\displaystyle=\left[\left\|\Gamma_{s}^{1/2}\Gamma_{0,s}^{-1}\Gamma_{s}^{1/2}\right\|_{op}\left\|\Gamma_{s}^{1/2}\Gamma_{1,s}^{-1}\Gamma_{s}^{1/2}\right\|_{op}\right]^{1/2}=o_{P}(n)$

and

$\displaystyle\left\|1_{n_{0,s}}(\alpha_{0,s}+\overline{X}_{s}^{\top}-\hat{\tau}_{0}(s))+e_{\aleph_{0,s}}(0)\right\|_{2}=O_{P}(1).$

This leads to the desired result in this step that
$\hat{\beta}_{1,s}^{\top}\Gamma_{s}\hat{\beta}_{0,s}/n_{s}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}cov(\phi_{i}(1),\phi_{i}(0)|S_{i}=s)$.

Step 3: Limit of $\hat{\Sigma}_{\mathcal{W}}$. For
$\hat{\Sigma}_{\mathcal{W}}$, by the proof of Theorem 3.1, we have
$\hat{\tau}_{1,s}-\hat{\tau}_{0,s}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\mathbb{E}(Y(1)-Y(0)|S=s)$ and
$\hat{\tau}^{adj}\stackrel{{\scriptstyle p}}{{\longrightarrow}}\tau$. By
Assumption 1, we also have $\hat{p}_{s}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}p_{s}$, which implies the desired result.

##  Appendix G Proof of Theorem 5.1

Step 1: Asymptotic normality and the expression of $\Omega$. The asymptotic
normality of $\hat{\tau}^{adj}$ has already been established by Ye et al.
(2022) (our $\hat{\tau}^{adj}$ is just their $\hat{\theta}_{A}$) when $k$ is
fixed. In the following, we just sketch the proof of the joint distribution of
$(\hat{\tau}^{adj},\hat{\tau}^{unadj})$ when $k$ is fixed. Following Remark 3,
we have

$\displaystyle\hat{\tau}^{adj}$
$\displaystyle=\frac{1}{n}\sum_{i\in[n]}\frac{A_{i}(Y_{i}-X_{i}^{\top}\hat{\beta}_{1,S_{i}})}{\hat{\pi}_{S_{i}}}-\frac{1}{n}\sum_{i\in[n]}\frac{(1-A_{i})(Y_{i}-X_{i}^{\top}\hat{\beta}_{0,S_{i}})}{1-\hat{\pi}_{S_{i}}}+\frac{1}{n}\sum_{i\in[n]}X_{i}^{\top}(\hat{\beta}_{1,S_{i}}-\hat{\beta}_{0,S_{i}})$
$\displaystyle=\frac{1}{n}\sum_{i\in[n]}\frac{A_{i}(Y_{i}(1)-X_{i}^{\top}\beta_{1,S_{i}}^{*})}{\hat{\pi}_{S_{i}}}-\frac{1}{n}\sum_{i\in[n]}\frac{(1-A_{i})(Y_{i}(0)-X_{i}^{\top}\beta_{0,S_{i}}^{*})}{1-\hat{\pi}_{S_{i}}}+\frac{1}{n}\sum_{i\in[n]}X_{i}^{\top}(\beta_{1,S_{i}}^{*}-\beta_{0,S_{i}}^{*})$
$\displaystyle+\sum_{s\in\mathcal{S}}\hat{p}_{s}\left(\frac{1}{n_{1,s}}\sum_{i\in\aleph_{1,s}}X_{i}-\frac{1}{n_{s}}\sum_{i\in\aleph_{s}}X_{i}\right)^{\top}(\beta_{1,s}^{*}-\hat{\beta}_{1,s})$
$\displaystyle-\sum_{s\in\mathcal{S}}\hat{p}_{s}\left(\frac{1}{n_{0,s}}\sum_{i\in\aleph_{0,s}}X_{i}-\frac{1}{n_{s}}\sum_{i\in\aleph_{s}}X_{i}\right)^{\top}(\beta_{0,s}^{*}-\hat{\beta}_{0,s})$
$\displaystyle=\frac{1}{n}\sum_{i\in[n]}\frac{A_{i}(Y_{i}(1)-X_{i}^{\top}\beta_{1,S_{i}}^{*})}{\hat{\pi}_{S_{i}}}-\frac{1}{n}\sum_{i\in[n]}\frac{(1-A_{i})(Y_{i}(0)-X_{i}^{\top}\beta_{0,S_{i}}^{*})}{1-\hat{\pi}_{S_{i}}}$
$\displaystyle+\frac{1}{n}\sum_{i\in[n]}X_{i}^{\top}(\beta_{1,S_{i}}^{*}-\beta_{0,S_{i}}^{*})+o_{P}(n^{-1/2}),$

where the last equality holds because

$\displaystyle\max_{a=0,1,s\in\mathcal{S}}\left\|\frac{1}{n_{a,s}}\sum_{i\in\aleph_{a,s}}X_{i}-\frac{1}{n_{s}}\sum_{i\in\aleph_{s}}X_{i}\right\|_{2}=O_{P}(n^{-1/2})\quad\text{and}\quad\max_{a=0,1,s\in\mathcal{S}}||\hat{\beta}_{a,s}-\beta_{a,s}^{*}||_{2}=o_{P}(1).$

Let $\delta_{Y,a,s,i}=Y_{i}(a)-\mathbb{E}(Y_{i}(a)|S_{i}=s)$ and
$\delta_{X,s,i}=X_{i}-\mathbb{E}(X_{i}|S_{i}=s)$. Then, we have

$\displaystyle\begin{pmatrix}\sqrt{n}(\hat{\tau}^{adj}-\tau)\\\
\sqrt{n}(\hat{\tau}^{unadj}-\tau)\end{pmatrix}$
$\displaystyle=\sum_{s\in\mathcal{S}}\hat{p}_{s}\sqrt{n}\begin{pmatrix}\frac{1}{n_{1,s}}\sum_{i\in\aleph_{1,s}}(\delta_{Y,1,s,i}-(1-\hat{\pi}_{s})\delta_{X,s,i}^{\top}\beta_{1,s}^{*}-\hat{\pi}_{s}\delta_{X,s,i}^{\top}\beta_{0,s}^{*})\\\
\frac{1}{n_{1,s}}\sum_{i\in\aleph_{1,s}}\delta_{Y,1,s,i}\end{pmatrix}$