Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
On the consistency of bootstrap for matching estimators
{5pt}
{5pt}
{5pt}
{5pt}
\hypersetup{colorlinks,breaklinks,urlcolor=blue,linkcolor=blue}
abstractIn a landmark paper, abadie2008failure showed that the naive bootstrap is inconsistent when applied to nearest neighbor matching estimators of the average treatment effect with a fixed number of matches. Since then, this finding has inspired numerous efforts to address the inconsistency issue, typically by employing alternative bootstrap methods. In contrast, this paper shows that the naive bootstrap is provably consistent for the original matching estimator, provided that the number of matches, $M$, diverges. The bootstrap inconsistency identified by abadie2008failure thus arises solely from the use of a fixed $M$.
{\bf Keywords:} matching estimators, naive bootstrap, diverging number of matches.
Introduction
Nearest neighbor (NN) matching methods, which impute missing potential outcomes by averaging outcomes from matches in the opposite treatment group, are among the most widely used tools for causal inference stuart2010matching, imbens2024causal. These methods are praised as “intuitive and easy to explain” imbens2024causal and, unlike competing approaches heckman1997matching, heckman1998matching, hirano2003efficient, chernozhukov2018double, hinge on an easy-to-interpret parameter—the number of matches—making them user-friendly to tune.
However, these advantages are offset by a significant limitation: Efron’s naive bootstrap efron1979bootstrap fails to capture the stochastic behavior of NN matching estimators when the number of matches, $M$, is {\it fixed}. Specifically, the naive bootstrap produces inconsistent estimates of both the asymptotic distribution and the variance of the estimators. As a result, NN matching estimators fall into an interestingly rare class of statistics that are asymptotically normal, $\sqrt{n}$ -consistent, but bootstrap-inconsistent. This phenomenon was rigorously established through a series of seminal papers by Abadie and Imbens abadie2006large, abadie2008failure, abadie2011bias and echoed in parallel work on Chatterjee’s rank correlation chatterjee2020new, lin2022limit, lin2023failure.
The bootstrap inconsistency of NN matching estimators has spurred extensive research on alternative resampling methods. Proposed solutions include the wild bootstrap, which resamples the linear component of the estimator otsu2017bootstrap, adusumilli2018bootstrap, bodory2024nonparametric, the block bootstrap, which resamples matched sets instead of individual observations abadie2022robust, the $m$-out-of-$n$ bootstrap walsh2023nearest, and the bag of little bootstraps kosko2024fast, among others. Some empirical evidence suggests that bootstrap-based methods improve coverage rates over asymptotic approximations; see, for example, bodory2020finite and references therein.
This paper contributes to this line of research by revisiting the original setting of abadie2008failure. Building on recent insights from lin2023estimation and he2024propensity, which connect NN matching to augmented inverse probability weighting estimators and highlight the importance of allowing $M\to\infty$ for improved efficiency, we show that the counterexample in abadie2008failure no longer holds when $M$ grows with the sample size. Specifically, we demonstrate that once $M$ is allowed to diverge, the naive bootstrap yields consistent distributional estimates for NN matching estimators. Thus, the bootstrap inconsistency identified in abadie2008failure arises entirely from the practical constraint of fixing $M$.
Setup and methods
Matching methods
We consider the standard potential outcome causal model for a binary treatment, assuming $n$ realizations, $\{(X_i, D_i, Y_i(0), Y_i(1))\}_{i=1}^n$, drawn from a quadruple
\[
(X, D, Y(0), Y(1)).
\]
Here $Y(0),Y(1) \in \bR$ are the potential outcomes under control and treatment, respectively neyman1923applications, rubin1974estimating, the treatment indicator $D \in \{0,1\}$ denotes whether the subject received the treatment, and $X \in \cX \subset \bR^d$ corresponds to pretreatment covariates. For each subject $i\in\llbracket n\rrbracket:=\{1,2,\ldots,n\}$, we observe only the triple
\[
O_i := (X_i, D_i, Y_i = Y_i(D_i)).
\]
Write $\mX=[X_1,\ldots,X_n]\in\bR^{d\times n}$, $\mD=(D_1,\ldots, D_n)^\top\in\bR^n$, and $\mY=(Y_1,\ldots,Y_n)^\top\in\bR^n$. The goal of this paper is to estimate the population average treatment effect (ATE),
align*[align* omitted — 51 chars of source]
using only the observed data $(\mX,\mD,\mY)$.
Abadie and Imbens abadie2006large, abadie2011bias introduced a (bias-corrected) matching approach to estimating $\tau$. Their method imputes the missing potential outcome $Y_i(1-D_i)$ for each subject $i$ by averaging the outcomes $Y_j(D_j)$'s from subjects with opposite treatment status, i.e., those whose covariates are closest to $X_i$ under the Euclidean distance $\|\cdot\|$ with $D_j=1-D_i$. The resulting estimator is then adjusted by subtracting a bias estimate.
In detail, for each subject $i\in\llbracket n\rrbracket$, let $\mathcal{J}_M(i)$ denote the index set of $M$-NNs of $X_i$ among subjects with $D_j=1-D_i$; formally,
\[
\mathcal{J}_M(i) := \Big\{ j\in\llbracket n\rrbracket: D_j=1-D_i, \sum_{\ell=1, D_\ell =1- D_i}^n \ind \Big( \lVert X_\ell - X_i \rVert < \lVert X_j - X_i \rVert\Big) < M\Big\}.
\]
The missing potential outcome $Y_i(1-D_i)$ is then imputed by
\[
\hat Y_i(1-D_i) := \frac{1}{M}\sum_{j\in\cJ_M(i)}Y_j,
\]
along with the convention that $\hat Y_i(D_i)=Y_i$. The {\it bias-uncorrected} $M$-NN matching estimator of $\tau$ is accordingly given by
\[
\hat\tau_M := \frac{1}{n}\sum_{i=1}^n\Big\{\hat Y_i(1)-\hat Y_i(0) \Big\},
\]
as studied in abadie2006large.
However, as noted in abadie2006large, unless $d=1$, the estimator $\hat\tau_M$ suffers from an asymptotically non-negligible bias of the form
\[
B_M := \frac{1}{n}\sum_{i = 1}^n\frac{2D_i - 1}{M}\sum_{j \in \cJ_M(i)}\Big\{\mu_{1 - D_i}(X_i) - \mu_{1- D_i}(X_j)\Big\},
\]
where $\mu_0(\cdot)$ and $\mu_1(\cdot)$ are the outcome conditional means
\[
\mu_0(x) := {\mathrm E}[Y \,|\, X=x, D=0]~~ {\rm and} ~~\mu_1(x) := {\mathrm E}[Y \,|\, X=x, D=1],
\]
respectively. To address this bias, abadie2011bias proposed a simple bias correction. Let $\hat \mu_0(\cdot)$ and $\hat \mu_1(\cdot)$ be estimators of $\mu_0(\cdot)$ and $\mu_1(\cdot)$, respectively. The bias correction term is defined as
\[
\hat B_M := \frac{1}{n}\sum_{i=1}^n\frac{2D_i-1}{M}\sum_{j\in\cJ_M(i)}\Big\{\hat\mu_{1-D_i}(X_i)-\hat\mu_{1-D_i}(X_j)\Big\}.
\]
The final bias-corrected $M$-NN matching estimator is
\[
\hat\tau_M^{\rm bc}:=\hat\tau_M-\hat B_M.
\]
The naive bootstrap
This paper focuses on Efron’s naive bootstrap for inferring $\hat\tau_M^{\rm bc}$. Specifically, we consider the bootstrap sample defined by abadie2008failure,
\[
\Big\{O_i^*=(X_i^*,D_i^*,Y_i^*)\Big\}_{i=1}^n,
\]
where observations are independently resampled {\it within} the treated and control groups separately. The bootstrap sample thus consists of $n_0$ and $n_1$ independent draws with replacement from $\{O_i\}_{i=1,D_i=0}^n$ and $\{O_i\}_{i=1,D_i=1}^n$, respectively.
Next, we implement the bootstrap procedure to infer the matching estimator. Following the same steps as in Section (ref), we first define the bootstrap matching set as
\[
\tilde{\mathcal{J}}_M^*(i) := \Big\{ j\in\llbracket n\rrbracket: D_j^*=1-D_i^*, \sum_{\ell=1, D_\ell^* =1- D_i^*}^n \ind \Big( \lVert X_\ell^* - X_i^* \rVert < \lVert X_j^* - X_i^* \rVert\Big) < M\Big\}.
\]
Since ties can occur in the bootstrap sample, $\tilde{\mathcal{J}}_M^*(i)$ may contain more than $M$ elements. Similar to lin2023failure, we break ties {\it arbitrarily} to obtain a subset
\[
\mathcal{J}_M^*(i)\subset \tilde{\mathcal{J}}_M^*(i),~~~\text{with cardinality }~\Big|\mathcal{J}_M^*(i)\Big|=M
\]
for each $i\in\llbracket n\rrbracket$. The corresponding imputations are then
align*[align* omitted — 370 chars of source]
The bias-uncorrected bootstrap matching estimator is accordingly given by
\[
\hat\tau^*_M:=\frac{1}{n}\sum_{i=1}^n\Big\{\hat Y_i^*(1)-\hat Y_i^*(0)\Big\},
\]
and the bias-corrected bootstrap matching estimator is
\[
\hat\tau_M^{*, {\rm bc}} := \hat\tau^*_M-\hat B_M^*,
\]
where
\[
\hat B_M^* := \frac{1}{n}\sum_{i=1}^n\frac{2D_i^*-1}{M}\sum_{j\in\cJ^*_M(i)}\Big\{\hat\mu_{1-D_i^*}(X_i^*)-\hat\mu_{1-D_i^*}(X_j^*)\Big\}.
\]
Bootstrap theory
This section gives the bootstrap consistency guarantee for $\hat\tau_M^{*, {\rm bc}}$ provided that $M$ grows to infinity at an appropriate rate with the sample size.
Setup
Let us first rigorously document the two sources of randomness: one from the original data generating scheme and the other from the bootstrap. Efron's bootstrap corresponds to some random weights $\mW_0$ and $\mW_1$ that account for the numbers of times each observation in $\{O_i\}_{i=1,D_i=1}^n$ and $\{O_i\}_{i=1,D_i=0}^n$ appears in the bootstrap sample. Specifically, we have
align*[align* omitted — 162 chars of source]
to be multinomially distributed. Let $\mW=(W_1,\ldots,W_n)^\top$ be the concatenated vector of $\mW_0$ and $\mW_1$ according to the order of $O_i$'s.
In the following, we use $(\Omega_O, \cA_O, {\mathrm P}_O )$ to denote the probability space from which the observed data $\{O_i\}_{i=1}^n$ arises and $( \Omega_W, \cA_W, {\mathrm P}_W )$ to denote the probability space for the bootstrap weights. Let the probability measure on the product measurable space be ${\mathrm P}_{OW}$,
so that
align*[align* omitted — 291 chars of source]
which we call the {\it bootstrap space}; throughout this paper, we make the assumption that the bootstrap weights are independent of the data, so that ${\mathrm P}_{OW} = {\mathrm P}_O^\infty \times {\mathrm P}_W$. Denote the corresponding expectation with respect to ${\mathrm P}_{OW}$ by ${\mathrm E}_{OW}$, similarly for the notations ${\mathrm E}_{W|O}$, ${\mathrm E}_{O}$, and ${\mathrm E}_{W}$. For simplicity, if no confusion is possible, we shorthand ${\mathrm P}_O^\infty, {\mathrm P}_W$ as ${\mathrm P}$, and ${\mathrm E}_{O}, {\mathrm E}_{W}$ as ${\mathrm E}$.
Lastly, given a real-valued random quantity $R_n$ defined over the bootstrap space, we say $R_n$ is of an order $o_{{\mathrm P}_W}(1)$ in ${\mathrm P}_O$-probability if, for any $\epsilon, \delta>0$
align*[align* omitted — 139 chars of source]
Theory
In the following, for any two sequences of positive real numbers $\{a_n\}$ and $\{b_n\}$, we write $a_n\lesssim b_n$, or equivalently $a_n=O(b_n)$, if there exists a universal constant $C>0$ such that $a_n/b_n\leq C$ for all sufficiently large $n$. We write $a_n=o(b_n)$ if $a_n/b_n\to 0$ as $n\to\infty$. We use $o_{\mathrm P}(\cdot)$ and $O_{\mathrm P}(\cdot)$ to represent the stochastic versions of $o(\cdot)$ and $O(\cdot)$, respectively.
Our theory is based on the following two sets of assumptions.
assumption\phantomsection
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
• For almost all $x \in \cX$, $D$ is independent of $(Y(0),Y(1))$ conditional on $X=x$, and there exists some constant $\eta > 0$ such that
\[
\eta < e(x):={\mathrm P}(D=1 \,|\, X=x) < 1-\eta.
\]
• $\{(X_i,D_i,Y_i(0),Y_i(1))\}_{i=1}^n$ are independent copies of $(X,D,Y(0),Y(1))$.
• $\sigma_\omega^2(x):={\mathrm E} [U^2_\omega \,|\, X=x] ={\mathrm E}[(Y(\omega)-\mu_{\omega}(X))^2\,|\, X=x]$ is uniformly bounded by a constant $C_{\sigma}>0$ for almost all $x \in \cX$ and $\omega \in \{0,1\}$.
• ${\mathrm E} [\mu^2_\omega]$ is bounded for $\omega \in \{0,1\}$.
\end{enumerate}
assumption\phantomsection
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
• For any $\omega \in \{0,1\}$, $\max_{t \in \Lambda_{\floor{d/2}{}+1}} \lVert \partial^t {\mu}_\omega \rVert_{\infty}$ is bounded, where $\|\cdot\|_{\infty}$ stands for the function supremum norm, $\Lambda_k$ denotes the set of all $d$-dimension vectors of nonnegative integers that sums to $k$, and $\floor{\cdot}{}$ stands for the floor function.
• For any $\omega \in \{0,1\}$, the outcome regression estimator $\hat{\mu}_\omega(\cdot)$ satisfies
\begin{align*}
\max_{t \in \Lambda_{\floor{d/2}+1}} \big\lVert \partial^t \hat{\mu}_\omega \big\rVert_{\infty} = O_{{\mathrm P}_O}(1)
\end{align*}
and
\begin{align*}
\max_{t \in \Lambda_{\ell}} \big\lVert \partial^t \hat{\mu}_\omega - \partial^t {\mu}_\omega \big\rVert_{\infty}= O_{{\mathrm P}_O} ( n^{-\gamma_{\ell}}), for all \ell \in \llbracket \floor{d/2}\rrbracket
\end{align*}
with some constants $\{\gamma_\ell \}_{\ell=1}^{\floor{d/2}{}}$ satisfying $\gamma_\ell > \frac{1}{2} - \frac{\ell}{d}$ for all $\ell \in \llbracket \floor{d/2}{}\rrbracket$.
\end{enumerate}
Assumption (ref) corresponds to Assumption 4.1 in lin2023estimation, and is used for proving the asymptotic normality of $\hat\tau_M$. Assumption (ref) is needed for justifying bias correction. It aligns with Assumptions 4.4–4.5 in lin2023estimation, except that Assumption 4.4(i)(ii) has been omitted. This relaxation is made possible by a new linearization result, which we will present in Theorem (ref) ahead. It eliminates the need for a Lindeberg-Feller-type argument used in the proof of Lemma C.1 in lin2023estimation.
With the above two assumptions, the following theorem is the main result in lin2023estimation.
theorem[lin2023estimation] Assume that Assumption (ref) holds true and either $({\mathrm P}_{X\,|\, D=0}, {\mathrm P}_{X\,|\, D=1})$ or $({\mathrm P}_{X\,|\, D=1}, {\mathrm P}_{X\,|\, D=0})$ satisfies Assumption (ref) ahead.
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
• If $M \to \infty$ and $M \log n/n \to 0$ as $n \to \infty$, we then have
\[
\sup_{t \in \bR} \Big\lvert {\mathrm P}\big\{\sqrt{n}(\hat\tau_M-B_M- \tau)\le t\big\} -{\mathrm P}\big(N(0, \sigma^2) \le t\big) \Big\rvert = o(1),
\]
where
\[
\sigma^2:={\mathrm E}\Big[(\mu_1(X)-\mu_0(X)-\tau)^2 \Big]+{\mathrm E}\Big[\frac{\sigma_1^2(X)}{e(X)}+\frac{\sigma_0^2(X)}{1-e(X)}\Big]
\]
is the semiparametric efficiency lower bound hahn1998role.
• If further assuming Assumption (ref) and $M/n^\gamma\to 0$ with
\[
\gamma := \Big\{ \min_{\ell \in \llbracket \lfloor d/2 \rfloor\rrbracket } \Big[1-\Big(\frac{1}{2} - \gamma_\ell \Big)\frac{d}{\ell}\Big]\Big\} \wedge \Big[1- \frac{d/2}{\lfloor d/2 \rfloor + 1}\Big],
\]
we then have
\[
\sup_{t \in \bR} \Big\lvert {\mathrm P}\big(\sqrt{n}(\hat\tau_M^{\rm bc}- \tau)\le t\big) -{\mathrm P}\big(N(0, \sigma^2) \le t\big) \Big\rvert = o(1).
\]
\end{enumerate}
Our first contribution to the literature is a new linear representation result for the matching estimator $\hat\tau_M$.
theorem[Linear representation for $\hat\tau_M$]
Assume Assumptions (ref) and either $({\mathrm P}_{X\,|\, D=0}, {\mathrm P}_{X\,|\, D=1})$ or $({\mathrm P}_{X\,|\, D=1}, {\mathrm P}_{X\,|\, D=0})$ satisfies Assumption (ref) ahead. If $M \to \infty$ and $M \log n/n \to 0$ as $n \to \infty$, we then have
\[
\sqrt{n} (\hat\tau_M-B_M-\tau) = \frac{1}{\sqrt{n}} \sum_{i=1}^n \big\{\chi_0(O_i)-\tau\big\} + o_{{\mathrm P}_O}(1),
\]
where
\begin{align*}
\chi_0(O_i) := \mu_1(X_i) - \mu_0(X_i) + \Big\{\frac{D_i}{e(X_i)} - \frac{1-D_i}{1-e(X_i)}\Big\}\big\{Y_i - \mu_{D_i}(X_i) \big\}.
\end{align*}
Building on the insights from lin2023estimation, Theorem (ref) offers a new perspective on matching estimators. Unlike the proof of Lemma C.1 in lin2023estimation, which establishes a central limit theorem for $\hat\tau_M$ using an implicit conditioning argument, Theorem (ref) provides a straightforward asymptotically linear representation for $\hat\tau_M$. This representation simplifies the subsequent analyses in both the original and bootstrap spaces. Moreover, it connects to the linear expansion initially discovered by abadie2006large and later employed by otsu2017bootstrap. The differences between our form and theirs, on the other hand, are self-explanatory.
With Theorem (ref), we are now ready to introduce the main theorem of this paper.
theorem[Bootstrap distribution consistency]
Assume that Assumption (ref) holds true and either $({\mathrm P}_{X\,|\, D=0}, {\mathrm P}_{X\,|\, D=1})$ or $({\mathrm P}_{X\,|\, D=1}, {\mathrm P}_{X\,|\, D=0})$ satisfies Assumption (ref) ahead. If $M \to \infty$, $M \log n/n \to 0$, and $\log(n/M)/M \to 0$ as $n \to \infty$, we then have
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
• {\rm the linearization of $\hat\tau_M^*$:}
\[
\sqrt{n} (\hat \tau_M^* - B_M^*-\tau) = \frac{1}{\sqrt{n}} \sum_{i=1}^n W_i \big\{\chi_0(O_i)-\tau\big\} + o_{{\mathrm P}_{W}}(1) \text{ in }{\mathrm P}_O\text{-probability},
\]
where
\[
B_M^* := \frac{1}{n}\sum_{i = 1}^n\frac{2D_i^* - 1}{M}\sum_{j \in \cJ_M^*(i)}\Big\{\mu_{1 - D_i^*}(X_i^*) - \mu_{1- D_i^*}(X_j^*)\Big\};
\]
• further assuming Assumption (ref) and $M/n^\gamma\to 0$,
\[
\sqrt{n} (\hat \tau_M^{*, {\rm bc}}-\hat\tau_M^{\rm bc}) = \frac{1}{\sqrt{n}} \sum_{i=1}^n (W_i-1) \chi_0(O_i) + o_{{\mathrm P}_{W}}(1) \text{ in }{\mathrm P}_O\text{-probability,} \addtocounter{equation}{1}\tag{\theequation} \label{eq:matching,BTGn}
\]
so that
\[
\sup_{t \in \bR} \Big\lvert {\mathrm P}_{\mW|\mO}\big(\sqrt{n} (\hat\tau_M^{*,{\rm bc}} - \hat\tau_M^{\rm bc}) \le t \big) - {\mathrm P}\big(N(0, \sigma^2) \le t\big) \Big\rvert = o_{{\mathrm P}_O}(1)
\]
and
\[
\sup_{t \in \bR} \Big\lvert {\mathrm P}_{\mW|\mO}\big(\sqrt{n} (\hat\tau_M^{*, {\rm bc}} - \hat\tau_M^{\rm bc}) \le t \big) - {\mathrm P}\big(\sqrt{n}(\hat\tau_M^{\rm bc}- \tau) \le t\big) \Big\rvert = o_{{\mathrm P}_O}(1).
\]
\end{enumerate}
A key component of the proof of Theorem (ref) is demonstrating that the NN-matching-based density ratio estimator, discussed in lin2023estimation, remains a consistent estimator of the density ratio in the bootstrap space. This result, which will be presented in Section (ref) ahead, represents one of the main technical contributions of this paper. Of note, the bootstrap validity is obtained under nearly the same conditions as Theorem (ref), which slightly weakens the requirement of lin2023estimation.
In the special case of $d=1$, abadie2006large noted that there is no need of bias correction and the above theorem could be further simplified. The following corollary shows that their note is also correct when $M$ is allowed to vary with the sample size.
corollary[Bootstrap distribution consistency, the special case of $d=1$] Assume that Assumption (ref) holds true and either $({\mathrm P}_{X\,|\, D=0}, {\mathrm P}_{X\,|\, D=1})$ or $({\mathrm P}_{X\,|\, D=1}, {\mathrm P}_{X\,|\, D=0})$ satisfies Assumption (ref) ahead. Further assume that $\mu_0(x)$ and $\mu_1(x)$ are Lipschitz on $\cX$ and $d=1$. Then, if $M \to \infty$, $M^2/n \to 0$, and $\log(n/M)/M \to 0$ as $n \to \infty$, we have
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
• $\sqrt{n} (\hat\tau_M-\tau) = \frac{1}{\sqrt{n}} \sum_{i=1}^n (\chi_0(O_i)-\tau) + o_{{\mathrm P}_O}(1)$;
• $\sqrt{n} (\hat \tau_M^* -\tau) = \frac{1}{\sqrt{n}} \sum_{i=1}^n W_i(\chi_0(O_i)-\tau) + o_{{\mathrm P}_{W}}(1)$ in ${\mathrm P}_O$-probability;
• $
\sqrt{n} (\hat \tau_M^{*}-\hat\tau_M) = \frac{1}{\sqrt{n}} \sum_{i=1}^n (W_i-1) \chi_0(O_i) + o_{{\mathrm P}_{W}}(1)$
in ${\mathrm P}_O$-probability, so that
\[
\sup_{t \in \bR} \Big\lvert {\mathrm P}_{\mW|\mO}\big(\sqrt{n} (\hat\tau_M^{*} - \hat\tau_M) \le t \big) - {\mathrm P}\big(N(0, \sigma^2) \le t\big) \Big\rvert = o_{{\mathrm P}_O}(1)
\]
and
\[
\sup_{t \in \bR} \Big\lvert {\mathrm P}_{\mW|\mO}\big(\sqrt{n} (\hat\tau_M^{*} - \hat\tau_M) \le t \big) - {\mathrm P}\big(\sqrt{n}(\hat\tau_M^{\rm bc}- \tau) \le t\big) \Big\rvert = o_{{\mathrm P}_O}(1).
\]
\end{enumerate}
Density ratio estimation in the bootstrap space
Like lin2023estimation, for giving the proof of bootstrap consistency, we have to first discuss a bootstrap counterpart of NN matching-based density ratio estimator that receives similar attention in lin2023estimation.
With a slight abuse of notation, consider $X,Z$ to be two general random vectors in $\bR^d$ that are defined on the same probability space with laws $\nu_0$ and $\nu_1$, respectively. Furthermore, assume that both $\nu_0$ and $\nu_1$ are absolutely continuous with respect to the Lebesgue measure $\lambda$ on $\bR^d$ that is equipped with the Euclidean norm $\lVert \cdot \rVert$. Denote the Lebesgue densities of $\nu_0$ and $\nu_1$ by $f_0$ and $f_1$, respectively. Lastly, assume that $\nu_1$ is absolutely continuous with respect to $\nu_0$ and denote the corresponding density ratio by
\[
r(x) := f_1(x)/f_0(x),~ \text{ with the convention that } 0/0=0.
\]
Suppose we have the original sample comprised of $N_0$ copies $\{X_i\}_{i=1}^{N_0}$ of $X$ and $N_1$ copies $\{Z_i\}_{i=1}^{N_1}$ of $Z$, and assume that they are mutually independent. Note that, by the setup, there is almost surely no ties in the observations.
definition[NN matching, Definition 2.1 in lin2023estimation]
For any $x, z \in \bR^d$ and $M \in \llbracket N_0\rrbracket$,
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
• let $\cX_{(M)}(\cdot): \bR^d \to \{X_i\}_{i=1}^{N_0}$
be the map
that returns the value of the input $z$'s $M$-th NN in the original sample $\{X_i\}_{i=1}^{N_0}$, i.e., the value of $x \in \{X_i\}_{i=1}^{N_0}$ such that
\begin{align*}
\sum_{i=1}^{N_0} \ind\big(\lVert X_i - z \rVert \le \lVert x - z \rVert\big) = M;
\end{align*}
• let $K_M(\cdot): \bR^d\to \llbracket N_1\rrbracket$ be the map that returns the number of matched times of $x$, i.e.,
\begin{align*}
K_M(x) = K_M\big(x;\{X_i\}_{i=1}^{N_0},\{Z_j\}_{j=1}^{N_1}\big) := \sum_{j=1}^{N_1} \ind\big(\lVert x - Z_j \rVert \le \lVert \cX_{(M)}(Z_j) - Z_j \rVert\big);
\end{align*}
• let $A_M(\cdot): \bR^d\to \cB(\bR^d)$ be the corresponding map from $\bR^d$ to the class of all Borel sets in $\bR^d$ so that
\begin{align*}
A_M(x) = A_M\big(x;\{X_i\}_{i=1}^{N_0}\big) := \big\{z \in \bR^d : \lVert x - z \rVert \le \lVert \cX_{(M)}(z) - z \rVert\big\}
\end{align*}
returns the catchment area of $x$ in the setting of Definition (ref)(ref);
\end{enumerate}
Now, suppose we resample $\{X_i\}_{i=1}^{N_0}$ with replacement $N_0$ times according to the naive bootstrap weights $\mW_X \sim {\rm Mult}(N_0;N_0^{-1}, \ldots, N_0^{-1})$. Here each $W_{X_i}$ refers to the number of times that $X_i$ is sampled. Similarly, we resample $\{Z_i\}_{i=1}^{N_1}$ with replacement $N_1$ times according to the naive bootstrap weights $\mW_Z \sim {\rm Mult}(N_1;N_1^{-1}, \ldots, N_1^{-1})$, and let $W_{Z_j}$ refer to the number of times that $Z_j$ is sampled.
The bootstrap sample then yields the following counterparts of Definition (ref).
definition[bootstrap NN matching]
For any $x, z \in \bR^d$ and $M \in \llbracket N_0\rrbracket$,
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
• let $\cX_{(M)}^*(\cdot): \bR^d \to \{X_i^*\}_{i=1}^{N_0}$
be the map
that returns the value of the input $z$'s $M$-th NN in the bootstrap sample $\{X_i^*\}_{i=1}^{N_0}$, i.e., the value of $x \in \{X_i^*\}_{i=1}^{N_0}$ such that
\begin{align*}
\sum_{i=1}^{N_0} \ind\big(\lVert X_i^* - z \rVert < \lVert x - z \rVert\big) < M and \sum_{i=1}^{N_0} \ind\big(\lVert X_i^* - z \rVert \le \lVert x - z \rVert\big) \ge M;
\end{align*}
• let $K_M^{*,+}(\cdot), K_M^{*,-}(\cdot): \bR^d\to \llbracket N_1\rrbracket$ be the two bootstrap matched times of $x$ that account for ties in different ways, i.e.,
\begin{align*}
K_M^{*,+}(x) &= K_M^{*,+}\big(x;\{X_i^*\}_{i=1}^{N_0},\{Z_j^*\}_{j=1}^{N_1}\big) := \sum_{j=1}^{N_1} \ind\big(\lVert x-Z_j^*\rVert \le \lVert \cX_{M}^*(Z_j^*)-Z_j^*\rVert\big), \\
K_M^{*,-}(x) &= K_M^{*,-}\big(x;\{X_i^*\}_{i=1}^{N_0},\{Z_j^*\}_{j=1}^{N_1}\big) := \sum_{j=1}^{N_1} \ind\big(\lVert x-Z_j^*\rVert < \lVert \cX_{M}^*(Z_j^*)-Z_j^*\rVert\big);
\end{align*}
• let $A_M^{*,+}(\cdot), A_M^{*,-}(\cdot): \bR^d\to \cB(\bR^d)$ be the corresponding maps from $\bR^d$ to the class of all Borel sets in $\bR^d$ so that
\begin{align*}
A_M^{*,+}(x) &= A_M^{*,+}\big(x;\{X_i^*\}_{i=1}^{N_0}\big) := \big\{z \in \bR^d : \lVert x - z \rVert \le \lVert \cX_{(M)}^*(z) - z \rVert\big\} \\
{\rm and} A_M^{*,-}(x) &= A_M^{*,-}\big(x;\{X_i^*\}_{i=1}^{N_0}\big) := \big\{z \in \bR^d : \lVert x - z \rVert < \lVert \cX_{(M)}^*(z) - z \rVert\big\}
\end{align*}
return the catchment areas of $x$.
\end{enumerate}
We then give the bootstrap density ratio estimators, $\hat{r}_M^{*,+}(x), \hat{r}_M^{*,-}(x)$, as
align*[align* omitted — 323 chars of source]
Note the following relation that is immediate from the above definitions:
align*[align* omitted — 299 chars of source]
For any fixed point $x\in\bR^d$, we first obtain the following lemma that characterizes the bootstrap catchment areas $A_M^{*,+}(x)$ and $A_M^{*,-}(x)$. This is analogous to lin2023estimation.
lemma[Asymptotic $L^p$ moments of bootstrap catchment areas's $\nu_1$-measure]
Assuming $M\log N_0/N_0 \to 0$ as $N_0 \to \infty$, then
\begin{align*}
\lim_{N_0\to\infty} \frac{N_0}{M} {\mathrm E}_{OW}\big[\nu_1\big(A_M^{*,+}(x)\big)\big] = \lim_{N_0\to\infty} \frac{N_0}{M} {\mathrm E}_{OW}\big[\nu_1\big(A_M^{*,-}(x)\big)\big] = r(x) \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
holds for $\nu_0$-almost all $x$. If further assuming that $M \to \infty$ and $\log (N_0/M)/M \to 0$ as $N_0 \to \infty$, then for any positive integer $p$,
\begin{align*}
\lim_{N_0\to\infty} \Big(\frac{N_0}{M}\Big)^p {\mathrm E}_{OW}\big[\nu_1^p\big(A_M^{*,+}(x)\big)\big] = \lim_{N_0\to\infty} \Big(\frac{N_0}{M}\Big)^p {\mathrm E}_{OW}\big[\nu_1^p\big(A_M^{*,-}(x)\big)\big] = \big[r(x)\big]^p \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
holds for $\nu_0$-almost all $x$.
Next, we give a result that is analogous to lin2023estimation.
theorem[Pointwise consistency in the bootstrap space] Assume $M\log N_0/N_0 \to 0$.
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
• (Asymptotic unbiasedness) For $\nu_0$-almost all $x$, we have
\begin{align*}
\lim_{N_0\to\infty} {\mathrm E}_{OW} \big[\hat{r}_M^{*,+}(x)\big] = \lim_{N_0\to\infty} {\mathrm E}_{OW} \big[\hat{r}_M^{*,-}(x)\big] = r(x). \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
• (Pointwise $L^p$ consistency) Let $p$ be any positive integer and assume $MN_1/N_0 \to \infty$, $M \to \infty$, and $\log (N_0/M)/M \to 0$ as $N_0 \to \infty$. Then for $\nu_0$-almost all $x$, we have
\begin{align*}
\lim_{N_0\to\infty} {\mathrm E}_{OW} \big[ \lvert \hat{r}_M^{*,+}(x) - r(x) \rvert^p \big] = \lim_{N_0\to\infty} {\mathrm E}_{OW} \big[ \lvert \hat{r}_M^{*,-}(x) - r(x) \rvert^p \big] = 0. \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
\end{enumerate}
In the following, let $S_0$ and $S_1$ be the supports of $\nu_0$ and $\nu_1$, respectively. In addition, for any set $S\subset\bR^d$, let ${\rm diam}(S)$ denote its diameter in $(\bR^d,\|\cdot\|)$, i.e.,
\[
{\rm diam}(S):=\sup_{x,y\in S}\|x-y\|.
\]
We then obtain an $L^p$ risk consistency result in the bootstrap space under the following assumption. It is identical to Assumption B.1 in lin2023estimation.
assumption\phantomsection
Assume the following hold.
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
• $\nu_0, \nu_1$ are two probability measures on $\bR^d$, both are absolutely continuous with respect to $\lambda$, and $\nu_1$ is absolutely continuous with respect to $\nu_0$.
• There exists a constant $R>0$ such that ${\rm diam}(S_0) \le R$.
• There exist two constants $f_L,f_U>0$ such that for any $x \in S_0$ and $z \in S_1$, $f_L \le f_0(x) \le f_U$ and $f_1(z) \le f_U$.
• There exists a constant $a \in (0,1)$ such that for any $\delta \in (0,{\rm diam}(S_0)]$ and $z \in S_1$,
\[
\lambda(B_{z,\delta} \cap S_0) \ge a \lambda(B_{z,\delta}),
\]
where $B_{z,\delta}$ represents the closed ball in $\bR^d$ with center at $z$ and radius $\delta$.
\end{enumerate}
The following is the main result in this section; it is analogous to lin2023estimation and is the key to the proof of Theorem (ref).
theorem[$L^p$ risk consistency in the bootstrap space]
Assume the pair of $\nu_0, \nu_1$ satisfies Assumption (ref). Let $p$ be any positive integer. Assume further that $M \to \infty$, $M\log N_0/N_0 \to 0$, $MN_1/N_0 \to \infty$, and $\log (N_0/M)/M \to 0$ as $N_0 \to \infty$. We then have
\[
\lim_{N_0\to\infty} {\mathrm E}_{OW} \Big[ \int_{\bR^d} \Big\lvert \hat{r}_M^{*,+}(x) - r(x) \Big\rvert^p f_0(x) {\mathrm d} x \Big] = \lim_{N_0\to\infty} {\mathrm E}_{OW} \Big[ \int_{\bR^d} \Big\lvert \hat{r}_M^{*,-}(x) - r(x) \Big\rvert^p f_0(x) {\mathrm d} x \Big] = 0.
\]
Proofs of main results
In the following, we shorthand $K_M(X_i)$ as $K_M(i)$, $K_M^{*,+}(X_i)$ as $K_M^{*,+}(i)$, and $K_M^{*,-}(X_i)$ as $K_M^{*,-}(i)$.
Proof of Theorem (ref)
proof[Proof of Theorem (ref)]
Same as the decomposition in abadie2011bias, we write
\begin{align*}
\hat{\tau}_{M}
&= \frac{1}{n} \sum_{i=1}^n \Big[\mu_1(X_i) - \mu_0(X_i)\Big] + \frac{1}{n} \sum_{i=1}^n (2D_i-1) \Big(1 + \frac{K_M(i)}{M}\Big) \big(Y_i - \mu_{D_i}(X_i) \big) \\
& + \frac{1}{n}\sum_{i=1}^n (2D_i-1) \Big[\frac{1}{M}\sum_{j \in \mathcal{J}_M(i)} \Big\{ \mu_{1-D_i}(X_i)-\mu_{1-D_i}(X_{j})\Big\} \Big]\\
& =: \bar{\tau}(\mX) + E_M + B_M.
\end{align*}
It remains to show
\begin{align*}
R_M & := \sqrt{n}\Big(\bar{\tau}(\mX) + E_M - \frac{1}{n} \sum_{i=1}^n \chi_0(O_i) \Big) \\
& =\frac{1}{\sqrt{n}} \sum_{i=1}^n \Big[ (2D_i-1)\Big(1+\frac{K_M(i)}{M}\Big) - \Big(\frac{D_i}{e(X_i)} - \frac{1-D_i}{1-e(X_i)} \Big)\Big] \big(Y_i - \mu_{D_i}(X_i)\big)\\
&= o_{{\mathrm P}_O}(1).
\end{align*}
We only show that
\begin{align*}
R_{M,1}:=&\frac{1}{\sqrt{n}} \sum_{D_i=1} \Big[ (2D_i-1)\Big(1+\frac{K_M(i)}{M}\Big) - \Big(\frac{D_i}{e(X_i)} - \frac{1-D_i}{1-e(X_i)} \Big)\Big] \big(Y_i - \mu_{D_i}(X_i)\big)\\
=&\frac{1}{\sqrt{n}} \sum_{D_i=1} \Big[ 1+\frac{K_M(i)}{M} - \frac{1}{e(X_i)} \Big] \big(Y_i - \mu_{1}(X_i)\big)\\
=&o_{{\mathrm P}_O}(1),
\end{align*}
and those units with $D_i=0$ can be treated in a similar way.
Noticing that ${\mathrm E}[R_{M,1} \,|\, \mX, \mD ]=0$, one has
\begin{align*}
{\mathrm E}\big[R_{M,1}^2\big] &= \Var\big[R_{M,1}\big] = {\mathrm E}\big[\Var[R_{M,1} \,|\, \mX, \mD] \big] \\
&= {\mathrm E}\Big[\frac{1}{n} \sum_{i=1, D_i=1}^n \Big(1+\frac{K_M(i)}{M} - \frac{1}{e(X_i)} \Big)^2 \sigma_{1}^2(X_i) \Big] \le C_{\sigma} {\mathrm E}\Big[ 1+\frac{K_M(i)}{M} - \frac{1}{e(X_i)} \Big]^2 = o(1)
\end{align*}
since all cross-terms have zero conditional means and the last equality here is obtained by applying Theorem B.2 in lin2023estimation. This completes the proof.
Proof of Theorem (ref)
To prove Theorem (ref), we first present a lemma that gives an alternative representation of the bootstrap estimator defined in Section (ref). This lemma corresponds to Equation (3.1) and Lemma 3.1 in lin2023estimation.
lemmaThe bias-uncorrected bootstrap matching estimator $\hat{\tau}_{M}^{*}$ can be written as
\begin{align*}
\hat{\tau}_{M}^{*} := \frac{1}{n} \sum_{i=1}^n W_i (2D_i-1) \Big(1 + \frac{K^*_M(i)}{M}\Big) Y_i.
\end{align*}
Similarly, the bias-corrected bootstrap matching estimator $\hat{\tau}_{M}^{*,{\rm bc}}$ can be written as
\begin{align*}
\hat{\tau}_{M}^{*,{\rm bc}} = \frac{1}{n} \sum_{i=1}^n W_i \big[\hat{\mu}_1(X_i) - \hat{\mu}_0(X_i)\big] + \frac{1}{n} \sum_{i=1}^n W_i (2D_i-1) \Big(1 + \frac{K^*_M(i)}{M}\Big) \big( Y_i - \hat{\mu}_{D_i}(X_i) \big).
\end{align*}
Here, we write $W_i$ as the bootstrap weight for $O_i$ and $K_M^*(i)$ is some nonnegative integer that lies between $K_M^{*,-}(i)$ and $K_M^{*,+}(i)$ for all $i \in \llbracket n\rrbracket$.
In addition, conditioning on $\mD$, if the assumptions in Theorem (ref) are satisfied by either $({\mathrm P}_{X\,|\, D=0}, {\mathrm P}_{X\,|\, D=1})$ or $({\mathrm P}_{X\,|\, D=1}, {\mathrm P}_{X\,|\, D=0})$, then for any positive integer $p$, we have
align*[align* omitted — 532 chars of source]
They jointly imply
align*[align* omitted — 340 chars of source]
Before proving Theorem (ref), we give one more lemma for bounding the moments of
\[
U_{M,i}^* := \cX_{(M)}^*(X_i^*) - X_i^*.
\]
This corresponds to Lemma C.2 in lin2023estimation.
lemmaLet $p$ be any positive integer and suppose $M \to \infty$, and $\log (n/M)/M \to 0$ as $n \to \infty$. There then exists a constant $C_p>0$ only depending on $p$ such that, for any $i \in \llbracket n\rrbracket$,
\begin{align*}
\Big(\frac{n_{1-D_i}}{M}\Big)^{\frac{p}{d}} {\mathrm E}_{OW} \Big[\Big\| U_{M,i}^* \Big\|^p \,\Big{|}\, \mD\Big] \le C_p \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
holds true with probability $1$.
Now we are in place to prove the main theorem.
proof[Proof of Theorem (ref)]
We divide our proof into two parts.
{\bf Part I.}
In the first part, we prove Theorem (ref)(ref). To this end, we decompose $\hat{\tau}_M^*$ as
\begin{align*}
\hat{\tau}_{M}^{*}
&= \frac{1}{n} \sum_{i=1}^n W_i \Big[\mu_1(X_i) - \mu_0(X_i)\Big] + \frac{1}{n} \sum_{i=1}^n W_i (2D_i-1) \Big(1 + \frac{K^*_M(i)}{M}\Big) \big(Y_i - \mu_{D_i}(X_i) \big) \\
& + \frac{1}{n}\sum_{i=1}^n (2D_i^*-1) \Big[\frac{1}{M}\sum_{j \in \mathcal{J}_M^*(i)} \Big\{\mu_{1-D_i^*}(X_i^*)-\mu_{1-D_i^*}(X_{j}^*)\Big\}\Big]\\
& =: \bar{\tau}(\mX)^* + E_M^* + B_M^*.
\end{align*}
It then suffices to show
\begin{align*}
R_M^* & := \sqrt{n}\Big(\bar{\tau}(\mX)^* + E_M^* - \frac{1}{n} \sum_{i=1}^n W_i \chi_0(O_i) \Big) \\
& =\frac{1}{\sqrt{n}} \sum_{i=1}^n W_i \Big[ (2D_i-1)\Big(1+\frac{K_M^*(i)}{M}\Big) - \Big(\frac{D_i}{e(X_i)} - \frac{1-D_i}{1-e(X_i)} \Big)\Big] \big(Y_i - \mu_{D_i}(X_i)\big)\\
&= o_{{\mathrm P}_{OW}}(1). \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
Similar to the proof of Theorem (ref), we only have to show that the sum of treated units, $R_{M,1}^*$, is $o_{{\mathrm P}_{OW}}(1)$, and the sum of those units with $D_i=0$ can be analyzed similarly.
Observe that ${\mathrm E}_{OW}[R_{M,1}^* \,|\, \mW, \mX, \mD ]=0$, and the independence relations given by the bootstrap design are
\begin{align*}
Y_i & is independent of Y_j \,|\, (X_i,D_i), (X_j, D_j), \\
Y_i & is independent of \mW \,|\, (X_i,D_i), (X_j, D_j), Y_j, \\
Y_i & is independent of Y_j, \mW \,|\, (X_i,D_i), (X_j, D_j), \\
Y_i & is independent of Y_j \,|\, \mW, (X_i,D_i), (X_j, D_j),
\end{align*}
for any $i\neq j$ with $D_i=D_j=1$. This implies that
\begin{align*}
{\mathrm E}_{OW} \big[R_{M,1}^{*2}\big] &= \Var_{OW} \big[R_{M,1}^*\big] \\
&= {\mathrm E}_{OW} \big[\Var_{OW} [R_{M,1}^* \,|\, \mW, \mX, \mD] \big] \\
&= {\mathrm E}_{OW} \Big[\frac{1}{n} \sum_{i=1, D_i=1}^n W_{i}^2 \Big(1+\frac{K_M^*(i)}{M} - \frac{1}{e(X_i)} \Big)^2 \sigma_{1}^2(X_i) \Big] \\
&\le C_{\sigma} {\mathrm E}_{OW} \Big[ W_{i} \Big(1+\frac{K_M^*(i)}{M} - \frac{1}{e(X_i)}\Big) \Big]^2 \\
&= o(1)
\end{align*}
since all cross-terms of $R_{M,1}^{*2}$ have conditional zero means so that we can apply (ref) combined with Hölder's inequality, and use the fact that $W_{i}$'s have bounded moments. This completes the first part of our proof.
{\bf Part II.} We now prove Theorem (ref)(ref). This is separated to three more parts.
{\bf Part II.1.} We first show that
\[
\sqrt{n} (\hat \tau_M^{*, {\rm bc}}-\hat\tau_M^{\rm bc}) = \frac{1}{\sqrt{n}} \sum_{i=1}^n (W_i-1) \chi_0(O_i) + o_{{\mathrm P}_{OW}}(1).
\]
We have, under conditions on the accuracy of $\hat \mu_\omega$ and the smoothness condition on $\mu_\omega$ for $\omega \in \{0,1\}$,
\begin{align*}
\hat \tau_M^{*, {\rm bc}} = \bar{\tau}(\mX)^* + E_M^* + B_M^* - \hat B_M^*,
\end{align*}
and we can establish
\begin{align*}
\sqrt{n} \big(B_M^* - \hat{B}_M^* \big) = o_{{\mathrm P}_{OW}}(1).
\end{align*}
To see this, write
\begin{align*}
& \big\lvert B_M^* - \hat{B}_M^* \big\rvert \\
=& \Big\lvert \frac{1}{n}\sum_{i=1}^n (2D_i^*-1) \Big[\frac{1}{M}\sum_{j \in \mathcal{J}_M^*(i)} \Big\{ \mu_{1-D_i^*}(X_i^*)-\mu_{1-D_i^*}(X_{j}^*) - \hat{\mu}_{1-D_i^*}(X_i^*) + \hat{\mu}_{1-D_i^*}(X_{j}^*)\Big\} \Big] \Big\rvert\\
\le& \max_{i \in \llbracket n\rrbracket} \frac{1}{M}\sum_{j \in \mathcal{J}_M^*(i)} \Big\lvert \mu_{1-D_i^*}(X_i^*)-\mu_{1-D_i^*}(X_{j}^*) - \hat{\mu}_{1-D_i^*}(X_i^*) + \hat{\mu}_{1-D_i^*}(X_{j}^*) \Big\rvert\\
\le& \max_{i \in \llbracket n\rrbracket, j \in \mathcal{J}_M^*(i), \omega \in \{0,1\}} \Big\lvert \mu_\omega(X_i^*)-\mu_\omega(X_{j}^*) - \hat{\mu}_\omega(X_i^*) + \hat{\mu}_\omega(X_{j}^*) \Big\rvert. \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
Let $k$ equal to $\floor{d/2}{}+1$ and denote $X_j^* - X_i^*$ by $U_{j,i}^*$ for all $i,j \in \llbracket n\rrbracket$. Employing Taylor expansion of $\mu_\omega(\cdot)$ around $X_i^*$ to the $k$-th order for $\omega \in \{0,1\}$, one obtains
\begin{align*}
\Big\lvert \mu_\omega (X_{j}^*) - \mu_{\omega} (X_i^*) - \sum_{\ell = 1}^{k-1} \frac{1}{\ell!} \sum_{t \in \Lambda_\ell} \partial^t \mu_{\omega} (X_i^*) U^{*t}_{j,i} \Big\rvert \le \max_{t \in \Lambda_k} \lVert \partial^t \mu_{\omega} \rVert_\infty \frac{1}{k!} \sum_{t \in \Lambda_k} \lVert U_{j,i}^* \rVert^k, \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
and similarly
\begin{align*}
\Big\lvert \hat{\mu}_\omega (X_{j}^*) - \hat{\mu}_{\omega} (X_i^*) - \sum_{\ell = 1}^{k-1} \frac{1}{\ell!} \sum_{t \in \Lambda_\ell} \partial^t \hat{\mu}_{\omega} (X_i^*) U^{*t}_{j,i} \Big\rvert \le \max_{t \in \Lambda_k} \lVert \partial^t \hat{\mu}_{\omega} \rVert_\infty \frac{1}{k!} \sum_{t \in \Lambda_k} \lVert U_{j,i}^* \rVert^k, \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
with the derivative terms satisfying
\begin{align*}
\Big\lvert \sum_{\ell = 1}^{k-1} \frac{1}{\ell!} \sum_{t \in \Lambda_\ell} (\partial^t \hat{\mu}_{\omega} (X_i^*) - \partial^t \mu_{\omega} (X_i^*)) U^{*t}_{j,i} \Big\rvert \le \sum_{\ell = 1}^{k-1} \max_{t \in \Lambda_\ell} \lVert \partial^t \hat{\mu}_{\omega} - \partial^t \mu_{\omega} \rVert_\infty \frac{1}{\ell!} \sum_{t \in \Lambda_\ell} \lVert U_{j,i}^* \rVert^\ell. \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
Combining (ref)-(ref), from Assumption (ref) one obtains
\begin{align*}
\big\lvert B_M^* - \hat{B}_M^* \big\rvert & \le \Big( \max_{\omega \in \{0,1\}} \max_{ t \in \Lambda_k} \lVert \partial^t \mu_{\omega} \rVert_\infty + \max_{\omega \in \{0,1\}} \max_{t \in \Lambda_k} \lVert \partial^t \hat{\mu}_{\omega} \rVert_\infty \Big) \max_{i \in \llbracket n\rrbracket} \lVert U_{M,i}^* \rVert^k \\
& + \sum_{\ell = 1}^{k-1} \Big( \max_{\omega \in \{0,1\}} \max_{t \in \Lambda_\ell} \lVert \partial^t \hat{\mu}_{\omega} - \partial^t \mu_{\omega} \rVert_\infty \Big) \max_{i \in \llbracket n\rrbracket} \lVert U_{M,i}^* \rVert^\ell \\
& \lesssim \big(1+O_{{\mathrm P}_{O}}(1) \big) \max_{i \in \llbracket n\rrbracket} \lVert U_{M,i}^* \rVert^k + \max_{\ell \in \llbracket k-1\rrbracket} O_{{\mathrm P}_{O}}(n^{-\gamma_\ell}) \max_{i \in \llbracket n\rrbracket} \lVert U_{M,i}^* \rVert^\ell.
\end{align*}
Using Lemma (ref) and Markov inequality, for any $\ell \in \llbracket k\rrbracket$, one has
\[
\max_{i \in \llbracket n\rrbracket} \lVert U_{M,i}^* \rVert^\ell = O_{{\mathrm P}_{OW}} ( (M/n)^{\ell/d}).
\]
Then, by the order of $M$, we obtain
\begin{align*}
\big\lvert B_M^* - \hat{B}_M^* \big\rvert &\lesssim O_{{\mathrm P}_{OW}} \Big(\Big(\frac{M}{n}\Big)^{\frac{k}{d}}\Big) + \max_{\ell \in \llbracket k-1\rrbracket} O_{{\mathrm P}_{OW}}\Big( n^{-\gamma_\ell} \Big(\frac{M}{n}\Big)^{\frac{\ell}{d}}\Big) \\
&= o_{{\mathrm P}_{OW}}\big( n^{(\gamma-1)\frac{k}{d}}\big) + \max_{\ell \in \llbracket k-1\rrbracket} o_{{\mathrm P}_{OW}}\big( n^{(\gamma-1)\frac{\ell}{d}-\gamma_{\ell}}\big) \\
&= o_{{\mathrm P}_{OW}}\big( n^{-\frac{1}{2}} \big)
\end{align*}
as desired.
{\bf Part II.2.} We then move on to the weak convergence part.
Writing
\[
\hat \sigma^2 := \frac{1}{n}\sum_{i=1}^n (\chi_0(O_i) - \bar \chi_0)^2,~~\text{with }\bar\chi_0:=\frac1n\sum_{i=1}^n\chi_0(O_i),
\]
it is enough to show
\begin{align*}
\sup_{t \in \bR} \Big\lvert {\mathrm P}_{\mW|\mO}\Big( \frac{1}{\sqrt{n}} \sum_{i=1}^n (W_i-1) \chi_0(O_i) \le t \Big) - \Phi\Big(\frac{t}{\hat \sigma}\Big) \Big\rvert +
\sup_{t \in \bR} \Big\lvert \Phi\Big(\frac{t}{\hat \sigma}\Big) - \Phi\Big(\frac{t}{\sigma}\Big) \Big\rvert = o(1) \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
for almost all $\{O_i\}_{i=1}^n$. For the second term in (ref) to be $o(1)$ almost surely, it follows from $\hat \sigma$ converges to $\sigma$ almost surely (i.e., $\hat\sigma\stackrel{\sf a.s.}{\to}\sigma$) by the strong law of large numbers and continuous mapping theorem, along with the fact that $\Phi(\cdot)$, the cumulative distribution function for the standard normal distribution, is uniformly continuous.
For the first term in (ref), we leverage the Berry-Essen bound durrett2019probability conditioning on the original data to obtain
\begin{align*}
\sup_{t \in \bR} \Big\lvert {\mathrm P}_{\mW|\mO}\Big( \frac{1}{\sqrt{n}} \sum_{i=1}^n (W_i-1) \chi_0(O_i) \le t \Big) - \Phi\Big(\frac{t}{\hat \sigma}\Big) \Big\rvert \le \frac{4}{5\sqrt{n}} \frac{{\mathrm E}_{\mW |\mO}\big[\lvert \chi_0(O_i^*) - \bar \chi_0 \rvert^3 \big] }{ \Var_{\mW |\mO}[\chi_0(O_i^*)]^{3/2} }.
\end{align*}
Furthermore we can bound
\begin{align*}
\frac{4}{5\sqrt{n}} \frac{{\mathrm E}_{\mW |\mO}\big[\lvert \chi_0(O_i^*) - \bar \chi_0 \rvert^3 \big] }{ \Var_{\mW |\mO}[\chi_0(O_i^*)]^{3/2} } &= \frac{4}{5\sqrt{n}}
\frac{ n^{-1} \sum_{i=1}^n \lvert \chi_0(O_i) - \bar \chi_0 \rvert^3 }{ \hat \sigma^3 } \\
& \le \frac{16}{5n^{3/2} \hat \sigma^3 } \Big[ \sum_{i=1}^n \big\lvert \chi_0(O_i) - {\mathrm E}[\chi_0(O_i)] \big\rvert^3 + n \big\lvert {\mathrm E}[\chi_0(O_i)] - \bar \chi_0 \big\rvert^3 \Big].
\end{align*}
Again, by the strong law of large numbers, one has $\bar \chi_0 \stackrel{\sf a.s.}{\longrightarrow} {\mathrm E}[\chi_0(O_i)]$ and $\sigma>0$. The second term above then satisfies
\[
\big\lvert {\mathrm E}[\chi_0(O_i)] - \bar \chi_0 \big\rvert^3 / \sqrt{n} \hat \sigma^3 \stackrel{\sf a.s.}{\longrightarrow} 0.
\]
Moreover, Lemma (ref) implies the first term here also converges to zero almost surely, as we can take $p=2/3$ and note ${\mathrm E}[ \chi_0(O_i)^2] < \infty$. This establishes (ref).
{\bf Part II.3.} Lastly, leveraging Lemma (ref), we have both $o_{{\mathrm P}_O}(1)$ and $o_{{\mathrm P}_{OW}}(1)$ above are of order $o_{{\mathrm P}_W}(1)$ in ${\mathrm P}_O$-probability. Accordingly, Theorem (ref)(ref) follows from Theorem (ref) and Polya theorem lehmann1999elements.
Proof of Corollary (ref)
proof[Proof of Corollary (ref)]
It suffices to show that, when $d=1$,
\begin{align*}
n^{1/2} \lvert B_M \rvert = O_{{\mathrm P}_O}\big(M/n^{1/2}\big) and n^{1/2} \lvert B_M^* \rvert = O_{{\mathrm P}_{OW}}\big(M/n^{1/2}\big). \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
The first part in (ref) corresponds to Corollary 4.1 in lin2023estimation and we include a proof here for facilitating the bootstrap part argument.
For any $d$, using Hölder's inequality implies
\begin{align*}
\Big(\frac{n}{M}\Big)^{2/d} {\mathrm E}[ B_M ]^2 &\le \Big(\frac{n}{M}\Big)^{2/d} \max_{i \in \llbracket n\rrbracket, j \in \cJ_M(i), \omega \in \{0,1\}}{\mathrm E} \big[ \lvert \mu_{\omega}(X_i)-\mu_{\omega}(X_{j}) \rvert^2 \big] \\
& \le \Big(\frac{n}{M}\Big)^{2/d} C_1 \max_{i \in \llbracket n\rrbracket, \omega \in \{0,1\}}{\mathrm E} \big[ \lVert U_{M,i} \rVert^2 \big], \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
where $C_1$ is some constant given by the Lipschitz condition on $\mu_0$ and $\mu_1$, and $U_{M,i}:= \cX_{(M)}(X_i)-X_i$. Since by Lemma C.2 in lin2023estimation we know that (ref) is uniformly bounded over $i$, it then holds true that
\[
B_M = O_{{\mathrm P}_O}((M/n)^{1/d}).
\]
For the second part of (ref), we similarly write
\begin{align*}
\Big(\frac{n}{M}\Big)^{2/d} {\mathrm E}[ B_M^* ]^2 &\le \Big(\frac{n}{M}\Big)^{2/d} \max_{i \in \llbracket n\rrbracket, j \in \cJ_M^*(i), \omega \in \{0,1\}}{\mathrm E} \big[ \lvert \mu_{\omega}(X_i^*)-\mu_{\omega}(X_{j}^*) \rvert^2 \big] \\
& \le \Big(\frac{n}{M}\Big)^{2/d} C_1 \max_{i \in \llbracket n\rrbracket, \omega \in \{0,1\}}{\mathrm E} \big[ \lVert U_{M,i}^* \rVert^2 \big]. \addtocounter{equation}{1}\tag{\theequation}
\end{align*}
Leveraging Lemma (ref), we can also bound (ref) and then $B_M^* = O_{{\mathrm P}_{OW}}((M/n)^{1/d})$. This completes our proof in combination with Lemma (ref).