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.
154,563 characters
Local Regression Distribution EstimatorsSupplemental Appendix
\title{
Local Regression Distribution Estimators\thanks{Cattaneo gratefully acknowledges financial support from the National Science Foundation through grant SES-1947805, and Jansson gratefully acknowledges financial support from the National Science Foundation through grant SES-1947662 and the research support of CREATES.}\bigskip\\ Supplemental Appendix\bigskip}
\author{
Matias D. Cattaneo\thanks{Department of Operations Research and Financial Engineering, Princeton University.}
\and
Michael Jansson\thanks{Department of Economics, UC Berkeley and CREATES.}
\and
Xinwei Ma\thanks{Department of Economics, UC San Diego.}
}
\maketitle
\begin{abstract}
\noindent This Supplemental Appendix contains general theoretical results encompassing those discussed in the main paper, includes proofs of those general results, and discusses additional methodological and technical results.
\end{abstract}
\vfill
\thispagestyle{empty}
\clearpage
\tableofcontents
\clearpage
\onehalfspacing
\pagestyle{plain}
\section{Setup}\label{section:setup}
Suppose $x_1,x_2,\cdots,x_n$ is a random sample from a univariate distribution with cumulative distribution function $F(\cdot)$. Also assume the distribution function admits a (sufficiently accurate) linear-in-parameters local approximation near an evaluation point $\mathsf{x}$:
\begin{align*}
\varrho(h,\mathsf{x}) := \sup_{|x-\mathsf{x}|\leq h}\left|F(x) - R(x-\mathsf{x})^\prime \theta(\mathsf{x})\right|\ \text{is small for $h$ small},
\end{align*}
where $R(\cdot)$ is a known basis function. The parameter $\theta(\mathsf{x})$ can be estimated by the following local $L^2$ method:
\begin{align}\label{eq:local L2 estimator}
\hat\theta_{G} &= \operatorname*{argmin}_{\theta} \int_{\mathcal{X}}\left( \hat{F}(u) - R(u-\mathsf{x})^\prime\theta \right)^2 \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right)\mathrm{d} G(u),\qquad \hat{F}(u) = \frac{1}{n}\sum_{i=1}^n \mathds{1}(x_i\leq u),
\end{align}
where $K(\cdot)$ is a kernel function, $\mathcal{X}$ is the support of $F(\cdot)$, and $G(\cdot)$ is a known weighting function to be specified later. The local $L^2$ estimator \eqref{eq:local L2 estimator} is closely related to another estimator, which is constructed by local regression:
\begin{align}\label{eq:local regression estimator}
\hat\theta &= \operatorname*{argmin}_{\theta} \sum_{i=1}^n \left( \hat{F}(x_i) - R(x_i-\mathsf{x})^\prime\theta \right)^2 \frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right).
\end{align}
The local regression estimator can be equivalently expressed as $\hat\theta_{\hat{F}}$, meaning that it can be viewed as a special case of the local $L^2$ estimator, with $G(\cdot)$ in \eqref{eq:local L2 estimator} replaced by the empirical distribution function $\hat{F}(\cdot)$.
For future reference, we first discuss some of the notation we use in the main paper and this Supplemental Appendix (SA). For a function $g(\cdot)$, we denote its $j$-th derivative as $g^{(j)}(\cdot)$. For simplicity, we also use the ``dot'' notation to denote the first derivative: $\dot{g}(\cdot) = g^{(1)}(\cdot)$. Assume $g(\cdot)$ is suitably smooth on $[\mathsf{x}-\delta,\mathsf{x})\cup(\mathsf{x},\mathsf{x}+\delta]$ for some $\delta>0$, but not necessarily continuous or differentiable at $\mathsf{x}$. If $g(\cdot)$ and its one-sided derivatives can be continuously extended to $\mathsf{x}$ from the two sides, we adopt the following special notation:
\begin{align*}
g^{(j)}_u &= \mathds{1}(u<0)g^{(j)}(\mathsf{x}-) + \mathds{1}(u\geq 0)g^{(j)}(\mathsf{x}+).
\end{align*}
With $j=0$, the above is simply $g_u = \mathds{1}(u<0)g(\mathsf{x}-)+\mathds{1}(u\geq 0)g(\mathsf{x}+)$. Also for $j=1$, we use $\dot{g}_u = g^{(1)}_u$. Convergence in probability and in distribution are denoted by $\overset{\mathbb{P}}{\to}$ and $\rightsquigarrow$, respectively, and limits are taken with respect to the sample size $n$ going to infinity unless otherwise specified. We use $|\cdot|$ to denote the Euclidean norm.
The following matrices will feature in asymptotic expansions of our estimators:
\begin{align*}
\Gamma_{h,\mathsf{x}} &= \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(u)^\prime K\left(u\right) g(\mathsf{x} + hv)\mathrm{d} u = \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(u)^\prime K\left(u\right) g_u\mathrm{d} u + O(h) = \Gamma_{1h,\mathsf{x}} + O(h),
\end{align*}
and
\begin{align*}
\Sigma_{h,\mathsf{x}} &= \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(v)^\prime\Big[ F(\mathsf{x}+h(u\wedge v)) - F(\mathsf{x}+hu)F(\mathsf{x}+hv) \Big] K\left(u\right)K\left(v\right) g(\mathsf{x}+hu)g(\mathsf{x}+hv)\mathrm{d} u\mathrm{d} v\\
&= F(\mathsf{x})(1-F(\mathsf{x}))\left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u) K(u) g_u\mathrm{d} u\right)\left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u) K(u) g_u\mathrm{d} u\right)^\prime\\
&\quad + h\iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(v)^\prime K(u)K(v)\Big[ - F(\mathsf{x})(uf_u+vf_v)g_ug_v + F(\mathsf{x})(1-F(\mathsf{x}))(u\dot{g}_ug_v+v\dot{g}_vg_u)\Big] \mathrm{d} u\mathrm{d} v\\
&\quad + h\iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(v)^\prime K(u)K(v)(u\wedge v)f_{u\wedge v}g_ug_v \mathrm{d} u\mathrm{d} v + O(h^2)\\
&= \Sigma_{1h,\mathsf{x}} + h\Sigma_{2h,\mathsf{x}} + O(h^2).
\end{align*}
Now we list the main assumptions.
\begin{assumption}\label{assumption:dgp pointwise}
$x_{1},\dots,x_{n}$ is a random sample from a distribution $F(\cdot)$ supported on $\mathcal{X}\subseteq\mathbb{R}$, and $\mathsf{x}\in\mathcal{X}$.
(i) For some $\delta>0$, $F(\cdot)$ is absolutely continuous on $[\mathsf{x}-\delta,\mathsf{x}+\delta]$ with a density $f(\cdot)$ admitting constants $f(\mathsf{x}-)$, $f(\mathsf{x}+)$, $\dot{f}(\mathsf{x}-)$, and $\dot{f}(\mathsf{x}+)$, such that
\begin{align*}
\sup_{u\in[-\delta,0)} \frac{f(\mathsf{x}+u) - f(\mathsf{x}-) - u\dot{f}(\mathsf{x}-)}{u^2} + \sup_{u\in(0,\delta]} \frac{f(\mathsf{x}+u) - f(\mathsf{x}+) - u\dot{f}(\mathsf{x}+)}{u^2} < \infty.
\end{align*}
(ii) $K(\cdot)$ is nonnegative, symmetric, and continuous on its support $[-1,1]$, and integrates to 1.
(iii) $R(\cdot)$ is locally bounded, and there exists a positive-definite diagonal matrix $\Upsilon_{h}$ for each $h>0$, such that $\Upsilon_{h}R(u)=R(u/h)$
(iv) For all $h$ sufficiently small, the minimum eigenvalues of $\Gamma_{h,\mathsf{x}}$ and $h^{-1}\Sigma_{h,\mathsf{x}}$ are bounded away from zero.
\qed
\end{assumption}
\begin{assumption}\label{assumption:design pointwise}
For some $\delta>0$, $G(\cdot)$ is absolutely continuous on $[\mathsf{x}-\delta,\mathsf{x}+\delta]$ with a derivative $g(\cdot)\geq 0$ admitting constants $g(\mathsf{x}-)$, $g(\mathsf{x}+)$, $\dot{g}(\mathsf{x}-)$, and $\dot{g}(\mathsf{x}+)$, such that
\begin{align*}
\sup_{u\in[-\delta,0)} \frac{g(\mathsf{x}+u) - g(\mathsf{x}-) - u\dot{g}(\mathsf{x}-)}{u^2} + \sup_{u\in(0,\delta]} \frac{g(\mathsf{x}+u) - g(\mathsf{x}+) - u\dot{g}(\mathsf{x}+)}{u^2} < \infty.
\end{align*}
\vskip-2em\qed
\end{assumption}
\begin{example}[Local Polynomial Estimator]\label{example:local polynomial density estimator}
Before closing this section, we briefly introduce the local polynomial estimator of \cite*{Cattaneo-Jansson-Ma_2020_JASA}, which is a special case of our local regression distribution estimator. The local polynomial estimator employs the following polynomial basis:
\[ R(u) = \Big(1,\ u,\ \frac{1}{2}u^2,\ \cdots,\ \frac{1}{p!}u^p\Big)^\prime,\]
for some $p\in \mathbb{N}$. As a result, it estimates the distribution function, the density function, and derivatives thereof. To be precise,
\[\theta(\mathsf{x}) = \Big( F(\mathsf{x}),\ f(\mathsf{x}),\ f^{(1)}(\mathsf{x}),\ \cdots,\ f^{(p-1)}(\mathsf{x}) \Big)^\prime.\]
With $R(\cdot)$ being a polynomial basis, it is straightforward to characterize the approximation bias $\varrho(h,\mathsf{x})$. Assuming the distribution function $F(\cdot)$ is at least $p+1$ times continuously differentiable in a neighborhood of $\mathsf{x}$, one can employ a Taylor expansion argument and show that $\varrho(h,\mathsf{x}) = O(h^{p+1})$. We will revisit this local polynomial estimator below as a leading example when we discuss pointwise and uniform asymptotic properties of our local distribution estimator.
\qed
\end{example}
\section{Pointwise Distribution Theory}\label{section:pointwise distribution theory}
We discuss pointwise (i.e., for a fixed evaluation point $\mathsf{x}\in\mathcal{X}$) large-sample properties of the local $L^2$ estimator \eqref{eq:local L2 estimator}, and that of the local regression estimator \eqref{eq:local regression estimator}. For ease of exposition, we suppress the dependence on the evaluation point $\mathsf{x}$ whenever possible.
\subsection{Local $L^2$ Distribution Estimation}\label{subsection:local projection: pointwise distribution theory}
With simple algebra, the local $L^2$ estimator in \eqref{eq:local L2 estimator} takes the following form
\begin{align*}
\hat\theta_G &= \left( \int_{\mathcal{X}} R(u-\mathsf{x})R(u-\mathsf{x})^\prime \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} G(u) \right)^{-1}\left( \int_{\mathcal{X}} R(u-\mathsf{x})\hat{F}(u) \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} G(u) \right).
\end{align*}
We can further simplify the above. First note that the ``denominator'' can be rewritten as
\begin{align*}
&\ \int_{\mathcal{X}} R(u-\mathsf{x})R(u-\mathsf{x})^\prime \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} G(u)\\
&= \Upsilon_h^{-1}\left(\int_{\mathcal{X}} \Upsilon_h R(u-\mathsf{x})R(u-\mathsf{x})^\prime \Upsilon_h \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) g(u)\mathrm{d} u\right) \Upsilon_h^{-1}= \Upsilon_h^{-1}\Gamma_h \Upsilon_h^{-1}.
\end{align*}
The same technique can be applied to the ``numerator'', which leads to
\begin{align}
\nonumber &\ \hat\theta_G - \theta = \Upsilon_h\Gamma_h^{-1} \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\hat{F}(\mathsf{x} + hu) K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u\right) - \theta\\
\label{eq:approximation bias, projection}=& \Upsilon_h\Gamma_h^{-1}\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[F(\mathsf{x} + hu) - \theta^\prime R(u)\Upsilon_h^{-1}\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u\\
\label{eq:linear variance, projection} &\qquad +\Upsilon_h \frac{1}{n}\sum_{i=1}^n\Gamma_h^{-1} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u.
\end{align}
The above provides a further expansion of the local $L^2$ estimator into a term that contributes as bias, and another term that contributes asymptotically to the variance.
The large-sample properties of the local $L^2$ estimator \eqref{eq:local L2 estimator} are as follows.
\begin{thm}[Local $L^2$ Distribution Estimation: Asymptotic Normality]\label{thm:local projection: asymptotic normality}
Assume Assumptions \ref{assumption:dgp pointwise} and \ref{assumption:design pointwise} hold, and that $h\to 0$, $nh\to \infty$ and $n\varrho(h)^2/h\to 0$. Then\\
(i) \eqref{eq:approximation bias, projection} satisfies
\begin{align*}
\left|\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[F(\mathsf{x} + hu) - \theta^\prime R(u)\Upsilon_h^{-1}\Big] K\left(u\right) \mathrm{d} u\right| = O(\varrho(h)).
\end{align*}
(ii) \eqref{eq:linear variance, projection} satisfies
\begin{align*}
\mathbb{V}\left[ \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u \right] &= \Sigma_h,
\end{align*}
and
\begin{align*}
\Sigma_h^{-1/2} \left(\frac{1}{\sqrt{n}}\sum_{i=1}^n\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u\right) \rightsquigarrow \mathcal{N}(0,I).
\end{align*}
(iii) The local $L^2$ estimator is asymptotically normally distributed
\begin{align*}
\sqrt{n}\left(\Gamma_h^{-1}\Sigma_h\Gamma_h^{-1}\right)^{-1/2}\Upsilon_h^{-1} (\hat\theta_G - \theta) \rightsquigarrow \mathcal{N}(0, I).
\end{align*}
\vskip-2em\qed
\end{thm}
For valid inference, one needs to construct standard errors. To start, note that $\Gamma_{h}$ is known, and hence we only need to estimate $\Sigma_{h}$. Consider the following:
\begin{align}
\nonumber\hat\Sigma_{h} &= \frac{1}{n}\sum_{i=1}^n\iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(v)^\prime\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - \hat{F}(\mathsf{x} + hu)\Big]\Big[\mathds{1}(x_i\leq \mathsf{x} + hv) - \hat{F}(\mathsf{x} + hv)\Big] \\
\label{eq:local projection: standard error pointwise}&\qquad \qquad \qquad \qquad K(u)K(v) g(\mathsf{x}+hu)g(\mathsf{x}+hv)\mathrm{d} u\mathrm{d} v,
\end{align}
where $\hat{F}(\cdot)$ is the empirical distribution function. The following theorem shows that standard errors constructed using $\hat{\Sigma}_h$ are consistent.
\begin{thm}[Local $L^2$ Distribution Estimation: Standard Errors]\label{thm:local projection: standard error pointwise}
Assume Assumptions \ref{assumption:dgp pointwise} and \ref{assumption:design pointwise} hold, and that $h\to 0$ and $nh\to \infty$. Let $c$ be a nonzero vector of suitable dimension, then
\begin{align*}
\left|\frac{c^\prime \hat{\Sigma}_h c}{c^\prime \Sigma_{h} c} - 1\right| &= O_{\mathbb{P}}\left( \sqrt{\frac{1}{nh}} \right).
\end{align*}
If, in addition that $n\varrho(h)^2/h\to 0$, then
\begin{align*}
\frac{c^\prime(\hat\theta_G - \theta)}{\sqrt{c^\prime \Upsilon_{h}\Gamma_h^{-1}\hat\Sigma_{h} \Gamma_h^{-1}\Upsilon_{h} c/n}} \rightsquigarrow \mathcal{N}(0,1).
\end{align*}
\vskip-2em\qed
\end{thm}
\subsection{Local Regression Distribution Estimation}\label{subsection:local regression: pointwise distribution theory}
The local regression distribution estimator \eqref{eq:local regression estimator} can be understood as a special case of the local $L^2$ estimator by setting $G=\hat{F}$ (i.e., using the empirical distribution as the design). However, the empirical measure $\hat{F}$ is not smooth, so that large-sample properties of the local regression estimator cannot be deduced directly from Theorem \ref{thm:local projection: asymptotic normality}. In this subsection, we will show that estimates obtained by the two approaches, \eqref{eq:local L2 estimator} and \eqref{eq:local regression estimator}, are asymptotically equivalent under suitable regularity conditions. To be precise, we establish the equivalence of the local regression distribution estimator, $\hat\theta$, and the (infeasible) local $L^2$ distribution estimator, $\hat{\theta}_F$ (i.e., using $F$ as the design weighting in \eqref{eq:local L2 estimator}). As before, we suppress the dependence on the evaluation point $\mathsf{x}$.
First, the local regression estimator can be written as
\begin{align*}
\hat\theta - \theta &= \left( \frac{1}{n}\sum_{i=1}^n R(x_i-\mathsf{x})R(x_i-\mathsf{x})^\prime \frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right)\right)^{-1}\\
&\qquad \qquad \left( \frac{1}{n}\sum_{i=1}^n R(x_i-\mathsf{x})\Big[ \hat{F}(x_i) -R(x_i-\mathsf{x})^\prime\theta \Big] \frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right)\right)\\
&= \Upsilon_h\hat{\Gamma}_h^{-1}\Gamma_h\Gamma_h^{-1}\left( \frac{1}{n}\sum_{i=1}^n \Upsilon_hR(x_i-\mathsf{x})\Big[ \hat{F}(x_i) -R(x_i-\mathsf{x})^\prime\theta \Big] \frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right)\right),
\end{align*}
where
\begin{align*}
\hat{\Gamma}_h &= \frac{1}{n}\sum_{i=1}^n \Upsilon_hR(x_i-\mathsf{x})R(x_i-\mathsf{x})^\prime\Upsilon_h \frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right),
\end{align*}
and $\Gamma_h$ is defined as before with $G=F$.
To proceed, we further expand as follows
\begin{align}
\nonumber&\ \frac{1}{n}\sum_{i=1}^n \Upsilon_hR(x_i-\mathsf{x})\Big[ \hat{F}(x_i) -R(x_i-\mathsf{x})^\prime\theta \Big] \frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right)\\
\nonumber=&\ \frac{1}{n^2}\sum_{i,j=1,i\neq j}^n \Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -F(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right)\\
\label{eq:leave-in bias}&\ + \frac{1}{n^2}\sum_{j=1}^n \Upsilon_hR(x_j-\mathsf{x})\Big[ 1 -F(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right)\\
\label{eq:approximation bias}&\ + \frac{1}{n}\sum_{j=1}^n \Upsilon_hR(x_j-\mathsf{x})\Big[ F(x_j) -R(x_j-\mathsf{x})^\prime\theta \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right).
\end{align}
The last two terms correspond to the leave-in bias and the approximation bias, respectively. We further decompose the first term with conditional expectation:
\begin{align}
\nonumber&\ \frac{1}{n^2}\sum_{i,j=1,i\neq j}^n \Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -F(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right)\\
\nonumber=&\ \frac{1}{n^2}\sum_{i,j=1,i\neq j}^n \mathbb{E}\left[\left.\Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -F(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right)\right| x_i\right]\\
\nonumber&\ + \frac{1}{n^2}\sum_{i,j=1,i\neq j}^n \Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -F(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right) \\
\nonumber&\qquad \qquad \qquad - \mathbb{E}\left[\left.\Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -F(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right)\right| x_i\right]\\
\label{eq:linear variance}=&\ \frac{n-1}{n^2}\sum_{i=1}^n \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) f(\mathsf{x} + hu) \mathrm{d} u\\
\nonumber&\ + \frac{1}{n^2}\sum_{i,j=1,i\neq j}^n \Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -F(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right) \\
\label{eq:quadratic variance}&\qquad \qquad \qquad- \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) f(\mathsf{x} + hu)\mathrm{d} u.
\end{align}
The following theorem studies the large-sample properties of each term in the above decomposition, and shows that the local regression distribution estimator is asymptotically equivalent to the local $L^2$ estimator by setting $G=F$, and hence it is asymptotically normally distributed.
\begin{thm}[Local Regression Distribution Estimation: Asymptotic Normality]\label{thm:local regression: asymptotic normality}
Assume Assumption \ref{assumption:dgp pointwise} holds, and that $h\to 0$, $nh^2\to \infty$ and $n\varrho(h)^2/h\to 0$. Then\\
(i) $\hat{\Gamma}_h$ satisfies
\begin{align*}
\left|\hat{\Gamma}_h - \Gamma_h\right| = O_{\mathbb{P}}\left( \sqrt{\frac{1}{nh}} \right).
\end{align*}
(ii) \eqref{eq:leave-in bias} and \eqref{eq:approximation bias} satisfy
\begin{align*}
\text{\eqref{eq:leave-in bias}} &= O_{\mathbb{P}}\left(\frac{1}{n}\right),\qquad \text{\eqref{eq:approximation bias}}= O_{\mathbb{P}}(\varrho(h)).
\end{align*}
(iii) \eqref{eq:quadratic variance} satisfies
\begin{align*}
\text{\eqref{eq:quadratic variance}} &= O_{\mathbb{P}}\left( \sqrt{\frac{1}{n^2h}} \right).
\end{align*}
(iv) The local regression distribution estimator \eqref{eq:local regression estimator} satisfies
\begin{align*}
\sqrt{n}\left(\Gamma_h^{-1}\Sigma_h\Gamma_h^{-1}\right)^{-1/2}\Upsilon_h^{-1} (\hat\theta - \theta) &= \sqrt{n}\left(\Gamma_h^{-1}\Sigma_h\Gamma_h^{-1}\right)^{-1/2}\Upsilon_h^{-1} (\hat\theta_F - \theta) + o_{\mathbb{P}}(1)
\rightsquigarrow \mathcal{N}(0,I).
\end{align*}
\vskip-2.5em\qed
\end{thm}
We now discuss how to construct standard errors in the local regression framework. Note that $\Gamma_{h}$ can be estimated by $\hat{\Gamma}_h$, whose properties have already been studied in Theorem \ref{thm:local regression: asymptotic normality}(i). To estimate $\Sigma_{h}$, we propose the following
\begin{align*}
\hat\Sigma_{h} &= \frac{1}{n}\sum_{i=1}^n \left[\frac{1}{n}\sum_{j=1}^n \Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -\hat{F}(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right)\right]\\
&\qquad \qquad \qquad \cdot\left[\frac{1}{n}\sum_{j=1}^n \Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -\hat{F}(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right)\right]^\prime.
\end{align*}
where $\hat{F}(\cdot)$ is the empirical distribution function. The following theorem shows that standard errors constructed using $\hat{\Sigma}_h$ are consistent.
\begin{thm}[Local Regression Distribution Estimation: Standard Errors]\label{thm:local regression: standard error pointwise}
Assume Assumption \ref{assumption:dgp pointwise} holds. In addition, assume $h\to 0$ and $nh^2\to \infty$. Let $c$ be a nonzero vector of suitable dimension. Then
\begin{align*}
\left|\frac{c^\prime \hat\Gamma_h^{-1} \hat{\Sigma}_h \hat\Gamma_h^{-1} c}{c^\prime \Gamma_h^{-1} \Sigma_{h} \Gamma_h^{-1} c} - 1\right| &= O_{\mathbb{P}}\left( \sqrt{\frac{1}{nh^2}} \right).
\end{align*}
If, in addition that $n\varrho(h)^2/h\to 0$, one has
\begin{align*}
\frac{c^\prime(\hat\theta - \theta)}{\sqrt{c^\prime \Upsilon_{h}\hat{\Gamma}_h^{-1}\hat\Sigma_{h} \hat{\Gamma}_h^{-1}\Upsilon_{h} c/n}} \rightsquigarrow \mathcal{N}(0,1).
\end{align*}
\vskip-2em\qed
\end{thm}
\section{Efficiency}\label{section:Efficiency}
For ease of presentation, we focus on the (infeasible) local $L^2$ distribution estimator $\hat\theta_F$,
\begin{align}\label{eq:local polynomial projection estimator}
\hat\theta_{F} &= \operatorname*{argmin}_{\theta} \int_{\mathcal{X}}\left( \hat{F}(u) - R(u-\mathsf{x})^\prime\theta \right)^2 \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right)\mathrm{d} F(u),
\end{align}
but all the results in this section are applicable to the local regression distribution estimator $\hat\theta$, as we showed earlier that it is asymptotically equivalent to $\hat\theta_F$. In addition, we consider a specific basis:
\begin{align}\label{eq:polynomial basis}
R(u) &= \left(1,\ P(u)^\prime,\ Q(u)\right)^\prime,
\end{align}
where $P(u)$ is a polynomial basis of order $p$:
\begin{align*}
P(u) = \left(u,\ \frac{1}{2}u^2,\ \cdots,\ \frac{1}{p!}u^p\right)^\prime,
\end{align*}
and $Q(u)$ is a scalar function, and hence is a ``redundant regressor.'' Without $Q(\cdot)$, the above reduces to the local polynomial estimator of \cite*{Cattaneo-Jansson-Ma_2020_JASA}. See Section \ref{section:setup} and Example \ref{example:local polynomial density estimator} for an introduction.
We consider additional regressors because they may help improve efficiency (i.e., reduce the asymptotic variance). Following Assumption \ref{assumption:dgp pointwise}, we assume there exists a scalar $\upsilon_h$ (depending on $h$) such that $\upsilon_hQ(u)=Q(u/h)$. Therefore, $\Upsilon_h$ is a diagonal matrix containing $1,h^{-1},h^{-2},\cdots,h^{-p},\upsilon_h$. As we consider a (local) polynomial basis, it is natural to impose smoothness assumptions on $F(\cdot)$. In particular,
\begin{assumption}\label{assumption:smoothness}
For some $\delta>0$, $F(\cdot)$ is $(p+1)$-times continuously differentiable in $\mathcal{X}\cap[\mathsf{x}-\delta,\mathsf{x}+\delta]$ for some $p\geq 1$, and $G(\cdot)$ is twice continuously differentiable in $\mathcal{X}\cap[\mathsf{x}-\delta,\mathsf{x}+\delta]$.
\qed
\end{assumption}
Under the above assumption, the approximation error satisfies $\varrho(h)=O(h^{p+1})$, and the parameter $\theta$ can be partitioned into the following:
\begin{align*}
\theta = \Big(\theta_{1},\ \theta_{P}^\prime,\ \theta_{Q} \Big)^\prime = \Big( F(\mathsf{x}),\quad f(\mathsf{x}),\ \cdots,\ f^{(p-1)}(\mathsf{x}),\quad 0 \Big)^\prime.
\end{align*}
We first state a corollary, which specializes Theorem \ref{thm:local projection: asymptotic normality} to the polynomial basis \eqref{eq:polynomial basis}.
\begin{coro}[Local Polynomial $L^2$ Distribution Estimation: Asymptotic Normality]\label{coro:asy normal loc pol projection estimator}
Assume Assumptions \ref{assumption:dgp pointwise} and \ref{assumption:smoothness} hold, and that $h\to 0$, $nh\to \infty$, and $n\varrho(h)^2/h\to 0$. Then the local polynomial $L^2$ distribution estimator in \eqref{eq:local polynomial projection estimator} satisfies
\begin{align*}
\sqrt{n}\left(\Gamma_h^{-1}\Sigma_h\Gamma_h^{-1}\right)^{-1/2}\Upsilon_h^{-1} (\hat\theta_F - \theta) \rightsquigarrow \mathcal{N}(0, I).
\end{align*}
\vskip-2em\qed
\end{coro}
\subsection{Effect of Orthogonalization}
To start, consider the following (sequentially) orthogonalized basis:
\begin{align}\label{eq:orthogonalized polynomial basis}
R^\perp(u) &= \left(1,\ P^\perp(u)^\prime,\ Q^\perp(u)\right)^\prime,
\end{align}
where
\begin{align*}
P^\perp(u) &= P^\perp(u) - \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u) P(u) \mathrm{d} u,\\
Q^\perp(u) &= Q(u) - \Big(1, P(v)^\prime\Big)\left( \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(v)\Big(1, P(v)^\prime\Big)^\prime\Big(1, P(v)^\prime\Big) \mathrm{d} v \right)^{-1}\left( \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(v)\Big(1, P(v)^\prime\Big)^\prime Q(v) \mathrm{d} v \right).
\end{align*}
The above transformation can be represented by the following:
\begin{align*}
R^\perp(u) &= \Lambda_h^\prime R(u),
\end{align*}
where $\Lambda_h$ is a nonsingular upper triangular matrix. (Note that the matrix $\Lambda_h$ depends on the bandwidth only because we would like to handle both interior and boundary evaluation points. If, for example, we fix the evaluation point to be in the interior of the support of the data, then $\Lambda_h$ is a fixed matrix and no longer depends on $h$. Alternatively, one could also use the notation ``$\Lambda_\mathsf{x}$'' to denote such dependence.) Now consider the following orthogonalized local polynomial $L^2$ estimator
\begin{align}\label{eq:local orthogonal polynomial projection estimator}
\hat\theta_{F}^\perp &= \operatorname*{argmin}_{\theta} \int_{\mathcal{X}}\left( \hat{F}(u) - \Lambda_h^\prime R(u-\mathsf{x})^\prime\theta \right)^2 \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right)\mathrm{d} F(u).
\end{align}
To discuss its properties, we partition the estimator and the target parameter as
\begin{align*}
\hat\theta_F^\perp &= \Big( \hat\theta_{1,F}^\perp,\ (\hat\theta_{P,F}^\perp)^\prime,\ \hat\theta_{Q,F}^\perp \Big)^\prime,
\end{align*}
where $\hat\theta_{1,F}^\perp$ is the first element of $\hat\theta_F^\perp$ and $\hat\theta_{Q,F}^\perp$ is the last element of $\hat\theta_F^\perp$. Similarly, we can partition the target parameter,
\begin{align*}
\theta^\perp &= \Lambda_h^{-1}\theta = \Big( \theta^\perp_{1},\ (\theta^\perp_{P})^\prime,\ \theta^\perp_{Q} \Big)^\prime,
\end{align*}
so that $\theta^\perp_{1}$ is the first element of $\Lambda_h^{-1}\theta$ and $\theta^\perp_{Q}$ is the last element of $\Lambda_h^{-1}\theta$. As $\theta_{Q}=0$, simple least squares algebra implies
\begin{align*}
\theta^\perp &= \Big( \theta^\perp_{1},\ \theta_{P}^\prime,\ 0 \Big)^\prime = \Big( \theta^\perp_{1},\ f(\mathsf{x}),\ f^{(1)}(\mathsf{x}),\ \cdots,\ f^{(p-1)}(\mathsf{x}),\ 0 \Big)^\prime.
\end{align*}
Note that, in general, $\theta^\perp_{1}\neq \theta_{1}$, meaning that after orthogonalization, the intercept of the local polynomial estimator no longer estimates the distribution function $F(\mathsf{x})$.
The following corollary gives the large-sample properties of the orthogonalized local polynomial estimator, excluding the intercept.
\begin{coro}[Orthogonalized Local Polynomial $L^2$ Distribution Estimation: Asymptotic Normality]\label{coro:asy normal ortho loc pol projection estimator}
Assume Assumptions \ref{assumption:dgp pointwise} and \ref{assumption:smoothness} hold, and that $h\to 0$, $nh\to \infty$, and $n\varrho(h)^2/h\to 0$. Then the orthogonalized local polynomial $L^2$ distribution estimator in \eqref{eq:local orthogonal polynomial projection estimator} satisfies
\begin{align*}
\begin{bmatrix}
(\Gamma_{P,h}^\perp)^{-1} \Sigma_{PP,h}^\perp (\Gamma_{P,h}^\perp)^{-1} & (\Gamma_{P,h}^\perp)^{-1} \Sigma_{PQ,h}^\perp (\Gamma_{Q,h}^\perp)^{-1} \\
(\Gamma_{Q,h}^\perp)^{-1} \Sigma_{QP,h}^\perp (\Gamma_{P,h}^\perp)^{-1} & (\Gamma_{Q,h}^\perp)^{-1} \Sigma_{QQ,h}^\perp (\Gamma_{Q,h}^\perp)^{-1}
\end{bmatrix}^{-1/2} \sqrt{\frac{n}{hf(\mathsf{x})}}\Upsilon_{-1, h}^{-1}\begin{bmatrix}
\hat\theta^\perp_{P,F} - \theta_{P} \\
\hat\theta^\perp_{Q,F}
\end{bmatrix} \rightsquigarrow \mathcal{N}(0, I),
\end{align*}
where
\begin{align*}
&\ \Gamma_{P,h}^\perp = \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} P^\perp(u)P^\perp(u)^\prime K(u)\mathrm{d} u,
\quad \Gamma_{Q,h}^\perp = \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} Q^\perp(u)^2 K(u)\mathrm{d} u,\\
&\ \Sigma_{PP,h}^\perp = \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)K(v)P^\perp(u)P^\perp(v)^\prime (u\wedge v) \mathrm{d} u\mathrm{d} v,\\
&\ \Sigma_{QQ,h}^\perp = \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)K(v)Q^\perp(u)Q^\perp(v) (u\wedge v) \mathrm{d} u\mathrm{d} v,\\
&\ \Sigma_{PQ,h}^\perp = (\Sigma_{QP,h}^\perp)^\prime = \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)K(v)P^\perp(u)Q^\perp(v) (u\wedge v) \mathrm{d} u\mathrm{d} v,
\end{align*}
and $\Upsilon_{-1, h}$ is a diagonal matrix containing $h^{-1},h^{-2},\cdots,h^{-p},\upsilon_h$.
\qed
\end{coro}
\subsection{Optimal $Q$}
Now we discuss the optimal choice of $Q$, which minimizes the asymptotic variance of the minimum distance estimator. Recall from the main paper that, with orthogonalized basis, the minimum distance estimator of ${f}^{(\ell)}(\mathsf{x})$, for $0\leq \ell\leq p-1$, has an asymptotic variance
\begin{align*}
f(\mathsf{x})\Big[e_\ell^\prime(\Gamma_{P,h}^\perp)^{-1} \Sigma_{PP,h}^\perp (\Gamma_{P,h}^\perp)^{-1}e_\ell - e_\ell^\prime(\Gamma_{P,h}^\perp)^{-1} \Sigma_{PQ,h}^\perp
(\Sigma_{QQ,h}^{\perp} )^{-1}
\Sigma_{QP,h}^\perp (\Gamma_{P,h}^\perp)^{-1}e_\ell\Big],
\end{align*}
where $e_\ell$ is the $(\ell+1)$-th standard basis vector. In subsequent analysis, we drop the multiplicative factor $f(\mathsf{x})$.
Let $p_\ell(u)$ be defined as
\begin{align*}
p_\ell(u) &= e_\ell^\prime(\Gamma_{P,h}^\perp)^{-1} P^\perp(u),
\end{align*}
then the objective is to maximize
\begin{align*}
\left( \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)K(v)Q^\perp(u)Q^\perp(v)(u\wedge v) \mathrm{d} u\mathrm{d} v \right)^{-1}\left( \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)K(v)p_\ell(u)Q^\perp(v)(u\wedge v) \mathrm{d} u\mathrm{d} v \right)^2.
\end{align*}
Alternatively, we would like to solve (recall that $Q(u)$ is a scaler function)
\begin{align*}
\text{maximize}&\ \frac{\left(\iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)K(v)p_\ell(u)q(v)(u\wedge v) \mathrm{d} u\mathrm{d} v\right)^2}{\iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)K(v)q(u)q(v)(u\wedge v) \mathrm{d} u\mathrm{d} v},\quad
\text{subject to}\ \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)q(u)(1, P(u)^\prime)\mathrm{d} u = 0.
\end{align*}
To proceed, define the following transformation for a function $g(\cdot)$:
\begin{align*}
\mathcal{H}(g)(u) &= \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} \mathds{1}(v\geq u)K(v)g(v)\mathrm{d} v.
\end{align*}
This transformation satisfies two important properties, which are summarized in the following lemma.
\begin{lem}[$\mathcal{H}$-transformation]\label{lem:survival transformation}\
(i) If $g_1(\cdot)$ and $g_2(\cdot)$ are bounded, and that either $\int_{\frac{\mathcal{X}-\mathsf{x}}{h}}K(u)g_1(u)\mathrm{d} u$ or $\int_{\frac{\mathcal{X}-\mathsf{x}}{h}}K(u)g_2(u)\mathrm{d} u$ is zero, then
\begin{align*}
\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \mathcal{H}(g_1)(u)\mathcal{H}(g_2)(u) \mathrm{d} u &= \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)K(v)g_1(u)g_1(v)(u\wedge v)\mathrm{d} u\mathrm{d} v.
\end{align*}
(ii) If $g_1(\cdot)$ and $g_2(\cdot)$ are bounded, $g_2(\cdot)$ is continuously differentiable with a bounded derivative, and that either $\int_{\frac{\mathcal{X}-\mathsf{x}}{h}}K(u)g_1(u)\mathrm{d} u$ or $\int_{\frac{\mathcal{X}-\mathsf{x}}{h}}K(u)g_2(u)\mathrm{d} u$ is zero, then
\begin{align*}
\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \mathcal{H}(g_1)(u)\dot{g}_2(u) \mathrm{d} u &= \int_{\frac{\mathcal{X}-\mathsf{x}}{h}}K(u)g_1(u) g_2(u) \mathrm{d} u.
\end{align*}
\vskip-2.3em\qed
\end{lem}
With the previous lemma, we can rewrite the maximization problem as
\begin{align}
\nonumber\text{maximize}&\ \qquad \frac{\left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \mathcal{H}(p_\ell)(u)\mathcal{H}(q)(u)\mathrm{d} u\right)^2}{\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \mathcal{H}(q)(u)^2 \mathrm{d} u}\\
\label{eq:maximization}\text{subject to}&\ \qquad \int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \dot{P}(u)\mathcal{H}(q)(u)\mathrm{d} u = 0,\qquad \mathcal{H}(q)\left( \frac{\inf \mathcal{X}-\mathsf{x}}{h}\vee (-1) \right) =0.
\end{align}
\begin{thm}[Variance Bound of the Minimum Distance Estimator]\label{thm:variance bound}
An upper bound of the maximization problem \eqref{eq:maximization} is
\begin{align*}
e_\ell^\prime(\Gamma_{P,h}^\perp)^{-1} \Sigma_{PP,h}^\perp (\Gamma_{P,h}^\perp)^{-1}e_\ell - e_\ell^\prime \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \dot{P}(u)\dot{P}(u)^\prime \mathrm{d} u\right)^{-1} e_{\ell}.
\end{align*}
Therefore, the asymptotic variance of the minimum distance estimator is bounded below by
\begin{align*}
f(\mathsf{x})e_\ell^\prime \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \dot{P}(u)\dot{P}(u)^\prime \mathrm{d} u\right)^{-1} e_{\ell},
\end{align*}
where $\dot{P}(u)=(1,\ u,\ u^2/2,\ u^3/3!,\ \cdots,\ u^{p-1}/(p-1)!)^\prime$.
\qed
\end{thm}
\begin{example}[Local Linear/Quadratic Minimum Distance Density Estimation]\label{example:p=1 v=1 interior}
Consider a simple example where $\ell=0$ and $P(u)=u$, which means we focus on the asymptotic variance of the estimated density in a local linear regression. Also assume we employ a uniform kernel: $K(u)=\frac{1}{2}\mathds{1}(|u|\leq 1)$, and that the integration region is $\frac{\mathcal{X}-\mathsf{x}}{h}=\mathbb{R}$ (i.e., $\mathsf{x}$ is an interior evaluation point). Note that this example also applies to local quadratic regressions, as $u$ and $u^2$ are orthogonal for interior evaluation points.
Taking $P(u) = u$, the variance bound in Theorem \ref{thm:variance bound} is easily found to be
\begin{align*}
f(\mathsf{x})\left(\int_{-1}^1 \dot{P}(u)\dot{P}(u)^\prime \mathrm{d} u\right)^{-1} = f(\mathsf{x})\frac{1}{2}.
\end{align*}
We now calculate the asymptotic variance of the minimum distance estimator. To be specific, we choose $Q(u) = u^{2j+1}$, which is a higher-order polynomial function. With tedious calculation, one can show that the minimum distance estimator has the following asymptotic variance
\begin{align*}
\mathrm{Asy}\mathbb{V}[\hat{f}_{\mathtt{MD}}(\mathsf{x})] &= f(\mathsf{x})\frac{11+4j}{20+8j},
\end{align*}
which asymptotes to $f(\mathsf{x})/2$ as $j\to \infty$. As a result, it is possible to achieve the maximum amount of efficiency gain by including one higher-order polynomial and using our minimum distance estimator.
In Figure \ref{fig:equi kerenl p=1 v=1}, we plot the equivalent kernel of the local linear minimum distance density estimator using a uniform kernel. Without the redundant regressor, it is equivalent to the kernel density estimator using the Epanechnikov kernel. As $j$ gets larger, however, the equivalent kernel of the minimum distance estimator becomes closer to the uniform kernel, which is why, as $j\to \infty$, the minimum distance estimator has an asymptotic variance the same as the kernel density estimator using the uniform kernel.
\qed
\end{example}
\begin{figure}[!t]
\centering
\includegraphics[width=0.45\textwidth]{SA-Tables-Figures/densityP1uniform.pdf}
\caption{Equivalent Kernel of the Local Linear Minimum Distance Density Estimator.} \label{fig:equi kerenl p=1 v=1}
\begin{flushleft}
\footnotesize
\textit{Notes}: The basis function $R(u)$ consists of an intercept, a linear term $u$ (i.e., local linear regression), and an odd higher-order polynomial term $u^{2j+1}$ for $j=1, 2, \cdots,30$. Without the higher-order polynomial regressor, the local linear density estimator using the uniform kernel is equivalent to the kernel density estimator using the Epanechnikov kernel (black line). Including a higher-order redundant regressor leads to an equivalent kernel that approaches the uniform kernel as $j$ tends to infinity (red).
\end{flushleft}
\end{figure}
\begin{example}[Local Cubic Minimum Distance Estimation]\label{example:p=3 v=1 interior}
We adopt the same setting in Example \ref{example:p=1 v=1 interior}, i.e., local polynomial density estimation with the uniform kernel at an interior evaluation point. The difference is that we now consider a local cubic regression: $P(u) = (u, \frac{1}{2}u^2, \frac{1}{3!}u^3)^\prime$.
As before, the variance bound in Theorem \ref{thm:variance bound} is easily found to be
\begin{align*}
f(\mathsf{x})\left(\int_{-1}^1 \dot{P}(u)\dot{P}(u)^\prime \mathrm{d} u\right)^{-1} = f(\mathsf{x})
\begin{bmatrix}
\frac{9}{8} & 0 & -\frac{15}{4} \\
0 & \frac{3}{2} & 0 \\
-\frac{15}{4} & 0 & \frac{45}{2} \\
\end{bmatrix}.
\end{align*}
Again, we compute the asymptotic variance of our minimum distance estimator. Note, however, that now we have both odd and even order polynomials in our basis $P(u)$, therefore we include two higher-order polynomials, that is, we set $Q(u)=(u^{2j}, u^{2j+1})^\prime$. The asymptotic variance of our minimum distance estimator is
\begin{align*}
\mathrm{Asy}\mathbb{V}\begin{bmatrix}
\hat{f}_{\mathtt{MD}}(\mathsf{x})\\
\hat{f}_{\mathtt{MD}}^{(1)}(\mathsf{x})\\
\hat{f}_{\mathtt{MD}}^{(2)}(\mathsf{x})
\end{bmatrix} &= f(\mathsf{x})\begin{bmatrix}
\frac{9 (4 j+15)}{16 (2 j+7)} & 0 & -\frac{15 (4 j+17)}{8 (2 j+7)} \\
0 & \frac{12 j+39}{8 j+20} & 0 \\
-\frac{15 (4 j+17)}{8 (2 j+7)} & 0 & \frac{45 (4 j+19)}{8 j+28} \\
\end{bmatrix},
\end{align*}
which, again, asymptotes to the variance bound as $j\to \infty$. See also Table \ref{table:density variance} for the efficiency gain of employing the minimum distance technique.
\qed
\end{example}
\begin{table}[!t]
\centering
\caption{Variance Comparison.}\label{table:density variance}
\subfloat[Density ${f}(\mathsf{x})$]{\resizebox{0.55\columnwidth}{!}{
\begin{tabular}{lrrrr}
\hline\hline
\multicolumn{1}{l}{}&\multicolumn{1}{c}{$p=1$}&\multicolumn{1}{c}{$p=2$}&\multicolumn{1}{c}{$p=3$}&\multicolumn{1}{c}{$p=4$}\tabularnewline
\hline
{\bfseries Kernel Function}&&&&\tabularnewline
~~Uniform&$0.600$&$0.600$&$1.250$&$1.250$\tabularnewline
~~Triangular&$0.743$&$0.743$&$1.452$&$1.452$\tabularnewline
~~Epanechnikov&$0.714$&$0.714$&$1.407$&$1.407$\tabularnewline
\hline
{\bfseries MD Variance Bound}&$0.500$&$0.500$&$1.125$&$1.125$\tabularnewline
\hline
\end{tabular}
}}
\vskip1em
\subfloat[Density Derivative ${f}^{(1)}(\mathsf{x})$]{\resizebox{0.55\columnwidth}{!}{
\begin{tabular}{lrrrr}
\hline\hline
\multicolumn{1}{l}{}&\multicolumn{1}{c}{$p=2$}&\multicolumn{1}{c}{$p=3$}&\multicolumn{1}{c}{$p=4$}&\multicolumn{1}{c}{$p=5$}\tabularnewline
\hline
{\bfseries Kernel Function}&&&&\tabularnewline
~~Uniform&$2.143$&$2.143$&$11.932$&$11.932$\tabularnewline
~~Triangular&$3.498$&$3.498$&$17.353$&$17.353$\tabularnewline
~~Epanechnikov&$3.182$&$3.182$&$15.970$&$15.970$\tabularnewline
\hline
{\bfseries MD Variance Bound}&$1.500$&$1.500$&$ 9.375$&$ 9.375$\tabularnewline
\hline
\end{tabular}
}}
\begin{flushleft}
\footnotesize
\textit{Notes}: Panel (a) compares asymptotic variance of the local polynomial density estimator of \cite*{Cattaneo-Jansson-Ma_2020_JASA} for different polynomial orders ($p=1$, $2$, $3$, and $4$) and different kernel functions (uniform, triangular and Epanechnikov). Also shown are the variance bound of the minimum distance estimator (MD Variance Bound), calculated according to Theorem \ref{thm:variance bound}. Panel(b) provides the same information for the estimated density derivative. All comparisons assume an interior evaluation point $\mathsf{x}$.
\end{flushleft}
\end{table}
\begin{example}[Local $p=5$ Minimum Distance Estimation]\label{example:p=5 v=1 interior}
We consider the same setting in Example \ref{example:p=1 v=1 interior} and \ref{example:p=3 v=1 interior}, but with $p=5$: $P(u) = (u, \frac{1}{2}u^2,\cdots, \frac{1}{5!}u^5)^\prime$.
The variance bound in Theorem \ref{thm:variance bound} is
\begin{align*}
f(\mathsf{x})\left(\int_{-1}^1 \dot{P}(u)\dot{P}(u)^\prime \mathrm{d} u\right)^{-1} = f(\mathsf{x})
\begin{bmatrix}
\frac{225}{128} & 0 & -\frac{525}{32} & 0 & \frac{2835}{16} \\
0 & \frac{75}{8} & 0 & -\frac{315}{4} & 0 \\
-\frac{525}{32} & 0 & \frac{2205}{8} & 0 & -\frac{14175}{4} \\
0 & -\frac{315}{4} & 0 & \frac{1575}{2} & 0 \\
\frac{2835}{16} & 0 & -\frac{14175}{4} & 0 & \frac{99225}{2} \\
\end{bmatrix}.
\end{align*}
Again, we include two higher order polynomials: $Q(u)=(u^{2j}, u^{2j+1})^\prime$. The asymptotic variance of our minimum distance estimator is
\begin{align*}
\mathrm{Asy}\mathbb{V}\begin{bmatrix}
\hat{f}_{\mathtt{MD}}(\mathsf{x})\\
\hat{f}_{\mathtt{MD}}^{(1)}(\mathsf{x})\\
\hat{f}_{\mathtt{MD}}^{(2)}(\mathsf{x})\\
\hat{f}_{\mathtt{MD}}^{(3)}(\mathsf{x})\\
\hat{f}_{\mathtt{MD}}^{(4)}(\mathsf{x})
\end{bmatrix} &= f(\mathsf{x})\begin{bmatrix}
\frac{225 (4 j+19)}{256 (2 j+9)} & 0 & -\frac{525 (4 j+21)}{64 (2 j+9)} & 0 & \frac{2835 (4 j+23)}{32 (2 j+9)} \\
0 & \frac{75 (4 j+17)}{16 (2 j+7)} & 0 & -\frac{315 (4 j+19)}{8 (2 j+7)} & 0 \\
-\frac{525 (4 j+21)}{64 (2 j+9)} & 0 & \frac{2205 (4 j+23)}{16 (2 j+9)} & 0 & -\frac{14175 (4 j+25)}{8 (2 j+9)} \\
0 & -\frac{315 (4 j+19)}{8 (2 j+7)} & 0 & \frac{1575 (4 j+21)}{8 j+28} & 0 \\
\frac{2835 (4 j+23)}{32 (2 j+9)} & 0 & -\frac{14175 (4 j+25)}{8 (2 j+9)} & 0 & \frac{99225 (4 j+27)}{8 j+36} \\
\end{bmatrix},
\end{align*}
which converges to the variance bound as $j\to \infty$. See also Table \ref{table:density variance} for the efficiency gain of employing the minimum distance technique.
\qed
\end{example}
Before closing this section, we make several remarks on the variance bound derived in Theorem \ref{thm:variance bound}, as well as to what extent it is achievable.
\begin{remark}[Achievability of the Variance Bound]
The previous two examples suggest that the variance bound derived in Theorem \ref{thm:variance bound} can be achieved by employing a minimum distance estimator with two additional regressors, one higher-order even polynomial and one higher-order odd polynomial. With analytic calculation, we verify that this is indeed the case for $p\leq 10$ when a uniform kernel function is used.
\qed
\end{remark}
\begin{remark}[Optimality of the Variance Bound]
\cite{Granovsky-Muller_1991_ISR} discuss the problem of finding the optimal kernel function for kernel-type estimators. To be precise, consider the following
\begin{align*}
\frac{1}{nh^{\ell+1}}\sum_{i=1}^n\phi_{\ell,k}\left( \frac{x_i-\mathsf{x}}{h} \right),
\end{align*}
where $\phi_{\ell,k}(u)$ is a function satisfying
\begin{align*}
\int_{-1}^1 u^{j}\phi_{\ell,k}(u)\mathrm{d} u &= \begin{cases}
0 & 0\leq j< k,\ j\neq \ell\\
\ell! & j=\ell
\end{cases},\qquad \int_{-1}^1 u^{k}\phi_{\ell,k}(u)\mathrm{d} u\neq 0.
\end{align*}
Then it is easy to see that, with a Taylor expansion argument,
\begin{align*}
\mathbb{E}\left[\frac{1}{nh^{\ell+1}}\sum_{i=1}^n\phi_{\ell,k}\left( \frac{x_i-\mathsf{x}}{h} \right)\right] &= \frac{1}{h^{\ell+1}}\int_{-1}^1 \phi_{\ell,k}\left( \frac{u-\mathsf{x}}{h} \right)f(u)\mathrm{d} u\\
&= \frac{1}{h^{\ell}}\int_{-1}^1 \phi_{\ell,k}\left( u \right)f(\mathsf{x}+hu)\mathrm{d} u\\
&= \frac{1}{h^{\ell}}\int_{-1}^1 \phi_{\ell,k}\left( u \right)\left[\sum_{j=0}^{k-1} \frac{(hu)^j}{j!}f^{(j)}(\mathsf{x}) + u^kO(h^k) \right]\mathrm{d} u\\
&= f^{(\ell)}(\mathsf{x}) + O(h^{k-\ell}).
\end{align*}
That is, the kernel $\phi_{\ell,k}(u)$ facilitates estimating the $\ell$-th derivative of the density function with a leading bias of order $h^{k-\ell}$. Asymptotic variance of this kernel-type estimator is easily found to be
\begin{align*}
\mathrm{Asy}\mathbb{V}\left[\frac{1}{nh^{\ell+1}}\sum_{i=1}^n\phi_{\ell,k}\left( \frac{x_i-\mathsf{x}}{h} \right)\right] &= f(\mathsf{x})\int_{-1}^1 \phi_{\ell,k}(u)^2\mathrm{d} u.
\end{align*}
\cite{Granovsky-Muller_1991_ISR} provide the exact form of the kernel function $\phi_{\ell,k}(u)$ that minimizes the asymptotic variance subject to the order of the leading bias.
Take $\ell=0$ and $k=2$, $\phi_{\ell,k}(u)$ takes the following form:
\begin{align*}
\phi_{\ell,k}(u) &= \frac{1}{2}\mathds{1}(|u|\leq 1),
\end{align*}
which is the uniform kernel and minimizes variance among all second order kernels for density estimation. As illustrated in Example \ref{example:p=1 v=1 interior}, our variance bound matches $f(\mathsf{x})\int_{-1}^1 \phi_{\ell,k}(u)^2\mathrm{d} u$.
Now take $\ell=1$ and $k=3$. This will give an estimator for the density derivative $f^{(1)}(\mathsf{x})$ with a leading bias of order $O(h^2)$. The optimal choice of $\phi_{\ell,k}(u)$ is
\begin{align*}
\phi_{\ell,k}(u) &= \frac{3}{2}u\mathds{1}(|u|\leq 1).
\end{align*}
to match the order of bias, we consider the minimum distance estimator with $p=3$. Again, the variance bound in Theorem \ref{thm:variance bound} matches $f(\mathsf{x})\int_{-1}^1 \phi_{\ell,k}(u)^2\mathrm{d} u$.
As a final illustration, take $\ell=1$ and $k=5$, which gives an estimator for the density derivative $f^{(1)}(\mathsf{x})$ with a leading bias of order $O(h^4)$. The optimal choice of $\phi_{\ell,k}(u)$ is
\begin{align*}
\phi_{\ell,k}(u) &= \left(\frac{75}{8}u - \frac{105}{8}u^3\right)\mathds{1}(|u|\leq 1).
\end{align*}
It is easy to see that $f(\mathsf{x})\int_{-1}^1 \phi_{\ell,k}(u)^2\mathrm{d} u = 75f(\mathsf{x})/8$. To match the bias order, we take $p=5$ for our minimum distance estimator. The variance bound is $75f(\mathsf{x})/8$, which is the same as $f(\mathsf{x})\int_{-1}^1 \phi_{\ell,k}(u)^2\mathrm{d} u$.
With analytic calculations, we verify that the variance bound stated in Theorem \ref{thm:variance bound} is the same as the minimum variance found in \cite{Granovsky-Muller_1991_ISR}. Together with the previous remark, we reach a much stronger conclusion: including two higher-order polynomials in our minimum distance estimator can help achieve the variance bound in Theorem \ref{thm:variance bound}, which, in turn, is the smallest variance any kernel-type estimator can achieve (given a specific leading bias order).
\qed
\end{remark}
\begin{remark}[Another Density Estimator Which Achieves the Variance Bound]
The following estimator achieves the bound of Theorem \ref{thm:variance bound}, although it does not belong to the class of estimators we consider in this paper.
\begin{align*}
\hat{\theta}_{\mathtt{ND}} &= \left( \int_{\mathcal{X}}\dot{P}(u-\mathsf{x})\dot{P}(u-\mathsf{x})^\prime \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} u \right)^{-1}\left(\frac{1}{n}\sum_{i=1}^n \dot{P}(x_i-\mathsf{x})\frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right)\right),
\end{align*}
where $\dot{P}(u)=(1,u,u^2/2,\cdots,u^{p-1}/(p-1)!)^\prime$ is the $(p-1)$-th order polynomial basis. The subscript represents ``numerical derivative,'' because the above estimator can be understood as
\begin{align*}
\hat{\theta}_{\mathtt{ND}} &= \left( \int_{\mathcal{X}}\dot{P}(u-\mathsf{x})\dot{P}(u-\mathsf{x})^\prime \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} u \right)^{-1}\left(\int_{\mathcal{X}} \dot{P}(u-\mathsf{x})\frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right)\frac{\mathrm{d} \hat{F}(u)}{\mathrm{d} u}\mathrm{d} u\right)\\
&= \operatorname*{argmin}_{\theta}\int_{\mathcal{X}} \left( \frac{\mathrm{d} \hat{F}(u)}{\mathrm{d} u} - \dot{P}(u-\mathsf{x})^\prime \theta \right)^2 \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right)\mathrm{d} u,
\end{align*}
where the derivative $\mathrm{d} \hat{F}(u)/\mathrm{d} u$ is interpreted in the sense of generalized functions. From the above, it is clear that this estimator requires the knowledge of the boundary position (that is, the knowledge of $\mathcal{X}$).
With straightforward calculations, this estimator has a leading bias
\begin{align*}
\mathbb{E}[\hat{\theta}_{\mathtt{ND}}] &= \left( \int_{\mathcal{X}}\dot{P}(u-\mathsf{x})\dot{P}(u-\mathsf{x})^\prime \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} u \right)^{-1}\mathbb{E}\left[ \dot{P}(x_i-\mathsf{x})\frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right) \right]\\
&= \theta + h^{p}\Upsilon_{h}f^{(p)}(\mathsf{x})\left( \int_{\frac{\mathcal{X}-\mathsf{x}}{h}}\dot{P}\left(u\right)\dot{P}\left(u\right)^\prime K\left(u\right) \mathrm{d} u \right)^{-1} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}}\dot{P}\left(u\right)u^pK\left(u\right) \mathrm{d} u + o(h^{p}\Upsilon_{h}),
\end{align*}
where $\Upsilon_{h}$ is a diagonal matrix containing $1$, $h^{-1}$, $\cdots$, $h^{-(p-1)}$. Its leading variance is also easy to establish:
\begin{align*}
\mathbb{V}[\hat{\theta}_{\mathtt{ND}}] &= \frac{1}{nh}\Upsilon_{h}f(\mathsf{x})\left( \int_{\frac{\mathcal{X}-\mathsf{x}}{h}}\dot{P}\left(u\right)\dot{P}\left(u\right)^\prime K\left(u\right) \mathrm{d} u \right)^{-1}
\left( \int_{\frac{\mathcal{X}-\mathsf{x}}{h}}\dot{P}\left(u\right)\dot{P}\left(u\right)^\prime K\left(u\right)^2 \mathrm{d} u \right)\\
&\qquad \qquad \qquad\qquad \cdot\left( \int_{\frac{\mathcal{X}-\mathsf{x}}{h}}\dot{P}\left(u\right)\dot{P}\left(u\right)^\prime K\left(u\right) \mathrm{d} u \right)^{-1}\Upsilon_{h}\\
&+ o\left(\frac{1}{nh}\Upsilon_{h}^2\right).
\end{align*}
To reach the efficiency bound in Theorem \ref{thm:variance bound}, it suffices to set $K(\cdot)$ to be the uniform kernel. Section 5.1.1 in \cite{Loader_2006_book_Book} also discussed this estimator, although it seems its efficiency property has not been realized in the literature.
\qed
\end{remark}
\section{Uniform Distribution Theory}\label{section:Uniform Distribution Theory}
We establish distribution approximation for $\{\hat\theta_G(\mathsf{x}),\mathsf{x}\in\mathcal{I}\}$ and $\{\hat\theta(\mathsf{x}),\mathsf{x}\in\mathcal{I}\}$, which can be viewed as processes indexed by the evaluation point $\mathsf{x}$ in some set $\mathcal{I}\subseteq \mathcal{X}$. Recall the definition of $\Gamma_{h,\mathsf{x}}$ and $\Sigma_{h,\mathsf{x}}$ from Section \ref{section:setup}, and we define $\Omega_{h,\mathsf{x}} = \Gamma_{h,\mathsf{x}}^{-1}\Sigma_{h,\mathsf{x}}\Gamma_{h,\mathsf{x}}^{-1}$.
We first study the following (infeasible) centered and Studentized process:
\begin{align}\label{eq:centered and studentized process}
\mathfrak{T}_G(\mathsf{x}) &= \frac{1}{\sqrt{n}}\sum_{i=1}^n\frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{x}}^{-1} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}},\ \mathsf{x} \in \mathcal{I},
\end{align}
where we consider linear combinations through a (known) vector $c_{h,\mathsf{x}}$, which can depend on the sample size through the bandwidth $h$, and can depend on the evaluation point. Again, we use the subscript $G$ to denote the local $L^2$ approach with $G$ being the design distribution. To economize notation, let
\begin{align*}
\mathcal{K}_{h,\mathsf{x}}(x) &= \frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{x}}^{-1} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}},
\end{align*}
then we can conveniently rewrite \eqref{eq:centered and studentized process} as
\begin{align*}
\mathfrak{T}_G(\mathsf{x}) &= \frac{1}{\sqrt{n}}\sum_{i=1}^n \mathcal{K}_{h,\mathsf{x}}(x_i),
\end{align*}
and hence the centered and Studentized process $\mathfrak{T}_G(\cdot)$ takes a kernel form. The difference compared to standard kernel density estimators, however, is that the (equivalent) kernel in our case changes with the evaluation point, which is why our estimator is able to adapt to boundary points automatically. From the pointwise distribution theory developed in Section \ref{section:pointwise distribution theory}, the process $\mathfrak{T}_G(\mathsf{x})$ has variance
\begin{align*}
\mathbb{V}\left[ \mathfrak{T}_G(\mathsf{x}) \right] = \mathbb{E}\left[ \mathcal{K}_{h,\mathsf{x}}(x_i)^2 \right] = 1.
\end{align*}
We can also compute the covariance as
\begin{align*}
\mathbb{C}\mathrm{ov}\left[ \mathfrak{T}_G(\mathsf{x}),\mathfrak{T}_G(\mathsf{y}) \right] = \mathbb{E}\left[\mathcal{K}_{h,\mathsf{x}}(x_i)\mathcal{K}_{h,\mathsf{y}}(x_i)\right] = \frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x},\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}} + O(h),
\end{align*}
where $\Omega_{h,\mathsf{x},\mathsf{y}} = \Gamma_{h,\mathsf{x}}^{-1}\Sigma_{h,\mathsf{x},\mathsf{y}}\Gamma_{h,\mathsf{y}}^{-1}$, and
\begin{align*}
\Sigma_{h,\mathsf{x},\mathsf{y}} &= \int_{\frac{\mathcal{X}-\mathsf{y}}{h}}\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(v)^\prime \Big[ F((\mathsf{x}+hu)\wedge (\mathsf{y}+hv)) - F(\mathsf{x}+hu)F(\mathsf{y}+hv) \Big] \\
&\qquad \qquad \qquad \qquad K(u)K(v)g(\mathsf{x}+hu)g(\mathsf{y}+hv)\mathrm{d} u\mathrm{d} v.
\end{align*}
Of course one can further expand the above, but this is unnecessary for our purpose.
For future reference, let
\begin{align*}
r_1(\varepsilon,h)&=\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I},|\mathsf{x}-\mathsf{y}|\leq \varepsilon}\left|c^\prime_{h,\mathsf{x}}\Upsilon_{h} - c^\prime_{h,\mathsf{y}}\Upsilon_{h}\right|,\qquad r_2(h) = \sup_{\mathsf{x}\in\mathcal{I}} \frac{1}{|c_{h,\mathsf{x}}^\prime\Upsilon_{h}|}.
\end{align*}
\begin{remark}[On the Order of $r_1(\varepsilon,h)$, $r_2(h)$ and $\sup_{\mathsf{x}\in\mathcal{I}}\varrho(h,\mathsf{x})$]\label{remark: order of some errors}
In general, it is not possible to give precise orders of the quantities introduced above. In this remark, we consider the local polynomial estimator of \cite*{Cattaneo-Jansson-Ma_2020_JASA} (see Section \ref{section:setup} for an introduction). The local polynomial estimator employs a polynomial basis, and hence estimates the density function and higher-order derivatives by (it also estimates the distribution function)
\begin{align*}
\hat{F}^{(\ell)}(\mathsf{x}) &= e_\ell^\prime \hat{\theta}(\mathsf{x}),
\end{align*}
where $e_\ell$ is the $(\ell+1)$-th standard basis vector. As a result, $c_{h,\mathsf{x}} = e_\ell$, which does not depend on the evaluation point. For the scaling matrix $\Upsilon_{h}$, we note that it is diagonal with elements $1, h^{-1},\cdots, h^{-p}$, and hence it does not depend on the evaluation point either. Therefore, we conclude that, for density (and higher-order) derivative estimation using the local polynomial estimator, $r_1(\varepsilon,h)$ is identically zero. Similarly, we have that $r_2(h) = h^{\ell}$. Finally, given the discussion in Section \ref{section:setup}, the bias term generally has order $\sup_{\mathsf{x}\in\mathcal{I}}\varrho(h,\mathsf{x}) = h^{p+1}$ for the local polynomial density estimator.
The above discussion restricts to the local polynomial density estimator, but more can be said about $r_2(h)$. We will argue that, in general, one should expect $r_2(h) = O(1)$. Recall that the leading variance of $c_{h,\mathsf{x}}^\prime \hat{\theta}(\mathsf{x})$ and $c_{h,\mathsf{x}}^\prime \hat{\theta}_G(\mathsf{x})$ is $\frac{1}{n}c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}$, and that the maximum eigenvalue of $\Omega_{h,\mathsf{x}}$ is bounded. Therefore, the variance has order $O(1/(nr_2(h)^2))$. In general, we do not expect the variance to shrink faster than $1/n$, which is why $r_2(h)$ is usually bounded. In fact, for most interesting cases, $c_{h,\mathsf{x}}^\prime \hat{\theta}(\mathsf{x})$ and $c_{h,\mathsf{x}}^\prime \hat{\theta}_G(\mathsf{x})$ will be ``nonparametric'' estimators in the sense that they estimate local features of the distribution function. If this is the case, we may even argue that $r_2(h)$ will be vanishing as the bandwidth shrinks.
\qed
\end{remark}
We also make some additional assumptions.
\begin{assumption}\label{assumption:uniform}
Let $\mathcal{I}$ be a compact interval.
(i) The density function is twice continuously differentiable and bounded away from zero in $\mathcal{I}$.
(ii) There exists some $\delta>0$ and compactly supported kernel functions $K^{\dag}(\cdot)$ and $\{K^{\ddag,d}(\cdot) \}_{d\leq \delta}$, such that (ii.1) $\sup_{u\in\mathbb{R}}| K^{\dag}(u) |, \sup_{d\leq \delta,u\in\mathbb{R}} | K^{\ddag,d}(u) |<\infty$; (ii.2) the support of $K^{\ddag,d}(\cdot)$ has Lebesgue measure bounded by $Cd$, where $C$ is independent of $d$; and (ii.3)
for all $u$ and $v$ such that $|u-v|\leq \delta$,
\begin{align*}
|K(u)-K(v)| \leq |u-v|\cdot K^{\dag}(u) + K^{\ddag,|u-v|}(u).
\end{align*}
(iii) The basis function $R(\cdot)$ is Lipschitz continuous in $[-1,1]$.
(iv) For all $h$ sufficiently small, the minimum eigenvalues of $\Gamma_{h,\mathsf{x}}$ and $h^{-1}\Sigma_{h,\mathsf{x}}$ are bounded away from zero uniformly for $\mathsf{x}\in\mathcal{I}$.
(v) $h\to 0$ and $nh/\log n\to\infty$ as $n\to \infty$.
(vi) For some $C_1>0$ and $C_2,\ C_3\geq0$,
\begin{align*}
r_1(\varepsilon,h) = O\left(\varepsilon^{C_1}h^{-C_2}\right),\qquad
r_2(h)= O \left(h^{C_3}\right).
\end{align*}
In addition,
\begin{align*}
\frac{\sup_{\mathsf{x}\in\mathcal{I}}|c_{h,\mathsf{x}}^\prime\Upsilon_{h}|}{\inf_{\mathsf{x}\in\mathcal{I}}|c_{h,\mathsf{x}}^\prime\Upsilon_{h}|} = O(1).
\end{align*}
\vskip-2em\qed
\end{assumption}
\begin{assumption}\label{assumption:design uniform}
The design density function $g(\cdot)$ is twice continuously differentiable and is bounded away from zero in $\mathcal{I}$.
\qed
\end{assumption}
For any $h>0$ (and fixed $n$), we can define a centered Gaussian process, $\{\mathfrak{B}_G(\mathsf{x}):\mathsf{x}\in\mathcal{I}\}$, which has the same variance-covariance structure as the process $\mathfrak{T}_G(\cdot)$. The following lemma shows that it is possible to construct such a process, and that $\mathfrak{T}_G(\cdot)$ and $\mathfrak{B}_G(\cdot)$ are ``close in distribution.''
\begin{thm}[Strong Approximation]\label{thm:strong approximation}
Assume Assumptions \ref{assumption:dgp pointwise}, \ref{assumption:design pointwise}, \ref{assumption:uniform} and \ref{assumption:design uniform} hold.
Then on a possibly enlarged probability space there exist two processes, $\{\tilde{\mathfrak{T}}_G(\mathsf{x}):\mathsf{x}\in\mathcal{I}\}$ and $\{\mathfrak{B}_G(\mathsf{x}):\mathsf{x}\in\mathcal{I}\}$, such that (i) $\tilde{\mathfrak{T}}_G(\cdot)$ has the same distribution as $\mathfrak{T}_G(\cdot)$; (ii) $\mathfrak{B}_G(\cdot)$ is a Gaussian process with the same covariance structure as $\mathfrak{T}_G(\cdot)$; and (iii)
\begin{align*}
\mathbb{P}\left[ \sup_{\mathsf{x}\in\mathcal{I}}\left| \tilde{\mathfrak{T}}_G(\mathsf{x}) - \mathfrak{B}_G(\mathsf{x}) \right| > \frac{C_4(u+C_5\log n)}{\sqrt{nh}} \right] \leq C_5e^{-C_5u},
\end{align*}
where $C_5$ is some constant that does not depend on $h$ or $n$.
\qed
\end{thm}
Next we consider the continuity property of the implied (equivalent) kernel of the process $\mathfrak{T}_G(\cdot)$, which will help control the complexity of the Gaussian process $\mathfrak{B}_G(\cdot)$. To be precise, define the pseudo-metric $\sigma_G(\mathsf{x},\mathsf{y})$ as
\begin{align*}
\sigma_G(\mathsf{x},\mathsf{y}) &= \sqrt{\mathbb{V}\left[ \mathfrak{T}_G(\mathsf{x}) - \mathfrak{T}_G(\mathsf{y}) \right]} = \sqrt{\mathbb{E}\left[ (\mathcal{K}_{h,\mathsf{x}}(x_i) - \mathcal{K}_{h,\mathsf{y}}(x_i))^2 \right]},
\end{align*}
we would like to provide an upper bound of $\sigma_G(\mathsf{x},\mathsf{y})$ in terms of $|\mathsf{x}-\mathsf{y}|$ (at least for all $\mathsf{x}$ and $\mathsf{y}$ such that $|\mathsf{x}-\mathsf{y}|$ is small enough).
\begin{lem}[VC-type Property]\label{lem:VC-type}
Assume Assumptions \ref{assumption:dgp pointwise}, \ref{assumption:design pointwise}, \ref{assumption:uniform} and \ref{assumption:design uniform} hold. Then for all $\mathsf{x},\mathsf{y}\in\mathcal{I}$ such that $|\mathsf{x}-\mathsf{y}|=\varepsilon\leq h$,
\begin{align*}
\sigma_G(\mathsf{x},\mathsf{y}) = O\left( \frac{1}{\sqrt{h}}\frac{\varepsilon}{h} + \frac{1}{\sqrt{h}}r_1(\varepsilon,h)r_2(h) + \frac{1}{h}r_1(\varepsilon,h)^2r_2(h)^2 \right).
\end{align*}
Therefore,
\begin{align*}
\mathbb{E}\left[ \sup_{\mathsf{x}\in\mathcal{I}}|\mathfrak{B}_G(\mathsf{x})| \right] = O\left(\sqrt{\log n}\right),\qquad \text{and}\qquad \mathbb{E}\left[ \sup_{\mathsf{x}\in\mathcal{I}}|\mathfrak{T}_G(\mathsf{x})| \right] = O\left(\sqrt{\log n} \right).
\end{align*}
\vskip-2em\qed
\end{lem}
\subsection{Local $L^2$ Distribution Estimation}
We first discuss the covariance estimator. For the local $L^2$ distribution estimator, let $\hat{\Omega}_{h,\mathsf{x},\mathsf{y}} = \Gamma_{h,\mathsf{x}}^{-1}\hat\Sigma_{h,\mathsf{x},\mathsf{y}}\Gamma_{h,\mathsf{y}}^{-1}$ with $\hat\Sigma_{h,\mathsf{x},\mathsf{y}}$ given by
\begin{align*} \hat\Sigma_{h,\mathsf{x},\mathsf{y}} &= \frac{1}{n}\sum_{i=1}^n\int_{\frac{\mathcal{X}-\mathsf{y}}{h}}\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(v)^\prime\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - \hat{F}(\mathsf{x} + hu)\Big]\Big[\mathds{1}(x_i\leq \mathsf{y} + hv) - \hat{F}(\mathsf{y} + hv)\Big] \\
&\ \qquad \qquad \qquad K(u)K(v) g(\mathsf{x}+hu)g(\mathsf{y}+hv)\mathrm{d} u\mathrm{d} v.
\end{align*}
The next lemma characterizes the convergence rate of $\hat{\Omega}_{h,\mathsf{x},\mathsf{y}}$.
\begin{lem}[Local $L^2$ Distribution Estimation: Covariance Estimation] \label{lemma:standard error, local projection uniform}
Assume Assumptions \ref{assumption:dgp pointwise}, \ref{assumption:design pointwise}, \ref{assumption:uniform} and \ref{assumption:design uniform} hold, and that $nh^2/\log n\to \infty$. Then
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}}\left|\frac{c^\prime_{h,\mathsf{x}}\Upsilon_h(\hat{\Omega}_{h,\mathsf{x},\mathsf{y}}-{\Omega}_{h,\mathsf{x},\mathsf{y}})\Upsilon_h c_{h,\mathsf{y}} }{\sqrt{c^\prime_{h,\mathsf{x}}\Upsilon_h{\Omega}_{h,\mathsf{x}}\Upsilon_h c_{h,\mathsf{x}}}\sqrt{c^\prime_{h,\mathsf{y}}\Upsilon_h{\Omega}_{h,\mathsf{y}}\Upsilon_h c_{h,\mathsf{y}}}}\right| = O_{\mathbb{P}}\left( \sqrt{\frac{\log n}{nh^2}} \right)
\end{align*}
\vskip-2em\qed
\end{lem}
We now consider the estimator $c^\prime_{h,\mathsf{x}}\hat\theta_G(\mathsf{x})$. From \eqref{eq:approximation bias, projection} and \eqref{eq:linear variance, projection}, one has
\begin{align}
\nonumber T_G(\mathsf{x}) &= \frac{\sqrt{n}c_{h,\mathsf{x}}^\prime\Big(\hat\theta_{G}(\mathsf{x}) - \theta(\mathsf{x})\Big)}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\\
\label{eq:local projection uniform term 1}&= \sqrt{n}\frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Gamma_h^{-1}\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[F(\mathsf{x} + hu) - \theta^\prime R(u)\Upsilon_h^{-1}\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\\
\label{eq:local projection uniform term 2}&\quad+ \frac{1}{\sqrt{n}}\sum_{i=1}^n\frac{c^\prime_{h,\mathsf{x}}\Upsilon_{h}\Gamma_{h,\mathsf{x}}^{-1} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}.
\end{align}
In the following lemma, we analyze the two terms in the above decomposition.
\begin{lem}\label{lem:local projection uniform term 1}
Assume Assumptions \ref{assumption:dgp pointwise}, \ref{assumption:design pointwise}, \ref{assumption:uniform} and \ref{assumption:design uniform} hold, and that $nh^2/\log n\to \infty$. Then
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}\Big|\text{\eqref{eq:local projection uniform term 1}}\Big| &= O_{\mathbb{P}}\left( \sqrt{\frac{n}{h}}\sup_{\mathsf{x}\in\mathcal{I}}\varrho(h,\mathsf{x}) \right),\qquad \sup_{\mathsf{x}\in\mathcal{I}}\Big|\text{\eqref{eq:local projection uniform term 2}} - \mathfrak{T}_G(\mathsf{x})\Big| = O_{\mathbb{P}}\left( \frac{\log n}{\sqrt{nh^2}} \right).
\end{align*}
\vskip-2em\qed
\end{lem}
Now we state the main result on uniform distributional approximation.
\begin{thm}[Local $L^2$ Distribution Estimation: Uniform Distributional Approximation]\label{thm:strong approximation, local projection}
Assume Assumptions \ref{assumption:dgp pointwise}, \ref{assumption:design pointwise}, \ref{assumption:uniform} and \ref{assumption:design uniform} hold, and that $nh^2/\log n\to \infty$. Then on a possibly enlarged probability space there exist two processes, $\{\tilde{\mathfrak{T}}_G(\mathsf{x}):\mathsf{x}\in\mathcal{I}\}$ and $\{\mathfrak{B}_G(\mathsf{x}):\mathsf{x}\in\mathcal{I}\}$, such that (i) $\tilde{\mathfrak{T}}_G(\cdot)$ has the same distribution as $\mathfrak{T}_G(\cdot)$; (ii) $\mathfrak{B}_G(\cdot)$ is a Gaussian process with the same covariance structure as $\mathfrak{T}_G(\cdot)$; and (iii)
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}\Big| T_G(\mathsf{x}) - \mathfrak{T}_G(\mathsf{x}) \Big| + \sup_{\mathsf{x}\in\mathcal{I}}\Big| \tilde{\mathfrak{T}}_G(\mathsf{x}) - \mathfrak{B}_G(\mathsf{x}) \Big| = O_{\mathbb{P}}\left(\frac{\log n}{\sqrt{nh^2}}+\sqrt{\frac{n}{h}}\sup_{\mathsf{x}\in\mathcal{I}}\varrho(h,\mathsf{x}) \right).
\end{align*}
\vskip-2em\qed
\end{thm}
The following theorem shows that a feasible approximation to the process $\mathfrak{B}_G(\cdot)$ can be achieved by simulating a Gaussian process with covariance estimated from the data. In the following, we use $\mathbb{P}^\star$, $\mathbb{E}^\star$ and $\mathbb{C}\mathrm{ov}^\star$ to denote the probability, expectation and covariance operator conditioning on the data $X_n=(x_1,x_2,\dots,x_n)'$.
\begin{thm}[Local $L^2$ Distribution Estimation: Feasible Distributional Approximation]\label{thm:feasible uniform approximation local projection}
Assume Assumptions \ref{assumption:dgp pointwise}, \ref{assumption:design pointwise}, \ref{assumption:uniform} and \ref{assumption:design uniform} hold, and that $nh^2/\log n\to \infty$. Then conditional on the data there exists a centered Gaussian process $\hat{\mathfrak{B}}_G(\cdot)$ with covariance
\begin{align*}
\mathbb{C}\mathrm{ov}^\star\left[ \hat{\mathfrak{B}}_G(\mathsf{x}),\hat{\mathfrak{B}}_G(\mathsf{y}) \right] &= \frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x},\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}},
\end{align*}
such that
\begin{align*}
\sup_{u\in \mathbb{R}}\left| \mathbb{P}\Big[ \sup_{\mathsf{x}\in\mathcal{I}}|\mathfrak{B}_G(\mathsf{x})| \leq u \Big]
-
\mathbb{P}^\star\Big[\sup_{\mathsf{x}\in\mathcal{I}}|\hat{\mathfrak{B}}_G(\mathsf{x})| \leq u \Big] \right| = O_{\mathbb{P}}\left(\left(\frac{\log^{5} n}{nh^2}\right)^\frac{1}{4} \right).
\end{align*}
\vskip-2em\qed
\end{thm}
\begin{remark}[On the Remainders in Theorems \ref{thm:strong approximation, local projection} and \ref{thm:feasible uniform approximation local projection}]\label{remark:bounding errors of strong approximation}
Recall that the local polynomial density estimator employs a polynomial basis, which implies that $\sup_{\mathsf{x}\in\mathcal{I}}\varrho(h,\mathsf{x}) = h^{p+1}$, where $p$ is the highest polynomial order. Then the error in Theorem \ref{thm:strong approximation, local projection} reduces to
\begin{align*}
\sqrt{nh^{2p+1}} + \frac{\log n}{\sqrt{nh^2}}.
\end{align*}
Therefore, a sufficient set of conditions for both errors to be negligible is $nh^{2p+1}\to 0$ and $nh^2/\log^5 n\to \infty$.
\qed
\end{remark}
\subsection{Local Regression Distribution Estimation}
Now we consider the local regression estimator $\{\hat\theta(\mathsf{x}),\mathsf{x}\in\mathcal{I}\}$. As before, we first discuss the construction of the covariance $\Omega_{h,\mathsf{x},\mathsf{y}}$. Let $\hat{\Omega}_{h,\mathsf{x},\mathsf{y}} = \hat{\Gamma}_{h,\mathsf{x}}^{-1}\hat{\Sigma}_{h,\mathsf{x},\mathsf{y}}\hat{\Gamma}_{h,\mathsf{y}}^{-1}$. Construction of $\hat{\Gamma}_{h,\mathsf{x}}$ is given in Section \ref{subsection:local regression: pointwise distribution theory}. The following lemma shows that $\hat{\Gamma}_{h,\mathsf{x}}$ is uniformly consistent.
\begin{lem}[Uniform Consistency of $\hat{\Gamma}_{h,\mathsf{x}}$]\label{lem:uniform consistency of Gamma}
Assume Assumptions \ref{assumption:dgp pointwise} and \ref{assumption:uniform} hold. Then
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}\left|\hat{\Gamma}_{h,\mathsf{x}}-\Gamma_{h,\mathsf{x}} \right| = O_{\mathbb{P}}\left(\sqrt{\frac{\log n}{nh}}\right).
\end{align*}
\vskip-2em\qed
\end{lem}
Construction of $\hat{\Sigma}_{h,\mathsf{x},\mathsf{y}}$ also mimics that in Section \ref{subsection:local regression: pointwise distribution theory}. To be precise, we let
\begin{align*}
\hat{\Sigma}_{h,\mathsf{x}, \mathsf{y}} &= \frac{1}{n}\sum_{i=1}^n \left[\frac{1}{n}\sum_{j=1}^n \Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -\hat{F}(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right)\right]\\
&\qquad \qquad \qquad \left[\frac{1}{n}\sum_{j=1}^n \Upsilon_hR(x_j-\mathsf{y})\Big[ \mathds{1}(x_i\leq x_j) -\hat{F}(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{y}}{h}\right)\right]^\prime.
\end{align*}
where $\hat{F}(\cdot)$ remains to be the empirical distribution function. The following result justifies consistency of $\hat{\Omega}_{h,\mathsf{x},\mathsf{y}}$.
\begin{lem}[Local Regression Distribution Estimation: Covariance Estimation] \label{lemma:standard error, local regression uniform}
Assume Assumptions \ref{assumption:dgp pointwise} and \ref{assumption:uniform} hold, and that $nh^2/\log n\to \infty$. Then
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}}\left|\frac{c^\prime_{h,\mathsf{x}}\Upsilon_h(\hat{\Omega}_{h,\mathsf{x},\mathsf{y}}-{\Omega}_{h,\mathsf{x},\mathsf{y}})\Upsilon_h c_{h,\mathsf{y}} }{\sqrt{c^\prime_{h,\mathsf{x}}\Upsilon_h{\Omega}_{h,\mathsf{x}}\Upsilon_h c_{h,\mathsf{x}}}\sqrt{c^\prime_{h,\mathsf{y}}\Upsilon_h{\Omega}_{h,\mathsf{y}}\Upsilon_h c_{h,\mathsf{y}}}}\right| = O_{\mathbb{P}}\left( \sqrt{\frac{\log n}{nh^2}} \right)
\end{align*}
\vskip-2em\qed
\end{lem}
The following is an expansion of $T(\cdot)$.
\begin{align}
\nonumber T(\mathsf{x}) &= \frac{\sqrt{n}c_{h,\mathsf{x}}^\prime\Big(\hat\theta(\mathsf{x}) - \theta(\mathsf{x})\Big)}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\\
\label{eq:local regression uniform term 1}&= \frac{1}{n\sqrt{n}}\sum_{i=1}^n \frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Gamma}_{h,\mathsf{x}}^{-1}\Upsilon_{h}R(x_i-\mathsf{x})[1-F(x_i)]\frac{1}{h}K(\frac{x_i-\mathsf{x}}{h})}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\\
\label{eq:local regression uniform term 2}&+ \frac{1}{\sqrt{n}}\sum_{i=1}^n \frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Gamma}_{h,\mathsf{x}}^{-1}\Upsilon_{h}R(x_i-\mathsf{x})[F(x_i)-\theta(\mathsf{x})^\prime R(x_i-\mathsf{x})]\frac{1}{h}K(\frac{x_i-\mathsf{x}}{h})}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\\
\nonumber &+ \frac{1}{n\sqrt{n}}\sum_{i,j=1,i\neq j}^n \frac{1}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}} \Bigg\{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Gamma}_{h,\mathsf{x}}^{-1}\Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -F(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right) \\
\label{eq:local regression uniform term 3}&\qquad \qquad \qquad\qquad- \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Gamma}_{h,\mathsf{x}}^{-1}R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) f(\mathsf{x} + hu)\mathrm{d} u\Bigg\}\\
\label{eq:local regression uniform term 4}&+ \frac{n-1}{n\sqrt{n}}\sum_{i=1}^n\frac{c^\prime_{h,\mathsf{x}}\Upsilon_{h}\hat{\Gamma}_{h,\mathsf{x}}^{-1} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) f(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}.
\end{align}
\begin{lem}\label{lem:local regression uniform term 1}
Assume Assumptions \ref{assumption:dgp pointwise} and \ref{assumption:uniform} hold, and that $nh^2/\log n\to \infty$. Then
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}\Big|\text{\eqref{eq:local regression uniform term 1}}\Big| &= O_{\mathbb{P}}\left( \frac{1}{\sqrt{nh}}\right),\quad
\sup_{\mathsf{x}\in\mathcal{I}}\Big|\text{\eqref{eq:local regression uniform term 2}}\Big| = O_{\mathbb{P}}\Big( \sqrt{\frac{n}{h}}\sup_{\mathsf{x}\in\mathcal{I}}\varrho(h,\mathsf{x})\Big),\quad
\sup_{\mathsf{x}\in\mathcal{I}}\Big|\text{\eqref{eq:local regression uniform term 3}}\Big| = O_{\mathbb{P}}\left( \frac{\log n}{\sqrt{nh^2}}\right).
\end{align*}
\vskip-2em\qed
\end{lem}
\begin{lem}\label{lem:local regression uniform term 2}
Assume Assumptions \ref{assumption:dgp pointwise} and \ref{assumption:uniform} hold, and that $nh^2/\log n\to \infty$. Then
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}\Big|\text{\eqref{eq:local regression uniform term 4}} - \mathfrak{T}_F(\mathsf{x})\Big| &= O_{\mathbb{P}}\left( \frac{\log n}{\sqrt{nh^2}}\right).
\end{align*}
\vskip-2em\qed
\end{lem}
Finally we have the following result on uniform distributional approximation for the local regression distribution estimator, as well as a feasible approximation by simulating from a Gaussian process with estimated covariance.
\begin{thm}[Local Regression Distribution Estimation: Uniform Distributional Approximation]\label{thm:strong approximation, local regression}
Assume Assumptions \ref{assumption:dgp pointwise} and \ref{assumption:uniform} hold, and that $nh^2/\log n\to \infty$. Then on a possibly enlarged probability space there exist two processes, $\{\tilde{\mathfrak{T}}_F(\mathsf{x}):\mathsf{x}\in\mathcal{I}\}$ and $\{\mathfrak{B}_F(\mathsf{x}):\mathsf{x}\in\mathcal{I}\}$, such that (i) $\tilde{\mathfrak{T}}_F(\cdot)$ has the same distribution as $\mathfrak{T}_F(\cdot)$; (ii) $\mathfrak{B}_F(\cdot)$ is a Gaussian process with the same covariance structure as $\mathfrak{T}_F(\cdot)$; and (iii)
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}\Big| T(\mathsf{x}) - \mathfrak{T}_F(\mathsf{x}) \Big| + \sup_{\mathsf{x}\in\mathcal{I}}\Big| \tilde{\mathfrak{T}}_F(\mathsf{x}) - \mathfrak{B}_F(\mathsf{x}) \Big| = O_{\mathbb{P}}\left( \frac{\log n}{\sqrt{nh^2}}+\sqrt{\frac{n}{h}}\sup_{\mathsf{x}\in\mathcal{I}}\varrho(h,\mathsf{x}) \right).
\end{align*}
\vskip-2em\qed
\end{thm}
\begin{thm}[Local Regression Distribution Estimation: Feasible Distributional Approximation]\label{thm:feasible uniform approximation local regression}
Assume Assumptions \ref{assumption:dgp pointwise} and \ref{assumption:uniform} hold, and that $nh^2/\log n\to \infty$. Then conditional on the data there exists a centered Gaussian process $\hat{\mathfrak{B}}_F(\cdot)$ with covariance
\begin{align*}
\mathbb{C}\mathrm{ov}^\star\left[ \hat{\mathfrak{B}}_F(\mathsf{x}),\hat{\mathfrak{B}}_F(\mathsf{y}) \right] &= \frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x},\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}},
\end{align*}
such that
\begin{align*}
\sup_{u\in \mathbb{R}}\left| \mathbb{P}\Big[ \sup_{\mathsf{x}\in\mathcal{I}}|\mathfrak{B}_F(\mathsf{x})| \leq u \Big]
-
\mathbb{P}^\star\Big[\sup_{\mathsf{x}\in\mathcal{I}}|\hat{\mathfrak{B}}_F(\mathsf{x})| \leq u \Big] \right| = O_{\mathbb{P}}\left(\left(\frac{\log^{5} n}{nh^2}\right)^{\frac{1}{4}} \right).
\end{align*}
\vskip-2em\qed
\end{thm}
\clearpage
{\footnotesize
\section{Proofs}
\subsection{Proof of Theorem \ref{thm:local projection: asymptotic normality}}
\subsubsection*{Part (i)}
The bias term can be bounded by
\begin{align*}
\left|\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[F(\mathsf{x} + hu) - \theta^\prime R(u)\Upsilon_h^{-1}\Big] K\left(u\right) \mathrm{d} u\right| &\leq \sup_{u\in[-1,1]}\Big|F(\mathsf{x} + hu) - \theta^\prime R(u)\Upsilon_h^{-1}\Big|\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} |R(u)| K\left(u\right) \mathrm{d} u\\
&= \varrho(h)\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} |R(u)| K\left(u\right) \mathrm{d} u.
\end{align*}
\subsubsection*{Part (ii)}
The variance can be found as
\begin{align*}
&\ \mathbb{V}\left[ \frac{1}{\sqrt{n}}\sum_{i=1}^n \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u \right]\\
=&\ \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(v)^\prime K(u)K(v) \Big[ F(\mathsf{x}+h(u\wedge v)) - F(\mathsf{x}+hu)F(\mathsf{x}+hv) \Big] g(\mathsf{x}+hu)g(\mathsf{x}+hv)\mathrm{d} u\mathrm{d} v.
\end{align*}
To establish asymptotic normality, we verify the Lyapunov condition with a fourth moment calculation. Take $c$ to be a nonzero vector of conformable dimension, and we employ the Cramer-Wold device:
\begin{align*}
\frac{1}{n}\left(c^\prime \Sigma_h c\right)^{-2}\mathbb{E}\left[ \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} c^\prime R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u\right)^4\right].
\end{align*}
If $c^\prime \Sigma_h c$ is bounded away from zero as the bandwidth decreases, the above will have order $n^{-1}$, as $K(\cdot)$ is bounded and compactly supported and $R(\cdot)$ is locally bounded. Therefore, the Lyapunov condition holds in this case. The more challenging case is when $c^\prime \Sigma_h c$ is of order $h$. In this case, it implies
\begin{align*}
F(\mathsf{x})(1-F(\mathsf{x})) \left|\iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} c^\prime R(u) K(u)g_u \mathrm{d} u\right|^2 = O(h).
\end{align*}
Now consider the fourth moment. The leading term is
\begin{align*}
F(\mathsf{x})(1-F(\mathsf{x}))(3F(\mathsf{x})^2-3F(\mathsf{x})+1)\left|\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} c^\prime R(u) K(u) g(\mathsf{x}+hu)\mathrm{d} u\right|^4 = O(h),
\end{align*}
meaning that for the Lyapunov condition to hold, we need the requirement that $nh\to \infty$.
\subsubsection*{Part (iii)}
This follows immediately from Part (i) and (ii).
\subsection{Proof of Theorem \ref{thm:local projection: standard error pointwise}}
To study the property of $\hat\Sigma_{h}$, we make the following decomposition:
\begin{align*}
\tag{I}\hat\Sigma_{h} &= \frac{1}{n}\sum_{i=1}^n\iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(v)^\prime\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - {F}(\mathsf{x} + hu)\Big]\Big[\mathds{1}(x_i\leq \mathsf{x} + hv) - {F}(\mathsf{x} + hv)\Big] K(u)K(v) g(\mathsf{x}+hu)g(\mathsf{x}+hv)\mathrm{d} u\mathrm{d} v\\
\tag{II}&- \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(v)^\prime\Big[\hat{F}( \mathsf{x} + hu) - {F}(\mathsf{x} + hu)\Big]\Big[\hat{F}(\mathsf{x} + hv) - {F}(\mathsf{x} + hv)\Big] K(u)K(v) g(\mathsf{x}+hu)g(\mathsf{x}+hv)\mathrm{d} u\mathrm{d} v.
\end{align*}
First, it is obvious that term (II) is of order $O_{\mathbb{P}}(1/n)$. Term (I) requires more delicate analysis. Let $c$ be a vector of unit length and suitable dimension, and define
\begin{align*}
c_i &= \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} c^\prime R(u)R(v)^\prime c\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - {F}(\mathsf{x} + hu)\Big]\Big[\mathds{1}(x_i\leq \mathsf{x} + hv) - {F}(\mathsf{x} + hv)\Big] K(u)K(v) g(\mathsf{x}+hu)g(\mathsf{x}+hv)\mathrm{d} u\mathrm{d} v.
\end{align*}
Then
\begin{align*}
c^\prime \text{(I)} c &= \mathbb{E}[c^\prime \text{(I)} c] + O_{\mathbb{P}}\left(\sqrt{\mathbb{V}[c^\prime \text{(I)} c]}\right) = \mathbb{E}[c_i] + O_{\mathbb{P}}\left(\sqrt{\frac{1}{n}\Big(\mathbb{E}[c_i^2] - (\mathbb{E}[c_i])^2 \Big)}\right),
\end{align*}
which implies that
\begin{align*}
\frac{c^\prime \text{(I)} c}{\mathbb{E}[c^\prime \text{(I)} c]} - 1 &= O_{\mathbb{P}}\left(\sqrt{\frac{1}{n}\left(\frac{\mathbb{E}[c_i^2]}{(\mathbb{E}[c_i])^2} - 1 \right)}\right).
\end{align*}
With the same argument used in the proof of Theorem \ref{thm:local projection: asymptotic normality}, one can show that
\begin{align*}
\frac{\mathbb{E}[c_i^2]}{(\mathbb{E}[c_i])^2} = O\left(\frac{1}{h}\right),
\end{align*}
which implies
\begin{align*}
\frac{c^\prime \text{(I)} c}{c^\prime\Sigma_{h} c} - 1 &= O_{\mathbb{P}}\left( \sqrt{\frac{1}{nh}}\right).
\end{align*}
\subsection{Proof of Theorem \ref{thm:local regression: asymptotic normality}}
\subsubsection*{Part (i)}
For the ``denominator,'' its variance is bounded by
\begin{align*}
&\ \left|\mathbb{V}\left[\frac{1}{n}\sum_{i=1}^n\Upsilon_hR(x_i-\mathsf{x})R(x_i-\mathsf{x})^\prime \Upsilon_h \frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right)\right]\right|
\leq \frac{1}{n}\mathbb{E}\left[ \left|\Upsilon_hR(x_i-\mathsf{x})R(x_i-\mathsf{x})^\prime \Upsilon_h\right|^2 \frac{1}{h^2}K\left(\frac{x_i-\mathsf{x}}{h}\right)^2 \right]\\
=&\ \frac{1}{n}\int_{\mathcal{X}} \left|\Upsilon_hR(u-\mathsf{x})R(u-\mathsf{x})^\prime \Upsilon_h\right|^2 \frac{1}{h^2}K\left(\frac{u-\mathsf{x}}{h}\right)^2 f(u)\mathrm{d} u
= \frac{1}{nh}\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} \left|R(u)R(u)^\prime \right|^2 K\left(u\right)^2 f(\mathsf{x}+hu)\mathrm{d} u\\
=&\ O\left(\frac{1}{nh}\right).
\end{align*}
Therefore, under the assumption that $h\to 0$ and $nh\to \infty$, we have
\begin{align*}
\left|\frac{1}{n}\sum_{i=1}^n\Upsilon_hR(x_i-\mathsf{x})R(x_i-\mathsf{x})^\prime \Upsilon_h \frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right)\ -\ \Gamma_h\right| = O_{\mathbb{P}}\left(\sqrt{\frac{1}{nh}} \right),
\end{align*}
which further implies that
\begin{align*}
\hat\theta - \theta &= \Upsilon_h\Gamma_h^{-1}\left( \frac{1}{n}\sum_{i=1}^n \Upsilon_hR(x_i-\mathsf{x})\Big[ \hat{F}(x_i) -R(x_i-\mathsf{x})^\prime\theta_0 \Big] \frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right)\right)\big(1 + o_{\mathbb{P}}(1)\big).
\end{align*}
\subsubsection*{Part (ii)}
The order of the leave-in bias is clearly $1/n$. For the approximation bias \eqref{eq:approximation bias}, we obtained its mean in the proof of Theorem \ref{thm:local projection: asymptotic normality} by setting $G=F$, which has an order of $\varrho(h)$. The approximation bias has a variance of order
\begin{align*}
&\ \left|\mathbb{V}\left[\frac{1}{n}\sum_{j=1}^n \Upsilon_hR(x_j-\mathsf{x})\Big[ F(x_j) -R(x_j-\mathsf{x})^\prime\theta_0 \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right)\right]\right|\\
\leq&\ \frac{1}{n}\mathbb{E}\left[ \left|\Upsilon_hR(x_j-\mathsf{x})\Big[ F(x_j) -R(x_j-\mathsf{x})^\prime\theta_0 \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right)\right|^2\right]\\
=&\ \frac{1}{n}\int_{\mathcal{X}}\left|\Upsilon_hR(u-\mathsf{x})\Big[ F(u) -R(u-\mathsf{x})^\prime\theta_0 \Big]\right|^2 \frac{1}{h^2}K\left(\frac{u-\mathsf{x}}{h}\right)^2f(u)\mathrm{d} u\\
=&\ \frac{1}{nh}\int_{\frac{\mathcal{X}-\mathsf{x}}{h}}\left|R(u)\Big[ F(\mathsf{x} + hu) -R(u-\mathsf{x})^\prime\theta_0 \Big]\right|^2 K(u)^2 f(\mathsf{x} + hu)\mathrm{d} u\\
\leq&\ \frac{1}{nh}\varrho(h)^2\int_{\frac{\mathcal{X}-\mathsf{x}}{h}}|R(u)|^2 K(u)^2 f(\mathsf{x} + hu)\mathrm{d} u
=\ O\left(\frac{\rho(h)^2}{nh}\right).
\end{align*}
Therefore,
\begin{align*}
\text{\eqref{eq:approximation bias}} &= O_{\mathbb{P}}\left( \varrho(h) +\varrho(h)\sqrt{\frac{1}{nh}} \right) = O_{\mathbb{P}}(\varrho(h)),
\end{align*}
provided that $nh\to \infty$.
\subsubsection*{Part (iii)}
We compute the variance of the U-statistic \eqref{eq:quadratic variance}. For simplicity, define
\begin{align*}
u_{ij} &= \Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -F(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right) - \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) f(\mathsf{x} + hu)\mathrm{d} u,
\end{align*}
which satisfies $\mathbb{E}[u_{ij}]=\mathbb{E}[u_{ij}|x_i]=\mathbb{E}[u_{ij}|x_j]=0$. Therefore
\begin{align*}
\mathbb{V}\left[ \text{\eqref{eq:quadratic variance}} \right] =& \frac{1}{n^4}\sum_{i,j=1,i\neq j}^n\sum_{i',j'=1,i\neq j}^n \mathbb{E}\left[u_{ij}u_{i'j'}^\prime\right]
= \frac{1}{n^4}\sum_{i,j=1,i\neq j}^n \mathbb{E}\left[u_{ij}u_{ij}^\prime\right] + \mathbb{E}\left[u_{ij}u_{ji}^\prime\right],
\end{align*}
meaning that
\begin{align*}
\text{\eqref{eq:quadratic variance}} &= O_{\mathbb{P}}\left( \sqrt{\frac{1}{n^2h}} \right).
\end{align*}
\subsubsection*{Part (iv)}
This follows immediately from Part (i)--(iii) and Theorem \ref{thm:local projection: asymptotic normality}.
\subsection{Proof of Theorem \ref{thm:local regression: standard error pointwise}}
We first decompose $\hat{\Sigma}_{h}$ into two terms,
\begin{align*}
\text{(I)} &= \frac{1}{n^3}\sum_{i,j,k=1}^n \Upsilon_{h}R_jR_k^\prime \Upsilon_{h} W_jW_k \Big(\mathds{1}(x_i\leq x_j) - F(x_j)\Big) \Big(\mathds{1}(x_i\leq x_k) - F(x_k)\Big)\\
\text{(II)} &= -\frac{1}{n^2}\sum_{j,k=1}^n \Upsilon_{h}R_jR_k^\prime \Upsilon_{h} W_jW_k \Big(\hat{F}(x_j)-{F}(x_j)\Big) \Big(\hat{F}(x_k) - {F}(x_k)\Big),
\end{align*}
where we use $R_i = R(x_i-\mathsf{x})$ and $W_i = K((x_i-\mathsf{x})/h)/h$ to conserve space.
(II) satisfies
\begin{align*}
\left|\text{(II)}\right| \leq \sup_{x}|\hat{F}(x) - F(x)|^2\frac{1}{n^2}\sum_{j,k=1}^n \left| \Upsilon_{h}R_jR_k^\prime \Upsilon_{h} W_jW_k\right|.
\end{align*}
It is obvious that
\begin{align*}
\sup_{x}|\hat{F}(x) - F(x)|^2 = O_{\mathbb{P}}\left( \frac{1}{n}\right).
\end{align*}
As for the second part, we have
\begin{align*}
\frac{1}{n^2}\sum_{j,k=1}^n\mathbb{E}\left[ \left| \Upsilon_{h}R_jR_k^\prime \Upsilon_{h} W_jW_k\right| \right]
=& \frac{n-1}{n}\mathbb{E}\Big[ \left| \Upsilon_{h}R_jR_k^\prime \Upsilon_{h} W_jW_k\right|\ \Big| j\neq k \Big] + \frac{1}{n}\mathbb{E}\Big[ \left| \Upsilon_{h}R_kR_k^\prime \Upsilon_{h} W_kW_k\right| \Big]\\
=& O_{\mathbb{P}}\left( 1 + \frac{1}{nh}\right) = O_{\mathbb{P}}\left(1\right),
\end{align*}
which holds as long as $nh\to \infty$. Then it further implies that
\begin{align*}
\text{(II)} = O_{\mathbb{P}}\left(\frac{1}{n}\right).
\end{align*}
To analyze (I), we further expand this term into ``diagonal'' and ``non-diagonal'' sums:
\begin{align*}
\tag{I.1}\text{(I)} &= \frac{1}{n^3}\sum_{\substack{i,j,k=1 \\ \text{distinct}}}^n \Upsilon_{h}R_jR_k^\prime \Upsilon_{h} W_jW_k \Big(\mathds{1}(x_i\leq x_j) - F(x_j)\Big) \Big(\mathds{1}(x_i\leq x_k) - F(x_k)\Big)\\
\tag{I.2}&\quad + \frac{1}{n^3}\sum_{\substack{i,k=1 \\ \text{distinct}}}^n \Upsilon_{h}R_iR_k^\prime \Upsilon_{h} W_iW_k \Big(\mathds{1}(x_i\leq x_i) - F(x_i)\Big) \Big(\mathds{1}(x_i\leq x_k) - F(x_k)\Big)\\
\tag{I.3}&\quad + \frac{1}{n^3}\sum_{\substack{i,j=1 \\ \text{distinct}}}^n \Upsilon_{h}R_jR_i^\prime \Upsilon_{h} W_jW_i \Big(\mathds{1}(x_i\leq x_j) - F(x_j)\Big) \Big(\mathds{1}(x_i\leq x_i) - F(x_i)\Big)\\
\tag{I.4}&\quad + \frac{1}{n^3}\sum_{\substack{i,j=1 \\ \text{distinct}}}^n \Upsilon_{h}R_jR_j^\prime \Upsilon_{h} W_jW_j \Big(\mathds{1}(x_i\leq x_j) - F(x_j)\Big) \Big(\mathds{1}(x_i\leq x_j) - F(x_j)\Big)\\
\tag{I.5}&\quad + \frac{1}{n^3}\sum_{i} \Upsilon_{h}R_iR_i^\prime \Upsilon_{h} W_iW_i \Big(\mathds{1}(x_i\leq x_i) - F(x_i)\Big) \Big(\mathds{1}(x_i\leq x_i) - F(x_i)\Big).
\end{align*}
By calculating the expectation of the absolute value of the summands above, it is straightforward to show
\begin{align*}
\text{(I.2)} = O_{\mathbb{P}}\left(\frac{1}{n}\right),\quad
\text{(I.3)} = O_{\mathbb{P}}\left(\frac{1}{n}\right),\quad
\text{(I.4)} = O_{\mathbb{P}}\left(\frac{1}{nh}\right),\quad
\text{(I.5)} = O_{\mathbb{P}}\left(\frac{1}{n^2h}\right).
\end{align*}
Therefore, we have
\begin{align*}
\hat{\Sigma}_{h} &= \text{(I.1)} + O_{\mathbb{P}}\left(\frac{1}{nh}\right)
= \frac{1}{n^3}\sum_{\substack{i,j,k=1 \\ \text{distinct}}}^n \Upsilon_{h}R_jR_k^\prime \Upsilon_{h} W_jW_k \Big(\mathds{1}(x_i\leq x_j) - F(x_j)\Big) \Big(\mathds{1}(x_i\leq x_k) - F(x_k)\Big) + O_{\mathbb{P}}\left(\frac{1}{nh}\right).
\end{align*}
To proceed, define
\begin{align*}
u_{ij} &= \Upsilon_{h}R_j W_j \Big(\mathds{1}(x_i\leq x_j) - F(x_j)\Big)\quad \text{and}\quad \bar{u}_i = \mathbb{E}[u_{ij}|x_i;i\neq j].
\end{align*}
Then we can further decompose (I.1) into
\begin{align*}
\text{(I.1)} &= \frac{1}{n^3}\sum_{\substack{i,j,k=1 \\ \text{distinct}}}^n u_{ij}u_{ik}^\prime
= \frac{1}{n^3}\sum_{\substack{i,j,k=1 \\ \text{distinct}}}^n \mathbb{E}[u_{ij}u_{ik}^\prime|x_i]
\quad + \frac{1}{n^3}\sum_{\substack{i,j,k=1 \\ \text{distinct}}}^n \Big(u_{ij}u_{ik}^\prime - \mathbb{E}[u_{ij}u_{ik}^\prime|x_i]\Big)\\
&= \underbrace{\frac{(n-1)(n-2)}{n^3}\sum_{i=1}^n \bar{u}_i\bar{u}_i^\prime}_{\textstyle \text{(I.1.1)}}
\ +\ \underbrace{\frac{1}{n^3}\sum_{\substack{i,j,k=1 \\ \text{distinct}}}^n \Big(u_{ij}u_{ik}^\prime - \bar{u}_i\bar{u}_i^\prime\Big)}_{\textstyle \text{(I.1.2)}}.
\end{align*}
We have already analyzed (I.1.1) in Theorem \ref{thm:local projection: standard error pointwise}, which suggests
\begin{align*}
\text{(I.1.1)}&= \Sigma_{h} + O_{\mathbb{P}}\left( \frac{1}{\sqrt{n}} \right).
\end{align*}
Now we study (I.1.2), which satisfies
\begin{align*}
\text{(I.1.2)} &= \underbrace{\frac{n-2}{n^3}\sum_{\substack{i,j=1 \\ \text{distinct}}}^n \Big(u_{ij} - \bar{u}_i\Big)\bar{u}_i^\prime}_{\textstyle \text{(I.1.2.1)}}
\ +\ \underbrace{\frac{n-2}{n^3}\sum_{\substack{i,j=1 \\ \text{distinct}}}^n \bar{u}_i\Big(u_{ij} - \bar{u}_i\Big)^\prime}_{\textstyle \text{(I.1.2.2)}}
\ +\ \underbrace{\frac{1}{n^3}\sum_{\substack{i,j,k=1 \\ \text{distinct}}}^n \Big(u_{ij} - \bar{u}_i\Big)\Big(u_{ik} - \bar{u}_i\Big)^\prime}_{\textstyle \text{(I.1.2.3)}}.
\end{align*}
With variance calculation, it is easy to see that
\begin{align*}
\text{(I.1.2.3)} &= O_{\mathbb{P}}\left( \frac{1}{nh} \right).
\end{align*}
Therefore we have
\begin{align*}
\frac{c^\prime (\hat\Sigma_{h} - \Sigma_{h}) c}{c^\prime \Sigma_{h} c} &= O_{\mathbb{P}}\left( \frac{1}{\sqrt{nh^2}} \right) + 2\frac{c^\prime \text{(I.1.2.1)} c}{c^\prime \Sigma_{h} c},
\end{align*}
since (I.1.2.1) and (I.1.2.2) are transpose of each other. To close the proof, we calculate the variance of the last term in the above.
\begin{align*}
\mathbb{V}\left[ \frac{c^\prime \text{(I.1.2.1)} c}{c^\prime \Sigma_{h} c} \right] &= \frac{1}{(c^\prime \Sigma_{h} c)^2}\frac{(n-2)^2}{n^6}\mathbb{E}\left[\sum_{\substack{i,j=1 \\ \text{distinct}}}^n\sum_{\substack{i',j'=1 \\ \text{distinct}}}^n c^\prime\Big(u_{ij} - \bar{u}_i\Big)\bar{u}_i^\prime c c^\prime\Big(u_{i'j'} - \bar{u}_{i'}\Big)\bar{u}_{i'}^\prime c \right]\\
&= \frac{1}{(c^\prime \Sigma_{h} c)^2}\frac{(n-2)^2}{n^6}\mathbb{E}\left[\sum_{\substack{i,j,i'=1 \\ \text{distinct}}}^n c^\prime u_{ij}\bar{u}_i^\prime c c^\prime u_{i'j}\bar{u}_{i'}^\prime c \right] + \text{higher order terms}.
\end{align*}
The expectation is further given by (note that $i$, $j$ and $i'$ are assumed to be distinct indices)
\begin{align*}
&\ \mathbb{E}\left[ c^\prime u_{ij}\bar{u}_i^\prime c c^\prime u_{i'j}\bar{u}_{i'}^\prime c \right]\\
=&\ \mathbb{E}\iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} W_j^2\left[c^\prime \Upsilon_h R_jR(u) c c^\prime \Upsilon_h R_jR(v) c\right] K(u) K(v) \\
&\qquad \left[ F(x_j\wedge (\mathsf{x} + hu)) - F(x_j)F(\mathsf{x} + hu) \right]\left[ F(x_j\wedge (\mathsf{x} + hv)) - F(x_j)F(\mathsf{x} + hv) \right]f(\mathsf{x} + hu)f(\mathsf{x} + hv)\mathrm{d} u \mathrm{d} v\\
=&\ \frac{1}{h}\iiint_{\frac{\mathcal{X}-\mathsf{x}}{h}} \left[c^\prime R(w)R(u) c c^\prime R(w)R(v) c\right] K(u) K(v) K(w)^2 \\
&\qquad \left[ F(\mathsf{x} + h(w\wedge u)) - F(\mathsf{x} + hw)F(\mathsf{x} + hu) \right]\left[ F(\mathsf{x} + h(w\wedge v)) - F(\mathsf{x} + hw)F(\mathsf{x} + hv) \right]f(\mathsf{x} + hw)f(\mathsf{x} + hu)f(\mathsf{x} + hv)\mathrm{d} w\mathrm{d} u \mathrm{d} v\\
=&\ \frac{1}{h}F(\mathsf{x})^2(1-F(\mathsf{x}))^2\iiint_{\frac{\mathcal{X}-\mathsf{x}}{h}} \left[c^\prime R(w)R(u) c c^\prime R(w)R(v) c\right] K(u) K(v) K(w)^2 f_wf_uf_v\mathrm{d} w\mathrm{d} u \mathrm{d} v + \text{higher-order terms}.
\end{align*}
If $c^\prime \Sigma_{h} c = O(1)$, then the above will have order $h$, which means
\begin{align*}
\mathbb{V}\left[ \frac{c^\prime \text{(I.1.2.1)} c}{c^\prime \Sigma_{h} c} \right] &= O\left(\frac{1}{nh}\right).
\end{align*}
If $c^\prime \Sigma_{h} c = O(h)$, however, $\mathbb{E}\left[ c^\prime u_{ij}\bar{u}_i^\prime c c^\prime u_{i'j}\bar{u}_{i'}^\prime c \right]$ will be $O(1)$, which will imply that
\begin{align*}
\mathbb{V}\left[ \frac{c^\prime \text{(I.1.2.1)} c}{c^\prime \Sigma_{h} c} \right] &= O\left(\frac{1}{nh^2}\right).
\end{align*}
As a result, we have
\begin{align*}
\frac{c^\prime (\hat\Sigma_{h} - \Sigma_{h}) c}{c^\prime \Sigma_{h} c} &= O_{\mathbb{P}}\left( \frac{1}{\sqrt{nh^2}} \right).
\end{align*}
Now consider
\begin{align*}
\frac{c^\prime \hat\Gamma_{h}^{-1} \hat\Sigma_{h} \hat\Gamma_{h}^{-1} c}{c^\prime \Gamma_{h}^{-1} \Sigma_{h} \Gamma_{h}^{-1} c} - 1 &= \frac{c^\prime \hat\Gamma_{h}^{-1} (\hat\Sigma_{h} - \Sigma_h) \hat\Gamma_{h}^{-1} c}{c^\prime \Gamma_{h}^{-1} \Sigma_{h} \Gamma_{h}^{-1} c}
+ \frac{c^\prime (\hat\Gamma_{h}^{-1} -\Gamma_{h}^{-1}) \Sigma_h \hat\Gamma_{h}^{-1} c}{c^\prime \Gamma_{h}^{-1} \Sigma_{h} \Gamma_{h}^{-1} c}
+ \frac{c^\prime (\hat\Gamma_{h}^{-1} -\Gamma_{h}^{-1}) \Sigma_h \Gamma_{h}^{-1} c}{c^\prime \Gamma_{h}^{-1} \Sigma_{h} \Gamma_{h}^{-1} c}\\
&= \frac{c^\prime \hat\Gamma_{h}^{-1} (\hat\Sigma_{h} - \Sigma_h) \hat\Gamma_{h}^{-1} c}{c^\prime \Gamma_{h}^{-1} \Sigma_{h} \Gamma_{h}^{-1} c}
+ 2\frac{c^\prime (\hat\Gamma_{h}^{-1} -\Gamma_{h}^{-1}) \Sigma_h \Gamma_{h}^{-1} c}{c^\prime \Gamma_{h}^{-1} \Sigma_{h} \Gamma_{h}^{-1} c}
+ \frac{c^\prime (\hat\Gamma_{h}^{-1} -\Gamma_{h}^{-1}) \Sigma_h (\hat\Gamma_{h}^{-1} -\Gamma_{h}^{-1}) c}{c^\prime \Gamma_{h}^{-1} \Sigma_{h} \Gamma_{h}^{-1} c}
.
\end{align*}
From the analysis of $\hat\Sigma_h$, we have
\begin{align*}
\frac{c^\prime \hat\Gamma_{h}^{-1} (\hat\Sigma_{h} - \Sigma_h) \hat\Gamma_{h}^{-1} c}{c^\prime \Gamma_{h}^{-1} \Sigma_{h} \Gamma_{h}^{-1} c} = O_{\mathbb{P}}\left( \frac{1}{\sqrt{nh^2}} \right).
\end{align*}
For the second term, we have
\begin{align*}
\left|\frac{c^\prime (\hat\Gamma_{h}^{-1} -\Gamma_{h}^{-1}) \Sigma_h \hat\Gamma_{h}^{-1} c}{c^\prime \Gamma_{h}^{-1} \Sigma_{h} \Gamma_{h}^{-1} c}\right|
&\leq \frac{|c^\prime (\hat\Gamma_{h}^{-1} -\Gamma_{h}^{-1}) \Sigma_h^{1/2}|\cdot|c^\prime \Gamma_{h}^{-1} \Sigma_h^{1/2}|}{|c^\prime \Gamma_{h}^{-1}\Sigma_h^{1/2} |^2}\\
&= \frac{|c^\prime (\hat\Gamma_{h}^{-1} -\Gamma_{h}^{-1}) \Sigma_h^{1/2}|}{|c^\prime \Gamma_{h}^{-1}\Sigma_h^{1/2} |} = O_{\mathbb{P}}\left( \sqrt{\frac{1}{nh^2}} \right).
\end{align*}
The third term has order
\begin{align*}
\frac{c^\prime (\hat\Gamma_{h}^{-1} -\Gamma_{h}^{-1}) \Sigma_h (\hat\Gamma_{h}^{-1} -\Gamma_{h}^{-1}) c}{c^\prime \Gamma_{h}^{-1} \Sigma_{h} \Gamma_{h}^{-1} c} = O_{\mathbb{P}}\left( \frac{1}{nh^2} \right).
\end{align*}
\subsection{Proof of Corollary \ref{coro:asy normal loc pol projection estimator}}
This follows directly from Theorem \ref{thm:local projection: asymptotic normality}.
\subsection{Proof of Corollary \ref{coro:asy normal ortho loc pol projection estimator}}
To understand \eqref{eq:local orthogonal polynomial projection estimator}, note that
\begin{align*}
\hat\theta_{F}^\perp &= \left( \int_{\mathcal{X}} \Lambda_h^\prime R(u-\mathsf{x})R(u-\mathsf{x})^\prime\Lambda_h \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} F(u) \right)^{-1}\left( \int_{\mathcal{X}} \Lambda_h^\prime R(u-\mathsf{x})\hat{F}(u) \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} F(u) \right)\\
&= \Lambda_h^{-1}\left( \int_{\mathcal{X}} R(u-\mathsf{x})R(u-\mathsf{x})^\prime \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} F(u) \right)^{-1}\left( \int_{\mathcal{X}} R(u-\mathsf{x})\hat{F}(u) \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} F(u) \right),
\end{align*}
which means $\hat\theta_{F}^\perp = \Lambda_h^{-1}\hat\theta_{F}$. Then we have (up to an approximation bias term)
\begin{align*}
&\ \hat\theta_{F}^\perp - \Lambda_h^{-1}\theta_0 = \Lambda_h^{-1}(\hat\theta_{F} - \theta_0) \\
=&\ \Lambda_h^{-1}\left( \int_{\mathcal{X}} R(u-\mathsf{x})R(u-\mathsf{x})^\prime \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} F(u) \right)^{-1}\left( \int_{\mathcal{X}} R(u-\mathsf{x})\Big(\hat{F}(u) - F(u)\Big) \frac{1}{h}K\left(\frac{u-\mathsf{x}}{h}\right) \mathrm{d} F(u) \right)\\
=&\ \Lambda_h^{-1}\Upsilon_h\left( \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)R(u)^\prime K\left(u\right) f(\mathsf{x}+hu)\mathrm{d} u \right)^{-1}\left( \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big(\hat{F}(\mathsf{x}+hu) - F(\mathsf{x}+hu) \Big) K\left(u\right) f(\mathsf{x}+hu)\mathrm{d} u \right)\\
=&\ \Lambda_h^{-1}\Upsilon_h\Lambda_h\left( \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R^\perp(u)R^\perp(u)^{ \prime} K\left(u\right) f(\mathsf{x}+hu)\mathrm{d} u \right)^{-1}\left( \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R^\perp(u)\Big(\hat{F}(\mathsf{x}+hu) - F(\mathsf{x}+hu) \Big) K\left(u\right) f(\mathsf{x}+hu)\mathrm{d} u \right).
\end{align*}
We first discuss the transformed parameter vector $\Lambda_h^{-1}\theta_0$. By construction, the matrix $\Lambda_h$ takes the following form:
\begin{align*}
\Lambda_h &= \begin{bmatrix}
1 & c_{1,2} & c_{1,3} & \cdots &c_{1,p+2} \\
0 & 1 & 0 & \cdots &c_{2,p+2} \\
0 & 0 & 1 & \cdots &c_{2,p+2} \\
\vdots & \vdots & \vdots & \ddots &\vdots \\
0 & 0 & 0 & \cdots & 1
\end{bmatrix}
\end{align*}
where $c_{i,j}$ are some constants (possibly depending on $h$). Therefore, the above matrix differs from the identity matrix only in its first row and in the last column. This observation also holds for $\Lambda_h^{-1}$. Since the last component of $\theta_0$ is zero (because the extra regressor $Q_h(\cdot)$ is redundant), we conclude that, except for the first element, $\Lambda_h\theta$ and $\theta$ are identical. More specifically, let $I_{-1}$ be the identity matrix excluding the first row:
\begin{align*}
I_{-1} &= \begin{bmatrix}
0 & 1 & 0 & 0 & \cdots & 0\\
0 & 0 & 1 & 0 & \cdots & 0\\
0 & 0 & 0 & 1 & \cdots & 0\\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & 0 & \cdots & 1\\
\end{bmatrix},
\end{align*}
which is used to extract all elements of a vector except for the first one, then by Theorem \ref{thm:local projection: asymptotic normality},
\begin{align*}
\sqrt{n}\left(I_{-1}(\Lambda_h^{-1}\Upsilon_h\Lambda_h)(\Gamma_h^\perp)^{-1}\Sigma_h^\perp(\Gamma_h^\perp)^{-1}(\Lambda_h^{-1}\Upsilon_h\Lambda_h)^\prime I_{-1}^\prime\right)^{-1/2} \begin{bmatrix}
\hat\theta^\perp_{P,F} - \theta_{P} \\
\hat\theta^\perp_{Q,F}
\end{bmatrix} \rightsquigarrow \mathcal{N}(0, I),
\end{align*}
where $\theta^\perp_{P,F}$ contains the second to the $p+1$-th element of $\theta^\perp_{F}$, and $\theta^\perp_{Q,F}$ is the last element.
Now we discuss the covariance matrix in the above display. Due to orthogonalization, $\Gamma_h^\perp$ is block diagonal. To be precise,
\begin{align*}
\Gamma_h^\perp = f(\mathsf{x})
\begin{bmatrix}
\Gamma_{1,h}^\perp & 0 & 0\\
& \Gamma_{P,h}^\perp & 0\\
0 & 0 & \Gamma_{Q,h}^\perp
\end{bmatrix},\ \Gamma_{1,h}^\perp = \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)\mathrm{d} u,\ \Gamma_{P,h}^\perp = \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} P^\perp(u)P^\perp(u)^\prime K(u)\mathrm{d} u,\ \Gamma_{Q,h}^\perp = \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} Q^\perp(u)^2 K(u)\mathrm{d} u.
\end{align*}
Finally, using the structure of $\Lambda_h$ and $\Upsilon_h$, we have
\begin{align*}
I_{-1}(\Lambda_h^{-1}\Upsilon_h\Lambda_h)(\Gamma_h^\perp)^{-1} &= I_{-1}\Upsilon_h(\Gamma_h^\perp)^{-1}.
\end{align*}
The form of $\Sigma_h^\perp$ is quite involved, but with some algebra, and using the fact that the basis $R(\cdot)$ (or $R^\perp(\cdot)$) includes a constant and polynomials, one can show the following:
\begin{align*}
(\Lambda_h^{-1}\Upsilon_h\Lambda_h)(\Gamma_h^\perp)^{-1}\Sigma_h^\perp(\Gamma_h^\perp)^{-1}(\Lambda_h^{-1}\Upsilon_h\Lambda_h)^\prime = hf(\mathsf{x})\Upsilon_{-1, h}(\Gamma_{-1,h}^\perp)^{-1}\Sigma_{-1,h}^\perp(\Gamma_{-1,h}^\perp)^{-1}\Upsilon_{-1, h},
\end{align*}
where $\Upsilon_{-1, h}$, $\Gamma_{-1,h}^\perp$ and $\Sigma_{-1,h}^\perp$ are obtained by excluding the first row and the first column of $\Upsilon_{h}$, $\Gamma_{h}^\perp$ and $\Sigma_{h}^\perp$, respectively:
\begin{align*}
\Upsilon_{-1, h} &= \begin{bmatrix}
h^{-1} & 0 & 0 & \cdots & 0\\
0 & h^{-2} & 0 & \cdots & 0\\
0 & 0 & h^{-3} & \cdots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & \cdots & \upsilon_h\\
\end{bmatrix},\ \Gamma_{-1,h}^\perp = f(\mathsf{x})\begin{bmatrix}
\Gamma_{P,h}^\perp & 0\\
0 & \Gamma_{Q,h}^\perp
\end{bmatrix},\ \Sigma_{-1,h}^\perp = f(\mathsf{x})^3\begin{bmatrix}
\Sigma_{PP,h}^\perp & \Sigma_{PQ,h}^\perp \\
\Sigma_{QP,h}^\perp & \Sigma_{QQ,h}^\perp
\end{bmatrix},
\end{align*}
and
\begin{align*}
\Sigma_{PP,h}^\perp &= \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)K(v)P^\perp(u)P^\perp(v)^\prime (u\wedge v) \mathrm{d} u\mathrm{d} v,\qquad \Sigma_{QQ,h}^\perp &= \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)K(v)Q^\perp(u)Q^\perp(v) (u\wedge v) \mathrm{d} u\mathrm{d} v\\
\Sigma_{PQ,h}^\perp &= (\Sigma_{QP,h}^\perp)^\prime = \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)K(v)P^\perp(u)Q^\perp(v) (u\wedge v) \mathrm{d} u\mathrm{d} v.
\end{align*}
\subsection{Proof of Lemma \ref{lem:survival transformation} }
\subsubsection*{Part (i)}
To start,
\begin{align*}
\int_{\frac{\mathcal{X}-\mathsf{x}}{h} \cap [-1,1]} \mathcal{H}(g_1)(u)\mathcal{H}(g_2)(u) \mathrm{d} u &= \int_{\frac{\mathcal{X}-\mathsf{x}}{h} \cap [-1,1]} \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} \mathds{1}(v_1\geq u)K(v_1)g(v_1)\mathrm{d} v_1\right)\left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} \mathds{1}(v_2\geq u)K(v_2)g(v_2)\mathrm{d} v_2\right) \mathrm{d} u\\
&= \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}}K(v_1)K(v_2)g(v_1)g(v_2)\left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h} \cap [-1,1]} \mathds{1}(v_1\geq u) \mathds{1}(v_2\geq u) \mathrm{d} u\right)\mathrm{d} v_1\mathrm{d} v_2\\
&= \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}}K(v_1)K(v_2)g(v_1)g(v_2)
\left[(v_1\wedge v_2)\wedge \left( \frac{\overline{x}-\mathsf{x}}{h}\wedge 1 \right) - \left( \frac{\underline{x}-\mathsf{x}}{h}\vee (-1) \right)\right]
\mathrm{d} v_1\mathrm{d} v_2\\
&= \iint_{\frac{\mathcal{X}-\mathsf{x}}{h}}K(v_1)K(v_2)g(v_1)g(v_2)(v_1\wedge v_2)\mathrm{d} v_1\mathrm{d} v_2,
\end{align*}
where to show the last equality, we used the fact that $v_1\leq \frac{\overline{x}-\mathsf{x}}{h}\wedge 1$ and $v_2\leq \frac{\overline{x}-\mathsf{x}}{h}\wedge 1$ for the outer double integral.
\subsubsection*{Part (ii)}
For this part,
\begin{align*}
\int_{\frac{\mathcal{X}-\mathsf{x}}{h} \cap [-1,1]} \mathcal{H}(g_1)(u)\dot{g}_2(u) \mathrm{d} u &= \int_{\frac{\mathcal{X}-\mathsf{x}}{h} \cap [-1,1]} \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} \mathds{1}(v\geq u)K(v)g_1(v)\mathrm{d} v\right)\dot{g}_2(u) \mathrm{d} u\\
&= \int_{\frac{\mathcal{X}-\mathsf{x}}{h}}K(v)g_1(v) \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h} \cap [-1,1]}\mathds{1}(v\geq u) \dot{g}_2(u) \mathrm{d} u \right)\mathrm{d} v\\
&= \int_{\frac{\mathcal{X}-\mathsf{x}}{h}}K(v)g_1(v) \left[g_2\left(v\wedge \frac{\overline{x}-\mathsf{x}}{h}\wedge 1 \right) - g_2\left( \frac{\underline{x}-\mathsf{x}}{h}\vee (-1) \right) \right] \mathrm{d} v\\
&= \int_{\frac{\mathcal{X}-\mathsf{x}}{h}}K(v)g_1(v) g_2(v) \mathrm{d} v.
\end{align*}
Again, to show the last equality, we used the fact that $v\leq \frac{\overline{x}-\mathsf{x}}{h}\wedge 1$ for the outer integral.
\subsection{Proof of Theorem \ref{thm:variance bound}}
To find a bound of the maximization problem, we note that for any $c\in\mathbb{R}^{p-1}$, one has
\begin{align*}
\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \mathcal{H}(p_\ell)(u)\mathcal{H}(q)(u)\mathrm{d} u = \int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \Big[\mathcal{H}(p_\ell)(u)+c^\prime \dot{P}(u)\Big]\mathcal{H}(q)(u)\mathrm{d} u,
\end{align*}
due to the constraint. Therefore, an upper bound of the objective function is (due to the Cauchy-Schwartz inequality)
\begin{align*}
&\ \inf_{c} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \Big[\mathcal{H}(p_\ell)(u)+c^\prime \dot{P}(u)\Big]^2\mathrm{d} u\\
=&\ \inf_{c} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \Big[\mathcal{H}(p_\ell)(u)^2 + 2c^\prime \dot{P}(u)\mathcal{H}(p_\ell)(u) + c^\prime \dot{P}(u)\dot{P}(u)^\prime c\Big]\mathrm{d} u\\
=&\ \int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]}\mathcal{H}(p_\ell)(u)^2\mathrm{d} u + \inf_{c} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \Big[2c^\prime \dot{P}(u)\mathcal{H}(p_\ell)(u) + c^\prime \dot{P}(u)\dot{P}(u)^\prime c\Big]\mathrm{d} u\\
=&\ \int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]}\mathcal{H}(p_\ell)(u)^2\mathrm{d} u + \inf_{c}\left[ 2c^\prime \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)P(u)p_{\ell}(u)\mathrm{d} u\right) + c^\prime \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \dot{P}(u)\dot{P}(u)^\prime \mathrm{d} u\right)c \right],
\end{align*}
which is minimized by setting
\begin{align*}
c &= -\left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \dot{P}(u)\dot{P}(u)^\prime \mathrm{d} u\right)^{-1} \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)P(u)p_{\ell}(u)\mathrm{d} u\right).
\end{align*}
As a result, an upper bound of \eqref{eq:maximization} is
\begin{align*}
&\ \int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]}\mathcal{H}(p_\ell)(u)^2\mathrm{d} u - \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)P(u)p_{\ell}(u)\mathrm{d} u\right)^\prime \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \dot{P}(u)\dot{P}(u)^\prime \mathrm{d} u\right)^{-1} \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)P(u)p_{\ell}(u)\mathrm{d} u\right).
\end{align*}
We may further simplify the above. First,
\begin{align*}
\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]}\mathcal{H}(p_\ell)(u)^2\mathrm{d} u &= e_\ell^\prime(\Gamma_{P,h}^\perp)^{-1} \Sigma_{PP,h}^\perp (\Gamma_{P,h}^\perp)^{-1}e_\ell.
\end{align*}
Second, note that
\begin{align*}
\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)P(u)p_{\ell}(u)\mathrm{d} u &= \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)P(u)P^\perp(u)^\prime \mathrm{d} u\right) (\Gamma_{P,h}^\perp)^{-1}e_\ell
= \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} K(u)P^\perp(u)P^\perp(u)^\prime \mathrm{d} u\right) (\Gamma_{P,h}^\perp)^{-1}e_\ell
= e_\ell.
\end{align*}
As a result, an upper bound of \eqref{eq:maximization} is
\begin{align*}
&\ e_\ell^\prime(\Gamma_{P,h}^\perp)^{-1} \Sigma_{PP,h}^\perp (\Gamma_{P,h}^\perp)^{-1}e_\ell - e_\ell^\prime \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \dot{P}(u)\dot{P}(u)^\prime \mathrm{d} u\right)^{-1} e_{\ell}\\
&= e_\ell^\prime\left[ (\Gamma_{P,h}^\perp)^{-1} \Sigma_{PP,h}^\perp (\Gamma_{P,h}^\perp)^{-1} - \left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}\cap[-1,1]} \dot{P}(u)\dot{P}(u)^\prime \mathrm{d} u\right)^{-1} \right]e_\ell.
\end{align*}
\subsection{Additional Preliminary Lemmas}
\begin{lem}\label{lem:uniform convergence of iid sum}
Assume $\{u_{i,h}(a):\ a\in A\subset\mathbb{R}^d\}$ are independent across $i$, and $\mathbb{E}[u_{i,h}(a)]=0$ for all $a\in A$ and all $h > 0$. In addition, assume for each $\varepsilon>0$ there exists $\{u_{i,h,\varepsilon}(a): a\in A\}$, such that
\begin{align*}
|a-b|\leq \varepsilon\quad \Rightarrow\quad |u_{i,h}(a)-u_{i,h}(b)|\leq u_{i,h,\varepsilon}(a).
\end{align*}
Define
\begin{alignat*}{2}
C_1 &= \sup_{a\in A}\max_{1\leq i\leq n}\mathbb{V}[u_{i,h}(a)],\qquad C_2 &&= \sup_{a\in A}\max_{1\leq i\leq n}|u_{i,h}(a)|\\
C_{1,\varepsilon} &= \sup_{a\in A}\max_{1\leq i\leq n}\mathbb{V}[u_{i,h,\varepsilon}(a)],\quad C_{2,\varepsilon} &&= \sup_{a\in A}\max_{1\leq i\leq n}|u_{i,h,\varepsilon}(a) - \mathbb{E}[u_{i,h,\varepsilon}(a)]|,\quad C_{3,\varepsilon}=\sup_{a\in A}\max_{1\leq i\leq n}\mathbb{E}[|u_{i,h,\varepsilon}(a)|].
\end{alignat*}
Then
\begin{align*}
\sup_{a\in A}\left|\frac{1}{n}\sum_{i=1}^n u_{i,h}(a)\right| = O_{\mathbb{P}}\left( \gamma + \gamma_{\varepsilon} + C_{3,\varepsilon}\right),
\end{align*}
where $\gamma$ and $\gamma_{\varepsilon}$ are any sequences satisfying
\begin{align*}
&\frac{\gamma^2n}{(C_{1} + \frac{1}{3} \gamma C_{2})\log N(\varepsilon, A, |\cdot|)}\quad \text{and} \quad
\frac{\gamma_{\varepsilon}^2n}{(C_{1,\varepsilon} + \frac{1}{3} \gamma_\varepsilon C_{2,\varepsilon})\log N(\varepsilon, A, |\cdot|)}\quad \text{are bounded from below},
\end{align*}
and $N(\varepsilon, A, |\cdot|)$ is the covering number of $A$.
\qed
\end{lem}
\begin{remark}
Provided that $u_{i,h}(\cdot)$ is reasonably smooth, one can always choose $\varepsilon$ (as a function of $n$ and $h$) small enough, and the leading order will be given by $\gamma$ (and hence is determined by $C_1$ and $C_2$). \qed
\end{remark}
\noindent\textbf{Proof.} Let $A_\varepsilon$ be an $\varepsilon$-covering of $A$, then
\begin{align*}
&\ \sup_{a\in A} \left|\frac{1}{n}\sum_{i=1}^n u_{i,h}(a) \right| \leq \sup_{a\in A_\varepsilon} \left|\frac{1}{n}\sum_{i=1}^n u_{i,h}(a) \right|
+\sup_{a\in A_\varepsilon,b\in A,|a-b|\leq \varepsilon}\left|\frac{1}{n}\sum_{i=1}^n u_{i,h}(a)-u_{i,h}(b)\right|.
\end{align*}
Next we apply the union bound and Bernstein's inequality:
\begin{align*}
\mathbb{P}\left[ \sup_{a\in A_\varepsilon} \left|\frac{1}{n}\sum_{i=1}^n u_{i,h}(a) \right| \geq \gamma u\right] &\leq N(\varepsilon, A, |\cdot|) \sup_{a\in A}\mathbb{P}\left[ \left|\frac{1}{n}\sum_{i=1}^n u_{i,h}(a) \right| \geq \gamma u\right]\\
&\leq 2N(\varepsilon, A, |\cdot|)\exp\left\{ -\frac{1}{2}\frac{\gamma^2nu^2}{ C_{1} + \frac{1}{3} \gamma C_{2} u } \right\}\\
&= 2\exp\left\{ -\frac{1}{2}\frac{\gamma^2nu^2}{ C_{1} + \frac{1}{3} \gamma C_{2} u } + \log N(\varepsilon, A, |\cdot|) \right\}.
\end{align*}
Now take $u$ sufficiently large, then the above is further bounded by:
\begin{align*}
\mathbb{P}\left[ \sup_{a\in A_\varepsilon} \left|\frac{1}{n}\sum_{i=1}^n u_{i,h}(a) \right| \geq \gamma u\right] &\leq 2\exp\left\{ -\log N(\varepsilon, A, |\cdot|)\left[\frac{1}{2}\frac{1}{\log N(\varepsilon, A, |\cdot|)}\frac{\gamma^2n}{ C_{1} + \frac{1}{3} \gamma C_{2} }u - 1 \right] \right\} ,
\end{align*}
which tends to zero if $\log N(\varepsilon, A, |\cdot|)\to\infty$ and
\begin{align*}
\frac{\gamma^2n}{(C_{1} + \frac{1}{3} \gamma C_{2})\log N(\varepsilon, A, |\cdot|)}\ \text{is bounded from below},
\end{align*}
in which case we have
\begin{align*}
\sup_{a\in A_\varepsilon} \left|\frac{1}{n}\sum_{i=1}^n u_{i,h}(a) \right| = O_{\mathbb{P}}\left( \gamma \right).
\end{align*}
We can apply the same technique to the other term, and obtain
\begin{align*}
\sup_{a\in A_\varepsilon,b\in A,|a-b|\leq \varepsilon}\left|\frac{1}{n}\sum_{i=1}^n u_{i,h}(a)-u_{i,h}(b)\right| = O_{\mathbb{P}}\left( \gamma_{\varepsilon} \right),
\end{align*}
where $\gamma_{\varepsilon}$ is any sequence satisfying
\begin{align*}
\frac{\gamma_{\varepsilon}^2n}{(C_{1,\varepsilon} + \frac{1}{3} \gamma_\varepsilon C_{2,\varepsilon})\log N(\varepsilon, A, |\cdot|)}\ \text{is bounded from below}.
\end{align*}
\qed
\begin{lem}[Corollary 5.1 in \citealt{chernozhukov2019central}]\label{lem:gaussian comparison}
Let $\mathbf{z}_1,\mathbf{z}_2\in\mathbb{R}^{\ell_n}$ be two mean-zero Gaussian random vectors with covariance matrices $\boldsymbol{\Omega}_{1}$ and $\boldsymbol{\Omega}_{2}$, respectively. Further assume that the diagonal elements in $\boldsymbol{\Omega}_{1}$ are all one. Then
\begin{align*}
\sup_{\substack{A\subseteq \mathbb{R}^{\ell_n} \\ A\ \text{rectangular}}} \left| \mathbb{P}\left[ \mathbf{z}_1 \in A \right]
-
\mathbb{P}\left[ \mathbf{z}_2 \in A \right] \right| \leq C\sqrt{\Vert \boldsymbol{\Omega}_{1} - \boldsymbol{\Omega}_{2} \Vert_{\infty}}\log \ell_n,
\end{align*}
where $\Vert \cdot\Vert_\infty$ denotes the supremum norm, and $C$ is an absolute constant.
\qed
\end{lem}
\begin{lem}[Equation (3.5) in \citealt*{Gine-Latala-Zinn_2000_Ustat}]\label{lem:ustatistic concentration inequality}
For a degenerate and decoupled second order U-statistic, $\sum_{i,j=1,i\neq j}^n h_{ij}(x_i,\tilde{x}_j)$, the following holds:
\begin{align*}
\mathbb{P}\left[ \left|\sum_{i,j,i\neq j}^n u_{ij}(x_i,\tilde{x}_j)\right| > t \right] \leq C\exp\left\{ -\frac{1}{C}\min\left[ \frac{t}{D},\ \left(\frac{t}{B}\right)^{\frac{2}{3}},\ \left( \frac{t}{A} \right)^{\frac{1}{2}} \right] \right\},
\end{align*}
where $C$ is some universal constant, and $A$, $B$ and $D$ are any constants satisfying
\begin{align*}
A &\geq \max_{1\leq i,j\leq n}\sup_{u,v}| u_{ij}(u,v) |\\
B^2 &\geq \max_{1\leq i,j\leq n}\left[ \sup_{v}\left|\sum_{i=1}^n \mathbb{E} u_{ij}(x_i,v)^2\right| ,\ \sup_{u}\left|\sum_{j=1}^n \mathbb{E} u_{ij}(u,\tilde{x}_j)^2\right| \right]\\
D^2 &\geq \sum_{i,j=1,i\neq j}^n \mathbb{E} u_{ij}(x_i,\tilde{x}_j)^2.
\end{align*}
where $\{{x}_i, 1\leq i\leq n\}$ are independent random variables, and $\{\tilde{x}_i, 1\leq i\leq n\}$ is an independent copy of $\{{x}_i, 1\leq i\leq n\}$. \qed
\end{lem}
\begin{remark}
To apply the above lemma, an additional decoupling step is usually needed. Fortunately, the decoupling step only introduces an extra constant, but will not affect the order of the tail probability bound. Formally,
\begin{lem}[\citealt{DeLaPena-MontgomerySmith_1995_AoP}]
Consider the setting of Lemma \ref{lem:ustatistic concentration inequality}. Then
\begin{align*}
\mathbb{P}\left[ \left|\sum_{i,j,i\neq j}^n u_{ij}(x_i,{x}_j)\right| > t \right] \leq C\cdot \mathbb{P}\left[ C\left|\sum_{i,j,i\neq j}^n u_{ij}(x_i,\tilde{x}_j)\right| > t \right],
\end{align*}
where $C$ is a universal constant.\qed
\end{lem}
As a result, we will apply Lemma \ref{lem:ustatistic concentration inequality} without explicitly mentioning the decoupling step or the extra constant it introduces.\qed
\end{remark}
\subsection{Proof of Theorem \ref{thm:strong approximation}}
To bound the distance between the two processes, $\tilde{\mathfrak{T}}_G(\cdot)$ and $\mathfrak{B}_G(\cdot)$, we employ the proof strategy of \cite*{Gine-Koltchinskii-Sakhanenko_2004_PTRF}. Recall that $F$ denotes the distribution of $x_i$, and we define
\begin{align*}
\mathcal{K}_{h,\mathsf{x}}\circ F^{-1}(x) = \mathcal{K}_{h,\mathsf{x}}(F^{-1}(x)).
\end{align*}
Take $v<v'$ in $[0,1]$, we have
\begin{align*}
&\ \left|\mathcal{K}_{h,\mathsf{x}}\circ F^{-1}(v) - \mathcal{K}_{h,\mathsf{x}}\circ F^{-1}(v')\right|\\
&= \left|\frac{ \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{x}}^{-1}R(u)\Big[\mathds{1}(F^{-1}(v)\leq \mathsf{x} + hu) - \mathds{1}(F^{-1}(v')\leq \mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\right|\\
&\leq \frac{ \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} \left|c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{x}}^{-1}R(u)\right|\Big[\mathds{1}(F^{-1}(v)\leq \mathsf{x} + hu) - \mathds{1}(F^{-1}(v')\leq \mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}.
\end{align*}
Therefore, the function $\mathcal{K}_{h,\mathsf{x}}\circ F^{-1}(\cdot)$ has a total variation bounded by
\begin{align*}
&\ \frac{ \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} \left|c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{x}}^{-1}R(u)\right|\Big[\mathds{1}(F^{-1}(0)\leq \mathsf{x} + hu) - \mathds{1}(F^{-1}(1)\leq \mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\\
&= \frac{ \int_{-1}^1 \left|c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{x}}^{-1}R(u)\right| K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\leq C_4\frac{1}{\sqrt{h}}.
\end{align*}
It is well-known that functions of bounded variation can be approximated (pointwise) by convex combination of indicator functions of half intervals. To be more precise,
\begin{align*}
\Big\{\mathcal{K}_{h,\mathsf{x}}\circ F^{-1}(\cdot):\ \mathsf{x}\in\mathcal{I}\Big\}\subset C_4\frac{1}{\sqrt{h}}\overline{\mathrm{conv}}\Big\{ \pm\mathds{1}(\cdot\leq t), \pm\mathds{1}(\cdot\geq t)\Big\}.
\end{align*}
Following (2.3) and (2.4) of \cite*{Gine-Koltchinskii-Sakhanenko_2004_PTRF}, we have
\begin{align*}
\mathbb{P}\left[ \sup_{\mathsf{x}\in\mathcal{I}}\left| \tilde{\mathfrak{T}}_G(\mathsf{x}) - \mathfrak{B}_G(\mathsf{x}) \right| > \frac{C_4(u+C_5\log n)}{\sqrt{nh}} \right] \leq C_5e^{-C_5u},
\end{align*}
where $C_5$ is some universal constant.
\subsection{Proof of Lemma \ref{lem:VC-type}}
Take $|\mathsf{x}-\mathsf{y}|\leq \varepsilon$ to be some small number, then
\begin{align*}
&\ \mathcal{K}_{h,\mathsf{x}}(x) - \mathcal{K}_{h,\mathsf{y}}(x)
= \frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{x}}^{-1} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\\
&\qquad\qquad\qquad\qquad\qquad\qquad - \frac{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{y}}^{-1} \int_{\frac{\mathcal{X}-\mathsf{y}}{h}} R(u)\Big[\mathds{1}(x\leq \mathsf{y} + hu) - F(\mathsf{y} + hu)\Big] K\left(u\right) g(\mathsf{y}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}}\\
=&\ \left(\frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{x}}^{-1}}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}} - \frac{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{y}}^{-1}}{\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}}\right)\left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u\right)\\
+&\left( \frac{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{y}}^{-1}}{\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}}\right)
\left(\frac{1}{h}\int_{\mathcal{X}} \Big[R\left(\frac{u-\mathsf{x}}{h}\right)K\left(\frac{u-\mathsf{x}}{h}\right) - R\left(\frac{u-\mathsf{y}}{h}\right)K\left(\frac{u-\mathsf{y}}{h}\right)\Big]\Big[\mathds{1}(x\leq u) - F(u)\Big] g(u)\mathrm{d} u\right)\\
\tag{I}=&\ \frac{1}{\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}}\left(c_{h,\mathsf{x}}^\prime\Upsilon_{h} - c_{h,\mathsf{y}}^\prime\Upsilon_{h}\right)\Gamma_{h,\mathsf{x}}^{-1}\left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u\right)\\
\tag{II}+&\ \frac{1}{\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}}c_{h,\mathsf{y}}^\prime\Upsilon_{h}\left(\Gamma_{h,\mathsf{x}}^{-1} - \Gamma_{h,\mathsf{y}}^{-1}\right)\left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u\right)\\
\tag{III}+& \left(\frac{1}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}} - \frac{1}{\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}}\right)
c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{x}}^{-1}\left(\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u\right)\\
\tag{IV}+&\left( \frac{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Gamma_{h,\mathsf{y}}^{-1}}{\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}}\right)
\left(\frac{1}{h}\int_{\mathcal{X}} \Big[R\left(\frac{u-\mathsf{x}}{h}\right)K\left(\frac{u-\mathsf{x}}{h}\right) - R\left(\frac{u-\mathsf{y}}{h}\right)K\left(\frac{u-\mathsf{y}}{h}\right)\Big]\Big[\mathds{1}(x\leq u) - F(u)\Big] g(u)\mathrm{d} u\right).
\end{align*}
For term (I), its variance (replace the placeholder $x$ by $x_i$) is
\begin{align*}
\mathbb{V}[\text{(I)}] &= \frac{1}{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}\left(c_{h,\mathsf{x}}^\prime\Upsilon_{h} - c_{h,\mathsf{y}}^\prime\Upsilon_{h}\right)\Omega_{h,\mathsf{x}}\left(c_{h,\mathsf{x}}^\prime\Upsilon_{h} - c_{h,\mathsf{y}}^\prime\Upsilon_{h}\right)^\prime = O\left(\frac{1}{h}r_1(\varepsilon,h)^2r_2(h)^2\right).
\end{align*}
Term (II) has variance
\begin{align*}
\mathbb{V}[\text{(II)}] &= \frac{1}{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}c_{h,\mathsf{y}}^\prime\Upsilon_{h}\left(\Gamma_{h,\mathsf{x}}^{-1} - \Gamma_{h,\mathsf{y}}^{-1}\right)\Sigma_{h,\mathsf{x}}\left(\Gamma_{h,\mathsf{x}}^{-1} - \Gamma_{h,\mathsf{y}}^{-1}\right)^\prime \left(c_{h,\mathsf{y}}^\prime\Upsilon_{h}\right)^\prime
= O \left(\frac{1}{h}\left(\frac{\varepsilon}{h}\wedge 1\right)^2\right),
\end{align*}
where the order $\frac{\varepsilon}{h}\wedge 1$ comes from the difference $\Gamma_{h,\mathsf{x}}^{-1} - \Gamma_{h,\mathsf{y}}^{-1}$.
Next for term (III), we have
\begin{align*}
\mathbb{V}[\text{(III)}] &= \left(\frac{1}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}} - \frac{1}{\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}}\right)^2c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}\\
&= \left(1 - \sqrt{1+\frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}-c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}}\right)^2\\
&\asymp \left(\frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}-c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}\right)^2\\
&= \left(\frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}(\Omega_{h,\mathsf{x}}-\Omega_{h,\mathsf{y}})\Upsilon_{h}c_{h,\mathsf{x}}}{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}
+\frac{(c_{h,\mathsf{x}}^\prime\Upsilon_{h}-c_{h,\mathsf{y}}^\prime\Upsilon_{h})\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{x}}+(c_{h,\mathsf{x}}^\prime\Upsilon_{h}-c_{h,\mathsf{y}}^\prime\Upsilon_{h})\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}\right)^2.
\end{align*}
The first term has bound
\begin{align*}
\frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}(\Omega_{h,\mathsf{x}}-\Omega_{h,\mathsf{y}})\Upsilon_{h}c_{h,\mathsf{x}}}{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}} = O \left(\frac{\varepsilon}{h}\right).
\end{align*}
The third term has bound
\begin{align*}
\frac{(c_{h,\mathsf{x}}^\prime\Upsilon_{h}-c_{h,\mathsf{y}}^\prime\Upsilon_{h})\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}
&\precsim \frac{|(c_{h,\mathsf{x}}^\prime\Upsilon_{h}-c_{h,\mathsf{y}}^\prime\Upsilon_{h})\Omega_{h,\mathsf{y}}^{1/2}|}{\sqrt{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}}
= O \left(\frac{1}{\sqrt{h}}r_1(\varepsilon,h)r_2(h)\right).
\end{align*}
Finally, the second term can be bounded as
\begin{align*}
\frac{(c_{h,\mathsf{x}}^\prime\Upsilon_{h}-c_{h,\mathsf{y}}^\prime\Upsilon_{h})\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{x}}}{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}} &= \frac{(c_{h,\mathsf{x}}^\prime\Upsilon_{h}-c_{h,\mathsf{y}}^\prime\Upsilon_{h})\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}+(c_{h,\mathsf{x}}^\prime\Upsilon_{h}-c_{h,\mathsf{y}}^\prime\Upsilon_{h})\Omega_{h,\mathsf{y}}(c_{h,\mathsf{x}}^\prime\Upsilon_{h}-c_{h,\mathsf{y}}^\prime\Upsilon_{h})^\prime}{c_{h,\mathsf{y}}^\prime\Upsilon_{h}\Omega_{h,\mathsf{y}}\Upsilon_{h}c_{h,\mathsf{y}}}\\
&= O \left(\frac{1}{\sqrt{h}}r_1(\varepsilon,h)r_2(h) + \frac{1}{h}r_1(\varepsilon,h)^2r_2(h)^2\right).
\end{align*}
Overall, we have that
\begin{align*}
\mathbb{V}[\text{(III)}] = O \left(\frac{\varepsilon^2}{h^2} + \frac{1}{h}r_1(\varepsilon,h)^2r_2(h)^2 + \frac{1}{h^2}r_1(\varepsilon,h)^4r_2(h)^4\right).
\end{align*}
Given our assumptions on the basis function and on the kernel function, it is obvious that term (IV) has variance
\begin{align*}
\mathbb{V}[\text{(IV)}] &= O \left(\frac{1}{h}\left(\frac{\varepsilon}{h}\wedge 1\right)^2\right).
\end{align*}
The bound on $\mathbb{E}[\sup_{\mathsf{x}\in\mathcal{I}}|\mathfrak{B}_G(\mathsf{x})|]$ can be found by standard entropy calculation, and the bound on $\mathbb{E}[\sup_{\mathsf{x}\in\mathcal{I}}|\mathfrak{T}_G(\mathsf{x})|]$ is obtained by the following fact
\begin{align*}
\mathbb{E}\left[\sup_{\mathsf{x}\in\mathcal{I}}|\mathfrak{T}_G(\mathsf{x})|\right] \leq \mathbb{E}\left[\sup_{\mathsf{x}\in\mathcal{I}}|\mathfrak{B}_G(\mathsf{x})|\right] + \mathbb{E}\left[\sup_{\mathsf{x}\in\mathcal{I}}|\tilde{\mathfrak{T}}_G(\mathsf{x})-B_G(\mathsf{x})|\right],
\end{align*}
and that
\begin{align*}
\mathbb{E}\left[\sup_{\mathsf{x}\in\mathcal{I}}|\tilde{\mathfrak{T}}_G(\mathsf{x})-\mathfrak{B}_G(\mathsf{x})|\right] &= \int_0^\infty \mathbb{P}\left[\sup_{\mathsf{x}\in\mathcal{I}}|\tilde{\mathfrak{T}}_G(\mathsf{x})-\mathfrak{B}_G(\mathsf{x})|>u\right]\mathrm{d} u = O\left(\frac{\log n}{\sqrt{nh}}\right) = o(\sqrt{\log n}),
\end{align*}
which follows from Theorem \ref{thm:strong approximation} and our assumption that $\log n/(nh)\to 0$.
\subsection{Proof of Lemma \ref{lemma:standard error, local projection uniform}}
We adopt the following decomposition (the integration is always on $\frac{\mathcal{X}-\mathsf{y}}{h}\times\frac{\mathcal{X}-\mathsf{x}}{h}$, unless otherwise specified):
\begin{align*}
\tag{I}&\frac{1}{n}\sum_{i=1}^n\iint R(u)R(v)^\prime\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - {F}(\mathsf{x} + hu)\Big]\Big[\mathds{1}(x_i\leq \mathsf{y} + hv) - {F}(\mathsf{y} + hv)\Big] K(u)K(v) g(\mathsf{x}+hu)g(\mathsf{y}+hv)\mathrm{d} u\mathrm{d} v\\
\tag{II}&- \iint R(u)R(v)^\prime\Big[\hat{F}(\mathsf{x} + hu) - {F}(\mathsf{x} + hu)\Big]\Big[\hat{F}(\mathsf{y} + hv) - {F}(\mathsf{y} + hv)\Big] K(u)K(v) g(\mathsf{x}+hu)g(\mathsf{y}+hv)\mathrm{d} u\mathrm{d} v.
\end{align*}
By the uniform convergence of the empirical distribution function, we have that
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}}|\text{(II)}| =O_{\mathbb{P}}\left(\frac{1}{n}\right).
\end{align*}
From the definition of $\Sigma_{h,\mathsf{x},\mathsf{y}}$, we know that
\begin{align*}
\mathbb{E}[\text{(I)}] &= \Sigma_{h,\mathsf{x},\mathsf{y}}.
\end{align*}
As (I) is a sum of bounded terms, we can apply Lemma \ref{lem:uniform convergence of iid sum} and easily show that
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}}\left|\text{(I)} - \Sigma_{h,\mathsf{x},\mathsf{y}}\right| + O_{\mathbb{P}}\left( \sqrt{\frac{\log n}{n}}\right).
\end{align*}
\subsection{Proof of Lemma \ref{lem:local projection uniform term 1}}
We rewrite \eqref{eq:local projection uniform term 1} as
\begin{align*}
|\text{\eqref{eq:local projection uniform term 1}}|
&= \left|\sqrt{n}\sqrt{\frac{{c_{h,\mathsf{x}}^\prime\Upsilon_{h}{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}{{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}}\frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}{\Gamma}_h^{-1}\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[F(\mathsf{x} + hu) - \theta^\prime R(u)\Upsilon_h^{-1}\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\right|\\
&\leq \sqrt{\frac{n}{h}}\left[\sup_{\mathsf{x}\in\mathcal{I}}\sqrt{\frac{{c_{h,\mathsf{x}}^\prime\Upsilon_{h}{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}{{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}}\right]\left[\sup_{\mathsf{x}\in\mathcal{I}} \left|\int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[F(\mathsf{x} + hu) - \theta^\prime R(u)\Upsilon_h^{-1}\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u\right| \right]\\
&= O_{\mathbb{P}}\left( \sqrt{\frac{n}{h}}\sup_{\mathsf{x}\in\mathcal{I}}\varrho(h,\mathsf{x}) \right),
\end{align*}
where the final bound holds uniformly for $\mathsf{x}\in\mathcal{I}$.
Next, we expand term \eqref{eq:local projection uniform term 2} as
\begin{align*}
\text{\eqref{eq:local projection uniform term 2}}
&= \frac{1}{\sqrt{n}}\sum_{i=1}^n\frac{c^\prime_{h,\mathsf{x}}\Upsilon_{h}{\Gamma}_{h,\mathsf{x}}^{-1} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\\
&\qquad+ \frac{1}{\sqrt{n}}\sum_{i=1}^n\left[ 1- \sqrt{\frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}} \right]\frac{c^\prime_{h,\mathsf{x}}\Upsilon_{h}{\Gamma}_{h,\mathsf{x}}^{-1} \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) g(\mathsf{x}+hu)\mathrm{d} u}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\\
&= \mathfrak{T}_G(\mathsf{x}) + \underbrace{\left[ 1- \sqrt{\frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\ \right]\mathfrak{T}_G(\mathsf{x})}_{\textstyle \text{(I)}}.
\end{align*}
Term (I) can be easily bounded by
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}|\text{(I)}| =O_{\mathbb{P}}\left( \left(\sqrt{\frac{\log n}{nh^2}}\right)\mathbb{E}\left[\sup_{\mathsf{x}\in\mathcal{I}}|\mathfrak{T}_G(\mathsf{x})|\right] \right) = O_{\mathbb{P}}\left( \frac{\log n}{\sqrt{nh^2}} \right).
\end{align*}
\subsection{Proof of Theorem \ref{thm:strong approximation, local projection}}
The claim follows from Theorem \ref{thm:strong approximation} and previous lemmas.
\subsection{Proof of Theorem \ref{thm:feasible uniform approximation local projection}}
Let $\mathcal{I}_\varepsilon$ be an $\varepsilon$-covering (with respect to the Euclidean metric) of $\mathcal{I}$, and assume $\varepsilon\leq h$. Then the process ${\mathfrak{B}}_G(\cdot)$ can be decomposed into:
\begin{align*}
{\mathfrak{B}}_G(\mathsf{x}) &= {\mathfrak{B}}_G(\Pi_{\mathcal{I}_\varepsilon}(\mathsf{x})) + {\mathfrak{B}}_G(\mathsf{x}) - {\mathfrak{B}}_G(\Pi_{\mathcal{I}_\varepsilon}(\mathsf{x})),
\end{align*}
where $\Pi_{\mathcal{I}_\varepsilon}:\mathcal{I}\to \mathcal{I}_\varepsilon$ is a mapping satisfying:
\begin{align*}
\Pi_{\mathcal{I}_\varepsilon}(\mathsf{x}) &= \operatorname*{argmin}_{\mathsf{y}\in \mathcal{I}_\varepsilon}|\mathsf{y}-\mathsf{x}|.
\end{align*}
We first study the properties of ${\mathfrak{B}}_G(\mathsf{x}) - {\mathfrak{B}}_G(\Pi_{\mathcal{I}_\varepsilon}(\mathsf{x}))$. With standard entropy calculation, one has:
\begin{align*}
\mathbb{E}\left[\sup_{\mathsf{x}\in\mathcal{I}} |{\mathfrak{B}}_G(\mathsf{x}) - {\mathfrak{B}}_G(\Pi_{\mathcal{I}_\varepsilon}(\mathsf{x}))|\right] &\leq \mathbb{E}\left[ \sup_{\mathsf{x},\mathsf{y}\in \mathcal{I},|\mathsf{x}-\mathsf{y}|\leq \varepsilon} |{\mathfrak{B}}_G(\mathsf{x})-{\mathfrak{B}}_G(\mathsf{y})| \right]
\leq \mathbb{E}\left[ \sup_{\mathsf{x},\mathsf{y}\in \mathcal{I},\sigma(\mathsf{x},\mathsf{y})\leq \delta(\varepsilon)} |{\mathfrak{B}}_G(\mathsf{x})-{\mathfrak{B}}_G(\mathsf{y})| \right]\\
&\precsim \int_0^{\delta(\varepsilon)} \sqrt{\log N(\lambda, \mathcal{I}, \sigma_G)}\mathrm{d} \lambda,
\end{align*}
where
\begin{align*}
\delta(\varepsilon) &= C\left(\frac{1}{\sqrt{h}}\frac{\varepsilon}{h} + \frac{1}{\sqrt{h}}r_1(\varepsilon,h)r_2(h) + \frac{1}{h}r_1(\varepsilon,h)^2r_2(h)^2\right),
\end{align*}
for some $C>0$ that does not depend on $\varepsilon$ and $h$, and $N(\lambda, \mathcal{I}, \sigma_G)$ is the covering number of $\mathcal{I}$ measured by the pseudo metric $\sigma_G(\cdot,\cdot)$, which satisfies
\begin{align*}
N(\lambda, \mathcal{I}, \sigma_G) \precsim \frac{1}{\delta^{-1}(\lambda)}.
\end{align*}
Therefore, we have
\begin{align*}
\tag{I}\mathbb{E}\left[\sup_{\mathsf{x}\in\mathcal{I}} |{\mathfrak{B}}_G(\mathsf{x}) - {\mathfrak{B}}_G(\Pi_{\mathcal{I}_\varepsilon}(\mathsf{x}))|\right] &\precsim \left(\frac{1}{\sqrt{h}}\frac{\varepsilon}{h} + \frac{1}{\sqrt{h}}r_1(\varepsilon,h)r_2(h) + \frac{1}{h}r_1(\varepsilon,h)^2r_2(h)^2\right)\sqrt{\log n}.
\end{align*}
A similar bound holds for the process $\hat{\mathfrak{B}}_G(\cdot)$ due to the uniform consistency of the covariance estimator.
Now consider the discretized version of ${\mathfrak{B}}_G(\cdot)$ and $\hat{\mathfrak{B}}_G(\cdot)$. By applying Lemmas \ref{lemma:standard error, local projection uniform} and \ref{lem:gaussian comparison}, we directly obtain the following bound:
\begin{align*}
\tag{II}\sup_{ A\ \text{rectangular}} \left| \mathbb{P}\Big[ \Big\{\mathfrak{B}_G(\Pi_{\mathcal{I}_\varepsilon}(\mathsf{x})),\mathsf{x}\in\mathcal{I}\Big\} \in A \Big]
-
\mathbb{P}^\star\Big[ \Big\{\hat{\mathfrak{B}}_G(\Pi_{\mathcal{I}_\varepsilon}(\mathsf{x})),\mathsf{x}\in\mathcal{I}\Big\} \in A \Big] \right| = O_{\mathbb{P}}\left( \left(\frac{\log n}{nh^2}\right)^{\frac{1}{4}}\log \frac{1}{\varepsilon} \right).
\end{align*}
As $\varepsilon$ appears in (I) polynomially but only logarithmically in (II), it is possible to choose $\varepsilon$ sufficiently small so that the discretization error becomes negligible. Therefore,
\begin{align*}
\sup_{u\in \mathbb{R}}\left| \mathbb{P}\Big[ \sup_{\mathsf{x}\in\mathcal{I}}|\mathfrak{B}_G(\mathsf{x})| \leq u \Big]
-
\mathbb{P}^\star\Big[\sup_{\mathsf{x}\in\mathcal{I}}|\hat{\mathfrak{B}}_G(\mathsf{x})| \leq u \Big] \right| = O_{\mathbb{P}}\left(\frac{\log^{\frac{5}{4}} n}{(nh^2)^{\frac{1}{4}}} \right).
\end{align*}
\subsection{Proof of Lemma \ref{lem:uniform consistency of Gamma}}
We apply Lemma \ref{lem:uniform convergence of iid sum}. For simplicity, assume $R(\cdot)$ is scalar, and let
\begin{align*}
u_{i,h}(\mathsf{x}) = R\left( \frac{x_i-\mathsf{x}}{h} \right)^2\frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right) - \Gamma_{h,\mathsf{x}}.
\end{align*}
Then it is easy to see that
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}\max_{1\leq i\leq n}\mathbb{V}[u_{i,h}(\mathsf{x})] = O(h^{-1}), \qquad \sup_{\mathsf{x}\in\mathcal{I}}\max_{1\leq i\leq n}|u_{i,h}(\mathsf{x})| = O(h^{-1}).
\end{align*}
Let $|\mathsf{x} - \mathsf{y}|\leq \varepsilon\leq h$, we also have
\begin{align*}
\left|u_{i,h}(\mathsf{x}) - u_{i,h}(\mathsf{y})\right| &\leq \left|R\left( \frac{x_i-\mathsf{x}}{h} \right)^2\frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right) - R\left( \frac{x_i-\mathsf{y}}{h} \right)^2\frac{1}{h}K\left(\frac{x_i-\mathsf{y}}{h}\right)\right| + \left|\Gamma_{h,\mathsf{x}} - \Gamma_{h,\mathsf{y}} \right|\\
&\leq \left|R\left( \frac{x_i-\mathsf{x}}{h} \right)^2-R\left( \frac{x_i-\mathsf{y}}{h} \right)^2\right|\frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right) +
R\left( \frac{x_i-\mathsf{y}}{h} \right)^2\frac{1}{h}\left|K\left(\frac{x_i-\mathsf{x}}{h}\right) - K\left(\frac{x_i-\mathsf{y}}{h}\right)\right| + \left|\Gamma_{h,\mathsf{x}} - \Gamma_{h,\mathsf{y}} \right|\\
&\leq M\left[\frac{\varepsilon}{h}\frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right) + \frac{\varepsilon}{h}\frac{1}{h}K^{\dag}\left(\frac{x_i-\mathsf{x}}{h}\right) + \frac{1}{h}K^{\ddag}\left(\frac{x_i-\mathsf{x}}{h}\right) + \frac{\varepsilon}{h}\right].
\end{align*}
where $M$ is some constant that does not depend on $n$, $h$ or $\varepsilon$. Then it is easy to see that
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}\max_{1\leq i\leq n}\mathbb{V}[u_{i,h,\varepsilon}(\mathsf{x})] = O\left(\frac{\varepsilon}{h^2}\right), \qquad \sup_{\mathsf{x}\in\mathcal{I}}\max_{1\leq i\leq n}|u_{i,h,\varepsilon}(\mathsf{x}) - \mathbb{E}[u_{i,h,\varepsilon}(\mathsf{x})]| = O(h^{-1}),\quad \sup_{\mathsf{x}\in\mathcal{I}}\max_{1\leq i\leq n}\mathbb{E}[|u_{i,h,\varepsilon}(\mathsf{x})|] = O\left(\frac{\varepsilon}{h}\right).
\end{align*}
Now take $\varepsilon = \sqrt{h\log n/n}$, then $\log N(\varepsilon, \mathcal{I}, |\cdot|)=O(\log n)$. Lemma \ref{lem:uniform convergence of iid sum} implies that
\begin{align*}
\sup_{\mathsf{x}\in \mathcal{I}}\left|\frac{1}{n}\sum_{i=1}^n R\left( \frac{x_i-\mathsf{x}}{h} \right)^2\frac{1}{h}K\left(\frac{x_i-\mathsf{x}}{h}\right) - \Gamma_{h,\mathsf{x}}\right| = O_{\mathbb{P}}\left( \sqrt{\frac{\log n}{nh}} \right).
\end{align*}
\subsection{Proof of Lemma \ref{lemma:standard error, local regression uniform}}
Let $R_i(\mathsf{x}) = R(x_i-\mathsf{x})$ and $W_i(\mathsf{x}) = K((x_i-\mathsf{x})/h)/h$, then we split $\hat{\Sigma}_{h,\mathsf{x},\mathsf{y}}$ into two terms,
\begin{align*}
\text{(I)} &= \frac{1}{n^3}\sum_{i,j,k} \Upsilon_{h}R_j(\mathsf{x})R_k(\mathsf{y})^\prime \Upsilon_{h} W_j(\mathsf{x})W_k(\mathsf{y}) \Big(\mathds{1}(x_i\leq x_j) - F(x_j)\Big) \Big(\mathds{1}(x_i\leq x_k) - F(x_k)\Big)\\
\text{(II)} &= -\frac{1}{n^2}\sum_{j,k} \Upsilon_{h}R_j(\mathsf{x})R_k(\mathsf{y})^\prime \Upsilon_{h} W_j(\mathsf{x})W_k(\mathsf{y}) \Big(\hat{F}(x_j)-{F}(x_j)\Big) \Big(\hat{F}(x_k) - {F}(x_k)\Big).
\end{align*}
(II) satisfies
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}}\left|\text{(II)}\right| \leq \sup_{x}|\hat{F}(x) - F(x)|^2\left(\sup_{\mathsf{x}\in \mathcal{I}}\frac{1}{n}\sum_{j} \left| \Upsilon_{h}R_j(\mathsf{x})W_j(\mathsf{x})\right|\right)^2.
\end{align*}
It is obvious that
\begin{align*}
\sup_{x}|\hat{F}(x) - F(x)|^2 = O_{\mathbb{P}}\left( \frac{1}{n}\right).
\end{align*}
As for the second part, one can employ the same technique used to prove Lemma \ref{lem:uniform consistency of Gamma} and show that
\begin{align*}
\sup_{\mathsf{x}\in \mathcal{I}}\frac{1}{n}\sum_{j} \left| \Upsilon_{h}R_j(\mathsf{x})W_j(\mathsf{x})\right| = O_{\mathbb{P}}(1),
\end{align*}
implying that
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}}\left|\text{(II)}\right| = O_{\mathbb{P}}\left(\frac{1}{n}\right).
\end{align*}
For (I), we first define
\begin{align*}
u_{ij}(\mathsf{x}) &= \Upsilon_{h}R_j(\mathsf{x}) W_j(\mathsf{x}) \Big(\mathds{1}(x_i\leq x_j) - F(x_j)\Big),
\end{align*}
and
\begin{align*}
\bar{u}_i(\mathsf{x}) &= \mathbb{E}[u_{ij}(\mathsf{x})|x_i;i\neq j],\qquad \hat{u}_i(\mathsf{x}) = \frac{1}{n}\sum_j u_{ij}(\mathsf{x}).
\end{align*}
Then
\begin{align*}
\text{(I)} &= \frac{1}{n}\sum_i\left( \frac{1}{n}\sum_j u_{ij}(\mathsf{x}) \right)\left( \frac{1}{n}\sum_j u_{ij}(\mathsf{y}) \right)^\prime = \frac{1}{n}\sum_i \hat{u}_i(\mathsf{x})\hat{u}_i(\mathsf{y})^\prime\\
&= \frac{1}{n}\sum_i \bar{u}_i(\mathsf{x})\bar{u}_i(\mathsf{y})^\prime + \frac{1}{n}\sum_i\left(\hat{u}_i(\mathsf{x}) - \bar{u}_i(\mathsf{x})\right)\hat{u}_i(\mathsf{y})^\prime + \frac{1}{n}\sum_i\bar{u}_i(\mathsf{x})\left(\hat{u}_i(\mathsf{y}) - \bar{u}_i(\mathsf{y})\right)^\prime\\
&= \underbrace{\frac{1}{n}\sum_i \bar{u}_i(\mathsf{x})\bar{u}_i(\mathsf{y})^\prime}_{\text{(I.1)}} + \underbrace{\frac{1}{n}\sum_i\left(\hat{u}_i(\mathsf{x}) - \bar{u}_i(\mathsf{x})\right)\bar{u}_i(\mathsf{y})^\prime}_{\text{(I.2)}} + \underbrace{\frac{1}{n}\sum_i\bar{u}_i(\mathsf{x})\left(\hat{u}_i(\mathsf{y}) - \bar{u}_i(\mathsf{y})\right)^\prime}_{\text{(I.3)}}\\
&\qquad\qquad\qquad + \underbrace{\frac{1}{n}\sum_i\left(\hat{u}_i(\mathsf{x}) - \bar{u}_i(\mathsf{x})\right)\left(\hat{u}_i(\mathsf{y}) - \bar{u}_i(\mathsf{y})\right)^\prime}_{\text{(I.4)}}.
\end{align*}
Term (I.1) has been analyzed in Lemma \ref{lemma:standard error, local projection uniform}, which satisfies
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}}\left|\text{(I.1)} - \Sigma_{h,\mathsf{x},\mathsf{y}}\right| = O_{\mathbb{P}}\left( \sqrt{\frac{\log n}{n}} \right).
\end{align*}
Term (I.2) has expansion:
\begin{align*}
\text{(I.2)} &= \frac{1}{n^2}\sum_{i,j}\left(u_{ij}(\mathsf{x}) - \bar{u}_i(\mathsf{x})\right)\bar{u}_i(\mathsf{y})^\prime
= \underbrace{\frac{1}{n^2}\sum_{\substack{i,j\\ \text{distinct}}}\left(u_{ij}(\mathsf{x}) - \bar{u}_i(\mathsf{x})\right)\bar{u}_i(\mathsf{y})^\prime}_{\text{(I.2.1)}}
+ \underbrace{\frac{1}{n^2}\sum_{i}\left(u_{ii}(\mathsf{x}) - \bar{u}_i(\mathsf{x})\right)\bar{u}_i(\mathsf{y})^\prime}_{\text{(I.2.2)}}.
\end{align*}
By the same technique of Lemma \ref{lem:uniform consistency of Gamma}, one can show that
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}}|\text{(I.2.2)}| = O_{\mathbb{P}}\left(\frac{1}{n}\right).
\end{align*}
We need a further decomposition to make (I.2.1) a degenerate U-statistic:
\begin{align*}
\text{(I.2.1)} &= \underbrace{\frac{n-1}{n^2}\sum_{j}\mathbb{E}\left[\left.\left(u_{ij}(\mathsf{x}) - \bar{u}_i(\mathsf{x})\right)\bar{u}_i(\mathsf{y})^\prime\right|x_j\right]}_{\text{(I.2.1.1)}}\\
&\qquad + \underbrace{\frac{1}{n^2}\sum_{\substack{i,j\\ \text{distinct}}}\Bigg\{\left(u_{ij}(\mathsf{x}) - \bar{u}_i(\mathsf{x})\right)\bar{u}_i(\mathsf{y})^\prime - \mathbb{E}\left[\left.\left(u_{ij}(\mathsf{x}) - \bar{u}_i(\mathsf{x})\right)\bar{u}_i(\mathsf{y})^\prime\right|x_j\right]\Bigg\}}_{\text{(I.2.1.2)}}.
\end{align*}
(I.2.1) has zero mean. By discretizing $\mathcal{I}$ and apply Bernstein's inequality, one can show that the (I.2.1.1) has order $O_{\mathbb{P}}\left(\sqrt{\log n/n}\right)$.
For (I.2.1.2), we first discretize $\mathcal{I}$ and then apply a Bernstein-type inequality (Lemma \ref{lem:ustatistic concentration inequality}) for degenerate U-statistics, which gives an order
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}}\left|\text{(I.2.1.2)} \right| = O_{\mathbb{P}}\left( \frac{\log n}{\sqrt{n^2h}} \right).
\end{align*}
Overall, we have
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}}|\text{(I.2)}| = O_{\mathbb{P}}\left( \frac{1}{n} + \sqrt{\frac{\log n}{n}} + \frac{\log n}{\sqrt{n^2h}} \right) = O_{\mathbb{P}}\left( \sqrt{\frac{\log n}{n}} \right),
\end{align*}
and the same bound applies to (I.3).
For (I.4), one can show that
\begin{align*}
\sup_{\mathsf{x}\in \mathcal{I}}\sup_{x\in\mathcal{X}} \left| \frac{1}{n}\sum_j \Upsilon_{h}R_j(\mathsf{x}) W_j(\mathsf{x}) \Big(\mathds{1}(x\leq x_j) - F(x_j)\Big) - \mathbb{E}\left[\Upsilon_{h}R_j(\mathsf{x}) W_j(\mathsf{x}) \Big(\mathds{1}(x\leq x_j) - F(x_j)\Big)\right] \right| = O_{\mathbb{P}}\left( \sqrt{\frac{\log n}{nh}} \right),
\end{align*}
which means
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}} |\text{(I.4)}| = O_{\mathbb{P}}\left( \frac{\log n}{nh} \right) = O_{\mathbb{P}}\left( \sqrt{\frac{\log n}{n}} \right),
\end{align*}
under our assumption that $\log n/(nh^2)\to 0$.
As a result, we have
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}} \left| \hat\Sigma_{h,\mathsf{x}, \mathsf{y}} - \Sigma_{h,\mathsf{x}, \mathsf{y}} \right| = O_{\mathbb{P}}\left( \sqrt{\frac{\log n}{n}} \right).
\end{align*}
Now take $c$ to be a generic vector. Then we have
\begin{align*}
\frac{c^\prime_{h,\mathsf{x}}\Upsilon_h(\hat{\Omega}_{h,\mathsf{x},\mathsf{y}}-{\Omega}_{h,\mathsf{x},\mathsf{y}})\Upsilon_h c_{h,\mathsf{y}} }{\sqrt{c^\prime_{h,\mathsf{x}}\Upsilon_h{\Omega}_{h,\mathsf{x}}\Upsilon_h c_{h,\mathsf{x}}}\sqrt{c^\prime_{h,\mathsf{y}}\Upsilon_h{\Omega}_{h,\mathsf{y}}\Upsilon_h c_{h,\mathsf{y}}}}
&= \frac{c^\prime_{h,\mathsf{x}}\Upsilon_h \hat\Gamma_{h,\mathsf{x}}^{-1} (\hat\Sigma_{h,\mathsf{x},\mathsf{y}} - \Sigma_{h,\mathsf{x},\mathsf{y}}) \hat\Gamma_{h,\mathsf{y}}^{-1} \Upsilon_h c_{h,\mathsf{y}}}{\sqrt{c^\prime_{h,\mathsf{x}}\Upsilon_h{\Omega}_{h,\mathsf{x}}\Upsilon_h c_{h,\mathsf{x}}}\sqrt{c^\prime_{h,\mathsf{y}}\Upsilon_h{\Omega}_{h,\mathsf{y}}\Upsilon_h c_{h,\mathsf{y}}}} \\
&\quad+ \frac{c^\prime_{h,\mathsf{x}}\Upsilon_h (\hat\Gamma_{h,\mathsf{x}}^{-1} -\Gamma_{h,\mathsf{x}}^{-1}) \Sigma_{h,\mathsf{x},\mathsf{y}} \hat\Gamma_{h,\mathsf{y}}^{-1} \Upsilon_h c_{h,\mathsf{y}}}{\sqrt{c^\prime_{h,\mathsf{x}}\Upsilon_h{\Omega}_{h,\mathsf{x}}\Upsilon_h c_{h,\mathsf{x}}}\sqrt{c^\prime_{h,\mathsf{y}}\Upsilon_h{\Omega}_{h,\mathsf{y}}\Upsilon_h c_{h,\mathsf{y}}}} \\
&\quad+ \frac{c^\prime_{h,\mathsf{x}}\Upsilon_h \Gamma_{h,\mathsf{x}}^{-1} \Sigma_{h,\mathsf{x},\mathsf{y}}(\hat\Gamma_{h,\mathsf{y}}^{-1} -\Gamma_{h,\mathsf{y}}^{-1}) \Upsilon_h c_{h,\mathsf{y}}}{\sqrt{c^\prime_{h,\mathsf{x}}\Upsilon_h{\Omega}_{h,\mathsf{x}}\Upsilon_h c_{h,\mathsf{x}}}\sqrt{c^\prime_{h,\mathsf{y}}\Upsilon_h{\Omega}_{h,\mathsf{y}}\Upsilon_h c_{h,\mathsf{y}}}}
.
\end{align*}
From the analysis of $\hat\Sigma_{h,\mathsf{x},\mathsf{y}}$, we have
\begin{align*}
\sup_{\mathsf{x},\mathsf{y}\in \mathcal{I}}\left|\frac{c^\prime_{h,\mathsf{x}}\Upsilon_h \hat\Gamma_{h,\mathsf{x}}^{-1} (\hat\Sigma_{h,\mathsf{x},\mathsf{y}} - \Sigma_{h,\mathsf{x},\mathsf{y}}) \hat\Gamma_{h,\mathsf{y}}^{-1} \Upsilon_h c_{h,\mathsf{y}}}{\sqrt{c^\prime_{h,\mathsf{x}}\Upsilon_h{\Omega}_{h,\mathsf{x}}\Upsilon_h c_{h,\mathsf{x}}}\sqrt{c^\prime_{h,\mathsf{y}}\Upsilon_h{\Omega}_{h,\mathsf{y}}\Upsilon_h c_{h,\mathsf{y}}}}\right| = O_{\mathbb{P}}\left( \sqrt{\frac{\log n}{nh^2}} \right).
\end{align*}
For the second term, we have
\begin{align*}
\left|\frac{c^\prime_{h,\mathsf{x}}\Upsilon_h (\hat\Gamma_{h,\mathsf{x}}^{-1} -\Gamma_{h,\mathsf{x}}^{-1}) \Sigma_{h,\mathsf{x},\mathsf{y}} \hat\Gamma_{h,\mathsf{y}}^{-1} \Upsilon_h c_{h,\mathsf{y}}}{\sqrt{c^\prime_{h,\mathsf{x}}\Upsilon_h{\Omega}_{h,\mathsf{x}}\Upsilon_h c_{h,\mathsf{x}}}\sqrt{c^\prime_{h,\mathsf{y}}\Upsilon_h{\Omega}_{h,\mathsf{y}}\Upsilon_h c_{h,\mathsf{y}}}}\right|
&\leq \frac{|c^\prime_{h,\mathsf{x}}\Upsilon_h (\hat\Gamma_{h,\mathsf{x}}^{-1} -\Gamma_{h,\mathsf{x}}^{-1}) \Sigma_{h,\mathsf{x}}^{1/2}|\cdot|c^\prime_{h,\mathsf{y}}\Upsilon_h \hat{\Gamma}_{h,\mathsf{y}}^{-1} \Sigma_{h,\mathsf{y}}^{1/2}|}{\sqrt{c^\prime_{h,\mathsf{x}}\Upsilon_h{\Omega}_{h,\mathsf{x}}\Upsilon_h c_{h,\mathsf{x}}}\sqrt{c^\prime_{h,\mathsf{y}}\Upsilon_h{\Omega}_{h,\mathsf{y}}\Upsilon_h c_{h,\mathsf{y}}}}\\
&= O_{\mathbb{P}}\left( \sqrt{\frac{\log n}{nh^2}} \right).
\end{align*}
The same bound holds for the third term.
\subsection{Proof of Lemma \ref{lem:local regression uniform term 1}}
We decompose \eqref{eq:local regression uniform term 1} as
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}|\text{\eqref{eq:local regression uniform term 1}}| &\leq
\frac{1}{\sqrt{n}}
\underbrace{\left[\sup_{\mathsf{x}\in\mathcal{I}}\left|\frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Gamma}_{h,\mathsf{x}}^{-1}}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}}\right|\right]}_{\text{(I)}}
\underbrace{\left[\sup_{\mathsf{x}\in\mathcal{I}}\left|\frac{1}{n}\sum_{i=1}^n {R((x_i-\mathsf{x})/h)[1-F(x_i)]\frac{1}{h}K(\frac{x_i-\mathsf{x}}{h})} \right|\right]}_{\text{(II)}}.
\end{align*}
As both $\hat\Gamma_{h,\mathsf{x}}$ and $c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}$ are uniformly consistent, term (I) has order
\begin{align*}
\text{(I)} = O_{\mathbb{P}}\left(\sqrt{\frac{1}{h}}\right).
\end{align*}
For (II), we can employ the same technique used to prove Lemma \ref{lem:uniform consistency of Gamma} and show that
\begin{align*}
\text{(II)} &= O_{\mathbb{P}}\left( 1 + \sqrt{\frac{\log n}{nh}} \right) = O_{\mathbb{P}}(1),
\end{align*}
where the leading order in the above represents the mean of ${R((x_i-\mathsf{x})/h)[1-F(x_i)]\frac{1}{h}K(\frac{x_i-\mathsf{x}}{h})}$.
Next, term \eqref{eq:local regression uniform term 2} is bounded by
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}|\eqref{eq:local regression uniform term 2}| \leq
\sqrt{n}
\underbrace{\left[ \sup_{\mathsf{x}\in\mathcal{I}}\left| \frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Gamma}_{h,\mathsf{x}}^{-1}}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}} \right| \right]}_{\text{(I)}}
\underbrace{\left[\sup_{\mathsf{x}\in\mathcal{I}}\left|\frac{1}{n}\sum_{i=1}^n {R((x_i-\mathsf{x})/h)[F(x_i)-\theta(\mathsf{x})^\prime R(x_i-\mathsf{x})]\frac{1}{h}K(\frac{x_i-\mathsf{x}}{h})}\right|\right]}_{\text{(II)}}.
\end{align*}
Employing the same argument used to prove Lemma \ref{lem:local regression uniform term 1}, we have
\begin{align*}
\text{(I)} = O_{\mathbb{P}}\left(\sqrt{\frac{1}{h}}\right).
\end{align*}
To bound term (II), recall that $K(\cdot)$ is supported on $[-1,1]$, meaning that
\begin{align*}
&\ \sup_{\mathsf{x}\in\mathcal{I}}\left|\frac{1}{n}\sum_{i=1}^n {R((x_i-\mathsf{x})/h)[F(x_i)-\theta(\mathsf{x})^\prime R(x_i-\mathsf{x})]\frac{1}{h}K(\frac{x_i-\mathsf{x}}{h})}\right|\\
&= \sup_{\mathsf{x}\in\mathcal{I}}\left|\frac{1}{n}\sum_{i=1}^n {R((x_i-\mathsf{x})/h)[F(x_i)-\theta(\mathsf{x})^\prime R(x_i-\mathsf{x})]\mathds{1}(|x_i-\mathsf{x}|\leq h)\frac{1}{h}K(\frac{x_i-\mathsf{x}}{h})}\right|\\
&\leq \underbrace{\left[\sup_{\mathsf{x}\in\mathcal{I}}\frac{1}{n}\sum_{i=1}^n \left|{R((x_i-\mathsf{x})/h)\frac{1}{h}K(\frac{x_i-\mathsf{x}}{h})}\right|\right]}_{\text{(II.1)}} \underbrace{\left[\sup_{\mathsf{x}\in\mathcal{I}}\sup_{u\in [\mathsf{x}-h,\mathsf{x}+h]}\left|\Big[F(u)-\theta(\mathsf{x})^\prime R(u-\mathsf{x})\Big]\right|\right]}_{\text{(II.2)}}.
\end{align*}
Term (II.2) has the bound $\sup_{\mathsf{x}\in\mathcal{I}}\varrho(h,\mathsf{x})$. Term (II.1) can be bounded by mean and variance calculations and adopting the proof of Lemma \ref{lem:uniform consistency of Gamma}, which leads to
\begin{align*}
\text{(II.1)} = O_{\mathbb{P}}\left(1 + \sqrt{\frac{\log n}{nh}}\right) = O_{\mathbb{P}}(1).
\end{align*}
To show the last conclusion, define the following:
\begin{align*}
u_{ij}(\mathsf{x}) &= \Upsilon_hR(x_j-\mathsf{x})\Big[ \mathds{1}(x_i\leq x_j) -F(x_j) \Big] \frac{1}{h}K\left(\frac{x_j-\mathsf{x}}{h}\right) - \int_{\frac{\mathcal{X}-\mathsf{x}}{h}} R(u)\Big[\mathds{1}(x_i\leq \mathsf{x} + hu) - F(\mathsf{x} + hu)\Big] K\left(u\right) f(\mathsf{x} + hu)\mathrm{d} u,
\end{align*}
then $n^{-2}\sum_{i,j=1, i\neq j}^n u_{ij}(\mathsf{x})$ is a degenerate U-statistic. We rewrite \eqref{eq:local regression uniform term 3} as
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}|\text{\eqref{eq:local regression uniform term 3}}| \leq
\sqrt{n}
\underbrace{\left[ \sup_{\mathsf{x}\in\mathcal{I}}\left| \frac{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Gamma}_{h,\mathsf{x}}^{-1}}{\sqrt{c_{h,\mathsf{x}}^\prime\Upsilon_{h}\hat{\Omega}_{h,\mathsf{x}}\Upsilon_{h}c_{h,\mathsf{x}}}} \right| \right]}_{\text{(I)}}
\underbrace{\left[ \sup_{\mathsf{x}\in\mathcal{I}}\left| \frac{1}{n^2}\sum_{i,j=1,i\neq j}^n u_{ij} \right| \right]}_{\text{(II)}}.
\end{align*}
As before, we have
\begin{align*}
\text{(I)} = O_{\mathbb{P}}\left(\sqrt{\frac{1}{h}}\right).
\end{align*}
Now we consider (II). Let $\mathcal{I}_\varepsilon$ be an $\frac{\varepsilon}{2}$-covering of $\mathcal{I}$, we have
\begin{align*}
\sup_{\mathsf{x}\in\mathcal{I}}\left|\frac{1}{n^2}\sum_{i,j=1, i\neq j}^n u_{ij}(\mathsf{x})\right| &\leq \underbrace{\max_{\mathsf{x}\in\mathcal{I}_\varepsilon}\left|\frac{1}{n^2}\sum_{i,j=1, i\neq j}^n u_{ij}(\mathsf{x})\right|}_{\text{(II.1)}}+ \underbrace{\max_{\mathsf{x}\in\mathcal{I}_\varepsilon,\mathsf{y}\in\mathcal{I},|\mathsf{x}-\mathsf{y}|\leq \varepsilon}\left|\frac{1}{n^2}\sum_{i,j=1, i\neq j}^n \Big(u_{ij}(\mathsf{x})-u_{ij}(\mathsf{y})\Big)\right|}_{\text{(II.2)}}.
\end{align*}
We rely on the concentration inequality in Lemma \ref{lem:ustatistic concentration inequality} for degenerate second order U-statistics. By our assumptions, $A$ can be chosen to be $C_1h^{-1}$ where $C_1$ is some constant that is independent of $\mathsf{x}$. Similarly, $B$ can be chosen to be $C_2\sqrt{n}h^{-1}$ for some constant $C_2$ which is independent of $\mathsf{x}$, and $D$ can be chosen as $C_3nh^{-1/2}$ for some $C_3$ independent of $\mathsf{x}$. Therefore, by setting $\eta=K\log n/\sqrt{n^2h}$ for some large constant $K$, we have
\begin{align*}
\mathbb{P}\left[ \text{(II.1)}\geq \eta \right] &\leq C\frac{1}{\varepsilon}\max_{\mathsf{x}\in\mathcal{I}_\varepsilon}\mathbb{P}\left[ \left|\sum_{i,j=1,i\neq j}^n u_{ij}(\mathsf{x})\right|\geq n^2\eta \right]\\
&\leq C\frac{1}{\varepsilon} \exp\left\{ -\frac{1}{C}\min\left[ \frac{n^2h^{1/2}\eta}{nc_3},\ \left(\frac{n^2h\eta}{n^{1/2}c_2}\right)^{\frac{2}{3}},\ \left( \frac{n^2h\eta}{c_1} \right)^{\frac{1}{2}} \right] \right\}\\
&= C\frac{1}{\varepsilon} \exp\left\{ -\frac{1}{C}\min\left[ \frac{K \log n}{c_3},\ \left(\frac{K\sqrt{nh}\log n}{c_2}\right)^{\frac{2}{3}},\ \left( \frac{K\sqrt{n^2h}\log n}{c_1} \right)^{\frac{1}{2}} \right] \right\}.
\end{align*}
As $\varepsilon$ is at most polynomial in $n$, the above tends to zero for all $K$ large enough, which implies
\begin{align*}
\text{(II.1)} = O_{\mathbb{P}}\left( \frac{\log n}{\sqrt{n^2h}} \right).
\end{align*}
With tedious but still straightforward calculations, it can be shown that
\begin{align*}
\text{(II.2)} &= O_{\mathbb{P}}\left( \frac{\varepsilon}{h} + \frac{\log n}{\sqrt{n^2h}} + \frac{\varepsilon}{h}\frac{\log n}{\sqrt{n^2h}} \right),
\end{align*}
and to match the rates, let $\varepsilon = h\log n/\sqrt{n^2h}$.
\subsection{Proof of Lemma \ref{lem:local regression uniform term 2}}
The proof resembles that of of Lemma \ref{lem:local projection uniform term 1}.
\subsection{Proof of Theorem \ref{thm:strong approximation, local regression}}
The proof resembles that of Theorem \ref{thm:strong approximation, local projection}.
\subsection{Proof of Theorem \ref{thm:feasible uniform approximation local regression}}
The proof resembles that of Theorem \ref{thm:feasible uniform approximation local projection}.
}
\singlespacing
\bibliographystyle{econometrica}
\bibliography{Cattaneo-Jansson-Ma_2021_JoE}
\clearpage
\clearpage