# Inference on the Best Policies with Many Covariates

Waverly Wei Division of Biostatistics, University of California, Berkeley.Co-
first authors in alphabetical order Yuqing Zhou 22footnotemark: 2 Department
of Statistics, University of Michigan Zeyu Zheng Department of Industrial
Engineering and Operations Research, University of California, Berkeley
Jingshen Wang 11footnotemark: 1 Correspondence<EMAIL_ADDRESS>
###### Abstract

Understanding the impact of the most effective policies or treatments on a
response variable of interest is desirable in many empirical works in
economics, statistics and other disciplines. Due to the widespread winner’s
curse phenomenon, conventional statistical inference assuming that the top
policies are chosen independent of the random sample may lead to overly
optimistic evaluations of the best policies. In recent years, given the
increased availability of large datasets, such an issue can be further
complicated when researchers include many covariates to estimate the policy or
treatment effects in an attempt to control for potential confounders. In this
manuscript, to simultaneously address the above-mentioned issues, we propose a
resampling-based procedure that not only lifts the winner’s curse in
evaluating the best policies observed in a random sample, but also is robust
to the presence of many covariates. The proposed inference procedure yields
accurate point estimates and valid frequentist confidence intervals that
achieve the exact nominal level as the sample size goes to infinity for
multiple best policy effect sizes. We illustrate the finite-sample performance
of our approach through Monte Carlo experiments and two empirical studies,
evaluating the most effective policies in charitable giving and the most
beneficial group of workers in the National Supported Work program.

JEL codes: C12; C13
Keywords: Winner’s curse; High dimensional data; Linear regression; Order
statistics.

###  1 Introduction

####  1.1 Motivation and our contribution

Many empirical work requires an understanding of the impact of the most
effective policies or treatments on a relevant response variable of interest.
For instance, in randomized (factorial) experiments with multiple treatments,
researchers may be interested in the most effective policies (combinations).
In online platforms, decision makers may be interested in the top five
advertising strategies. In financial portfolio management, managers might want
to learn about the best-performing strategies among many alternatives. In
practice, after different policy effect sizes are estimated from a random
sample, researchers may naturally look into those policies with the largest
effect sizes. Accurately measuring the performance of top policies allows
policy makers to deliver better-informed decisions for forecasting the effects
of future policy implementations.

Nevertheless, given the well-recognized “winner’s curse” phenomenon, there can
be considerable uncertainties concerning if the top policies with large
estimated effect sizes are indeed effective in the population (see Section 1.2
for a literature review). In fact, due to the winner’s curse phenomenon,
literature documents that the estimated effect sizes of the best-performing
policies without additional adjustments tend to be overly optimistic,
rendering under-covered confidence intervals [41, 7]. In this manuscript, we
refer to the optimistic bias introduced by the winner’s curse phenomenon as
the winner’s curse bias. To mitigate this bias issue, we focus on the problem
of constructing accurate point estimates and valid confidence intervals for
the true effect sizes of the (observed) best policies. By the best policies,
we refer to a user-supplied number of policies that have the largest
(estimated) effects among a set of candidate policies (see Section 2.1 for a
concrete problem setup), as we would expect that in practice researchers might
want to focus on a few top policies of interest.

Other than the winner’s curse phenomenon discussed above, an additional
consideration gains prominence in the evaluation of the most effective
policies. Since policy (or intervention) variables are often not exogenous,
researchers may adopt observational methods to estimate their effects. In
recent years, given the increased availability of large datasets with rich
covariate information, one commonly adopted approach in empirical works is to
assume that the policy variables are exogenous after controlling for a
sufficiently large set of factors or covariates. Such a consideration
demandingly requires empirical researchers to estimate the policy effects in
the presence of many covariates.

To simultaneously address the above-mentioned issues, in this article, we
propose a procedure that not only is robust to the presence of many
covariates, but also provides accurate point estimates and valid frequentist
confidence intervals for multiple best policy effect sizes. By many
covariates, we allow the number of covariates $q_{n}$ to diverge with the
sample size $n$ as long as $\limsup_{n\rightarrow\infty}q_{n}/n<1$. Note that
this does not rule out the cases where $q_{n}$ is fixed or $q_{n}=o(n)$. In
other words, our inferential method remains valid when the dimension of the
covariates $q_{n}$ is fixed or $q_{n}=o(n)$. Our proposed confidence intervals
are built upon resampling methods, and we demonstrate that they achieve exact
nominal coverage as the sample size goes to infinity under fairly moderate
assumptions. Our empirical evidence shows that conventional estimates ignoring
the winner’s curse issue are substantially upward biased, while our
corrections reduce the winner’s curse bias and increase coverage. As far as we
know, valid statistical inferential tools on multiple best policies that lift
the winner’s curse while incorporating possibly many covariates have been
lacking, and the contribution of our work is to bridge this gap and help
policy makers deliver well-informed decisions in practice.

We illustrate our method with two empirical applications. In the first case
study, we use the charitable giving data from [38] to evaluate the best
pricing policies that motivate donors to give. Our results suggest that simple
methods without adjusting for the winner’s curse bias could be potentially
overly optimistic in identifying the most effective polices. After accounting
for the winner’s curse bias, we do not find sufficient evidence to support
that the second best pricing policy–asking the donor to give 25% more than
his/her highest historical donation–is effective, implying that asking for a
more “expensive” donation may not encourage donors to give. We nevertheless
note that given our calibration only marginally reduces the effect size of the
second best policy, the above conclusion might not warrant a different
economic interpretation. In the second case study, we evaluate the
effectiveness of the national supported work (NSW) program in different groups
of workers. The NSW program is a job training program designed to prepare
disadvantaged workers for employment, and it has been investigated in various
studies [19, 40]. We apply the proposed approach to evaluate the performance
of the NSW program on the most-affected subgroups of workers observed in the
dataset. Our study results potentially suggest that married black workers
might benefit from the NSW program with an average increase of $4,410 for
their annual income.

####  1.2 Connection to the existing literature

One fundamental trend that drives the motivation of the methodology developed
in this manuscript is the increasing availability of massive datasets and the
associated increasing dimensionality. Such a trend brings scientists
opportunities to deliver better-informed policies but, at the same time,
presents challenges in developing econometric and statistical tools; see [28],
[9], [25], [12] for example. A recent book [27] provides a thorough discussion
of analytical methods that aim to address such challenges. Specifically, the
increasing data availability brings challenges and also opportunities to
better understand various policies whose effects can be inferred from data.
Along this line, our manuscript aims at providing understating for policies
that are estimated and selected to be the most effective from a pool of
policies.

The winner’s curse phenomenon and its related issues have been widely
recognized in economics, statistics, and data science at large. Seminal works
by [30, 29] point out that spurious discoveries can easily arise when target
parameters are selected through data mining and statistical machine learning
algorithms. Recent work by [7] considers performing conditional and
unconditional inference on observed best policy and [6] extends the work to
more general ranking problems, which is still different from our goal in
conducting unconditional inference on multiple top policies. Moreover, while
the conditional approaches in [7] and [6] produce optimal confidence intervals
for the observed policy effects, their point estimates and confidence
intervals can be conservative when they are applied unconditionally. [20]
considers a method to handle the winner’s curse bias with Tweedie’s formula
concerning the empirical Bayes theory. [41] consider a plug-in correction of
the winner’s curse bias and propose to construct confidence interval based on
bootstrapping in the context of A/B testing, but the proposed method lacks
theoretical justifications. In clinical trials for evaluating the largest
observed treatment effect in multiple subpopulations, [33] propose a
bootstrap-based confidence interval that achieves the exact nominal level as
the sample size goes to infinity, though generalizing their method to make
inference on several top policies might not be straightforward, especially in
the presence of many covariates.

Our manuscript builds upon the literature on linear regression models with
many or high dimensional covariates; see [35], [44], [43], [3], [22], [11],
[10], [37] and the reference therein. In particular, [44] has established the
asymptotic normality results for any contrasts of the ordinary least squares
(OLS) coefficient vector estimator, when the dimension of the covariates
divided by the sample size vanishes asymptotically. More recently, [12] have
shown that a small subset of the OLS estimators for the regression
coefficients are asymptotically normal without restricting the dimension of
the covariates to be a vanishing fraction of the sample size. Moreover, [12]
have proposed a robust covariance matrix estimator for the subset of the the
OLS estimator under fairly general conditions. [37] has proposed an
alternative covariance matrix estimator that can deal with designs with even
large number of covariates under additional assumptions (Assumption 4 in the
current manuscript).

Making inference on the best-performing policies is related to the literature
on constructing confidence intervals for extrema parameters with bootstrap;
see [5], [24], [52], [15], [17] and the reference therein. Given the
asymptotic distributions of extrema parameter estimators are often not normal,
bootstrap-based methods can face serious difficulties when used to replicate
the distribution of extrema of parameter estimators [44, 45]. While
subsampling could overcome this issue faced by the classical bootstrap, it can
exhibit very poor finite-sample performance because of the noise introduced by
the vanishing subsample size. Different from our goal in constructing
confidence intervals that achieve the exact nominal level, Hall and Miller
[34] and [15] propose to construct conservative bootstrap confidence intervals
for extrema of parameters. In our current problem setup with many covariates,
the problem becomes even more acute as [23] show through a mix of simulation
and theoretical analyses that the bootstrap is fraught with problems in
moderate high dimensions. In the context of meta-analyses, [17] propose an
approach to make inference on ordered fixed study-specific parameters when
different parameters are estimated independently from multiple studies.

Our method also contributes to the rapidly growing literature on program
evaluations; see [31], [9], [39], [8], [2], [14], [26], [48] among many
others. Under our asymptotic regime where the number of covariates $q_{n}$
grows with the sample size $n$, the Neyman orthogonalization based approaches
often need to work with models with sparse regression coefficients [9]. Rather
than imposing such a sparsity assumption, our approach estimates the policy
effects with regression adjustments without requiring the regression
coefficients to be sparse. Because our approach only requires a consistent
covariance matrix estimation for different policy effect estimators, we expect
that the proposed framework on evaluating the best policies can be generalized
when different policy effects are estimated with other off-shelf methods and
we relegate such extensions for future work.

Notation.  We work with triangular array data
$\\{\omega_{i,n}:i=1,\ldots,n;n=1,2,\ldots\\}$ where for each $n$,
$\\{\omega_{i,n}:i=1,\ldots,n\\}$ is defined on the probability space
$(\Omega,\mathcal{S},P_{n})$. All parameters that characterize the
distribution of $\\{\omega_{i,n}:i=1,\ldots,n\\}$ are implicitly indexed by
$P_{n}$ and thus by $n$. We write vectors and matrices in bold font, and use
regular font for univariate variables and constants.

###  2 Model setup and methodology

####  2.1 Problem setup and a revisit to the winner’s curse phenomenon

Suppose we have a random sample
$\\{(y_{i,n},\bm{x}_{i,n}^{\prime},\bm{w}_{i,n}^{\prime})^{\prime}\\}_{i=1}^{n}$,
we pose the problem in the framework of a linear regression model under
heteroscedasticity

$\displaystyle
y_{i,n}=\bm{x}_{i,n}^{\prime}\bm{\beta}+\bm{w}_{i,n}^{\prime}\bm{\gamma}_{n}+u_{i,n},\quad
i=1,\ldots,n,$ (1)

where $y_{i,n}$ is the outcome variable, $\bm{x}_{i,n}\in\mathbb{R}^{d}$ are
the treatment or policy variables of interest,
$\bm{w}_{i,n}\in\mathbb{R}^{q_{n}}$ contains the confounding factors,
$u_{i,n}$ is an unobserved error term, and the coefficient vector
$\bm{\beta}=(\beta_{1},\ldots,\beta_{d})^{\prime}$ contains the treatment
effect of $\bm{x}_{i,n}$ on the outcome $y_{i,n}$. We allow the linear model
(1) to hold approximately by allowing
$\mathbb{E}[u_{i,n}|\\{\bm{x}_{i,n}\\}_{i=1}^{n},\\{\bm{w}_{i,n}\\}_{i=1}^{n}]\neq
0$. We are also in a scenario where $\bm{w}_{i,n}$ is high-dimensional, in the
sense that $q_{n}$ can be a vanishing fraction of the sample size $n$ as long
as $\limsup_{n\rightarrow\infty}q_{n}/n<1$. To simplify notations, we drop
subscript $n$ in univariate random variables in the rest of the manuscript.
That is, for example, we denote $\gamma_{j,n}$ as $\gamma_{j}$.

We write the ordered values of $\beta_{1},\ldots,\beta_{d}$ as
$\beta_{(1)}\geq\ldots\geq\beta_{(d)}$. We adopt the ordinary least-squares
(OLS) estimator $\widehat{\bm{\beta}}$ (see Remark 2 for other possible
estimates) to estimate $\bm{\beta}$ and write the order statistics of
$\widehat{\bm{\beta}}$ as
$\widehat{\beta}_{(1)}\geq\ldots\geq\widehat{\beta}_{(d)}$. Because
researchers in practice might hope to focus on a few top policies, given that
$d_{0}$ is a user-supplied positive integer, our goal is to construct accurate
point estimates and valid confidence intervals for two sets of quantities:

  1. (1)

the best policy effect sizes in the population:
$\beta_{(1)},\ldots,\beta_{(d_{0})}$,

  2. (2)

the observed best policy effect sizes: $\beta_{\widehat{j}}$, where
$\widehat{j}=\sum_{k=1}^{d}k\cdot\mathds{1}(\widehat{\beta}_{k}=\widehat{\beta}_{(j)})$,
for $j=1,\ldots,d_{0}$.

The first set of quantities characterizes the effects of the top $d_{0}$
policies in the population and are thus fixed parameters. The second set of
quantities describes the true effect sizes of the best performing policies
observed in the random sample, and these quantities are thus “data-dependent
parameters.” Both sets of quantities can be of interest in different empirical
applications [18, 16, 46], and our proposed procedure can be used to deliver
valid statistical inference on both quantities (Theorem 1 and Corollary 1).

###### Remark 1 (Ties in the estimated policy effects)

The second set of parameters is well defined if the observed policies do not
have exact ties in the sense that
$\widehat{\beta}_{(1)}>\ldots>\widehat{\beta}_{(d)}$. When the policies effect
estimators solve to the interior points of the feasible parameter space, it is
likely that no exact ties appear in the random sample. On the other hand,
there can exist scenarios where, for example, $d_{0}$ is set as 2 but there
are multiple policy effect sizes that tie at rank 2. In this case, one may
choose instead a data-dependent
$\widehat{d}_{0}=\max\\{k:\widehat{\beta}_{(k)}>\widehat{\beta}_{(2)}-C_{1}\cdot
n^{-0.25}\\}$. This new random $d_{0}$ will asymptotically be able to
incorporate all the effect sizes that actually are equal to the true effect
size associated with $\widehat{\beta}_{(2)}$. In this way, the limiting value
of $\widehat{d}_{0}$ will not necessarily be 2, but can be a larger number
than 2 to incorporate “very close” effect sizes with the rank-2 effect size.
We also provide some related discussions in Remark 4.

###### Remark 2 (Other possible estimators of $\bm{\beta}$)

In the presence of many covariates when $q_{n}$ is potentially large
($\limsup_{n\rightarrow\infty}q_{n}/n\rightarrow 1$ in our asymptotic regime)
without assuming the coefficient $\mathbf{\gamma}_{n}$ to be sparse, we adopt
the OLS estimator to estimate $\bm{\beta}$, because the OLS estimator has been
thoroughly studied in the existing literature and enjoys favorable theoretical
guarantees. In high dimensions when $q_{n}\gg n$, other estimators of
$\bm{\beta}$ that incorporate model selection procedures can be adopted as
well. Our procedure can produce valid statistical inference as long as the
covariance matrix of $\widehat{\bm{\beta}}$ can be consistently estimated. For
example, under the sparsity assumption on $\bm{\gamma}_{n}$ documented in the
literature [27], we may adopt the covariance matrix estimator from the de-
sparsified Lasso procedure [32, 53].

To fully realize the challenges on delivering valid statistical inference on
these two sets of parameters in our current problem setup, we revisit the
winner’s curse phenomenon. When first discussed in common-value auctions, the
winner’s curse refers to the bidding behavior where bidders systematically
overbid, resulting in an expected loss [13]. In our context of policy
evaluations, the winner’s curse refers to the issue that the observed best
policies have the tendency to over-estimate the best policies in the
population. We would thus often expect that neither
$\mathbb{E}[\widehat{\beta}_{(j)}-\beta_{(j)}]$ nor
$\mathbb{E}[\widehat{\beta}_{(j)}-\beta_{\widehat{j}}]$ is close to zero, and
the resulting confidence interval may fail to reach the nominal level. Such an
issue becomes even more acute as we have many covariates $\bm{w}_{i,n}$
entering the inferential process.

We next illustrate the winner’s curse issue through Example 1 with a simple
simulation study, where we observe substantial winner’s curse bias and under-
covered confidence intervals for the top polices. In particular, Figure 1(b)
demonstrates that coverage probabilities are worsened when a larger number of
covariates are incorporated for estimating $\bm{\beta}$. It is worth pointing
out that when $d=3$, $\widehat{\beta}_{(2)}$ is the median policy effect.
Thus, the estimation bias is around 0, and the true standard deviation is much
smaller than the estimated standard deviation, resulting in a confidence
interval with close to 100% coverage. When $d$ increases, the coverage
probability gradually drops due to a larger estimation bias and inaccurately
estimated standard deviation.

###### Example 1 (A simulation study demonstrating the winner’s curse
phenomenon with many covariates)

We generate 1000 Monte Carlo samples following the setup in Model (1). We
generate $\bm{x}_{i,n}\sim\mathcal{N}(0,\bm{\Sigma})$ with
$\Sigma_{jk}=0.5^{|j-k|}$ for $j,k=1,\ldots,d$,
$\bm{w}_{i,n}=\mathds{1}(\widetilde{\bm{w}}_{i,n}\geq\Phi^{-1}(0.98))$ with
$\widetilde{\bm{w}}_{i,n}\sim\mathcal{N}(0,\bm{I}_{q_{n}})$, where
$\bm{I}_{q_{n}}$ is a $q_{n}$-dimensional identity matrix. We consider the
case where no policy is effective (so that $\bm{\beta}=0$,
$\beta_{\widehat{1}}=\beta_{\widehat{2}}=0$) and $\gamma_{j}=1/j$, for
$j=1,\ldots,q_{n}$. We report the asymptotic bias of the conventional
estimator (i.e., $\sqrt{n}\cdot\mathbb{E}[\widehat{\beta}_{(j)}-\beta_{(j)}]$)
as well as the coverage probability of confidence intervals constructed based
on normal approximation with the Eicker-White [21, 51] covariance matrix
estimator defined in Eq (6).

(a) Winner’s curse bias ($q_{n}=141$)

(b) Coverage probability

Figure 1: Demonstration of the winner’s curse phenomenon following the
simulation setup in Example 1. The maximum Monte Carlo standard error for the
asymptotic bias is 0.88. Panel (a) captures the asymptotic winner’s curse bias
when $q_{n}=141$; Panel (b) captures the coverage probability when
$q_{n}\in\\{141,631\\}$ and the nominal level is 0.95.

####  2.2 Methodology

Our method starts with the ordinary least-squares (OLS) estimator of
$\bm{\beta}$, that is

$\displaystyle\widehat{\bm{\beta}}=\Big{(}\sum_{i=1}^{n}\widehat{\bm{v}}_{i,n}\widehat{\bm{v}}_{i,n}^{\prime}\Big{)}^{-1}\Big{(}\sum_{i=1}^{n}\widehat{\bm{v}}_{i,n}y_{i,n}\Big{)},$

where $\widehat{\bm{v}}_{i,n}=\sum_{j=1}^{n}(\bm{M}_{n})_{i,j}\bm{x}_{j,n}$,
and
$(\bm{M}_{n})_{i,j}\triangleq\mathds{1}(i=j)-\bm{w}^{\prime}_{i,n}\Big{(}\sum_{k=1}^{n}\bm{w}_{k,n}\bm{w}_{k,n}^{\prime}\Big{)}^{-1}\bm{w}_{j,n}$.
As we focus on the case when $q_{n}$ can be a non-vanishing fraction of $n$,
$n\rightarrow\infty$, we adopt the robust covariance matrix estimator proposed
in [37]. We try to follow the author’s notation as closely as possible:

$\displaystyle\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}\triangleq\widehat{\bm{\Gamma}}_{n}^{-1}\widehat{\Sigma}^{\texttt{KJ}}_{n}\widehat{\bm{\Gamma}}_{n}^{-1},$

where

$\displaystyle\widehat{\bm{\Gamma}}_{n}=\frac{1}{n}\sum_{i=1}^{n}\widehat{\bm{v}}_{i,n}\widehat{\bm{v}}_{i,n}^{\prime},\quad\widehat{\Sigma}^{\texttt{KJ}}_{n}\triangleq\frac{1}{n}\sum_{i=1}^{n}\widehat{\bm{v}}_{i,n}\widehat{\bm{v}}_{i,n}^{\prime}y_{i,n}\acute{u}_{i,n},$

where $\acute{u}_{i,n}=\frac{\widehat{u}_{i,n}}{(\bm{M}_{n})_{i,i}}$,
$\widehat{u}_{i,n}=\sum_{j=1}^{n}(\bm{M}_{n})_{i,j}(y_{j,n}-\bm{x}_{j,n}^{\prime}\widehat{\bm{\beta}})$,
for $i=1,\ldots,n$. Such an estimator is well-defined as long as
$\min_{i}(\bm{M}_{n})_{i,i}>0$. If $\min_{i}(\bm{M}_{n})_{i,i}=0$, it means
that the auxiliary regression produces a perfect prediction. So the
observation does not carry information on $\bm{\beta}$ and can be ignored.

As shown in Example 1, the estimated top policy effect sizes with
$\widehat{\beta}_{(1)},\ldots,\widehat{\beta}_{(d_{0})}$ are often biased
upward for our target parameters due to the winner’s curse phenomenon.
Inspired by the procedure proposed by [17]111Note that there is a typo in [17]
for the definition of the near tie set. Although the near tie
$\mathcal{H}_{(j)}$ in their manuscript was originally defined as
$\mathcal{H}_{(j)}=\big{\\{}k:|\beta_{k}-\beta_{(j)}|=O(n^{-\frac{1}{2}}),\
k=1,\ldots,d\big{\\}}$, their proof goes through when the near tie set is
defined with
$\mathcal{H}_{(j)}=\big{\\{}k:|\beta_{k}-\beta_{(j)}|=o(n^{-\frac{1}{2}}),\
k=1,\ldots,d\big{\\}}$.  for meta-analyses, we generate replicates of
$\widehat{\bm{\beta}}$ from a multivariate normal distribution

$\displaystyle\widehat{\bm{\bm{\beta}}}^{*}\big{|}\\{(y_{i,n},\bm{x}_{i,n}^{\prime},\bm{w}_{i,n}^{\prime})^{\prime}\\}_{i=1}^{n}\sim\mathcal{N}(\widehat{\bm{\beta}},\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}/n),\quad\text{
where
}\widehat{\bm{\bm{\beta}}}^{*}=(\widehat{\beta}^{*}_{1},\ldots,\widehat{\beta}^{*}_{d})^{\prime},$
(2)

and we denote the ordered values of the vector $\widehat{\bm{\bm{\beta}}}^{*}$
as $\widehat{\beta}^{*}_{(1)}\geq\ldots\geq\widehat{\beta}^{*}_{(d)}$. Note
that the above description of $\widehat{\bm{\bm{\beta}}}^{*}$ differs from
some previous work on bootstrapping insofar we have suppressed the role of
“multiplier variables,” and we have defined $\widehat{\bm{\bm{\beta}}}^{*}$ as
a sample from
$\mathcal{N}(\widehat{\bm{\beta}},\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}/n)$.
Different from [17] that requires different estimators to be estimated from
independent studies with non-overlapping random samples, our approach relaxes
such an requirement and allows
$\widehat{\beta}_{1},\ldots,\widehat{\beta}_{d}$ to be correlated.

Next, given properly chosen $b_{\texttt{L}}$ and $b_{\texttt{R}}$ so that
$b_{\texttt{L}}-b_{\texttt{R}}=O(n^{-\delta})$ with $\delta\in(0,\frac{1}{2})$
(see Supplementary Materials Section C.1 for their data-adaptive choices, and
robustness to different choices of tuning parameters in Supplementary
Materials, Section C.2), we estimate a “near tie” set that captures policies
that have similar effect sizes to the $j$-th largest policy:

$\displaystyle\widehat{\mathcal{H}}_{(j)}=\big{\\{}k:\widehat{\beta}^{*}_{{{(j)}}}-b_{\texttt{L}}\leq\widehat{\beta}^{*}_{k}\leq\widehat{\beta}^{*}_{{{(j)}}}+b_{\texttt{R}},\
k=1,\ldots,d\big{\\}}.$

We then record the averages of
$\widehat{\beta}^{*}_{1},\ldots,\widehat{\beta}^{*}_{d}$ and of
$\widehat{\beta}_{1},\ldots,\widehat{\beta}_{d}$ in the estimated tie set
$\widehat{\mathcal{H}}_{(j)}$ as

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\beta}_{(j)}^{*}=\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}^{*}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}},\text{
and
}\widetilde{\beta}_{(j)}=\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|},$
(3)

where $|\widehat{\mathcal{H}}_{(j)}|$ denotes the cardinality of the set
$\widehat{\mathcal{H}}_{(j)}$.

Finally, we apply the above resampling procedure to construct point estimates
and confidence intervals for $\beta_{(j)}$ as well as $\beta_{\widehat{j}}$
(as defined in Section 2.1 and in Eq (5)), $j=1,\ldots,d_{0}$. Specifically,
for confidence interval construction, we generate $B$ independent samples of
$\widetilde{\beta}^{*}_{(j)}$ as in Eq (3), and then define
$\widehat{q}_{(j)}(\alpha/2)$ to be the empirical $\alpha/2$-quantile of the
$B\geq 1$ samples (and similarly for $\widehat{q}_{(j)}(1-\alpha/2)$), leading
to a level-$\alpha$ confidence interval for $\beta_{(j)}$ with

$\displaystyle\big{[}\widehat{q}_{(j)}(\alpha/2),\
\widehat{q}_{(j)}(1-\alpha/2)\big{]},\quad j=1,\ldots,d_{0}.$

Corollary 1 demonstrates that the above confidence interval also serves as an
asymptotically exact level-$\alpha$ prediction interval for
$\beta_{\widehat{j}}$. For point estimates, we may either use
$\widetilde{\beta}_{(j)}$ or the averaged resampled statistics
$\widetilde{\beta}_{(j)}^{*}$ to estimate $\beta_{(j)}$ and
$\beta_{\widehat{j}}$.

###  3 Theoretical investigation

####  3.1 Notations and assumptions

Before discussing the theoretical results in detail, we revisit and introduce
some notations and assumptions adopted in the manuscript. We denote the sample
$\\{(y_{i,n},\bm{x}_{i,n}^{\prime},\bm{w}_{i,n}^{\prime})^{\prime}\\}_{i=1}^{n}$
as $\\{\bm{z}_{i,n}\\}_{i=1}^{n}$. Recall $u_{i,n}$ is the random error in the
considered linear model (1), we define

$\displaystyle\varepsilon_{i,n}=u_{i,n}-\mathbb{E}[u_{i,n}|\\{\bm{w}_{i,n}\\}_{i=1}^{n},\\{\bm{x}_{i,n}\\}_{i=1}^{n}],\quad\bm{v}_{i,n}=\bm{x}_{i,n}-\mathbb{E}[\bm{x}_{i,n}|\\{\bm{w}_{i,n}\\}_{i=1}^{n}],\quad
i=1,\ldots,n.$ (4)

Let
$e_{i,n}=\mathbb{E}[u_{i,n}|\\{\bm{w}_{i,n}\\}_{i=1}^{n},\\{\bm{x}_{i,n}\\}_{i=1}^{n}]$,
we further denote

$\displaystyle\sigma_{i,n}^{2}=\mathbb{E}[\varepsilon_{i,n}^{2}|\\{\bm{w}_{i,n}\\}_{i=1}^{n},\\{\bm{x}_{i,n}\\}_{i=1}^{n}],\quad\widetilde{\bm{v}}_{i,n}=\sum_{j=1}^{n}(\bm{M}_{n})_{i,j}\bm{v}_{j,n},$
$\displaystyle\rho_{n}^{1}=\frac{\sum_{i=1}^{n}\mathbb{E}\big{[}e_{i,n}^{2}\big{]}}{n},\quad\rho_{n}^{2}=\frac{\sum_{i=1}^{n}\mathbb{E}\Big{[}\mathbb{E}\big{(}e_{i,n}|\\{\bm{w}_{i,n}\\}_{i=1}^{n}\big{)}^{2}\Big{]}}{n},$
$\displaystyle\bm{Q}_{i,n}=\mathbb{E}\Big{[}\bm{x}_{i,n}-\big{(}\sum_{j=1}^{n}\mathbb{E}[\bm{w}_{j,n}\bm{w}_{j,n}^{\prime}]\big{)}^{-1}\sum_{j=1}^{n}\mathbb{E}[\bm{w}_{j,n}\bm{x}^{\prime}_{j,n}]\Big{|}\\{\bm{w}_{i,n}\\}_{i=1}^{n}\Big{]},$
$\displaystyle\widetilde{\bm{Q}}_{i,n}=\sum_{j=1}^{n}(\bm{M}_{n})_{i,j}\bm{Q}_{i,n}.$

For a policy $j$, we define the near tie set in the population as:

$\displaystyle\mathcal{H}_{(j)}=\big{\\{}k:|\beta_{k}-\beta_{(j)}|=o(n^{-\frac{1}{2}}),\
k=1,\ldots,d\big{\\}}.$

Next, let $\widehat{\bm{e}}_{j}$ denote a $d$-dimensional (sparse) vector with

$\displaystyle\widehat{\bm{e}}_{j}=(\widehat{e}_{j,1},\ldots,\widehat{e}_{j,d}),\quad\widehat{e}_{jk}=\frac{\mathds{1}(k\in\widehat{\mathcal{H}}_{(j)})}{|\widehat{\mathcal{H}}_{(j)}|},\quad
k=1,\ldots,d.$

We will use the notation $\mathbb{P}(\cdot|\\{\bm{z}_{i,n}\\}_{i=1}^{n})$ to
refer to the probability that is conditional on the random variables
$\\{\bm{z}_{i,n}\\}_{i=1}^{n}$.

We make following assumptions throughout this section. Note that Assumptions
1-4 listed below largely follow the assumptions in [12] and [37], we list
these assumptions along with their interpretations to present a full picture
for our readers.

###### Assumption 1 (Sampling)

The errors $\varepsilon_{i,n}$ are uncorrelated across $i$ conditional on
$\\{\bm{x}_{i,n}\\}_{i=1}^{n}$ and $\\{\bm{w}_{i,n}\\}_{i=1}^{n}$. Let
$\\{N_{1},\ldots,N_{G_{n}}\\}$ represents a partition of $\\{1,\ldots,n\\}$
with $\underset{1\leq g\leq G_{n}}{\max}{|N_{g}|}=O(1)$ such that
$\\{(\varepsilon_{i,n},\bm{v}_{i,n}),i\in N_{g}\\}$ (defined in (4)) are
independent across $g$ conditional on $\\{\bm{w}_{i,n}\\}_{i=1}^{n}$.

Assumption 1 generalizes the classical independent and identically distributed
(i.i.d.) setting to allow for repeated measurements or group structures in the
observed data. For example, Assumption 1 allows the observed data to form
clusters of finite sample sizes, and within-cluster dependency is allowed as
long as the observations between clusters are independent.

###### Assumption 2 (Design)

The dimension of the covariates $\bm{w}_{i,n}$ satisfies that
$\limsup_{n\rightarrow\infty}q_{n}/n<1$. The minimum eigenvalue of the matrix
$\sum_{i=1}^{n}\bm{w}_{i,n}\bm{w}_{i,n}^{\prime}$ is bounded away from 0 with
probability approaching one, that is

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}\lambda_{\min}\big{(}\sum_{i=1}^{n}\bm{w}_{i,n}\bm{w}_{i,n}^{\prime}\big{)}>0\Big{)}=1.$

Lastly,

$\displaystyle\underset{1\leq i\leq
n}{\max}\Big{\\{}\mathbb{E}[\varepsilon_{i,n}^{4}|\\{\bm{w}_{i,n}\\}_{i=1}^{n},\\{\bm{x}_{i,n}\\}_{i=1}^{n}],\
\frac{1}{\sigma_{i,n}^{2}},$
$\displaystyle\quad\quad\quad\quad\mathbb{E}[\bm{v}_{i,n}^{4}|\\{\bm{w}_{i,n}\\}_{i=1}^{n}],\
1/\lambda_{\min}\big{(}\frac{\sum_{i=1}^{n}\mathbb{E}\big{[}\widetilde{\bm{v}}_{i,n}\widetilde{\bm{v}}_{i,n}^{\prime}|\\{\bm{w}_{i,n}\\}_{i=1}^{n}\big{]}}{n}\big{)}\Big{\\}}=O_{p}(1).$

Assumption 2 contains three conditions. The first condition allows the
dimension of the covariates $\bm{w}_{i,n}$ to grow at the sample rate as the
sample size $n$. The second condition requires the matrix
$\sum_{i=1}^{n}\bm{w}_{i,n}\bm{w}_{i,n}^{\prime}$ to be full rank, which is
necessary otherwise the OLS estimator would not be able to calculate the
matrix $\bm{M}_{n}$. Furthermore, as noted in [12], such an assumption can be
imposed by dropping any covariates in $\bm{w}_{i,n}$ that are collinear. The
third condition contains conventional moment conditions for the covariates and
heteroscedasticity.

###### Assumption 3 (Linear model approximation)

$\sum_{i=1}^{n}\mathbb{E}[||\bm{Q}_{i,n}||^{2}]/n=O(1)$,
$\rho_{n}^{1}+n(\rho_{n}^{1}-\rho_{n}^{2})+\rho_{n}^{1}\cdot\sum_{i=1}^{n}\mathbb{E}[||\bm{Q}_{i,n}||^{2}]=o(1)$,
and $\underset{1\leq i\leq
n}{\max}||\widehat{\bm{v}}_{i,n}||/\sqrt{n}=o_{p}(1)$, $n\rho_{n}^{1}=O(1)$.

Assumption 3 mainly characterizes the difference between the mean squares of
the conditional errors $\rho_{n}^{1}$ and the projection $\rho_{n}^{2}$ into
the covariate space $\\{\bm{w}_{i,n}\\}$’s. The characterization of this
difference involves $\sum_{i=1}^{n}\mathbb{E}[||\bm{Q}_{i,n}||^{2}]$ where
$\bm{Q}_{i,n}$ describes the deviation of $\bm{x}_{i,n}$ from its population
linear projection. Residuals of this linear projection, represented by
$\widehat{\bm{v}}_{i,n}$’s, are assumed to satisfy a negligibility condition
after a maximization over all $i$’s. This negligibility condition regularizes
the distributional connection between $\bm{x}_{i,n}$’s and $\bm{w}_{i,n}$’s.
We note that if the mean squares of $\bm{x}_{i,n}$’s are bounded and that an
exogeneity condition $e_{i,n}=0$ holds for all $i$ and $n$, then the linear
model approximation assumption naturally holds. Otherwise, if the exogeneity
condition does not hold, Assumption 3 requires a small-bias condition
$n\rho_{n}^{1}=O(1)$.

###### Assumption 4 (Variance estimation)

$\lim_{n\rightarrow\infty}\mathbb{P}(\min_{i}\big{(}\bm{M}_{n}\big{)}_{i,i}>0)=1$,

$\displaystyle\mathbb{P}\Big{(}\min_{i}\big{(}\bm{M}_{n}\big{)}_{i,i}>0\Big{)}=O_{p}(1),\quad\frac{\sum_{i=1}^{n}||\widetilde{\bm{Q}}_{i,n}||^{4}}{n}=O_{p}(1),$

and $\max_{i}||\mu_{i,n}||/\sqrt{n}=o_{p}(1)$ with
$\mu_{i,n}=\mathbb{E}\big{[}y_{i,n}|\\{\bm{x}_{i,n}\\}_{i=1}^{n},\\{\bm{w}_{i,n}\\}_{i=1}^{n}\big{]}$.

Assumption 4 has two major parts. The first part regularizes the diagonal
elements $(\bm{M}_{n}\big{)}_{i,i}$’s, essentially requiring the smallest
diagonal element to be consistently bounded away from zero when $n$ tends to
infinity. Even though it is difficult to provide broadly general primitives to
validate this assumption, Assumption 2 of [12], Assumption 4 of [37], and the
discussions therein provide sufficient conditions for this assumption to hold.
The second part regularizes $\mu_{i,n}$’s and $\widetilde{\mathbf{Q}}_{i,n}$’s
in order to control the variance of $y_{i,n}$’s and the variance of
$\mathbb{E}(\bm{v}_{i,n}|\\{\bm{w}_{i,n}\\}_{i=1}^{n})$’s.

###### Assumption 5 (Policy effect sizes)

For $\delta\in(0,\frac{1}{2})$, the asymptotic distance between the effects of
policy $k\not\in\mathcal{H}_{(j)}$ and $j\in\mathcal{H}_{(j)}$ diverges as
$n\rightarrow\infty$:

$\displaystyle
n^{\delta}\cdot\min_{k\not\in\mathcal{H}_{(j)}}\big{|}\beta_{(j)}-\beta_{k}\big{|}\rightarrow\infty,\text{
as }n\rightarrow\infty,\quad j=1,\ldots,d.$

Assumption 5 requires that any policies outside the near tie set
$\mathcal{H}_{(j)}$ have effect sizes sufficiently different from the ones in
$\mathcal{H}_{(j)}$. In fixed dimensions when $q_{n}$ does not grow with $n$,
the underlying policy effect sizes $\beta_{1},\ldots,\beta_{d}$ are constant
with respect to the sample size $n$. The near tie set reduces to a “precise”
tie set $\mathcal{H}_{(j)}=\big{\\{}k:\beta_{k}=\beta_{(j)},\
k=1,\ldots,d\big{\\}}$, suggesting that
$\min_{k\not\in\mathcal{H}_{(j)}}\big{|}\beta_{(j)}-\beta_{k}\big{|}$ is a
positive constant bounded away from zero. In such a case, Assumption 5 is
automatically satisfied.

####  3.2 Properties of the proposed estimator

For the proposed estimator, we show that the following theorem holds:

###### Theorem 1

Under Assumptions 1-5, for any $t\in\mathbb{R}$, for the resampled statistics,
the following holds

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\left(\frac{\sqrt{n}\big{(}\widetilde{\beta}_{(j)}^{*}-\widetilde{\beta}_{(j)}\big{)}}{(\widehat{\bm{e}}_{j}^{\prime}\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}\widehat{\bm{e}}_{j})^{\frac{1}{2}}}\leq
t\Big{|}\\{(y_{i,n},\bm{x}_{i,n}^{\prime},\bm{w}_{i,n}^{\prime})^{\prime}\\}_{i=1}^{n}\right)=\Phi(t).$

For the original statistics, it holds that

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\left(\frac{\sqrt{n}\big{(}\widetilde{\beta}_{(j)}-{\beta}_{(j)}\big{)}}{(\widehat{\bm{e}}_{j}^{\prime}\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}\widehat{\bm{e}}_{j})^{\frac{1}{2}}}\leq
t\right)=\Phi(t).$

Furthermore, we have that
$\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}\mathbb{P}\big{(}\widetilde{\beta}_{(j)}^{*}\leq\beta_{(j)}|\\{\bm{z}_{i,n}\\}_{i=1}^{n}\big{)}\leq
s\Big{)}=s$.

Theorem 1 confirms that our proposed confidence interval for $\beta_{(j)}$
achieves exact $1-\alpha$ coverage probability as the sample size goes to
infinity when $B$ is sufficiently large, which distinguishes the proposed
inference procedure from simultaneous methods. Furthermore, Theorem 1 says
that $\widetilde{\beta}_{(j)}$ is a root-$n$ consistent estimator of
$\beta_{(j)}$, in the sense that $\forall\varepsilon>0$, there exists $M>0$
such that
$\mathbb{P}\big{(}|\sqrt{n}(\widetilde{\beta}_{(j)}-{\beta}_{(j)})|>M\big{)}\leq\varepsilon$,
for $n\geq 1$.

As for the observed best policies, recall that we denote the observed $j$-th
largest policy as

$\displaystyle\widehat{j}=\sum_{k=1}^{d}k\cdot\mathds{1}(\widehat{\beta}_{k}=\widehat{\beta}_{(j)}).$
(5)

The following corollary suggests that the proposed confidence interval for
$\beta_{(j)}$ can also serve as an exact prediction interval for
$\beta_{\widehat{j}}$. Therefore, the proposed procedure in Section 2.2 can
also be used to make inference on the observed top policies in a random
sample:

###### Corollary 1

Under Assumptions 1-5, we have that
$\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}\mathbb{P}\big{(}\widetilde{\beta}_{(j)}^{*}\leq\beta_{\widehat{j}}|\\{\bm{z}_{i,n}\\}_{i=1}^{n}\big{)}\leq
s\Big{)}=s$. Furthermore, $\widetilde{\beta}_{(j)}$ is a “root-$n$ consistent”
estimator of the data-dependent parameter $\beta_{\widehat{j}}$ in the sense
that $\forall\varepsilon>0$, there exists $M>0$ such that
$\mathbb{P}\big{(}|\sqrt{n}(\widetilde{\beta}_{(j)}-{\beta}_{\widehat{j}})|>M\big{)}\leq\varepsilon$,
for $n\geq 1$.

###### Remark 3 (Regression models with fixed effects)

The proposed resampling-based approach can be used to calibrate multiple best
policies when fixed effects are introduced in linear regression models (see
[49] for comprehensive discussion). This suggests that our approach not only
applies to independently sampled data, but also remains valid when there are
repeated-measurements present in the data. These may include short panel data,
and datasets in which, for example, two individuals have sampled from each
household. To conserve space in the main manuscript, we have leave the
detailed discussion in the Supplementary Materials (Section D).

###### Remark 4 (Data dependent choice of $d_{0}$)

In addition to a deterministic choice of $d_{0}$, another practically relevant
scenario is a data dependent choice of $d_{0}$. An example of such a data
dependent choice is $\widehat{d}_{0}=\max\\{k:\widehat{\beta}_{(k)}>C\\}$,
where $C$ is a user-specified threshold for the effect size. A relatively
complicated situation is that $C$ coincides with some of the policy sizes in
$\beta_{1},\beta_{2},\cdots,\beta_{d}$. In this situation, it is possible that
no matter how large $n$ is, $\widehat{d}_{0}$ does not converge to a
deterministic value but instead to a non-degenerate random variable. For the
purpose of separation, we may adjust
$\widehat{d}_{0}=\max\\{k:\widehat{\beta}_{(k)}>C\\}$ to be
$\widehat{d}^{\prime}_{0}=\max\\{k:\widehat{\beta}_{(k)}>C+C_{1}\cdot
n^{-0.25}\\}$, where $C_{1}$ is a constant that does not depend on $n$. The
choice of $-0.25$ is tunable and may be of independent interest. By this new
choice of $\widehat{d}^{\prime}_{0}$, the policy effects that exactly equal
$C$ will be eliminated almost surely when $n$ tends to infinity. This
elimination exactly matches the target to select all the policy sizes that are
larger than $C.$ In the limit of $n$ tending to infinity,
$\max\\{k:\widehat{\beta}_{(k)}>C+C_{1}\cdot n^{-0.25}\\}$ will converge
almost surely to a set that contains all effect sizes larger than $C$.
Therefore, the large-sample theory results for a pre-specified deterministic
integer would still hold by plugging in $\widehat{d}^{\prime}_{0}$.

###  4 Simulation studies

####  4.1 Simulation design

We generate i.i.d. Monte Carlo samples of
$\\{(y_{i,n},\bm{x}_{i,n}^{\prime},\bm{w}_{i,n}^{\prime})\\}_{i=1}^{n}$ from
the model

$\displaystyle
y_{i,n}=\bm{x}_{i,n}^{\prime}\bm{\beta}+\bm{w}_{i,n}^{\prime}\gamma_{n}+\varepsilon_{i,n},\quad
i=1,\ldots,n.$

We consider various data generating processes (DGP) for different choices of
the policy variable $\bm{x}_{i,n}$, the covariates $\bm{w}_{i,n}$ and the
random noise $\varepsilon_{i,n}$.  The first DGP follows a similar setup
taking from [37] and [12], where we generate many (sparse) dummy variables
entering the estimation of $\bm{\beta}$. We generate
$\bm{x}_{i,n}\sim\mathcal{N}(0,\Sigma)$ with $\Sigma_{jk}=0.5^{|j-k|}$ for
$j,k=1,\ldots,d$,
$\bm{w}_{i,n}=\mathds{1}(\widetilde{\bm{w}}_{i,n}\geq\Phi^{-1}(0.98))$ with
$\widetilde{\bm{w}}_{i,n}\sim\mathcal{N}(0,\bm{I}_{q_{n}})$ and
$\bm{I}_{q_{n}}$ is a $q_{n}$-dimensional identity matrix, and
$\varepsilon_{i,n}\sim\mathcal{N}(0,1)$. The second DGP considers a case with
dummy policy random variables and normal covariates, where we generate
$\bm{x}_{i,n}=\mathds{1}(\widetilde{\bm{x}}_{i,n}>0)$ with
$\widetilde{\bm{x}}_{i,n}\sim\mathcal{N}(0,\Sigma)$,
$\bm{w}_{i,n}\sim\mathcal{N}(0,\bm{I}_{q_{n}})$ and
$\varepsilon_{i,n}\sim\mathcal{N}(0,1)$. In the Supplementary Materials, we
have further included DGPs with more realistic error terms beyond normal
distribution, including error terms with asymmetric and bimodal distributions.
For most of the DGPs, we investigate both homoscedastic as well as
heteroscedastic models. See Supplementary Materials Section C for detailed
description and simulation results.

As for the coefficients, we consider three DGPs that vary in $\bm{\beta}$ and
$\bm{\gamma}_{n}$. The first DGP considers the case in which no policy is
effective (meaning that $\bm{\beta}=0$), and the coefficient $\gamma_{j}=1/j$,
for $j=1,\ldots,q_{n}$. We refer to this case as the “homogeneity” case since
$\beta_{j}$’s take the same value zero. The second and the third DGPs consider
cases where policy effects are generated from
$\beta_{j}=\Phi^{-1}\big{(}\frac{j}{d+1}\big{)}$ for $j=1,\ldots,d$, and the
coefficients are either $\bm{\gamma}_{n}=0$ or $\gamma_{j}=1/j$, for
$j=1,\ldots,q_{n}$. We refer to this case as the “heterogeneity(1)” case and
“heterogeneity(2)” case, respectively, since different policies have
heterogeneous effects.

We set the sample size $n\in\\{700,2000\\}$ to mimic the sample size in our
case studies, the number of policies $d\in\\{5,10\\}$, and the dimension of
the covariates $q_{n}$ from $q_{n}\in\\{1,141,281,421,561,631\\}$. All
statistics reported below are computed based on over 1,000 Monte Carlo
replications. To avoid redundancy, we present the results for $n=700$ and
$d=5$ in the main manuscript, and rests are provided in the Supplementary
Materials (Section C).

To demonstrate the robustness of the adopted covariance matrix estimator, we
compare our proposal with three alternative covariance matrix estimators. The
first one we compare with is the covariance matrix estimator proposed by [12]:
$\widehat{\bm{\Omega}}^{\texttt{HCK}}_{n}=\widehat{\bm{\Gamma}}_{n}^{-1}\widehat{\Sigma}^{\texttt{HCK}}_{n}\widehat{\bm{\Gamma}}_{n}^{-1}$,
where
$\widehat{\Sigma}^{\texttt{HCK}}_{n}\triangleq\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}\kappa_{ij,n}^{\texttt{HCK}}\widehat{\bm{v}}_{i,n}\widehat{\bm{v}}_{i,n}^{\prime}\widehat{u}_{j,n}^{2}$,
$\widehat{u}_{j,n}=\sum_{k=1}^{n}(\bm{M}_{n})_{j,k}(y_{k,n}-\bm{x}_{k,n}^{\prime}\widehat{\bm{\beta}})$,
and

$\displaystyle\kappa_{n}^{\texttt{HCK}}=\begin{pmatrix}M_{11,n}^{2}&\cdots&M_{1n,n}^{2}\\\
\vdots&\ddots&\vdots\\\ M_{n1,n}^{2}&\cdots&M_{nn,n}^{2}\\\
\end{pmatrix}^{-1}$ $\displaystyle=(\bm{M}_{n}\odot\bm{M}_{n})^{-1},$

with $\odot$ denoting the Hadamard product. The estimator
$\widehat{\Sigma}^{\texttt{HCK}}_{n}$ is well-defined whenever
$(\bm{M}_{n}\odot\bm{M}_{n})$ is invertible. We use the acronym “HCK” to
denote this estimator in the following parts. The second one we compare with
is the classical Eicker-White covariance matrix estimator [21, 51] of the
form:

$\displaystyle\widehat{\bm{\Omega}}^{\texttt{EW}}_{n}=\widehat{\bm{\Gamma}}_{n}^{-1}\widehat{\Sigma}^{\texttt{EW}}_{n}\widehat{\bm{\Gamma}}_{n}^{-1},$
(6)

where
$\widehat{\Sigma}^{\texttt{EW}}_{n}\triangleq\frac{1}{n}\sum_{i=1}^{n}\widehat{\bm{v}}_{i,n}\widehat{\bm{v}}_{i,n}^{\prime}\widehat{u}_{i,n}^{2}$
and
$\widehat{u}_{i,n}=\sum_{j=1}^{n}(\bm{M}_{n})_{i,j}(y_{j,n}-\bm{x}_{j,n}^{\prime}\widehat{\bm{\beta}})$.
We use the acronym “EW” to denote this estimator in our simulation results
section. Huber-Eicker-White standard error is also known as the HC0 standard
error, where HC stands for “heteroskedasticity robust.” The last covariance
matrix estimator we adopted is a variant of the HC0 estimator:

$\displaystyle\widehat{\bm{\Omega}}^{\texttt{HC3}}_{n}=\widehat{\bm{\Gamma}}_{n}^{-1}\widehat{\Sigma}^{\texttt{HC3}}_{n}\widehat{\bm{\Gamma}}_{n}^{-1},\quad\text{
where
}\widehat{\Sigma}^{\texttt{HC3}}_{n}\triangleq\frac{1}{n}\sum_{i=1}^{n}\widehat{\bm{v}}_{i,n}\widehat{\bm{v}}_{i,n}^{\prime}\frac{\widehat{u}_{i,n}^{2}}{(\bm{M}_{n})_{i,j}^{2}}.$
(7)

The above estimator upward reweights regression residuals, and we use the
acronym “HC3” to denote this estimator in our simulation results section.

####  4.2 Simulation results

We summarize our main takeaways from the simulation results presented in Table
1-3, where we have compared our proposed approach (“Proposed + KJ”) in Section
2.2 with four other methods. “Proposed + EW”, “Proposed + HC3”, and “Proposed
+ HCK” refer to methods adjusting for the winner’s curse bias but use
$\widehat{\bm{\Omega}}_{n}^{\texttt{EW}}$,
$\widehat{\bm{\Omega}}_{n}^{\texttt{HC3}}$, and
$\widehat{\bm{\Omega}}_{n}^{\texttt{HCK}}$, respectively, to estimate the
covariance matrix of $\bm{\beta}$. “No adjustment+KJ” refers to the approach
with no adjustment for the winner’s curse bias and adopts the robust
covariance matrix estimator proposed by [37] to make inference on the best
policies. We present the coverage probabilities and $\sqrt{n}$-scaled biases
for the top two policies in the population, i.e., $\beta_{(1)}$ and
$\beta_{(2)}$. As the simulation results are rather similar for the observed
top two policies in the random sample, i.e., $\beta_{\widehat{1}}$,
$\beta_{\widehat{2}}$, we present these results in the Supplementary Materials
(Section C.4).

Our simulation results confirm our theoretical results presented in Theorem 1.
When no policy is effective, our proposed method not only successfully
suppresses the winner’s curse bias for the top two policies but also attains
near nominal coverage (Table 2). Similar pattern can also be observed when top
policies are effective (i.e., $\beta_{j}$’s are heterogeneous, and Table 1 in
particular). In nearly all designs and for a range of considered values of
$q_{n}$, our proposal yields close to nominal confidence interval, though some
under coverage is observed for large values of $q_{n}$. The method with no
adjustment is obviously biased upward due to the winner’s curse phenomenon,
thus it provides under-covered confidence intervals and point estimates with
rather large biases. In all considered cases, both the EW-based method and the
HC3-based method tend to lose coverage when $q_{n}\geq 141$, and the HCK-based
method tends to produce under-covered confidence interval whenever $q_{n}\geq
561$. In moderately high dimensions so that $q_{n}/n$ is approximately one
half, the proposed method with the HCK variance estimator has comparable
performances with our approach.

Table 1:  Simulation results ($d=5,\text{heterogeneity},\beta_{(1)}$)

$\beta_{j}=\Phi^{-1}\big{(}\frac{j}{d+1}\big{)},\quad\bm{\gamma}_{n}=0$, j =
1,…, d
$\bm{x}_{i,n}\sim\mathcal{N}(0,\Sigma),\quad\bm{w}_{i,n}=\mathds{1}(\widetilde{\bm{w}}_{i,n}\geq\Phi^{-1}(0.98))$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.97(0.01) 0.96(0.01) 0.96(0.01) 0.96(0.01) 0.97(0.01) $\sqrt{n}$Bias
-0.04(0.06) -0.03(0.04) -0.03(0.04) -0.04(0.05) 0.05(0.06) $q_{n}=141$ Cover
0.96(0.01) 0.96(0.01) 0.94(0.01) 0.95(0.01) 0.95(0.01) $\sqrt{n}$Bias
-0.04(0.05) -0.04(0.04) -0.04(0.04) 0.06(0.06) 0.06(0.07) $q_{n}=281$ Cover
0.96(0.01) 0.95(0.01) 0.82(0.02) 0.80(0.01) 0.94(0.01) $\sqrt{n}$Bias
-0.05(0.06) -0.06(0.05) -0.06(0.03) -0.10(0.07) -0.08(0.08) $q_{n}=421$ Cover
0.95(0.02) 0.94(0.01) 0.79(0.01) 0.76(0.01) 0.78(0.01) $\sqrt{n}$Bias
-0.05(0.05) -0.06(0.05) -0.07(0.05) -0.12(0.09) 0.11(0.09) $q_{n}=561$ Cover
0.95(0.01) 0.92(0.01) 0.65(0.02) 0.63(0.01) 0.68(0.01) $\sqrt{n}$Bias
-0.07(0.07) -0.09(0.07) -0.17(0.10) -0.20(0.12) 0.15(0.13) $q_{n}=631^{*}$
Cover 0.93(0.01) 0.91(0.01) 0.51(0.02) 0.48(0.01) 0.55(0.01) $\sqrt{n}$Bias
-0.17(0.08) -0.19(0.10) -0.28(0.11) -0.35(0.22) -0.26(0.13)
$\bm{x}_{i,n}=\mathds{1}(\widetilde{\bm{x}}_{i,n}>0),\quad\bm{w}_{i,n}\sim\mathcal{N}(0,I)$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed + EW No adjustment+KJ $q_{n}=1$
Cover 0.97(0.01) 0.97(0.01) 0.95(0.01) 0.96(0.01) 0.97(0.01) $\sqrt{n}$Bias
-0.02(0.07) -0.02(0.04) -0.01(0.03) -0.07(0.11) -0.05(0.09) $q_{n}=141$ Cover
0.96(0.01) 0.95(0.01) 0.94(0.01) 0.94(0.01) 0.96(0.01) $\sqrt{n}$Bias
-0.02(0.03) -0.02(0.02) -0.03(0.02) 0.11(0.12) -0.06(0.12) $q_{n}=281$ Cover
0.95(0.01) 0.94(0.01) 0.87(0.01) 0.85(0.01) 0.95(0.01) $\sqrt{n}$Bias
-0.03(0.04) -0.03(0.03) -0.04(0.02) 0.14(0.12) -0.08(0.13) $q_{n}=421$ Cover
0.95(0.01) 0.94(0.01) 0.78(0.01) 0.76(0.01) 0.75(0.01) $\sqrt{n}$Bias
-0.03(0.03) -0.05(0.04) -0.08(0.05) -0.19(0.17) 0.19(0.14) $q_{n}=561$ Cover
0.95(0.01) 0.92(0.01) 0.63(0.02) 0.61(0.01) 0.63(0.01) $\sqrt{n}$Bias
-0.04(0.04) -0.08(0.06) -0.19(0.10) -0.30(0.26) -0.24(0.22) $q_{n}=631$ Cover
0.94(0.01) 0.91(0.01) 0.49(0.02) 0.45(0.01) 0.68(0.01) $\sqrt{n}$Bias
-0.06(0.07) -0.11(0.09) -0.23(0.13) -0.42(0.20) 0.39(0.29)

  * •

Note: “Cover” is the empirical coverage of the 95% confidence interval for
$\beta_{(1)}$ and “ $\sqrt{n}$Bias ” captures the root-$n$ scaled Monte Carlo
bias for estimating $\beta_{(1)}$. “ * ” indicates that
$\widehat{\bm{\Omega}}_{n}^{\texttt{KJ}}$ is not positive semi-definite in
some Monte Carlo samples.

Table 2:  Simulation results ($d=5,\text{homogeneity},\beta_{(1)}$)

$\bm{\beta}=0,\quad\bm{\gamma}_{j}=1/j$
$\bm{x}_{i,n}\sim\mathcal{N}(0,\Sigma),\quad\bm{w}_{i,n}=\mathds{1}(\widetilde{\bm{w}}_{i,n}\geq\Phi^{-1}(0.98))$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.97(0.01) 0.96(0.01) 0.96(0.01) 0.93(0.02) 0.90(0.01) $\sqrt{n}$Bias
0.02(0.03) 0.02(0.03) 0.03(0.03) 0.03(0.03) 1.64(0.04) $q_{n}=141$ Cover
0.96(0.01) 0.96(0.01) 0.96(0.01) 0.89(0.02) 0.88(0.01) $\sqrt{n}$Bias
0.03(0.04) 0.04(0.04) 0.03(0.04) 0.12(0.04) 1.78(0.04) $q_{n}=281$ Cover
0.96(0.01) 0.94(0.01) 0.90(0.02) 0.85(0.01) 0.83(0.01) $\sqrt{n}$Bias
0.03(0.04) 0.04(0.04) 0.05(0.04) 0.22(0.03) 2.03(0.05) $q_{n}=421$ Cover
0.95(0.01) 0.93(0.01) 0.82(0.02) 0.79(0.01) 0.74(0.02) $\sqrt{n}$Bias
0.05(0.05) 0.18(0.05) 0.24(0.06) 0.36(0.03) 2.63(0.06) $q_{n}=561$ Cover
0.95(0.01) 0.93(0.01) 0.67(0.02) 0.73(0.01) 0.63(0.02) $\sqrt{n}$Bias
0.08(0.09) 0.51(0.05) 0.74(0.06) 0.44(0.04) 3.74(0.09) $q_{n}=631^{*}$ Cover
0.93(0.01) 0.89(0.01) 0.53(0.02) 0.50(0.01) 0.45(0.02) $\sqrt{n}$Bias
0.18(0.09) 1.21(0.09) 1.84(0.11) 2.42(0.06) 5.10(0.12)
$\bm{x}_{i,n}=\mathds{1}(\widetilde{\bm{x}}_{i,n}>0),\quad\bm{w}_{i,n}\sim\mathcal{N}(0,I)$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.96(0.01) 0.96(0.01) 0.96(0.01) 0.94(0.01) 0.90(0.01) $\sqrt{n}$Bias
0.03(0.04) 0.04(0.05) 0.05(0.05) 0.07(0.07) 2.75(0.06) $q_{n}=141$ Cover
0.96(0.01) 0.96(0.01) 0.93(0.01) 0.90(0.01) 0.83(0.01) $\sqrt{n}$Bias
0.05(0.05) 0.05(0.06) 0.17(0.07) 0.31(0.07) 3.29(0.08) $q_{n}=281$ Cover
0.95(0.01) 0.95(0.01) 0.90(0.01) 0.88(0.01) 0.75(0.02) $\sqrt{n}$Bias
0.07(0.08) 0.07(0.07) 0.31(0.08) 0.54(0.05) 3.59(0.08) $q_{n}=421$ Cover
0.95(0.01) 0.95(0.01) 0.85(0.01) 0.86(0.01) 0.65(0.02) $\sqrt{n}$Bias
0.04(0.04) 0.10(0.11) 0.73(0.13) 0.64(0.06) 4.58(0.11) $q_{n}=561$ Cover
0.93(0.01) 0.90(0.02) 0.59(0.02) 0.61(0.01) 0.53(0.02) $\sqrt{n}$Bias
0.13(0.07) 0.19(0.13) 2.00(0.13) 1.73(0.08) 5.90(0.13) $q_{n}=631$ Cover
0.90(0.01) 0.78(0.02) 0.38(0.02) 0.30(0.01) 0.33(0.02) $\sqrt{n}$Bias
0.50(0.12) 2.47(0.16) 5.16(0.19) 7.68(0.18) 6.51(0.19)

  * •

Note: “Cover” is the empirical coverage of the 95% confidence interval for
$\beta_{(1)}$ and “ $\sqrt{n}$Bias ” captures the root-$n$ scaled Monte Carlo
bias for estimating $\beta_{(1)}$. “ * ” indicates that
$\widehat{\bm{\Omega}}_{n}^{\texttt{KJ}}$ is not positive semi-definite in
some Monte Carlo samples.

Table 3:  Simulation results ($d=5,\text{heterogeneity},\beta_{(2)}$)

$\beta_{j}=\Phi^{-1}\big{(}\frac{j}{d+1}\big{)},\quad\bm{\gamma}_{n}=0$
$\bm{x}_{i,n}\sim\mathcal{N}(0,\Sigma),\quad\bm{w}_{i,n}=\mathds{1}(\widetilde{\bm{w}}_{i,n}\geq\Phi^{-1}(0.98))$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.97(0.01) 0.96(0.01) 0.96(0.01) 0.96(0.01) 0.96(0.01) $\sqrt{n}$Bias
0.02(0.06) 0.02(0.06) 0.01(0.06) -0.05(0.07) -0.04(0.07) $q_{n}=141$ Cover
0.97(0.01) 0.96(0.01) 0.93(0.01) 0.94(0.01) 0.94(0.01) $\sqrt{n}$Bias
-0.02(0.04) -0.02(0.04) -0.04(0.03) -0.06(0.06) -0.07(0.08) $q_{n}=281$ Cover
0.95(0.01) 0.94(0.01) 0.88(0.02) 0.89(0.01) 0.85(0.01) $\sqrt{n}$Bias
-0.03(0.04) -0.03(0.03) -0.07(0.03) -0.11(0.10) 0.19(0.11) $q_{n}=421$ Cover
0.95(0.01) 0.94(0.02) 0.81(0.02) 0.80(0.01) 0.77(0.01) $\sqrt{n}$Bias
-0.03(0.03) -0.03(0.04) -0.10(0.06) -0.15(0.12) -0.23(0.15) $q_{n}=561$ Cover
0.95(0.01) 0.94(0.02) 0.67(0.02) 0.65(0.01) 0.63(0.01) $\sqrt{n}$Bias
0.03(0.03) 0.05(0.06) 0.12(0.08) -0.17(0.13) -0.26(0.18) $q_{n}=631^{*}$ Cover
0.94(0.01) 0.93(0.01) 0.56(0.02) 0.53(0.01) 0.50(0.01) $\sqrt{n}$Bias
-0.07(0.07) -0.18(0.06) -0.23(0.08) -0.26(0.17) 0.47(0.22)
$\bm{x}_{i,n}=\mathds{1}(\widetilde{\bm{x}}_{i,n}>0),\quad\bm{w}_{i,n}\sim\mathcal{N}(0,I)$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.98(0.01) 0.98(0.01) 0.98(0.01) 0.96(0.01) 0.98(0.01) $\sqrt{n}$Bias
-0.08(0.12) -0.09(0.12) -0.10(0.13) 0.09(0.12) -0.07(0.10) $q_{n}=141$ Cover
0.97(0.01) 0.97(0.01) 0.95(0.01) 0.95(0.01) 0.97(0.01) $\sqrt{n}$Bias
-0.09(0.12) -0.10(0.13) -0.12(0.13) 0.09(0.13) -0.08(0.10) $q_{n}=281$ Cover
0.97(0.01) 0.97(0.01) 0.90(0.01) 0.87(0.01) 0.96(0.01) $\sqrt{n}$Bias
-0.11(0.14) -0.10(0.14) -0.15(0.14) -0.18(0.14) -0.10(0.11) $q_{n}=421$ Cover
0.96(0.01) 0.95(0.02) 0.80(0.02) 0.75(0.01) 0.94(0.01) $\sqrt{n}$Bias
0.14(0.16) -0.16(0.18) -0.20(0.18) -0.22(0.17) -0.15(0.15) $q_{n}=561$ Cover
0.96(0.01) 0.94(0.02) 0.63(0.02) 0.60(0.01) 0.92(0.01) $\sqrt{n}$Bias
0.14(0.18) 0.20(0.23) -0.24(0.22) -0.30(0.23) 0.19(0.20) $q_{n}=631$ Cover
0.94(0.01) 0.93(0.02) 0.58(0.02) 0.55(0.01) 0.52(0.01) $\sqrt{n}$Bias
0.15(0.20) 0.24(0.26) 0.28(0.25) -0.35(0.13) 0.56(0.24)

  * •

Note: “Cover” is the empirical coverage of the 95% confidence interval for
$\beta_{(2)}$ and “ $\sqrt{n}$Bias ” captures the root-$n$ scaled Monte Carlo
bias for estimating $\beta_{(2)}$. “ * ” indicates that
$\widehat{\bm{\Omega}}_{n}^{\texttt{KJ}}$ is not positive semi-definite in
some Monte Carlo samples.

###  5 Case studies

####  5.1 Case study I: Charitable giving

In the past half century, charitable giving by individuals in the United
States has grown and it has contributed to more than two percent of the annual
GDP since 1998 [42]. Charitable giving is often driven by altruism, while as
suggested by many field experiments, improper policies adopted by the demand
side–fundraisers–may impair the supply side’s (individual donors) motivation
of giving [4]. Therefore, to effectively attract resources from individual
donors, fundraisers need to properly design donation incentives. One of the
donation incentives is matching grant which means that a matching donor
pledges to match any donation from other donors with certain ratio and up to
some threshold. As the price elasticity of matching donation may differ from
other donation incentives, we hope to carefully investigate different pricing
policies in a matching donation and study if the observed top two performing
policies are indeed effective.

We work with the charitable giving data in [38]. [38] conduct a field
experiment that explores the price elasticity in a matching donation. The
field experiment involves $50,083$ previous donors to a political charity.
Individuals are randomly assigned to two groups: treatment ($n=33,396$) and
control ($n=16,687$). In the control group, individuals receive a standard
letter with no matching details. In the treatment group, each potential donor
receives a letter with three strategies: (1) match ratio, (2) match size, and
(3) ask amount. Within each strategy, individuals are randomly assigned to a
sub-policy detailed below.

For the match ratio strategy, there are three sub-policies: (1) 1:1 (the
matching donor contributes the same amount as the individual donor), (2) 2:1
(the matching donor contributes twice as many as the individual donor), (3)
3:1 (the matching donor contributes three times as many as the individual
donor). For the match size strategy, there are four sub-policies with
different pledge amounts: (1) $25,000, (2) $50,000, (3) $100,000, and (4)
unstated amount. For the ask amount strategy, individual donors are asked to
give same amount, 25% more or 50% more than their largest past donation.

In our study, we focus on the treatment “ask amount” with three pricing
policies, and we study the subpopulation ($n=7,938$) of unmarried males living
in red counties or red states. Red county (state) refers to a county (state)
in which residents predominantly vote for the Republican Party. The outcome of
interest is the donation amount. We have adjusted $q_{n}=1,049$ covariates
including the donors’ demographic information (26 variables), census
information (27 variables), and their two-way interaction terms. Our results
are summarized in Table 4.

Method | Policies(Ask amount) | Est (95% CI) |  $p$-value |
---|---|---|---|---
Uncalibrated | Same | 0.67 (0.09, 1.25) | 0.023* |
| 25% more | 0.66 (0.01, 1.31) | 0.046* |
| 50% more | 0.33 (-0.21, 0.86) | 0.235 |
Calibrated | Same | 0.63 (0.10, 1.20) | 0.025* |
| 25% more | 0.56(-0.01, 1.07) | 0.052 |
Table 4: Estimated treatment effects (Est), 95% confidence intervals (95% CI),
and two-sided $p$-values for the three “ask amount” policies. “Uncalibrated”
refers to the study results obtained without any adjustment, and the
confidence intervals are constructed based on normal approximation with the
estimated covariance matrix $\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}$.
“Calibrated” refers to our proposed methodology. The computational time is
$741$ seconds on a Lenovo NeXtScale nx360m5 node (24 cores per node) equipped
with Intel Xeon Haswell processor.

Results in Table 4 suggest that, without any calibration, asking the donor
either to give the same amount or to give 25% more than their highest past
donation seems to be the best policies that significantly increase the
donation amount. Specifically, our results from running a simple linear
regression model suggest that asking the individual donor to give the same
amount of their largest past donation appears to be the most effective pricing
policy, and it on average raises $0.67 (95% CI = (0.09, 1.25), $p$-value =
0.023) per donor. Asking the individual donor to give 25% more than their
largest past donation is the second most effective policy, with an increased
donation by $0.66 (95% CI = (0.01, 1.31), $p$-value = 0.046) per donor.

Because we pick the most effective policies from a random sample, these
estimates are potentially subject to the winner’s curse bias. We thus apply
the proposed method to carefully examine these seemly effective policies.
After calibrating for the winner’s curse bias, we confirm that the asking for
the same amount policy remains as the most effective policy, though with a
slightly smaller estimated effect size (Est = $0.63, 95% CI = (0.10, 1.20),
$p$-value = 0.025). This result is moderately aligned with the analysis in
[38], whose results suggest that donors from red states or red counties are
more willing to contribute, partially because the collaborating charity is
politically oriented. However, for the effect of the second best policy–asking
to donate 25% more than past donation–is shifted downward, and it no longer
has significant impact in promoting the donation amount (Est = $0.56, 95% CI =
(-0.01, 1.07), $p$-value = 0.052). This result might be partially explained by
the observation that donors are more motivated by a lower “price” of donation
[50]. In sum, our analyses suggest that the best pricing policy of charitable
giving for unmarried males living in the Republican Party dominated voting
regions could be asking for the same amount as their highest previous
donation, and asking for more donations may not incentivize the donors to
give. Though given the obtained p-value before and after calibration for the
second best policy is rather close to the 5 percent threshold, we note that
such a conclusion should also be viewed with caution.

####  5.2 Case study II: National supported work (NSW) program

In this case study, we revisit a dataset from the National Supported Work
(NSW) program. The NSW program is a labor training program implemented in
1970’s that provides work experience to disadvantaged workers. Our proposed
method can also be used to evaluate if the job training program is indeed
beneficial for certain groups of workers. To do so, the structural component
$\bm{x}_{i,n}$ in the model (1) would include variables representing the
interactions between the treatment variable (the job training program) and
different subgroup indicator variables of interest.

We use the field experiment dataset adopted in [19] ($n=455$), in which 185
workers are in the treatment group and 260 workers are in the control group.
This dataset consists of a treatment indicator variable, an outcome variable
measured by the participant post-treatment earnings in 1978, and eight
baseline variables ( including age, years of education, an indicator for high
school degree, indicators for Black and Hispanic, marital status, and pre-
treatment earnings in 1974 and 1975). We further add three sets of additional
covariates following the setup in [31]: (1) $\mathds{1}$(1974 earnings =0) and
$\mathds{1}$(1975 earnings =0); (2) all first-order interactions; (3) all
polynomials up to the 2nd-order. The final dataset includes $51$ covariates.
We aim to investigate the effectiveness of the NSW program in four groups of
workers: (1) married Black workers, (2) unmarried Black workers, (3) married
Non-Black workers, and (4) unmarried Non-Black workers. The summarized results
are shown in Table 5.

Method | Subgroups | Est (95% CI) ($$10^{3}$) |  $p$-value |
---|---|---|---|---
Uncalibrated | Black, married | 4.35 (0.89, 7.81) | 0.014* |
| Black, unmarried | 1.10(-0.55, 2.75) | 0.190 |
| Non-Black, married | 1.33(-6.63, 9.29) | 0.743 |
| Non-Black, unmarried | 1.40(-2.61, 5.40) | 0.494 |
Calibrated | Black, married | 4.41(1.74, 8.50) | 0.009* |
Table 5: Estimated treatment effects (Est), 95% confidence intervals (CI), in
units $$10^{3}$/year, and two-sided $p$-values for the four subgroups in the
NSW study ($n=445$, $q_{n}=51$). “Uncalibrated” refers to the study results
obtained without any adjustment, and the confidence intervals are constructed
based on normal approximation with the estimated covariance matrix
$\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}$. “Calibrated” refers to our proposed
methodology. The computational time is $122$ seconds on a Lenovo NeXtScale
nx360m5 node (24 cores per node) equipped with Intel Xeon Haswell processor.

Table 5 demonstrates that without adjusting for the winner’s curse bias,
married Black workers (estimated treatment effect = 4.35, 95% CI = (0.89,
7.81), $p$-value = 0.014, in units $\$10^{3}$) seem to benefit from the
program the most. After accounting for the winner’s curse bias issue, our
approach potentially confirms that the treatment effect of the NSW program for
the married Black workers is still significant, and the calibrated treatment
effect remains roughly the same ( Est= 4.41, 95% CI = (1.74, 8.50), $p$-value
= 0.009, in units $\$10^{3}$/year).

The dataset collected from the NSW program has been frequently analyzed in the
past decade, and our results are largely in-line with current understandings
gathered in past studies. For example, although not focusing on the same
groups of workers, [36] suggest that married and unemployed Black workers with
some college education have increased their post-treatment earnings for about
38%. [19] show that the job training program yields positive treatment effect
on the overall Black participants. In this case study, our approach may help
to confirm the seemly effective subgroup observed in a random sample while
providing a statistically justified estimate accounting for the winner’s curse
bias.

###  6 Concluding remarks

In this article, we have introduced an approach to evaluate multiple best
policies based on resampling in the context of a linear model with many
covariates. While our approach is numerically reliable and theoretically
grounded, it is worthwhile to generalize our framework so that the policy
effects can be estimated with other off-shelf methods that are, for example,
robust to the high-dimensional confounders or to the presence of interference
and noncompliance. Our current theoretical analysis suggests that our proposed
approach can be readily extended as long as the covariance matrix between
different policies can be consistently estimated. It is thus desirable for us
to provide a general framework to broaden future applications for other
disciplines in general.

Acknowledgement. The research of Jingshen Wang is supported in part by the
National Science Foundation (DMS 2015325) and the National Institute of Health
(R01MH125746).

### References

  * [1]
  * Abadie and Cattaneo [2018] Abadie, A. and Cattaneo, M. D. (2018). “Econometric methods for program evaluation,” Annual Review of Economics, 10, 465–503.
  * Anatolyev [2012] Anatolyev, S. (2012). “Inference in regression models with many regressors,” Journal of Econometrics, 170(2), 368–382.
  * Andreoni and Miller [2002] Andreoni, J. and Miller, J. (2002). “Giving according to GARP: An experimental test of the consistency of preferences for altruism,” Econometrica, 70(2), 737–753.
  * Andrews [2000] Andrews, D. W. (2000). “Inconsistency of the bootstrap when a parameter is on the boundary of the parameter space,” Econometrica, 399–405.
  * Andrews et al. [2022+] Andrews, I., Bowen, D., Kitagawa, T., and McCloskey, A. (2022+). “Inference for Losers,” American Economic Association Papers and Proceedings.
  * Andrews et al. [2019] Andrews, I., Kitagawa, T., and McCloskey, A. (2019). “Inference on winners,”Technical report, National Bureau of Economic Research.
  * Athey and Imbens [2017] Athey, S. and Imbens, G. W. (2017). “The state of applied econometrics: Causality and policy evaluation,” Journal of Economic Perspectives, 31(2), 3–32.
  * Belloni et al. [2014] Belloni, A., Chernozhukov, V., and Hansen, C. (2014). “Inference on treatment effects after selection among high-dimensional controls,” The Review of Economic Studies, 81(2), 608–650.
  * Cattaneo et al. [2019] Cattaneo, M. D., Jansson, M., and Ma, X. (2019). “Two-step estimation and inference with possibly many included covariates,” The Review of Economic Studies, 86(3), 1095–1122.
  * Cattaneo et al. [2018a] Cattaneo, M. D., Jansson, M., and Newey, W. K. (2018a). “Alternative asymptotics and the partially linear model with many regressors,” Econometric Theory, 34(2), 277–301.
  * Cattaneo et al. [2018b] Cattaneo, M. D., Jansson, M., and Newey, W. K. (2018b). “Inference in linear regression models with many covariates and heteroscedasticity,” Journal of the American Statistical Association, 113(523), 1350–1361.
  * Charness and Levin [2009] Charness, G. and Levin, D. (2009). “The origin of the winner’s curse: a laboratory study,” American Economic Journal: Microeconomics, 1(1), 207–36.
  * Chernozhukov et al. [2018] Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018). “Double/debiased machine learning for treatment and structural parameters.”
  * Chernozhukov et al. [2013a] Chernozhukov, V., Chetverikov, D., and Kato, K. (2013a). “Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors,” The Annals of Statistics, 41(6), 2786–2819.
  * Chernozhukov et al. [2013b] Chernozhukov, V., Lee, S., and Rosen, A. M. (2013b). “Intersection bounds: estimation and inference,” Econometrica, 81(2), 667–737.
  * Claggett et al. [2014] Claggett, B., Xie, M., and Tian, L. (2014). “Meta-analysis with fixed, unknown, study-specific parameters,” Journal of the American Statistical Association, 109(508), 1660–1671.
  * Dawid [1994] Dawid, A. (1994). “Selection paradoxes of Bayesian inference,” In Multivariate analysis and its applications: Institute of Mathematical Statistics, 211–220.
  * Dehejia and Wahba [2002] Dehejia, R. H. and Wahba, S. (2002). “Propensity score-matching methods for nonexperimental causal studies,” Review of Economics and statistics, 84(1), 151–161.
  * Efron [2011] Efron, B. (2011). “Tweedie’s formula and selection bias,” Journal of the American Statistical Association, 106(496), 1602–1614.
  * Eicker [1963] Eicker, F. (1963). “Asymptotic normality and consistency of the least squares estimators for families of linear regressions,” The annals of mathematical statistics, 447–456.
  * El Karoui et al. [2013] El Karoui, N., Bean, D., Bickel, P. J., Lim, C., and Yu, B. (2013). “On robust regression with high-dimensional predictors,” Proceedings of the National Academy of Sciences, 110(36), 14557–14562.
  * El Karoui and Purdom [2018] El Karoui, N. and Purdom, E. (2018). “Can we trust the bootstrap in high-dimensions? The case of linear models,” The Journal of Machine Learning Research, 19(1), 170–235.
  * Fan et al. [2007] Fan, J., Hall, P., and Yao, Q. (2007). “To how many simultaneous hypothesis tests can normal, student’s t or bootstrap calibration be applied?” Journal of the American Statistical Association, 102(480), 1282–1288.
  * Fan et al. [2014] Fan, J., Han, F., and Liu, H. (2014). “Challenges of big data analysis,” National science review, 1(2), 293–314.
  * Fan et al. [2021] Fan, J., Imai, K., Lee, I., Liu, H., Ning, Y., and Yang, X. (2021). “Optimal covariate balancing conditions in propensity score estimation,” Journal of Business & Economic Statistics, 1–14.
  * Fan et al. [2020] Fan, J., Li, R., Zhang, C.-H., and Zou, H. (2020). Statistical foundations of data science: Chapman and Hall/CRC.
  * Fan et al. [2011] Fan, J., Lv, J., and Qi, L. (2011). “Sparse high-dimensional models in economics,” Annu. Rev. Econ., 3(1), 291–317.
  * Fan et al. [2018] Fan, J., Shao, Q.-M., and Zhou, W.-X. (2018). “Are discoveries spurious? Distributions of maximum spurious correlations and their applications,” Annals of statistics, 46(3), 989.
  * Fan and Zhou [2016] Fan, J. and Zhou, W.-X. (2016). “Guarding against spurious discoveries in high dimensions,” The Journal of Machine Learning Research, 17(1), 7123–7156.
  * Farrell [2015] Farrell, M. H. (2015). “Robust inference on average treatment effects with possibly more covariates than observations,” Journal of Econometrics, 189(1), 1–23.
  * Geer et al. [2014] Van de Geer, S., Bühlmann, P., Ritov, Y., and Dezeure, R. (2014). “On asymptotically optimal confidence regions and tests for high-dimensional models,” The Annals of Statistics, 42(3), 1166–1202.
  * Guo and He [2020] Guo, X. and He, X. (2020). “Inference on selected subgroups in clinical trials,” Journal of the American Statistical Association, 1–19.
  * Hall and Miller [2010] Hall, P. and Miller, H. (2010). “Bootstrap confidence intervals and hypothesis tests for extrema of parameters,” Biometrika, 97(4), 881–892.
  * Huber [1973] Huber, P. J. (1973). “Robust regression: asymptotics, conjectures and Monte Carlo,” The annals of statistics, 799–821.
  * Imai and Ratkovic [2013] Imai, K. and Ratkovic, M. (2013). “Estimating treatment effect heterogeneity in randomized program evaluation,” The Annals of Applied Statistics, 7(1), 443–470.
  * Jochmans [2020] Jochmans, K. (2020). “Heteroscedasticity-robust inference in linear regression models with many covariates,” Journal of the American Statistical Association, 1–10.
  * Karlan and List [2007] Karlan, D. and List, J. A. (2007). “Does price matter in charitable giving? Evidence from a large-scale natural field experiment,” American Economic Review, 97(5), 1774–1793.
  * Kueck et al. [2017] Kueck, J., Luo, Y., Spindler, M., and Wang, Z. (2017). “Estimation and Inference of Treatment Effects with $L\\_2$-Boosting in High-Dimensional Settings,” arXiv preprint arXiv:1801.00364.
  * LaLonde [1986] LaLonde, R. J. (1986). “Evaluating the econometric evaluations of training programs with experimental data,” The American economic review, 604–620.
  * Lee and Shen [2018] Lee, M. R. and Shen, M. (2018). “Winner’s curse: Bias estimation for total effects of features in online controlled experiments,” In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 491–499.
  * List [2011] List, J. A. (2011). “The market for charitable giving,” Journal of Economic Perspectives, 25(2), 157–80.
  * Long and Ervin [2000] Long, J. S. and Ervin, L. H. (2000). “Using heteroscedasticity consistent standard errors in the linear regression model,” The American Statistician, 54(3), 217–224.
  * Mammen [1993] Mammen, E. (1993). “Bootstrap and wild bootstrap for high dimensional linear models,” The annals of statistics, 255–285.
  * Mammen [2012] Mammen, E. (2012). When does bootstrap work?: asymptotic results and simulations, 77: Springer Science & Business Media.
  * Rai [2018] Rai, Y. (2018). “Statistical inference for treatment assignment policies,” Unpublished Manuscript.
  * Stock and Watson [2008] Stock, J. H. and Watson, M. W. (2008). “Heteroskedasticity-robust standard errors for fixed effects panel data regression,” Econometrica, 76(1), 155–174.
  * Vazquez-Bare [2021] Vazquez-Bare, G. (2021). “Identification and estimation of spillover effects in randomized experiments,” Journal of Econometrics, Forthcoming.
  * Verdier [2020] Verdier, V. (2020). “Estimation and inference for linear models with two-way fixed effects and sparsely matched data,” Review of Economics and Statistics, 102(1), 1–16.
  * Warwick [2003] Warwick, M. (2003). Testing, Testing 1, 2, 3: Raise More Money with Direct Mail Tests: John Wiley & Sons.
  * White [1980] White, H. (1980). “A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity,” Econometrica: journal of the Econometric Society, 817–838.
  * Xie et al. [2009] Xie, M., Singh, K., and Zhang, C.-H. (2009). “Confidence intervals for population ranks in the presence of ties and near ties,” Journal of the American Statistical Association, 104(486), 775–788.
  * Zhang and Zhang [2014] Zhang, C.-H. and Zhang, S. S. (2014). “Confidence intervals for low dimensional parameters in high dimensional linear models,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1), 217–242.

## Appendix

###  Appendix A Theorem 1

####  A.1 Review of notations and assumptions

We denote the sample
$\\{(y_{i,n},\bm{x}_{i,n}^{\prime},\bm{w}_{i,n}^{\prime})^{\prime}\\}_{i=1}^{n}$
as $\\{\bm{z}_{i,n}\\}_{i=1}^{n}$. Recall $u_{i,n}$ is the random error in the
considered linear model:

$\displaystyle
y_{i,n}=\bm{x}_{i,n}^{\prime}\bm{\beta}+\bm{w}_{i,n}^{\prime}\bm{\gamma}_{n}+u_{i,n},\quad
i=1,\ldots,n,$ (A.1)

we define

$\displaystyle\varepsilon_{i,n}=u_{i,n}-\mathbb{E}[u_{i,n}|\\{\bm{w}_{i,n}\\}_{i=1}^{n},\\{\bm{x}_{i,n}\\}_{i=1}^{n}],\quad\bm{v}_{i,n}=\bm{x}_{i,n}-\mathbb{E}[\bm{x}_{i,n}|\\{\bm{w}_{i,n}\\}_{i=1}^{n}],\quad
i=1,\ldots,n.$ (A.2)

Let
$e_{i,n}=\mathbb{E}[u_{i,n}|\\{\bm{w}_{i,n}\\}_{i=1}^{n},\\{\bm{x}_{i,n}\\}_{i=1}^{n}]$,
we further denote

$\displaystyle\widehat{u}_{i}=\sum_{j=1}^{n}(\bm{M}_{n})_{i,j}(y_{j,n}-\bm{x}_{j,n}^{\prime}\widehat{\bm{\beta}}),\quad\widehat{\bm{v}}_{i,n}=\sum_{i=1}^{n}(\bm{M}_{n})_{i,j}\bm{x}_{j,n},$
$\displaystyle(\bm{M}_{n})_{i,j}=\mathds{1}(i=j)-\bm{w}^{\prime}_{i,n}\Big{(}\sum_{k=1}^{n}\bm{w}_{k,n}\bm{w}_{k,n}^{\prime}\Big{)}^{-1}\bm{w}_{j,n},$
$\displaystyle\sigma_{i,n}^{2}=\mathbb{E}[\varepsilon_{i,n}^{2}|\\{\bm{w}_{i,n}\\}_{i=1}^{n},\\{\bm{x}_{i,n}\\}_{i=1}^{n}],\quad\widetilde{\bm{v}}_{i,n}=\sum_{j=1}^{n}(\bm{M}_{n})_{i,j}\bm{v}_{j,n},$
$\displaystyle\rho_{n}^{1}=\frac{\sum_{i=1}^{n}\mathbb{E}\big{[}e_{i,n}^{2}\big{]}}{n},\quad\rho_{n}^{2}=\frac{\sum_{i=1}^{n}\mathbb{E}\Big{[}\mathbb{E}\big{(}e_{i,n}|\\{\bm{w}_{i,n}\\}_{i=1}^{n}\big{)}^{2}\Big{]}}{n},$
$\displaystyle\bm{Q}_{i,n}=\mathbb{E}\Big{[}\bm{x}_{i,n}-\big{(}\sum_{j=1}^{n}\mathbb{E}[\bm{w}_{j,n}\bm{w}_{j,n}^{\prime}]\big{)}^{-1}\sum_{j=1}^{n}\mathbb{E}[\bm{w}_{j,n}\bm{x}^{\prime}_{j,n}]\Big{|}\\{\bm{w}_{i,n}\\}_{i=1}^{n}\Big{]}.$

We will use the notation $\mathbb{P}(\cdot|\\{\bm{z}_{i,n}\\}_{i=1}^{n})$ to
refer to the probability that is conditional on the random variables
$\\{\bm{z}_{i,n}\\}_{i=1}^{n}$.

For a policy $j$, recall our definition of the near tie set:

$\displaystyle\mathcal{H}_{(j)}=\big{\\{}k:|\beta_{k}-\beta_{(j)}|=o(n^{-\frac{1}{2}}),\
k=1,\ldots,d\big{\\}}.$

This suggests that $\forall k\in\mathcal{H}_{(j)}$, there exist a sequence
$\delta_{n}\rightarrow 0$ as $n\rightarrow 0$, such that

$\displaystyle\beta_{k}=\beta_{(j)}+n^{-\frac{1}{2}}\cdot\delta_{n},\quad\forall
k\in\mathcal{H}_{(j)}.$

Next, let $\widehat{\bm{e}}_{j}$ denote a $d$-dimensional (sparse) vector
based on the estimated tie set $\widehat{\mathcal{H}}_{(j)}$ with

$\displaystyle\widehat{\bm{e}}_{j}=(\widehat{e}_{j,1},\ldots,\widehat{e}_{j,d}),\quad\widehat{e}_{jk}=\frac{\mathds{1}(k\in\widehat{\mathcal{H}}_{(j)})}{|\widehat{\mathcal{H}}_{(j)}|},\quad
k=1,\ldots,d,$

and define a $d-$dimensional sparse index vector based on the true near-tie
set ${\mathcal{H}}_{(j)}$ as

$\displaystyle{\bm{e}}_{j}=({e}_{j,1},\ldots,{e}_{j,d}),\quad{e}_{jk}=\frac{\mathds{1}(k\in{\mathcal{H}}_{(j)})}{|{\mathcal{H}}_{(j)}|},\quad
k=1,\ldots,d.$ (A.3)

We make following assumptions throughout this section:

###### Assumption 6 (Sampling)

The errors $\varepsilon_{i,n}$ are uncorrelated across $i$ conditional on
$\\{\bm{x}_{i,n}\\}_{i=1}^{n}$ and $\\{\bm{w}_{i,n}\\}_{i=1}^{n}$. Let
$\\{N_{1},\ldots,N_{G_{n}}\\}$ represents a partition of $\\{1,\ldots,n\\}$
with $\underset{1\leq g\leq G_{n}}{\max}{|N_{g}|}=O(1)$ such that
$\\{(\varepsilon_{i,n},\bm{v}_{i,n}),i\in N_{g}\\}$ (defined in (A.2)) are
independent across $g$ conditional on $\\{\bm{w}_{i,n}\\}_{i=1}^{n}$.

###### Assumption 7 (Design)

The dimension of the covariates $\bm{w}_{i,n}$ satisfies that
$\limsup_{n\rightarrow\infty}q_{n}/n<1$. The minimum eigenvalue of of the
matrix $\sum_{i=1}^{n}\bm{w}_{i,n}\bm{w}_{i,n}^{\prime}$ is bounded away from
0 with probability approaching one, that is

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}\lambda_{\min}\big{(}\sum_{i=1}^{n}\bm{w}_{i,n}\bm{w}_{i,n}^{\prime}\big{)}>0\big{)}\Big{)}=1.$

Lastly,

$\displaystyle\underset{1\leq i\leq
n}{\max}\Big{\\{}\mathbb{E}[\varepsilon_{i,n}^{4}|\\{\bm{w}_{i,n}\\}_{i=1}^{n},\\{\bm{x}_{i,n}\\}_{i=1}^{n}],\
\frac{1}{\sigma_{i,n}^{2}},$
$\displaystyle\quad\quad\quad\quad\mathbb{E}[\bm{v}_{i,n}^{4}|\\{\bm{w}_{i,n}\\}_{i=1}^{n}],\
1/\lambda_{\min}\big{(}\frac{\sum_{i=1}^{n}\mathbb{E}\big{[}\widetilde{\bm{v}}_{i,n}\widetilde{\bm{v}}_{i,n}^{\prime}|\\{\bm{w}_{i,n}\\}_{i=1}^{n}\big{]}}{n}\big{)}\Big{\\}}=O_{p}(1).$

###### Assumption 8 (Linear model approximation)

$\sum_{i=1}^{n}\mathbb{E}[||Q_{i,n}||^{2}]/n=O(1)$,
$\rho_{n}^{1}+n(\rho_{n}^{1}-\rho_{n}^{2})+\rho_{n}^{2}\cdot\sum_{i=1}^{n}\mathbb{E}[||Q_{i,n}||^{2}]=o(1)$,
and $\underset{1\leq i\leq
n}{\max}||\widehat{\bm{v}}_{i,n}||/\sqrt{n}=o_{p}(1)$, $n\rho_{n}^{1}=O(1)$.

###### Assumption 9 (Variance estimation)

$\lim_{n\rightarrow\infty}\mathbb{P}(\min_{i}\big{(}\bm{M}_{n}\big{)}_{i,i}>0)=1$,

$\displaystyle\mathbb{P}\Big{(}\min_{i}\big{(}\bm{M}_{n}\big{)}_{i,i}>0\Big{)}=O_{p}(1),\quad\frac{\sum_{i=1}^{n}||\widetilde{\bm{Q}}_{i,n}||^{4}}{n}=O_{p}(1),$

and $\max_{i}||\mu_{i,n}||/\sqrt{n}=o_{p}(1)$ with
$\mu_{i,n}=\mathbb{E}\big{[}y_{i,n}|\\{\mathbf{x}_{i,n}\\}_{i=1}^{n},\\{\mathbf{w}_{i,n}\\}_{i=1}^{n}\big{]}$.

###### Assumption 10 (Policy effect sizes)

For $\delta\in(0,\frac{1}{2})$, the asymptotic distance between the effects of
policy $k$ and $j$ diverges as $n\rightarrow\infty$:

$\displaystyle
n^{\delta}\cdot\min_{k\not\in\mathcal{H}_{(j)}}\big{|}\beta_{(j)}-\beta_{k}\big{|}\rightarrow\infty,\text{
as }n\rightarrow\infty,\quad j=1,\ldots,d.$

####  A.2 Proof of Theorem 1

In this section, we show the following theorem holds:

###### Theorem 2

Under Assumptions 6-10, for any $t\in\mathbb{R}$, for the resampled
statistics, the following holds precisely

$\displaystyle\mathbb{P}\left(\frac{\sqrt{n}\big{(}\widetilde{\beta}_{(j)}^{*}-\widetilde{\beta}_{(j)}\big{)}}{(\widehat{\bm{e}}_{j}^{\prime}\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}\widehat{\bm{e}}_{j})^{\frac{1}{2}}}\leq
t\Big{|}\\{(y_{i,n},\bm{x}_{i,n}^{\prime},\bm{w}_{i,n}^{\prime})^{\prime}\\}_{i=1}^{n}\right)=\Phi(t).$

For the original statistics, it holds that

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\left(\frac{\sqrt{n}\big{(}\widetilde{\beta}_{(j)}-{\beta}_{(j)}\big{)}}{(\widehat{\bm{e}}_{j}^{\prime}\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}\widehat{\bm{e}}_{j})^{\frac{1}{2}}}\leq
t\right)=\Phi(t).$

In addition, we show that

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}\mathbb{P}\big{(}\widetilde{\beta}_{(j)}^{*}\leq\beta_{(j)}|\\{\bm{z}_{i,n}\\}_{i=1}^{n}\big{)}\leq
s\Big{)}=s.$

###### Proof.

Our proof of Theorem 1 entails the following steps:

Step 1.

Under Assumptions 6-8, [37] has shown the following holds

$\displaystyle(\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n})^{-\frac{1}{2}}\sqrt{n}(\widehat{\bm{\beta}}-\bm{\beta})\leadsto
N(0,\bm{I}_{d}),$

where $\bm{I}_{d}$ is a $d$-dimensional identity matrix. Therefore, following
the definition of $\bm{e}_{j}$ in Eq (A.3), we have

$\displaystyle({\bm{e}}_{j}^{\prime}\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)}\leadsto\mathcal{N}(0,1).$
(A.4)

Step 2.

Because of Lemma 1, we have

$\lim_{n\rightarrow\infty}\mathbb{P}(\widehat{\mathcal{H}}_{(j)}\neq\mathcal{H}_{(j)})=0.$

Combing this with (A.4), we have

$\displaystyle\Phi(t)$
$\displaystyle=\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}({\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)}\leq
t\Big{)}$
$\displaystyle=\lim_{n\rightarrow\infty}\Bigg{[}\mathbb{P}\Big{(}({\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)}\leq
t\Big{|}\widehat{\mathcal{H}}_{(j)}=\mathcal{H}_{(j)}\Big{)}\mathbb{P}(\widehat{\mathcal{H}}_{(j)}=\mathcal{H}_{(j)})+o_{\mathbb{P}}(1)\Bigg{]}$
$\displaystyle=\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}({\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)}\leq
t\Big{|}\widehat{\mathcal{H}}_{(j)}=\mathcal{H}_{(j)}\Big{)},$

in which $o_{\mathbb{P}}(1)$ is lower bounded by zero and upper bounded by
$\mathbb{P}(\widehat{\mathcal{H}}_{(j)}\neq\mathcal{H}_{(j)})$, which tends to
zero when $n$ tends to infinity. We use this same $o_{\mathbb{P}}(1)$ notion
throughout this proof. Now we have

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}(\widehat{\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}\widehat{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}\big{)}\leq
t\Big{|}\widehat{\mathcal{H}}_{(j)}=\mathcal{H}_{(j)}\Big{)}=\Phi(t).$

Next, we have

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}(\widehat{\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}\widehat{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}\big{)}\leq
t\Big{)}$ $\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\Bigg{[}\mathbb{P}\Big{(}(\widehat{\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}\widehat{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}\big{)}\leq
t\Big{|}\widehat{\mathcal{H}}_{(j)}=\mathcal{H}_{(j)}\Big{)}\Big{)}\mathbb{P}(\widehat{\mathcal{H}}_{(j)}=\mathcal{H}_{(j)})+o_{\mathbb{P}}(1)\Bigg{]}$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}(\widehat{\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}\widehat{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}\big{)}\leq
t\Big{|}\widehat{\mathcal{H}}_{(j)}=\mathcal{H}_{(j)}\Big{)}$ $\displaystyle=$
$\displaystyle\,\Phi(t)$

The following holds precisely following the definition of the resampling
procedure:

$\displaystyle\mathbb{P}\Big{(}({\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}^{*}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)}\leq
t\big{|}\\{\bm{z}_{i,n}\\}_{i=1}^{n}\Big{)}=\Phi(t).$ (A.5)

We now show that

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}(\widehat{\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}\widehat{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}^{*}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}\big{)}\leq
t\big{|}\\{\bm{z}_{i,n}\\}_{i=1}^{n}\Big{)}$ $\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\Bigg{[}\mathbb{P}\Big{(}(\widehat{\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}\widehat{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}^{*}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|}\big{)}\leq
t\big{|}\\{\bm{z}_{i,n}\\}_{i=1}^{n},\widehat{\mathcal{H}}_{(j)}=\mathcal{H}_{(j)}\Big{)}\mathbb{P}(\widehat{\mathcal{H}}_{(j)}=\mathcal{H}_{(j)})+o_{\mathbb{P}}(1)\Bigg{]}$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\Bigg{[}\mathbb{P}\Big{(}({\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}^{*}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)}\leq
t\big{|}\\{\bm{z}_{i,n}\\}_{i=1}^{n},\widehat{\mathcal{H}}_{(j)}=\mathcal{H}_{(j)}\Big{)}\mathbb{P}(\widehat{\mathcal{H}}_{(j)}=\mathcal{H}_{(j)})+o_{\mathbb{P}}(1)\Bigg{]}$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}({\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}^{*}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)}\leq
t\big{|}\\{\bm{z}_{i,n}\\}_{i=1}^{n}\Big{)}$ $\displaystyle=$
$\displaystyle\,\Phi(t).$

Recall our definition in the main manuscript

$\displaystyle\widetilde{\beta}_{(j)}^{*}=\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}^{*}_{k}}{|\widehat{\mathcal{H}}_{(j)}|},\text{
and
}\widetilde{\beta}_{(j)}=\frac{\sum_{k\in\widehat{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|\widehat{\mathcal{H}}_{(j)}|},$
(A.6)

we thus have reached the conclusion presented in the theorem:

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}(\widehat{\bm{e}}_{j}^{\prime}\widehat{\Sigma}^{\texttt{KJ}}_{n}\widehat{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\widetilde{\beta}_{(j)}^{*}-\widetilde{\beta}_{(j)}\big{)}\leq
t\big{|}\\{\bm{z}_{i,n}\\}_{i=1}^{n}\Big{)}=\Phi(t).$

Step 3.

Lastly, to prove the bootstrap consistency, we show that

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}\mathbb{P}\big{(}\widetilde{\beta}_{(j)}^{*}\leq\beta_{(j)}|\\{\bm{z}_{i,n}\\}_{i=1}^{n}\big{)}\leq
s\Big{)}=s.$

Note that

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\big{(}\widetilde{\beta}_{(j)}^{*}\leq\beta_{(j)}|\\{\bm{z}_{i,n}\\}_{i=1}^{n}\big{)}$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}({\bm{e}}_{j}^{\prime}\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}^{*}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)}$
$\displaystyle\leq$
$\displaystyle({\bm{e}}_{j}^{\prime}\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)}\big{|}\\{\bm{z}_{i,n}\\}_{i=1}^{n}\Big{)}$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\Phi(({\bm{e}}_{j}^{\prime}\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)}).$

Therefore, by the bounded convergence theorem, we have

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}\mathbb{P}\big{(}\widetilde{\beta}_{(j)}^{*}\leq\beta_{(j)}|\\{\bm{z}_{i,n}\\}_{i=1}^{n}\big{)}\leq
s\Big{)}$ $\displaystyle=$ $\displaystyle\
\mathbb{P}(\lim_{n\rightarrow\infty}\Phi(({\bm{e}}_{j}^{\prime}\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)})\leq
s)$ $\displaystyle=$ $\displaystyle\
\mathbb{P}(\lim_{n\rightarrow\infty}({\bm{e}}_{j}^{\prime}\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}{\bm{e}}_{j})^{-\frac{1}{2}}\cdot\sqrt{n}\big{(}\frac{\sum_{k\in{\mathcal{H}}_{(j)}}{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}-\frac{\sum_{k\in{\mathcal{H}}_{(j)}}\widehat{\beta}_{k}}{|{\mathcal{H}}_{(j)}|}\big{)}\leq\Phi^{-1}(s))$
$\displaystyle=$ $\displaystyle\ \mathbb{P}(N(0,1)\leq\Phi^{-1}(s))$
$\displaystyle=$ $\displaystyle\ s.$

∎

###  Appendix B Lemmas and corollary

####  B.1 Lemma 1

###### Lemma 1

Suppose $w_{k,(j)}=\mathds{1}(k\in\widehat{\mathcal{H}}_{(j)})$, for
$j,k=1,\ldots,d$, under Assumptions 6-10, we have the following argument holds
$\forall\varepsilon>0$,

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\big{(}|w_{k,(j)}-\mathds{1}(k\in\mathcal{H}_{(j)})|>\varepsilon\big{)}=0.$

###### Proof.

We start with reviewing and introducing some notations to pave the way for a
clear proof. Recall that $\widetilde{\beta}_{(j)}^{*}$ is the $j$-th largest
effect size for the resampled statistics
$\widehat{\beta}_{(1)},\ldots,\widehat{\beta}_{(d)}$, suppose
$\widetilde{\beta}_{(j)}^{*}$ is resampled statistics from the normal
distribution centered at $\widehat{\beta}_{\widecheck{j}}$, that is

$\displaystyle\widetilde{\beta}_{(j)}^{*}\big{|}\\{\bm{z}_{i,n}\\}_{i=1}^{n}\sim\mathcal{N}(\widehat{\beta}_{\widecheck{j}},(\widehat{\Sigma}^{\texttt{KJ}}_{n})_{\widecheck{j}\widecheck{j}}),\quad\widecheck{j}=\sum_{k=1}^{d}k\cdot\mathds{1}(\widehat{\beta}_{k}^{*}=\widetilde{\beta}_{(j)}^{*}),$

where $(\widehat{\Sigma}^{\texttt{KJ}}_{n})_{\widecheck{j}\widecheck{j}}$ is
the $\widecheck{j}$th component in the diagonal of the matrix
$\widehat{\Sigma}^{\texttt{KJ}}_{n}$.

Recall we define the near tie $\mathcal{H}_{(j)}$ set as

$\displaystyle\mathcal{H}_{(j)}=\big{\\{}k:|\beta_{k}-\beta_{(j)}|=o(n^{-\frac{1}{2}}),\
k=1,\ldots,d\big{\\}}.$

We further define two sets of policies that have effect sizes lower/larger
than the policies in the set $\mathcal{H}_{(j)}$:

$\displaystyle\mathcal{H}_{(j)}^{L}=\big{\\{}k:\
\beta_{k}<\min_{m\in\mathcal{H}_{(j)}}\\{\beta_{m}\\}\
k=1,\ldots,d\big{\\}},\quad\mathcal{H}_{(j)}^{U}=\big{\\{}k:\
\beta_{k}>\max_{m\in\mathcal{H}_{(j)}}\\{\beta_{m}\\}\ k=1,\ldots,d\big{\\}}.$

As for the estimated near tie set, we have for any
$j\in\widehat{\mathcal{H}}_{(j)}$ that

$\displaystyle-
b_{L}\leq\widehat{\beta}^{*}_{k}-\widetilde{\beta}^{*}_{{{(j)}}}\leq
b_{R},\quad\text{with }|b_{R}-b_{L}|=O(n^{-\delta}),$

where $\delta\in(0,0.5)$. Thus, there exists a positive constant $C$ such that

$\displaystyle\frac{|\widehat{\beta}^{*}_{k}-\widetilde{\beta}^{*}_{{{(j)}}}|}{n^{-\delta}}<C,\quad\forall
j\in\widehat{\mathcal{H}}_{(j)}.$

Our proof is composed of the following three steps:

Step 1.

We first show that the policy with the $j$th largest policy effect size in the
resampled statistics falls into the set $\mathcal{H}_{(j)}$ with high
probability, that is

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}\widecheck{j}\in\mathcal{H}_{(j)}\Big{)}=1.$
(B.1)

Because
$\widehat{\beta}_{\widecheck{j}}^{*}\in[\min_{j\in\mathcal{H}_{(j)}}\widehat{\beta}_{j}^{*},\max_{j\in\mathcal{H}_{(j)}}\widehat{\beta}_{j}^{*}]$
by definition, coupled with the fact that

$\displaystyle\Big{\\{}\max_{k\in\mathcal{H}_{(j)}^{L}}\widehat{\beta}_{k}^{*}<\min_{j\in\mathcal{H}_{(j)}}\widehat{\beta}_{j}^{*}\leq\max_{j\in\mathcal{H}_{(j)}}\widehat{\beta}_{j}^{*}<\min_{k\in\mathcal{H}_{(j)}^{U}}\widehat{\beta}_{k}^{*}\Big{\\}}\subset\Big{(}\widecheck{j}\in\mathcal{H}_{(j)}\Big{)},$

it is suffice to show

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\left(\max_{k\in\mathcal{H}_{(j)}^{L}}\widehat{\beta}_{k}^{*}<\min_{j\in\mathcal{H}_{(j)}}\widehat{\beta}_{j}^{*}\leq\max_{j\in\mathcal{H}_{(j)}}\widehat{\beta}_{j}^{*}<\min_{k\in\mathcal{H}_{(j)}^{U}}\widehat{\beta}_{k}^{*}\right)=1.$

Under Assumption 10, by Lemma 2, for any $k\in\mathcal{H}_{(j)}^{L}$ and
$m\in\mathcal{H}_{(j)}$, we have

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\left(\widehat{\beta}_{k}^{*}<\widehat{\beta}_{m}^{*}\right)=1.$

Similarly, for any $k\in\mathcal{H}_{(j)}^{U}$ and $m\in\mathcal{H}_{(j)}$, we
have

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\left(\widehat{\beta}_{m}^{*}<\widehat{\beta}_{k}^{*}\right)=1.$

Step 2.

We then show, for $k\not\in\mathcal{H}_{(j)}$

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\big{(}w_{k,(j)}>\varepsilon,k\not\in\mathcal{H}_{(j)}\big{)}=0.$

For any $\varepsilon>0$ and $k\not\in\mathcal{H}_{(j)}$, we have the following
holds

$\displaystyle\mathbb{P}\big{(}w_{k,(j)}>\varepsilon\big{)}=$ $\displaystyle\
\mathbb{P}\big{(}\mathds{1}(k\in\widehat{\mathcal{H}}_{(j)})>\varepsilon\big{)}$
$\displaystyle=$ $\displaystyle\
\mathbb{P}\big{(}\mathds{1}(k\in\widehat{\mathcal{H}}_{(j)})>\varepsilon|k\in\widehat{\mathcal{H}}_{(j)}\big{)}\cdot\mathbb{P}\big{(}k\in\widehat{\mathcal{H}}_{(j)}\big{)}$
$\displaystyle\
+\mathbb{P}\big{(}\mathds{1}(k\in\widehat{\mathcal{H}}_{(j)})>\varepsilon|k\not\in\widehat{\mathcal{H}}_{(j)}\big{)}\cdot\mathbb{P}\big{(}k\not\in\widehat{\mathcal{H}}_{(j)}\big{)}$
$\displaystyle\leq$ $\displaystyle\
\mathbb{P}\big{(}k\in\widehat{\mathcal{H}}_{(j)}\big{)}$
$\displaystyle\overset{\text{Def}}{=}$ $\displaystyle\
\mathbb{P}\Big{(}\frac{|\widehat{\beta}^{*}_{k}-\widehat{\beta}^{*}_{{{(j)}}}|}{n^{-\delta}}<C\Big{)}$
$\displaystyle=$ $\displaystyle\
\mathbb{P}\Big{(}\frac{|\widehat{\beta}^{*}_{k}-\widehat{\beta}^{*}_{\widecheck{j}}|}{n^{-\delta}}<C\Big{)}$
$\displaystyle=$ $\displaystyle\
\mathbb{P}\Big{(}|\widehat{\beta}^{*}_{k}-\widehat{\beta}^{*}_{\widecheck{j}}|<n^{-\delta}\cdot
C\Big{)}$ $\displaystyle=$ $\displaystyle\
\mathbb{P}\Big{(}|(\widehat{\beta}^{*}_{k}-\beta_{k})-(\widehat{\beta}^{*}_{\widecheck{j}}-\beta_{\widecheck{j}})+(\beta_{k}-\beta_{\widecheck{j}})|<n^{-\delta}\cdot
C\Big{)}$ $\displaystyle\leq$ $\displaystyle\
\mathbb{P}\Big{(}n^{\delta}|\beta_{k}-\beta_{\widecheck{j}}|-n^{\delta}|\widehat{\beta}^{*}_{k}-\beta_{k}|-n^{\delta}|\widehat{\beta}^{*}_{\widecheck{j}}-\beta_{\widecheck{j}}|<C\Big{)}$
$\displaystyle=$ $\displaystyle\
\mathbb{P}\Big{(}n^{\delta}|\widehat{\beta}^{*}_{k}-\beta_{k}|+n^{\delta}|\widehat{\beta}^{*}_{\widecheck{j}}-\beta_{\widecheck{j}}|>n^{\delta}|\beta_{k}-\beta_{\widecheck{j}}|+C\Big{)}$

By definition, for $k\not\in\mathcal{H}_{(j)}$

$\displaystyle\mathbb{P}\big{(}n^{\delta}|\beta_{k}-\beta_{\widecheck{j}}|<C\big{)}\leq$
$\displaystyle\
\mathbb{P}\big{(}n^{\delta}|\beta_{k}-\beta_{\widecheck{j}}|<C,\widecheck{j}\in\mathcal{H}_{(j)}\big{)}+\mathbb{P}\big{(}\widecheck{j}\not\in\mathcal{H}_{(j)}\big{)}$
$\displaystyle\leq$ $\displaystyle\ \max_{j\in\mathcal{H}_{j}}\
\mathbb{P}\big{(}n^{\delta}|\beta_{k}-\beta_{{j}}|<C,{j}\in\mathcal{H}_{(j)}\big{)}+\mathbb{P}\big{(}\widecheck{j}\not\in\mathcal{H}_{(j)}\big{)}.$

Under Assumption 10, Lemma 2 and the conclusion in Eq (B.1) in Step 1 suggest
that by letting $n\rightarrow\infty$ on both side, we have the above
probability converges to zero. Based on above derivation, we have shown that
$\lim_{n\rightarrow\infty}\mathbb{P}\big{(}w_{k,(j)}>\varepsilon,k\not\in\mathcal{H}_{(j)}\big{)}=0$,
for $k\not\in\mathcal{H}_{(j)}$.

Step 3.

We are left to prove that for all $k\in\mathcal{H}_{(j)}$, the following holds
$\forall\varepsilon>0$:

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\big{(}|w_{k,(j)}-1|>\varepsilon\big{)}=0.$

Following similar arguments, for a positive constant $C$, we have for
$k,j\in\mathcal{H}_{(j)}$ the following statement holds

$\displaystyle\mathbb{P}\big{(}|w_{k,(j)}-1|>\varepsilon\big{)}=$
$\displaystyle\
\mathbb{P}\big{(}|\mathds{1}(k\in\widehat{\mathcal{H}}_{(j)})-1|>\varepsilon\big{)}$
$\displaystyle=$ $\displaystyle\
\mathbb{P}\big{(}|\mathds{1}(k\in\widehat{\mathcal{H}}_{(j)})-1|>\varepsilon|k\in\widehat{\mathcal{H}}_{(j)}\big{)}\cdot\mathbb{P}\big{(}k\in\widehat{\mathcal{H}}_{(j)}\big{)}$
$\displaystyle\
+\mathbb{P}\big{(}|\mathds{1}(k\in\widehat{\mathcal{H}}_{(j)})-1|>\varepsilon|k\not\in\widehat{\mathcal{H}}_{(j)}\big{)}\cdot\mathbb{P}\big{(}k\not\in\widehat{\mathcal{H}}_{(j)}\big{)}$
$\displaystyle\leq$ $\displaystyle\
\mathbb{P}\big{(}k\not\in\widehat{\mathcal{H}}_{(j)}\big{)}$
$\displaystyle\overset{\text{Def}}{=}$ $\displaystyle\
\mathbb{P}\Big{(}\frac{|\widehat{\beta}^{*}_{k}-\widehat{\beta}^{*}_{{{(j)}}}|}{n^{-\delta}}\geq
C\Big{)}$ $\displaystyle=$ $\displaystyle\
\mathbb{P}\Big{(}\frac{|\widehat{\beta}^{*}_{k}-\widehat{\beta}^{*}_{\widecheck{j}}|}{n^{-\delta}}\geq
C\Big{)}$ $\displaystyle=$ $\displaystyle\
\mathbb{P}\Big{(}|(\widehat{\beta}^{*}_{k}-\beta_{k})-(\widehat{\beta}^{*}_{\widecheck{j}}-\beta_{\widecheck{j}})+(\beta_{k}-\beta_{\widecheck{j}})|\geq
n^{-\delta}\cdot C\Big{)}$ $\displaystyle\leq$ $\displaystyle\
\mathbb{P}\Big{(}|\widehat{\beta}^{*}_{k}-\beta_{k}|+|\widehat{\beta}^{*}_{\widecheck{j}}-\beta_{\widecheck{j}}|+|\beta_{k}-\beta_{\widecheck{j}}|\geq
n^{-\delta}\cdot C\Big{)}$ $\displaystyle\leq$ $\displaystyle\
\mathbb{P}\Big{(}n^{\frac{1}{2}}|\widehat{\beta}^{*}_{k}-\beta_{k}|+n^{\frac{1}{2}}|\widehat{\beta}^{*}_{\widecheck{j}}-\beta_{\widecheck{j}}|+n^{\frac{1}{2}}|\beta_{k}-\beta_{\widecheck{j}}|\geq
n^{\frac{1}{2}-\delta}\cdot C\Big{)}.$

By definition of the near-tie set, for $k\in\mathcal{H}_{(j)}$, we have

$\displaystyle\mathbb{P}\big{(}n^{\frac{1}{2}}|\beta_{k}-\beta_{\widecheck{j}}|<C\big{)}\leq$
$\displaystyle\
\mathbb{P}\big{(}n^{\delta}|\beta_{k}-\beta_{\widecheck{j}}|<C,\widecheck{j}\in\mathcal{H}_{(j)}\big{)}+\mathbb{P}\big{(}\widecheck{j}\not\in\mathcal{H}_{(j)}\big{)}$
$\displaystyle\leq$ $\displaystyle\ \max_{j\in\mathcal{H}_{j}}\
\mathbb{P}\big{(}n^{\delta}|\beta_{k}-\beta_{{j}}|<C,{j}\in\mathcal{H}_{(j)}\big{)}+\mathbb{P}\big{(}\widecheck{j}\not\in\mathcal{H}_{(j)}\big{)}.$

Again, under Assumption 10, Lemma 2 and the conclusion in Eq (B.1) we have
derived in Step 1, by letting $n\rightarrow\infty$ on both side, we have the
above probability converges to 1.

∎

####  B.2 Lemma 2

###### Lemma 2

Under Assumption 10, we show that for all $k\in\\{1,\ldots,d\\}$, any positive
constant $C$ and $\delta\in(0,\frac{1}{2})$, the following statement holds

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\big{(}|\widehat{\beta}^{*}_{k}-\beta_{k}|\geq
n^{-\delta}\cdot C\big{)}=0.$

###### Proof.

Note that

$\displaystyle\sqrt{n}(\widehat{\beta}^{*}_{k}-\beta_{k})$
$\displaystyle=\sqrt{n}(\widehat{\beta}_{k}-\beta_{k})+\mathcal{N}\big{(}0,(\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n})_{k,k}\big{)}.$
(B.2)

Because $\sqrt{n}(\widehat{\beta}_{k}-\beta_{k})$ converges in distribution to
a finite-value random variable and $\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}$
converges in probability to a finite-value matrix when $n$ tends to infinity
[37], for any given $\varepsilon>0$, there exists an $M$, such that

$\displaystyle\mathbb{P}(\sqrt{n}|\widehat{\beta}^{*}_{k}-\beta_{k}|>M)<\varepsilon.$
(B.3)

Then, for any $n$ such that

$n>\big{(}\frac{M}{C}\big{)}^{\frac{1}{\frac{1}{2}-\delta}}$

we have that

$\mathbb{P}\big{(}|\widehat{\beta}^{*}_{k}-\beta_{k}|\geq n^{-\delta}\cdot
C\big{)}<\varepsilon$

and therefore

$\limsup_{n\rightarrow\infty}\mathbb{P}\big{(}|\widehat{\beta}^{*}_{k}-\beta_{k}|\geq
n^{-\delta}\cdot C\big{)}<\varepsilon.$

Note that the above inequality holds for arbitrary $\varepsilon>0$. Therefore,
we have

$\limsup_{n\rightarrow\infty}\mathbb{P}\big{(}|\widehat{\beta}^{*}_{k}-\beta_{k}|\geq
n^{-\delta}\cdot C\big{)}=0,$

completing the proof. ∎

###### Lemma 3

Denote the selected policy as

$\displaystyle\widehat{j}=\sum_{k=1}^{d}k\cdot\mathds{1}(\widehat{\beta}_{k}=\widehat{\beta}_{(j)}),$

we show that

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\big{(}\widehat{j}\in\mathcal{H}_{(j)}\big{)}=1.$

###### Proof.

This is a direct result from Step 1 and Step 2 in the proof for Theorem 1. ∎

####  B.3 Corollary 1

###### Corollary 2

Under Assumptions 6-10, we have that
$\lim_{n\rightarrow\infty}\mathbb{P}\Big{(}\mathbb{P}\big{(}\widetilde{\beta}_{(j)}^{*}\leq\beta_{\widehat{j}}|\\{\bm{z}_{i,n}\\}_{i=1}^{n}\big{)}\leq
s\Big{)}=s$.

###### Proof.

Because of the consistency in Lemma 3, we have

$\lim_{n\rightarrow\infty}\mathbb{P}(\beta_{\widehat{j}}=\beta_{(j)})=1.$

Therefore,

$\displaystyle\lim_{n\rightarrow\infty}|\mathbb{P}\big{(}\widetilde{\beta}_{(j)}^{*}\leq\beta_{(j)}|\\{\bm{z}_{i,n}\\}_{i=1}^{n}\big{)}-\mathbb{P}\big{(}\widetilde{\beta}_{(j)}^{*}\leq\beta_{(j)}|\\{\bm{z}_{i,n}\\}_{i=1}^{n}\big{)}|\leq\lim_{n\rightarrow\infty}\mathbb{P}(\beta_{(j)}\neq\widehat{\beta}_{(j)}|\\{\bm{z}_{i,n}\\}_{i=1}^{n})=0.$

The result then follows by applying Step 3 in the Proof of Theorem 1. ∎

###  Appendix C Additional simulation and empirical results

####  C.1 Practical implementation

In this section, we discuss the choice of tuning parameters (including $B$,
$\delta$, $b_{\texttt{L}}$ and $b_{\texttt{R}}$) of the proposed method in
Section 2.2. For the number of repetitions for our resampling procedure, we
recommend using $B=2,000$ as a good balance between computational load and
statistical inference accuracy.

For the tuning pair $(b_{\texttt{R}},b_{\texttt{L}})$, from our theoretical
analysis, we need to ensure that the distance between $b_{\texttt{R}}$ and
$b_{\texttt{L}}$ is of the order $n^{\delta}$ with $\delta\in(0,0.5)$ to
guarantee the statistical validity of our proposed procedure. To achieve this
goal, for the policy $\beta_{(j)}$, we adopt the tuning pair of the form

$\displaystyle b_{\texttt{L}}^{j}=n^{-\delta}\cdot
s_{\widehat{j}}^{\delta}\cdot c_{\texttt{L}}^{j},\quad
b_{\texttt{R}}^{j}=n^{-\delta}\cdot s_{\widehat{j}}^{\delta}\cdot
c_{\texttt{R}}^{j},$

where $s_{\widehat{j}}$ is the $\widehat{j}$th element in the diagonal of the
estimated covariance matrix $\widehat{\bm{\Omega}}^{\texttt{KJ}}_{n}$,
$c_{\texttt{L}}^{j}$ and $c_{\texttt{R}}^{j}$ are positive constants.

The constants $c_{\texttt{L}}^{j}$ and $c_{\texttt{R}}^{j}$ can significantly
impact the performance of the proposed approach in finite samples. In the
extreme cases, on the one hand, if both $c_{\texttt{L}}^{j}$ and
$c_{\texttt{R}}^{j}$ are overly large, the estimated near tie set might
include more policies than necessary and our approach is only valid if all
true policy effects are closely ties. On the other hand, if
$c_{\texttt{L}}^{j}$ and $c_{\texttt{R}}^{j}$ are both closer to zero, our
approach reduces to a standard parametric bootstrap approach, which is
problematic in the presence of tied policy effects. To present a robust
algorithm in finite samples, we thus adopt the following “double-bootstrap”
method as discussed in [17] (note that double-bootstrapped statistics are
labelled with double-star superscripts):

  1. 1.

For $j=1,...,d$, set
$\beta_{j}^{*}=\Delta\cdot\frac{\sum_{j=1}^{d}\widehat{\beta}_{j}}{d}+(1-\Delta)\cdot\widehat{\beta}_{j}$,
where

$\Delta={\min}\Big{\\{}1,\frac{\sum_{j=1}^{d}s_{\widehat{j}}}{n\sum_{j=1}^{d}(\widehat{\beta}_{j}-\bar{\widehat{\beta}})^{2}}\times
n^{0.05}\Big{\\}}.$

  2. 2.

For every candidate pair $(c_{\texttt{L}},c_{\texttt{R}})$ such that
$c_{\texttt{L}}\in\mathcal{C}_{\texttt{L}}$ and
$c_{\texttt{R}}\in\mathcal{C}_{\texttt{R}}$, do

    1. (a)

For $t\leftarrow 1$ to $T$, do

      1. i.

Generate
$\widehat{\bm{\beta}}^{*}=(\widehat{\beta}_{1}^{*},...,\widehat{\beta}_{d}^{*})$
from $\mathcal{N}(\bm{\beta}^{*},\widehat{\Omega}_{n}^{\texttt{KJ}}/n)$, where
$\bm{\beta}^{*}=(\beta_{1}^{*},\ldots,\beta^{*}_{d})^{\prime}$, and denote the
ordered values in $\bm{\beta}^{*}$ as
$\beta_{(1)}^{*}\geq\ldots\geq\beta_{(d)}^{*}$.

      2. ii.

For $r\leftarrow 1$ to $R$, do

        1. A.

Generate double bootstrap statistics
$\widehat{\bm{\beta}}^{**}\triangleq(\widehat{\beta}_{1}^{**},...,\widehat{\beta}_{d}^{**})^{\prime}$
from
$\mathcal{N}(\widehat{\bm{\beta}}^{*},\widehat{\Omega}_{n}^{\texttt{KJ}}/n)$,
and denote the ordered values of $\widehat{\bm{\beta}}^{**}$ as
$\widehat{\beta}_{(1)}^{**}\geq...\geq\widehat{\beta}_{(d)}^{**}$.

        2. B.

Record $w_{k,(j)}^{**}=\mathds{1}\\{-c_{\texttt{L}}\cdot n^{-\delta}\cdot
s_{\widehat{j}}^{\delta}\leq(\widehat{\beta}_{k}^{**}-\beta_{(j)}^{**})\leq
c_{\texttt{R}}\cdot n^{-\delta}\cdot s_{\widehat{j}}^{\delta}\\}$ and
$\widetilde{\beta}^{**}_{(j)}=\sum_{k=1}^{d}w_{k,(j)}^{**}\widehat{\beta}_{j}^{**}/\sum_{k=1}^{d}w_{k,(j)}^{**}$,
for $j=1,\ldots,d$.

      3. iii.

Calculate
${B}_{j,t}(c_{\texttt{L}},c_{\texttt{R}})=\frac{1}{R}\sum_{r=1}^{R}\mathds{1}\big{(}\widetilde{\beta}^{**,r}_{(j)}\leq\beta^{*,r}_{(j)}\big{)}$,
for $j=1,\ldots,d$.

  3. 3.

Record the loss function

$\displaystyle{L}_{j}(c_{\texttt{L}},c_{\texttt{R}})=\frac{1}{T}\sum_{t=1}^{T}\Big{(}{B}_{j,(t)}(c_{\texttt{L}},c_{\texttt{R}})-\frac{t}{T+1}\Big{)}^{2},$
(C.1)

where ${B}_{j,(t)}(c_{\texttt{L}},c_{\texttt{R}})$ is the $t$-th smallest
statistics in
${B}_{j,1}(c_{\texttt{L}},c_{\texttt{R}}),\ldots,{B}_{j,T}(c_{\texttt{L}},c_{\texttt{R}})$.

  4. 4.

Choose the pair $(c_{\texttt{L}}^{j},c_{\texttt{R}}^{j})$ for inferring
$\beta_{(j)}$ and $\beta_{\widehat{j}}$ that minimizes
${L}_{j}(c_{\texttt{L}},c_{\texttt{R}})$, that is

$\displaystyle(c_{\texttt{L}}^{j},c_{\texttt{R}}^{j})=\underset{(c_{\texttt{L}},c_{\texttt{R}})\in\mathcal{C}}{\min}{L}_{j}(c_{\texttt{L}},c_{\texttt{R}}),\quad
j=1,\ldots,d_{0}.$

Note that we only use the above procedure to choose the tuning parameters
$c_{\texttt{L}}^{j}$ and $c_{\texttt{R}}^{j}$, meaning that we do not use the
resampled statistics in Step 1 to carry out inference on $\beta_{(j)}$. In
Step 1, $\Delta$ is adopted to stabilize the performance of the tuning
parameter selection in finite samples, and $\Delta$ only takes a close-to-zero
value whenever limited variation is found between policy effect estimates.

Following Theorem 1, we know that
$\mathbb{P}\big{(}\widetilde{\beta}_{(j)}^{*}\leq\beta_{(j)}|\\{\bm{z}_{i,n}\\}_{i=1}^{n}\big{)}$
roughly follows $\text{Unif}(0,1)$ when the sample size $n$ is large. Given a
desirable tuning pair $(c_{\texttt{L}},c_{\texttt{R}})$, we would thus expect
that
${B}_{j,(1)}(c_{\texttt{L}},c_{\texttt{R}}),...,{B}_{j,(T)}(c_{\texttt{L}},c_{\texttt{R}})$
share a similar distribution with the ordered statistics of i.i.d.
$\text{Unif}(0,1)$ random variables. The loss function defined in Eq (C.1)
measures the average of squared differences between
${B}_{j,(t)}(c_{\texttt{L}},c_{\texttt{R}})$ and the expected value of the
order statistics of the $\text{Unif}(0,1)$ random variables. Given the
rational above, we would expect that the optimal tuning pair
$(c_{\texttt{L}}^{j},c_{\texttt{R}}^{j})$ minimize such a loss.

We further comment on several implementation details. Our numerical results
suggest using $R=200$ and $T=40$ can provide reasonable choice of the tuning
parameters in finite samples. In addition, when the loss function
${L}_{j}(c_{\texttt{L}},c_{\texttt{R}})$ do not fluctuate substantially over
all considered pairs $(c_{\texttt{L}},c_{\texttt{R}})$. In this case, let
$\gamma$ denote the 97.5th percentile of
${L}_{U}=\frac{1}{T}\sum_{t=1}^{T}({U}_{(b)}-\frac{t}{T+1})^{2}$ and
${U}_{(1)},...,{U}_{(T)}$ are ordered observations from $\text{Unif}(0,1)$
distribution, we choose $(\bar{c}_{L}^{j},\bar{c}_{R}^{j})$ which is the mean
of all plausible pairs such that
${L}_{j}(c_{\texttt{L}},c_{\texttt{R}})<\gamma$. Lastly, as for the candidate
region of $c_{\texttt{L}}$ and $c_{\texttt{R}}$, we first consider selecting
$c_{L}$ from 0 to
$\frac{2(\widehat{\beta}_{(1)}-\widehat{\beta}_{(j)})n^{\delta}}{s_{\widehat{j}}^{\delta}}$
and $c_{R}$ from 0 to
$\frac{2(\widehat{\beta}_{(j)}-\widehat{\beta}_{(d)})n^{\delta}}{s_{\widehat{j}}\delta}$.
Then based on the values of $L_{j}(c_{\texttt{L}},c_{\texttt{R}})$ for
different tuning pairs, we may choose to expand or shrink the candidate region
to make our algorithm more efficient.

####  C.2 Robustness to different tuning choice

We summarizing our simulation results with different choices of the tuning
parameter $\delta\in\\{0.05,0.15,0.25\\}$, $R\in\\{200,500\\}$ and
$T\in\\{40,100,200\\}$. To avoid redundancy, we showcase the results with
$\bm{\beta}=0$ and $\beta_{j}=\Phi^{-1}\big{(}\frac{j}{d+1}\big{)},\
j=1,\ldots,d$ while $q_{n}$ takes value $141$ or $561$. We report the coverage
probabilities and asymptotic biases for estimating the top two policies (i.e.,
$d_{0}=2$) $\beta_{(1)}$ and $\beta_{(2)}$. Supplementary Materials Table C.1
summarizes the simulation results with different choice of $\delta$ and fixed
$R=200$ and $T=40$. There, we observe that the performance of our method is
overall robust to the choice of different $\delta$ in a variety of settings.
Though when $q_{n}$ is large and no policy is effective, smaller $\delta$
likely leads to under-covered confidence intervals. Supplementary Materials
Table C.2 summarizes the simulation results under $R\in\\{200,500\\}$ and
$T\in\\{40,100,200\\}$, while fixing $\delta=0.25$. Our results demonstrate
that when $R$ or $T$ increases, the coverage probabilities are slightly
increased and biases are marginally reduced. Overall, we observe that the
proposed method is not very sensitive to the choice of various tuning
parameters $\delta$, $T$, and $R$. To guide readers for the selection of
tuning parameters to reach an optimal accuracy and computational efficiency
trade-off, we further provide the computational time under various choices of
$T$ and $R$ in the Supplementary Materials Table C.3. In practice, to reduce
computational cost while maintaining valid statistical inference, we adopt the
following tuning set in the rest of the numerical studies: $R=200$, $T=100$,
and $\delta=0.25$.

Table C.1: Coverage probability and asymptotic bias with different choices of
$\delta$

No policy is effective, $\beta_{(1)}=\beta_{(2)}=0$ $q_{n}=141$ $q_{n}=561$
$\delta$=0.05 $\delta$=0.15 $\delta$=0.25 $\delta=0.05$ $\delta$=0.15
$\delta$=0.25 $\beta_{(1)}$ Cover 0.97(0.00) 0.95(0.01) 0.96(0.01) 0.93(0.01)
0.92(0.01) 0.96(0.01) $\sqrt{n}$Bias 0.01(0.02) 0.03(0.04) -0.02(0.02)
0.03(0.01) 0.03(0.01) -0.01(0.01) $\beta_{(2)}$ Cover 0.98(0.00) 0.94(0.01)
0.97(0.01) 0.93 (0.01) 0.93(0.01) 0.96(0.01) $\sqrt{n}$Bias 0.00(0.01)
0.02(0.01) 0.01(0.01) 0.03(0.01) 0.02(0.01) 0.01(0.01) Top two policies are
effective, $\beta_{(1)}=0.97,\beta_{(2)}=0.43$ $q_{n}=141$ $q_{n}=561$
$\delta$=0.05 $\delta$=0.15 $\delta$=0.25 $\delta=0.05$ $\delta$=0.15
$\delta$=0.25 $\beta_{(1)}$ Cover 0.94(0.01) 0.98(0.00) 0.95(0.01) 0.93(0.01)
0.95(0.01) 0.94(0.01) $\sqrt{n}$Bias 0.05(0.05) -0.03(0.04) -0.04(0.05)
-0.07(0.06) -0.06(0.07) -0.07(0.07) $\beta_{(2)}$ Cover 0.94(0.01) 0.97(0.01)
0.95(0.01) 0.95(0.01) 0.96(0.01) 0.95(0.01) $\sqrt{n}$Bias 0.08(0.08)
0.05(0.06) 0.06(0.06) 0.08(0.08) 0.07(0.08) 0.08(0.09)

  * •

Note: We fix $R=200$ and $T=40$. “Cover” is the empirical coverage of the 95%
confidence interval for $\beta_{(j)}$ and “ $\sqrt{n}$Bias ” captures the
root-$n$ scaled Monte Carlo bias for estimating $\beta_{(j)}$. Monte Carlo
standard errors are provided in the parenthesis.

####  C.3 Computational time with different tuning parameters

In this section, we summarize computational time with respect to various
choices of $R$, $T$, and $n$ to make the computational costs transparent for
readers. Here, we fix $\delta$ at $0.25$ and select $20$ candidate tuning
parameters. Table C.3 demonstrates that the computational costs are largely
determined by $T$ and $R$. When both $T$ and $R$ reach $500$, running our
method once takes approximately one hour. For simulation study with multiple
iterations, we recommend setting $R=200$ and $T\leq 200$ to achieve a
reasonable trade-off between accuracy and computational efficiency.

Table C.2: Coverage probability and asymptotic bias with different choices of
$T$ and $R$

No policy is effective, $\beta_{(1)}=\beta_{(2)}=0$ $q_{n}=141$ $q_{n}=561$
$T=40$ $T=100$ $T=200$ $T=40$ $T=100$ $T=200$ $\beta_{(1)}$ $R=200$ Cover
0.96(0.01) 0.96(0.01) 0.96(0.01) 0.96(0.01) 0.96(0.02) 0.96(0.02)
$\sqrt{n}$Bias -0.02(0.02) 0.03(0.06) 0.02(0.06) -0.01(0.01) -0.01(0.12)
0.01(0.12) $R=500$ Cover 0.96(0.01) 0.95(0.01) 0.96(0.01) 0.97(0.01)
0.96(0.01) 0.96(0.01) $\sqrt{n}$Bias 0.01(0.02) -0.02(0.04) 0.01(0.02)
0.01(0.01) 0.01(0.01) 0.01(0.01) $\beta_{(2)}$ $R=200$ Cover 0.96(0.01)
0.96(0.01) 0.97(0.01) 0.96(0.01) 0.96(0.01) 0.97(0.01) $\sqrt{n}$Bias
0.01(0.01) -0.01(0.01) 0.01(0.01) 0.01(0.01) 0.01(0.01) -0.01(0.01) $R=500$
Cover 0.95(0.01) 0.95(0.01) 0.96(0.01) 0.96(0.01) 0.95(0.01) 0.96(0.01)
$\sqrt{n}$Bias 0.01(0.01) 0.01(0.01) 0.01(0.01) 0.01(0.01) 0.01(0.01)
0.00(0.01) Top two policies are effective, $\beta_{(1)}=0.97,\beta_{(2)}=0.43$
$q_{n}=141$ $q_{n}=561$ $T=40$ $T=100$ $T=200$ $T=40$ $T=100$ $T=200$
$\beta_{(1)}$ $R=200$ Cover 0.95(0.01) 0.95(0.01) 0.96(0.01) 0.94(0.01)
0.95(0.01) 0.95(0.01) $\sqrt{n}$Bias -0.04(0.05) 0.02(0.08) 0.02(0.08)
-0.07(0.07) 0.07(0.16) 0.06(0.18) $R=500$ Cover 0.96(0.01) 0.96(0.01)
0.96(0.01) 0.94(0.01) 0.95(0.01) 0.95(0.01) $\sqrt{n}$Bias 0.02(0.04)
0.01(0.04) 0.01(0.05) 0.05(0.07) 0.05(0.07) 0.04(0.07) $\beta_{(2)}$ $R=200$
Cover 0.95(0.01) 0.95(0.01) 0.95(0.01) 0.95(0.01) 0.95(0.01) 0.96(0.01)
$\sqrt{n}$Bias 0.06(0.06) 0.05(0.06) 0.05(0.06) 0.08(0.09) 0.07(0.08)
0.07(0.09) $R=500$ Cover 0.96(0.01) 0.96(0.01) 0.96(0.01) 0.95(0.01)
0.96(0.01) 0.96(0.01) $\sqrt{n}$Bias 0.06(0.08) 0.05(0.06) 0.05(0.06)
0.07(0.08) 0.07(0.08) 0.06(0.09)

  * •

Note: We fix $\delta=0.25$. “Cover” is the empirical coverage of the 95%
confidence interval for $\beta_{(j)}$ and “ $\sqrt{n}$Bias ” captures the
root-$n$ scaled Monte Carlo bias for estimating $\beta_{(j)}$. Monte Carlo
standard errors are provided in the parenthesis.

Computational time with respect to various $n$, $q_{n}$, $T$, and $R$
$q_{n}=141$ ($s\times 10^{3}$) $q_{n}=561$ ($s\times 10^{3}$) $T=40$ $T=200$
$T=500$ $T=40$ $T=200$ $T=500$ $n=500$ $R=200$ 0.10 0.48 1.27 0.11 0.48 1.35
$R=500$ 0.23 1.20 3.31 0.24 1.38 3.50 $n=2000$ $R=200$ 0.10 0.50 1.29 0.12
0.51 1.40 $R=500$ 0.26 1.23 3.44 0.27 1.39 3.98 $n=5000$ $R=200$ 0.11 0.52
1.32 0.14 0.55 1.51 $R=500$ 0.27 1.25 3.50 0.30 1.40 4.05

Table C.3:  The unit: 1,000 seconds. We fix $\delta=0.25$ and set $20$
candidate tuning pairs for $(c_{\texttt{L}},c_{\texttt{R}})$. The simulations
are performed on a Lenovo NeXtScale nx360m5 node (24 cores per node) equipped
with Intel Xeon Haswell processor. The core frequency is 2.3 Ghz and supports
16 floating-point operations per clock period.

####  C.4 Simulation results: $\beta_{\widehat{j}}$

The simulation results presented in Table C.4 and C.5 help us confirm our
theoretical analyses in Corollary 1, and we observe similar trends compared to
the results in the main manuscript.

Table C.4:  Simulation results
($d=5,\text{heterogeneity},\beta_{\widehat{1}}$)

$\beta_{j}=\Phi^{-1}\big{(}\frac{j}{d+1}\big{)},\quad\bm{\gamma}_{n}=0,\quad
j=1,\ldots,d$
$\bm{x}_{i,n}\sim\mathcal{N}(0,\Sigma),\quad\bm{w}_{i,n}=\mathds{1}(\widetilde{\bm{w}}_{i,n}\geq\Phi^{-1}(0.98))$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.96(0.01) 0.96(0.01) 0.95(0.01) 0.96(0.01) 0.95(0.01) $\sqrt{n}$Bias
0.04(0.06) 0.04(0.05) 0.04(0.05) -0.04(0.05) 0.04(0.06) $q_{n}=141$ Cover
0.95(0.01) 0.95(0.01) 0.93(0.01) 0.95(0.01) 0.95(0.01) $\sqrt{n}$Bias
0.07(0.07) 0.06(0.07) 0.07(0.06) 0.06(0.06) 0.05(0.07) $q_{n}=281$ Cover
0.94(0.01) 0.95(0.01) 0.84(0.01) 0.82(0.01) 0.94(0.01) $\sqrt{n}$Bias
-0.09(0.08) -0.07(0.07) -0.10(0.08) -0.11(0.07) -0.07(0.08) $q_{n}=421$ Cover
0.94(0.01) 0.91(0.01) 0.76(0.02) 0.75(0.01) 0.93(0.01) $\sqrt{n}$Bias
-0.09(0.10) -0.10(0.09) -0.15(0.09) -0.16(0.09) -0.10(0.09) $q_{n}=561$ Cover
0.94(0.01) 0.90(0.01) 0.67(0.02) 0.65(0.01) 0.78(0.01) $\sqrt{n}$Bias
-0.15(0.14) -0.12(0.10) -0.17(0.23) -0.25(0.12) 0.15(0.11) $q_{n}=631^{*}$
Cover 0.92(0.01) 0.89(0.02) 0.45(0.02) 0.42(0.01) 0.68(0.01) $\sqrt{n}$Bias
-0.19(0.18) -0.22(0.13) -0.35(0.29) -0.54(0.22) 0.28(0.18)
$\bm{x}_{i,n}=\mathds{1}(\widetilde{\bm{x}}_{i,n}>0),\quad\bm{w}_{i,n}\sim\mathcal{N}(0,I)$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.97(0.01) 0.97(0.01) 0.95(0.01) 0.95(0.01) 0.96(0.01) $\sqrt{n}$Bias
-0.05(0.10) -0.06(0.09) -0.06(0.10) -0.10(0.11) -0.04(0.09) $q_{n}=141$ Cover
0.97(0.01) 0.95(0.01) 0.94(0.01) 0.93(0.01) 0.95(0.01) $\sqrt{n}$Bias
-0.06(0.11) -0.08(0.12) -0.07(0.11) 0.13(0.12) 0.09(0.12) $q_{n}=281$ Cover
0.96(0.01) 0.94(0.01) 0.86(0.02) 0.85(0.01) 0.95(0.01) $\sqrt{n}$Bias
-0.09(0.13) -0.10(0.13) -0.10(0.13) -0.15(0.12) -0.09(0.13) $q_{n}=421$ Cover
0.94(0.01) 0.93(0.01) 0.75(0.02) 0.72(0.01) 0.93(0.01) $\sqrt{n}$Bias
0.11(0.17) -0.12(0.13) 0.18(0.17) -0.20(0.17) 0.14(0.14) $q_{n}=561$ Cover
0.94(0.01) 0.90(0.01) 0.51(0.02) 0.48(0.01) 0.92(0.01) $\sqrt{n}$Bias
-0.15(0.22) -0.21(0.20) -0.25(0.23) -0.46(0.26) -0.21(0.20) $q_{n}=631$ Cover
0.91(0.01) 0.90(0.01) 0.48(0.02) 0.45(0.01) 0.80(0.01) $\sqrt{n}$Bias
-0.21(0.20) -0.23(0.22) 0.41(0.30) -0.53(0.20) 0.35(0.22)

  * •

Note: “Cover” is the empirical coverage of the 95% confidence interval for
$\beta_{\widehat{1}}$ and “ $\sqrt{n}$Bias ” captures the root-$n$ scaled
Monte Carlo bias for estimating $\beta_{\widehat{1}}$. “ * ” indicates that
$\widehat{\bm{\Omega}}_{n}^{\texttt{KJ}}$ is not positive semi-definite in
some Monte Carlo samples.

Table C.5:  Simulation results
($d=5,\text{heterogeneity},\beta_{\widehat{2}}$)

$\bm{\beta}_{j}=\Phi^{-1}\big{(}\frac{j}{d+1}\big{)},\quad\bm{\gamma}_{n}=0,\quad
j=1,\ldots,d$
$\bm{x}_{i,n}\sim\mathcal{N}(0,\Sigma),\quad\bm{w}_{i,n}=\mathds{1}(\widetilde{\bm{w}}_{i,n}\geq\Phi^{-1}(0.98))$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.96(0.01) 0.97(0.01) 0.97(0.01) 0.96(0.01) 0.97(0.01) $\sqrt{n}$Bias
-0.04(0.06) -0.03(0.06) -0.03(0.06) -0.04(0.07) -0.04(0.06) $q_{n}=141$ Cover
0.96(0.01) 0.96(0.02) 0.88(0.02) 0.90(0.01) 0.94(0.01) $\sqrt{n}$Bias
-0.05(0.08) -0.05(0.08) -0.10(0.08) -0.09(0.06) -0.07(0.07) $q_{n}=281$ Cover
0.95(0.01) 0.94(0.02) 0.86(0.02) 0.84(0.01) 0.91(0.02) $\sqrt{n}$Bias
0.07(0.09) 0.07(0.09) 0.13(0.08) -0.15(0.10) 0.12(0.10) $q_{n}=421$ Cover
0.94(0.01) 0.93(0.02) 0.77(0.02) 0.72(0.02) 0.71(0.02) $\sqrt{n}$Bias
-0.10(0.13) -0.12(0.13) -0.15(0.11) -0.17(0.13) -0.19(0.17) $q_{n}=561$ Cover
0.94(0.01) 0.92(0.02) 0.65(0.02) 0.60(0.01) 0.69(0.02) $\sqrt{n}$Bias
-0.16(0.17) -0.18(0.16) -0.18(0.15) 0.20(0.13) 0.35(0.22) $q_{n}=631^{*}$
Cover 0.93(0.01) 0.92(0.02) 0.44(0.02) 0.42(0.01) 0.50(0.02) $\sqrt{n}$Bias
-0.18(0.17) -0.23(0.21) -0.45(0.19) -0.49(0.17) 0.48(0.30)
$\bm{x}_{i,n}=\mathds{1}(\widetilde{\bm{x}}_{i,n}>0),\quad\bm{w}_{i,n}\sim\mathcal{N}(0,I)$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.96(0.01) 0.98(0.01) 0.96(0.01) 0.96(0.01) 0.97(0.01) $\sqrt{n}$Bias
-0.08(0.13) -0.07(0.13) -0.10(0.12) 0.09(0.12) -0.10(0.10) $q_{n}=141$ Cover
0.96(0.01) 0.96(0.01) 0.94(0.01) 0.95(0.01) 0.95(0.01) $\sqrt{n}$Bias
0.09(0.14) 0.10(0.14) 0.13(0.14) 0.12(0.13) 0.10(0.14) $q_{n}=281$ Cover
0.95(0.01) 0.95(0.01) 0.92(0.01) 0.90(0.01) 0.95(0.01) $\sqrt{n}$Bias
-0.11(0.16) -0.11(0.15) 0.17(0.15) 0.19(0.14) 0.13(0.17) $q_{n}=421$ Cover
0.95(0.01) 0.92(0.02) 0.82(0.02) 0.78(0.01) 0.94(0.01) $\sqrt{n}$Bias
0.18(0.20) -0.20(0.20) -0.25(0.20) -0.32(0.17) -0.20(0.20) $q_{n}=561$ Cover
0.94(0.01) 0.92(0.02) 0.67(0.02) 0.64(0.01) 0.73(0.01) $\sqrt{n}$Bias
-0.26(0.24) 0.22(0.24) -0.41(0.23) -0.48(0.23) 0.27(0.25) $q_{n}=631$ Cover
0.92(0.01) 0.87(0.02) 0.52(0.02) 0.48(0.01) 0.50(0.01) $\sqrt{n}$Bias
0.30(0.27) 0.41(0.29) 0.50(0.28) -0.65(0.13) 0.55(0.28)

  * •

Note: “Cover” is the empirical coverage of the 95% confidence interval for
$\beta_{\widehat{2}}$ and “ $\sqrt{n}$Bias ” captures the root-$n$ scaled
Monte Carlo bias for estimating $\beta_{\widehat{2}}$. “ * ” indicates that
$\widehat{\bm{\Omega}}_{n}^{\texttt{KJ}}$ is not positive semi-definite in
some Monte Carlo samples.

####  C.5 Simulation results: $d=10$

This section provides an additional set of simulation results when $d=10$,
which is larger than the setting ($d=5$) adopted in the main manuscript. We
investigate the performance of the five methods for estimating $\beta_{(2)}$,
$\beta_{(5)}$, and $\beta_{(10)}$. Table C.6 \- C.9 demonstrate that without
adjustment, the coverage probabilities for $\beta_{(5)}$ fall below $80\%$
when $q_{n}\geq 281$, while our proposed method reaches nominal level coverage
regardless the ranking of $\beta_{j}$. “Proposed + EW”, “Proposed + HCK”, and
“Proposed + HC3” show similar trends compared to those when $d=5$.

Table C.6:  Simulation results ($d=10,\text{heterogeneity},\beta_{(2)}$)

$\beta_{j}=\Phi^{-1}\big{(}\frac{j}{d+1}\big{)},\quad\bm{\gamma}_{n}=0$, j =
1, …, d
$\bm{x}_{i,n}\sim\mathcal{N}(0,\Sigma),\quad\bm{w}_{i,n}=\mathds{1}(\widetilde{\bm{w}}_{i,n}\geq\Phi^{-1}(0.98))$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.96(0.01) 0.96(0.01) 0.96(0.01) 0.97(0.01) 0.94(0.01) $\sqrt{n}$Bias
-0.05(0.07) -0.04(0.07) -0.05(0.06) -0.06(0.07) 0.06(0.06) $q_{n}=141$ Cover
0.96(0.01) 0.95(0.01) 0.96(0.01) 0.96(0.01) 0.94(0.01) $\sqrt{n}$Bias
-0.06(0.07) -0.07(0.07) 0.07(0.08) -0.07(0.07) -0.06(0.06) $q_{n}=281$ Cover
0.95(0.01) 0.95(0.01) 0.92(0.01) 0.89(0.01) 0.92(0.01) $\sqrt{n}$Bias
-0.07(0.09) -0.08(0.09) 0.15(0.10) -0.18(0.08) 0.20(0.08) $q_{n}=421$ Cover
0.96(0.01) 0.94(0.01) 0.88(0.01) 0.80(0.02) 0.90(0.01) $\sqrt{n}$Bias
-0.08(0.11) -0.10(0.11) 0.20(0.11) -0.25(0.10) -0.25(0.09) $q_{n}=561$ Cover
0.95(0.01) 0.94(0.01) 0.75(0.01) 0.67(0.02) 0.89(0.01) $\sqrt{n}$Bias
0.14(0.15) 0.15(0.15) -0.27(0.16) -0.31(0.13) -0.28(0.12) $q_{n}=631^{*}$
Cover 0.93(0.01) 0.91(0.01) 0.65(0.01) 0.57(0.02) 0.88(0.01) $\sqrt{n}$Bias
0.22(0.20) 0.27(0.19) -0.33(0.20) 0.35(0.17) 0.38(0.16)
$\bm{x}_{i,n}=\mathds{1}(\widetilde{\bm{x}}_{i,n}>0),\quad\bm{w}_{i,n}\sim\mathcal{N}(0,I)$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.97(0.01) 0.96(0.01) 0.96(0.01) 0.96(0.01) 0.95(0.01) $\sqrt{n}$Bias
-0.05(0.11) -0.07(0.11) 0.06(0.12) -0.08(0.11) -0.09(0.09) $q_{n}=141$ Cover
0.96(0.01) 0.95(0.01) 0.94(0.01) 0.93(0.01) 0.94(0.01) $\sqrt{n}$Bias
0.07(0.12) 0.10(0.13) 0.13(0.13) 0.15(0.13) 0.10(0.11) $q_{n}=281$ Cover
0.95(0.01) 0.95(0.01) 0.91(0.01) 0.88(0.01) 0.93(0.01) $\sqrt{n}$Bias
0.14(0.15) 0.16(0.16) 0.20(0.13) -0.28(0.15) -0.29(0.13) $q_{n}=421$ Cover
0.94(0.01) 0.93 (0.01) 0.85(0.01) 0.76(0.01) 0.92(0.01) $\sqrt{n}$Bias
0.20(0.19) 0.22(0.18) 0.25(0.17) 0.37(0.18) 0.32(0.16) $q_{n}=561$ Cover
0.93(0.01) 0.93(0.01) 0.77(0.01) 0.67(0.02) 0.90(0.01) $\sqrt{n}$Bias
0.22(0.23) 0.25(0.23) -0.30(0.19) 0.41(0.22) 0.44(0.19) $q_{n}=631$ Cover
0.90(0.01) 0.88(0.01) 0.67(0.01) 0.60(0.02) 0.88(0.01) $\sqrt{n}$Bias
0.31(0.25) 0.45(0.26) -0.48(0.13) 0.55(0.25) 0.61(0.23)

  * •

Note: “Cover” is the empirical coverage of the 95% confidence interval for
$\beta_{(2)}$ and “ $\sqrt{n}$Bias ” captures the root-$n$ scaled Monte Carlo
bias for estimating $\beta_{(2)}$. “∗” indicates that
$\widehat{\bm{\Omega}}_{n}^{\texttt{KJ}}$ is not positive semi-definite in
some Monte Carlo samples.

Table C.7:  Simulation results ($d=10,\text{homogeneity},\beta_{(5)}$)

$\bm{\beta}=0,\quad\bm{\gamma}_{j}=1/j$, $j=1,\ldots,q_{n}$
$\bm{x}_{i,n}\sim\mathcal{N}(0,\Sigma),\quad\bm{w}_{i,n}=\mathds{1}(\widetilde{\bm{w}}_{i,n}\geq\Phi^{-1}(0.98))$
Proposed+KJ Proposed+HCK Proposed+HC3 Proposed+EW No adjustment+KJ $q_{n}=1$
Cover 0.95(0.01) 0.95(0.01) 0.95(0.01) 0.93(0.01) 0.92(0.01) $\sqrt{n}$Bias
0.01(0.01) 0.01(0.01) 0.01(0.01) 0.01(0.01) 0.04(0.02) $q_{n}=141$ Cover
0.94(0.01) 0.94(0.01) 0.94(0.01) 0.93(0.01) 0.82(0.01) $\sqrt{n}$Bias
-0.01(0.01) -0.01(0.01) 0.02(0.02) 0.02(0.01) -0.17(0.02) $q_{n}=281$ Cover
0.94(0.01) 0.94(0.01) 0.94(0.01) 0.92(0.01) 0.80(0.01) $\sqrt{n}$Bias
0.02(0.02) 0.02(0.02) 0.02(0.02) 0.05(0.02) -0.18(0.02) $q_{n}=421$ Cover