The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
178,740 characters
Estimation and Inference in Boundary Discontinuity Designs: Location-Based Methods Supplemental Appendix
\title{Estimation and Inference in Boundary Discontinuity Designs: Location-Based Methods\\ Supplemental Appendix\bigskip}
\author{Matias D. Cattaneo\thanks{Department of Operations Research and Financial Engineering, Princeton University.} \and
Rocio Titiunik\thanks{Department of Politics, Princeton University.} \and
Ruiqi (Rae) Yu\thanks{Department of Operations Research and Financial Engineering, Princeton University.}
}
\maketitle
\begin{abstract}
This supplemental appendix presents more general theoretical results encompassing those reported in the paper, their theoretical proofs, and other technical results. In particular, it presents a new strong approximation result for residual-based empirical processes leveraging and extending ideas from \cite{Cattaneo-Yu_2025_AOS}.
\end{abstract}
\textit{Keywords}: regression discontinuity, treatment effects estimation, causal inference.
\thispagestyle{empty}
\clearpage
\onehalfspacing
\setcounter{page}{1}
\pagestyle{plain}
\setcounter{tocdepth}{2}
\setcounter{secnumdepth}{4}
\tableofcontents
\clearpage
\section{Setup}\label{section: Setup}
This supplemental appendix considers a generalized version of the problem studied in the main paper: the location variable $\mathbf{X}_i$ is $d$-dimensional with $d\geq1$ and support $\mathcal{X}\subseteq\mathbb{R}^d$, and the boundary region $\mathcal{B}$ is a low dimensional manifold with ``effective dimension'' $d-1$. The special case considered in the paper is $d = 2$, that is, $\mathbf{X}_i$ is bivariate and $\mathcal{B}$ is a one-dimensional (boundary) curve.
Assumption 1 from the paper is generalized to the following.
\begin{assumption}[Data Generating Process]\label{sa-assump: DGP}
Let $t\in\{0,1\}$.
\begin{enumerate}[label=\normalfont(\roman*),noitemsep,leftmargin=*]
\item $(Y_1(t), \mathbf{X}_1^\top)^\top,\ldots, (Y_n(t), \mathbf{X}_n^\top)^\top$ are independent and identically distributed random vectors with $\mathcal{X} = \prod_{l = 1}^d [a_l, b_l]$ for $-\infty < a_l < b_l < \infty$ for $l = 1,\cdots,d$.
\item The distribution of $\mathbf{X}_i$ has a Lebesgue density $f_X(\mathbf{x})$ that is continuous and bounded away from zero on $\mathcal{X}$.
\item $\mu_t(\mathbf{x}) = \mathbb{E}[Y_i(t)| \mathbf{X}_i = \mathbf{x}]$ is $(p+1)$-times continuously differentiable on $\mathcal{X}$.
\item $\sigma^2_t(\mathbf{x}) = \mathbb{V}[Y_i(t)|\mathbf{X}_i = \mathbf{x}]$ is bounded away from zero and continuous on $\mathcal{X}$.
\item $\sup_{\mathbf{x} \in \mathcal{X}}\mathbb{E}[|Y_i(t)|^{2+v}|\mathbf{X}_i = \mathbf{x}] < \infty$ for some $v \geq 2$.
\end{enumerate}
\end{assumption}
We partition $\mathcal{X}$ into two areas, $\mathcal{A}_t\subset\mathbb{R}^d$ with $t\in\{0,1\}$, which represent the control and treatment regions, respectively. That is, $\mathcal{X} = \mathcal{A}_0 \cup \mathcal{A}_1$, where $\mathcal{A}_0$ and $\mathcal{A}_1$ are two disjoint regions in $\mathbb{R}^d$. The observed outcome is $Y_i = \mathds{1}(\mathbf{X}_i \in \mathcal{A}_0) Y_i(0) + \mathds{1}(\mathbf{X}_i \in \mathcal{A}_1) Y_i(1)$. The ``boundary'' now becomes $\mathcal{B} = \mathtt{bd}(\mathcal{A}_0) \cap \mathtt{bd}(\mathcal{A}_1)$ denotes the boundary determined by the assignment regions $\mathcal{A}_t$, $t\in\{0,1\}$, where $\mathtt{bd}(\mathcal{A}_t)$ denotes the topological boundary of $\mathcal{A}_t$. As in the paper, we assume that $\mathcal{B}$ belongs to $\mathcal{A}_1$, that is, $\mathcal{B} \subseteq \mathcal{A}_1$ and $\mathcal{B} \cup \mathcal{A}_0 = \emptyset$.
The multidimensional generalization of the three causal parameters studied in the paper are:
\begin{enumerate}
\item \emph{Boundary average treatment effect curve} (BATEC):
\begin{align*}
\tau(\mathbf{x}) = \mathbb{E}[Y_i(1) - Y_i(0)|\mathbf{X}_i = \mathbf{x}], \qquad \mathbf{x}\in\mathcal{B}\subseteq\mathbb{R}^{d-1}.
\end{align*}
\item \emph{Weighted Boundary average treatment effect} (WBATE):
\begin{align*}
\tau_{\mathtt{WBATE}} = \frac{\int_{\mathcal{B}} \tau(\mathbf{x}) w(\mathbf{x}) d \mathfrak{H}^{d-1}(\mathbf{x})}
{\int_{\mathcal{B}} w(\mathbf{x}) d \mathfrak{H}^{d-1}(\mathbf{x})}.
\end{align*}
\item \emph{Largest Boundary average treatment effect} (LBATE):
\begin{align*}
\tau(\mathbf{x}) = \sup_{\mathbf{x}\in\mathcal{B}} \tau(\mathbf{x}).
\end{align*}
\end{enumerate}
More generally, this supplemental appendix also considers the derivatives of the BATEC parameter:
\begin{align*}
\tau^{(\boldsymbol{\nu})}(\mathbf{x}) = \mu^{(\boldsymbol{\nu})}_1(\mathbf{x}) - \mu^{(\boldsymbol{\nu})}_0(\mathbf{x}), \qquad \mathbf{x} \in \mathcal{B},
\end{align*}
where, using standard multi-index notation, $\boldsymbol{\nu}=(\nu_1,\ldots,\nu_d)^\top\in\mathbb{N}_0$ with $|\boldsymbol{\nu}| = \nu_1+\ldots+\nu_d \leq p$ and $\mu^{(\boldsymbol{\nu})}_t(\mathbf{x}) = \partial^{\nu_1}\cdots\partial^{\nu_d} \mu_t(\mathbf{x})$ for $t\in\{0,1\}$.
The treatment effect estimator process along the boundary (submanifold) is
\begin{align*}
\Big(\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x})
= \widehat{\mu}^{(\boldsymbol{\nu})}_1(\mathbf{x}) - \widehat{\mu}^{(\boldsymbol{\nu})}_0(\mathbf{x}) : \mathbf{x}\in\mathcal{B} \Big),
\end{align*}
where $\widehat{\mu}_t^{(\boldsymbol{\nu})}(\mathbf{x}) = \mathbf{e}_{1 + \boldsymbol{\nu}}^\top \widehat{\boldsymbol{\beta}}_t(\mathbf{x})$ for $t\in\{0,1\}$ with
\begin{align*}
\widehat{\boldsymbol{\beta}}_t(\mathbf{x})
= \operatorname*{arg\,min}_{\boldsymbol{\beta} \in \mathbb{R}^{\mathfrak{p}_p+1}}
\mathbb{E}_n \Big[ \big(Y_i - \mathbf{r}_p(\mathbf{X}_i - \mathbf{x})^{\top}\boldsymbol{\beta} \big)^2 K_h(\mathbf{X}_i - \mathbf{x})\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) \Big], \qquad \mathbf{x} \in \mathcal{B},
\end{align*}
with $\mathfrak{p}_p = \frac{(d + p)!}{d! p !}$, $\mathbf{r}_p(\mathbf{u})$ denotes the $p$th order polynomial expansion of the $d$-variate vector $\mathbf{u}=(u_1,\cdots, u_d)^\top$, $K_h(\mathbf{u})=K(u_1/h,\cdots, u_d/h)/h^d$ for a $d$-variate kernel function $K(\cdot)$ and a bandwidth parameter $h$.
We impose the following assumption on the $d$-variate kernel function and assignment boundary (submanifold) $\mathcal{B}$.
\begin{assumption}[Kernel and Boundary]\label{sa-assump: Kernel and Boundary} Let $t \in \{0,1\}$.
\begin{enumerate}[label=\normalfont(\roman*)]
\item $\mathcal{B}$ is compact $(d-1)$-rectifiable, with $\mathfrak{H}^{d-1}(\mathcal{B})$ positive and finite.
\item $K: \mathbb{R}^d \to [0,\infty)$ is compact supported and Lipschitz continuous or $K(\mathbf{u}) = \mathds{1}(\mathbf{u} \in [-1,1]^d)$.
\item There exists a set $U \subseteq \mathbb{R}^d$, such that $K(\mathbf{u}) \geq \kappa > 0$ for all $\mathbf{u} \in U$, $\lambda_{\min} (\int_U \mathbf{r}_p(\mathbf{z}) \mathbf{r}_p(\mathbf{z})^{\top} d \mathbf{z}) > 0$, and $\liminf_{h \downarrow 0}\inf_{\mathbf{x} \in \mathcal{B}} \int_{U} K(\mathbf{u}) \mathds{1}(\mathbf{x} + h \mathbf{u} \in \mathcal{A}_t) d \mathbf{u} \gtrsim 1$.
\end{enumerate}
\end{assumption}
Note that in case $d = 2$, if we assume $\mathcal{B}$ is a rectifiable curve, then Assumption~\ref{sa-assump: Kernel and Boundary} (i) holds.
Under the assumptions imposed,
\begin{align*}
\widehat{\boldsymbol{\beta}}_{t}(\mathbf{x}) = \mathbf{H}^{-1} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} \mathbb{E}_n \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{x}) Y_i \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) \Big],
\end{align*}
where $\mathbf{H} = \operatorname*{diag}((h^{|\mathbf{v}|})_{0 \leq |\mathbf{v}| \leq p})$ with $\mathbf{v}$ running through all $\frac{d+\mathfrak{p}}{d!\mathfrak{p}!}$ multi-indices such that $|\mathbf{v}| \leq p$, and
\begin{align*}
\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}} = \mathbb{E}_n \Big[ \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)^{\top} K_h(\mathbf{X}_i - \mathbf{x}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big],
\end{align*}
where its population analogue is
\begin{align*}
\boldsymbol{\Gamma}_{t, \mathbf{x}} = \mathbb{E} \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)^{\top} K_h(\mathbf{X}_i - \mathbf{x}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) \Big].
\end{align*}
Note that $\left\lVert\mathbf{e}_{1 + \boldsymbol{\nu}}^\top \mathbf{H}^{-1}\right\rVert_2 = \left\lVert\mathbf{e}_{1 + \boldsymbol{\nu}}^\top \mathbf{H}^{-1}\right\rVert_{\infty} = h^{-|\boldsymbol{\nu}|}$. In addition, define
\begin{align*}
\mathbf{Q}_{t,\mathbf{x}} = \mathbb{E}_n \left[\mathbf{r}_p \left(\frac{\mathbf{X}_i - \mathbf{x}}{h}\right) K_h(\mathbf{X}_i - \mathbf{x}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) u_i\right],
\end{align*}
where $u_i = Y_i - [\mathds{1}(\mathbf{X}_i \in \mathcal{A}_0) \mu_0(\mathbf{X}_i) + \mathds{1}(\mathbf{X}_i \in \mathcal{A}_1) \mu_1(\mathbf{X}_i)] = Y_i - \mathbb{E}[Y_i|\mathbf{X}_i]$.
For $\mathbf{x}_1, \mathbf{x}_2 \in \mathcal{B}$ and $t \in \{0,1\}$, we introduce the following quantities:
\begin{align*}
\widehat{\boldsymbol{\Sigma}}_{t,\mathbf{x}_1,\mathbf{x}_2}
&= h^d \mathbb{E}_n \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}_1}{h}\Big) \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}_2}{h}\Big)^{\top} K_h(\mathbf{X}_i - \mathbf{x}_1) K_h(\mathbf{X}_i - \mathbf{x}_2) \widehat{\varepsilon}_i(\mathbf{x}_1) \widehat{\varepsilon}_i(\mathbf{x}_2) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big],\\
\boldsymbol{\Sigma}_{t, \mathbf{x}_1,\mathbf{x}_2}
&= h^d \mathbb{E} \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}_1}{h}\Big) \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}_2}{h}\Big)^{\top} K_h(\mathbf{X}_i - \mathbf{x}_1) K_h(\mathbf{X}_i - \mathbf{x}_2) \sigma_t^2(\mathbf{X}_i) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big],\\
\widehat{\Omega}_{t,\mathbf{x}_1,\mathbf{x}_2}^{(\boldsymbol{\nu})}
&= \frac{1}{n h^{d + 2 |\boldsymbol{\nu}|}} \mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}_1}^{-1} \widehat{\boldsymbol{\Sigma}}_{t,\mathbf{x}_1,\mathbf{x}_2} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}_2}^{-1} \mathbf{e}_{1 + \boldsymbol{\nu}},
\qquad\qquad \widehat{\Omega}_{\mathbf{x}_1,\mathbf{x}_2}^{(\boldsymbol{\nu})} = \widehat{\Omega}_{0,\mathbf{x}_1,\mathbf{x}_2}^{(\boldsymbol{\nu})} + \widehat{\Omega}_{1,\mathbf{x}_1, \mathbf{x}_2}^{(\boldsymbol{\nu})},\\
\Omega_{t,\mathbf{x}_1,\mathbf{x}_2}^{(\boldsymbol{\nu})}
&= \frac{1}{n h^{d + 2 |\boldsymbol{\nu}|}} \mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \boldsymbol{\Gamma}_{t,\mathbf{x}_1}^{-1} \boldsymbol{\Sigma}_{t,\mathbf{x}_1,\mathbf{x}_2} \boldsymbol{\Gamma}_{t,\mathbf{x}_2}^{-1} \mathbf{e}_{1 + \boldsymbol{\nu}},
\qquad\qquad \Omega_{\mathbf{x}_1,\mathbf{x}_2}^{(\boldsymbol{\nu})} = \Omega_{0,\mathbf{x}_1,\mathbf{x}_2}^{(\boldsymbol{\nu})} + \Omega_{1,\mathbf{x}_1,\mathbf{x}_2}^{(\boldsymbol{\nu})},
\end{align*}
where $\widehat{\varepsilon}_{i}(\mathbf{x}) = Y_i - \mathbf{r}_p(\mathbf{X}_i - \mathbf{x})^\top[\mathds{1}(\mathbf{X}_i \in \mathcal{A}_0) \widehat{\boldsymbol{\beta}}_0(\mathbf{x}) + \mathds{1}(\mathbf{X}_i \in \mathcal{A}_1) \widehat{\boldsymbol{\beta}}_1(\mathbf{x})]$.
Finally, to ensure that various weighted integral functionals on submanifolds (over the assignment boundary $\mathcal{B}$) are well-defined, we impose the following conditions on the weight function.
\begin{assumption}[Weight Function and Boundary]\label{sa-assump: Weight Function and Boundary}
Let $w:\mathcal{B} \mapsto \mathbb{R}$ with $\sup_{\mathbf{x}\in\mathcal{B}} |w(\mathbf{x})| < \infty$, $\inf_{\mathbf{x}\in\mathcal{B}} |w(\mathbf{x})| >0$, and $\int_{\mathcal{B}} |w(\mathbf{x})| d\mathfrak{H}^{d-1}(\mathbf{x}) < \infty$.
\end{assumption}
\subsection{Notation and Definitions}\label{sa-sec: notations and defns}
For textbook references on empirical process, see \cite{van-der-Vaart-Wellner_1996_Book}, \cite{dudley2014uniform}, and \cite{Gine-Nickl_2016_Book}. For textbook reference on geometric measure theory, see \cite{simon1984lectures}, \cite{federer2014geometric}, and \cite{folland2002advanced}.
\begin{enumerate}[label=(\roman*)]
\item \textit{Multi-index Notations}. For a multi-index $\mathbf{u} = (u_1, \ldots, u_d) \in \mathbb{N}^d$, denote $|\mathbf{u}| = \sum_{i = 1}^d u_d$, $\mathbf{u}! = \Pi_{i = 1}^d u_d$. Denote $\mathbf{r}_p(\mathbf{u}) = (1, u_1, \ldots, u_d, u_1^2, \ldots, u_d^2, \ldots, u_1^p, \ldots, u_d^p)$, that is, all monomials $u_1^{\alpha_1} \cdots u_d^{\alpha_d}$ such that $\alpha_i \in \mathbb{N}$ and $\sum_{i = 1}^d \alpha_i \leq p$. Define $\mathbf{e}_{1 + \boldsymbol{\nu}}$ to be the $p_d = \frac{(d + p)!}{d! p!}$-dimensional vector such that $\mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \mathbf{r}_p(\mathbf{u}) = \mathbf{u}^{\boldsymbol{\nu}}$ for all $\mathbf{u} \in \mathbb{R}^d$ and $|\boldsymbol{\nu}|\leq p$.
\item \textit{Norms}. For a vector $\mathbf{v} \in \mathbb{R}^k$, $\left\lVert\mathbf{v}\right\rVert = (\sum_{i = 1}^k \mathbf{v}_i^2)^{1/2}$, $\lVert \mathbf{v} \rVert_{\infty} = \max_{1 \leq i \leq k}|\mathbf{v}_i|$. For a matrix $A \in \mathbb{R}^{m \times n}$, $\left\lVertA\right\rVert_p = \sup_{\left\lVert\mathbf{x}\right\rVert_p = 1} \left\lVertA\mathbf{x}\right\rVert_p$, $p \in \mathbb{N} \cup \{\infty\}$, and $\lambda_{\min}(A)$ denotes its minimum eigenvalue. For a function $f$ on a metric space $(S, d)$, $\lVert f \rVert_{\infty} = \sup_{\mathbf{x} \in \mathcal{X}} |f(\mathbf{x})|$. For a probability measure $Q$ on $(\mathcal{S}, \mathscr{S})$ and $p \geq 1$, define $\left\lVertf\right\rVert_{Q,p} = (\int_{\mathcal{S}} |f|^p d Q)^{1/p}$. For a set $E \subseteq \mathbb{R}^d$, denote by $\mathfrak{m}(E)$ the Lebesgue measure of $E$.
\item \textit{Empirical Process}. We use standard empirical process notations: $\mathbb{E}_n[g(\mathbf{v}_i)] = \frac{1}{n} \sum_{i = 1}^n g(\mathbf{v}_i)$ and $\mathbb{G}_n[g(\mathbf{v}_i)] = \frac{1}{\sqrt{n}} \sum_{i = 1}^n (g(\mathbf{v}_i) - \mathbb{E}[g(\mathbf{v}_i)])$. Let $(\mathcal{S},d)$ be a semi-metric space. The covering number $N(\mathcal{S}, d, \varepsilon)$ is the minimal number of balls $B_s(\varepsilon) =\{t: d(t,s) < \varepsilon\}$ needed to cover $\mathcal{S}$. A \emph{$\mathbb{P}$-Brownian bridge} is a mean-zero Gaussian random function $W_n(f), f \in L_2(\mathcal{X}, \mathbb{P})$ with the covariance $\mathbb{E}[W_{\mathbb{P}}(f)W_{\mathbb{P}}(g)] = \mathbb{P}(fg) - \mathbb{P}(f)\mathbb{P}(g)$, for $f,g \in L_2(\mathcal{X},\mathbb{P})$. A class $\mathcal{F} \subseteq L_2(\mathcal{X}, \mathbb{P})$ is \emph{$\mathbb{P}$-pregaussian} if there is a version of $\mathbb{P}$-Brownian bridge $W_{\mathbb{P}}$ such that $W_{\mathbb{P}} \in C(\mathcal{F}; \rho_{\mathbb{P}})$ almost surely, where $\rho_{\mathbb{P}}$ is the semi-metric on $L_2(\mathcal{X},\mathbb{P})$ is defined by $\rho_{\mathbb{P}}(f, g) = (\|f - g\|_{\mathbb{P},2}^2 - (\int f \, d\mathbb{P} - \int g \, d\mathbb{P})^2)^{1/2}$, for $f, g \in L_2(\mathcal{X},\mathbb{P})$.
\item \textit{Geometric Measure Theory}. For a set $E \subseteq \mathcal{X}$, the \emph{De Giorgi perimeter of $E$ related to $\mathcal{X}$} is $\mathcal{L}(E) = \mathtt{TV}_{\{\mathds{1}_{E}\},\mathcal{X}}$. For $d \in \mathbb{N}$ and $0 \leq m \leq d$, the $m$-dimensional Hausdorff (outer) measure is given by $\mathfrak{H}^m(A) = \lim_{\delta \downarrow 0}\mathfrak{H}^m_{\delta}(A)$, $A \subseteq \mathbb{R}^d$, where for each $\delta > 0$, $\mathfrak{H}^m_{\delta}(A)$ is defined by taking $\mathfrak{H}^m_{\delta}(\emptyset) = 0$, and for any non-empty $A \subseteq \mathbb{R}^d$, $\mathfrak{H}^m_{\delta}(A) = \frac{\pi^{m/2}}{\Gamma(m/2+1)} \inf \sum_{j = 1}^{\infty} (\operatorname{diam}(C_j)/2)^m$, and the infimum is taken over all countable collections $C_1, C_2, \cdots$ of subsets of $\mathbb{R}^d$ such that $\operatorname{diam}(C_j) < \delta$ and $A \subseteq \cup_{j = 1}^{\infty}C_j$. Integration against $\mathfrak{H}^m$ is defined via Carathéodory's Theorem following the classical measure-theoretic literature. The Hausdorff dimension $\dim_{\mathfrak{H}}(A)$ of $A$ is defined by $\dim_{\mathfrak{H}}(A) = \inf\{t \geq 0: \mathfrak{H}^t(A) = 0\}$. A set $A \subseteq \mathbb{R}^d$ is said to be $k$-rectifiable if $A$ is of Hausdorff dimension $k$, and there exist a countable collection $\{f_i\}$ of continuously differentiable maps $f_i: \mathbb{R}^k \to \mathbb{R}^d$ such that $\mathfrak{H}^{k}(E \setminus \cup_{i = 0}^{\infty} f_i(\mathbb{R}^k)) = 0$. $B$ is a \emph{rectifiable curve} if there exists a Lipschitz continuous function $\gamma:[0,1] \to \mathbb{R}$ such that $B=\gamma([0,1])$. We define the curve length function of $B$ to be $\mathfrak{L}({B}) = \sup_{\pi \in \Pi} s(\pi, \gamma)$, where $\Pi = \left\{(t_0, t_1, \ldots, t_N): N \in \mathbb{N}, 0 \leq t_0 < t_1 < \ldots \leq t_N \leq 1\right\}$ and $s(\pi,\gamma) = \sum_{i = 0}^{N}\left\lVert\gamma(t_{i}) - \gamma(t_{i+1})\right\rVert_2$ for $\pi = (t_0, t_1, \ldots, t_N)$.
\item \textit{Bounds and Asymptotics}. For reals sequences $a_n = o(b_n)$ if $\limsup \frac{|a_n|}{|b_n|} = 0$, $a_n \lesssim b_n$ if there exists some constant $C$ and $N > 0$ such that $n > N$ implies $|a_n| \leq C |b_n|$. For sequences of random variables $a_n = o_{\mathbb{P}}(b_n)$ if $\operatorname{plim}_{n \to \infty}\frac{a_n}{b_n} = 0, |a_n| \lesssim_{\mathbb{P}} |b_n|$ if $\limsup_{M \to \infty} \limsup_{n \to \infty} \mathbb{P}[|\frac{a_n}{b_n}| \geq M] = 0$.
\end{enumerate}
All limits are taken such that $h\to0$ as $n\to\infty$. Most of our results hold with $h$ fixed but small enough, but we do not make this distinction explicit to avoid overly-complex statements.
\subsection{Mapping Between Paper and Supplement}
The results in the paper are special cases of the results in this supplemental appendix as follows.
\begin{itemize}
\item Theorem 1 in the paper corresponds to Theorem \ref{sa-thm: Convergence Rates} with $d=2$.
\item Theorem 2 in the paper corresponds to Theorem \ref{sa-thm: MSE} with $d=2$.
\item Theorem 3 in the paper corresponds to Theorems \ref{sa-thm: Confidence Intervals} and \ref{sa-thm: Confidence Bands} with $d=2$.
\item Theorem 4 in the paper corresponds to Theorem \ref{sa-thm: Integral: MSE Expansion} with $d=2$.
\item Theorem 5 in the paper corresponds to Theorem \ref{sa-thm: Integral: Asymptotic Normality} with $d=2$.
\item Theorem 6 in the paper corresponds to Theorem \ref{sa-thm: Confidence Bands for max} with $d=2$.
\end{itemize}
\section{Preliminary Lemmas}\label{sa-sec:technicals}
Let $\mathbf{X} = (\mathbf{X}_1^{\top}, \cdots, \mathbf{X}_n^{\top})$, and recall that $t \in \{0,1\}$.
\begin{lem}[Invertibility]\label{sa-lem:invert}
Suppose Assumption~\ref{sa-assump: DGP}(i,ii) and Assumption~\ref{sa-assump: Kernel and Boundary} hold. Then for $t = 0,1$,
$$\liminf_{h \to 0} \inf_{\mathbf{x} \in \mathcal{B}}\lambda_{\min}(\boldsymbol{\Gamma}_{t,\mathbf{x}}) > 0.$$
\end{lem}
\begin{lem}[Gram]\label{sa-lem: gram}
Suppose Assumption~\ref{sa-assump: DGP}(i,ii) and Assumption~\ref{sa-assump: Kernel and Boundary} hold. If $\frac{\log(1/h)}{n h^d} = o(1)$, then
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Gamma}}_{t, \mathbf{x}} - \boldsymbol{\Gamma}_{t, \mathbf{x}}\big\| \lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{n h^d}},
\qquad
\sup_{\mathbf{x} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Gamma}}_{t, \mathbf{x}}^{-1} - \boldsymbol{\Gamma}_{t, \mathbf{x}}^{-1}\big\| \lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{n h^d}},
\end{align*}
and if further $h = o(1)$, then $1 \lesssim_{\mathbb{P}} \inf_{\mathbf{x} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Gamma}}_{t, \mathbf{x}}\big\| \lesssim \sup_{\mathbf{x} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Gamma}}_{t, \mathbf{x}}\big\| \lesssim_{\mathbb{P}} 1$.
\end{lem}
\begin{lem}[Stochastic Linear Approximation]\label{sa-lem: Q}
Suppose Assumption~\ref{sa-assump: DGP}(i,ii,iv,v) and Assumption~\ref{sa-assump: Kernel and Boundary} hold. Suppose $\frac{\log(1/h)}{n h^d} = o(1)$, then
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big|\mathbf{Q}_{t,\mathbf{x}} \big| \lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{1+v}{2+v}}h^d},
\end{align*}
and if further $h = o(1)$,
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{\mu}_t^{(\boldsymbol{\nu})}(\mathbf{x}) - \mathbb{E}\big[\widehat{\mu}_t^{(\boldsymbol{\nu})}(\mathbf{x})\big|\mathbf{X}\big] - \mathbf{e}_{1+ \boldsymbol{\nu}}^{\top}\mathbf{H}^{-1}\boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1}\mathbf{Q}_{t,\mathbf{x}}\big| \lesssim_{\mathbb{P}} h^{-|\boldsymbol{\nu}|} \sqrt{\frac{\log(1/h)}{n h^d}} \bigg(\sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{1+v}{2+v}}h^d} \bigg).
\end{align*}
\end{lem}
\begin{lem}[Covariance]\label{sa-lem: covariance}
Suppose Assumptions \ref{sa-assump: DGP} and \ref{sa-assump: Kernel and Boundary} hold. If $\frac{\log(1/h)}{n h^d} = o(1)$, then
\begin{align*}
\sup_{\mathbf{x}_1,\mathbf{x}_2 \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Sigma}}_{t, \mathbf{x}_1,\mathbf{x}_2} - \boldsymbol{\Sigma}_{t, \mathbf{x}_1,\mathbf{x}_2}\big\|
\lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{v}{2+v}}h^d} + h^{p+1},
\end{align*}
\begin{align*}
\sup_{\mathbf{x}_1,\mathbf{x}_2 \in \mathcal{B}} \big|\widehat{\Omega}_{\mathbf{x}_1,\mathbf{x}_2}^{(\boldsymbol{\nu})} - \Omega_{\mathbf{x}_1,\mathbf{x}_2}^{(\boldsymbol{\nu})}\big|
\lesssim_{\mathbb{P}} (n h^{d + 2|\boldsymbol{\nu}|})^{-1}\bigg(\sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{v}{2+v}}h^d}+ h^{p+1}\bigg),
\end{align*}
and
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big|(\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2} - (\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2} \big|
\lesssim_{\mathbb{P}} \sqrt{n h^{d + 2 |\boldsymbol{\nu}|}}\bigg(\sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{v}{2+v}}h^d} + h^{p+1}\bigg).
\end{align*}
\end{lem}
\begin{lem}[Bias]\label{sa-lem: bias}
Suppose Assumption~\ref{sa-assump: DGP}(i,ii,iii) and Assumption~\ref{sa-assump: Kernel and Boundary} hold. If $\frac{\log(1/h)}{n h^d} = o(1)$ and $h = o(1)$, then
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big|\mathbb{E}[\widehat{\mu}_t^{(\boldsymbol{\nu})}(\mathbf{x})|\mathbf{X}] - \mu_t^{(\boldsymbol{\nu})}(\mathbf{x})\big| \lesssim_{\mathbb{P}} h^{p+1-|\boldsymbol{\nu}|},
\end{align*}
implying
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big|\mathbb{E}[\widehat{\mu}_t^{(\boldsymbol{\nu})}(\mathbf{x})|\mathbf{X}] - \mu_t^{(\boldsymbol{\nu})}(\mathbf{x}) - h^{p+1-|\boldsymbol{\nu}|} \widehat{B}_{t,\mathbf{x}}^{(\boldsymbol{\nu})}\big| = o_{\mathbb{P}}(h^{p+1-|\boldsymbol{\nu}|}),
\end{align*}
with $\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{B}_{t,\mathbf{x}}^{(\boldsymbol{\nu})} - B_{t,\mathbf{x}}^{(\boldsymbol{\nu})} \big| \lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{n h^d}}$, and hence $\sup_{\mathbf{x} \in \mathcal{B}}|\widehat{B}_{t,\mathbf{x}}^{(\boldsymbol{\nu})}| \lesssim_{\mathbb{P}} 1$.
\end{lem}
\section{Boundary Average Treatment Effect Curve}
\subsection{Point Estimation and MSE Expansions}
\begin{thm}[Convergence Rates]\label{sa-thm: Convergence Rates}
Suppose Assumptions \ref{sa-assump: DGP} and \ref{sa-assump: Kernel and Boundary} hold. If $\frac{\log(1/h)}{n h^d} = o(1)$ and $h = o(1)$, then
\begin{align*}
\big|\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) - \tau^{(\boldsymbol{\nu})}(\mathbf{x}) \big|
\lesssim_{\mathbb{P}} h^{-|\boldsymbol{\nu}|} \Big(\frac{1}{\sqrt{n h^d}} + \frac{1}{n^{\frac{1+v}{2+v}}h^d} + h^{p+1}\Big)
\end{align*}
for $\mathbf{x} \in \mathcal{B}$, and
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) - \tau^{(\boldsymbol{\nu})}(\mathbf{x})\big|
\lesssim_{\mathbb{P}} h^{-|\boldsymbol{\nu}|} \Big(\sqrt{\frac{\log(1/h)}{ n h^d}} + \frac{\log(1/h)}{n^{\frac{1+v}{2+v}}h^d} + h^{p+1}\Big).
\end{align*}
\end{thm}
The conditional mean-squared error (MSE) is
\begin{align*}
& \operatorname{MSE}_{\boldsymbol{\nu}}(\mathbf{x}) = \mathbb{E} \Big[(\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) - \tau^{(\boldsymbol{\nu})}(\mathbf{x}))^2 \Big| \mathbf{X} \Big]
\end{align*}
for $\mathbf{x}\in \mathcal{B}$, and the conditional integrated MSE (IMSE) is
\begin{align*}
& \operatorname{IMSE}_{\boldsymbol{\nu}} = \int_{\mathcal{B}} \operatorname{MSE}_{\boldsymbol{\nu}}(\mathbf{x}) w(\mathbf{x}) d \mathfrak{H}^{d-1}(\mathbf{x}),
\end{align*}
where $w(\mathbf{x})$ satisfies Assumption \ref{sa-assump: Weight Function and Boundary}. To state the MSE expansions, we introduce some more notation for the leading bias and variance:
\begin{align*}
B_{\mathbf{x}}^{(\boldsymbol{\nu})} = B_{1,\mathbf{x}}^{(\boldsymbol{\nu})} - B_{0,\mathbf{x}}^{(\boldsymbol{\nu})},
\qquad
B_{t,\mathbf{x}}^{(\boldsymbol{\nu})}
= \mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \sum_{|\boldsymbol{\omega}| = p + 1} \frac{\mu_t^{(\boldsymbol{\omega})}(\mathbf{x})}{\boldsymbol{\omega}!}
\mathbb{E} \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) \Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)^{\boldsymbol{\omega}} K_h(\mathbf{X}_i - \mathbf{x}) \Big],
\end{align*}
where
\begin{align*}
V_{\mathbf{x}}^{(\boldsymbol{\nu})} = V_{0,\mathbf{x}}^{(\boldsymbol{\nu})} + V_{1,\mathbf{x}}^{(\boldsymbol{\nu})},
\qquad
V_{t,\mathbf{x}}^{(\boldsymbol{\nu})} = \mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \boldsymbol{\Sigma}_{t,\mathbf{x},\mathbf{x}} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1}\mathbf{e}_{1 + \boldsymbol{\nu}} = nh^{d + 2|\boldsymbol{\nu}|} \Omega_{t,\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}.
\end{align*}
\begin{thm}[MSE Expansions]\label{sa-thm: MSE}
Suppose Assumptions \ref{sa-assump: DGP}, \ref{sa-assump: Kernel and Boundary} and \ref{sa-assump: Weight Function and Boundary} hold. If $\frac{\log(1/h)}{n h^d} = o(1)$ and $h = o(1)$, then
\begin{align*}
\operatorname{MSE}_{\boldsymbol{\nu}}(\mathbf{x})
= \big(h^{p+1-|\boldsymbol{\nu}|} B_{\mathbf{x}}^{(\boldsymbol{\nu})}\big)^2
+ \frac{1}{nh^{d + 2|\boldsymbol{\nu}|}} V_{\mathbf{x}}^{(\boldsymbol{\nu})}
+ o_{\mathbb{P}} \big(h^{2p+2-2|\boldsymbol{\nu}|} + n^{-1} h^{-d - 2|\boldsymbol{\nu}|}\big)
\end{align*}
for $\mathbf{x} \in \mathcal{B}$, and
\begin{align*}
\operatorname{IMSE}_{\boldsymbol{\nu}}
= \int_{\mathcal{B}} \Big[ \big( h^{p+1-|\boldsymbol{\nu}} B_{\mathbf{x}}^{(\boldsymbol{\nu})}\big)^2
+ \frac{1}{nh^{d + 2|\boldsymbol{\nu}|}} V_{\mathbf{x}}^{(\boldsymbol{\nu})} \Big] w(\mathbf{x}) d\mathfrak{H}^{d-1}(\mathbf{x})
+ o_{\mathbb{P}} \big(h^{2p+2-2|\boldsymbol{\nu}|} + n^{-1} h^{-d - 2|\boldsymbol{\nu}|}\big).
\end{align*}
\end{thm}
Theorem \ref{sa-thm: MSE} can be used to develop (feasible) bandwidth selectors. If $\widehat{B}_{\mathbf{x}}^{(\boldsymbol{\nu})} \neq 0$, the asymptotic MSE-optimal bandwidth is
\begin{align*}
h_{\operatorname{MSE},\boldsymbol{\nu},p}(\mathbf{x})
= \left(\frac{(d + 2|\boldsymbol{\nu}|) V_{\mathbf{x}}^{(\boldsymbol{\nu})}}
{(2p+2 - 2 |\boldsymbol{\nu}|) (B_{\mathbf{x}}^{(\boldsymbol{\nu})})^2} \frac{1}{n}\right)^{\frac{1}{2p+d+2}}
\end{align*}
for $\mathbf{x} \in \mathcal{B}$. Similarly, if $\int_{\mathcal{B}} (B_{\mathbf{x}}^{(\boldsymbol{\nu})})^2 w(\mathbf{x})d H^{d-1}(\mathbf{x}) \neq 0$, the asymptotic IMSE-optimal bandwidth is
\begin{align*}
h_{\operatorname{IMSE},\boldsymbol{\nu},p}
= \left(\frac{(d + 2 |\boldsymbol{\nu}|) \int_{\mathcal{B}} V_{\mathbf{x}}^{(\boldsymbol{\nu})} w(\mathbf{x})d\mathfrak{H}^{d-1}(\mathbf{x})}
{(2p+2 - 2|\boldsymbol{\nu}|) \int_{\mathcal{B}} (B_{\mathbf{x}}^{(\boldsymbol{\nu})})^2 w(\mathbf{x}) d\mathfrak{H}^{d-1}(\mathbf{x})} \frac{1}{n}\right)^{\frac{1}{2p+d+2}}.
\end{align*}
In practice, the the unknown bias and variance quantities can be replaced with (consistent) estimators thereof. For example, $\widehat{B}_{\mathbf{x}}^{(\boldsymbol{\nu})} = \widehat{B}_{1,\mathbf{x}}^{(\boldsymbol{\nu})} - \widehat{B}_{0,\mathbf{x}}^{(\boldsymbol{\nu})}$ with
\begin{align*}
\widehat{B}_{t,\mathbf{x}}^{(\boldsymbol{\nu})}
= \mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} \sum_{|\boldsymbol{\omega}| = p + 1} \frac{\mu_t^{(\boldsymbol{\omega})}(\mathbf{x})}{\boldsymbol{\omega}!}
\mathbb{E}_n \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) \Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)^{\boldsymbol{\omega}} K_h(\mathbf{X}_i - \mathbf{x}) \Big],
\qquad \\
\end{align*}
where the unknown functions $\mu_t^{(\boldsymbol{\omega})}(\mathbf{x})$ can be estimated using higher-order local polynomial estimators, and $\widehat{V}_{\mathbf{x}}^{(\boldsymbol{\nu})} = \widehat{V}_{0,\mathbf{x}}^{(\boldsymbol{\nu})} + \widehat{V}_{1,\mathbf{x}}^{(\boldsymbol{\nu})}$ with
\begin{align*}
\widehat{V}_{t,\mathbf{x}}^{(\boldsymbol{\nu})} = \mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} \widehat{\boldsymbol{\Sigma}}_{t,\mathbf{x},\mathbf{x}} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1}\mathbf{e}_{1 + \boldsymbol{\nu}},
\end{align*}
which corresponds to a standard variance estimator (which is also used for asymptotic inference as discussed below).
Finally, notice that the pointwise convergence rate and MSE expansion can be obtained under the slightly weaker side rate condition $n h^d\to\infty$. We do not make this distinction explicit to simplify the exposition.
\subsection{Distributional Approximation and Inference}\label{sa-sec: distributional approximation and inference}
Let $\mathbf{W} = ((\mathbf{X}_1^{\top},Y_1), \cdots, (\mathbf{X}_n^{\top},Y_n))$, and recall that $t \in \{0,1\}$. For $|\boldsymbol{\nu}| \leq p$, define the feasible $t$-statistic
\begin{align*}
\widehat{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})
= \frac{\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) - \tau^{(\boldsymbol{\nu})}(\mathbf{x})}{\sqrt{\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}}},
\qquad \mathbf{x} \in \mathcal{B}.
\end{align*}
The associated $100(1-\alpha)\%$ confidence interval estimator is
\begin{align*}
\widehat{\operatorname{I}}^{(\boldsymbol{\nu})}_{\alpha}(\mathbf{x})
= \bigg[\;\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) - \mathcal{q}_{\alpha} \sqrt{\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}}
\; , \;
\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) + \mathcal{q}_{\alpha} \sqrt{\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}}\;\bigg],
\end{align*}
where $\mathcal{q}_{\alpha}$ denotes an appropriate quantile depending on the desired confidence level $\alpha\in(0,1)$, and coverage objective (pointwise vs. uniform over $\mathcal{B}$). The following theorem establishes pointwise asymptotic normality and validity of confidence intervals. Let $\Phi(\cdot)$ be the cumulative distribution function of a standard univariate Gaussian random variable.
\begin{thm}[Confidence Intervals]\label{sa-thm: Confidence Intervals}
Suppose Assumptions \ref{sa-assump: DGP} and \ref{sa-assump: Kernel and Boundary} hold. If $n h^d \to \infty$ and $n h^d h^{2(p+1)} \to 0$, then
\begin{align*}
\sup_{u \in \mathbb{R}} \Big|\mathbb{P} \big(\widehat{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) \leq u \big) - \Phi(u)\Big| = o(1),
\qquad \mathbf{x} \in \mathcal{B},
\end{align*}
and
\begin{align*}
\mathbb{P} \big(\tau^{(\boldsymbol{\nu})}(\mathbf{x}) \in \widehat{\operatorname{I}}_{\alpha}^{(\boldsymbol{\nu})}(\mathbf{x}) \big) = 1 - \alpha + o(1),
\qquad \mathbf{x} \in \mathcal{B},
\end{align*}
provided that $\mathcal{q}_{\alpha} = \inf \{c > 0: \mathbb{P}( |\widehat{Z}| \geq c | \mathbf{W}) \leq \alpha \}$ with $\widehat{Z}|\mathbf{W} \thicksim \mathsf{Normal}(0, \widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})$.
\end{thm}
For uniform inference, we rely on a new strong approximation result established in Section \ref{sa-sec: Gaussian Strong Approximation}. First, we simplify the statistic $\widehat{\operatorname{T}}^{(\boldsymbol{\nu})}$, which is not directly a sum of independent random variables. Let
\begin{align*}
\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})
= \mathbb{E}_n \Big[(\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2} \mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \Big[\mathds{1}(\mathbf{X}_i \in \mathcal{A}_1) \boldsymbol{\Gamma}_{1, \mathbf{x}}^{-1} - \mathds{1}(\mathbf{X}_i \in \mathcal{A}_0) \boldsymbol{\Gamma}_{0, \mathbf{x}}^{-1}\Big] \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{x}) u_i \Big],
\end{align*}
where recall that $u_i = Y_i - \sum_{t \in \{0,1\}}\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) \mu_t(\mathbf{X}_i) = \mathbb{E}[Y_i | \mathbf{X}_i]$.
\begin{thm}[Stochastic Linearization]\label{sa-thm: stochastic linearization}
Suppose Assumptions \ref{sa-assump: DGP} and \ref{sa-assump: Kernel and Boundary} hold. If $\frac{\log(1/h)}{n^{\frac{v}{2+v}}h^d} = o(1)$ and $h = o(1)$, then
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \Big|\widehat{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) - \overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) \Big|
\lesssim_{\mathbb{P}} h^{p+1}\sqrt{n h^d} + \sqrt{\log(1/h)} \Big(\sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{v}{2+v}}h^d}\Big).
\end{align*}
\end{thm}
We can now exploit the linear structure of $(\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}): \mathbf{x} \in \mathcal{B})$, that is, an average of i.n.i.d. random vectors. Define the following functions indexed by $\mathbf{x} \in \mathcal{B}$:
\begin{align*}
g_{\mathbf{x}}(\mathbf{u}) = \mathds{1}(\mathbf{u} \in \mathcal{A}_1) \mathscr{K}^{(\boldsymbol{\nu})}_1(\mathbf{u};\mathbf{x}) - \mathds{1}(\mathbf{u} \in \mathcal{A}_0) \mathscr{K}^{(\boldsymbol{\nu})}_0(\mathbf{u};\mathbf{x}),
\qquad \mathbf{u} \in \mathcal{X},
\end{align*}
and
\begin{align*}
\mathscr{K}^{(\boldsymbol{\nu})}_t(\mathbf{u};\mathbf{x}) = n^{-1/2} (\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2}\mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \mathbf{r}_p \left(\frac{\mathbf{u} - \mathbf{x}}{h}\right) K_h(\mathbf{u} - \mathbf{x}), \qquad \mathbf{u} \in \mathcal{X},\quad t \in \{0,1\}.
\end{align*}
Define the associated class of functions $\mathcal{G} = \{g_{\mathbf{x}}: \mathbf{x} \in \mathcal{B}\}$ and $\mathcal{R} = \{\operatorname{Id}\}$, where $\operatorname{Id}(x) = x$, for all $x \in \mathbb{R}$. Then, the \emph{residual-based empirical process} is
\begin{align*}
R_n(g,r) = n^{-1/2}\sum_{i = 1}^n \Big[g(\mathbf{X}_i)r(Y_i) - g(\mathbf{X}_i) \mathbb{E}[r(Y_i)|\mathbf{X}_i] \Big], \qquad g \in \mathcal{G}, r \in \mathcal{R},
\end{align*}
and therefore
\begin{align*}
\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) = R_n(g_{\mathbf{x}},\operatorname{Id}), \qquad \mathbf{x} \in \mathcal{B}.
\end{align*}
Leveraging ideas in \cite{Cattaneo-Yu_2025_AOS}, Theorem~\ref{sa-lem: sa thm} gives a new Gaussian strong approximation that can be applied to our current setup. Specifically, our new theorem allows for polynomial moment bound on the conditional distribution of $Y_i|\mathbf{X}_i$.
\begin{thm}[Gaussian Strong Approximation: $\overline{\operatorname{T}}^{(\boldsymbol{\nu})}$]\label{sa-thm: Gaussian Strong Approximation: Tstat}
Suppose Assumptions \ref{sa-assump: DGP} and \ref{sa-assump: Kernel and Boundary} hold, and that there exists a constant $C > 0$ such that for $t \in \{0,1\}$ and for any $\mathbf{x} \in \mathcal{B}$, the De Giorgi perimeter of the set $E_{t,\mathbf{x}} = \{\mathbf{y} \in \mathcal{A}_t: (\mathbf{y} - \mathbf{x})/h \in \operatorname{Supp}(K)\}$ satisfies $\mathcal{L}(E_{t,\mathbf{x}}) \leq C h^{d-1}$. If $\liminf_{n \to \infty} \frac{\log h}{\log n} > - \infty$ and $n h^d \to \infty$ as $n \to \infty$, then (on a possibly enlarged probability space) there exists a mean-zero Gaussian process $Z^{(\boldsymbol{\nu})}$ indexed by $\mathcal{B}$ with almost surely continuous sample path such that
\begin{align*}
\mathbb{E}\Big[\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) - Z^{(\boldsymbol{\nu})}(\mathbf{x}) \big| \Big]
\lesssim (\log n)^{\frac{3}{2}} \bigg(\frac{1}{n h^d}\bigg)^{\frac{1}{2d+2}\cdot \frac{v}{v + 2}} + \log (n) \left(\frac{1}{n^{\frac{v}{2+v}}h^d}\right)^{\frac{1}{2}},
\end{align*}
where $\lesssim$ is up to a universal constant, and $Z^{(\boldsymbol{\nu})}$ has the same covariance structure as $\overline{\operatorname{T}}^{(\boldsymbol{\nu})}$; that is, $\mathbb{C}\mathrm{ov}[\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}_1), \overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}_2)] = \mathbb{C}\mathrm{ov}[Z^{(\boldsymbol{\nu})}(\mathbf{x}_1), Z^{(\boldsymbol{\nu})}(\mathbf{x}_2)]$ for all $\mathbf{x}_1, \mathbf{x}_2 \in \mathcal{B}$.
\end{thm}
Theorem \ref{sa-thm: Gaussian Strong Approximation: Tstat} can be used to construct confidence bands for $(\tau^{(\boldsymbol{\nu})}(\mathbf{x}):\mathbf{x}\in\mathcal{B})$. Let $(\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x}):\mathbf{x} \in \mathcal{B})$ be a (conditionally on $\mathbf{W}$) mean-zero Gaussian process with feasible (conditional) covariance function
\begin{align*}
\mathbb{C}\mathrm{ov} \Big[\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x}_1),\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x}_2) \Big|\mathbf{W} \Big]
= \frac{\widehat{\Omega}_{\mathbf{x}_1,\mathbf{x}_2}^{(\boldsymbol{\nu})}}
{\sqrt{\widehat{\Omega}_{\mathbf{x}_1,\mathbf{x}_1}^{(\boldsymbol{\nu})} \widehat{\Omega}_{\mathbf{x}_2,\mathbf{x}_2}^{(\boldsymbol{\nu})}}},
\qquad \mathbf{x}_1, \mathbf{x}_2 \in \mathcal{B}.
\end{align*}
\begin{thm}[Confidence Bands]\label{sa-thm: Confidence Bands}
Suppose the assumptions and conditions in Theorem \ref{sa-thm: Gaussian Strong Approximation: Tstat} hold. If $\liminf_{n \to \infty} \frac{\log h}{\log n} > - \infty$, $\frac{(\log n)^{3}}{n^{\frac{v}{2+v}}h^d} = o(1)$ and $h^{p+1}\sqrt{n h^d} = o(1)$, then
\begin{align*}
\sup_{u \in \mathbb{R}} \Big|\mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u \Big)
- \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u \Big| \mathbf{W} \Big) \Big| = o_{\mathbb{P}}(1)
\end{align*}
and
\begin{align*}
\mathbb{P}\Big[\tau^{(\boldsymbol{\nu})}(\mathbf{x}) \in \widehat{\operatorname{I}}_{\alpha}^{(\boldsymbol{\nu})}(\mathbf{x}), \text{ for all } \mathbf{x} \in \mathcal{B} \Big] = 1 - \alpha + o(1),
\end{align*}
provided that $\mathcal{q}_{\alpha} = \inf \big\{c > 0: \mathbb{P} \big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x})\big|\geq c \big| \mathbf{W} \big) \leq \alpha \big\}$.
\end{thm}
\section{Weighted Boundary Average Treatment Effect}
Without loss of generality, we set $\int_{\mathcal{B}} w(\mathbf{b}) d\mathfrak{H}^{d-1} (\mathbf{b}) = 1$, and the parameter of interest is the (weighted) average treatment effect along the boundary:
\begin{align*}
\tau_{\mathtt{WBATE}} = \int_{\mathcal{B}} \tau(\mathbf{b}) w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b}),
\end{align*}
where the weight function $w:\mathcal{X} \mapsto \mathbb{R}$ satisfies Assumption \ref{sa-assump: Weight Function and Boundary}.
The (weighted) boundary average treatment effect estimator along the boundary is
\begin{align*}
\widehat{\tau}_{\mathtt{WBATE}} = \int_\mathcal{B} \widehat{\tau}(\mathbf{b}) w(\mathbf{b}) \, d \mathfrak{H}^{d-1}(\mathbf{b}),
\end{align*}
Our first lemma in this section studies the conditional bias of $\tau_{\mathtt{WBATE}}$. Let
\begin{align*}
B_{\mathtt{WBATE}} = B_{1,\mathtt{WBATE}} - B_{0,\mathtt{WBATE}},
\qquad
B_{t,\mathtt{WBATE}} = \int_{\mathcal{B}} B^{(\mathbf{0})}_{t,\mathbf{b}} w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b}),
\end{align*}
for $t\in\{0,1\}$.
\begin{lem}[Bias: WBATE]\label{sa-lem: Integral: Bias}
Suppose Assumption~\ref{sa-assump: DGP}(i)-(iii), \ref{sa-assump: Kernel and Boundary} and \ref{sa-assump: Weight Function and Boundary} hold. If $\frac{\log(1/h)}{n h^d} = o(1)$ and $h = o(1)$, then
\begin{align*}
\mathbb{E}[\widehat{\tau}_{\mathtt{WBATE}}|\mathbf{X}] - \tau_{\mathtt{WBATE}} = h^{p+1} B_{\mathtt{WBATE}} + o_{\mathbb{P}}(h^{p+1}).
\end{align*}
\end{lem}
The next lemma studies the conditional variance of $\tau_{\mathtt{WBATE}}$, and a plug-in estimator thereof. Let
\begin{align*}
\Omega_{\mathtt{WBATE}} = \Omega_{1,\mathtt{WBATE}} + \Omega_{0,\mathtt{WBATE}},
\qquad
\Omega_{t,\mathtt{WBATE}} = \int_{\mathcal{B}} \int_{\mathcal{B}} \Omega^{(\mathbf{0})}_{t,\mathbf{b}_1,\mathbf{b}_2} w(\mathbf{b}_1) w(\mathbf{b}_2) d \mathfrak{H}^{d-1}(\mathbf{b}_1) d \mathfrak{H}^{d-1}(\mathbf{b}_2)
\end{align*}
and
\begin{align*}
\widehat{\Omega}_{\mathtt{WBATE}} = \widehat{\Omega}_{1,\mathtt{WBATE}} + \widehat{\Omega}_{0,\mathtt{WBATE}},
\qquad
\widehat{\Omega}_{t,\mathtt{WBATE}} = \int_{\mathcal{B}} \int_{\mathcal{B}} \widehat{\Omega}^{(\mathbf{0})}_{t,\mathbf{b}_1, \mathbf{b}_2} w(\mathbf{b}_1) w(\mathbf{b}_2) d \mathfrak{H}^{d-1}(\mathbf{b}_1) d \mathfrak{H}^{d-1}(\mathbf{b}_2),
\end{align*}
for $t\in\{0,1\}$.
\begin{lem}[Variance: WBATE]\label{sa-lem: Integral: Variance}
Suppose Assumptions \ref{sa-assump: DGP}, \ref{sa-assump: Kernel and Boundary} and \ref{sa-assump: Weight Function and Boundary} hold. If $\frac{\log(1/h)}{n h^d} = o(1)$ and $h = o(1)$, then
\begin{align*}
\mathbb{V}[\widehat{\tau}_{\mathtt{WBATE}}|\mathbf{X}]
= \Omega_{\mathtt{WBATE}} + O_{\mathbb{P}} \Big(h^{d-1}\frac{\log(1/h)^{1/2}}{(n h^d)^{3/2}}\Big)
= \Omega_{\mathtt{WBATE}} + o_{\mathbb{P}}((n h)^{-1}),
\end{align*}
where
\begin{align*}
(n h)^{-1} \lesssim \Omega_{\mathtt{WBATE}} \lesssim (n h)^{-1}.
\end{align*}
If, in addition, $\frac{\log(1/h)}{n^{\frac{v}{2+v}}h^d} = o(1)$, then
\begin{align*}
\mathbb{V}[\widehat{\tau}_{\mathtt{WBATE}}|\mathbf{X}] & = \widehat{\Omega}_{\mathtt{WBATE}} + o_{\mathbb{P}}((n h)^{-1}).
\end{align*}
\end{lem}
\begin{thm}[MSE Expansion: WBATE]\label{sa-thm: Integral: MSE Expansion}
Suppose Assumptions \ref{sa-assump: DGP}, \ref{sa-assump: Kernel and Boundary} and \ref{sa-assump: Weight Function and Boundary} hold. If $\frac{\log(1/h)}{n^{\frac{v}{2+v}}h^d} = o(1)$ and $h = o(1)$, then
\begin{align*}
\mathbb{E} [(\widehat{\tau}_{\mathtt{WBATE}} - \tau_{\mathtt{WBATE}})^2|\mathbf{X}]
= \Omega_{\mathtt{WBATE}} + h^{2p+2} B_{\mathtt{WBATE}}^2 + o_{\mathbb{P}}((n h)^{-1}) + o_{\mathbb{P}}(h^{2p+2}).
\end{align*}
\end{thm}
MSE-optimal bandwidth selection follows directly from Theorem \ref{sa-thm: Integral: MSE Expansion}.
For inference, we consider the feasible $t$-statistics
\begin{align*}
\widehat{\operatorname{T}}_{\mathtt{WBATE}} = \frac{\widehat{\tau}_{\mathtt{WBATE}} - \tau_{\mathtt{WBATE}}}{\sqrt{\widehat{\Omega}_{\mathtt{WBATE}}}}.
\end{align*}
\begin{thm}[Asymptotic Normality: WBATE]\label{sa-thm: Integral: Asymptotic Normality}
Suppose Assumptions \ref{sa-assump: DGP}, \ref{sa-assump: Kernel and Boundary} and \ref{sa-assump: Weight Function and Boundary} hold. If $\frac{\log(1/h)}{n^{\frac{v}{2+v}}h^d} = o(1)$ and $n h^{2p+3} = o(1)$, then
\begin{align*}
\sup_{u \in \mathbb{R}} \big|\mathbb{P} (\widehat{\operatorname{T}}_{\mathtt{WBATE}} \leq u ) - \Phi(u) \big| = o(1).
\end{align*}
\end{thm}
\section{Largest Boundary Average Treatment Effect}
Consider the maximum treatment effect over the boundary, defined by
\begin{align*}
\tau_{\mathtt{LBATE}} = \sup_{\mathbf{b} \in \mathcal{B}} \tau(\mathbf{b}).
\end{align*}
\begin{thm}[Convergence Rate: LBATE]\label{sa-thm: Convergence Rate for max}
Suppose Assumptions \ref{sa-assump: DGP} and \ref{sa-assump: Kernel and Boundary} hold. If $\frac{\log(1/h)}{n h^d} = o(1)$ and $h = o(1)$, then
\begin{align*}
\big|\widehat{\tau}_{\mathtt{LBATE}} - \tau_{\mathtt{LBATE}} \big|
\lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{ n h^d}} + \frac{\log(1/h)}{n^{\frac{1+v}{2+v}}h^d} + h^{p+1}.
\end{align*}
\end{thm}
Recall from Section~\ref{sa-sec: distributional approximation and inference} that $(\widehat{Z}(\mathbf{x}):\mathbf{x} \in \mathcal{B})$ is a (conditionally on $\mathbf{W}$) mean-zero Gaussian process with feasible (conditional) covariance function
\begin{align*}
\mathbb{C}\mathrm{ov} \Big[\widehat{Z}(\mathbf{x}_1),\widehat{Z}(\mathbf{x}_2) \Big|\mathbf{W} \Big]
= \frac{\widehat{\Omega}_{\mathbf{x}_1,\mathbf{x}_2}}
{\sqrt{\widehat{\Omega}_{\mathbf{x}_1,\mathbf{x}_1} \widehat{\Omega}_{\mathbf{x}_2,\mathbf{x}_2}}},
\qquad \mathbf{x}_1, \mathbf{x}_2 \in \mathcal{B}.
\end{align*}
Consider the confidence interval given by
\begin{align*}
\widehat{\operatorname{I}}_{\alpha,\mathtt{LBATE}} = \bigg[\sup_{\mathbf{b} \in \mathcal{B}} \Big(\widehat{\tau}(\mathbf{b}) - \mathcal{q}_{\alpha} \sqrt{\widehat{\Omega}_{\mathbf{b},\mathbf{b}}}\Big), \; \sup_{\mathbf{b} \in \mathcal{B}} \Big(\widehat{\tau}(\mathbf{b}) + \mathcal{q}_{\alpha} \sqrt{\widehat{\Omega}_{\mathbf{b},\mathbf{b}}}\Big)\bigg],
\end{align*}
where $\mathcal{q}_{\alpha} = \inf \big\{c > 0: \mathbb{P} \big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{Z}(\mathbf{x})\big|\geq c \big| \mathbf{W} \big) \leq \alpha \big\}$.
\begin{thm}[Confidence Interval: LBATE]\label{sa-thm: Confidence Bands for max}
Suppose the assumptions and conditions in Theorem \ref{sa-thm: Gaussian Strong Approximation: Tstat} hold. If $\liminf_{n \to \infty} \frac{\log h}{\log n} > - \infty$, $\frac{(\log n)^{3}}{n^{\frac{v}{2+v}}h^d} = o(1)$ and $h^{p+1}\sqrt{n h^d} = o(1)$, then
\begin{align*}
\mathbb{P}\Big[\tau_{\mathtt{LBATE}} \in \widehat{\operatorname{I}}_{\alpha,\mathtt{LBATE}} \Big] \geq 1 - \alpha + o(1).
\end{align*}
\end{thm}
\section{Gaussian Strong Approximation}\label{sa-sec: Gaussian Strong Approximation}
We present a Gaussian strong approximation theorem, which is the key technical tool behind Theorem~\ref{sa-thm: Gaussian Strong Approximation: Tstat}. The theorem builds on and generalizes the results in \cite{Cattaneo-Yu_2025_AOS}. Consider the \emph{residual-based empirical process} given by
\begin{align*}
R_n[g,r] = \frac{1}{\sqrt{n}} \sum_{i = 1}^n \Big[g(\mathbf{x}_i) r(y_i) - \mathbb{E}[g(\mathbf{x}_i) r(y_i)|\mathbf{x}_i] \Big], \qquad g \in \mathcal{G}, r \in \mathcal{R},
\end{align*}
where $\mathcal{G}$ and $\mathcal{R}$ are classes of functions satisfying certain regularity conditions.
\subsection{Definitions for Function Spaces}
Let $\mathcal{F}$ be a class of measurable functions from a probability space $(\mathbb{R}^q, \mathcal{B}(\mathbb{R}^q), \mathbb{P})$ to $\mathbb{R}$. We introduce several definitions that capture properties of $\mathcal{F}$.
\begin{enumerate}[label=(\roman*)]
\item $\mathcal{F}$ is pointwise measurable if it contains a countable subset $\mathcal{G}$ such that for any $f \in \mathcal{F}$, there exists a sequence $(g_m:m\geq1) \subseteq \mathcal{G}$ such that $\lim_{m \to \infty} g_m(\mathbf{u}) = f(\mathbf{u})$ for all $\mathbf{u} \in \mathbb{R}^q$.
\item Let $\operatorname{Supp}(\mathcal{F}) = \cup_{f \in \mathcal{F}}\operatorname{Supp}(f)$. A probability measure $\mathbb{Q}_\mathcal{F}$ on $(\mathbb{R}^q,\mathcal{B}(\mathbb{R}^q))$ is a surrogate measure for $\mathbb{P}$ with respect to $\mathcal{F}$ if
\begin{enumerate}[label=(\roman*)]
\item $\mathbb{Q}_\mathcal{F}$ agrees with $\mathbb{P}$ on $\operatorname{Supp}(\mathbb{P}) \cap \operatorname{Supp}(\mathcal{F})$.
\item $\mathbb{Q}_\mathcal{F}(\operatorname{Supp}(\mathcal{F}) \setminus \operatorname{Supp}(\mathbb{P})) = 0$.
\end{enumerate}
Let $\mathcal{Q}_\mathcal{F}=\operatorname{Supp}(\mathbb{Q}_\mathcal{F})$.
\item For $q=1$ and an interval $\mathcal{I}\subseteq\mathbb{R}$, the pointwise total variation of $\mathcal{F}$ over $\mathcal{I}$ is
\begin{align*}
\mathtt{pTV}_{\mathcal{F},\mathcal{I}} = \sup_{f \in \mathcal{F}} \sup_{P\geq1}\sup_{\mathcal{P}_P \in \mathcal{I}} \sum_{i = 1}^{P-1}|f(a_{i+1}) - f(a_i)|,
\end{align*}
where $\mathcal{P}_P=\{(a_1,\dots,a_P):a_1 \leq \cdots \leq a_P\}$ denotes the collection of all partitions of $\mathcal{I}$.
\item For a non-empty $\mathcal{C} \subseteq \mathbb{R}^q$, the total variation of $\mathcal{F}$ over $\mathcal{C}$ is
\begin{align*}
\mathtt{TV}_{\mathcal{F}, \mathcal{C}} = \inf_{\mathcal{U} \in \mathcal{O}(\mathcal{C})}\sup_{f \in \mathcal{F}} \sup_{\phi \in \mathscr{D}_{q}(\mathcal{U})} \int_{\mathbb{R}^q} f(\mathbf{u})\operatorname{div}(\phi)(\mathbf{u}) d \mathbf{u} / \lVert \left\lVert\phi\right\rVert_2 \rVert_{\infty},
\end{align*}
where $\mathcal{O}(\mathcal{C})$ denotes the collection of all open sets that contains $\mathcal{C}$, and $\mathscr{D}_{q}(\mathcal{U})$ denotes the space of infinitely differentiable functions from $\mathbb{R}^q$ to $\mathbb{R}^q$ with compact support contained in $\mathcal{U}$.
\item For a non-empty $\mathcal{C} \subseteq \mathbb{R}^q$, the local total variation constant of $\mathcal{F}$ over $\mathcal{C}$, is a positive number $\mathtt{K}_{\mathcal{F},\mathcal{C}}$ such that for any cube $\mathcal{D} \subseteq \mathbb{R}^q$ with edges of length $\ell$ parallel to the coordinate axises,
\begin{align*}
\mathtt{TV}_{\mathcal{F}, \mathcal{D} \cap \mathcal{C}} \leq \mathtt{K}_{\mathcal{F}, \mathcal{C}} \ell^{d-1}.
\end{align*}
\item For a non-empty $\mathcal{C} \subseteq \mathbb{R}^q$, the envelopes of $\mathcal{F}$ over $\mathcal{C}$ are
\begin{align*}
\mathtt{M}_{\mathcal{F},\mathcal{C}} = \sup_{\mathbf{u} \in \mathcal{C} }M_{\mathcal{F},\mathcal{C}}(\mathbf{u}),
\qquad M_{\mathcal{F},\mathcal{C}}(\mathbf{u}) = \sup_{f \in \mathcal{F}}|f(\mathbf{u})|,
\qquad \mathbf{u} \in \mathcal{C}.
\end{align*}
\item For a non-empty $\mathcal{C} \subseteq \mathbb{R}^q$, the Lipschitz constant of $\mathcal{F}$ over $\mathcal{C}$ is
\begin{align*}
\mathtt{L}_{\mathcal{F},\mathcal{C}} = \sup_{f \in \mathcal{F}}\sup_{\mathbf{u}_1, \mathbf{u}_2 \in \mathcal{C}} \frac{|f(\mathbf{u}_1) - f(\mathbf{u}_2)|}{\|\mathbf{u}_1 - \mathbf{u}_2\|_\infty}.
\end{align*}
\item For a non-empty $\mathcal{C} \subseteq \mathbb{R}^q$, the $L_1$ bound of $\mathcal{F}$ over $\mathcal{C}$ is
\begin{align*}
\mathtt{E}_{\mathcal{F},\mathcal{C}} = \sup_{f \in \mathcal{F}} \int_{\mathcal{C}} |f| d\mathbb{P}.
\end{align*}
\item For a non-empty $\mathcal{C} \subseteq \mathbb{R}^q$, the uniform covering number of $\mathcal{F}$ with envelope $M_{\mathcal{F},\mathcal{C}}$ over $\mathcal{C}$ is
\begin{align*}
\mathtt{N}_{\mathcal{F},\mathcal{C}}(\delta,M_{\mathcal{F},\mathcal{C}}) = \sup_{\mu} N(\mathcal{F},\left\lVert\cdot\right\rVert_{\mu,2},\delta \left\lVertM_{\mathcal{F},\mathcal{C}}\right\rVert_{\mu,2}),
\qquad \delta \in (0, \infty),
\end{align*}
where the supremum is taken over all finite discrete measures on $(\mathcal{C}, \mathcal{B}(\mathcal{C}))$. We assume that $M_{\mathcal{F},\mathcal{C}}(\mathbf{u})$ is finite for every $\mathbf{u} \in \mathcal{C}$.
\item For a non-empty $\mathcal{C} \subseteq \mathbb{R}^q$, the uniform entropy integral of $\mathcal{F}$ with envelope $M_{\mathcal{F},\mathcal{C}}$ over $\mathcal{C}$ is
\begin{align*}
J_\mathcal{C}(\delta, \mathcal{F}, M_{\mathcal{F},\mathcal{C}}) = \int_0^{\delta} \sqrt{1 + \log \mathtt{N}_{\mathcal{F},\mathcal{C}}(\varepsilon,M_{\mathcal{F},\mathcal{C}})} d \varepsilon,
\end{align*}
where it is assumed that $M_{\mathcal{F},\mathcal{C}}(\mathbf{u})$ is finite for every $\mathbf{u} \in \mathcal{C}$.
\item For a non-empty $\mathcal{C} \subseteq \mathbb{R}^q$, $\mathcal{F}$ is a VC-type class with envelope $M_{\mathcal{F},\mathcal{C}}$ over $\mathcal{C}$ if (i) $M_{\mathcal{F},\mathcal{C}}$ is measurable and $M_{\mathcal{F},\mathcal{C}}(\mathbf{u})$ is finite for every $\mathbf{u} \in \mathcal{C}$, and (ii) there exist $\mathtt{c}_{\mathcal{F},\mathcal{C}}>0$ and $\mathtt{d}_{\mathcal{F},\mathcal{C}}>0$ such that
\begin{align*}
\mathtt{N}_{\mathcal{F},\mathcal{C}}(\varepsilon,M_{\mathcal{F},\mathcal{C}}) \leq \mathtt{c}_{\mathcal{F},\mathcal{C}} \varepsilon^{-\mathtt{d}_{\mathcal{F},\mathcal{C}}},
\qquad \varepsilon\in(0,1).
\end{align*}
\end{enumerate}
If a surrogate measure $\mathbb{Q}_\mathcal{F}$ for $\mathbb{P}$ with respect to $\mathcal{F}$ has been assumed, and it is clear from the context, we drop the dependence on $\mathcal{C} = \mathcal{Q}_{\mathcal{F}}$ for all quantities in the previous definitions. That is, to save notation, we set $\mathtt{TV}_{\mathcal{F}}=\mathtt{TV}_{\mathcal{F},\mathcal{Q}_{\mathcal{F}}}$, $\mathtt{K}_{\mathcal{F}}=\mathtt{K}_{\mathcal{F},\mathcal{Q}_{\mathcal{F}}}$, $\mathtt{M}_{\mathcal{F}}=\mathtt{M}_{\mathcal{F},\mathcal{Q}_{\mathcal{F}}}$, $M_{\mathcal{F}}(\mathbf{u})=M_{\mathcal{F},\mathcal{Q}_{\mathcal{F}}}(\mathbf{u})$, $\mathtt{L}_{\mathcal{F}}=\mathtt{L}_{\mathcal{F},\mathcal{Q}_{\mathcal{F}}}$, and so on, whenever there is no confusion.
\subsection{Residual-based Empirical Process}
The following theorem generalizes \citet[Theorem 2]{Cattaneo-Yu_2025_AOS} by requiring only bounded polynomial moments for $y_i$ conditional on $\mathbf{x}_i$.
\begin{thm}[Strong Approximation for Residual-based Empirical Processes]\label{sa-lem: sa thm}
Suppose $(\mathbf{z}_i=(\mathbf{x}_i, y_i): 1 \leq i \leq n)$ are i.i.d. random vectors taking values in $(\mathbb{R}^{d+1}, \mathcal{B}(\mathbb{R}^{d+1}))$ with common law $\mathbb{P}_Z$, where $\mathbf{x}_i$ has distribution $\mathbb{P}_X$ supported on $\mathcal{X}\subseteq\mathbb{R}^d$, $y_i$ has distribution $\mathbb{P}_Y$ supported on $\mathcal{Y}\subseteq\mathbb{R}$, $\sup_{\mathbf{x} \in \mathcal{X}}\mathbb{E}[|y_i|^{2 + v}|\mathbf{x}_i = \mathbf{x}] \leq 2$ for some $v > 0$, and the following conditions hold:
\begin{enumerate}[label=\emph{(\roman*)}]
\item $\mathcal{G}$ is a real-valued pointwise measurable class of functions on $(\mathbb{R}^d, \mathcal{B}(\mathbb{R}^d), \mathbb{P}_X)$.
\item There exists a surrogate measure $\mathbb{Q}_\mathcal{G}$ for $\mathbb{P}_X$ with respect to $\mathcal{G}$ such that $\mathbb{Q}_\mathcal{G} = \mathfrak{m} \circ \phi_\mathcal{G}$, where the \textit{normalizing transformation} $\phi_{\mathcal{G}}: \mathcal{Q}_\mathcal{G} \mapsto [0,1]^d$ is a diffeomorphism.
\item $\mathcal{G}$ is a VC-type class with envelope $\mathtt{M}_{\mathcal{G}}$ over $\mathcal{Q}_\mathcal{G}$ with $\mathtt{c}_{\mathcal{G}} \geq e$ and $\mathtt{d}_{\mathcal{G}} \geq 1$.
\item $\mathcal{R}$ is a real-valued pointwise measurable class of functions on $(\mathbb{R}, \mathsf{Borel}(\mathbb{R}),\mathbb{P}_Y)$.
\item $\mathcal{R}$ is a VC-type class with envelope $M_{\mathcal{R},\mathcal{Y}}$ over $\mathcal{Y}$ with $\mathtt{c}_{\mathcal{R},\mathcal{Y}}\geq e$ and $\mathtt{d}_{\mathcal{R},\mathcal{Y}}\geq 1$, where $M_{\mathcal{R},\mathcal{Y}}(y) + \mathtt{pTV}_{\mathcal{R},(-|y|,|y|)} \leq \mathtt{v} (1 + |y|)$ for all $y \in \mathcal{Y}$, for some $\mathtt{v}>0$.
\item There exists a constant $\mathtt{k}$ such that $|\log_2 \mathtt{E}_{\mathcal{G}}| + |\log_2 \mathtt{TV}| + |\log_2 \mathtt{M}_{\mathcal{G}}| \leq \mathtt{k} \log_2 n$, where the constant $\mathtt{TV} = \max \{\mathtt{TV}_{\mathcal{G}}, \mathtt{TV}_{\mathcal{G} \times \mathscr{U}_{\mathcal{R}},\mathcal{Q}_\mathcal{G}}\}$ with $\mathscr{U}_{\mathcal{R}} = \{\theta(\cdot,r, \tau) : r \in \mathcal{R}, \tau \in (0,\infty]\}$, and $\theta(\mathbf{x},r,\tau) = \mathbb{E}[r(y_i) \mathds{1}(|y_i| \leq \tau)|\mathbf{x}_i = \mathbf{x}]$ for $\mathbf{x} \in \mathcal{X}$.
\end{enumerate}
Define the residual based empirical process
\begin{align*}
R_n(g,r) = \frac{1}{\sqrt{n}} \sum_{i = 1}^n g(\mathbf{x}_i)(r(y_i) - \mathbb{E}[r(y_i)|\mathbf{x}_i]), \qquad g \in \mathcal{G}, r \in R.
\end{align*}
Then (1) on a possibly enlarged probability space, there exists a sequence of mean-zero Gaussian processes $(Z_n^R(g,r): g\in\mathcal{G}, r \in \mathcal{R})$ with almost sure continuous trajectories such that:
\begin{itemize}
\item $\mathbb{E}[R_n(g_1, r_1) R_n(g_2, r_2)] = \mathbb{E}[Z^R_n(g_1, r_1) Z^R_n(g_2, r_2)]$ for all $(g_1, r_1), (g_2, r_2) \in \mathcal{G} \times \mathcal{R}$, and
\item $\mathbb{E}\big[\left\lVertR_n - Z_n^R\right\rVert_{\mathcal{G} \times \mathcal{R}}\big] \leq C \mathtt{v} \mathtt{d} \log(\mathtt{c} n) \rho_n$,
\end{itemize}
with
\begin{align*}
\rho_n = \sqrt{\mathtt{d} \log (\mathtt{c} n)} \; \mathtt{r}_n^{\frac{v}{v +2}}(\sqrt{\mathtt{M}_{\mathcal{G}}\mathtt{E}_{\mathcal{G}}})^{\frac{2}{v+2}} + \mathtt{M}_{\mathcal{G}} n^{-\frac{v/2}{2+v}} + \frac{\mathtt{M}_{\mathcal{G}}}{\sqrt{n}} \Big(\frac{\sqrt{\mathtt{M}_{\mathcal{G}} \mathtt{E}_{\mathcal{G}}}}{\mathtt{r}_n}\Big)^{\frac{2}{v+2}},
\end{align*}
where $C$ is a positive universal constant, $\mathtt{c} = \mathtt{c}_{\mathcal{G}} + \mathtt{c}_{\mathcal{R},\mathcal{Y}} + \mathtt{k}$, $\mathtt{d} = \mathtt{d}_{\mathcal{G}} \mathtt{d}_{\mathcal{R},\mathcal{Y}} \mathtt{k}$, and
\begin{gather*}
\mathtt{r}_n = \min\Big\{\frac{(\mathtt{c}_1^d \mathtt{M}_{\mathcal{G}}^{d+1} \mathtt{TV}^d \mathtt{E}_{\mathcal{G}})^{1/(2d+2)} }{n^{1/(2d+2)}}, \frac{(\mathtt{c}_1^{d/2} \mathtt{c}_2^{d/2} \mathtt{M}_{\mathcal{G}} \mathtt{TV}^{d/2} \mathtt{E}_{\mathcal{G}} \mathtt{L}^{d/2})^{1/(d+2)}}{n^{1/(d+2)}} \Big\}, \\
\mathtt{c}_1 = d \sup_{\mathbf{x} \in \mathcal{Q}_{\mathcal{G}}} \prod_{j = 1}^{d-1} \sigma_j(\nabla \phi_{\mathcal{G}}(\mathbf{x})), \qquad
\mathtt{c}_2 = \sup_{\mathbf{x} \in \mathcal{Q}_{\mathcal{G}}} \frac{1}{\sigma_{d}(\nabla \phi_{\mathcal{G}}(\mathbf{x}))}, \qquad
\mathtt{L} = \max\{\mathtt{L}_{\mathcal{G}}, \mathtt{L}_{\mathcal{G} \times \mathscr{U}_{\mathcal{R}},\mathcal{Q}_{\mathcal{R}}}\};
\end{gather*}
and (2) if $\mathcal{R}$ is a singleton, then we can replace $\mathtt{TV}$ and $\mathtt{L}$ in the previous conditions and statements by $\mathtt{TV}_{\text{sing}} = \max\{\mathtt{TV}_{\mathcal{G}}, \mathtt{TV}_{\mathcal{G} \times \mathcal{V}_{\mathcal{R}}, \mathcal{Q}_{\mathcal{G}}}\}$, and $\mathtt{L}_{\text{sing}} = \max\{\mathtt{L}_{\mathcal{G}}, \mathtt{L}_{\mathcal{G} \times \mathscr{\mathcal{V}}_{\mathcal{R}},\mathcal{Q}_{\mathcal{R}}}\}$, respectively, with $\mathcal{V}_{\mathcal{R}} = \{\theta(\cdot,r) : r \in \mathcal{R}\}$, and $\theta(\mathbf{x},r) = \mathbb{E}[r(y_i) |\mathbf{x}_i = \mathbf{x}]$ for $\mathbf{x} \in \mathcal{X}$.
\end{thm}
\begin{remark}
The class $\mathscr{U}_{\mathcal R}$ comprises truncated conditional means at all truncation levels. Its Lipschitz and total-variation constants can be bounded, for example, if $f(y \mid \mathbf{x})$ is Lipschitz in $\mathbf{x}$ uniformly over $(\mathbf{x},y)$ in the support of $(\mathbf{x}_i,y_i)$. When $\mathcal R$ is a singleton, it suffices to assume regularity only for $\mathscr{V}_{\mathcal R}$, the class containing the (untruncated) conditional mean functions, which is easily justified.
\end{remark}
\section{Proofs}\label{sa-sec: Proofs}
\subsection{Proof of Lemma~\ref{sa-lem:invert}}
Assumption~\ref{sa-assump: DGP} (ii) implies
\begin{align*}
\boldsymbol{\Gamma}_{t,\mathbf{x}} & = \mathbb{E} \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)^{\top} K_h(\mathbf{X}_i - \mathbf{x}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) \Big] \\
& = \int_{\mathcal{A}_t} \mathbf{r}_p \Big(\frac{\mathbf{u} - \mathbf{x}}{h}\Big) \mathbf{r}_p \Big(\frac{\mathbf{u} - \mathbf{x}}{h}\Big)^{\top} K_h(\mathbf{u} - \mathbf{x}) f(\mathbf{u}) d \mathbf{u} \\
& = f(\mathbf{x})\int_{\mathcal{A}_t} \mathbf{r}_p \Big(\frac{\mathbf{u} - \mathbf{x}}{h}\Big) \mathbf{r}_p \Big(\frac{\mathbf{u} - \mathbf{x}}{h}\Big)^{\top} K_h(\mathbf{u} - \mathbf{x}) d \mathbf{u} + o(1),
\end{align*}
where in the last line we have used $\int_{\mathcal{A}_t} (\frac{\mathbf{u} - \mathbf{x}}{h})^{\mathbf{v}} K_h(\mathbf{u} - \mathbf{x}) d \mathbf{u} = O(1)$ for any multi-index $\mathbf{v}$ from standard change of variable argument.
\begin{center}
\textbf{I. Polynomial Representation of Minimum Eigenvalue}
\end{center}
For simplicity, call
\begin{align*}
\mathbf{S}_{t,\mathbf{x}} = \lim_{h \to 0} \mathbf{S}_{t,\mathbf{x}}(h), \qquad \mathbf{S}_{t,\mathbf{x}}(h) = \int_{\mathcal{A}_t} \mathbf{r}_p \Big(\frac{\mathbf{u} - \mathbf{x}}{h}\Big) \mathbf{r}_p \Big(\frac{\mathbf{u} - \mathbf{x}}{h}\Big)^{\top} K_h(\mathbf{u} - \mathbf{x}) d \mathbf{u}.
\end{align*}
A change of variable gives
\begin{align*}
\mathbf{S}_{t,\mathbf{x}}(h) = \int \mathbf{r}_p(\mathbf{z}) \mathbf{r}_p(\mathbf{z})^{\top} K(\mathbf{z}) \mathds{1}(\mathbf{x} + h \mathbf{z} \in \mathcal{A}_t) d \mathbf{z}.
\end{align*}
Let $\mathbf{a} \in \mathbb{R}^{\mathfrak{p}_p}$, where $\mathfrak{p}_p = \frac{(d + p)!}{d! p !}$.
Then the equivalent representation of minimum eigenvalue gives
\begin{align}\label{sa-eq: min eigenvalue}
\nonumber \lambda_{\min}(\mathbf{S}_{t,\mathbf{x}}(h))
& = \min_{\left\lVert\mathbf{a}\right\rVert = 1} \int (\mathbf{a}^{\top} \mathbf{r}_p(\mathbf{z}))^2 K(\mathbf{z}) \mathds{1}(\mathbf{x} + h \mathbf{z} \in \mathcal{A}_t) d \mathbf{z} \\
& \geq \kappa \min_{\left\lVert\mathbf{a}\right\rVert = 1} \int_{U} (\mathbf{a}^{\top} \mathbf{r}_p(\mathbf{z}))^2 \mathds{1}(\mathbf{x} + h \mathbf{z} \in \mathcal{A}_t) d \mathbf{z},
\end{align}
where in the last line we have used $K(\mathbf{u}) \geq \kappa$ for all $u \in U$.
\begin{center}
\textbf{II. Mass Retaining Ratio in Treatment/Control Region}
\end{center}
Denote $E_h(\mathbf{x},t) = \{\mathbf{z} \in U: \mathbf{x} + h \mathbf{z} \in \mathcal{A}_t \}$. Assumption~\ref{sa-assump: Kernel and Boundary} (ii) implies there is some upper bound $\Lambda > 0$ of $K(\cdot)$. Hence for $c_0 = 1/2 \; \liminf_{h \downarrow 0}\inf_{\mathbf{x} \in \mathcal{B}} \int_{U} K(\mathbf{u}) \mathds{1}(\mathbf{x} + h \mathbf{u} \in \mathcal{A}_t) d \mathbf{u}$, we have
\begin{align*}
\Lambda \mathfrak{m}(E_h(\mathbf{x},t))
& \geq \int_{U} K(\mathbf{u}) \mathds{1}(\mathbf{x} + h \mathbf{u} \in \mathcal{A}_t)
\geq c_0
\end{align*}
for small enough $h$, which implies
\begin{align}\label{sa-eq: mass relation}
\mathfrak{m}(E_h(\mathbf{x},t)) \geq \alpha \mathfrak{m}(U), \qquad \alpha = \frac{c_0}{\Lambda \mathfrak{m}(U)}.
\end{align}
\begin{center}
\textbf{III. $L_2$ Integral of Polynomials in Full v.s. Treatment/Control Regions}
\end{center}
Consider $S = \{f \in \mathcal{P}_{\mathfrak{p}}: \int_U f(\mathbf{u})^2 d \mathbf{u} = 1\}$, where $\mathcal{P}_{\mathfrak{p}}$ is the collection of all $\mathfrak{p}$-order polynomials. Let $(\phi_j, 1 \leq j \leq \mathfrak{p})$ be a set of orthonormal basis of $(\mathcal{P}_{\mathfrak{p}}, \lVert \cdot \rVert_{L_2})$. Then $T(\mathbf{a}) = \sum_{j = 1}^{\mathfrak{p}} a_j \phi_j$ is an isometry. Since $T(S) = \{\mathbf{a} \in \mathbb{R}^{\mathfrak{p}}: \left\lVert\mathbf{a}\right\rVert = 1\}$ is compact, $S$ is also compact in $(\mathcal{P}_{\mathfrak{p}}, \lVert \cdot \rVert_{L_2})$. Since $\mathcal{P}_{\mathfrak{p}}$ is $\mathfrak{p}$-dimensional, equivalent of norms implies that $S$ is also compact in $(\mathcal{P}_{\mathfrak{p}}, \lVert \cdot \rVert_{L_{\infty}})$. Now consider
\begin{align*}
\Phi_q(\varepsilon) = \mathfrak{m}(\{\mathbf{u} \in U: |q(u)| < \varepsilon \}), \qquad q \in S, \varepsilon > 0,
\end{align*}
and
\begin{align*}
\psi(q) = \sup \Big\{\varepsilon > 0: \Phi_q(\varepsilon) \leq \frac{\alpha}{2} \mathfrak{m}(U) \Big\}.
\end{align*}
Since $\int_U q^2 = 1$ and $q$ is polynomial, $\lim_{\varepsilon \downarrow 0} \Phi_q(\varepsilon) = 0$ and $\Phi_q(\lVert q\rVert_{\infty}) = \mathfrak{m}(U)$. Continuity and Lipchitzness of $q \in S$ imply $\psi(q) > 0$ for all $q \in S$.
Next, we want to show $\psi$ is lower-semicontinous function on $(\mathcal{P}_{\mathfrak{p}}, \lVert \cdot \rVert_{L_{\infty}})$. Suppose $q_n \to q$ uniformly on $U$. For every $\varepsilon_0 \in (0, \psi(q))$, there exists $\eta > 0$ such that $\Phi_q(\varepsilon_0) \leq \frac{\alpha}{2} \mathfrak{m}(U) - \eta$. Continuity of polynomials and the fact that level sets of polynomials have zero Lebesgue measure imply $\mathds{1}_{\{|q_n| < \varepsilon_0\}}(\cdot) \to \mathds{1}_{\{|q| < \varepsilon_0\}}(\cdot)$ almost surely. By Dominated Convergence Theorem, $\Phi_{q_n}(\varepsilon_0) \to \Phi_q(\varepsilon_0)$. Hence for large enough $n$, $\Phi_{q_n}(\varepsilon_0) \leq \frac{\alpha}{2}\mathfrak{m}(U)$, which implies $\varepsilon_0 \leq \psi(q_n)$. This implies $\liminf_{n \to \infty} \psi(q_n) \geq \varepsilon_0$. Since $\varepsilon_0$ is arbitrary in $(0, \psi(q))$, we have $\liminf_{n \to \infty} \psi(q_n) \geq \psi(q)$.
Compactness of $S$ and lower-semicontinuity of $\psi$ implies $\psi$ attains its minimum on $S$. Since $\psi(q) > 0$ for all $q \in S$, we know $\varepsilon_* = \inf_{q \in S} \psi(q) > 0$. Then for every $q \in S$,
\begin{align*}
\int_{E_h(\mathbf{x},t)} q^2
& \geq \varepsilon_*^2 \; \mathfrak{m} \Big(E_h(\mathbf{x},t) \setminus \{|q| \leq \varepsilon_*\} \Big) \\
& \geq \varepsilon_*^2 \; \Big(\mathfrak{m}(E_h(\mathbf{x},t)) - \mathfrak{m}(\{|q| \leq \varepsilon_*\}) \Big) \\
& \geq \varepsilon_*^2 \; \frac{\alpha}{2} \mathfrak{m}(U).
\end{align*}
Scaling $q$ from $S$ gives
\begin{align}\label{sa-eq: integral relation}
\int_{E_h(\mathbf{x},t)} q^2 \geq \varepsilon_*^2 \; \frac{\alpha}{2} \int_{U} q^2, \qquad q \in \mathcal{P}_{\mathfrak{p}}.
\end{align}
\begin{center}
\textbf{IV. Lower Bound of Minimum Eigenvalue}
\end{center}
Equations~\eqref{sa-eq: min eigenvalue}, \eqref{sa-eq: mass relation} and \eqref{sa-eq: integral relation} together give for small enough $h$,
\begin{align*}
\inf_{\mathbf{x} \in \mathcal{B}}\lambda_{\min}(\mathbf{S}_{t,\mathbf{x}}(h))
& \geq \kappa \inf_{\mathbf{x} \in \mathcal{B}}\min_{\left\lVert\mathbf{a}\right\rVert = 1} \int_{E_h(\mathbf{x},t)} (\mathbf{a}^{\top} \mathbf{r}_p(\mathbf{z}))^2 d \mathbf{z}, \\
& \geq \kappa \varepsilon_*^2 \; \frac{\alpha}{2} \min_{\left\lVert\mathbf{a}\right\rVert = 1} \int_{U} (\mathbf{a}^{\top} \mathbf{r}_p(\mathbf{z}))^2 d \mathbf{z} \\
& \geq \kappa \varepsilon_*^2 \; \frac{\alpha}{2} \lambda_{\min} \Big(\int_U \mathbf{r}_p(\mathbf{z}) \mathbf{r}_p(\mathbf{z})^{\top} d \mathbf{z} \Big),
\end{align*}
which implies $\liminf_{h \to 0} \inf_{\mathbf{x} \in \mathcal{B}}\lambda_{\min}(\mathbf{S}_{t,\mathbf{x}}(h)) > 0$.
\subsection{Proof of Lemma~\ref{sa-lem: gram}}
Since $\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}$ is a finite dimensional matrix, it suffices to show the stated rate of convergence for each entry. Let $\mathbf{v}$ be a multi-index such that $|\mathbf{v}| \leq 2 p$. Define
\begin{align*}
g_n(\xi,\mathbf{x}) = \left(\frac{\xi - \mathbf{x}}{h}\right)^{\mathbf{v}} \frac{1}{h^d} K \left(\frac{\xi - \mathbf{x}}{h}\right)\mathds{1}(\xi \in \mathcal{A}_t), \qquad \xi \in \mathcal{X}, \mathbf{x} \in \mathcal{B}.
\end{align*}
and $\mathcal{F} = \{g_n(\cdot, \mathbf{x}): \mathbf{x} \in \mathcal{B}\}$. We will show $\mathcal{F}$ is a VC-type of class. In order to do this, we study the following quantities.
\medskip\textit{Constant Envelope Function}. We assume $K$ is continuous and has compact support, or $K = \mathds{1}(\cdot \in [-1,1]^d)$. Hence there exists a constant $C_1$ such that for all $l \in \mathcal{F}$, for all $\mathbf{x} \in \mathcal{B}$, $|l(\mathbf{x})| \leq C_1 h^{-d} = F$.
\medskip\textit{Diameter of $\mathcal{F}$ in $L_2$}. $ \sup_{l \in \mathcal{F}} \left\lVertl\right\rVert_{\mathbb{P},2}
= \sup_{\mathbf{x} \in \mathcal{B}}(\int_{\frac{\mathcal{A}_t - \mathbf{x}}{h}} \frac{1}{h^d}\mathbf{y}^{2\mathbf{v}}K(\mathbf{y})^2 f_X(\mathbf{x} + h \mathbf{y}) d \mathbf{y})^{1/2} \leq C_2 h^{-d/2}$ for some constant $C_2$. We can take $C_1$ large enough so that $\sigma = C_2 h^{-d/2} \leq F = C_1 h^{-d}$.
\medskip\textit{Ratio}. For some constant $C_3$, $\delta = \frac{\sigma}{F} = C_3 \sqrt{h^d}$.
\medskip\textit{Covering Numbers}. Case 1: $K$ is Lipschitz. Let $\mathbf{x},\mathbf{x}' \in \mathcal{B}$. Then, for a generic evaluation points $\mathbf{x} = (x_1,\ldots,x_d)^\top$ and $\mathbf{x}' = (x_1',\ldots,x_d')^\top$,
\begin{align*}
\sup_{\xi \in \mathcal{X}} \left|g_n(\xi,\mathbf{x}) - g_n(\xi, \mathbf{x}') \right|
& \leq \left|\Big(\frac{\xi_1 - x_1}{h}\Big)^{v_1} \cdots \Big(\frac{\xi - x_d}{h}\Big)^{v_d} - \Big(\frac{\xi_1 - x_1'}{h}\Big)^{v_1} \cdots \Big(\frac{\xi - x_d'}{h}\Big)^{v_d} \right| K_h(\xi - \mathbf{x}) \\
& \qquad + \Big(\frac{\xi_1 - x_1'}{h}\Big)^{v_1} \cdots \Big(\frac{\xi - x_d'}{h}\Big)^{v_d} \Big|K_h(\xi - \mathbf{x}) - K_h(\xi - \mathbf{x}') \Big| \\
& \lesssim h^{-d-1} \lVert \mathbf{x} - \mathbf{x}' \rVert_{\infty},
\end{align*}
since we have assumed that $K$ has compact support and is Lipschitz continuous. Hence, for any $\varepsilon \in (0,1]$ and for any finitely supported measure $Q$ and metric $\left\lVert\cdot\right\rVert_{Q,2}$ based on $L_2(Q)$,
\begin{align*}
N(\mathcal{F}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon \left\lVertF\right\rVert_{Q,2})
\leq N(\mathcal{X}, \lVert \cdot \rVert_{\infty}, \varepsilon \left\lVertF\right\rVert_{Q,2} h^{d+1})
\stackrel{(i)}{\lesssim} \left(\frac{\operatorname{diam}(\mathcal{X})}{ \varepsilon \left\lVertF\right\rVert_{Q,2} h^{d+1}}\right)^d
\lesssim \left(\frac{\operatorname{diam}(\mathcal{X})}{\varepsilon h}\right)^d,
\end{align*}
where in ($i$) we used the fact that $\varepsilon \left\lVertF\right\rVert_{Q,2} h^{d+1} \lesssim \varepsilon h \lesssim 1$.
Hence, $\mathcal{F}$ forms a VC-type class, and taking $A_1 = \operatorname{diam}(\mathcal{X})/h$ and $A_2 = d$, $\sup_Q N(\mathcal{F}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon \left\lVertF\right\rVert_{Q,2}) \lesssim (A_1/ \varepsilon)^{A_2}$, $\varepsilon \in (0,1]$, and where the supremum is over all finite discrete measure.
Case 2: $K = \mathds{1}(\cdot \in [-1,1]^d)$. Consider
\begin{align*}
m_n(\xi, \mathbf{x})
= \Big(\frac{\xi - \mathbf{x}}{h}\Big)^{\mathbf{v}} \frac{1}{h^d} \mathds{1}(\xi \in \mathcal{A}_t),
\qquad \xi, \mathbf{x} \in \mathcal{X},
\end{align*}
$\mathcal{M} = \{m_n(\cdot, \mathbf{x}): \mathbf{x} \in \mathcal{B}\}$ and the constant envelope function $M = C_4 h^{-|\mathbf{v}| - d}$, for some constant $C_4$ only depending on diameter of $\mathcal{X}$. The same argument as before shows that for any discrete measure $Q$, we have
\begin{align*}
N(\mathcal{M}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon \left\lVertM\right\rVert_{Q,2})
& \leq N(\mathcal{X}, \lVert \cdot \rVert_{\infty}, \varepsilon \left\lVertM\right\rVert_{Q,2} h^{d + |\mathbf{v}|+1})
\lesssim \Big(\frac{\operatorname{diam}(\mathcal{X})}{\varepsilon \left\lVertM\right\rVert_{Q,2} h^{d + |\mathbf{v}|+1}}\Big)^d
\lesssim \Big(\frac{\operatorname{diam}(\mathcal{X})}{\varepsilon h}\Big)^d.
\end{align*}
The class $\mathcal{G} = \{\mathds{1}(\cdot - \mathbf{x} \in [-1,1]^d): \mathbf{x} \in \mathcal{B}\}$ has VC dimension no greater than $2d$ \cite[Example 2.6.1]{van-der-Vaart-Wellner_1996_Book}, and by \citet[Theorem 2.6.4]{van-der-Vaart-Wellner_1996_Book}, for any discrete measure $Q$, $N(\mathcal{G}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon) \leq 2d (4 e)^{2d} \varepsilon^{-4d}$, $0 < \varepsilon \leq 1$. It then follows that for any discrete measure $Q$,
\begin{align*}
N(\mathcal{F}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon \left\lVertH\right\rVert_{Q,2})
\lesssim N(\mathscr{H}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon/2 \left\lVertH\right\rVert_{Q,2}) + N(\mathcal{G}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon/2)
\lesssim 2^d h^{-d} \varepsilon^{-d} + 2d (32 e)^{d} \varepsilon^{-4d}.
\end{align*}
Hence, taking $A_1 = (2^d h^{-d} + 2d (32 e)^d) h^{-|\mathbf{v}|}$ and $A_2 = 4d$, $\sup_Q N(\mathcal{F}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon \left\lVertF\right\rVert_{Q,2}) \lesssim (A_1/ \varepsilon)^{A_2}$, $\varepsilon \in (0,1]$, where the supremum is over all finite discrete measure.
\medskip\textit{Maximal Inequality}. By Corollary 5.1 in \cite{Chernozhukov-Chetverikov-Kato_2014b_AoS} for the empirical process on class $\mathcal{F}$,
\begin{align*}
\mathbb{E} \Big[\sup_{\mathbf{x} \in \mathcal{B}} \big|\mathbb{E}_n [g_n(\mathbf{X}_i, \mathbf{x})] - \mathbb{E}[g_n(\mathbf{X}_i, \mathbf{x})]\big| \Big]
& \lesssim \frac{\sigma}{\sqrt{n}}\sqrt{A_2\log(A_1/\delta)} + \frac{\lVert F \rVert_{\mathbb{P},2} A_2\log(A_1/\delta)}{n} \\
& \lesssim \sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n h^d},
\end{align*}
where $A_1, A_2, \sigma, F, \delta$ are all given previously. Assuming $\frac{\log(h^{-1})}{n h^d} \to 0$ as $n \to \infty$, we conclude that $\sup_{\mathbf{x} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}} - \boldsymbol{\Gamma}_{t,\mathbf{x}}\big\| \lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{n h^d}}$. Hence, $1 \lesssim_{\mathbb{P}} \inf_{\mathbf{x} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}\big\| \lesssim_{\mathbb{P}} \sup_{\mathbf{x} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}\big\| \lesssim_{\mathbb{P}} 1$. By Weyl's Theorem, $\sup_{\mathbf{x} \in \mathcal{B}}|\lambda_{\min}(\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}) - \lambda_{\min}(\boldsymbol{\Gamma}_{t,\mathbf{x}})| \leq \sup_{\mathbf{x} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}} - \boldsymbol{\Gamma}_{t,\mathbf{x}}\big\| \lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{n h^d }}$. Assuming that $\lambda_{\min}(\boldsymbol{\Gamma}_{t,\mathbf{x}}) \gtrsim 1$ (which we will verify in the last part of the proof), then we can lower the minimum eigenvalue by $\inf_{\mathbf{x} \in \mathcal{B}} \lambda_{\min}(\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}) \geq \inf_{\mathbf{x} \in \mathcal{B}}\lambda_{\min}(\boldsymbol{\Gamma}_{t,\mathbf{x}}) - \sup_{\mathbf{x} \in \mathcal{B}}|\lambda_{\min}(\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}) - \lambda_{\min}(\boldsymbol{\Gamma}_{t,\mathbf{x}})| \gtrsim_{\mathbb{P}} 1$. It follows that $\sup_{\mathbf{x} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1}\big\| \lesssim_{\mathbb{P}} 1$ and hence $\sup_{\mathbf{x} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} - \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1}\big\| \leq \sup_{\mathbf{x} \in \mathcal{B}} \big\|\boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1}\big\| \big\|\boldsymbol{\Gamma}_{t,\mathbf{x}} -\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}\big\| \big\|\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1}\big\| \lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{n h^{d}}}$.
\subsection{Proof of Lemma~\ref{sa-lem: Q}}
The proof is similar to the proof of Lemma~\ref{sa-lem: gram}. Let $\mathbf{v}$ be a multi-index such that $ 0 \leq |\mathbf{v}| \leq p$. Let
\begin{align*}
g_n(\xi,\mathbf{x}) = \Big(\frac{\xi - \mathbf{x}}{h}\Big)^{\mathbf{v}} K_h(\xi - \mathbf{x})\mathds{1}(\xi\in \mathcal{A}_t), \qquad \xi , \mathbf{x} \in \mathcal{X}.
\end{align*}
Define the class of functions $\mathcal{F} = \{(\xi, u) \in \mathcal{X} \times \mathbb{R} \mapsto g_n(\xi,\mathbf{x}): \mathbf{x} \in \mathcal{B}\}$.
\medskip\textit{Envelope Function}. Since $K$ is continuous on its compact support, there exists a constant $C_1 > 0$ such that $|g_n(\xi,\mathbf{x})u| \leq C_1 h^{-d} |u|$, for $\xi, \mathbf{x} \in \mathcal{X}$ and $ u \in \mathbb{R}$. We define the envelope function $F(\xi, u) = C_1 h^{-d}|u|$, for $\xi \in \mathcal{X}$ and $ u \in \mathbb{R}$. Moreover, by Assumption~\ref{sa-assump: DGP}(v), let $M = \max_{1 \leq i \leq n} F(\mathbf{X}_i, u_i)$, then
\begin{align*}
\mathbb{E}[M^2]^{1/2}
\lesssim h^{-d} \mathbb{E} \big[\max_{1 \leq i \leq n} |u_i|^2 \big]^{1/2}
\lesssim h^{-d} \mathbb{E} \big[\max_{1 \leq i \leq n} |u_i|^{2+v}\big]^{1/(2+v)}
\lesssim n^{1/(2+v)} h^{-d}.
\end{align*}
\medskip\textit{Diameter of $\mathcal{F}$ in $L_2$}. Recall we denote $u_i = Y_i - \mathbb{E}[Y_i|\mathbf{X}_i]$, then
\begin{align*}
\sup_{l \in \mathcal{F}} \mathbb{E}[l(\mathbf{X}_i, u_i)^2]^{1/2}
\leq \sup_{\xi \in \mathcal{X}} \mathbb{E}[u_i^2 | \mathbf{X}_i = \xi]^{1/2} \sup_{\xi \in \mathcal{X}} \mathbb{E}[g_n(\mathbf{X}_i, \xi)^2]^{1/2} \leq C_3 h^{-d/2}
= \sigma.
\end{align*}
\medskip\textit{Ratio}. We set $\delta = \frac{\sigma}{\left\lVertF\right\rVert_{\mathbb{P},2}} \lesssim h^{d/2}$.
\medskip\textit{Covering Numbers}. Case 1: $K$ is Lipschitz. Let $\mathbb{Q}$ be a finite distribution on $(\mathcal{X} \times \mathbb{R}, \mathcal{B}(\mathcal{X}) \otimes \mathsf{Borel}(\mathbb{R}))$. Let $\mathbf{x}, \mathbf{x}' \in \mathcal{X}$. In the proof of Lemma~\ref{sa-lem: gram}, we showed that $\sup_{\xi \in \mathcal{X}} \sup_{\mathbf{x}, \mathbf{x}' \in \mathcal{X} } \frac{|g_n(\xi,\mathbf{x}) - g_n(\xi, \mathbf{x}')|}{\lVert \mathbf{x} - \mathbf{x}' \rVert_{\infty}} \lesssim h^{-d-1}$. Hence,
\begin{align*}
\| g_n(\mathbf{X}_i, \mathbf{x}) u_i - g_n(\mathbf{X}_i, \mathbf{x}')u_i \|_{\mathbb{Q},2}
\leq \lVert g_n(\cdot, \mathbf{x}) - g_n(\cdot,\mathbf{x}') \rVert_{\infty} \left\lVertu_i\right\rVert_{\mathbb{Q},2}
\lesssim h^{-1} \left\lVertF\right\rVert_{\mathbb{Q},2} \lVert \mathbf{x}- \mathbf{x}' \rVert_{\infty}.
\end{align*}
It follows that $\sup_{Q} N(\mathcal{F}, \lVert \cdot \rVert_{Q,2}, \epsilon \| F \|_{Q,2}) \lesssim (\frac{\operatorname{diam}(\mathcal{X})}{\epsilon h}))^d$, where $\sup$ is over all finite probability distributions on $(\mathcal{X} \times \mathbb{R}, \mathcal{B}(\mathcal{X}) \otimes \mathsf{Borel}(\mathbb{R}))$. Letting $A_1 = \frac{\operatorname{diam}(\mathcal{X})}{h}$ and $ A_2 = d$, we conclude that
\begin{align*}
\sup_Q N(\mathcal{F}, \left\lVert\cdot\right\rVert_{Q,2}, \epsilon \left\lVertF\right\rVert_{Q,2}) \lesssim (A_1/\epsilon)^{A_2}, \qquad \epsilon \in (0,1].
\end{align*}
Case 2: $K$ is the uniform kernel. Let
\begin{align*}
m_n(\xi, \mathbf{x}) = \left(\frac{\xi - \mathbf{x}}{h}\right)^{\mathbf{v}} \frac{1}{h^d} \mathds{1}(\xi \in \mathcal{A}_t), \qquad \xi, \mathbf{x} \in \mathcal{X},
\end{align*}
with $\mathcal{M} = \{(\xi,u) \in \mathcal{X} \times \mathbb{R} \to m_n(\xi,\mathbf{x})u: \mathbf{x} \in \mathcal{B}\}$ and envelop function $M(\mathbf{x},u) = C_1 h^{-d-|\mathbf{v}|} |u|$, for a positive constant $C_1$ depending only on $K$. By similar arguments as Case 1 and the proof of Lemma~\ref{sa-lem: gram}, it follows that $\sup_Q N(\mathcal{M}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon \|M\|_{Q,2}) \lesssim \big(\frac{\operatorname{diam}(\mathcal{X})}{\varepsilon h}\big)^d$, where the supremum is taken over all finite discrete measures. Taking $\mathcal{G} = \{\mathds{1}(\cdot - \mathbf{x} \in [-1,1]^d): \mathbf{x} \in \mathcal{B}\}$, the proof of Lemma~\ref{sa-lem: gram} shows that
\begin{align*}
\sup_Q N(\mathcal{G}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon) \leq 2d (4 e)^{2d} \varepsilon^{-4d}, \qquad \varepsilon \in [0,1],
\end{align*}
where the supremum is taken over all finite discrete measures. Taking $A_1 = (2^d h^{-d} + 2d (32 e)^d) h^{-|\mathbf{v}|}$ and $A_2 = 4d$, we have
\begin{align*}
\sup_Q N(\mathcal{F}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon \left\lVertF\right\rVert_{Q,2}) \lesssim (A_1/ \varepsilon)^{A_2}, \qquad \varepsilon \in (0,1],
\end{align*}
the supremum is over all finite discrete measure.
\medskip\textit{Maximal Inequality}. By Corollary 5.1 in \cite{Chernozhukov-Chetverikov-Kato_2014b_AoS},
\begin{align*}
\mathbb{E} \Big[\sup_{\mathbf{x} \in \mathcal{X}} \big|\mathbb{E}_n [g_n(\mathbf{X}_i,\mathbf{x})u_i] \big| \Big]
& \lesssim \frac{\sigma}{\sqrt{n}}\sqrt{A_2\log(A_1/\delta)} + \frac{\lVert M \rVert_{\mathbb{P},2} A_2\log(A_1/\delta)}{n}\\
& \lesssim \sqrt{\frac{\log(1/h)}{n h^d}} +\frac{\log(1/h)}{n^{\frac{1+v}{2+v}} h^d}.
\end{align*}
Since $\mathbf{Q}_{t,\mathbf{x}}$ is finite-dimensional, entry-wise convergence implies convergence in norm with the same rate. Hence, $\sup_{\mathbf{x} \in \mathcal{X}} \big\|\mathbf{Q}_{t,\mathbf{x}}\big\| \lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{1+v}{2+v}}h^d}$. By Lemma~\ref{sa-lem: gram},
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{X}} \big|\widehat{\mu}_t^{(\boldsymbol{\nu})}(\mathbf{x}) - \mathbb{E} [\widehat{\mu}_t^{(\boldsymbol{\nu})}(\mathbf{x}) | \mathbf{X} ]
- \mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \mathbf{Q}_{t,\mathbf{x}} \big|
&= \sup_{\mathbf{x} \in \mathcal{X}} \big| \mathbf{e}_{1+\boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \big(\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} - \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \big) \mathbf{Q}_{t,\mathbf{x}}\big| \\
& \lesssim_{\mathbb{P}} h^{-|\boldsymbol{\nu}|} \sqrt{\frac{\log(1/h)}{n h^d}} \Big(\sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{1+v}{2+v}}h^d} \Big),
\end{align*}
and
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{X}} \big|\widehat{\mu}_t^{(\boldsymbol{\nu})}(\mathbf{x}) - \mathbb{E}[\widehat{\mu}_t^{(\boldsymbol{\nu})}(\mathbf{x}) | \mathbf{X} ] \big|
\lesssim_{\mathbb{P}} h^{-|\boldsymbol{\nu}|} \Big(\sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{1+v}{2+v}}h^d}\Big),
\end{align*}
which completes the proof.
\qed
\subsection{Proof of Lemma~\ref{sa-lem: covariance}}
Let $\eta_i(\mathbf{x}) = \sum_{t \in \{0,1\}} \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) (\mu_t(\mathbf{X}_i) - \widehat{\boldsymbol{\beta}}_t(\mathbf{x})^\top \mathbf{R}_p(\mathbf{X}_i - \mathbf{x}))$. Then, for all $\mathbf{x}, \mathbf{y} \in \mathcal{B}$, the difference between the estimated and true variance matrices is
\begin{align*}
\widehat{\boldsymbol{\Sigma}}_{t,\mathbf{x},\mathbf{y}} - \boldsymbol{\Sigma}_{t,\mathbf{x},\mathbf{y}} = \mathbf{M}_{1,\mathbf{x},\mathbf{y}} + \mathbf{M}_{2,\mathbf{x},\mathbf{y}} + \mathbf{M}_{3,\mathbf{x},\mathbf{y}} + \mathbf{M}_{4,\mathbf{x},\mathbf{y}}
\end{align*}
where
\begin{align*}
\mathbf{M}_{1,\mathbf{x},\mathbf{y}}
&= \mathbb{E}_n \Big[\mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big)^{\top} \frac{1}{h^d} K\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)K\Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big)\eta_i(\mathbf{x}) \eta_i(\mathbf{y}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big],\\
\mathbf{M}_{2,\mathbf{x},\mathbf{y}}
&= \mathbb{E}_n \Big[\mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) \mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big)^{\top} \frac{1}{h^d} K\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)K\Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big) (\eta_i(\mathbf{x}) + \eta_i(\mathbf{y})) u_i \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big],\\
\mathbf{M}_{3,\mathbf{x},\mathbf{y}}
&= \mathbb{E}_n \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big)^{\top}\frac{1}{h^d}K \Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)K \Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big)(u_i^2 - \sigma_t(\mathbf{X}_i)^2)\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big], \\
\mathbf{M}_{4,\mathbf{x},\mathbf{y}}
&= \mathbb{E}_n \Big[\mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)\mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big)^{\top}\frac{1}{h^d}K \Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)K\Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big)\sigma_t(\mathbf{X}_i)^2\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big]\\
&\qquad - \mathbb{E} \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big)^{\top}\frac{1}{h^d}K\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)K\Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big)\sigma_t(\mathbf{X}_i)^2 \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big].
\end{align*}
For $\mathbf{u}$ and $\mathbf{v}$ multi-indices, let $g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y}) = \frac{1}{h^d}(\frac{\mathbf{X}_i - \mathbf{x}}{h})^{\mathbf{u}}(\frac{\mathbf{X}_i - \mathbf{y}}{h})^{\mathbf{v}} K(\frac{\mathbf{X}_i - \mathbf{x}}{h})K(\frac{\mathbf{X}_i - \mathbf{y}}{h}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)$. Set
\begin{align*}
\mathfrak{R}_n = \sqrt{\frac{\log (1/h)}{n h^d}} + \frac{\log (1/h)}{n^{\frac{1+v}{2+v}}h^d}.
\end{align*}
First, we present a bound on $\max_{1 \leq i \leq n} |\eta_i(\mathbf{x})| \mathds{1}((\mathbf{X}_i - \mathbf{x})/h \in \operatorname{Supp}(K))$. By Lemma~\ref{sa-lem: bias} and Lemma~\ref{sa-lem: Q}, and multi-index $\boldsymbol{\nu}$ such that $|\boldsymbol{\nu}| \leq p$,
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}}|\mathbf{e}_{1 + \boldsymbol{\nu}}^\top \widehat{\mu}_t(\mathbf{x}) - \mathbf{e}_{1 + \boldsymbol{\nu}}^\top \mu_t(\mathbf{x})| \lesssim_{\mathbb{P}} h^{-|\boldsymbol{\nu}|} (h^{p+1} + \mathfrak{R}_n).
\end{align*}
Since $K$ is compactly supported, we have
\begin{align*}
\max_{1 \leq i \leq n} \Big|\sum_{t \in \{0,1\}}\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) (\widehat{\boldsymbol{\beta}}_t(\mathbf{x}) - \boldsymbol{\beta}_t(\mathbf{x}))^\top \mathbf{R}_p(\mathbf{X}_i - \mathbf{x}) \mathds{1}((\mathbf{X}_i - \mathbf{x})/h \in \operatorname{Supp}(K))\Big|
\lesssim_{\mathbb{P}} h^{p+1} + \mathfrak{R}_n.
\end{align*}
Since $\mu_t$ is $p+1$ times continuously differentiable,
\begin{align*}
\max_{1 \leq i \leq n} \Big|\sum_{t \in \{0,1\}}\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) (\mu_t(\mathbf{X}_i)
- \boldsymbol{\beta}_t(\mathbf{x})^\top \mathbf{R}_p(\mathbf{X}_i - \mathbf{x})) \mathds{1}((\mathbf{X}_i - \mathbf{x})/h \in \operatorname{Supp}(K))\Big|
\lesssim h^{p+1}.
\end{align*}
It follows that
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}}\max_{1 \leq i \leq n} |\eta_i(\mathbf{x})| \mathds{1}((\mathbf{X}_i - \mathbf{x})/h \in \operatorname{Supp}(K))
\lesssim_{\mathbb{P}} h^{p+1} + \mathfrak{R}_n.
\end{align*}
\medskip\textit{Term $\mathbf{M}_{1,\mathbf{x},\mathbf{y}}$}. From the proof for Lemma~\ref{sa-lem: gram}, $\sup_{\mathbf{x},\mathbf{y} \in \mathcal{X}} \big|\mathbb{E}_n[ g_n(\mathbf{X}_i; \mathbf{x},\mathbf{y})] - \mathbb{E}[g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y})]\big|\lesssim_{\mathbb{P}} \sqrt{\frac{\log (1/h)}{n h^d}}$. Moreover, $\sup_{\mathbf{x},\mathbf{y} \in \mathcal{X}}\big|\mathbb{E}[g_n(\mathbf{X}_i; \mathbf{x},\mathbf{y})]\big| \lesssim_{\mathbb{P}} 1$. Hence $\sup_{\mathbf{x},\mathbf{y} \in \mathcal{X}} \big|\mathbb{E}_n[g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y})]\big| \lesssim_{\mathbb{P}} 1$. Thus,
\begin{align*}
\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \big|\mathbb{E}_n[g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y}) \eta_i(\mathbf{x}) \eta_i(\mathbf{y})] \big|
& \leq \sup_{\mathbf{x} \in \mathcal{X}}\max_{1 \leq i \leq n} |\eta_i(\mathbf{x})| \mathds{1}((\mathbf{X}_i - \mathbf{x})/h \in \operatorname{Supp}(K))
\cdot \sup_{\mathbf{x}, \mathbf{y} \in \mathcal{X}}\big|\mathbb{E}_n[g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y})]\big| \\
& \lesssim_{\mathbb{P}} (h^{p+1} + \mathfrak{R}_n)^2,
\end{align*}
where we have used Theorem~\ref{sa-thm: Convergence Rates}, which does not depend on this lemma, for $\sup_{\mathbf{x} \in \mathcal{B}} \big| \widehat{\mu}_t(\mathbf{x}) - \mu_t(\mathbf{x})\big| \lesssim_{\mathbb{P}} h^{p+1} + \mathfrak{R}_n$. Finite dimensionality of $\mathbf{M}_{1,\mathbf{x},\mathbf{y}}$ then implies
\begin{align*}
\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}}\left\lVert\mathbf{M}_{1,\mathbf{x},\mathbf{y}}\right\rVert \lesssim_{\mathbb{P}} (h^{p+1} + \mathfrak{R}_n)^2.
\end{align*}
\medskip\textit{Term $\mathbf{M}_{2,\mathbf{x},\mathbf{y}}$}. From the proof of Lemma~\ref{sa-lem: Q}, $\sup_{\mathbf{x},\mathbf{y} \in \mathcal{X}} \big|\mathbb{E}_n[g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y})u_i] - \mathbb{E}[g_n(\mathbf{X}_i; \mathbf{x},\mathbf{y}) u_i]\big| \lesssim_{\mathbb{P}} \mathfrak{R}_n$. Moreover, $\sup_{\mathbf{x},\mathbf{y} \in \mathcal{X}}\big|\mathbb{E}[g_n(\mathbf{X}_i; \mathbf{x},\mathbf{y})u_i]big\| \lesssim_{\mathbb{P}} 1$. Hence, $\sup_{\mathbf{x},\mathbf{y} \in \mathcal{X}} \big|\mathbb{E}_n[g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y})u_i]\big| \lesssim_{\mathbb{P}} 1$. Thus,
\begin{align*}
\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \big|\mathbb{E}_n[g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y})(\eta_i(\mathbf{x}) + \eta_i(\mathbf{y})) u_i] \big|
& \leq \sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \big|\widehat{\mu}_t(\mathbf{x}) - \mu_t(\mathbf{x})\big| \sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \mathbb{E}_n[|g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y}) u_i|]
\lesssim_{\mathbb{P}} h^{p+1} + \mathfrak{R}_n,
\end{align*}
which implies that
\begin{align*}
\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}}\left\lVert\mathbf{M}_{2,\mathbf{x},\mathbf{y}}\right\rVert \lesssim_{\mathbb{P}} h^{p+1} + \mathfrak{R}_n.
\end{align*}
\medskip\textit{Term $\mathbf{M}_{3,\mathbf{x},\mathbf{y}}$}. Define $l_n(\cdot, \cdot;\mathbf{x},\mathbf{y}): \mathcal{X} \times \mathbb{R} \to \mathbb{R}$ as
\begin{align*}
l_n(\xi, \varepsilon;\mathbf{x}, \mathbf{y})
= \frac{1}{h^d} \Big(\frac{\xi - \mathbf{x}}{h}\Big)^{\mathbf{u}} \Big(\frac{\xi - \mathbf{y}}{h}\Big)^{\mathbf{v}} K\Big(\frac{\xi - \mathbf{x}}{h}\Big) K\Big(\frac{\xi - \mathbf{y}}{h}\Big) \mathds{1}(\xi \in \mathcal{A}_t) (\varepsilon^2 - \sigma_t^2(\xi)),
\end{align*}
and consider the function class $\mathcal{L} = \{l_n(\cdot,\cdot;\mathbf{x},\mathbf{y}): \mathbf{x},\mathbf{y} \in \mathcal{X}\}$. Let $L: \mathcal{X} \times \mathbb{R} \to \mathbb{R}$ be $L(\xi, \varepsilon) = \frac{c}{h^d}|\varepsilon^2 - \sigma_t^2(\xi)|$ with $c = \sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \big|\big(\frac{\xi - \mathbf{x}}{h}\big)^{\mathbf{u}} \big(\frac{\xi - \mathbf{y}}{h}\big)^{\mathbf{v}}K\big(\frac{\xi - \mathbf{x}}{h}\big) K\big(\frac{\xi - \mathbf{y}}{h}\big)\big|$. By similar argument as in the proof for Lemma~\ref{sa-lem: Q}, we can show $\mathcal{L}$ is a VC-type class such that $\mathbb{E} [l_n(\mathbf{X}_i,u_i; \mathbf{x}, \mathbf{y})] = 0$, for all $\mathbf{x},\mathbf{y} \in \mathcal{X}$,
\begin{align*}
\sup_{\mathbf{x},\mathbf{y} \in \mathcal{X}}\mathbb{E}[l_n(\mathbf{X}_i,\varepsilon;\mathbf{x},\mathbf{y})^2]^{\frac{1}{2}}
\lesssim \sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \mathbb{E}[g_n(\mathbf{X}_i, u_i; \mathbf{x}, \mathbf{y})^2]^{\frac{1}{2}} \sup_{\xi \in \mathcal{X}} \mathbb{V}[u_i^2|\mathbf{X}_i = \xi]
\lesssim h^{-d/2}
\end{align*}
and
\begin{align*}
\mathbb{E}\big[\max_{1 \leq i \leq n} L(\mathbf{X}_i,u_i)^2\big]^{\frac{1}{2}}
\lesssim h^{-d} \mathbb{E}\big[\max_{1 \leq i \leq n}u_i^4\big]^{1/2}
\lesssim h^{-d} \mathbb{E}\big[\max_{1 \leq i \leq n} u_i^{2+v}\big]^{\frac{2}{2+v}}
\lesssim h^{-d}n^{\frac{2}{2+v}}.
\end{align*}
Applying Corollary 5.1 in \cite{Chernozhukov-Chetverikov-Kato_2014b_AoS}, we obtain
\begin{align*}
\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \big|\mathbb{E}_n[l_n(\mathbf{X}_i,u_i;\mathbf{x},\mathbf{y})] \big|
\lesssim_{\mathbb{P}} \sqrt{\frac{\log (1/h)}{n h^d}} + \frac{\log (1/h)}{n^{\frac{v}{2+v}}h^d}
\end{align*}
and
\begin{align*}
\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \left\lVert\mathbf{M}_{3,\mathbf{x},\mathbf{y}}\right\rVert
\lesssim_{\mathbb{P}} \sqrt{\frac{\log (1/h)}{n h^d}} + \frac{\log (1/h)}{n^{\frac{v}{2+v}}h^d}.
\end{align*}
\medskip\textit{Term $\mathbf{M}_{4,\mathbf{x},\mathbf{y}}$}. Notice that $\{g_n(\cdot; \mathbf{x}, \mathbf{y})\sigma_t^2(\cdot): \mathbf{x}, \mathbf{y} \in \mathcal{B}\}$ is a VC-type of class with constant envelope function $C h^{-d}$ for some positive constant $C$, where $\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \sup_{\xi \in \mathcal{X}}|g_n(\xi; \mathbf{x}, \mathbf{y})\sigma^2(\xi)| \lesssim h^{-d}$ and $\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \mathbb{E}[g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y})^2 \sigma_t(\mathbf{X}_i)^2]^{\frac{1}{2}} \lesssim h^{-d/2}$. Then, similar to the proof of $\mathbf{M}_{1,\mathbf{x},\mathbf{y}}$, we conclude that
\begin{align*}
\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \big|\mathbb{E}_n[g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y})] -\mathbb{E}[g_n(\mathbf{X}_i; \mathbf{x}, \mathbf{y})] \big|
\lesssim \sqrt{\frac{\log (1/h)}{n h^d}}
\end{align*}
and
\begin{align*}
\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \left\lVert\mathbf{M}_{4,\mathbf{x},\mathbf{y}}\right\rVert
\lesssim \sqrt{\frac{\log (1/h)}{n h^d}}.
\end{align*}
\medskip\textit{Final result}. Combining the the upper bounds of the four terms,
\begin{align*}
\sup_{\mathbf{x},\mathbf{y} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Sigma}}_{1, \mathbf{x},\mathbf{y}} - \boldsymbol{\Sigma}_{1, \mathbf{x},\mathbf{y}} \big\|
\lesssim_{\mathbb{P}} h^{p+1} + \sqrt{\frac{\log (1/h)}{n h^d}} + \frac{\log (1/h)}{n^{\frac{v}{2+v}}h^d},
\end{align*}
which implies $\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \|\widehat{\boldsymbol{\Sigma}}_{1, \mathbf{x},\mathbf{y}}\| \lesssim_{\mathbb{P}} 1$. It follows that
\begin{align*}
\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} |\widehat{\Omega}_{1, \mathbf{x},\mathbf{y}}^{(\boldsymbol{\nu})} - \Omega_{1, \mathbf{x},\mathbf{y}}^{(\boldsymbol{\nu})}|
& \leq \frac{1}{n h^{d + 2|\boldsymbol{\nu}|}} \Big( \sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \big\|\widehat{\boldsymbol{\Gamma}}_{1, \mathbf{x}}^{-1} - \boldsymbol{\Gamma}_{1, \mathbf{x}}^{-1} \big\| \big\|\widehat{\boldsymbol{\Sigma}}_{1, \mathbf{x},\mathbf{y}}\big\| \big\|\widehat{\boldsymbol{\Gamma}}_{1, \mathbf{y}}^{-1}\big\|\\
&\qquad + \sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \big\|\boldsymbol{\Gamma}_{1, \mathbf{x}}^{-1}\big\| \big\|\widehat{\boldsymbol{\Sigma}}_{1, \mathbf{x},\mathbf{y}} - \boldsymbol{\Sigma}_{1,\mathbf{x},\mathbf{y}}\big\| \big\|\widehat{\boldsymbol{\Gamma}}_{1, \mathbf{y}}^{-1}\big\|\\
&\qquad + \sup_{\mathbf{x}, \mathbf{y} \in \mathcal{B}} \big\|\boldsymbol{\Gamma}_{1, \mathbf{x}}^{-1}\big\| \big\|\boldsymbol{\Sigma}_{1, \mathbf{x},\mathbf{y}}\big\| \big\|\widehat{\boldsymbol{\Gamma}}_{1, \mathbf{y}}^{-1} - \boldsymbol{\Gamma}_{1, \mathbf{y}}^{-1}\big\| \Big) \\
& \leq \frac{1}{n h^{d + 2|\boldsymbol{\nu}|}} \Big(h^{p+1} + \sqrt{\frac{\log (1/h)}{n h^d}} + \frac{\log (1/h)}{n^{\frac{v}{2+v}}h^d}\Big).
\end{align*}
By Assumption~\ref{sa-assump: DGP}(iv) and Assumption~\ref{sa-assump: Kernel and Boundary}(ii), $\inf_{\mathbf{x} \in \mathcal{B}}\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})} \gtrsim_{\mathbb{P}} (n h^{d + 2 |\boldsymbol{\nu}|})^{-1}$. Therefore, $\inf_{\mathbf{x} \in \mathcal{B}}\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})} \gtrsim (n h^{d + 2 |\boldsymbol{\nu}|})^{-1}$. Furthermore,
\begin{align*}
& \sup_{\mathbf{x} \in \mathcal{B}} \Big|\sqrt{\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}} - \sqrt{\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}}\Big|
\lesssim_{\mathbb{P}} \sup_{\mathbf{x} \in \mathcal{B}}\sqrt{n h^{d + 2|\boldsymbol{\nu}|}} \Big|\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})} - \Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}\Big|
\lesssim_{\mathbb{P}} \frac{1}{\sqrt{n h^{d + 2|\boldsymbol{\nu}|}}} \Big(h^{p+1} + \sqrt{\frac{\log (1/h)}{n h^d}} + \frac{\log (1/h)}{n^{\frac{v}{2+v}}h^d}\Big)
\end{align*}
and
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \left|\frac{h^{-|\boldsymbol{\nu}|}}{\sqrt{\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}}} - \frac{h^{-|\boldsymbol{\nu}|}}{\sqrt{\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}}}\right|
= h^{-|\boldsymbol{\nu}|}\sup_{\mathbf{x} \in \mathcal{B}} \left|\frac{\sqrt{\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}} - \sqrt{\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}}}{\sqrt{\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})} \Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}}} \right|
\lesssim_{\mathbb{P}} \sqrt{n h^d} \Big(h^{p+1} + \sqrt{\frac{\log (1/h)}{n h^d}} + \frac{\log (1/h)}{n^{\frac{v}{2+v}}h^d}\Big),
\end{align*}
which completes the proof.
\qed
\subsection{Proof of Lemma~\ref{sa-lem: bias}}
Define
\begin{align*}
\boldsymbol{\chi}_{t, \mathbf{x}}
= \mathbb{E}_n \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{x})\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\mathfrak{r}_t(\mathbf{X}_i;\mathbf{x}) \Big],
\qquad
\mathfrak{r}_t(\xi;\mathbf{x}) = \mu_t(\xi) - \sum_{0 \leq |\boldsymbol{\omega}| \leq p}\frac{\mu^{(\boldsymbol{\omega})}_t(\mathbf{x})}{\boldsymbol{\omega}!}(\xi - \mathbf{x})^{\boldsymbol{\omega}}.
\end{align*}
Since $\mu_t$ is $(p+1)$-times continuously differentiable, there exists $\boldsymbol{\alpha}_{\mathbf{x},\mathbf{X}_i,t} \in \mathbb{R}^{p+1}$ such that
\begin{align*}
\left\lVert\boldsymbol{\chi}_{t,\mathbf{x}}\right\rVert^2
& = \Big\lVert \frac{1}{n} \sum_{i=1}^n \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{x})\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)^{\top} (\mathbf{0}^{\top}, \boldsymbol{\alpha}_{\mathbf{x},\mathbf{X}_i,t}^{\top})^{\top} \Big \rVert^2 h^{2(p+1)} \\
& \leq \Big(\mathbb{E}_n \Big[ \Big \lVert \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{x})\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)^{\top} \Big \rVert^2 \Big] \Big) \Big( \mathbb{E}_n \big[\|\boldsymbol{\alpha}_{\mathbf{x},\mathbf{X}_i,t}\|^2\big] \Big) h^{2(p+1)},
\end{align*}
where $\sup_{\mathbf{x} \in \mathcal{B}} \max_{t \in \{0,1\}} \max_{1 \leq i \leq n}\left\lVert\boldsymbol{\alpha}_{\mathbf{x},\mathbf{X}_i,t}\right\rVert \lesssim 1$. Since $\frac{\log(1/h)}{n h^d} = o(1)$, the same argument as the proof of Lemma~\ref{sa-lem: gram} shows
\begin{align*}
\mathbb{E}_n \Big[\Big \lVert \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{x})\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)^{\top} \Big \rVert^2\Big] \lesssim_{\mathbb{P}} 1.
\end{align*}
It then follows from Lemma~\ref{sa-lem: gram} that
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big|\mathbb{E}[\widehat{\mu}_t^{(\boldsymbol{\nu})}(\mathbf{x})|\mathbf{X}] - \mu_t^{(\boldsymbol{\nu})}(\mathbf{x})\big|
= \sup_{\mathbf{x} \in \mathcal{B}} \big|\mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} \boldsymbol{\chi}_{t,\mathbf{x}}\big|
\lesssim_{\mathbb{P}} h^{p + 1 - |\boldsymbol{\nu}|}.
\end{align*}
Now also assume that $h = o(1)$. Then, for all $\mathbf{x} \in \mathcal{B}$ and $\xi \in \mathcal{X}$,
\begin{align*}
\mathds{1} \left(K_h(\xi - \mathbf{x}) \neq 0 \right)\left|\gamma_{\mathbf{v}}(\xi; \mathbf{x}) - \frac{|\mathbf{v}|}{\mathbf{v}!} \partial^{\mathbf{v}} \mu_t(\mathbf{x})\right| \leq \frac{|\mathbf{v}|}{\mathbf{v}!} \sup_{\left\lVert\mathbf{u} -\mathbf{u}'\right\rVert \leq h} \left| \partial^{\mathbf{v}}\mu_t(\mathbf{u}) - \partial^{\mathbf{v}}\mu_t(\mathbf{u}')\right| = M_n,
\end{align*}
where $\gamma_{\mathbf{v}}(\xi; \mathbf{x}) = \frac{|\mathbf{v}|}{\mathbf{v}!} \int_{0}^1 (1 - t)^{|\mathbf{v}|-1} \partial_{\mathbf{v}} \mu_t(\mathbf{x} + t(\xi - \mathbf{x})) dt$. By Assumption~\ref{sa-assump: DGP}(iii), $\partial^{\mathbf{v}} \mu_t$ is uniformly continuous on the compact set $\mathcal{X}$. This implies that when $h = o(1)$, $M_n = o(1)$. Letting
\begin{align*}
\widetilde{\boldsymbol{\chi}}_{t,\mathbf{x}}
= \mathbb{E}_n \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{x})\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) ( \sum_{|\mathbf{v}| = p+1} \frac{|\mathbf{v}|}{\mathbf{v}!} \partial^{\mathbf{v}}\mu_t(\mathbf{x})(\mathbf{X}_i - \mathbf{x})^{\mathbf{v}})\Big],
\end{align*}
we conclude that
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big\| \boldsymbol{\chi}_{t,\mathbf{x}} - \widetilde{\boldsymbol{\chi}}_{t,\mathbf{x}} \big\|
\lesssim M_n \sup_{\mathbf{x} \in \mathcal{B}} \Big\| \mathbb{E}_n \Big[\mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{x}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) \Big(\sum_{|\mathbf{v}| = p+1} \frac{|\mathbf{v}|}{\mathbf{v}!} \left|\mathbf{X}_i - \mathbf{x}\right|^{\mathbf{v}} \Big) \Big] \Big\|
= o_{\mathbb{P}}(h^{p+1}),
\end{align*}
where the last equality employs the same arguments as in the proof of Lemma~\ref{sa-lem: gram}. Hence,
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \Big|\mathbb{E}[\widehat{\mu}_t^{(\boldsymbol{\nu})}(\mathbf{x})|\mathbf{X}] - \mu_t^{(\boldsymbol{\nu})}(\mathbf{x}) - h^{p+1-|\boldsymbol{\nu}|} \widehat{B}_{t,\mathbf{x}}^{(\boldsymbol{\nu})}\Big|
= \sup_{\mathbf{x} \in \mathcal{B}} \Big|\mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} \boldsymbol{\chi}_{t,\mathbf{x}} - \mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} \widetilde{\boldsymbol{\chi}}_{t,\mathbf{x}}\Big|
= o_{\mathbb{P}}(h^{p+1-|\boldsymbol{\nu}|}).
\end{align*}
Using Lemma~\ref{sa-lem: gram} and the maximal inequality as in the proof of Lemma~\ref{sa-lem: gram}, we conclude that
\begin{align*}
\max_{t \in \{0,1\}} \sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{B}_{t,\mathbf{x}}^{(\boldsymbol{\nu})} - B_{t,\mathbf{x}}^{(\boldsymbol{\nu})}\big|
\lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{n h^d}}.
\end{align*}
Since $\max_{t \in \{0,1\}}\sup_{\mathbf{x} \in \mathcal{B}}|B_{t,\mathbf{x}}^{(\boldsymbol{\nu})}| \lesssim 1$, it follows that $\max_{t \in \{0,1\}} \sup_{\mathbf{x} \in \mathcal{B}}|\widehat{B}_{t,\mathbf{x}}^{(\boldsymbol{\nu})}| \lesssim_{\mathbb{P}} 1$.
\qed
\subsection{Proof of Theorem~\ref{sa-thm: Convergence Rates}}
The results follow from Lemma~\ref{sa-lem: bias} and Lemma~\ref{sa-lem: Q}.
\qed
\subsection{Proof of Theorem~\ref{sa-thm: MSE}}
For the conditional bias, by Lemma~\ref{sa-lem: bias},
\begin{align*}
& \sup_{\mathbf{x} \in \mathcal{B}} \big|\mathbb{E} \big[\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) - \tau^{(\boldsymbol{\nu})}(\mathbf{x}) \big| \mathbf{X} \big]^2 - (h^{p+1-|\boldsymbol{\nu}|}B_{\mathbf{x}}^{(\boldsymbol{\nu})})^2 \big| \\
& \leq \sup_{\mathbf{x} \in \mathcal{B}} \big|\mathbb{E} \big[\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) - \tau^{(\boldsymbol{\nu})}(\mathbf{x}) \big| \mathbf{X} \big] - h^{p+1-|\boldsymbol{\nu}|}B_{\mathbf{x}}^{(\boldsymbol{\nu})} \big| \cdot \sup_{\mathbf{x} \in \mathcal{B}} \big|\mathbb{E} \big[\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) - \tau^{(\boldsymbol{\nu})}(\mathbf{x}) \big| \mathbf{X} \big] + h^{p+1-|\boldsymbol{\nu}|}B_{\mathbf{x}}^{(\boldsymbol{\nu})} \big| \\
& = o_{\mathbb{P}}(h^{p+1-|\boldsymbol{\nu}|}).
\end{align*}
Since $\sup_{\mathbf{x} \in \mathcal{B}}|B_{t,\mathbf{x}}^{(\boldsymbol{\nu})} - B_{t,\mathbf{x}}^{(\boldsymbol{\nu})}| \lesssim_{\mathbb{P}} \sqrt{\frac{\log(1/h)}{n h^d}}$ from Lemma~\ref{sa-lem: bias},
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big|\mathbb{E} \big[\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) - \tau^{(\boldsymbol{\nu})}(\mathbf{x}) \big| \mathbf{X} \big]^2 - (h^{p+1-|\boldsymbol{\nu}|}B_{\mathbf{x}}^{(\boldsymbol{\nu})})^2 \big| = o_{\mathbb{P}}(h^{p+1-|\boldsymbol{\nu}|}).
\end{align*}
For the conditional variance, by Lemma~\ref{sa-lem: covariance},
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big|\mathbb{V}\big[\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) \big|\mathbf{X} \big] - (n h^{d+2|\boldsymbol{\nu}|})^{-1} V_{\mathbf{x}}^{(\boldsymbol{\nu})} \big|
= o_{\mathbb{P}}((n h^{d+2|\boldsymbol{\nu}|})^{-1}).
\end{align*}
The pointwise MSE expansion follows directly. For the IMSE expansion, notice that
\begin{align*}
& \Big|\operatorname{IMSE}_{\boldsymbol{\nu}} - \int_{\mathcal{B}} \big[(h^{p+1-|\boldsymbol{\nu}|} B_{\mathbf{x}}^{(\boldsymbol{\nu})})^2 + (n h^{d+2|\boldsymbol{\nu}|})^{-1} V_{\mathbf{x}}^{(\boldsymbol{\nu})}\big] w(\mathbf{x}) d\mathfrak{H}^{d-1}(\mathbf{x})\Big|\\
& \leq \int_{\mathcal{B}} |w(\mathbf{x})| d\mathfrak{H}^{d-1}(\mathbf{x}) \cdot \sup_{\mathbf{x} \in \mathcal{B}} \big|\operatorname{MSE}_{\boldsymbol{\nu}}(\mathbf{x}) - (h^{p+1-|\boldsymbol{\nu}|} B_{\mathbf{x}}^{(\boldsymbol{\nu})})^2 - (n h^{d+2|\boldsymbol{\nu}|})^{-1} V_{\mathbf{x}}^{(\boldsymbol{\nu})}\big|\\
& = o_{\mathbb{P}} \big(h^{2p+2-2|\boldsymbol{\nu}|} + (n h^{d+2|\boldsymbol{\nu}|})^{-1} \big),
\end{align*}
which completes the proof.
\qed
\subsection{Proof of Theorem~\ref{sa-thm: Confidence Intervals}}
We have $\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) = \sum_{i = 1}^n Z_i$ with
\begin{align*}
Z_i = \sum_{t \in \{0,1\}} n^{-1} (\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2} \mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{x}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) u_i,
\end{align*}
where $\mathbb{E}[Z_i] = 0$ and $\mathbb{V}[Z_i] = n^{-1}$. By the Berry-Essen Theorem,
\begin{align*}
\sup_{u \in \mathbb{R}} \Big|\mathbb{P} \big(\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) \leq u\big) - \Phi(u) \Big|
\lesssim B_n^{-1} \sum_{i = 1}^n \mathbb{E}[|Z_i|^3],
\end{align*}
where $B_n = \sum_{i = 1}^n \mathbb{V}[Z_i] = 1$. Moreover,
\begin{align*}
\sum_{i = 1}^n \mathbb{E}[|Z_i|^3]
& = n^{-3} (\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-3/2} \sum_{i = 1}^n \mathbb{E} \Big[\Big|\sum_{t \in \{0,1\}} \mathbf{e}_{1+\boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \mathbf{r}_p \left(\frac{\mathbf{X}_i - \mathbf{x}}{h}\right) K_h \left(\mathbf{X}_i - \mathbf{x}\right) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) u_i \Big|^3 \Big]\\
& \lesssim n^{-3} (\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-3/2} \sum_{i = 1}^n \mathbb{E} \Big[\Big|\sum_{t \in \{0,1\}}\mathbf{e}_{1+\boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \mathbf{r}_p \left(\frac{\mathbf{X}_i - \mathbf{x}}{h}\right) K_h \left(\mathbf{X}_i - \mathbf{x}\right) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big|^3 \Big]\\
& \lesssim n^{-2} h^{-|\boldsymbol{\nu}|-d} (\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-3/2} \mathbb{E} \Big[\Big|\sum_{t \in \{0,1\}}\mathbf{e}_{1+\boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \mathbf{r}_p \left(\frac{\mathbf{X}_i - \mathbf{x}}{h}\right) K_h \left(\mathbf{X}_i - \mathbf{x}\right) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big|^2 \Big]\\
& \lesssim n^{-1} h^{-|\boldsymbol{\nu}|-d} (\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2}\\
& \lesssim (n h^d)^{-1/2},
\end{align*}
where the second line uses Assumption~\ref{sa-assump: DGP}(v), the third line uses
\begin{align*}
\bigg|\sum_{t \in \{0,1\}}\mathbf{e}_{1+\boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \mathbf{r}_p \left(\frac{\mathbf{X}_i - \mathbf{x}}{h}\right) K_h \left(\mathbf{X}_i - \mathbf{x}\right) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\bigg| \lesssim h^{-|\boldsymbol{\nu}| - d}
\end{align*}
and the fourth line uses the definition of $\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}$.
Finally, although Lemma~\ref{sa-lem: gram} through Lemma~\ref{sa-lem: covariance} provide convergence results uniformly in $\mathbf{x}$, for pointwise results with fix $\mathbf{x} \in \mathcal{B}$, we can replace the class of functions in those proofs by one containing a \emph{singleton} (corresponding to the evaluation point $\mathbf{x}$). Thus, we obtain the following result:
\begin{align}\label{sa-eq: lin error 2d}
\Big|\widehat{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) - \overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\Big|
\lesssim_{\mathbb{P}} h^{p+1}\sqrt{n h^d} + 1/\sqrt{n h^d} + 1/(n^{\frac{v}{2+v}}h^d),
\end{align}
provided that $h^{p+1}\sqrt{n h^d} \to 0$ and $n^{\frac{v}{2+v}} h^d \to 0$.
The final results follow by weak convergence to a Gaussian distribution, and properties of the distribution function.
\qed
\subsection{Proof of Theorem~\ref{sa-thm: stochastic linearization}}
For all $\mathbf{x} \in \mathcal{B}$, we have $\widehat{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) = \overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) + G_1^{(\boldsymbol{\nu})}(\mathbf{x}) + G_2^{(\boldsymbol{\nu})}(\mathbf{x})$, where
\begin{align*}
G_1^{(\boldsymbol{\nu})}(\mathbf{x})
= \Big(\mathbb{E}\big[\widehat{\tau}^{(\boldsymbol{\nu})}(\mathbf{x}) \big| \mathbf{X} \big] - \tau^{(\boldsymbol{\nu})}(\mathbf{x})\Big) (\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2},
\end{align*}
and
\begin{align*}
G_2^{(\boldsymbol{\nu})}(\mathbf{x})
= \mathbf{e}_{1 + \boldsymbol{\nu}}^{\top}\mathbf{H}^{-1} \Big[\big(\widehat{\boldsymbol{\Gamma}}_{1, \mathbf{x}}^{-1} \mathbf{Q}_{1, \mathbf{x}}- \widehat{\boldsymbol{\Gamma}}_{0, \mathbf{x}}^{-1} \mathbf{Q}_{0, \mathbf{x}}\big) (\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-\frac{1}{2}} - \big(\boldsymbol{\Gamma}_{1, \mathbf{x}}^{-1} \mathbf{Q}_{1, \mathbf{x}}- \boldsymbol{\Gamma}_{0, \mathbf{x}}^{-1} \mathbf{Q}_{0, \mathbf{x}}\big)(\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-\frac{1}{2}}\Big].
\end{align*}
By Lemma~\ref{sa-lem: bias} and Lemma~\ref{sa-lem: covariance},
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \big|G_1^{(\boldsymbol{\nu})}(\mathbf{x})\big| \lesssim_{\mathbb{P}} h^{p+1-|\boldsymbol{\nu}|} (n h^{d + 2|\boldsymbol{\nu}|})^{1/2} \lesssim h^{p+1}\sqrt{n h^d}.
\end{align*}
By Lemma~\ref{sa-lem: gram}, Lemma~\ref{sa-lem: Q} and Lemma~\ref{sa-lem: covariance},
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \Big|e_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \big[\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} - \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1}\big] \mathbf{Q}_{t,\mathbf{x}} (\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2}\Big|
\lesssim_{\mathbb{P}} \sqrt{\log(1/h)} \bigg(\sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{1 + v}{2+v}}h^d}\bigg)
\end{align*}
and
\begin{align*}
& \sup_{\mathbf{x} \in \mathcal{B}} \Big|e_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \mathbf{Q}_{t,\mathbf{x}} \Big[(\widehat{\Omega}_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2} - (\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2} \Big] \Big| \\
& \lesssim_{\mathbb{P}} h^{-|\boldsymbol{\nu}|} \cdot \bigg(\sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{1+v}{2+v}}h^d} \bigg) \cdot \sqrt{n h^{d+ 2\boldsymbol{\nu}}} \bigg(\sqrt{\frac{\log(1/h)}{n h^d}} + \frac{\log(1/h)}{n^{\frac{v}{2+v}}h^d} + h^{p+1}\bigg) \\
& \lesssim \frac{\log(1/h)}{\sqrt{n h^d}} + \frac{(\log(1/h))^{3/2}}{n^{\frac{v}{2+v}}h^d}.
\end{align*}
The result now follows from combining the bounds above.
\qed
\subsection{Proof of Theorem~\ref{sa-thm: Gaussian Strong Approximation: Tstat}}
We verify the high-level conditions of Theorem~\ref{sa-lem: sa thm}. We will employ the following technical lemma.
\begin{lem}[VC Class to VC2 Class]\label{sa-lem: VC Class to VC2 Class}
Assume $\mathcal{F}$ is a VC class on a measure space $(\mathcal{X}, \mathcal{B})$: there exists an envelope function $F$ and positive constants $c(\mathcal{F}), d(\mathcal{F})$ such that for all $\varepsilon\in(0,1)$,
\begin{align*}
\sup_{Q}N(\mathcal{F}, \left\lVert\cdot\right\rVert_{Q,1}, \varepsilon \left\lVertF\right\rVert_{Q,1}) \leq c(\mathcal{F})\varepsilon^{-d(\mathcal{F})},
\end{align*}
where the supremum is taken over all finite discrete measures. Then, $\mathcal{F}$ is also VC2 class: for all $\varepsilon\in(0,1)$,
\begin{align*}
\sup_{Q} N(\mathcal{F}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon \left\lVertF\right\rVert_{Q,2}) \leq c(\mathcal{F}) (\varepsilon^2/2)^{-d(\mathcal{F})},
\end{align*}
where the supremum is taken over all finite discrete measures.
\end{lem}
\noindent\textbf{Proof of Lemma~\ref{sa-lem: VC Class to VC2 Class}}. Let $Q$ be a finite discrete probability measure. Let $f,g \in \mathcal{F}$. Then, $\int |f - g|^2 d Q \leq 2 \int |f - g| |F| d Q$. Define another probability measure $\tilde{Q}(c_k) = F(c_k) Q(c_k) / \left\lVertF\right\rVert_{Q,1}$ on the support of $Q$, denoted by $\{c_1, \ldots, c_k, \ldots\}$. Then,
\begin{align*}
\int |f - g|^2 d Q
\leq 2 \left\lVertF\right\rVert_{Q,1} \int |f - g| d \tilde{Q}
\leq 2 \left\lVertF\right\rVert_{Q,1} \left\lVertf - g\right\rVert_{\tilde{Q},1}.
\end{align*}
Hence, if we take an $\varepsilon^2/2$-net in $(\mathcal{F}, \left\lVert\cdot\right\rVert_{\tilde{Q},1})$ with cardinality no greater than $c(\mathcal{F}) \varepsilon^{- d (\mathcal{F})}$, then for any $f \in \mathcal{F}$, there exists a $g \in \mathcal{F}$ such that $\left\lVertf - g\right\rVert_{\tilde{Q},1} \leq \varepsilon^2/2 \left\lVertF\right\rVert_{\tilde{Q},1}$, and hence
\begin{align*}
\left\lVertf - g\right\rVert_{Q,2}^2 \leq 2 \varepsilon^2/2 \left\lVertF\right\rVert_{Q,1} \left\lVertF\right\rVert_{\tilde{Q},1} \leq \varepsilon^2 \left\lVertF\right\rVert_{Q,2}^2,
\end{align*}
which gives the result.
\qed
Without loss of generality, we assume $\mathcal{X} = [0,1]^d$, and $\mathcal{Q}_{\mathcal{F}_t} = \mathbb{P}_X$ is a valid surrogate measure for $\mathbb{P}_X$ with respect to $\mathcal{F}_t$, and $\phi_{\mathcal{F}_t} = \operatorname{Id}$ is a valid normalizing transformation (as in ). This implies the constants $\mathtt{c}_1$ and $\mathtt{c}_2$ from Theorem~\ref{sa-lem: sa thm} are all $1$.
Consider first the class of functions $\mathcal{F}_t = \{\mathscr{K}_t^{(\boldsymbol{\nu})}(\cdot; \mathbf{x}): \mathbf{x} \in \mathcal{B}\}$, for $t \in \{0,1\}$.
\medskip\textit{Envelope Function}. By Lemma~\ref{sa-lem: gram} and Lemma~\ref{sa-lem: covariance} and the fact that $\operatorname{Supp}(K)$ is compact,
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \sup_{\xi \in \mathcal{X}} \big|\mathscr{K}_t^{(\boldsymbol{\nu})}(\xi; \mathbf{x})\big|
\lesssim \frac{1}{\sqrt{n}h^{d + \boldsymbol{\nu}}} \sup_{\mathbf{x} \in \mathcal{B}}\big( \|\boldsymbol{\Gamma}_{1, \mathbf{x}}^{-1}\| + \|\boldsymbol{\Gamma}_{0, \mathbf{x}}^{-1}\|\big) \sup_{\mathbf{x} \in \mathcal{B}} \big|(\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{- \frac{1}{2}} \big|
\lesssim h^{-d/2}.
\end{align*}
Hence, there exists a constant $C_1 > 0$ such that $\mathtt{M}_{\mathcal{F}_t} = C_1 h^{-d/2}$ is a constant envelope function.
\medskip\textit{$L_1$ Bound}. We have $\mathtt{E}_{\mathcal{F}_t} = \sup_{\mathbf{x} \in \mathcal{B}} \mathbb{E}[|\mathscr{K}_t^{(\boldsymbol{\nu})}(\mathbf{X}_i; \mathbf{x})|] \lesssim h^{d/2}$.
\medskip\textit{Uniform Variation}. Case 1: $K$ is Lipschitz. By Assumption \ref{sa-assump: DGP}(iv) and Assumption \ref{sa-assump: Kernel and Boundary},
\begin{align*}
\mathtt{L}_{\mathcal{F}_t}
= \sup_{\mathbf{x} \in \mathcal{B}}\sup_{\xi, \xi' \in \mathcal{X}} \frac{|\mathscr{K}_t^{(\boldsymbol{\nu})}(\xi; \mathbf{x}) - \mathscr{K}_t^{(\boldsymbol{\nu})}(\xi'; \mathbf{x})|}{\lVert \xi - \xi' \rVert_{\infty}}
\lesssim h^{-d/2-1}.
\end{align*}
Each entry of $\boldsymbol{\Gamma}_{t,\mathbf{x}}$ and $\boldsymbol{\Sigma}_{t,\mathbf{x}}$ are of the form $\int (\frac{\xi - \mathbf{x}}{h})^{\mathbf{u} + \mathbf{v}} K_h(\xi - \mathbf{x})\mathds{1}(\xi \in \mathcal{A}_t)f(\xi)d \xi$ and $\int (\frac{\xi - \mathbf{x}}{h})^{\mathbf{u} + \mathbf{v}} K_h(\xi - \mathbf{x})\sigma_t(\xi)^2 \mathds{1}(\xi \in \mathcal{A}_t)d \xi$ for some multi-index $\mathbf{u}$ and $\mathbf{v}$, respectively. Hence, by Assumption~\ref{sa-assump: Kernel and Boundary}, each entry of $\boldsymbol{\Gamma}_{t,\mathbf{x}}$ and $\boldsymbol{\Sigma}_{t,\mathbf{x}}$ are $h^{-1}$-Lipschitz in $\mathbf{x}$. It follows that there exists a constant $C_2$ such that for all $\mathbf{x}, \mathbf{x}' \in \mathcal{B}$,
\begin{align*}
\big\|\boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} - \boldsymbol{\Gamma}_{t,\mathbf{x}'}^{-1}\big\|
\leq \|\boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1}\| \|\boldsymbol{\Gamma}_{t,\mathbf{x}} - \boldsymbol{\Gamma}_{t,\mathbf{x}}\| \|\boldsymbol{\Gamma}_{t,\mathbf{x}'}^{-1}\|
\leq C_2 h^{-1} \left\lVert\mathbf{x} - \mathbf{x}'\right\rVert.
\end{align*}
Also, by definition of $\Omega_{t,\mathbf{x}}$ and Assumption~\ref{sa-assump: Kernel and Boundary}(iv), there exists $C_3$ such that for all $\mathbf{x},\mathbf{x}' \in \mathcal{X}$,
\begin{align*}
\big|\Omega_{t,\mathbf{x}}^{(\boldsymbol{\nu})} - \Omega_{t,\mathbf{x}'}^{(\boldsymbol{\nu})}\big|
\leq C_3 (n h^{d + 2|\boldsymbol{\nu}| + 1})^{-1} \lVert \mathbf{x} - \mathbf{x}' \rVert_{\infty},
\end{align*}
and
\begin{align*}
\big|(\Omega_{t,\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2} - (\Omega_{t,\mathbf{x}'}^{(\boldsymbol{\nu})})^{-1/2}\big|
\leq \frac{1}{2} \inf_{\mathbf{z} \in \mathcal{X}} (\Omega_{t,\mathbf{z}}^{(\boldsymbol{\nu})})^{-3/2} \big|\Omega_{t,\mathbf{x}}^{(\boldsymbol{\nu})} - \Omega_{t,\mathbf{x}'}^{(\boldsymbol{\nu})}\big|
\leq \frac{1}{2}C_3 h^{-1} (n h^{d + 2 |\boldsymbol{\nu}|})^{1/2} \lVert \mathbf{x} - \mathbf{x}' \rVert_{\infty}.
\end{align*}
It then follows that we have a uniform Lipschitz property with respect to the point of evaluation:
\begin{align*}
\mathtt{l}_{\mathcal{F}_t} = \sup_{\xi \in \mathcal{X}}\sup_{\mathbf{x}, \mathbf{x}' \in \mathcal{B}} \frac{\left|\mathscr{K}_t^{(\boldsymbol{\nu})}(\xi; \mathbf{x}) - \mathscr{K}_t^{(\boldsymbol{\nu})}(\xi; \mathbf{x}') \right|}{\lVert \mathbf{x} - \mathbf{x}' \rVert_{\infty}} \lesssim h^{-d/2-1}.
\end{align*}
Let $\mathbf{x} \in \mathcal{B}$. Then, $\mathscr{K}_t^{(\boldsymbol{\nu})}(\cdot; \mathbf{x})$ is supported on $\mathbf{x} + \mathbf{c} [-h, h]^d$. Then,
\begin{align*}
\mathtt{TV}_{\mathcal{F}_t}
\lesssim \mathfrak{m}\big(\mathbf{c} [-h,h]^d \big) \mathtt{L}_{\mathcal{F}_t}
\lesssim h^{d/2-1}
\end{align*}
Case 2: $K = \mathds{1}(\cdot \in [-1,1]^d)$. Consider
\begin{align*}
\tilde{\mathscr{K}}^{(\boldsymbol{\nu})}_t(\mathbf{u};\mathbf{x})
= n^{-1/2} (\Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})})^{-1/2}\mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \boldsymbol{\Gamma}_{t,\mathbf{x}}^{-1} \mathbf{r}_p \Big(\frac{\mathbf{u} - \mathbf{x}}{h}\Big) h^{-d},
\qquad \mathbf{u} \in \mathcal{X},\quad t \in \{0,1\}.
\end{align*}
Then, $\mathscr{K}^{(\boldsymbol{\nu})}(\mathbf{u};\mathbf{x}) = \tilde{\mathscr{K}}^{(\boldsymbol{\nu})}(\mathbf{u};\mathbf{x}) \mathds{1}(\mathbf{u} - \mathbf{x} \in [-1,1]^d)$ for all $\mathbf{u} \in \mathcal{X}$ and $\mathbf{x} \in \mathcal{B}$, and we set $\tilde{\mathcal{F}}_t = \{\tilde{\mathscr{K}}^{(\boldsymbol{\nu})}(\cdot;\mathbf{x}): \mathbf{x} \in \mathcal{B}\}$, $t \in \{0,1\}$. Then, the argument above implies that $\mathtt{TV}_{\tilde{\mathcal{F}}_t}\lesssim \mathfrak{m}\left(\mathbf{c} [-h,h]^d \right) \mathtt{L}_{\mathcal{F}_t} \lesssim h^{d/2-1}$. Next, set $\mathscr{L} = \{\mathds{1}((\cdot - \mathbf{x})/h \in [-1,1]^d): \mathbf{x} \in \mathcal{B} \}$. Then, using a product rule, we have
\begin{align*}
\mathtt{TV}_{\mathcal{F}_t}
\leq \mathtt{TV}_{\tilde{\mathcal{F}}_t} \mathtt{M}_{\mathscr{L}} + \mathtt{M}_{\tilde{\mathcal{F}}_t} \mathtt{TV}_{\mathscr{L}}
\lesssim h^{d/2 - 1} \cdot 1 + h^{-d/2} h^{d-1}
\lesssim h^{d/2 - 1}.
\end{align*}
\medskip\textit{VC-type Class}. Case 1: $K$ is Lipschitz. We apply \citet[Lemma 7]{Cattaneo-Chandak-Jansson-Ma_2024_Bernoulli}. To make the notation consistent, define
\begin{align*}
f_{\mathbf{x}}(\cdot)
= \frac{1}{\sqrt{n \Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}}}\mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \boldsymbol{\Gamma}_t^{-1} \mathbf{r}_p \left(\cdot\right) K\left(\cdot\right),
\qquad \mathbf{x} \in \mathcal{B},
\end{align*}
and $\mathcal{H} = \{g_{\mathbf{x}} \left(\frac{\cdot - \mathbf{x}}{h}\right): \mathbf{x} \in \mathcal{B} \}$. Notice that $f_{\mathbf{x}}(\frac{\cdot - \mathbf{x}}{h}) = h^d \frac{1}{\sqrt{n \Omega_{\mathbf{x},\mathbf{x}}^{(\boldsymbol{\nu})}}}\mathbf{e}_{1 + \boldsymbol{\nu}}^{\top} \mathbf{H}^{-1} \boldsymbol{\Gamma}^{-1} \mathbf{r}_p(\frac{\cdot - \mathbf{x}}{h}) K_h(\cdot - \mathbf{x})$. Then, the following conditions in \citet[Lemma 7]{Cattaneo-Chandak-Jansson-Ma_2024_Bernoulli} hold (for $\mathbf{z},\mathbf{z}',\mathbf{z}'' \in \mathcal{X}$):
\begin{enumerate}[label=\normalfont(\roman*),noitemsep,leftmargin=*]
\item boundedness: $\sup_{\mathbf{z}} \sup_{\mathbf{z}'} \left|f_{\mathbf{z}}(\mathbf{z}') \right| \leq \mathbf{c}$,
\item compact support: $\operatorname{supp}(f_{\mathbf{z}}(\cdot)) \subseteq [-\mathbf{c}, \mathbf{c}]^d$,
\item Lipschitz continuity: $\sup_{\mathbf{z}} \left|f_{\mathbf{z}}(\mathbf{z}') - f_{\mathbf{z}}(\mathbf{z}'') \right| \leq \mathbf{c} |\mathbf{z}' - \mathbf{z}''|$ and $\sup_{\mathbf{z}}|f_{\mathbf{z}'}(\mathbf{z}) - f_{\mathbf{z}''}(\mathbf{z})| \leq \mathbf{c} h^{-1} |\mathbf{z}' - \mathbf{z}''|$,
\end{enumerate}
and therefore there exists a constant $\mathbf{c}'$ only depending on $\mathbf{c}$ and $d$ that for any $0 \leq \varepsilon \leq 1$,
\begin{align*}
\sup_{Q} N\big(\mathcal{H}, \left\lVert\cdot\right\rVert_{Q,1}, (2c+1)^{d+1}\varepsilon \big)
\leq \mathbf{c}' \varepsilon^{-d-2} + 1,
\end{align*}
where the supremum is taken over all finite discrete measures on $\mathcal{X} = [0,1]^d$. It then follows from Lemma~\ref{sa-lem: VC Class to VC2 Class} that with the constant envelope function $\mathtt{M}_{\mathcal{F}_t} = h^{-d/2}$, for any $0 \leq \varepsilon \leq 1$,
\begin{align*}
\sup_{Q} N\big(\mathcal{F}_t, \left\lVert\cdot\right\rVert_{Q,2}, (2c+1)^{d+1}\varepsilon \mathtt{M}_{\mathcal{F}_t}\big) \leq \mathbf{c}' 2^{2d+4}\varepsilon^{-2d-4} + 1,
\end{align*}
where the supremum is taken over all finite discrete measures.
\emph{Case 2: Suppose $K = \mathds{1}(\cdot \in [-1,1]^d)$.} Recall $\tilde{\mathcal{F}}_t$ and $\mathscr{L}$ defined in the analysis of \textit{uniform variation}. The same argument as before shows
\begin{align*}
\sup_{Q} N\big(\tilde{\mathcal{F}}_t, \left\lVert\cdot\right\rVert_{Q,2}, (2c+1)^{d+1}\varepsilon \mathtt{M}_{\tilde{\mathcal{F}}_t}\big)
\leq \mathbf{c}' 2^{2d+4} \varepsilon^{-2d-4} + 1, \qquad \varepsilon \in (0,1],
\end{align*}
where the supremum is taken over all finite discrete measures, and $\tilde{\mathcal{F}}_t = h^{-d/2}$. By \citet[Example 2.6.1]{van-der-Vaart-Wellner_1996_Book}, the class $\mathscr{L} = \{\mathds{1}((\cdot - \mathbf{x})/h \in [-1,1]^d): \mathbf{x} \in \mathcal{B}\}$ has VC dimension no greater than $2d$, and by \citet[Theorem 2.6.4]{van-der-Vaart-Wellner_1996_Book},
\begin{align*}
\sup_{Q}N(\mathscr{L}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon) \leq 2d (4 e)^{2d} \varepsilon^{-4d}, \qquad 0 < \varepsilon \leq 1,
\end{align*}
where the supremum is taken over all finite discrete measures on $\mathcal{X} = [0,1]^d$. Putting together, we have
\begin{align*}
\sup_{Q}N(\mathcal{F}_t, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon C_1 \mathtt{M}_{\tilde{\mathcal{F}}_t}) \leq C_2 \varepsilon^{-4d},
\end{align*}
where $C_1$, $C_2$ are constants only depending on $d$, and the supremum is taken over all finite discrete measures on $\mathcal{X} = [0,1]^d$.
\bigskip
Consider next the class of functions $\mathcal{G} = \{g_{\mathbf{x}}: \mathbf{x} \in \mathcal{B}\}$, where $g_{\mathbf{x}}(\mathbf{u}) = \mathds{1}(\mathbf{u}\in\mathcal{A}_1) \mathscr{K}_1^{(\boldsymbol{\nu})}(\mathbf{u};\mathbf{x}) - \mathds{1}(\mathbf{u}\in\mathcal{A}_0) \mathscr{K}_0^{(\boldsymbol{\nu})}(\mathbf{u};\mathbf{x})$. We have immediately that $\mathtt{M}_{\mathcal{G}} \lesssim h^{-d/2}$, $\mathtt{E}_{\mathcal{G}} \lesssim h^{d/2}$, and
\begin{align*}
\sup_{Q} N(\mathcal{G}, \left\lVert\cdot\right\rVert_{Q,2}, \varepsilon (2 c +1)^{d+1} \mathtt{M}_{\mathcal{G}}) \leq 2 \mathbf{c}' \varepsilon^{-4d-4} + 2,
\end{align*}
where the supremum is taken over all finite discrete measures.
\medskip\textit{Total Variation}. Observe that $\mathds{1}(\mathbf{u}\in\mathcal{A}_t) \mathscr{K}_t^{(\boldsymbol{\nu})}(\mathbf{u};\mathbf{x}) \neq 0$ implies $E_{t,\mathbf{x}} = \mathbf{u} \in \{\mathbf{y} \in \mathcal{A}_t: (\mathbf{y} - \mathbf{x})/h \in \operatorname{Supp}(K)\}$, and $\mathds{1}(\mathbf{u} \in \mathcal{A}_t) \mathscr{K}_t^{(\boldsymbol{\nu})}(\mathbf{u};\mathbf{x}) = \mathds{1}(\mathbf{u} \in E_{t,\mathbf{x}}) \mathscr{K}_t^{(\boldsymbol{\nu})}(\mathbf{u};\mathbf{x})$, for all $\mathbf{u} \in \mathcal{X}$. By the assumption that the De Giorgi perimeter of $E_{t,\mathbf{x}}$ satisfies $\mathscr{L}(E_{t,\mathbf{x}}) \leq C h^{d-1}$ and using $\mathtt{TV}_{\{gf\}} \leq \mathtt{M}_{\{g\}} \mathtt{TV}_{\{f\}} + \mathtt{M}_{\{f\}} \mathtt{TV}_{\{g\}}$ for any two functions $g$ and $f$, we have
\begin{align*}
\mathtt{TV}_{\mathcal{G}} = \sup_{\mathbf{x} \in \mathcal{B}}\mathtt{TV}_{\{g_{\mathbf{x}}\}} & \leq \sup_{\mathbf{x} \in \mathcal{B}} \sum_{t \in \{0,1\}} \mathtt{TV}_{\{\mathds{1}_{\mathcal{A}_t} \mathscr{K}_t^{(\boldsymbol{\nu})}(\cdot;\mathbf{x}) \}}
\leq \sup_{\mathbf{x} \in \mathcal{B}} \sum_{t \in \{0,1\}} \mathtt{TV}_{\{\mathscr{K}_t^{(\boldsymbol{\nu})}(\cdot;\mathbf{x}) \}} + \mathtt{M}_{\mathcal{F}_t} \mathtt{TV}_{\{\mathds{1}_{E_{t,\mathbf{x}}}\}}
\lesssim h^{d/2-1}.
\end{align*}
We completed the verification of all the high-level sufficient conditions of Theorem~\ref{sa-lem: sa thm}, which immediately give the result.
\qed
\subsection{Proof of Theorem~\ref{sa-thm: Confidence Bands}}
The proof is divided in three technical lemmas.
\begin{lem}[KS Distance Between $\overline{\operatorname{T}}^{(\boldsymbol{\nu})}$ and $Z^{(\boldsymbol{\nu})}$]\label{sa-lem: infeasible gaussian to bahadur representation}
Suppose the conditions of Theorem~\ref{sa-thm: Gaussian Strong Approximation: Tstat} hold. Then, for any multi-index $|\boldsymbol{\nu}| \leq p$,
\begin{align*}
\sup_{u \in \mathbb{R}} \Big|\mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u\Big)
- \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u\Big)\Big|
\lesssim \bigg((\log n)^{\frac{3}{2}}\Big(\frac{1}{n h^d}\Big)^{\frac{1}{d+2}\cdot \frac{v}{v+2}}
+ \log(n) \sqrt{\frac{1}{n^{\frac{v}{v+2}}h^d}} \bigg)^{1/2}.
\end{align*}
\end{lem}
\noindent\textbf{Proof of Lemma~\ref{sa-lem: infeasible gaussian to bahadur representation}.} Let $\mathfrak{R}_n = (\log n)^{\frac{3}{2}}(\frac{1}{n h^d})^{\frac{1}{d+2}\cdot \frac{v}{v+2}} + \log(n)\sqrt{\frac{1}{n^{\frac{v}{2+v}}h^d}}$, and $a_n$ positive sequence to be determined below. For any $u > 0$,
\begin{align*}
&\mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u \Big)\\
&\leq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq \sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) - Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| + u \Big) \\
&\leq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u + a_n \Big)
+ \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x}) - \overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| > a_n \Big)\\
&\leq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u \Big)
+ 4 a_n \Big(\mathbb{E} \Big[\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big|\Big] + 1\Big)
+ \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x}) - \overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| > a_n \Big) \\
&\leq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u \Big)
+ 4 a_n \Big(\mathbb{E} \Big[\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big|\Big] + 1\Big) + \frac{C \mathfrak{R}_n}{a_n},
\end{align*}
where in the fourth line we have used the Gaussian Anti-concentration Inequality in \cite[Theorem 2.1]{Chernozhukov-Chetverikov-Kato_2014a_AoS}, and in the last line we have used the tail bound in Theorem~\ref{sa-thm: Gaussian Strong Approximation: Tstat}. Similarly, for any $u > 0$, we have the lower bound
\begin{align*}
&\mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u \Big)\\
&\geq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u - \sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) - Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \Big) \\
&\geq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u - a_n \Big)
- \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x}) - \overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| > a_n \Big)\\
&\geq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u \Big)
- 4 a_n \Big(\mathbb{E} \Big[\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big|\Big] + 1\Big)
- \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x}) - \overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| > a_n \Big) \\
&\geq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u \Big)
- 4 a_n \Big(\mathbb{E} \Big[\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big|\Big] + 1\Big)
- \frac{C \mathfrak{R}_n}{a_n}.
\end{align*}
Notice that $Z^{(\boldsymbol{\nu})}(\mathbf{x}), \mathbf{x} \in \mathcal{B}$ is a mean-zero Gaussian process satisfying
\begin{align*}
\mathbb{E} \big[\big(Z^{(\boldsymbol{\nu})}(\mathbf{x}) - Z^{(\boldsymbol{\nu})}(\mathbf{y})\big)^2\big]^{\frac{1}{2}}
= \mathbb{E} \big[(\mathscr{K}(\mathbf{X}_i, \mathbf{x}) - \mathscr{K}(\mathbf{X}_i, \mathbf{y}))^2 \sigma(\mathbf{X}_i)^2\big]^{\frac{1}{2}}
\leq C' l_{n,2} \lVert \mathbf{x} - \mathbf{y} \rVert_{\infty},
\end{align*}
where $C'$ is a constant and $l_{n,2} \asymp h^{-1}$, and hence
\begin{align*}
\sup_{\mathbf{x} \in \mathcal{B}} \mathbb{E} \big[\big(Z^{(\boldsymbol{\nu})}(\mathbf{x}) - Z^{(\boldsymbol{\nu})}(\mathbf{y})\big)^2\big]^{\frac{1}{2}}
= \sup_{\mathbf{x} \in \mathcal{B}} \mathbb{E} \big[\mathscr{K}(\mathbf{X}_i, \mathbf{x})^2 \sigma^2(\mathbf{X}_i)\big] \lesssim 1.
\end{align*}
Then, by Corollary 2.2.8 in \cite{van-der-Vaart-Wellner_1996_Book}, we have $\mathbb{E} \big[\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big|\big] \lesssim 1$. Choosing $a_n \asymp \sqrt{\mathfrak{R}_n}$, the result now follows.
\qed
\begin{lem}[KS Distance Between $\widehat{\operatorname{T}}^{(\boldsymbol{\nu})}$ and $\overline{\operatorname{T}}^{(\boldsymbol{\nu})}$]\label{sa-lem: t stats to bahadur}
Suppose the conditions in Theorem~\ref{sa-thm: Gaussian Strong Approximation: Tstat} hold. Then, for any multi-index $|\boldsymbol{\nu}| \leq p$,
\begin{align*}
\sup_{u \in \mathbb{R}} \Big| \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u\Big)
- \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u \Big)\Big| = o(1).
\end{align*}
\end{lem}
\noindent\textbf{Proof of Lemma~\ref{sa-lem: t stats to bahadur}.} Let $\mathfrak{R}_n = (\log n)^{\frac{3}{2}}(\frac{1}{n h^d})^{\frac{1}{d+2}\cdot \frac{v}{v+2}} + \log(n)\sqrt{\frac{1}{n^{\frac{v}{2+v}}h^d}}$ and
\begin{align*}
a_n
= o\bigg(\sqrt{\log(1/h)}\Big(\sqrt{n^{-1} h^{-d} \log(1/h)} + \frac{\log(1/h)}{n^{\frac{v}{2+v}}h^d} \Big)
+ h^{p+1} \sqrt{n h^d}\bigg).
\end{align*}
Then, $\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) - \widehat{\operatorname{T}}(\mathbf{x})\big| = o_{\mathbb{P}}(a_n)$. Hence, for any $u > 0$,
\begin{align*}
& \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{\operatorname{T}}(\mathbf{x})\big| \leq u\Big)\\
& \leq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u + a_n\Big)
+ \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) - \widehat{\operatorname{T}}(\mathbf{x})\big| \geq a_n \Big)\\
& \leq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u + a_n\Big) + \sqrt{\mathfrak{R}_n} + o(1)\\
& \leq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u\Big)
+ 4 a_n \Big(\mathbb{E}\Big[\sup_{\mathbf{x} \in \mathcal{B}}\big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big|\Big] + 1 \Big)
+ \sqrt{\mathfrak{R}_n} + o(1) \\
& \leq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u\Big)
+ 4 a_n \Big(\mathbb{E}\Big[\sup_{\mathbf{x} \in \mathcal{B}}\big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big|\Big] + 1 \Big)
+ 2 \sqrt{\mathfrak{R}_n} + o(1),
\end{align*}
where the third line uses Lemma~\ref{sa-lem: infeasible gaussian to bahadur representation} and $\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) - \widehat{\operatorname{T}}(\mathbf{x})\big| = o_{\mathbb{P}}(a_n)$, the fourth line uses \cite[Theorem 2.1]{Chernozhukov-Chetverikov-Kato_2014a_AoS}, and the last line uses Lemma~\ref{sa-lem: infeasible gaussian to bahadur representation} again. Similarly,
\begin{align*}
& \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{\operatorname{T}}(\mathbf{x})\big| \leq u\Big)\\
& \geq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u - a_n\Big)
- \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x}) - \widehat{\operatorname{T}}(\mathbf{x})\big| \geq a_n \Big)\\
& \geq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u - a_n\big) - \sqrt{\mathfrak{R}_n} + o(1)\\
& \geq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u\Big)
- 4 a_n \Big(\mathbb{E}\Big[\sup_{\mathbf{x} \in \mathcal{B}}\big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big|\Big] + 1 \Big)
- \sqrt{\mathfrak{R}_n} + o(1) \\
& \geq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\overline{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u\Big)
- 4 a_n \Big(\mathbb{E}\Big[\sup_{\mathbf{x} \in \mathcal{B}}\big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big|\Big] + 1 \Big)
- 2 \sqrt{\mathfrak{R}_n} + o(1).
\end{align*}
From the proof of Lemma~\ref{sa-lem: infeasible gaussian to bahadur representation}, $\mathbb{E}\left[\sup_{\mathbf{x} \in \mathcal{B}}\left|Z^{(\boldsymbol{\nu})}(\mathbf{x})\right|\right] \lesssim 1$. Hence, the result follows.
\qed
\begin{lem}[KS Distance Between $Z^{(\boldsymbol{\nu})}$ and $\widehat{Z}^{(\boldsymbol{\nu})}$]\label{sa-lem: feasible gaussian to infeasible gaussian}
Suppose the conditions for Theorem~\ref{sa-thm: Gaussian Strong Approximation: Tstat} hold. Then, for any multi-index $|\boldsymbol{\nu}| \leq p$,
\begin{align*}
\sup_{\mathbf{u} \in \mathbb{R}} \Big|\mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq u \Big)
- \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x}) \big| \leq u \big| \mathbf{W} \Big) \Big|
\lesssim_{\mathbb{P}} \log(n) \Big(\sqrt{\frac{\log n}{n h^d}} + \frac{\log n}{n^{\frac{v}{2+v}}h^d} + h^{p+1}\Big)^{1/2}.
\end{align*}
\end{lem}
\noindent \textbf{Proof of Lemma~\ref{sa-lem: feasible gaussian to infeasible gaussian}.} First, using Lemma~\ref{sa-lem: covariance}, we provide an upper bound between covariance functions of the feasible Gaussian process and the infeasible Gaussian process. Letting $\boldsymbol{\Pi}_{\mathbf{x}, \mathbf{y}} = \Omega_{\mathbf{x}, \mathbf{y}} / \sqrt{\Omega_{\mathbf{x},\mathbf{x}} \Omega_{\mathbf{y}}}$ and $\widehat{\boldsymbol{\Pi}}_{\mathbf{x},\mathbf{y}} = \widehat{\Omega}_{\mathbf{x},\mathbf{y}} / \sqrt{\widehat{\Omega}_{\mathbf{x},\mathbf{x}} \widehat{\Omega}_{\mathbf{y}}}$,
\begin{align*}
\sup_{\mathbf{x},\mathbf{y} \in \mathcal{X}} \big|\boldsymbol{\Pi}_{\mathbf{x}, \mathbf{y}} - \widehat{\boldsymbol{\Pi}}_{\mathbf{x},\mathbf{y}}\big|
= \sup_{\mathbf{x}, \mathbf{y} \in \mathcal{X}}
\Big|\frac{\Omega_{\mathbf{x}, \mathbf{y}} - \widehat{\Omega}_{\mathbf{x}, \mathbf{y}}}{\sqrt{\Omega_{\mathbf{x},\mathbf{x}} \Omega_{\mathbf{y}}}}
+ \frac{\widehat{\Omega}_{\mathbf{x}, \mathbf{y}}}{\sqrt{\widehat{\Omega}_{\mathbf{x},\mathbf{x}} \widehat{\Omega}_{\mathbf{y}}}} \Big(\sqrt{\frac{\widehat{\Omega}_{\mathbf{x},\mathbf{x}} \widehat{\Omega}_{\mathbf{y}}}{\Omega_{\mathbf{x},\mathbf{x}} \Omega_{\mathbf{y}}}} - 1\Big)\Big|
\end{align*}
From Lemma~\ref{sa-lem: covariance} and the fact that $|\sqrt{x} - \sqrt{y}| \leq (x \wedge y)^{-1/2} |x - y|/2$ for $x, y > 0$,
\begin{align*}
\sup_{\mathbf{x},\mathbf{y} \in \mathcal{X}}
\frac{\big|\big(\widehat{\Omega}_{\mathbf{x},\mathbf{x}} \widehat{\Omega}_{\mathbf{y}}\big)^{1/2} - \big(\Omega_{\mathbf{x},\mathbf{x}} \Omega_{\mathbf{y}}\big)^{1/2}\big|}
{\big(\Omega_{\mathbf{x},\mathbf{x}} \Omega_{\mathbf{y}}\big)^{1/2}}
\lesssim \frac{\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{X}} \big|\widehat{\Omega}_{\mathbf{x},\mathbf{x}} \widehat{\Omega}_{\mathbf{y}} - \Omega_{\mathbf{x},\mathbf{x}} \Omega_{\mathbf{y}}\big|}
{\inf_{\mathbf{x},\mathbf{y}} \widehat{\Omega}_{\mathbf{x},\mathbf{x}} \widehat{\Omega}_{\mathbf{y}} \wedge \inf_{\mathbf{x},\mathbf{y}} \Omega_{\mathbf{x},\mathbf{x}} \Omega_{\mathbf{y}}}
\lesssim_{\mathbb{P}} h^{p+1} + \sqrt{\frac{\log n}{n h^d}} + \frac{\log n}{n^{\frac{v}{2+v}}h^d}
\end{align*}
and
\begin{align*}
\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{X}} \frac{\big|\Omega_{\mathbf{x},\mathbf{y}} - \widehat{\Omega}_{\mathbf{x},\mathbf{y}}\big|}{\sqrt{\Omega_{\mathbf{x},\mathbf{x}} \Omega_{\mathbf{y}}}}
\lesssim_{\mathbb{P}} h^{p+1} + \sqrt{\frac{\log n}{n h^d}} + \frac{\log n}{n^{\frac{v}{2+v}}h^d}.
\end{align*}
Therefore, letting $\mathfrak{R}_n = \sqrt{\frac{\log n}{n h^d}} + \frac{\log n}{n^{\frac{v}{2+v}}h^d}$, it follows that $\sup_{\mathbf{x}, \mathbf{y} \in \mathcal{X}} |\boldsymbol{\Pi}_{\mathbf{x}, \mathbf{y}} - \widehat{\boldsymbol{\Pi}}_{\mathbf{x},\mathbf{y}}| \lesssim_{\mathbb{P}} h^{p+1} + \mathfrak{R}_n$. Then, we bound the KS distance between the maximum of $Z_n$ and $\widehat{Z}^{(\boldsymbol{\nu})}$ on a $\delta_n$-net of $\mathcal{X}$, denoted by $\mathcal{X}_{\delta_n}$: for all $\mathbf{x} \in \mathcal{B}$, there exists $\mathbf{z} \in \mathcal{X}_{\delta_n}$ such that $\lVert \mathbf{x} - \mathbf{z} \rVert_{\infty} \leq \delta_n$. Since $\mathcal{X}$ is compact, we can assume $M : = \operatorname{Card}\left(\mathcal{X}_{\delta_n}\right) \lesssim \delta_n^{-d}$. Denote $\mathbf{Z}_n^{\delta_n}$ and $\widehat{\mathbf{Z}}_n^{\delta_n}$ to the process $Z_n$ and $\widehat{Z}^{(\boldsymbol{\nu})}$ restricted on $\mathcal{X}_{\delta_n}$, respectively. Then, by \cite[Theorem 2.1]{chernozhuokov2022improved},
\begin{align*}
\sup_{\mathbf{y} \in \mathbb{R}^{M}} \big|\mathbb{P}(\mathbf{Z}_n^{\delta_n}\leq \mathbf{y}) - \mathbb{P}(\widehat{\mathbf{Z}}_n^{\delta_n} \leq \mathbf{y} | \mathbf{W}) \big|
\lesssim \log(M) \sup_{\mathbf{x},\mathbf{y} \in \mathcal{X}} \big|\boldsymbol{\Pi}_{\mathbf{x}, \mathbf{y}} - \widehat{\boldsymbol{\Pi}}_{\mathbf{x},\mathbf{y}}\big|^{\frac{1}{2}}
\lesssim_{\mathbb{P}} \log(M) (\mathfrak{R}_n + h^{p+1})^{\frac{1}{2}},
\end{align*}
and hence
\begin{align*}
\sup_{x \in \mathbb{R}} \big|\mathbb{P}(\lVert \mathbf{Z}_n^{\delta_n} \rVert_{\infty} \leq x ) - \mathbb{P}(\lVert \widehat{\mathbf{Z}}_n^{\delta_n} \rVert_{\infty} \leq x | \mathbf{W})\big|
&\leq \sup_{x \in \mathbb{R}} \big|\mathbb{P}(- x\mathbf{1} \leq \mathbf{Z}_n^{\delta_n} \leq x \mathbf{1}) - \mathbb{P}(- x\mathbf{1} \leq \widehat{\mathbf{Z}}_n^{\delta_n} \leq x \mathbf{1} | \mathbf{W})\big| \\
& \lesssim_{\mathbb{P}} \log(M) (\mathfrak{R}_n + h^{p+1})^{\frac{1}{2}} =: \mathfrak{R}_M.
\end{align*}
Finally, we bound the KS distance on the whole $\mathcal{X}$ with the help of a sequence $a_n > 0$ to be determined. Let
\begin{align*}
\Psi_{\delta_n}(a_n)
= \mathbb{P} \Big(\sup_{\lVert \mathbf{x} - \mathbf{y} \rVert_{\infty} \leq \delta_n} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x}) - Z^{(\boldsymbol{\nu})}(\mathbf{y})\big| \geq a_n \Big)
\end{align*}
and
\begin{align*}
\widehat{\Psi}_{\delta_n}(a_n)
= \mathbb{P} \Big(\sup_{\lVert \mathbf{x} - \mathbf{y} \rVert_{\infty} \leq \delta_n} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x}) - \widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{y})\big| \geq a_n \Big| \mathbf{W} \Big).
\end{align*}
Then, for all $t > 0$,
\begin{align*}
& \mathbb{P} \Big( \sup_{\mathbf{x} \in \mathcal{B}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq t \Big) \\
& \leq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}_{\delta_n}} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x}) \big| \leq t + a_n \Big) + \Psi_{\delta_n}(a_n)\\
& \leq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}_{\delta_n}} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x}) \big| \leq t + a_n \Big| \mathbf{W}\Big)
+ \Psi_{\delta_n}(a_n) + \mathfrak{R}_M\\
& \leq \mathbb{P} \Big(\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x}) \big| \leq t + a_n \Big| \mathbf{W}\Big)
+ \Psi_{\delta_n}(a_n) + \widehat{\Psi}_{\delta_n} (a_n) + \mathfrak{R}_M \\
& \leq \mathbb{P} \Big( \sup_{\mathbf{x} \in \mathcal{B}} \big| \widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq t \big| \mathbf{W}\Big)
+ 4 a_n \Big(\mathbb{E} \Big[\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \Big| \mathbf{W}\Big] + 1\Big)
+ \Psi_{\delta_n}(a_n) + \widehat{\Psi}_{\delta_n}(a_n) + \mathfrak{R}_M.
\end{align*}
Similarly, for all $t > 0$,
\begin{align*}
\mathbb{P} \Big( \sup_{\mathbf{x} \in \mathcal{B}} \Big|Z^{(\boldsymbol{\nu})}(\mathbf{x})\Big| \leq t \Big)
& \geq \mathbb{P} \Big( \sup_{\mathbf{x} \in \mathcal{B}} \big| \widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq t \Big| \mathbf{W} \Big)
- 4 a_n \Big(\mathbb{E} \Big[\sup_{\mathbf{x} \in \mathcal{B}} \Big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x})\Big| \Big| \mathbf{W} \Big] + 1\Big)\\
&\qquad - \Psi_{\delta_n}(a_n) - \widehat{\Psi}_{\delta_n} (a_n) - \mathfrak{R}_M.
\end{align*}
Since $\mathfrak{R}_M$ depends on $\delta_n$ through $\log M \asymp \log (\delta_n^{-d})$, by choosing $\delta_n = n^{-s}$ for large enough $s$, the term $\mathfrak{R}_M$ will dominate the terms $\Psi_{\delta_n}(a_n)$ and $\widehat{\Psi}_{\delta_n}(a_n)$. More precisely, for any $\delta$,
\begin{align*}
& \sup_{\lVert \mathbf{x} - \mathbf{y} \rVert_{\infty} \leq \delta} \mathbb{E}\Big[\big(\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x}) - \widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{y})\big)^2 \Big| \mathbf{W} \Big]\\
&= \sup_{\lVert \mathbf{x} - \mathbf{y} \rVert_{\infty} \leq \delta}
\big(\widehat{\Omega}_{\mathbf{x},\mathbf{x}} \widehat{\Omega}_{\mathbf{y}}\big)^{-\frac{1}{2}} \Big(\frac{1}{n h^d}\Big)^2
\sum_{i=1}^n \widehat{\varepsilon}^2_i \mathds{1}(\mathbf{X}_i \in \mathcal{A}_1)\\
&\qquad\qquad \cdot \Big(\mathbf{e}_1^T \widehat{\boldsymbol{\Gamma}}_{1, \mathbf{x}}^{-1} \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) K\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)
- \mathbf{e}_1^T \widehat{\boldsymbol{\Gamma}}_{1, \mathbf{x}}^{-1} \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big) K\Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big)\Big)^2 \\
&\qquad + \sup_{\lVert \mathbf{x} - \mathbf{y} \rVert_{\infty} \leq \delta}
\big(\widehat{\Omega}_{\mathbf{x},\mathbf{x}} \widehat{\Omega}_{\mathbf{y}} \big)^{-\frac{1}{2}} \Big(\frac{1}{n h^d}\Big)^2
\sum_{i=1}^n \widehat{\varepsilon_i}^2 \mathds{1} \left(\mathbf{X}_i \in \mathcal{A}_0\right)\\
&\qquad\qquad \cdot \Big(\mathbf{e}_1^T \widehat{\boldsymbol{\Gamma}}_{0, \mathbf{x}}^{-1} \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) K\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) - \mathbf{e}_1^T \widehat{\boldsymbol{\Gamma}}_{0, \mathbf{x}}^{-1} \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big) K\Big(\frac{\mathbf{X}_i - \mathbf{y}}{h}\Big)\Big)^2 \\
& \lesssim_{\mathbb{P}} h^{-d-2} \delta^2,
\end{align*}
where the last line uses Lemma~\ref{sa-lem: covariance}, Lemma~\ref{sa-lem: gram}, and the almost sure bound on the Lipschitz constant from the proof of Theorem~\ref{sa-thm: Gaussian Strong Approximation: Tstat}, for some constant $C > 0$. Similarly, for any $\delta > 0$,
\begin{align*}
\sup_{\lVert \mathbf{x} - \mathbf{y} \rVert_{\infty} \leq \delta} \mathbb{E} \Big[\big(Z^{(\boldsymbol{\nu})}(\mathbf{x}) - Z^{(\boldsymbol{\nu})}(\mathbf{y})\big)^2 \Big]
= \sup_{\lVert \mathbf{x} - \mathbf{y} \rVert_{\infty} \leq \delta} \mathbb{E} \Big[\big(\mathscr{K}(\mathbf{X}_i, \mathbf{x}) - \mathscr{K}(\mathbf{X}_i, \mathbf{y})\big)^2 \varepsilon_i^2 \Big]
\leq C' h^{-2} \delta^2,
\end{align*}
Then, by \cite[Corollary 2.2.5]{van-der-Vaart-Wellner_1996_Book},
\begin{align*}
\mathbb{E} \Big[\sup_{\lVert \mathbf{x} - \mathbf{y} \rVert_{\infty}\leq \delta_n} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x}) - \widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{y})\big| \Big| \mathbf{W} \Big]
\lesssim_{\mathbb{P}} \int_{0}^{C h^{-d/2 - 1}\delta_n} \sqrt{d \log \left(\frac{1}{\varepsilon h^{d/2 + 1}}\right)} d \varepsilon
\lesssim \sqrt{\log n} h^{-d/2 - 1}\delta_n
\end{align*}
and
\begin{align*}
\mathbb{E} \Big[\sup_{\lVert \mathbf{x} - \mathbf{y} \rVert_{\infty}\leq \delta_n} \big|Z^{(\boldsymbol{\nu})}(\mathbf{x}) - Z^{(\boldsymbol{\nu})}(\mathbf{y})\big|\Big]
\lesssim \int_{0}^{C h^{- 1}\delta_n} \sqrt{d \log \left(\frac{1}{\varepsilon h}\right)} d \varepsilon
\lesssim \sqrt{\log n} h^{- 1}\delta_n.
\end{align*}
In addition, using the fact that $\mathbb{E} \big[\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x})\big|\big| \mathbf{W} \big] \lesssim 1$, and choosing $a_n \asymp (\sqrt{\log n} h^{-d/2 - 1} \delta_n)^{1/2}$ and $\delta_n \asymp n^{-s}$ for some large constant $s > 0$, we conclude that
\begin{align*}
& 4 a_n \Big(\mathbb{E} \Big[\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \Big| \mathbf{W}\Big] + 1\Big)
+ \Psi_{\delta_n}(a_n) + \widehat{\Psi}_{\delta_n}(a_n) + \mathfrak{R}_M \\
& \lesssim_{\mathbb{P}} \big(\sqrt{\log n} h^{-d/2 - 1} \delta_n \big)^{1/2} + d \log(\delta_n^{-1}) ( a_n + h^{p+1})^{1/2}
\lesssim d \log(n) (a_n + h^{p+1})^{1/2},
\end{align*}
and putting all the intermediate results together, the lemma follows.
\qed
\bigskip
The proof of Theorem \ref{sa-thm: Confidence Bands} now follows directly from Lemma~\ref{sa-lem: infeasible gaussian to bahadur representation}, Lemma~\ref{sa-lem: t stats to bahadur} and Lemma~\ref{sa-lem: feasible gaussian to infeasible gaussian}. Furthermore, by definition of $\widehat{\operatorname{I}}_{\alpha}^{(\boldsymbol{\nu})}(\mathbf{x})$,
\begin{align*}
\mathbb{P}\Big[\mu^{(\boldsymbol{\nu})}(\mathbf{x}) \in \widehat{\operatorname{I}}_{\alpha}^{(\boldsymbol{\nu})}(\mathbf{x}), \text{ for all } \mathbf{x} \in \mathcal{B} \Big]
& = \mathbb{P} \Big[\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{\operatorname{T}}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq \mathcal{q}_{\alpha}\Big] \\
& = \mathbb{P} \Big[\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x})\big|\leq \mathcal{q}_{\alpha}\Big] + o(1) \\
& = \mathbb{E} \Big[\mathbb{P} \Big[\sup_{\mathbf{x} \in \mathcal{B}} \big|\widehat{Z}^{(\boldsymbol{\nu})}(\mathbf{x})\big| \leq \mathcal{q}_{\alpha}\Big| \mathbf{W} \Big]\Big] + o(1) \\
& = 1 - \alpha + o(1),
\end{align*}
which completes the proof of the theorem.
\qed
\subsection{Proof of Lemma~\ref{sa-lem: Integral: Bias}}
Follows from Lemma~\ref{sa-lem: bias} and the assumption that $\int_{\mathcal{B}} |w(\mathbf{x})| d \mathfrak{H}^{d-1}(\mathbf{x}) <\infty$.
\qed
\subsection{Proof of Lemma~\ref{sa-lem: Integral: Variance}}
Since $\mathbb{V}[\widehat{\tau}_{\mathtt{WBATE}}|\mathbf{X}] = \mathbb{V}[\int_{\mathcal{B}}\widehat{\mu}_0(\mathbf{b}) w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b})|\mathbf{X}] + \mathbb{V}[\int_{\mathcal{B}}\widehat{\mu}_1(\mathbf{b}) w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b})|\mathbf{X}]$, it is enough to consider only one treatment assignment group $t \in \{0,1\}$. In addition,
\begin{align*}
\mathbb{V} \Big[\int_{\mathcal{B}}\widehat{\mu}_t(\mathbf{b}) w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b}) \Big| \mathbf{X} \Big]
= \int_{\mathcal{B}} \int_{\mathcal{B}} \mathbb{C}\mathrm{ov} \big[\widehat{\mu}_t(\mathbf{b}_1), \widehat{\mu}_t(\mathbf{b}_2) \big| \mathbf{X} \big] w(\mathbf{b}_1) w(\mathbf{b}_2) d \mathfrak{H}^{d-1}(\mathbf{b}_1) d \mathfrak{H}^{d-1}(\mathbf{b}_2)
\end{align*}
and
\begin{align*}
\Omega_{t,\mathtt{WBATE}}
= \int_{\mathcal{B}} \int_{\mathcal{B}} \Omega_{t,\mathbf{b}_1, \mathbf{b}_2} w(\mathbf{b}_1) w(\mathbf{b}_2) d \mathfrak{H}^{d-1}(\mathbf{b}_1) d \mathfrak{H}^{d-1}(\mathbf{b}_2).
\end{align*}
Proceeding as in the proof of Lemma \ref{sa-lem: covariance}, we have
\begin{align*}
\sup_{\mathbf{b}_1, \mathbf{b}_2 \in \mathcal{B}} \big| \mathbb{C}\mathrm{ov}[\widehat{\mu}_t(\mathbf{b}_1), \widehat{\mu}_t(\mathbf{b}_2)| \mathbf{X} ] - \Omega_{t, \mathbf{b}_1, \mathbf{b}_2} \big|
\lesssim_{\mathbb{P}} \frac{\log(1/h)^{1/2}}{(n h^d)^{3/2}}.
\end{align*}
Since $K$ is supported on a compact set, let $R\in(0,\infty)$ denote the diameter of the support, and define the ``effective domain'' $\mathcal{E}(h) = \{(\mathbf{x}, \mathbf{y}) \in \mathcal{B} \times \mathcal{B}: \lVert \mathbf{x} - \mathbf{y} \rVert \leq h R\}$. Since $\mathcal{B}$ is $(d-1)$ dimensional, we have $\nu_d(\mathcal{E}(h)) \lesssim h^{d-1}$, where $\nu_{d}$ is the product measure $\mathfrak{H}^{d-1} \times \mathfrak{H}^{d-1}$. Therefore,
\begin{align*}
&\Big| \mathbb{V} \Big[\int_{\mathcal{B}}\widehat{\mu}_t(\mathbf{b}) w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b}) \Big| \mathbf{X} \Big] - \Omega_t \Big|\\
&= \Big| \int_{\mathcal{B}} \int_{\mathcal{B}} \big( \mathbb{C}\mathrm{ov}[\widehat{\mu}_t(\mathbf{b}_1), \widehat{\mu}_t(\mathbf{b}_2)| \mathbf{X}] - \Omega_{t,\mathbf{b}_1, \mathbf{b}_2} \big) w(\mathbf{b}_1) w(\mathbf{b}_2) d \mathfrak{H}^{d-1}(\mathbf{b}_1) d \mathfrak{H}^{d-1}(\mathbf{b}_2)\Big|\\
& \lesssim \sup_{\mathbf{b}_1, \mathbf{b}_2 \in \mathcal{B}} \big| \mathbb{C}\mathrm{ov}[\widehat{\mu}_t(\mathbf{b}_1), \widehat{\mu}_t(\mathbf{b}_2)|\mathbf{X}] -
\Omega_{t, \mathbf{b}_1, \mathbf{b}_2}\big| \int_{\mathcal{B}} \int_{\mathcal{B}} \mathds{1}((\mathbf{b}_1, \mathbf{b}_2) \in \mathcal{E}(h)) w(\mathbf{b}_1) w(\mathbf{b}_2) d \mathfrak{H}^{d-1}(\mathbf{b}_1) d \mathfrak{H}^{d-1}(\mathbf{b}_2) \\
&\lesssim_{\mathbb{P}} h^{d-1}\frac{\log(1/h)^{1/2}}{(n h^d)^{3/2}} = o_{\mathbb{P}}((n h)^{-1}),
\end{align*}
because $\frac{\log(1/h)}{n h^d} = o(1)$. This proves the first claim. Next,
\begin{align*}
\Omega_{t,\mathtt{WBATE}}
& = \int_{\mathcal{B}} \int_{\mathcal{B}} \Omega_{t,\mathbf{b}_1, \mathbf{b}_2} w(\mathbf{b}_1) w(\mathbf{b}_2) d \mathfrak{H}^{d-1}(\mathbf{b}_1) d \mathfrak{H}^{d-1}(\mathbf{b}_2)\\
& \leq \sup_{\mathbf{b}_1, \mathbf{b}_2 \in \mathcal{B}} \big|\Omega_{t,\mathbf{b}_1, \mathbf{b}_2}\big| \int_{\mathcal{B}} \int_{\mathcal{B}} \mathds{1}((\mathbf{b}_1, \mathbf{b}_2) \in \mathcal{E}(h)) w(\mathbf{b}_1) w(\mathbf{b}_2) d \mathfrak{H}^{d-1}(\mathbf{b}_1) d \mathfrak{H}^{d-1}(\mathbf{b}_2) \\
& \lesssim (n h^d)^{-1} \nu_d(\mathcal{E}(h)) \lesssim (n h)^{-1},
\end{align*}
which verifies the upper bound. For the lower bound, let $\mathbf{b}_1 \in \mathcal{B}$ and $\mathbf{b}_2 = \mathbf{b}_1 + h \boldsymbol{\delta}$ for some vector $\boldsymbol{\delta}$ such that $\sup_{\mathbf{x} \in \mathcal{X}}K_h(\mathbf{x} - \mathbf{b}_1) K_h(\mathbf{x} - \mathbf{b}_2) > 0$. For multi-indexes $\mathbf{u}$ and $\mathbf{v}$, and using change of variables, a typical element of $\boldsymbol{\Sigma}_{t,\mathbf{b}_1, \mathbf{b}_2}$ is
\begin{align*}
& \mathbb{E} \Big[\Big(\frac{\mathbf{X}_i - \mathbf{b}_1}{h}\Big)^{\mathbf{u}} \Big(\frac{\mathbf{X}_i - \mathbf{b}_1 - \boldsymbol{\delta} h}{h}\Big)^{\mathbf{v}} \frac{1}{h^d} K \Big(\frac{\mathbf{X}_i - \mathbf{b}_1}{h}\Big) K \Big(\frac{\mathbf{X}_i - \mathbf{b}_1 - \boldsymbol{\delta} h}{h}\Big) \sigma_t^2(\mathbf{X}_i) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big] \\
& = \int_{\mathbf{b}_1 + h \mathcal{A}_t} \mathbf{s}^{\mathbf{u}} (\mathbf{s} - \boldsymbol{\delta})^{\mathbf{v}} K(\mathbf{s}) K(\mathbf{s} + \boldsymbol{\delta}) \sigma_t^2(\mathbf{b}_1 + h\mathbf{s}) f(\mathbf{s}) d \mathbf{s} \gtrsim 1,
\end{align*}
which implies that $|\Omega_{t,\mathbf{b}_1, \mathbf{b}_2}| \gtrsim (n h^d)^{-1}$ for $(\mathbf{b}_1, \mathbf{b}_2)$ on a set $\mathcal{E}'(h)$ such that $\nu_d(\mathcal{E}'(h)) \gtrsim h^{d-1}$. This verifies lower bound in the second claim.
The third and final claim of the lemma follows from Lemma~\ref{sa-lem: covariance} and the same analysis as above.
\subsection{Proof of Theorem~\ref{sa-thm: Integral: MSE Expansion}}
Follows from Lemma \ref{sa-lem: Integral: Bias} and Lemma \ref{sa-lem: Integral: Variance}.
\qed
\subsection{Proof of Theorem~\ref{sa-thm: Integral: Asymptotic Normality}}
Since $\widehat{\tau}_{\mathtt{WBATE}} - \tau_{\mathtt{WBATE}} = (\widehat{\mu}_{1,\mathtt{WBATE}} - \mu_{1,\mathtt{WBATE}}) - (\widehat{\mu}_{0,\mathtt{WBATE}} - \mu_{0,\mathtt{WBATE}})$, it is enough to start with only one treatment assignment group $t \in \{0,1\}$. Furthermore,
\begin{align*}
\widehat{\mu}_{t,\mathtt{WBATE}} - \mu_{t,\mathtt{WBATE}}
&= \int_{\mathcal{B}} (\widehat{\mu}_{1}(\mathbf{b}) - \mu_1(\mathbf{b}) ) w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b})\\
&= \int_{\mathcal{B}} \mathbf{e}_1^{\top} \boldsymbol{\Gamma}_{t,\mathbf{b}}^{-1} \mathbf{Q}_{t,\mathbf{b}} w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b})
+ \int_{\mathcal{B}} \mathbf{e}_1^{\top} (\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{b}}^{-1} - \boldsymbol{\Gamma}_{t,\mathbf{b}}^{-1}) \mathbf{Q}_{t,\mathbf{b}} w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b})
+ O_\mathbb{P}(h^{p+1})
\end{align*}
using Lemma~\ref{sa-lem: bias} to bound the approximation error.
For the second integral, let
\begin{align*}
\overline{\boldsymbol{\Sigma}}_{t,\mathbf{x}_1,\mathbf{x}_2}
= \mathbb{E}_n \Big[\mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{x}_1}{h}\Big) \mathbf{r}_p \Big( \frac{\mathbf{X}_i - \mathbf{x}_2}{h}\Big)^{\top} h^d K_h(\mathbf{X}_i - \mathbf{x}_1) K_h(\mathbf{X}_i - \mathbf{x}_2) \sigma_t^2(\mathbf{X}_i) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)\Big],
\end{align*}
and since $\overline{\boldsymbol{\Sigma}}_{t,\mathbf{b}_1, \mathbf{b}_2} = 0$ if $\mathbf{b}_1$ and $\mathbf{b}_2$ are farther away form each other than the diameter of $\operatorname{Supp}(K)$,
\begin{align*}
&\mathbb{E} \bigg[ \bigg(\int_{\mathcal{B}} \mathbf{e}_1^{\top} (\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{b}}^{-1} - \boldsymbol{\Gamma}_{t,\mathbf{b}}^{-1}) \mathbf{Q}_{t,\mathbf{b}} w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b}) \bigg)^2\bigg|\mathbf{X}\bigg]\\
&= \int_{\mathcal{B}} \int_{\mathcal{B}} \mathbf{e}_1^{\top} (\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{b}_1}^{-1} - \boldsymbol{\Gamma}_{t,\mathbf{b}_1}^{-1}) (n h^d)^{-1} \overline{\boldsymbol{\Sigma}}_{t,\mathbf{b}_1, \mathbf{b}_2} (\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{b}_2}^{-1} - \boldsymbol{\Gamma}_{t,\mathbf{b}_2}^{-1}) \mathbf{e}_1 w(\mathbf{b}_1) w(\mathbf{b}_2) d \mathfrak{H}^{d-1}(\mathbf{b}_1) d \mathfrak{H}^{d-1}(\mathbf{b}_2),\\
&\leq \sup_{\mathbf{b} \in \mathcal{B}} \|\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{b}}^{-1} - \boldsymbol{\Gamma}_{t,\mathbf{b}}^{-1}\|^2 \sup_{\mathbf{b}_1, \mathbf{b}_2 \in \mathcal{B}} \lVert \overline{\boldsymbol{\Sigma}}_{t,\mathbf{b}_1, \mathbf{b}_2} \rVert \sup_{\mathbf{b} \in \mathcal{B}} |w(\mathbf{b})| (n h^d)^{-1} \mathfrak{m}(\mathcal{E}(h)) \\
&\lesssim_{\mathbb{P}} (n h)^{-1},
\end{align*}
and hence $\int_{\mathcal{B}} \mathbf{e}_1^{\top} (\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{b}}^{-1} - \boldsymbol{\Gamma}_{t,\mathbf{b}}^{-1}) \mathbf{Q}_{t,\mathbf{b}} w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b}) = o_\mathbb{P}((n h)^{-1})$.
Next, using Lemma~\ref{sa-lem: Integral: Variance} and the previous results,
\begin{align*}
\widehat{\operatorname{T}}_{\mathtt{WBATE}} - \overline{\operatorname{T}}_{\mathtt{WBATE}}
= \big(\widehat{\Omega}_{\mathtt{WBATE}}^{-1/2} - \Omega_{\mathtt{WBATE}}^{-1/2}\big)
\int_{\mathcal{B}} \mathbf{e}_1^{\top} \boldsymbol{\Gamma}_{t,\mathbf{b}}^{-1} \mathbf{Q}_{t,\mathbf{b}} d \mathfrak{H}^{d-1}(\mathbf{b}) + o_{\mathbb{P}}(1)
= o_{\mathbb{P}}(1),
\end{align*}
where
\begin{align*}
\overline{\operatorname{T}}_{\mathtt{WBATE}} = \Omega_{\mathtt{WBATE}}^{-1/2} \int_{\mathcal{B}} \mathbf{e}_1^{\top} \boldsymbol{\Gamma}_{t,\mathbf{b}}^{-1} \mathbf{Q}_{t,\mathbf{b}} d \mathfrak{H}^{d-1}(\mathbf{b}).
\end{align*}
Finally, we apply the Berry-Esseen lemma to the statistic $\overline{\operatorname{T}}_{w} = \sum_{i = 1}^n Z_i$, where
\begin{align*}
Z_i
= n^{-1}\Omega_{\mathtt{WBATE}}^{-1/2} \int_{\mathcal{B}} \mathbf{e}_1^{\top} \boldsymbol{\Gamma}_{t,\mathbf{b}}^{-1} \mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{b}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{b}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) u_i w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b}),
\end{align*}
which satisfies $\mathbb{E}[Z_i] = 0$. The definition of $\Omega_{\mathtt{WBATE}}$ implies that $\sum_{i = 1}^n\mathbb{V}[Z_i] = \Omega_{\mathtt{WBATE}}^{-1/2} \Omega_{\mathtt{WBATE}} \Omega_{\mathtt{WBATE}}^{-1/2} = 1$. Hence, it remains to bound
\begin{align*}
\sum_{i = 1}^n \mathbb{E}[|Z_i|^3]
= n^{-3} \Omega_{\mathtt{WBATE}}^{-3/2} \sum_{i = 1}^n \mathbb{E} \bigg[\bigg| \int_{\mathcal{B}} \mathbf{e}_1^{\top} \boldsymbol{\Gamma}_{t,\mathbf{b}}^{-1} \mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{b}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{b}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) u_i w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b}) \bigg|^3\bigg].
\end{align*}
Let $R$ denote the diameter of the (compact) support of $K$, and define $\mathcal{E}(h) = \{(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3) \in \mathcal{B}^3: \lVert \mathbf{b}_i - \mathbf{b}_j\rVert \leq R, j = 1,2,3\}$. Since $\mathcal{B}$ is $d-1$ dimensional, $\mathfrak{m}(\mathcal{E}(h)) \lesssim h^{2d-2}$. Then,
\begin{align*}
& \mathbb{E} \bigg[\bigg| \int_{\mathcal{B}} \mathbf{e}_1^{\top} \boldsymbol{\Gamma}_{t,\mathbf{b}}^{-1} \mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{b}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{b}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) u_i w(\mathbf{b}) d \mathfrak{H}^{d-1}(\mathbf{b}) \bigg|^3\bigg] \\
& \leq \mathbb{E} \bigg[\int_{\mathbf{b}_1 \in \mathcal{B}} \int_{\mathbf{b}_2 \in \mathcal{B}} \int_{\mathbf{b}_3 \in \mathcal{B}} |G(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3)| w(\mathbf{b}_1) w(\mathbf{b}_2) w(\mathbf{b}_3) d \mathfrak{H}^{d-1}(\mathbf{b}_1) d \mathfrak{H}^{d-1}(\mathbf{b}_2) d \mathfrak{H}^{d-1}(\mathbf{b}_3)\bigg] \\
& \lesssim \mathfrak{m}(\mathcal{E}(h)) \sup_{\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3 \in \mathcal{B}} \mathbb{E}[|G(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3)|],
\end{align*}
where $G(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3) = g(\mathbf{X}_i, u_i, \mathbf{b}_1) g(\mathbf{X}_i, u_i, \mathbf{b}_2) g(\mathbf{X}_i, u_i, \mathbf{b}_3)$ with
\begin{align*}
g(\mathbf{X}_i, u_i, \mathbf{b})
= \mathbf{e}_1^{\top} \boldsymbol{\Gamma}_{t,\mathbf{b}}^{-1} \mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{b}}{h}\Big) K_h(\mathbf{X}_i - \mathbf{b}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) u_i.
\end{align*}
Proceeding as in the proof of Lemma~\ref{sa-lem: gram} and Lemma~\ref{sa-lem: Integral: Variance}, it can be shown that
\begin{align*}
\sup_{\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3 \in \mathcal{B}} \mathbb{E}[|G(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3)|] \lesssim h^{-2d}
\end{align*}
provided that $\frac{\log(1/n)}{n h^d} = o(1)$. Therefore, together with the rate of $\Omega_{\mathtt{WBATE}}$ from Lemma~\ref{sa-lem: Integral: Variance}, we have $\sum_{i = 1}^n \mathbb{E}[|Z_i^3|] \lesssim (n h)^{-1/2}$, and the result follows.
\qed
\subsection{Proof of Theorem~\ref{sa-thm: Convergence Rate for max}}
Follows by Theorem \ref{sa-thm: Convergence Rates} after noting that $ |\widehat{\tau}_{\mathtt{LBATE}} - \tau_{\mathtt{LBATE}} | \leq \sup_{\mathbf{x}\in\mathcal{B}} |\widehat{\tau}(\mathbf{x}) - \tau(\mathbf{x})|$.
\qed
\subsection{Proof of Theorem~\ref{sa-thm: Confidence Bands for max}}
Consider the event $E = \Big\{\sup_{\mathbf{b} \in \mathcal{B}} \frac{|\widehat{\tau}(\mathbf{b}) - \tau(\mathbf{b})|}{\widehat{\Omega}_{\mathbf{b}, \mathbf{b}}^{1/2}} \leq \mathcal{q}_{\alpha}\Big\}$. Theorem~\ref{sa-thm: Confidence Bands} implies that $\mathbb{P}(E) = 1 - \alpha + o(1)$. On the event $E$, we also have
\begin{align*}
\widehat{\tau}(\mathbf{b}) - \mathcal{q}_{\alpha} \widehat{\Omega}_{\mathbf{b},\mathbf{b}}^{1/2} \leq \tau(\mathbf{b})
\leq \widehat{\tau}(\mathbf{b}) + \mathcal{q}_{\alpha} \widehat{\Omega}_{\mathbf{b},\mathbf{b}}^{1/2}, \qquad \forall \; \mathbf{b} \in \mathcal{B},
\end{align*}
which implies
\begin{align*}
\sup_{\mathbf{b} \in \mathcal{B}} \widehat{\tau}(\mathbf{b}) - \mathcal{q}_{\alpha} \widehat{\Omega}_{\mathbf{b},\mathbf{b}}^{1/2}
\leq \sup_{\mathbf{b} \in \mathcal{B}} \tau(\mathbf{b})
\leq \sup_{\mathbf{b} \in \mathcal{B}} \widehat{\tau}(\mathbf{b}) + \mathcal{q}_{\alpha} \widehat{\Omega}_{\mathbf{b},\mathbf{b}}^{1/2}.
\end{align*}
The stated result then follows.
\qed
\subsection{Proof of Theorem~\ref{sa-lem: sa thm}}\label{sa-sec: proof of sa lemma}
We will use a truncation argument. Let $\kappa_n > 0$ be the level of truncation. For each $r \in \mathcal{R}$, define
\begin{align*}
\tilde{r}(y) = r(y) \mathds{1}(|y| \leq \kappa_n), \qquad y \in \mathbb{R},
\end{align*}
and define the class $\tilde{\mathcal{R}} = \{\tilde{r}: r \in \mathcal{R}\}$. For an overview of our argument, suppose $Z_n^R$ is some mean-zero Gaussian process indexed by $\mathcal{G} \times \mathcal{R} \cup \mathcal{G} \times \tilde{\mathcal{R}}$, whose existence will be shown below, then we can decompose by:
\begin{align*}
R_n(g,r) - Z_n^R(g,r)
& = \big[R_n(g,\tilde{r}) - Z_n^R(g,\tilde{r})\big] + \big[R_n(g,r) - R_n(g,\tilde{r})\big] + \big[Z_n^R(g,r) - Z_n^R(g,\tilde{r})\big].
\end{align*}
\subsubsection*{Part 1: Strong approximation for truncated residual empirical process.}
Observe that $\mathtt{M}_{\tilde{\mathcal{R}},\mathcal{Y}} \lesssim \kappa_n$ and $\mathtt{pTV}_{\tilde{\mathcal{R}},\mathcal{Y}} \lesssim \kappa_n$, and $\tilde{\mathcal{R}}$ is a VC-type class with envelope $M_{\tilde{\mathcal{R}},\mathcal{Y}} = M_{\mathcal{R},\mathcal{Y}} \mathds{1}(|\cdot| \leq \kappa_n)$ over $\mathcal{Y}$ with constants $\mathtt{c}_{\mathcal{R},\mathcal{Y}}$ and $\mathtt{d}_{\mathcal{R},\mathcal{Y}}$. Then, \citet[Theorem 2]{Cattaneo-Yu_2025_AOS} with $\mathtt{v} = \kappa_n$ and $\alpha = 0$ for the class of functions $\mathcal{G}$ and $\tilde{\mathcal{R}}$ implies on a possibly enlarged probability space, there exists a sequence of mean-zero Gaussian processes $(Z_n^R(g,r): (g,r)\in \mathcal{G}\times \tilde{\mathcal{R}})$ with almost sure continuous trajectories on $(\mathcal{G} \times \tilde{\mathcal{R}}, \rho_{\mathbb{P}})$ such that $\mathbb{E}[R_n(g_1, r_1) R_n(g_2, r_2)] = \mathbb{E}[Z^R_n(g_1, r_1) Z^R_n(g_2, r_2)]$ for all $(g_1, r_1), (g_2, r_2) \in \mathcal{G} \times \tilde{\mathcal{R}}$, and
\begin{align*}
& \mathbb{E}[\lVert R_n(g,\tilde{r}) - Z_n^R(g,\tilde{r}) \rVert_{\mathcal{G} \times \mathcal{R}}] \\
& \leq C_1 \mathtt{v} \kappa_n \bigg(\sqrt{d} \min\Big\{\frac{(\mathtt{c}_1^d \mathtt{M}_{\mathcal{G}}^{d+1} \mathtt{TV}^d \mathtt{E}_{\mathcal{G}})^{\frac{1}{2d+2}} }{n^{1/(2d+2)}},
\frac{(\mathtt{c}_1^{\frac{d}{2}} \mathtt{c}_2^{\frac{d}{2}}\mathtt{M}_{\mathcal{G}} \mathtt{TV}^{\frac{d}{2}} \mathtt{E}_{\mathcal{G}} \mathtt{L}^{\frac{d}{2}})^{\frac{1}{d+2}}}{n^{1/(d+2)}} \Big\} ((\mathtt{d} + \mathtt{k})\log (\mathtt{c} n))^{3/2}
+ \frac{(\mathtt{d} + \mathtt{k})\log (\mathtt{c} n)}{\sqrt{n}} \mathtt{M}_{\mathcal{G}} \bigg) \\
& = C_1 \mathtt{v} \kappa_n \bigg(\sqrt{d} \mathtt{r}_n ((\mathtt{d} + \mathtt{k})\log (\mathtt{c} n))^{\frac{3}{2}}
+ \frac{(\mathtt{d} + \mathtt{k})\log (\mathtt{c} n)}{\sqrt{n}} \mathtt{M}_{\mathcal{G}}\bigg),
\end{align*}
where $C_1$ is some positive universal constant. Notice that we use $\mathtt{TV} = \max \{\mathtt{TV}_{\mathcal{G}}, \mathtt{TV}_{\mathcal{G} \times \mathscr{U}_{\mathcal{R}},\mathcal{Q}_\mathcal{G}}\}$ as an upper bound for $\max \{\mathtt{TV}_{\mathcal{G}}, \mathtt{TV}_{\mathcal{G} \times \mathscr{V}_{\tilde{\mathcal{R}}},\mathcal{Q}_\mathcal{G}}\}$, and similarly $\mathtt{L}$ as an upper bound for $\max \{\mathtt{L}_{\mathcal{G}}, \mathtt{L}_{\mathcal{G} \times \mathscr{V}_{\tilde{\mathcal{R}}},\mathcal{Q}_\mathcal{G}}\}$.
In the special case that $\mathcal{R} = \{r_{\ast}\}$ is a singleton, take $\tilde{y}_i = r_{\ast}(y_i) \mathds{1}(|y_i| \leq \kappa_n)/(\mathtt{v} \kappa_n)$, then we have $\mathbb{E}[\exp(|\tilde{y}_i|)] \leq 2$. Also $\tilde{y}_i$ is supported on $\tilde{\mathcal{Y}} = [-1,1]$. Moreover,
\begin{align*}
\frac{1}{\mathtt{v} \kappa_n}R_n(g, \tilde{r}_{\ast}) = \frac{1}{n} \sum_{i = 1}^n g(\mathbf{x}_i) (\tilde{y}_i - \mathbb{E}[\tilde{y}_i]), \qquad g \in \mathcal{G}.
\end{align*}
In particular, the right hand side can be viewed as a residual empirical process based on sample $(\mathbf{x}_i, \tilde{y}_i), 1 \leq i \leq n$, indexed by $\mathcal{G} \times \{\operatorname{Id}\}$, where $\operatorname{Id}: \mathbb{R} \to \mathbb{R}$ is the identity function. Then we can apply \citet[Theorem 2]{Cattaneo-Yu_2025_AOS} with $\mathtt{v} = 1$ and $\alpha = 0$ on the latter empirical process to get the upper bound with $\mathtt{TV}$ and $\mathtt{L}$ replaced by $\mathtt{TV}_{\text{sing}}$ and $\mathtt{L}_{\text{sing}}$.
\subsubsection*{Part 2: Truncation error for the empirical process --- $\lVert R_n(g,r) - R_n(g, \widetilde{r})\rVert_{\mathcal{G} \times \mathcal{R}}$}
Consider the class of differences due to truncation, that is, $\Delta \mathcal{R} = \{r - \tilde{r}: r \in \mathcal{R}\}$. Our assumptions imply $\mathcal{G} \times \Delta \mathcal{R}$ is VC-type in the sense that for all $0 <\varepsilon < 1$,
\begin{align*}
\sup_{\mathcal{Q}} N(\mathcal{G} \times \Delta \mathcal{R}, \left\lVert\cdot\right\rVert_{\mathcal{Q},2}, \varepsilon \lVert \mathtt{M}_{\mathcal{G}} (M_{\mathcal{R},\mathcal{Y}} - M_{\tilde{\mathcal{R}},\mathcal{Y}}) \rVert_{\mathcal{Q},2}) \leq \mathtt{c}_{\mathcal{G}} \mathtt{c}_{\mathcal{R},\mathcal{Y}} (\varepsilon^2/4)^{-\mathtt{d}_{\mathcal{G}} - \mathtt{d}_{\mathcal{R},\mathcal{Y}}},
\end{align*}
where $\sup$ is over all finite discrete measure on $\mathbb{R}^{d+1}$, and $M_{\tilde{\mathcal{R}},\mathcal{Y}}(y) = M_{\mathcal{R},\mathcal{Y}}(y) \mathds{1}(|y| \leq \kappa_n)$. We can check that $\mathtt{M}_{\mathcal{G}} (M_{\mathcal{R},\mathcal{Y}} - M_{\tilde{\mathcal{R}},\mathcal{Y}})$ is an envelope function for $\mathcal{G} \times \Delta \mathcal{R}$, since all functions in $\Delta \mathcal{R}$ are evaluated to zero on $[-\kappa_n, \kappa_n]$. Denote $\mathbf{X} = (\mathbf{x}_i)_{1 \leq i \leq n}$,
\begin{align*}
\mathbb{E} \Big[\max_{1 \leq i \leq n}\mathtt{M}_{\mathcal{G}}^2 (M_{\mathcal{R},\mathcal{Y}}(y_i) - M_{\tilde{\mathcal{R}},\mathcal{Y}}(y_i))^2\Big|\mathbf{X}\Big]^{\frac{1}{2}}
& \lesssim \mathtt{M}_{\mathcal{G}} \mathbb{E} \Big[ \big( \max_{1 \leq i \leq n}M_{\mathcal{R},\mathcal{Y}}(y_i)\big)^2\Big| \mathbf{X} \Big]^{\frac{1}{2}} \lesssim \mathtt{M}_{\mathcal{G}} n^{\frac{1}{2+v}}, \\
\sup_{(g,r) \in \mathcal{G} \times \mathcal{R}} \mathbb{E}[g(\mathbf{x}_i)^2 r(y_i)^2 \mathds{1}(|y_i| \geq \kappa_n^{1/\alpha})]^{\frac{1}{2}}
& \lesssim \sup_{(g,r) \in \mathcal{G} \times \mathcal{R}} \mathbb{E} \left[g(\mathbf{x}_i)^2 \mathbb{E}[r(y_i)^{2+v}|\mathbf{x}_i]^{\frac{2}{2+v}} \mathbb{P}(|y_i| \geq \kappa_n| \mathbf{x}_i)^{\frac{v}{2+v}}\right] \\
& \lesssim \sqrt{\mathtt{M}_{\mathcal{G}} \mathtt{E}_{\mathcal{G}} \kappa_n}.
\end{align*}
By Jensen's inequality, we also have
\begin{align*}
& \mathbb{E} \Big[\max_{1 \leq i \leq n}\mathtt{M}_{\mathcal{G}}^2 (\mathbb{E}[M_{\mathcal{R},\mathcal{Y}}(y_i) - M_{\tilde{\mathcal{R}},\mathcal{Y}}(y_i)|\mathbf{x}_i])^2\Big|\mathbf{X}\Big]^{\frac{1}{2}} \lesssim \mathtt{M}_{\mathcal{G}} n^{\frac{1}{2+v}}, \\
& \sup_{(g,r) \in \mathcal{G} \times \mathcal{R}} \mathbb{E}[g(\mathbf{x}_i)^2 \mathbb{E}[r(y_i) - \tilde{r}(y_i)|\mathbf{x}_i]^2]^{\frac{1}{2}} \lesssim \sqrt{\mathtt{M}_{\mathcal{G}} \mathtt{E}_{\mathcal{G}} \kappa_n^{-v}}, \\
& \mathbb{E}[\mathtt{M}_{\mathcal{G}}^2(M_{\mathcal{R},\mathcal{Y}}(y_i) - M_{\tilde{\mathcal{R}},\mathcal{Y}}(y_i))^2]^{1/2} \lesssim \mathtt{M}_{\mathcal{G}} \kappa_n^{-v/2}.
\end{align*}
Denote $A = (\mathtt{c}_{\mathcal{G}} \mathtt{c}_{\mathcal{R}})^{\frac{1}{2 \mathtt{d}_{\mathcal{G}} + 2 \mathtt{d}_{\mathcal{R}}}}/4$ and $D = 2 \mathtt{d}_{\mathcal{G}} + 2 \mathtt{d}_{\mathcal{R}}$, \citet[Corollary 5.1]{Chernozhukov-Chetverikov-Kato_2014b_AoS} gives
\begin{align*}
\mathbb{E} \left[\lVert R_n(g, r) - R_n(g \widetilde{r})\rVert_{\mathcal{G} \times \mathcal{R}}\right] &
\lesssim \mathbb{E} \left[\sup_{g \in \mathcal{G}} \sup_{h \in \Delta \mathcal{R}} \frac{1}{\sqrt{n}} \sum_{i = 1}^n g(\mathbf{x}_i)(h(y_i) - \mathbb{E}[h(y_i)|\mathbf{x}_i])\right]\\
& \lesssim \sqrt{D \mathtt{M}_{\mathcal{G}} \mathtt{E}_{\mathcal{G}}\kappa_n^{-v} \log (A \sqrt{\mathtt{M}_{\mathcal{G}}/\mathtt{E}_{\mathcal{G}}})} + \frac{D \mathtt{M}_{\mathcal{G}} n^{\frac{1}{2+V}}}{\sqrt{n}} \log (A \sqrt{\mathtt{M}_{\mathcal{G}}/\mathtt{E}_{\mathcal{G}}})\\
& \lesssim \sqrt{D \log (A \sqrt{\mathtt{M}_{\mathcal{G}}/\mathtt{E}_{\mathcal{G}}})} \sqrt{\mathtt{M}_{\mathcal{G}} \mathtt{E}_{\mathcal{G}}} \kappa_n^{-v/2} + \frac{D \log (A \sqrt{\mathtt{M}_{\mathcal{G}}/\mathtt{E}_{\mathcal{G}}}) \mathtt{M}_{\mathcal{G}}}{\sqrt{n^{\frac{v}{2+v}}}}.
\end{align*}
\subsubsection*{Part 3: Truncation error for the Gaussian process --- $\lVert Z_n^R(g,r) - Z_n^R(g,\tilde{r}) \rVert_{\mathcal{G} \times \mathcal{R}}$}
Our assumptions imply $\mathcal{G} \times \tilde{\mathcal{R}} \cup \mathcal{G} \times \mathcal{R}$ is VC-type w.r.p envelope function $2 \mathtt{M}_{\mathcal{G}} \mathtt{M}_{\mathcal{R},\mathcal{Y}}$ in the sense that for all $0 <\varepsilon < 1$,
\begin{align*}
\sup_{Q} N(\mathcal{G} \times \mathcal{R} \cup \mathcal{G} \times \tilde{\mathcal{R}}, \left\lVert\cdot\right\rVert_{Q,2}, 2 \varepsilon \lVert \mathtt{M}_{\mathcal{G}} \mathtt{M}_{\mathcal{R},\mathcal{Y}} \rVert_{Q,2}) \leq \mathtt{c}_{\mathcal{G}} \mathtt{c}_{\mathcal{R}} (\varepsilon^2/4)^{-\mathtt{d}_{\mathcal{G}}-\mathtt{d}_{\mathcal{R}}},
\end{align*}
where $\sup$ is over all finite discrete measure on $\mathbb{R}^{d+1}$. Hence $\mathcal{G} \times \tilde{\mathcal{R}} \cup \mathcal{G} \times \mathcal{R}$ is pre-Gaussian, and on some probability space, there exists a mean-zero Gaussian process $\bar{Z}_n^R$ indexed by $\mathcal{F} = \mathcal{G} \times \tilde{\mathcal{R}} \cup \mathcal{G} \times \mathcal{R}$ with the same covariance structure as $R_n$, and has almost sure continuous path w.r.p the metric $\rho$, given by
\begin{align*}
\rho((g_1,r_1),(g_2,r_2)) = \mathbb{E}[(Z_n^R(g_1, r_1) - Z_n^R(g_2, r_2))^2]^{\frac{1}{2}} = \mathbb{E}[(R_n(g_1, r_1) - R_n(g_2, r_2))^2]^{\frac{1}{2}}, (g_1, r_1), (g_2, r_2) \in \mathcal{F}.
\end{align*}
Recall the definition of $\mathcal{G} \times \Delta\mathcal{R}$ in Part 2. Then, we have shown previously that
\begin{align*}
\sigma \equiv \sup_{f \in \mathcal{G} \times \Delta \mathcal{R}}\rho(f,f) \leq \sqrt{\mathtt{M}_{\mathcal{G}} \mathtt{E}_{\mathcal{G}} \kappa_n^{-v}},
\end{align*}
Our assumptions imply for all $0 < \varepsilon < 1$,
\begin{align*}
N(\mathcal{G} \times \mathcal{R} \cup \mathcal{G} \times \tilde{\mathcal{R}}, \rho, \rho(2 \varepsilon \mathtt{M}_{\mathcal{G}} M_{\mathcal{R},\mathcal{Y}}, 2 \varepsilon \lVert \mathtt{M}_{\mathcal{G}} M_{\mathcal{R},\mathcal{Y}})^{1/2}) \leq \mathtt{c}_{\mathcal{G}} \mathtt{c}_{\mathcal{R}} (\varepsilon^2/4)^{-\mathtt{d}_{\mathcal{G}}-\mathtt{d}_{\mathcal{R}}}
\end{align*}
Denote $A = (\mathtt{c}_{\mathcal{G}} \mathtt{c}_{\mathcal{R}})^{\frac{1}{2 \mathtt{d}_{\mathcal{G}} + 2 \mathtt{d}_{\mathcal{R}}}}/4$ and $D = 2 \mathtt{d}_{\mathcal{G}} + 2 \mathtt{d}_{\mathcal{R}}$. Then, by \citet[Corollary 2.2.8]{van-der-Vaart-Wellner_1996_Book}, choose any $(g_0, r_0) \in \mathcal{G} \times \mathcal{R}$, we have
\begin{align*}
\mathbb{E}\bigg[\lVert \bar{Z}_n^R(g,r) - \bar{Z}_n^R(g,\tilde{r}) \rVert_{\mathcal{G} \times \mathcal{R}}\bigg] & \lesssim \mathbb{E} \left[\left|\bar{Z}_n^R(g_0, r_0) - \bar{Z}_n^R(g_0, \tilde{r}_0)\right|\right] + \int_{0}^{\sigma} \sqrt{\log\Big( \mathtt{c}_{\mathcal{G}} \mathtt{c}_{\mathcal{R}}\big(\frac{\mathtt{M}_{\mathcal{G}}}{\varepsilon}\big)^{\mathtt{d}_{\mathcal{G}} + \mathtt{d}_{\mathcal{R}}}\Big)} d \varepsilon \\
& \leq \sqrt{D \log (A \sqrt{\mathtt{M}_{\mathcal{G}}/\mathtt{E}_{\mathcal{G}}})} \sqrt{\mathtt{M}_{\mathcal{G}} \mathtt{E}_{\mathcal{G}}} \kappa_n^{-v/2} \\
& \lesssim \sqrt{(\mathtt{d}_{\mathcal{G}} + \mathtt{d}_{\mathcal{R},\mathcal{Y}}) \log(\mathtt{c}_{\mathcal{G}} \mathtt{c}_{\mathcal{R},\mathcal{Y}} \mathtt{k} n)} \sqrt{\mathtt{M}_{\mathcal{G}} \mathtt{E}_{\mathcal{G}}} \kappa_n^{-v/2}.
\end{align*}
Since $(\bar{Z}_n^R(g,r): g \in \mathcal{G}, r \in \mathcal{R})$ has the same distribution as $(Z_n^R(g,r): g \in \mathcal{G}, r \in \mathcal{R})$, we know from Vorob'ev–Berkes–Philipp theorem \citep[Theorem 1.31]{dudley2014uniform} that $\bar{Z}_n^R$ can be constructed on the same probability space as $(\mathbf{x}_i,y_i)_{1 \leq i \leq n}$ and $Z_n^R$, such that $\bar{Z}_n^R$ and $Z_n^R$ coincide on $\mathcal{G} \times \mathcal{R}$. By an abuse of notation, call $\bar{Z}_n^R$ now $Z_n^R$, the outputted Gaussian process.
\subsubsection*{Part 4: Putting Together}
If follows from the definition of $\tilde{\mathcal{R}}$ and the previous three parts that if we choose $\kappa_n$ such that
\begin{align*}
\mathtt{r}_n \kappa_n \asymp \sqrt{\mathtt{M}_{\mathcal{G}} \mathtt{E}_{\mathcal{G}}} \kappa_n^{-v/2},
\end{align*}
then the approximation error can be bounded by
\begin{align*}
\mathbb{E} \big[\lVert R_n - Z_n^R \rVert_{\mathcal{G} \times \mathcal{R}} \big] & \lesssim (\mathtt{d}\log (\mathtt{c} n))^{3/2} \mathtt{r}_n^{\frac{v}{v +2}}(\sqrt{\mathtt{M}_{\mathcal{G}}\mathtt{E}_{\mathcal{G}}})^{\frac{2}{v+2}} + \mathtt{d}\log(\mathtt{c} n) \mathtt{M}_{\mathcal{G}} n^{-\frac{v/2}{2+v}} \\
& \qquad + \mathtt{d}\log(\mathtt{c} n) \mathtt{M}_{\mathcal{G}} n^{-1/2} \Big(\frac{\sqrt{\mathtt{M}_{\mathcal{G}} \mathtt{E}_{\mathcal{G}}}}{\mathtt{r}_n}\Big)^{\frac{2}{v+2}},
\end{align*}
where $\mathtt{d} = \mathtt{d}_{\mathcal{G}} + \mathtt{d}_{\mathcal{R},\mathcal{Y}} + \mathtt{k}$, and $\mathtt{c} = \mathtt{c}_{\mathcal{G}} \mathtt{c}_{\mathcal{R},\mathcal{Y}} \mathtt{k}$.
\qed
\bibliography{CTY_2025_BDD-Location--bib}
\bibliographystyle{plainnat}