Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
Simple Local Polynomial Density Estimators Supplemental Appendix
abstractThis Supplemental Appendix contains general theoretical results and their proofs, which encompass those discussed in the main paper, discusses additional methodological and technical results, and reports simulation evidence.
\thispagestyle{empty}
\setcounter{tocdepth}{2}
\thispagestyle{empty}
\onehalfspacing
\setcounter{page}{1}
\pagestyle{plain}
\pagestyle{plain}
Setup
We repeat the setup in the main paper for completeness. Recall that $\{ x_i \}_{1\leq i\leq n}$ is a random sample from the cumulative distribution function (hereafter CDF) $F$, supported on $\mathcal{X}=[{x}_{\mathtt{L}},{x}_{\mathtt{U}}]$. Note that it is possible to have ${x}_{\mathtt{L}}=-\infty$ and/or ${x}_{\mathtt{U}}=\infty$. We will assume both ${x}_{\mathtt{L}}$ and ${x}_{\mathtt{U}}$ are finite, to facilitate discussion on boundary estimation issues.
Define the empirical distribution function (hereafter EDF)
align*[align* omitted — 70 chars of source]
Note that in the main paper, we use $\hat{F}(\cdot)$ to denote the above EDF We avoid such notation in this Supplemental Appendix, and instead use $\tilde{F}(\cdot)$, because a (smoothed) CDF estimator can be obtained from our local polynomial approach.
Given $p\in\mathbb{N}$, our local polynomial distribution estimator is defined as
align*[align* omitted — 196 chars of source]
where $\mathbf{r}_p(u)=[1, u, u^2, \cdots, u^p]$ is a (one-dimensional) polynomial expansion; $K$ is a kernel function whose properties are to be specified later; $h=h_n$ is a bandwidth sequence. The estimator, $\hat{\boldsymbol{\beta}}_p(x)$, is motivated as a local Taylor series expansion, hence the target parameter is (i.e., the population counterpart, assuming exists)
align*[align* omitted — 139 chars of source]
Therefore, we also write
align*[align* omitted — 169 chars of source]
or equivalently, $\hat{F}_p^{(v)} = v! \mathbf{e}_v^\prime \hat{\boldsymbol{\beta}}_p(x)$, provided that $v\leq p$, and $\mathbf{e}_v$ is the $(v+1)$-th unit vector of $\mathbb{R}^{p+1}$. (The subscript $p$ is omitted in the main paper to economize notation.) We also use $f=F^{(1)}$ to denote the corresponding probability density function (hereafter PDF) for convenience.
The estimator has the following matrix form:
align*[align* omitted — 310 chars of source]
where
$\mathbf{K}_h$ is a diagonal matrix collecting $\{ h^{-1}K((x_i-x)/h) \}_{1\leq i\leq n}$, and $\mathbf{Y}$ is a column vector collecting $\{ \tilde{F}(x_i) \}_{1\leq i\leq n}$. We also adopt the convention $K_h(u) = h^{-1}K(u/h)$.
In this Supplemental Appendix, we use $n$ to denote sample size, and limits are taken with $n\to \infty$, unless otherwise specified. The standard Euclidean norm is denoted by $|\cdot|$, and other norms will be defined at their first appearances. Maximum and minimum of two real numbers $a$ and $b$ are denoted by $a\vee b$ and $a\wedge b$, respectively. For sequence of numbers (or random variables), $a_n\precsim b_n$ implies $\limsup_n |a_n/b_n|$ is finite, and $a_n\asymp b_n$ implies both directions. The notation $a_n\precsim_\mathbb{P} b_n$ is used to denote that $|a_n/b_n|$ is asymptotically tight: $\limsup_{\varepsilon\uparrow \infty}\limsup_{n}\mathbb{P}[|a_n/b_n|\geq \varepsilon]=0$. $a_n\asymp_{\mathbb{P}} b_n$ implies both $a_n\precsim_\mathbb{P} b_n$ and $b_n\precsim_\mathbb{P} a_n$. When $b_n$ is a sequence of nonnegative numbers, $a_n=O(b_n)$ is sometimes used for $a_n\precsim b_n$, so does $a_n=O_{\mathbb{P}}(b_n)$. For probabilistic convergence, we use $\to_\mathbb{P}$ for convergence in probability and $\rightsquigarrow$ for weak convergence (convergence in distribution). Standard normal distribution is denoted as $\mathcal{N}(0,1)$, with CDF $\Phi$ and PDF $\phi$. Throughout, we use $C$ to denote generic constants which do not depend on sample size. The exact value can change given the context.
Overview of Main Results
In this subsection, we give an overview of our results, including a (first order) mean squared error (hereafter MSE) expansion, and asymptotic normality. Fix some $v\geq 1$ and $p$, we have the following:
align*[align* omitted — 211 chars of source]
The previous result gives MSE expansion for derivative estimators, $1\leq v\leq p$, but not for $v=0$. With $v=0$, $\hat{F}_p(x)$ is essentially a smoothed EDF, which estimates the CDF $F(x)$. Since $F(x)$ is $\sqrt{n}$-estimable, one should be expected that the estimated distribution function will have very different properties compared to the estimated derivatives. Indeed, we have
align*[align* omitted — 189 chars of source]
There is another complication, however, when $x$ is in the boundary region. For a drifting sequence $x$ in the boundary region, the EDF $\tilde{F}(x)$ is “super-consistent” in the sense that it converges at rate $\sqrt{h/n}$. The reason is that when $x$ is near ${x}_{\mathtt{L}}$ or ${x}_{\mathtt{U}}$, $\tilde{F}(x)$ is essentially estimating 0 or 1, and the variance, $F(x)(1-F(x))$ vanishes asymptotically, giving rise to the additional factor $\sqrt{h}$. This is shared by our estimator: for $v=0$ and $x$ in the boundary region, the CDF estimator $\tilde{F}_p(x)$ is super-consistent, with $\mathcal{V}_{p,0}(x)\asymp h$.
Also note that for the MSE expansion, we provide not only the first order bias but also the second order bias. The second order bias will be used for bandwidth selection, since it is well-known that in some cases the first order bias can vanish. (More precisely, when $x$ is an interior evaluation point and $p-v$ is even. See, for example, Fan-Gijbels_1996_Book.)
The MSE expansion provides the rate of convergence of our estimator. The following shows that, under suitable regularity conditions, they are also asymptotically normal. Again first consider $v\geq 1$.
align*[align* omitted — 169 chars of source]
provided that the bandwidth is not too large, so that after scaling, the remaining bias does not feature in first-order asymptotics. For $v=0$, i.e. the smoothed EDF, we have
align*[align* omitted — 157 chars of source]
where we moved the variance $\mathcal{V}_{p,0}(x)$ as a scaling factor in the above display, to encompass the situation where $x$ lies near boundaries.
Some Matrices
In this subsection we collect some matrices which will be used throughout this Supplemental Appendix. They show up in asymptotic results as components of the (leading) bias and variance. Note that $x$ can be either a fixed point, or it can be a drifting sequence to capture the issue of estimation and inference in boundary regions. For the latter case, $x$ takes the form $x={x}_{\mathtt{L}}+ch$ or $x={x}_{\mathtt{U}}-ch$ for some $c\in[0,1)$.
Define
align*[align* omitted — 757 chars of source]
Later we will assume that the kernel function $K$ is supported on $[-1, 1]$, hence with a shrinking bandwidth sequence $h\downarrow 0$, the region of integration in the above display can be replaced by
center[center omitted — 277 chars of source]
Since we do not allow ${x}_{\mathtt{L}}={x}_{\mathtt{U}}$, no drifting sequence $x$ can be in both lower and upper boundary regions, at least in large samples.
Large Sample Properties
Assumptions
In this section we give assumptions, preliminary lemmas and our main results. Other assumptions specific to certain results will be given in corresponding sections.
Let $\mathcal{O}$ be a connected subset of $\mathbb{R}$ with nonempty interior, $\mathcal{C}^s(\mathcal{O})$ denotes functions that are at least $s$-times continuously differentiable in the interior of $\mathcal{O}$, and that the derivatives can be continuously extended to the boundary of $\mathcal{O}$.
Assumption[DGP]
$ $\\
$\{x_i\}_{1\leq i\leq n}$ is a random sample from distribution $F$, supported on $\mathcal{X}=[{x}_{\mathtt{L}},{x}_{\mathtt{U}}]$. Further, $F\in\mathcal{C}^{\alpha_x}(\mathcal{X})$ for some $\alpha_x\geq 1$, and $f(x)=F^{(1)}(x)>0$ for all $x\in\mathcal{X}$.
Assumption[Kernel]
$ $\\
The kernel function $K(\cdot)$ is nonnegative, symmetric, and belongs to $\mathcal{C}^0([-1,1])$. Further, it integrates to one: $\int_{\mathbb{R}}K(u)\mathrm{d} u=1$.
Preliminary Lemmas
We first consider the object $\mathbf{X}_h^{\prime}\mathbf{K}_h\mathbf{X}_h/n$
lemAssume Assumptions (ref) and (ref) hold, $h\to 0$ and $nh\to \infty$. Then
\begin{align*}
\frac{1}{n}\mathbf{X}_h^\prime\mathbf{K}_h\mathbf{X}_h = f(x)\mathbf{S}_{p,x} + o(1) + O_{\mathbb{P}}\left( 1/\sqrt{nh} \right).
\end{align*}
Lemma (ref) shows that the matrix $\mathbf{X}_h^{\prime}\mathbf{K}_h\mathbf{X}_h/n$ is asymptotically invertible. Also note that this result covers both interior and boundary evaluation point $x$, and depending on the nature of $x$, the exact form of $\mathbf{S}_{p,x}$ differs.
With simple algebra, one has
align*[align* omitted — 275 chars of source]
and the following gives a further decomposition of the “numerator.”
align*[align* omitted — 837 chars of source]
The first part represents the smoothing bias, and the second part can be analyzed as a sample average. The real challenge comes from the third term, which can have a nonnegligible (first order) contribution. We further decompose it as
align*[align* omitted — 462 chars of source]
As a result,
align*[align* omitted — 1,085 chars of source]
To provide intuition for the above decomposition, the smoothing bias is a typical feature of nonparametric estimators; leave-in bias arises since each observation is used twice, in constructing the EDF $\tilde{F}$ and as a design point (that is, $\tilde{F}$ has to be evaluated at $x_i$); and a second order U-statistic shows up because the “dependent variable,” $\mathbf{Y}$, is estimated, which leads to double summation.
We first analyze the bias terms.
lemAssume Assumptions (ref) and (ref) hold with $\alpha_x\geq p+1$, $h\to 0$ and $nh\to\infty$. Then
\begin{align*}
\hat{\mathbf{B}}_{\mathtt{S}} &= h^{p+1}\frac{F^{(p+1)}(x)f(x)}{(p+1)!}\mathbf{c}_{p,x} + o_{\mathbb{P}}(h^{p+1}),\qquad \hat{\mathbf{B}}_{\mathtt{LI}} = O_{\mathbb{P}}\left(n^{-1}\right).
\end{align*}
By imposing additional smoothness, it is also possible to characterize the next term in the smoothing bias, which has order $h^{p+2}$. We report the higher order bias in a later section as it is used for bandwidth selection.
Next we consider the “influence function” part, $\hat{\mathbf{L}}$. This term is crucial in the sense that (under suitable conditions so that $\hat{\mathbf{R}}$ becomes negligible) it determines the asymptotic variance of our estimator, and with correct scaling, it is asymptotically normally distributed.
lemAssume Assumptions (ref)and (ref) hold with $\alpha_x\geq 2$, $h\to 0$ and $nh\to\infty$. Define the scaling matrix
\begin{align*}
\mathbf{N}_x &= \begin{cases}
\mathrm{diag}\Big\{ 1,\qquad h^{-1/2},\ h^{-1/2},\ \cdots,\ h^{-1/2} \Big\} & $x$ interior,\\
\mathrm{diag}\Big\{ h^{-1/2},\ h^{-1/2},\ h^{-1/2},\ \cdots,\ h^{-1/2} \Big\} & $x$ boundary,\\
\end{cases}
\end{align*}
then
\begin{align*}
\sqrt{n}\mathbf{N}_x \Big[f(x)\mathbf{S}_{p,x}\Big]^{-1}\hat{\mathbf{L}} \rightsquigarrow \mathcal{N}(\mathbf{0},\ \mathsf{V}_{p,x}),
\end{align*}
with
\begin{align*}
\mathsf{V}_{p,x} = \begin{cases}
F(x)(1-F(x))\mathbf{e}_0\mathbf{e}_0^\prime + f(x) (\mathbf{I}-\mathbf{e}_0\mathbf{e}_0^\prime)\mathbf{S}_{p,x}^{-1}\boldsymbol{\Gamma}_{p,x}\mathbf{S}_{p,x}^{-1}(\mathbf{I}-\mathbf{e}_0\mathbf{e}_0^\prime) & $x$ interior\\
f(x) \left(\mathbf{S}_{p,x}^{-1}\boldsymbol{\Gamma}_{p,x}\mathbf{S}_{p,x}^{-1}+c\mathbf{e}_0\mathbf{e}_0^\prime\right) & x={x}_{\mathtt{L}}+ch\\
f(x) \left(\mathbf{S}_{p,x}^{-1}\boldsymbol{\Gamma}_{p,x}\mathbf{S}_{p,x}^{-1}+ c\mathbf{e}_0\mathbf{e}_0^\prime-(\mathbf{e}_1\mathbf{e}_0^\prime + \mathbf{e}_0\mathbf{e}_1^\prime)\right) & x={x}_{\mathtt{U}}-ch.
\end{cases}
\end{align*}
The scaling matrix depends on whether the evaluation point is located in the interior or boundary, which is a unique feature of our estimator. To see the intuition, consider an interior point $x$, and recall that the first element of $\hat{\boldsymbol{\beta}}_p(x)$ is the smoothed EDF, which is $\sqrt{n}$-estimable. Therefore, the property of $\hat{F}_p(x)$ is very different from those of the estimated density and higher order derivatives.
When $x$ is either in the lower or upper boundary region, $\hat{F}_p(x)$ essentially estimates 0 or 1, respectively, hence it is super-consistent in the sense that it converges even faster than $1/\sqrt{n}$. In this case, the leading $1/\sqrt{n}$-variance vanishes, and higher order residual noise dominates, which makes $\hat{F}_p(x)$ no longer independent of the estimated density and derivatives, justifying the formula of boundary evaluation points.
Finally we consider the second order U-statistic component.
lemAssume Assumptions (ref) and (ref) hold, $h\to 0$ and $nh\to\infty$. Then
\begin{align*}
\mathbb{V}[\hat{\mathbf{R}}] &= \frac{2}{n^2h}f(x)F(x)(1-F(x))\mathbf{T}_{p,x} + O(n^{-2}).
\end{align*}
In particular, when $x$ is in the boundary region, the above has order $O(n^{-2})$.
Main Results
In this section we provide two main results, one on asymptotic normality, and the other on standard error.
thm[Asymptotic Normality]
Assume Assumptions (ref) and (ref) hold with $\alpha_x\geq p+1$ for some integer $p\geq 0$. Further $h\to 0$, $nh^2\to \infty$ and $nh^{2p+1}=O(1)$. Then
\begin{align*}
\sqrt{nh^{2v-1}}\Big( \hat{F}_p^{(v)}(x) - F^{(v)}(x) - h^{p+1-v}\mathcal{B}_{p,v}(x) \Big) &\rightsquigarrow \mathcal{N}\Big(0,\ \mathcal{V}_{p,v}(x)\Big),\qquad 1\leq v\leq p,\\
\sqrt{\frac{n}{\mathcal{V}_{p,0}(x)}} \Big( \hat{F}_p(x) - F(x) - h^{p+1}\mathcal{B}_{p,0}(x) \Big) &\rightsquigarrow \mathcal{N}\Big(0,\ 1\Big).
\end{align*}
The constants are
\begin{align*}
\mathcal{B}_{p,v}(x) &= v! \frac{F^{(p+1)}(x)}{(p+1)!} \mathbf{e}_v^\prime \mathbf{S}_{p,x}^{-1} \mathbf{c}_{p,x},
\end{align*}
and
\begin{align*}
\mathcal{V}_{p,v}(x) &= \begin{cases}
(v!)^2 f(x) \mathbf{e}_v^\prime \mathbf{S}_{p,x}^{-1}\boldsymbol{\Gamma}_{p,x}\mathbf{S}_{p,x}^{-1} \mathbf{e}_v \qquad& 1\leq v\leq p\\
F(x)(1 - F(x)) & $v=0$, $x$ interior\\
hf(x) \left(\mathbf{e}_0^\prime\mathbf{S}_{p,x}^{-1}\boldsymbol{\Gamma}_{p,x}\mathbf{S}_{p,x}^{-1}\mathbf{e}_0+c\right) & $v=0$, $x={x}_{\mathtt{L}}+ch$ or ${x}_{\mathtt{U}}-ch$.
\end{cases}
\end{align*}
remark[On $nh^{2p+1}=O(1)$]
This condition ensures that higher order bias, after scaling, is asymptotically negligible. \qquad \qquad $\parallel$
remark[On $nh^2\to \infty$]
This condition ensures that the second order U-statistic, $\hat{\mathbf{R}}$, has smaller order compared to $\hat{\mathbf{L}}$. Note that this condition can be dropped for boundary $x$ or when the parameter of interest is the CDF $\hat{F}_p$. \qquad \qquad $\parallel$
Now we provide a standard error, which is also boundary adaptive. Given the formula in Theorem (ref), it is possible to estimate the asymptotic variance by plugging in unknown quantities regarding the data generating process. For example consider $\mathcal{V}_{p,1}(x)$ for the estimated density. Assume the researcher knows the location of the boundary ${x}_{\mathtt{L}}$ and ${x}_{\mathtt{U}}$, the matrices $\mathbf{S}_{p,x}$ and $\boldsymbol{\Gamma}_{p,x}$ can be constructed with numerical integration, since they are related to features of the kernel function, not the data generating process. The unknown density $f(x)$ can also be replaced by its estimate, as long as $p\geq 1$.
Another approach is to utilize the decomposition of the estimator, in particular the $\hat{\mathbf{L}}$ term. To introduce our variance estimator, we make the following definitions.
align*[align* omitted — 487 chars of source]
Following is the main result regarding variance estimation. It is automatic and fully-adaptive, in the sense that no knowledge about the boundary location is needed.
thm[Variance Estimation]\ \\
Assume Assumptions (ref) and (ref) hold with $\alpha_x\geq p+1$ for some integer $p\geq 0$. Further $h\to 0$, $nh^2\to \infty$ and $nh^{2p+1}=O(1)$. Then
\begin{align*}
\hat{\mathcal{V}}_{p,v}(x) \equiv (v!)^2 \mathbf{e}_v^\prime \mathbf{N}_x \hat{\mathbf{S}}_{p,x}^{-1}\hat{\boldsymbol{\Gamma}}_{p,x}\hat{\mathbf{S}}_{p,x}^{-1} \mathbf{N}_x \mathbf{e}_v\ \to_\mathbb{P}\ \mathcal{V}_{p,v}(x).
\end{align*}
Define the standard error as
\begin{align*}
\hat{\sigma}_{p,v}(x) &\equiv (v!) \sqrt{\frac{1}{nh^{2v}} \mathbf{e}_v^\prime \hat{\mathbf{S}}_{p,x}^{-1}\hat{\boldsymbol{\Gamma}}_{p,x}\hat{\mathbf{S}}_{p,x}^{-1} \mathbf{e}_v},
\end{align*}
then
\begin{align*}
\hat{\sigma}_{p,v}(x)^{-1} \Big( \hat{F}_p^{(v)}(x) - F^{(v)}(x) - h^{p+1-v}\mathcal{B}_{p,v}(x) \Big) \rightsquigarrow \mathcal{N}\Big( 0,\ 1\Big).
\end{align*}
Bandwidth Selection
In this section we consider the problem of constructing MSE-optimal bandwidth for our local polynomial regression-based distribution estimators. We focus exclusively on the case $v\geq 1$, hence the object of interest will be either the density function or derivatives thereof. Valid bandwidth choice for the distribution function $\hat{F}_p(x)$ is also an interesting topic, but difficulty arises since it is estimated at the parametric rate. We will briefly mention MSE expansion of the estimated CDF at the end.
For Density and Derivatives Estimates ($v\geq 1$)
Consider some $1\leq v\leq p$, the following lemma gives finer characterization of the bias.
lemAssume Assumptions (ref) and (ref) hold with $\alpha_x\geq p+2$, $h\to 0$ and $nh^3\to \infty$. Then the leading bias of $\hat{F}_p^{(v)}(x)$ is
\begin{align*}
h^{p+1-v}\mathcal{B}_{p,v}(x) &= h^{p+1-v}\left\{
\frac{F^{(p+1)}(x)}{(p+1)!}v!\mathbf{e}_v^\prime \mathbf{S}_{p,x}^{-1}\mathbf{c}_{p,x} +
h\left( \frac{F^{(p+2)}(x)}{(p+2)!} + \frac{F^{(p+1)}(x)}{(p+1)!}\frac{F^{(2)}(x)}{f(x)} \right)v!\mathbf{e}_v^\prime \mathbf{S}_{p,x}^{-1}\tilde{\mathbf{c}}_{p,x} \right\}.
\end{align*}
The above lemma is a refinement of Lemma (ref) and (ref), and characterizes the higher-order bias. To see its necessity, we note that when $p-v$ is even and $x$ is an interior evaluation point, the leading bias is zero. This is because $\mathbf{e}_v^\prime \mathbf{S}_{p,x}^{-1}\mathbf{c}_{p,x}$ is zero, which is explained in Fan-Gijbels_1996_Book. Except for rare cases such as $F^{(p+1)}(x)=0$ or $F^{(p+2)}(x)=0$, we have
center[center omitted — 233 chars of source]
Note that for boundary evaluation points, the leading bias never vanishes.
The leading variance is also characterized by Theorem (ref), and we reproduce it here:
align*[align* omitted — 186 chars of source]
The MSE-optimal bandwidth is defined as a minimizer of the following
align*[align* omitted — 133 chars of source]
Given the discussion we had earlier on the bias, it is easy to see that the MSE-optimal bandwidth has the following asymptotic order:
center[center omitted — 272 chars of source]
Again only the case where $p-v$ is even and $x$ is interior needs special attention.
There are two notions of bandwidth consistency. Let $h$ be some non-stochastic bandwidth sequence, and $\hat{h}$ be an estimated bandwidth. Then $\hat{h}$ is consistent in rate if $\hat{h}\asymp h$ (in most cases it is even true that $\hat{h}/h\to_\mathbb{P} C\in(0,\infty)$). And $\hat{h}$ is consistent in rate and constant if $\hat{h}/h\to_\mathbb{P} 1$.
To construct consistent bandwidth, either rate consistent or consistent in both rate and constant, we need estimates of both the bias and variance. The variance part is relatively easy, as we have already demonstrated in Theorem (ref):
align*[align* omitted — 96 chars of source]
where $\ell$ is some preliminary bandwidth used to construct $\hat{\sigma}_{p,v}(x)$.
To introduce our bias estimate, first assume there are consistent estimators for $F^{(p+1)}(x)$ and $F^{(p+2)}(x)$, denoted by $\hat{F}^{(p+1)}(x)$ and $\hat{F}^{(p+2)}(x)$. They can be obtained, for example, using our local polynomial regression-based approach, or can be constructed with some reference model (such as the normal distribution). The critical step is to obtain consistent estimators of the matrices, which are given in the following lemma.
lemAssume Assumptions (ref) and (ref) hold, $\ell\to 0$ and $n\ell\to \infty$. Then
\begin{align*}
\widehat{\mathbf{S}_{p,x}^{-1}\mathbf{c}_{p,x}}
&= \left(\frac{1}{n}\sum_i\mathbf{r}_p\left(\frac{x_i-x}{\ell}\right)\mathbf{r}_p\left(\frac{x_i-x}{\ell}\right)^\prime K_\ell(x_i-x) \right)^{-1}\left( \frac{1}{n}\sum_i \left(\frac{x_i-x}{\ell}\right)^{p+1}\mathbf{r}_p\left(\frac{x_i-x}{\ell}\right)K_\ell(x_i-x) \right)\\
&\qquad \qquad \qquad \to_\mathbb{P} \mathbf{S}_{p,x}^{-1}\mathbf{c}_{p,x},
\end{align*}
and
\begin{align*}
\widehat{\mathbf{S}_{p,x}^{-1}\tilde{\mathbf{c}}_{p,x}}
&= \left(\frac{1}{n}\sum_i\mathbf{r}_p\left(\frac{x_i-x}{\ell}\right)\mathbf{r}_p\left(\frac{x_i-x}{\ell}\right)^\prime K_\ell(x_i-x) \right)^{-1}\left( \frac{1}{n}\sum_i \left(\frac{x_i-x}{\ell}\right)^{p+2}\mathbf{r}_p\left(\frac{x_i-x}{\ell}\right)K_\ell(x_i-x) \right)\\
&\qquad \qquad \qquad \to_\mathbb{P} \mathbf{S}_{p,x}^{-1}\tilde{\mathbf{c}}_{p,x}.
\end{align*}
Note that we used different notation, $\ell$, as it corresponds to a preliminary bandwidth. Define
align*[align* omitted — 296 chars of source]
and assume that $\hat{\sigma}_{p,v}(x)$ is constructed using the preliminary bandwidth $\ell$. Then
align*[align* omitted — 156 chars of source]
We make some remarks here.
remark[Preliminary bandwidth $\ell$]
The optimization argument $h$ enters the RHS of the previous display in three places. First it is part of the variance component, by $1/h^{2v-1}$. Second it shows as a multiplicative factor of the bias component, $h^{2p-2v+2}$. Finally within the definition of $\hat{\mathcal{B}}_{p,v}(x)$, there is another multiplicative $h$, in front of the higher order bias.
The preliminary bandwidth $\ell$, serves a different role. It is used to estimate the variance and bias components. Of course one can use different preliminary bandwidths for $\hat{\sigma}_{p,v}(x)$, $\widehat{\mathbf{S}_{p,x}^{-1}\mathbf{c}_{p,x}}$ and $\widehat{\mathbf{S}_{p,x}^{-1}\tilde{\mathbf{c}}_{p,x}} $.
\qquad \qquad $\parallel$
remark[Consistent bias estimator]
The bias estimator we proposed, $h^{p-v+1}\hat{\mathcal{B}}_{p,v}(x)$, is consistent in rate for the true leading bias, but not necessarily in constant. Compare $\hat{\mathcal{B}}_{p,v}(x)$ and ${\mathcal{B}}_{p,v}(x)$, it is easily seen that the term involving $F^{(p+1)}(x)F^{(2)}(x)/f(x)$ is not captured. To capture this term, we need one additional nonparametric estimator for $F^{(2)}(x)$. This is indeed feasible, and one can employ our local polynomial regression-based estimator for this purpose.
\qquad \qquad $\parallel$
thm[Consistent bandwidth]
Let $1\leq v\leq p$. Assume the preliminary bandwidth $\ell$ is chosen such that $nh^{2v-1}\hat{\sigma}_{p,v}(x)^2/\mathcal{V}_{p,v}(x) \to_\mathbb{P} 1$, $\widehat{\mathbf{S}_{p,x}^{-1}\mathbf{c}_{p,x}}\to_\mathbb{P} {\mathbf{S}_{p,x}^{-1}\mathbf{c}_{p,x}}$, and $\widehat{\mathbf{S}_{p,x}^{-1}\tilde{\mathbf{c}}_{p,x}} \to_\mathbb{P} {\mathbf{S}_{p,x}^{-1}\tilde{\mathbf{c}}_{p,x}}$. Under the conditions of Lemma (ref) and Theorem (ref):
\begin{itemize}
• If either $x$ is in boundary regions or $p-v$ is odd, let $\hat{F}^{(p+1)}(x)$ be consistent for $F^{(p+1)}\neq 0$. Then
\begin{align*}
\frac{\hat{h}_{ p,v}(x)}{{h}_{ p,v}(x)} \to_\mathbb{P} 1.
\end{align*}
• If $x$ is in interior and $p-v$ is even, let $\hat{F}^{(p+2)}(x)$ be consistent for $F^{(p+2)}\neq 0$. Further assume $nh^3\to 0$ and ${h}_{ p,v}(x)$ is well-defined. Then
\begin{align*}
\frac{\hat{h}_{ p,v}(x)}{{h}_{ p,v}(x)} \to_\mathbb{P} C\in(0,\infty).
\end{align*}
\end{itemize}
For CDF Estimate ($v=0$)
In this subsection we mention briefly how to choose bandwidth for the CDF estimate, $\hat{F}_p^{(0)}(x)\equiv \hat{F}_p(x)$. We assume $x$ is in interior. Previous discussions on bias remains to apply:
align*[align* omitted — 319 chars of source]
which means the bias of $\hat{F}_p(x)$ has order $h^{p+1}$ if either $x$ is boundary or $p$ is odd, and $h^{p+2}$ otherwise. Difficulty arises since the CDF estimator has leading variance of order
align*[align* omitted — 88 chars of source]
which cannot be used for bandwidth selection, because the above is proportional to the bandwidth (i.e., there is no bias-variance trade-off).
The trick is to use a higher order variance term. Recall that the local polynomial regression-based estimator is essentially a second order U-statistic, which is then decomposed into two terms, a linear term $\hat{\mathbf{L}}$ and a quadratic term $\hat{\mathbf{R}}$, where the latter is a degenerate second-order U-statistic. The variance of the quadratic term $\hat{\mathbf{R}}$ has been ignored so far, as it is negligible compared to the variance of the linear term. For the CDF estimator, however, it is the variance of this quadratic term that leads to a bias-variance trade-off. The exact form of this variance is given in Lemma (ref). With this additional variance term included, we have (with some abuse of notation)
align*[align* omitted — 139 chars of source]
Provided $x$ is an interior point, the additional variance term increases as the bandwidth shrinks. As a result, a MSE-optimal bandwidth for $\hat{F}_p(x)$ is well-defined, and estimating this bandwidth is also straightforward.
center[center omitted — 248 chars of source]
What if $x$ is in a boundary region? Then the MSE-optimal bandwidth for $\hat{F}_p(x)$ is not well defined. The leading variance now takes the form $h/n + 1/n^2$, which is proportional to the bandwidth. (This is not surprising, since for boundary $x$ the CDF is known, and a very small bandwidth gives a super-consistent estimator.). Although MSE-optimal bandwidth for $\hat{F}_p(x)$ is not well-defined for boundary $x$, it is still feasible to minimize the empirical MSE. To see how this works, one first estimate the bias term and variance term with some preliminary bandwidth $\ell$, leading to $\hat{\mathcal{B}}_{p,0}(x)$ and $\hat{\mathcal{V}}_{p,0}(x)$. Then the MSE-optimal bandwidth can be constructed by minimizing the empirical MSE. Under regularity conditions, $\hat{\mathcal{B}}_{p,0}(x)$ will converge to some nonzero constant, while, if $x$ is boundary, $\hat{\mathcal{V}}_{p,0}(x)$ has order $\ell$, the same as the preliminary bandwidth. Then the MSE-optimal bandwidth constructed in this way will have the following order:
center[center omitted — 306 chars of source]
Note that the preliminary bandwidth enters the rate of $\hat{h}_{ p,0}(x)$ for boundary $x$, because it determines the rate at which the variance estimator $\hat{\mathcal{V}}_{p,0}(x)$ vanishes. Although this estimated bandwidth is not consistent for any well-defined object, it can be useful in practice, and it reflects the fact that for boundary $x$ it is appropriate to use bandwidth shrinks fast when the object of interest is the CDF
Application to Manipulation Testing
We devote this section to density discontinuity (manipulation) tests in regression discontinuity designs. Assume there is a natural (and known) partition of the support $\mathcal{X}=[{x}_{\mathtt{L}},{x}_{\mathtt{U}}] = [{x}_{\mathtt{L}},\bar{x})\cup [\bar{x}, {x}_{\mathtt{U}}] = \mathcal{X}_-\cup \mathcal{X}_+$, and the regularity conditions we imposed so far are satisfied on each of the partitions, $\mathcal{X}_-$ and $\mathcal{X}_+$. To be precise, assume the distribution $F$ is continuously differentiable to a certain order on each of the partitions, but the derivatives are not necessarily continuous across the cutoff $\bar{x}$. In this case consistent estimates of the densities (and derivatives thereof) require fitting local polynomials separately on each sides of $\bar{x}$. Alternatively, one can use the joint estimation framework introduced below.
Unrestricted Model
By an unrestricted model with cutoff $\bar{x}$, we consider the following polynomial basis $\mathbf{r}_p$
align*[align* omitted — 269 chars of source]
The following two vectors will arise later, which we give the definition here:
align*[align* omitted — 205 chars of source]
Also we define the vectors to extract the corresponding derivatives
align*[align* omitted — 194 chars of source]
With the above definition, the estimator at the cutoff is\footnote{The EDF is defined with the whole sample as before: $\tilde{F}(u)=n^{-1}\sum_i \mathds{1}[x_i\leq u]$.}
align*[align* omitted — 188 chars of source]
We assume the same bandwidth is used below and above the cutoff to avoid cumbersome notation. Generalizing to using different bandwidths is straightforward. Other notations (for example $\mathbf{X}$ and $\mathbf{X}_h$) are redefined similarly, with the scaling matrix $\mathbf{H}$ adjusted so that $\mathbf{H}^{-1}\mathbf{r}_p(u) = \mathbf{r}_p(h^{-1}u)$ is always true. we denote the estimates by
align*[align* omitted — 188 chars of source]
remark[Separate estimation]
An alternative implementation is to apply our local polynomial-based estimator separately to the two samples, one with observations below the cutoff, and the other with observations above the cutoff. To be precise, let $\tilde{F}_{-}(\cdot)$ and $\tilde{F}_{+}(\cdot)$ be the empirical distribution functions constructed by the two samples. That is,
\begin{align*}
\tilde{F}_{-}(x) &= \frac{1}{n_-}\sum_{i:\ x_i<\bar{x}}\mathds{1}[x_i\leq x],\qquad \tilde{F}_{+}(x) = \frac{1}{n_+}\sum_{i:\ x_i\geq\bar{x}}\mathds{1}[x_i\leq x],
\end{align*}
where $n_-$ and $n_+$ denote the size of the two samples, respectively. The the local polynomial approach, applied to $\tilde{F}_{-}(\cdot)$ and $\tilde{F}_{+}(\cdot)$ separately, will yield two sets of estimates, which we denote by $\hat{F}_{p,-}^{(v)}(\bar{x})$ and $\hat{F}_{p,+}^{(v)}(\bar{x})$. To see the relation between joint and separate estimations, we note the following (which can be easily seen using least squares algebra)
\begin{alignat*}{3}
v=0 &\qquad \hat{F}_{p,-}(\bar{x}) &&= \frac{n}{n_-}\hat{F}_{p}(\bar{x}-),\quad &&\hat{F}_{p,+}(\bar{x}) = \frac{n}{n_+}\hat{F}_{p}(\bar{x}+) - \frac{n_-}{n_+}\\
v\geq 1 &\qquad \hat{F}_{p,-}^{(v)}(\bar{x}) &&= \frac{n}{n_-}\hat{F}_{p}^{(v)}(\bar{x}-),\quad &&\hat{F}_{p,+}^{(v)}(\bar{x}) = \frac{n}{n_+}\hat{F}_{p}^{(v)}(\bar{x}+).
\end{alignat*}
The difference comes from the fact that by separate estimation, one obtains estimates of the conditional CDF and the derivatives.
\qquad \qquad $\parallel$
In the following lemmas, we will give asymptotic results for the joint estimation problem. Proofs are omitted.
lemLet Assumptions of Lemma (ref) hold separately on $\mathcal{X}_-$ and $\mathcal{X}_+$, then
\begin{align*}
\nonumber\frac{1}{n}\mathbf{X}_h'\mathbf{K}_h\mathbf{X}_h &= f(\bar{x}-)\mathbf{S}_{-,p} +f(\bar{x}+)\mathbf{S}_{+,p} + O\Big(h\Big) + O_{\mathbb{P}}\Big(1/\sqrt{nh}\Big),
\end{align*}
where
\begin{align*}
\mathbf{S}_{-,p} = \int_{-1}^0 \mathbf{r}_{-,p}(u)\mathbf{r}_{-,p}(u)'K(u) \mathrm{d} u,\qquad \mathbf{S}_{+,p} = \int_0^1 \mathbf{r}_{+,p}(u)\mathbf{r}_{+,p}(u)'K(u) \mathrm{d} u.
\end{align*}
Again we decompose the estimator into four terms, namely $\hat{\mathbf{B}}_{\mathtt{LI}}$, $\hat{\mathbf{B}}_{\mathtt{S}}$, $\hat{\mathbf{L}}$ and $\hat{\mathbf{R}}$.
lemLet Assumptions of Lemma (ref) hold separately on $\mathcal{X}_-$ and $\mathcal{X}_+$, then
\begin{align*}
\hat{\mathbf{B}}_{\mathtt{S}}&= h^{p+1} \left\{ \frac{F^{(p+1)}(\bar{x}-)f(\bar{x}-)}{(p+1)!}\mathbf{c}_{-,p} + \frac{F^{(p+1)}(\bar{x}+)f(x+)}{(p+1)!}\mathbf{c}_{+,p} \right\}+ o_{\mathbb{P}}(h^{p+1}),\qquad \hat{\mathbf{B}}_{\mathtt{LI}} = O_{\mathbb{P}}\left(\frac{1}{n}\right),
\end{align*}
where
\begin{align*}
\mathbf{c}_{-,p} = \int_{-1}^{0} u^{p+1}\mathbf{r}_{-,p}(u)K(u)\mathrm{d} u,\qquad \mathbf{c}_{+,p} = \int_{0}^{1} u^{p+1}\mathbf{r}_{+,p}(u)K(u)\mathrm{d} u.
\end{align*}
lemLet Assumptions of Lemma (ref) hold separately on $\mathcal{X}_-$ and $\mathcal{X}_+$, then
\begin{align*}
\mathbb{V}\left[ \sqrt{\frac{n}{h}}\Big(\mathbf{e}_{1,+}-\mathbf{e}_{1,-}\Big)'\Big(f(\bar{x}+)\mathbf{S}_{+,p} +f(\bar{x}-)\mathbf{S}_{-,p}\Big)^{-1}\hat{\mathbf{L}} \right]
&= f(\bar{x}-) \mathbf{e}_{1,-}'\mathbf{S}_{-,p}^{-1}\boldsymbol{\Gamma}_{-,p}\mathbf{S}_{-,p}^{-1}\mathbf{e}_{1,-} \\
&\qquad+ f(\bar{x}+) \mathbf{e}_{1,+}'\mathbf{S}_{+,p}^{-1}\boldsymbol{\Gamma}_{+,p}\mathbf{S}_{+,p}^{-1}\mathbf{e}_{1,+} + O(h),
\end{align*}
where
\begin{align*}
\boldsymbol{\Gamma}_{-,p} &= \iint_{[-1,0]^2} (u\wedge v)\mathbf{r}_{-,p}(u)\mathbf{r}_{-,p}(v)'K(u)K(v)\ \mathrm{d} u\mathrm{d} v,\quad \boldsymbol{\Gamma}_{+,p} = \iint_{[0,1]^2} (u\wedge v)\mathbf{r}_{+,p}(u)\mathbf{r}_{+,p}(v)'K(u)K(v)\ \mathrm{d} u\mathrm{d} v.
\end{align*}
Note that the above gives the asymptotic variance of the difference $\hat{f}(\bar{x}+) - \hat{f}(\bar{x}-)$, and the variance takes an additive form. This is not surprising, since the two density estimates, $\hat{f}(\bar{x}+)$ and $\hat{f}(\bar{x}-)$, rely on distinctive subsamples, meaning that they are asymptotically independent.
Finally the order of $\hat{\mathbf{R}}$ can also be established.
lemLet Assumptions of Lemma (ref) hold separately on $\mathcal{X}_-$ and $\mathcal{X}_+$, then
\begin{align*}
\hat{\mathbf{R}} = O_{\mathbb{P}}\left( \sqrt{\frac{1}{n^2h}} \right).
\end{align*}
Now we state the main result concerning the manipulation testing. Let $\hat{\mathbf{S}}_{p,\bar{x}}$ and $\hat{\boldsymbol{\Gamma}}_{p,\bar{x}}$ be constructed as in Section (ref), and
align*[align* omitted — 231 chars of source]
coroAssume Assumptions (ref) and (ref) hold separately on $\mathcal{X}_-$ and $\mathcal{X}_+$ with $\alpha_x\geq p+1$ for some integer $p\geq 1$. Further, $n \cdot h^2\to \infty$ and $n\cdot h^{2p+1}\to 0$. Then under the null hypothesis $\mathsf{H}_0:f(\bar{x}+) = f(\bar{x}-)$,
\begin{align*}
T_p(h) &= \frac{\hat{f}_p(\bar{x}+) - \hat{f}_p(\bar{x}-) }{\sqrt{\frac{1}{nh}\hat{\mathcal{V}}_{p,1}(\bar{x})}} \rightsquigarrow \mathcal{N}(0,1).
\end{align*}
As a result, under the alternative hypothesis $\mathsf{H}_1:f(\bar{x}+) \neq f(\bar{x}-)$,
\begin{align*}
\lim_{n\to\infty}\mathbb{P}[|T_p(h)| \geq \Phi_{1-\alpha/2}] = 1.
\end{align*}
Here $\Phi_{1-\alpha/2}$ is the $(1-\alpha/2)$-quantile of the standard normal distribution.
remark[Separate estimation]
Recall that it is possible to implement our local polynomial estimator separately for the two subsamples, below and above the cutoff $\bar{x}$. Let $\hat{f}_{p,-}(\bar{x})$ and $\hat{f}_{p,+}(\bar{x})$ be the two density estimates, and $\hat{\mathcal{V}}_{p,1,-}(\bar{x})$ and $\hat{\mathcal{V}}_{p,1,+}(\bar{x})$ be the associated variance estimates. Then the test statistic is equivalently:
\begin{align*}
T_p(h) &= \frac{\frac{n_+}{n}\hat{f}_{p,+}(\bar{x}) - \frac{n_-}{n}\hat{f}_{p,-}(\bar{x}) }{\sqrt{\frac{1}{nh}\left(\frac{n_+}{n}\hat{\mathcal{V}}_{p,1,+}(\bar{x})+\frac{n_-}{n}\hat{\mathcal{V}}_{p,1,-}(\bar{x})\right)}}.
\end{align*}
$\parallel$
Restricted Model
In the previous subsection, we gave a test procedure on the discontinuity of the density by estimating on the two sides of the cutoff separately. This procedure is flexible and requires minimum assumptions. There are ways, however, to improve the power of the test when the densities are estimated with additional assumptions on the smoothness of the CDF
In a restricted model, the polynomial basis is re-defined as
align*[align* omitted — 154 chars of source]
and the estimator in the fully restricted model is
align*[align* omitted — 365 chars of source]
Again the notations (for example $\mathbf{X}$ and $\mathbf{X}_h$) are redefined similarly, with the scaling matrix $\mathbf{H}$ adjusted to ensure $\mathbf{H}^{-1}\mathbf{r}_p(u) = \mathbf{r}_p(h^{-1}u)$. Here $\hat F_p(\bar{x})$ is the estimated CDF and $\frac{1}{2}\hat F_p^{(2)}(\bar{x})$, $\cdots$, $\frac{1}{p!}\hat F_p^{(p)}(\bar{x})$ are the estimated higher order derivatives, which we assume are all continuous at $\bar{x}$, while $\hat f_{p}(\bar{x}-)$ and $\hat f_{p}(\bar{x}+)$ are the estimated densities on the two sides of $\bar{x}$. Therefore we call the above model restricted, since it only allows discontinuity of the first derivative of $F$ (i.e. the density) but not the other derivatives.
With the modification of the polynomial basis, all other matrices in the previous subsection are redefined similarly, and
align*[align* omitted — 173 chars of source]
where the subscripts indicate the corresponding derivatives to extract. Moreover
align*[align* omitted — 187 chars of source]
lemLet Assumptions of Lemma (ref) hold with the exception that $f$ may be discontinuous across $\bar{x}$, then
\begin{align*}
\nonumber\frac{1}{n}\mathbf{X}_h'\mathbf{K}_h\mathbf{X}_h &= \left\{f(\bar{x}-)\mathbf{S}_{-,p} +f(\bar{x}+)\mathbf{S}_{+,p} \right\} + O\left(h\right) + O_{\mathbb{P}}(1/\sqrt{nh}),
\end{align*}
where
\begin{align*}
\mathbf{S}_{-,p} = \int_{-1}^0 \mathbf{r}_{-,p}(u)\mathbf{r}_{-,p}(u)'K(u) \mathrm{d} u,\qquad \mathbf{S}_{+,p} = \int_0^1 \mathbf{r}_{+,p}(u)\mathbf{r}_{+,p}(u)'K(u) \mathrm{d} u.
\end{align*}
Again we decompose the estimator into four terms, $\hat{\mathbf{B}}_{\mathtt{LI}}$, $\hat{\mathbf{B}}_{\mathtt{S}}$, $\hat{\mathbf{L}}$ and $\hat{\mathbf{R}}$, which correspond to leave-in bias, smoothing bias, linear variance and quadratic variance, respectively.
lemLet Assumptions of Lemma (ref) hold with the exception that $f$ may be discontinuous across $\bar{x}$, then
\begin{align}
\hat{\mathbf{B}}_{\mathtt{S}}&= h^{p+1} \left\{ \frac{F^{(p+1)}(\bar{x}-)f(\bar{x}-)}{(p+1)!}\mathbf{c}_{-,p} + \frac{F^{(p+1)}(\bar{x}+)f(\bar{x}+)}{(p+1)!}\mathbf{c}_{+,p} \right\}+o_{\mathbb{P}}(h^{p+1}),\qquad \hat{\mathbf{B}}_{\mathtt{LI}}=O_{\mathbb{P}}\left(\frac{1}{n}\right),
\end{align}
where
\begin{align*}
\mathbf{c}_{-,p} = \int_{-1}^{0} u^{p+1}\mathbf{r}_{-,p}(u)K(u)\mathrm{d} u,\qquad \mathbf{c}_{+,p} = \int_{0}^{1} u^{p+1}\mathbf{r}_{+,p}(u)K(u)\mathrm{d} u.
\end{align*}
lemLet Assumptions of Lemma (ref) hold with the exception that $f$ may be discontinuous across $\bar{x}$, then
\begin{align*}
&\mathbb{V}\left[ \sqrt{\frac{n}{h}}\Big(\mathbf{e}_{1,+}-\mathbf{e}_{1,-}\Big)'\Big(f(\bar{x}+)\mathbf{S}_{+,p} +f(\bar{x}-)\mathbf{S}_{-,p}\Big)^{-1}\hat{\mathbf{L}} \right]\\
&\qquad = (\mathbf{e}_{1,+}-\mathbf{e}_{1,-})'(f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})^{-1}
(f(\bar{x}+)^3 \boldsymbol{\Gamma}_{+,p}\\
&\qquad\qquad+ f(\bar{x}-)^3 \boldsymbol{\Psi}\boldsymbol{\Gamma}_{+,p}\boldsymbol{\Psi})
(f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})^{-1}(\mathbf{e}_{1,+}-\mathbf{e}_{1,-}) + O(h),
\end{align*}
where
\begin{align*}
\boldsymbol{\Gamma}_{-,p} &= \iint_{[-1,0]^2} (u\wedge v)\mathbf{r}_{-,p}(u)\mathbf{r}_{-,p}(v)'K(u)K(v)\ \mathrm{d} u\mathrm{d} v,\quad \boldsymbol{\Gamma}_{+,p} = \iint_{[0,1]^2} (u\wedge v)\mathbf{r}_{+,p}(u)\mathbf{r}_{+,p}(v)'K(u)K(v)\ \mathrm{d} u\mathrm{d} v.
\end{align*}
and
\begin{align*}
\boldsymbol{\Psi} = \begin{bmatrix}
(-1)^0 & & & & & & \\
& & (-1)^1& & & & \\
& (-1)^1& & & & & \\
& & & (-1)^2& & & \\
& & & &(-1)^3& & \\
& & & & & \ddots & \\
& & & & & & (-1)^p
\end{bmatrix}.
\end{align*}
Again we can show that the quadratic part is negligible.
lemLet Assumptions of Lemma (ref) hold with the exception that $f$ may not be continuous across $\bar{x}$, then
\begin{align*}
\hat{\mathbf{R}} = O_{\mathbb{P}}\left( \sqrt{\frac{1}{n^2h}} \right).
\end{align*}
Now we state the main result concerning the manipulation testing. Let $\hat{\mathbf{S}}_{p,\bar{x}}$ and $\hat{\boldsymbol{\Gamma}}_{p,\bar{x}}$ be constructed as in Section (ref), and
align*[align* omitted — 231 chars of source]
coroAssume Assumptions (ref) and (ref) hold separately on $\mathcal{X}_-$ and $\mathcal{X}_+$ with $\alpha_x\geq p+1$ for some integer $p\geq 1$. Further, $n \cdot h^2\to \infty$ and $n\cdot h^{2p+1}\to 0$. Then under the null hypothesis $\mathsf{H}_0:f(\bar{x}+) = f(\bar{x}-)$,
\begin{align*}
T_p(h) &= \frac{\hat{f}_p(\bar{x}+) - \hat{f}_p(\bar{x}-) }{\sqrt{\frac{1}{nh}\hat{\mathcal{V}}_{p,1}(\bar{x})}} \rightsquigarrow \mathcal{N}(0,1).
\end{align*}
As a result, under the alternative hypothesis $\mathsf{H}_1:f(\bar{x}+) \neq f(\bar{x}-)$,
\begin{align*}
\lim_{n\to\infty}\mathbb{P}[|T_p(h)| \geq \Phi_{1-\alpha/2}] = 1.
\end{align*}
Here $\Phi_{1-\alpha/2}$ is the $(1-\alpha/2)$-quantile of the standard normal distribution.
Other Standard Error Estimators
The standard error $\hat{\sigma}_{p,v}(x)$ (see Theorem (ref)) is fully automatic and adapts to both interior and boundary regions. In this section we consider two other ways to construct a standard error.
Plug-in Standard Error
Take $v\geq 1$. Then the asymptotic variance of $\hat{F}^{(v)}_p(x)$ takes the following form:
align*[align* omitted — 151 chars of source]
One way of constructing estimate of the above quantity is to plug-in a consistent estimator of $f(x)$, which is simply the estimated density. Hence we can use
align*[align* omitted — 165 chars of source]
The next question is how $\mathbf{S}_{p,x}$ and $\boldsymbol{\Gamma}_{p,x}$ should be constructed. Note that they are related to the kernel, evaluation point $x$ and the bandwidth $h$, but not the data generating process. Therefore the three matrices can be constructed by either analytical integration or numerical method.
Jackknife-based Standard Error
The standard error $\hat{\sigma}_{p,v}(x)$ is obtained by inspecting the asymptotic linear representation. It is fully automatic and adapts to both interior and boundaries. In this part, we present another standard error which resembles $\hat{\sigma}_{p,v}(x)$, albeit with a different motivation.
Recall that $\hat{\boldsymbol{\beta}}_p(x)$ is essentially a second order U-statistic, and the following expansion is justified:
align*[align* omitted — 727 chars of source]
where the remainder represents leave-in bias. Note that the above could be written as a U-statistic, and to apply the Hoeffding decomposition, define
align*[align* omitted — 307 chars of source]
which is symmetric in its two arguments. Then
align*[align* omitted — 413 chars of source]
Here $\mathbf{U}_{1}(x_i) = \mathbb{E}\left[\left. \mathbf{U}(x_i,x_j)\right| x_i\right]$. The second line in the above display is the analogue of $\hat{\mathbf{L}}$, which contributes to the leading variance, and the third line is negligible. The new standard error, we call the jackknife-based standard error, is given by the following:
align*[align* omitted — 225 chars of source]
with
align*[align* omitted — 376 chars of source]
and
align*[align* omitted — 324 chars of source]
The name jackknife comes from the fact that we use leave-one-out “estimator” for $\mathbf{U}_{1}(x_i)$: with $x_i$ fixed,
align*[align* omitted — 123 chars of source]
Under the same conditions specified in Theorem (ref), one can show that the jackknife-based standard error is consistent.
Simulation Study
DGP 1: Truncated Normal Distribution
In this subsection, we conduct simulation study based on truncated normal distribution. To be more specific, the underlying distribution of $x_i$ is the standard normal distribution truncated below at $-0.8$. Therefore,
align*[align* omitted — 87 chars of source]
and zero otherwise. Equivalently, $x_i$ has Lebesgue density $\Phi^{(1)}(x)/(1-\Phi(-0.8))$ on $[-0.8,\infty]$.
In this simulation study, the target parameter is the density function evaluated at various points. Note that both the variance and the bias of our estimator depend on the evaluation point, and in particular, the magnitude of the bias depends on higher order derivatives of the distribution function.
enumerate• Evaluation point. We estimate the density at $x\in\{ -0.8,\ -0.5,\ 0.5,\ 1.5 \}$. Note that $-0.8$ is the boundary point, where classical density estimators such as the kernel density estimator has high bias. The point $-0.5$, given our bandwidth choice, is fairly close to the boundary, hence should be understood as in the lower boundary region. The two points $0.5$ and $1.5$ are interior, but the curvature of the normal density is quite different at those two points, and we expect to see the estimators having different bias behaviors.
• Polynomial order. We consider $p\in\{2,3\}$. For density estimation using our estimators, $p=2$ should be the default choice, since it corresponds to estimating conditional mean with local linear regression. Such choice is also recommended by Fan-Gijbels_1996_Book, according to which one should always choose $p-s=2-1=1$ to be an odd number. We include $p=3$ for completeness.
• Kernel function. For local polynomial regression, the choice of kernel function is usually not very important. We use the triangular kernel $k(u)=(1-|u|)\vee 0$.
• Sample size. The sample size used consists of $n\in\{1000,\ 2000\}$. For most empirical studies employing nonparametric density estimation, the sample size is well above 1000, hence $n=2000$ is more representative.
Overall, we have $4\times 2\times 2 = 16$ designs, and for each design, we conduct 5000 Monte Carlo repetitions.
We consider a grid of bandwidth choices, which correspond to multiples of the MSE-optimal bandwidth, ranging from $0.1h_{\mathtt{MSE}}$ to $2h_{\mathtt{MSE}}$. We also consider the estimated bandwidth. The MSE-optimal bandwidth, $h_{\mathtt{MSE}}$, is chosen by minimizing the asymptotic mean squared error, using the true underlying distribution.
For each design, we report the empirical bias of the estimator, $\mathbb{E}[\hat{f}_p(x)-f(x)]$, under bias. And empirical standard deviations, $\mathbb{V}^{1/2}[\hat{f}_p(x)]$, and empirical root-MSE, under sd and $\sqrt{\text{mse}}$, respectively. For the standard errors constructed from the variance estimators, we report their empirical average under mean, which should be compared to sd. We also report the empirical rejection rate of t-statistics at 5% nominal level, under size. The t-statistic is $(\hat{f}_p(x)-\mathbb{E} \hat{f}_p(x))/\text{se}$, which is exactly centered, hence rejection rate thereof is a measure of accuracy of normal approximation.
DGP 2: Exponential Distribution
In this subsection, we conduct simulation study based on exponential distribution. To be more specific, the underlying distribution of $x_i$ is $F(x)=1-e^{-x}$. Equivalently, $x_i$ has Lebesgue density $e^{-x}$ for $x\geq 0$.
In this simulation study, the target parameter is the density function evaluated at various points. Note that both the variance and the bias of our estimator depend on the evaluation point, and in particular, the magnitude of the bias depends on higher order derivatives of the distribution function.
enumerate• Evaluation point. We estimate the density at $x\in\{ 0,\ 1,\ 1.5 \}$. Note that $0$ is the boundary point, where classical density estimators such as the kernel density estimator has high bias. The two points $1$ and $1.5$ are interior.
• Polynomial order. We consider $p\in\{2,3\}$. For density estimation using our estimators, $p=2$ should be the default choice, since it corresponds to estimating conditional mean with local linear regression. Such choice is also recommended by Fan-Gijbels_1996_Book, according to which one should always choose $p-s=2-1=1$ to be an odd number. We include $p=3$ for completeness.
• Kernel function. For local polynomial regression, the choice of kernel function is usually not very important. We use the triangular kernel $k(u)=(1-|u|)\vee 0$.
• Sample size. The sample size used consists of $n\in\{1000,\ 2000\}$. For most empirical studies employing nonparametric density estimation, the sample size is well above 1000, hence $n=2000$ is more representative.
Overall, we have $3\times 2\times 2 = 12$ designs, and for each design, we conduct 5000 Monte Carlo repetitions.
We consider a grid of bandwidth choices, which correspond to multiples of the MSE-optimal bandwidth, ranging from $0.1h_{\mathtt{MSE}}$ to $2h_{\mathtt{MSE}}$. We also consider the estimated bandwidth. The MSE-optimal bandwidth, $h_{\mathtt{MSE}}$, is chosen by minimizing the asymptotic mean squared error, using the true underlying distribution.
For each design, we report the empirical bias of the estimator, $\mathbb{E}[\hat{f}_p(x)-f(x)]$, under bias. And empirical standard deviations, $\mathbb{V}^{1/2}[\hat{f}_p(x)]$, and empirical root-MSE, under sd and $\sqrt{\text{mse}}$, respectively. For the standard errors constructed from the variance estimators, we report their empirical average under mean, which should be compared to sd. We also report the empirical rejection rate of t-statistics at 5% nominal level, under size. The t-statistic is $(\hat{f}_p(x)-\mathbb{E} \hat{f}_p(x))/\text{se}$, which is exactly centered, hence rejection rate thereof is a measure of accuracy of normal approximation.
\addcontentsline{toc}{section}{References}
Proof
{
\singlespacing
Proof of Lemma (ref)
proofA generic element of the matrix $\frac{1}{n}\mathbf{X}_h'\mathbf{K}_h\mathbf{X}_h$ takes the form:
\begin{align*}
\frac{1}{n}\sum_i \frac{1}{h}\left(\frac{x_i-x}{h}\right)^s K\left(\frac{x_i-x}{h}\right),\qquad 0\leq s\leq 2p.
\end{align*}
Then we compute the expectation:
\begin{align*}
&\ \mathbb{E}\left[ \frac{1}{n}\sum_i \frac{1}{h}\left(\frac{x_i-x}{h}\right)^s K\left(\frac{x_i-x}{h}\right) \right] = \mathbb{E}\left[ \frac{1}{h}\left(\frac{x_i-x}{h}\right)^s K\left(\frac{x_i-x}{h}\right) \right]\\
&= \int_{{x}_{\mathtt{L}}}^{{x}_{\mathtt{U}}} \frac{1}{h}\left(\frac{u-x}{h}\right)^s K\left(\frac{u-x}{h}\right)f(u)\mathrm{d} u= \int_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} v^s K\left(v\right)f(x+vh) \mathrm{d} v
= \int_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} v^s K\left(v\right)f(x+vh) \mathrm{d} v,
\end{align*}
hence for $x$ in the interior,
\begin{align*}
&\ \mathbb{E}\left[ \frac{1}{n}\sum_i \frac{1}{h}\left(\frac{x_i-x}{h}\right)^s K\left(\frac{x_i-x}{h}\right) \right]= f(x)\int_{\mathbb{R}} \mathbf{r}_p(v)\mathbf{r}_p(v)^\prime K(v)\mathrm{d} v + o(1),
\end{align*}
and for $x={x}_{\mathtt{L}} + ch$ with $c\in[0, 1]$,
\begin{align*}
&\ \mathbb{E}\left[ \frac{1}{n}\sum_i \frac{1}{h}\left(\frac{x_i-x}{h}\right)^s K\left(\frac{x_i-x}{h}\right) \right]= f({x}_{\mathtt{L}})\int_{-c}^\infty \mathbf{r}_p(v)\mathbf{r}_p(v)^\prime K(v)\mathrm{d} v + o(1),
\end{align*}
and for $x={x}_{\mathtt{U}} - ch$ with $c\in[0, 1]$,
\begin{align*}
&\ \mathbb{E}\left[ \frac{1}{n}\sum_i \frac{1}{h}\left(\frac{x_i-x}{h}\right)^s K\left(\frac{x_i-x}{h}\right) \right]= f({x}_{\mathtt{U}})\int_{-\infty}^c \mathbf{r}_p(v)\mathbf{r}_p(v)^\prime K(v)\mathrm{d} v + o(1),
\end{align*}
provided that $F\in\mathcal{C}^1$.
The variance satisfies
\begin{align*}
&\ \mathbb{V}\left[ \frac{1}{n}\sum_i \frac{1}{h}\left(\frac{x_i-x}{h}\right)^s K\left(\frac{x_i-x}{h}\right) \right] = \frac{1}{n}\mathbb{V}\left[\frac{1}{h}\left(\frac{x_i-x}{h}\right)^s K\left(\frac{x_i-x}{h}\right)\right]\\
&\leq \frac{1}{n}\mathbb{E}\left[ \frac{1}{h^2}\left(\frac{x_i-x}{h}\right)^{2s} K\left(\frac{x_i-x}{h}\right)^2 \right]
= O\left(\frac{1}{nh}\right),
\end{align*}
provided that $F\in\mathcal{C}^1$.
Proof of Lemma (ref)
proofFirst consider the smoothing bias. The leading term can be easily obtain by taking expectation together with Taylor expansion of $F$ to power $p+1$. The variance of this term has order $n^{-1}h^{-1}h^{2p+2}$, which gives the residual estimate $o_{\mathbb{P}}(h^{p+1})$ since it is assumed that $nh\to \infty$.
Next for the leave-in bias, note that it has expectation of order $n^{-1}$, and variance of order $n^{-3}h^{-1}$, hence overall this term of order $O_{\mathbb{P}}(n^{-1})$.
Proof of Lemma (ref)
proofWe first compute the variance. Note that
\begin{align*}
& \int_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} \mathbf{r}_p\left(u\right)\Big(\tilde{F}(x+hu) - F(x+hu) \Big)K(u)f(x+hu)\mathrm{d} u\\
&\qquad = \frac{1}{n}\int_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} \mathbf{r}_p\left(u\right) \Big(\mathds{1}[x_i\leq x+hu] - F(x+hu) \Big)K(u)f(x+hu)\mathrm{d} u,
\end{align*}
and
\begin{align*}
&\ \mathbb{V}\left[\int_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} \mathbf{r}_p\left(u\right) \Big(\mathds{1}[x_i\leq x+hu] - F(x+hu) \Big)K(u)f(x+hu)\mathrm{d} u \right]\\
&= \iint_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} \mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v)f(x+hu)f(x+hv)\\
& \qquad\qquad\qquad\qquad\times \left[\int_{\mathbb{R}}\left(\mathds{1} [t\leq x+hu]-F(x+hu)\right)\left(\mathds{1} [t\leq x+hv]-F(x+hv)\right) f(t)\mathrm{d} t\right] \mathrm{d} u\mathrm{d} v\\
\tag{I}
&= \iint_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} \mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v)f(x+hu)f(x+hv) \Big(F(x+h(u\wedge v)) - F(x+hu)F(x+hv)\Big) \mathrm{d} u\mathrm{d} v.
\end{align*}
We first consider the interior case, where the above reduces to:
\begin{align*}
&\ (I)_{interior} \\
&= \iint_{\mathbb{R}} \mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v)f(x)^2 \Big(F(x) - F(x)^2\Big) \mathrm{d} u\mathrm{d} v\\
&\qquad+ h\iint_{\mathbb{R}} (u\wedge v)\mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v) f(x)^3 \mathrm{d} u\mathrm{d} v\\
&\qquad- h\iint_{\mathbb{R}} (u+ v)\mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v)f(x)^3 F(x)\mathrm{d} u\mathrm{d} v\\
&\qquad + h\iint_{\mathbb{R}} (u+ v)\mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v)f(x)F^{(2)}(x)\Big(F(x) - F(x)^2\Big) \mathrm{d} u\mathrm{d} v + o(h)\\
&= f(x)^2 \Big(F(x) - F(x)^2\Big)\mathbf{S}_{p,x}\mathbf{e}_0\mathbf{e}_0^\prime\mathbf{S}_{p,x} \\
&\qquad- h f(x)^3 F(x) \mathbf{S}_{p,x}(\mathbf{e}_1\mathbf{e}_0^\prime + \mathbf{e}_0\mathbf{e}_1^\prime)\mathbf{S}_{p,x}\\
&\qquad+ h f(x)F^{(2)}(x)\Big(F(x) - F(x)^2\Big) \mathbf{S}_{p,x}(\mathbf{e}_1\mathbf{e}_0^\prime + \mathbf{e}_0\mathbf{e}_1^\prime)\mathbf{S}_{p,x}\\
&\qquad+ h f(x)^3 \boldsymbol{\Gamma}_{p,x} + o(h).
\end{align*}
For $x={x}_{\mathtt{L}}+hc$ with $c\in[0, 1)$ in the lower boundary region,
\begin{align*}
&\ (I)_{lower boundary} \\
&= h\iint_{\mathbb{R}} (u\wedge v+c)\mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v) f({x}_{\mathtt{L}})^3 \mathrm{d} u\mathrm{d} v + o(h)= h f({x}_{\mathtt{L}})^3 \left(\boldsymbol{\Gamma}_{p,x}+c\mathbf{S}_{p,x}\mathbf{e}_0\mathbf{e}_0^\prime\mathbf{S}_{p,x}\right) + o(h).
\end{align*}
Finally, we have
\begin{align*}
&\ (I)_{upper boundary} \\
&= h\iint_{\mathbb{R}} (u\wedge v -c)\mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v) f({x}_{\mathtt{U}})^3 \mathrm{d} u\mathrm{d} v- h\iint_{\mathbb{R}} (u+ v - 2c)\mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v)f({x}_{\mathtt{U}})^3 \mathrm{d} u\mathrm{d} v + o(h)\\
&= h f({x}_{\mathtt{U}})^2 f({x}_{\mathtt{U}}) \left(\boldsymbol{\Gamma}_{p,x}+ c\mathbf{S}_{p,x}\mathbf{e}_0\mathbf{e}_0^\prime\mathbf{S}_{p,x}-\mathbf{S}_{p,x}(\mathbf{e}_1\mathbf{e}_0^\prime + \mathbf{e}_0\mathbf{e}_1^\prime)\mathbf{S}_{p,x}\right)+o(h).
\end{align*}
With the above results, it is easy to verify the variance formula, provided that we can show the asymptotic normality.
We first consider the interior case, and verify the Lindeberg condition on the fourth moment. Let $\boldsymbol{\alpha}\in\mathbb{R}^{p+1}$ be an arbitrary nonzero vector, then
\begin{align*}
&\ \sum_i \mathbb{E}\left( \frac{1}{\sqrt{n}} \boldsymbol{\alpha}^\prime\mathbf{N}_x(f(x)\mathbf{S}_{p,x})^{-1}\int_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} \mathbf{r}_p\left(u\right) \Big(\mathds{1}[x_i\leq x+hu] - F(x+hu) \Big)K(u)f(x+hu)\mathrm{d} u \right)^4\\
&= \frac{1}{n}\mathbb{E}\left( \boldsymbol{\alpha}^\prime\mathbf{N}_x (f(x)\mathbf{S}_{p,x})^{-1}\int_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} \mathbf{r}_p\left(u\right) \Big(\mathds{1}[x_i\leq x+hu] - F(x+hu) \Big)K(u)f(x+hu)\mathrm{d} u \right)^4\\
&= \frac{1}{n}\iiiint_{\mathcal{A}} \prod_{j=1,2,3,4} \Big(\boldsymbol{\alpha}^\prime\mathbf{N}_x(f(x)\mathbf{S}_{p,x})^{-1}\mathbf{r}_p\left(u_j\right) K(u_j)\Big) f(x+hu_j) \\
&\qquad \left[\int_{\mathbb{R}} \prod_{j=1,2,3,4}\Big(\mathds{1} [t\leq x+hu_j]-F(x+hu_j)\Big) f(t)\mathrm{d} t\right] \mathrm{d} u_1 \mathrm{d} u_2 \mathrm{d} u_3 \mathrm{d} u_4\\
&\leq \frac{C}{n} \cdot \iiiint_{\mathcal{A}} \prod_{j=1,2,3,4} \Big(\boldsymbol{\alpha}^\prime\mathbf{N}_x(f(x)\mathbf{S}_{p,x})^{-1}\mathbf{r}_p\left(u_j\right) K(u_j)\Big) f(x)\mathrm{d} u_1 \mathrm{d} u_2 \mathrm{d} u_3 \mathrm{d} u_4 + O\left(\frac{1}{nh}\right),
\end{align*}
where $\mathcal{A}=[\frac{{x}_{\mathtt{L}}-x}{h},\frac{{x}_{\mathtt{U}}-x}{h}]^4\subset \mathbb{R}^4$. The first term in the above display is asymptotically negligible, since it is takes the form $C\cdot(\boldsymbol{\alpha}^\prime\mathbf{N}_x\mathbf{e}_0)^4/n$ where the constant $C$ depends on the DGP, and is finite. The order of the next term is $1/(nh)$, which comes from multiplying $n^{-1}$, $h^{-2}$ (from the scaling matrix $\mathbf{N}_x$), and $h$ (from linearization), hence is also negligible.
Under the assumption that $nh\to \infty$, the Lindeberg condition is verified for interior case. The same logic applies to the boundary case, whose proof is easier than the interior case, since the leading term in the calculation is identically zero for $x$ in either the lower or upper boundary.
Proof of Lemma (ref)
proofFor $\hat{\mathbf{R}}$, we rewrite it as a second order degenerate U-statistic:
\begin{align*}
\hat{\mathbf{R}} &= \frac{1}{n^2}\sum_{i,j;i< j} \hat{\mathbf{U}}_{ij},
\end{align*}
where
\begin{align*}
\hat{\mathbf{U}}_{ij} &= \mathbf{r}_p\left(\frac{x_i-x}{h}\right)\Big(\mathds{1}[x_j\leq x_i] - F(x_i)\Big)K_h(x_i-x)+ \mathbf{r}_p\left(\frac{x_j-x}{h}\right)\Big(\mathds{1}[x_i\leq x_j] - F(x_j)\Big)K_h(x_j-x)\\
&- \mathbb{E}\left[ \mathbf{r}_p\left(\frac{x_i-x}{h}\right)\Big(\mathds{1}[x_j\leq x_i] - F(x_i)\Big)K_h(x_i-x)\Big| x_j \right]- \mathbb{E}\left[ \mathbf{r}_p\left(\frac{x_j-x}{h}\right)\Big(\mathds{1}[x_i\leq x_j] - F(x_j)\Big)K_h(x_j-x)\Big| x_i \right].
\end{align*}
To compute the leading term, it suffices to consider
\begin{align*}
&\ 2\mathbb{E}\left[ \mathbf{r}_p\left(\frac{x_i-x}{h}\right)\mathbf{r}_p\left(\frac{x_i-x}{h}\right)^\prime \left(\mathds{1} [x_j\leq x_i]-F(x_i)\right)^2K_{h}(x_i-x)^2 \right]\\
&= 2\mathbb{E}\left[ \mathbf{r}_p\left(\frac{x_i-x}{h}\right)\mathbf{r}_p\left(\frac{x_i-x}{h}\right)^\prime \Big( F(x_i) - F(x_i)^2 \Big)K_{h}(x_i-x)^2 \right]\\
&= \frac{2}{h}\int_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} \mathbf{r}_p\left(v\right)\mathbf{r}_p\left(v\right)^\prime \Big( F(x+hv) - F(x+hv)^2 \Big) K(v)^2 f(x+ hv)\mathrm{d} v\\
&= \frac{2}{h}\int_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} \mathbf{r}_p\left(v\right)\mathbf{r}_p\left(v\right)^\prime \Big( F(x) - F(x)^2 \Big) K(v)^2 f(x)\mathrm{d} v + O(1)\\
&=_{interior} \frac{2}{h}f(x)\left[ F(x) - F(x)^2 \right] \mathbf{T}_{p,x} + O(1),\\
&=_{boundary} O(1),
\end{align*}
which closes the proof.
Proof of Theorem (ref)
proofThis follows from previous lemmas.
Proof of Theorem (ref)
proofFirst we note that the second half of the theorem follows from the first half and the asymptotic normality result of Theorem (ref), hence it suffices to prove the first half, i.e. the consistency of $\hat{\mathcal{V}}_{p,v}(x)$.
The analysis of this estimator is quite involved, since it takes the form of a third order V-statistic. Moreover, since the empirical d.f. $\tilde{F}$ is involved in the formula, a full expansion leads to a fifth order V-statistic. However, some simple tricks will greatly simplify the problem.
We first split $\hat{\boldsymbol{\Gamma}}_{p,x}$ into four terms, respectively
\begin{align*}
\hat{\boldsymbol{\Sigma}}_{p,x,1} &= \frac{1}{n^3}\sum_{i,j,k} \mathbf{r}_p\left( \frac{x_j-x}{h} \right)\mathbf{r}_p\left( \frac{x_k-x}{h} \right)^\prime K_h(x_j-x)K_h(x_k-x) \Big(\mathds{1}[x_i\leq x_j] - F(x_j)\Big) \Big(\mathds{1}[x_i\leq x_k] - F(x_k)\Big)\\
\hat{\boldsymbol{\Sigma}}_{p,x,2} &= \frac{1}{n^3}\sum_{i,j,k} \mathbf{r}_p\left( \frac{x_j-x}{h} \right)\mathbf{r}_p\left( \frac{x_k-x}{h} \right)^\prime K_h(x_j-x)K_h(x_k-x)\Big( F(x_j)-\tilde{F}(x_j) \Big)\Big( \mathds{1}[x_i\leq x_k]-\tilde{F}(x_k) \Big)\\
\hat{\boldsymbol{\Sigma}}_{p,x,3} &= \frac{1}{n^3}\sum_{i,j,k} \mathbf{r}_p\left( \frac{x_j-x}{h} \right)\mathbf{r}_p\left( \frac{x_k-x}{h} \right)^\prime K_h(x_j-x)K_h(x_k-x)\Big( \mathds{1}[x_i\leq x_j]-\tilde{F}(x_j) \Big)\Big( F(x_k)-\tilde{F}(x_k) \Big)\\
\hat{\boldsymbol{\Sigma}}_{p,x,4} &= \frac{1}{n^3}\sum_{i,j,k} \mathbf{r}_p\left( \frac{x_j-x}{h} \right)\mathbf{r}_p\left( \frac{x_k-x}{h} \right)^\prime K_h(x_j-x)K_h(x_k-x) \Big(F(x_j)-\tilde{F}(x_j)\Big) \Big(F(x_k) - \tilde{F}(x_k)\Big).
\end{align*}
Leaving $\hat{\boldsymbol{\Sigma}}_{p,x,1}$ for a while, since it is the key component in this variance estimator. We first consider $\mathbf{N}_x\hat{\mathbf{S}}_{p,x}^{-1}\hat{\boldsymbol{\Sigma}}_{p,x,4}\hat{\mathbf{S}}_{p,x}^{-1}\mathbf{N}_x$. By the uniform consistency of the empirical d.f., it can be shown easily that
\begin{align*}
\mathbf{N}_x\hat{\mathbf{S}}_{p,x}^{-1}\hat{\boldsymbol{\Sigma}}_{p,x,4}\hat{\mathbf{S}}_{p,x}^{-1}\mathbf{N}_x &= O_{\mathbb{P}}\left( (nh)^{-1} \right).
\end{align*}
Note that the extra $h^{-1}$ comes from the scaling matrix $\mathbf{N}_x$, but not the kernel function $K_h$. Next we consider $\mathbf{N}_x\hat{\mathbf{S}}_{p,x}^{-1}\hat{\boldsymbol{\Sigma}}_{p,x,2}\hat{\mathbf{S}}_{p,x}^{-1}\mathbf{N}_x$, which takes the following form (up to the negligible smoothing bias):
\begin{align*}
\mathbf{N}_x\hat{\mathbf{S}}_{p,x}^{-1}\hat{\boldsymbol{\Sigma}}_{p,x,2}\hat{\mathbf{S}}_{p,x}^{-1}\mathbf{N}_x =& \mathbf{N}_x\mathbf{H}(\boldsymbol{\beta}_p(x) - \hat{\boldsymbol{\beta}}_p(x)) \left( \frac{1}{n^2}\sum_{i,k} \mathbf{r}_p\left( \frac{x_k-x}{h} \right)^\prime K_h(x_k-x)\Big( \mathds{1}[x_i\leq x_k]-\tilde{F}(x_k) \Big) \right)\hat{\mathbf{S}}_{p,x}^{-1}\mathbf{N}_x\\
=& O_{\mathbb{P}}((nh)^{-1/2}) = o_{\mathbb{P}}(1),
\end{align*}
where the last line uses the asymptotic normality of $\hat{\boldsymbol{\beta}}_p(x)$. For $\hat{\boldsymbol{\Sigma}}_{p,x,1}$, we make the observation that it is possible to ignore all “diagonal” terms, meaning that
\begin{align*}
\hat{\boldsymbol{\Sigma}}_{p,x,1} &= \frac{1}{n^3}\sum_{\substack{i,j,k \\ distinct}} \mathbf{r}_p\left( \frac{x_j-x}{h} \right)\mathbf{r}_p\left( \frac{x_k-x}{h} \right)^\prime K_h(x_j-x)K_h(x_k-x) \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}}(h),
\end{align*}
under the assumption that $nh^2\to \infty$. As a surrogate, define
\begin{align*}
\mathbf{U}_{i,j,k} &= \mathbf{r}_p\left( \frac{x_j-x}{h} \right)\mathbf{r}_p\left( \frac{x_k-x}{h} \right)^\prime K_h(x_j-x)K_h(x_k-x)\Big(\mathds{1}[x_i\leq x_j] - F(x_j)\Big) \Big(\mathds{1}[x_i\leq x_k] - F(x_k)\Big),
\end{align*}
which means
\begin{align*}
\hat{\boldsymbol{\Sigma}}_{p,x,1} &= \frac{1}{n^3}\sum_{\substack{i,j,k \\ distinct}} \mathbf{U}_{i,j,k}.
\end{align*}
The critical step is to further decompose the above into
\begin{align*}
\tag{I}\hat{\boldsymbol{\Sigma}}_{p,x,1} &= \frac{1}{n^3}\sum_{\substack{i,j,k \\ distinct}} \mathbb{E}[\mathbf{U}_{i,j,k}|x_i] \\
\tag{II}&+ \frac{1}{n^3}\sum_{\substack{i,j,k \\ distinct}} \Big(\mathbf{U}_{i,j,k} - \mathbb{E}[\mathbf{U}_{i,j,k}|x_i,x_j]\Big)\\
\tag{III}&+ \frac{1}{n^3}\sum_{\substack{i,j,k \\ distinct}} \Big(\mathbb{E}[\mathbf{U}_{i,j,k}|x_i,x_j] - \mathbb{E}[\mathbf{U}_{i,j,k}|x_i]\Big).
\end{align*}
We already investigated the properties of term (I) in Lemma (ref), hence it remains to show that both (II) and (III) are $o(h)$, hence does not affect the estimation of asymptotic variance. We consider (II) as an example, and the analysis of (III) is similar. Since (II) has zero expectation, we consider its variance (for simplicity treat $\mathbf{U}$ as a scaler):
\begin{align*}
\mathbb{V}[(II)] &= \mathbb{E}\left[ \frac{1}{n^6}\sum_{\substack{i,j,k \\ \text{distinct}}}\sum_{\substack{i',j',k' \\ \text{distinct}}} \Big(\mathbf{U}_{i,j,k} - \mathbb{E}[\mathbf{U}_{i,j,k}|x_i,x_j]\Big)\Big(\mathbf{U}_{i,j,k} - \mathbb{E}[\mathbf{U}_{i',j',k'}|x_{i'},x_{j'}]\Big) \right].
\end{align*}
The expectation will be zero if the six indices are all distinct. Similarly, when there are only two indices among the six are equal, the expectation will be zero \textit{unless} $k=k'$, hence
\begin{align*}
\mathbb{V}[\text{(II)}] &= \mathbb{E}\left[ \frac{1}{n^6}\sum_{\substack{i,j,k \\ \text{distinct}}}\sum_{\substack{i',j',k' \\ \text{distinct}}} \Big(\mathbf{U}_{i,j,k} - \mathbb{E}[\mathbf{U}_{i,j,k}|x_i,x_j]\Big)\Big(\mathbf{U}_{i,j,k} - \mathbb{E}[\mathbf{U}_{i',j',k'}|x_{i'},x_{j'}]\Big) \right]\\
&= \mathbb{E}\left[ \frac{1}{n^6}\sum_{\substack{i,j,k,i'j' \\ \text{distinct}}} \Big(\mathbf{U}_{i,j,k} - \mathbb{E}[\mathbf{U}_{i,j,k}|x_i,x_j]\Big)\Big(\mathbf{U}_{i,j,k} - \mathbb{E}[\mathbf{U}_{i',j',k}|x_{i'},x_{j'}]\Big) \right]\\
&+ \cdots,
\end{align*}
where $\cdots$ represent cases where more than two indices among the six are equal. We can easily compute the order from the above as
\begin{align*}
\mathbb{V}[\text{(II)}] &= O(n^{-1}) + O((nh)^{-2}),
\end{align*}
which shows that
\begin{align*}
\text{(II)} &= O_{\mathbb{P}}(n^{-1/2} + (nh)^{-1}) = o_{\mathbb{P}}(h),
\end{align*}
which closes the proof.
Proof of Lemma (ref)
proofWe rely on Lemma (ref) and (ref) (note that whether the weights are estimated is irrelevant here), hence will not repeat arguments already established there. Instead, extra care will be given to ensure the characterization of higher order bias.
Consider the case where with enough smoothness on $G$, then the bias is characterized by
\begin{align*}
&\ h^{-v}v!\mathbf{e}_v^\prime\left[f(x)\mathbf{S}_{p,x} + hF^{(2)}(x)\tilde{\mathbf{S}}_{p,x} + o(h) + O_{\mathbb{P}}(1/\sqrt{nh}) \right]^{-1} \\
&\ \qquad \left[ h^{p+1}\frac{F^{(p+1)}(x)}{(p+1)!}f(x)\mathbf{c}_{p,x}+h^{p+2}\left[\frac{F^{(p+2)}(x)}{(p+2)!}f(x)+\frac{F^{(p+1)}(x)}{(p+1)!}F^{(2)}(x)\right]\tilde{\mathbf{c}}_{p,x} + o(h^{p+2})\right]\\
&= h^{-v}v!\mathbf{e}_v^\prime\left[ \frac{1}{f(x)}\mathbf{S}_{p,x}^{-1}-h\frac{F^{(2)}(x)}{[f(x)]^2}\mathbf{S}_{p,x}^{-1}\tilde \mathbf{S}_{p,x} \mathbf{S}_{p,x}^{-1}+O_{\mathbb{P}}\left(1/\sqrt{nh}\right) \right]\\
&\ \qquad \left[ h^{p+1}\frac{F^{(p+1)}(x)}{(p+1)!}f(x)\mathbf{c}_{p,x}+h^{p+2}\left[\frac{F^{(p+2)}(x)}{(p+2)!}f(x)+\frac{F^{(p+1)}(x)}{(p+1)!}F^{(2)}(x)\right]\tilde{\mathbf{c}}_{p,x} + o(h^{p+2})\right]\{1+o_{\mathbb{P}}(1)\},
\end{align*}
which gives the desired result. Here $\tilde{\mathbf{S}}_{p,x} = \int_{\frac{{x}_{\mathtt{L}}-x}{h}}^{\frac{{x}_{\mathtt{U}}-x}{h}} u\mathbf{r}_p(u)\mathbf{r}_p(u)^\prime k(u)\mathrm{d} u$.
And for the last line to hold, one needs the extra condition $nh^3\to \infty$ so that $O_{\mathbb{P}}\left(1/\sqrt{nh}\right)=o_{\mathbb{P}}(h)$. See Fan-Gijbels_1996_Book (Theorem 3.1, pp. 62).
Proof of Lemma (ref)
proofThe proof resembles that of Lemma (ref), and is omitted here.
Proof of Theorem (ref)
proofThe proof splits into two cases. We sketch one of them. Assume either $x$ is boundary or $p-v$ is odd, the MSE-optimal bandwidth is asymptotically equivalent to the following:
\begin{align*}
\frac{\tilde{h}_{ p,v}(x)}{{h}_{ p,v}(x)}\to 1,\qquad \tilde{h}_{ p,v}(x) = \left(\frac{1}{n} \frac{(2v-1)f(x)\mathbf{e}_v^\prime\mathbf{S}_{p,x}^{-1}\boldsymbol{\Gamma}_{p,x}\mathbf{S}_{p,x}^{-1}\mathbf{e}_v}{(2p-2v+2)(\frac{F^{(p+1)}(x)}{(p+1)!}\mathbf{e}_v^\prime \mathbf{S}_{p,x}^{-1}\mathbf{c}_{p,x} )^2} \right)^{\frac{1}{2p+1}},
\end{align*}
which is obtained by optimizing MSE ignoring the higher order bias term. With consistency of the preliminary estimates, it can be shown that
\begin{align*}
\hat{h}_{ p,v}(x) &= \left(\frac{1}{n} \frac{(2v-1)\hat{\sigma}_{p,v}(x)^2n\ell^{2v-1}}{(2p-2v+2)(v!\frac{\hat{F}^{(p+1)}(x)}{(p+1)!}\mathbf{e}_v^\prime \mathbf{S}_{p,x}^{-1}\mathbf{c}_{p,x} )^2} \right)^{\frac{1}{2p+1}}\{1+o_{\mathbb{P}}(1)\}.
\end{align*}
Apply the consistency assumption of the preliminary estimates again, one can easily show that $\hat{h}_{ p,v}(x)$ is consistent both in rate and constant.
A similar argument can be made for the other case, and is omitted here.
Proof of Lemma (ref)
proofThis resembles the proof of Lemma (ref), and we only perform the mean computation. To start,
\begin{align*}
&\ \mathbb{E}\left[ \frac{1}{n}\mathbf{X}_h' \mathbf{K}_h \mathbf{X}_h \right] = \mathbb{E}\left[ \mathbf{r}_p\left(\frac{x_i-\bar{x}}{h}\right)\mathbf{r}_p\left(\frac{x_i-\bar{x}}{h}\right)'\frac{1}{h}K\left(\frac{x_i-\bar{x}}{h}\right) \right]\\
&= \mathbb{E}\left[\left. \mathbf{r}_p\left(\frac{x_i-\bar{x}}{h}\right)\mathbf{r}_p\left(\frac{x_i-\bar{x}}{h}\right)'\frac{1}{h}K\left(\frac{x_i-\bar{x}}{h}\right) \right| x_i<\bar{x}\right]F(\bar{x})\\
&\qquad + \mathbb{E}\left[\left. \mathbf{r}_p\left(\frac{x_i-\bar{x}}{h}\right)\mathbf{r}_p\left(\frac{x_i-\bar{x}}{h}\right)'\frac{1}{h}K\left(\frac{x_i-\bar{x}}{h}\right) \right| x_i\geq\bar{x}\right](1-F(\bar{x})).
\end{align*}
Then by Lemma (ref), the first term takes the form:
\begin{align*}
&\ \mathbb{E}\left[\left. \mathbf{r}_p\left(\frac{x_i-\bar{x}}{h}\right)\mathbf{r}_p\left(\frac{x_i-\bar{x}}{h}\right)'\frac{1}{h}K\left(\frac{x_i-\bar{x}}{h}\right) \right| x_i<\bar{x}\right]F(\bar{x})\\
&= f(\bar{x}- |x_i<\bar{x})F(\bar{x}) \int_{-1}^0 \mathbf{r}_{-,p}(u)\mathbf{r}_{-,p}(u)'K(u)\mathrm{d} u + O(h),
\end{align*}
where $f(\bar{x}- |x_i<\bar{x})$ is the one-sided density of $x_i$ at the cutoff, conditional on $x_i<\bar{x}$. Alternatively, we can simplify by the fact that $f(\bar{x} |x_i<\bar{x})F(\bar{x}) = f(\bar{x} -)$. Similarly, one has
\begin{align*}
&\ \mathbb{E}\left[\left. \mathbf{r}_p\left(\frac{x_i-\bar{x}}{h}\right)\mathbf{r}_p\left(\frac{x_i-\bar{x}}{h}\right)'\frac{1}{h}K\left(\frac{x_i-\bar{x}}{h}\right) \right| x_i\geq \bar{x}\right](1-F(\bar{x}))\\
&= f(\bar{x}+ |x_i\geq\bar{x})(1-F(\bar{x})) \int_{0}^1 \mathbf{r}_{+,p}(u)\mathbf{r}_{+,p}(u)'K(u)\mathrm{d} u + O(h),
\end{align*}
and that $f(\bar{x}+ |x_i\geq\bar{x})(1-F(\bar{x})) = f(\bar{x}+ )$. The rest of the proof follows standard variance calculation, and is not repeated here.
Proof of Lemma (ref)
proofThis follows from Lemma (ref) by splitting the bias calculation for the two subsamples, below and above the cutoff $\bar{x}$.
Proof of Lemma (ref)
proofTo start,
\begin{align*}
&\ \int_{-1}^{1} \mathbf{r}_p\left(u\right)\Big(\tilde{F}(\bar{x}+hu) - F(\bar{x}+hu) \Big)K(u)f(\bar{x}+hu)\mathrm{d} u = \frac{1}{n}\int_{-1}^{1} \mathbf{r}_p\left(u\right) \Big(\mathds{1}[x_i\leq \bar{x}+hu] - F(\bar{x}+hu) \Big)K(u)f(\bar{x}+hu)\mathrm{d} u,
\end{align*}
and
\begin{align*}
&\ \mathbb{V}\left[\int_{-1}^{1} \mathbf{r}_p\left(u\right) \Big(\mathds{1}[x_i\leq \bar{x}+hu] - F(\bar{x}+hu) \Big)K(u)f(\bar{x}+hu)\mathrm{d} u \right]\\
&= \iint_{-1}^{1} \mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv)\\
&\qquad\qquad\qquad\qquad\times\left[\int_{\mathbb{R}}\left(\mathds{1} [t\leq \bar{x}+hu]-F(\bar{x}+hu)\right)\left(\mathds{1} [t\leq \bar{x}+hv]-F(\bar{x}+hv)\right) f(t)\mathrm{d} t\right] \mathrm{d} u\mathrm{d} v\\
\tag{I}
&= \iint_{-1}^{1} \mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv) \Big(F(\bar{x}+h(u\wedge v)) - F(\bar{x}+hu)F(\bar{x}+hv)\Big) \mathrm{d} u\mathrm{d} v.
\end{align*}
Now we split the integral of (I) into four regions.
\begin{align*}
(u<0,v<0)\ (I) &= \iint_{-1}^{0} \mathbf{r}_{-,p}\left(u\right)\mathbf{r}_{-,p}\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv) \Big(F(\bar{x}+h(u\wedge v)) - F(\bar{x}+hu)F(\bar{x}+hv)\Big) \mathrm{d} u\mathrm{d} v\\
&= f(\bar{x}-)^2 \Big(F(\bar{x}) - F(\bar{x})^2\Big)\mathbf{S}_{-,p}\mathbf{e}_{0,-}\mathbf{e}_{0,-}^\prime\mathbf{S}_{-,p} \\
&\qquad- h f(\bar{x}-)^3 F(\bar{x}) \mathbf{S}_{-,p}(\mathbf{e}_{1,-}\mathbf{e}_{0,-}^\prime + \mathbf{e}_{0,-}\mathbf{e}_{1,-}^\prime)\mathbf{S}_{-,p}\\
&\qquad+ h f(\bar{x}-)F^{(2)}(\bar{x}-)\Big(F(\bar{x}) - F(\bar{x})^2\Big) \mathbf{S}_{-,p}(\mathbf{e}_{1,-}\mathbf{e}_{0,-}^\prime + \mathbf{e}_{0,-}\mathbf{e}_{1,-}^\prime)\mathbf{S}_{-,p}\\
&\qquad+ h f(\bar{x}-)^3 \boldsymbol{\Gamma}_{-,p} + O(h^2),
\end{align*}
and
\begin{align*}
(u\geq 0,v\geq 0)\ (I) &= \iint_{0}^{1} \mathbf{r}_{+,p}\left(u\right)\mathbf{r}_{+,p}\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv) \Big(F(\bar{x}+h(u\wedge v)) - F(\bar{x}+hu)F(\bar{x}+hv)\Big) \mathrm{d} u\mathrm{d} v\\
&= f(\bar{x}+)^2 \Big(F(\bar{x}) - F(\bar{x})^2\Big)\mathbf{S}_{+,p}\mathbf{e}_{0,+}\mathbf{e}_{0,+}^\prime\mathbf{S}_{+,p} \\
&\qquad- h f(\bar{x}+)^3 F(\bar{x}) \mathbf{S}_{+,p}(\mathbf{e}_{1,+}\mathbf{e}_{0,+}^\prime + \mathbf{e}_{0,+}\mathbf{e}_{1,+}^\prime)\mathbf{S}_{+,p}\\
&\qquad+ h f(\bar{x}+)F^{(2)}(\bar{x}+)\Big(F(\bar{x}) - F(\bar{x})^2\Big) \mathbf{S}_{+,p}(\mathbf{e}_{1,+}\mathbf{e}_{0,+}^\prime + \mathbf{e}_{0,+}\mathbf{e}_{1,+}^\prime)\mathbf{S}_{+,p}\\
&\qquad+ h f(\bar{x}+)^3 \boldsymbol{\Gamma}_{+,p} + O(h^2),
\end{align*}
and
\begin{align*}
(u<0,v\geq 0)\ (I) &= \iint_{[-1,0]\times[0,1]} \mathbf{r}_{-,p}\left(u\right)\mathbf{r}_{+,p}\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv) F(\bar{x}+hu) \Big(1- F(\bar{x}+hv)\Big) \mathrm{d} u\mathrm{d} v\\
&= \left[\int_{-1}^0 \mathbf{r}_{-,p}\left(u\right) K(u)f(\bar{x}+hu)F(\bar{x}+hu) \mathrm{d} u\right]\left[\int_{0}^1 \mathbf{r}_{+,p}\left(v\right)^\prime K(v)f(\bar{x}+hv) \Big(1- F(\bar{x}+hv)\Big) \mathrm{d} v\right]\\
&= \Big[ f(\bar{x}-)F(\bar{x}) \mathbf{S}_{-,p}\mathbf{e}_{0,-} + h\Big( f(\bar{x}-)^2 + F^{(2)}(\bar{x}-)F(\bar{x}) \Big) \mathbf{S}_{-,p}\mathbf{e}_{1,-} + O(h^2)\Big]\\
&\qquad \qquad \Big[ f(\bar{x}+)(1-F(\bar{x})) \mathbf{S}_{+,p}\mathbf{e}_{0,+} + h\Big( -f(\bar{x}+)^2 + F^{(2)}(\bar{x}+)(1-F(\bar{x})) \Big) \mathbf{S}_{+,p}\mathbf{e}_{1,+} + O(h^2)\Big]',
\end{align*}
and
\begin{align*}
(u\geq 0,v< 0)\ (I) &= \iint_{[0,1]\times[-1,0]} \mathbf{r}_{-,p}\left(u\right)\mathbf{r}_{+,p}\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv) F(\bar{x}+hv) \Big(1- F(\bar{x}+hu)\Big) \mathrm{d} u\mathrm{d} v\\
&= \left[\int_{0}^1 \mathbf{r}_{+,p}\left(u\right) K(u)f(\bar{x}+hu)(1-F(\bar{x}+hu)) \mathrm{d} u\right]\left[\int_{-1}^{0} \mathbf{r}_{-,p}\left(v\right)^\prime K(v)f(\bar{x}+hv) F(\bar{x}+hv)\Big) \mathrm{d} v\right]\\
&= \Big[ f(\bar{x}+)(1-F(\bar{x})) \mathbf{S}_{+,p}\mathbf{e}_{0,+} + h\Big( -f(\bar{x}+)^2 + F^{(2)}(\bar{x}+)(1-F(\bar{x})) \Big) \mathbf{S}_{+,p}\mathbf{e}_{1,+} + O(h^2)\Big]\\
&\qquad \qquad \Big[ f(\bar{x}-)F(\bar{x}-) \mathbf{S}_{-,p}\mathbf{e}_{0,-} + h\Big( f(\bar{x}-)^2 + F^{(2)}(\bar{x}-)F(\bar{x}) \Big) \mathbf{S}_{-,p}\mathbf{e}_{1,-} + O(h^2)\Big]'.
\end{align*}
Let $\mathbf{S}_{-,p}^{-1}$ and $\mathbf{S}_{+,p}^{-1}$ be the Moore–Penrose inverse of $\mathbf{S}_{-,p}$ and $\mathbf{S}_{+,p}$, respectively. Then
\begin{align*}
&\ \mathbb{V}\left[ (\mathbf{e}_{1,+}-\mathbf{e}_{1,-})'\sqrt{\frac{n}{h}}(f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})^{-1}\hat{\mathbf{L}} \right]\\
&= f(\bar{x}-) \mathbf{e}_{1,-}'\mathbf{S}_{-,p}^{-1}\boldsymbol{\Gamma}_{-,p}\mathbf{S}_{-,p}^{-1}\mathbf{e}_{1,-} + f(\bar{x}+) \mathbf{e}_{1,+}'\mathbf{S}_{+,p}^{-1}\boldsymbol{\Gamma}_{+,p}\mathbf{S}_{+,p}^{-1}\mathbf{e}_{1,+} + O(h).
\end{align*}
Proof of Lemma (ref)
proofThis follows from Lemma (ref) by splitting the bias calculation for the two subsamples, below and above the cutoff $\bar{x}$.
Proof of Corollary (ref)
proofThis follows from the previous lemmas and verifying the Lindeberg condition. See also the proof of Lemma (ref), Theorem (ref) and Theorem (ref).
Proof of Lemma (ref)
proofThis follows from Lemma (ref) by splitting the bias calculation for the two subsamples, below and above the cutoff $\bar{x}$. See also the proof of Lemma (ref).
Proof of Lemma (ref)
proofThis follows from Lemma (ref) by splitting the bias calculation for the two subsamples, below and above the cutoff $\bar{x}$.
Proof of Lemma (ref)
proofTo start,
\begin{align*}
&\ \int_{-1}^{1} \mathbf{r}_p\left(u\right)\Big(\tilde{F}(\bar{x}+hu) - F(\bar{x}+hu) \Big)K(u)f(\bar{x}+hu)\mathrm{d} u = \frac{1}{n}\int_{-1}^{1} \mathbf{r}_p\left(u\right) \Big(\mathds{1}[x_i\leq \bar{x}+hu] - F(\bar{x}+hu) \Big)K(u)f(\bar{x}+hu)\mathrm{d} u,
\end{align*}
and
\begin{align*}
&\ \mathbb{V}\left[\int_{-1}^{1} \mathbf{r}_p\left(u\right) \Big(\mathds{1}[x_i\leq \bar{x}+hu] - F(\bar{x}+hu) \Big)K(u)f(\bar{x}+hu)\mathrm{d} u \right]\\
&= \iint_{-1}^{1} \mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv)\\
&\qquad\qquad\qquad\qquad\times\left[\int_{\mathbb{R}}\left(\mathds{1} [t\leq \bar{x}+hu]-F(\bar{x}+hu)\right)\left(\mathds{1} [t\leq \bar{x}+hv]-F(\bar{x}+hv)\right) f(t)\mathrm{d} t\right] \mathrm{d} u\mathrm{d} v\\
\tag{I}&= \iint_{-1}^{1} \mathbf{r}_p\left(u\right)\mathbf{r}_p\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv) \Big(F(\bar{x}+h(u\wedge v)) - F(\bar{x}+hu)F(\bar{x}+hv)\Big) \mathrm{d} u\mathrm{d} v.
\end{align*}
Now we split the integral of (I) into four regions.
\begin{align*}
(u<0,v<0)\ (I) &= \iint_{-1}^{0} \mathbf{r}_{-,p}\left(u\right)\mathbf{r}_{-,p}\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv) \Big(F(\bar{x}+h(u\wedge v)) - F(\bar{x}+hu)F(\bar{x}+hv)\Big) \mathrm{d} u\mathrm{d} v\\
&= f(\bar{x}-)^2 \Big(F(\bar{x}) - F(\bar{x})^2\Big)\mathbf{S}_{-,p}\mathbf{e}_{0}\mathbf{e}_{0}^\prime\mathbf{S}_{-,p} \\
&\qquad- h f(\bar{x}-)^3 F(\bar{x}) \mathbf{S}_{-,p}(\mathbf{e}_{1,-}\mathbf{e}_{0}^\prime + \mathbf{e}_{0}\mathbf{e}_{1,-}^\prime)\mathbf{S}_{-,p}\\
&\qquad+ h f(\bar{x}-)F^{(2)}(\bar{x})\Big(F(\bar{x}) - F(\bar{x})^2\Big) \mathbf{S}_{-,p}(\mathbf{e}_{1,-}\mathbf{e}_{0}^\prime + \mathbf{e}_{0}\mathbf{e}_{1,-}^\prime)\mathbf{S}_{-,p}\\
&\qquad+ h f(\bar{x}-)^3 \boldsymbol{\Gamma}_{-,p} + O(h^2),
\end{align*}
and
\begin{align*}
(u\geq 0,v\geq 0)\ (I) &= \iint_{0}^{1} \mathbf{r}_{+,p}\left(u\right)\mathbf{r}_{+,p}\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv) \Big(F(\bar{x}+h(u\wedge v)) - F(\bar{x}+hu)F(\bar{x}+hv)\Big) \mathrm{d} u\mathrm{d} v\\
&= f(\bar{x}+)^2 \Big(F(\bar{x}) - F(\bar{x})^2\Big)\mathbf{S}_{+,p}\mathbf{e}_{0}\mathbf{e}_{0}^\prime\mathbf{S}_{+,p} \\
&\qquad- h f(\bar{x}+)^3 F(\bar{x}) \mathbf{S}_{+,p}(\mathbf{e}_{1,+}\mathbf{e}_{0}^\prime + \mathbf{e}_{0}\mathbf{e}_{1,+}^\prime)\mathbf{S}_{+,p}\\
&\qquad+ h f(\bar{x}+)F^{(2)}(\bar{x})\Big(F(\bar{x}) - F(\bar{x})^2\Big) \mathbf{S}_{+,p}(\mathbf{e}_{1,+}\mathbf{e}_{0}^\prime + \mathbf{e}_{0}\mathbf{e}_{1,+}^\prime)\mathbf{S}_{+,p}\\
&\qquad+ h f(\bar{x}+)^3 \boldsymbol{\Gamma}_{+,p} + O(h^2),
\end{align*}
and
\begin{align*}
(u<0,v\geq 0)\ (I) &= \iint_{[-1,0]\times[0,1]} \mathbf{r}_{-,p}\left(u\right)\mathbf{r}_{+,p}\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv) F(\bar{x}+hu) \Big(1- F(\bar{x}+hv)\Big) \mathrm{d} u\mathrm{d} v\\
&= \left[\int_{-1}^0 \mathbf{r}_{-,p}\left(u\right) K(u)f(\bar{x}+hu)F(\bar{x}+hu) \mathrm{d} u\right]\left[\int_{0}^1 \mathbf{r}_{+,p}\left(v\right)^\prime K(v)f(\bar{x}+hv) \Big(1- F(\bar{x}+hv)\Big) \mathrm{d} v\right]\\
&= \Big[ f(\bar{x}-)F(\bar{x}) \mathbf{S}_{-,p}\mathbf{e}_{0} + h\Big( f(\bar{x}-)^2 + F^{(2)}(\bar{x})F(\bar{x}) \Big) \mathbf{S}_{-,p}\mathbf{e}_{1,-} + O(h^2)\Big]\\
&\qquad \qquad \Big[ f(\bar{x}+)(1-F(\bar{x})) \mathbf{S}_{+,p}\mathbf{e}_{0} + h\Big( -f(\bar{x}+)^2 + F^{(2)}(\bar{x})(1-F(\bar{x})) \Big) \mathbf{S}_{+,p}\mathbf{e}_{1,+} + O(h^2)\Big]',
\end{align*}
and
\begin{align*}
(u\geq 0,v< 0)\ (I) &= \iint_{[0,1]\times[-1,0]} \mathbf{r}_{-,p}\left(u\right)\mathbf{r}_{+,p}\left(v\right)^\prime K(u)K(v)f(\bar{x}+hu)f(\bar{x}+hv) F(\bar{x}+hv) \Big(1- F(\bar{x}+hu)\Big) \mathrm{d} u\mathrm{d} v\\
&= \left[\int_{0}^1 \mathbf{r}_{+,p}\left(u\right) K(u)f(\bar{x}+hu)(1-F(\bar{x}+hu)) \mathrm{d} u\right]\left[\int_{-1}^{0} \mathbf{r}_{-,p}\left(v\right)^\prime K(v)f(\bar{x}+hv) F(\bar{x}+hv)\Big) \mathrm{d} v\right]\\
&= \Big[ f(\bar{x}+)(1-F(\bar{x})) \mathbf{S}_{+,p}\mathbf{e}_{0} + h\Big( -f(\bar{x}+)^2 + F^{(2)}(\bar{x})(1-F(\bar{x})) \Big) \mathbf{S}_{+,p}\mathbf{e}_{1,+} + O(h^2)\Big]\\
&\qquad \qquad \Big[ f(\bar{x}-)F(\bar{x}) \mathbf{S}_{-,p}\mathbf{e}_{0} + h\Big( f(\bar{x}-)^2 + F^{(2)}(\bar{x})F(\bar{x}) \Big) \mathbf{S}_{-,p}\mathbf{e}_{1,-} + O(h^2)\Big]'.
\end{align*}
By collecting terms, one has
\begin{align*}
(I) &= \Big(f(\bar{x}+)\mathbf{S}_{+,p} + f(\bar{x}+)\mathbf{S}_{-,p}\Big)\mathbf{e}_{0}\mathbf{e}_{0}'\Big(f(\bar{x}+)\mathbf{S}_{+,p} + f(\bar{x}+)\mathbf{S}_{-,p}\Big)'\\
& \qquad- h f(\bar{x}-) F(\bar{x}) f(\bar{x}-)\mathbf{S}_{-,p}\mathbf{e}_{1,-}\mathbf{e}_{0}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+ h \frac{F^{(2)}(\bar{x})}{f(\bar{x}-)}F(\bar{x})(1 - F(\bar{x})) f(\bar{x}-)\mathbf{S}_{-,p}\mathbf{e}_{1,-}\mathbf{e}_{0}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad- h f(\bar{x}-) F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,-}^\prime f(\bar{x}-)\mathbf{S}_{-,p}\\
& \qquad+h \frac{F^{(2)}(\bar{x})}{f(\bar{x}-)}F(\bar{x}) (1-F(\bar{x})) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,-}^\prime f(\bar{x}-)\mathbf{S}_{-,p}\\
& \qquad- h f(\bar{x}+) F(\bar{x}) f(\bar{x}+)\mathbf{S}_{+,p}\mathbf{e}_{1,+}\mathbf{e}_{0}^\prime(f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+h \frac{F^{(2)}(\bar{x})}{f(\bar{x}+)}(1-F(\bar{x}))F(\bar{x}) f(\bar{x}+)\mathbf{S}_{+,p}\mathbf{e}_{1,+}\mathbf{e}_{0}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad- h f(\bar{x}+) F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,+}^\prime f(\bar{x}+)\mathbf{S}_{+,p}\\
& \qquad+h \frac{F^{(2)}(\bar{x})}{f(\bar{x}+)}F(\bar{x})(1-F(\bar{x})) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,+}^\prime f(\bar{x}+)\mathbf{S}_{+,p}\\
& \qquad+hf(\bar{x}-) f(\bar{x}-)\mathbf{S}_{-,p}\mathbf{e}_{1,-}\mathbf{e}_{0}^\prime f(\bar{x}+)\mathbf{S}_{+,p}\\
& \qquad+h f(\bar{x}-) f(\bar{x}+)\mathbf{S}_{+,p}\mathbf{e}_{0}\mathbf{e}_{1,-}^\prime f(\bar{x}-)\mathbf{S}_{-,p}\\
& \qquad+h (f(\bar{x}+)^3 \boldsymbol{\Gamma}_{+,p}+ f(\bar{x}-)^3 \boldsymbol{\Gamma}_{-,p}).
\end{align*}
Next, we note that
\begin{align*}
\mathbf{S}_{+,p}\mathbf{e}_{1,-} = \mathbf{S}_{-,p}\mathbf{e}_{1,+} = \mathbf{0},
\end{align*}
which implies
\begin{align*}
(I) &= \Big(f(\bar{x}+)\mathbf{S}_{+,p} + f(\bar{x}+)\mathbf{S}_{-,p}\Big)\mathbf{e}_{0}\mathbf{e}_{0}'\Big(f(\bar{x}+)\mathbf{S}_{+,p} + f(\bar{x}+)\mathbf{S}_{-,p}\Big)'\\
& \qquad- h f(\bar{x}-) F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{1,-}\mathbf{e}_{0}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+ h \frac{F^{(2)}(\bar{x})}{f(\bar{x}-)}F(\bar{x})(1 - F(\bar{x})) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{1,-}\mathbf{e}_{0}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad- h f(\bar{x}-) F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,-}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+h \frac{F^{(2)}(\bar{x})}{f(\bar{x}-)}F(\bar{x}) (1-F(\bar{x})) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,-}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad- h f(\bar{x}+) F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{1,+}\mathbf{e}_{0}^\prime(f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+h \frac{F^{(2)}(\bar{x})}{f(\bar{x}+)}(1-F(\bar{x}))F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{1,+}\mathbf{e}_{0}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad- h f(\bar{x}+) F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,+}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+h \frac{F^{(2)}(\bar{x})}{f(\bar{x}+)}F(\bar{x})(1-F(\bar{x})) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,+}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+hf(\bar{x}-) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{1,-}\mathbf{e}_{0}^\prime f(\bar{x}+)\mathbf{S}_{+,p}\\
& \qquad+h f(\bar{x}-) f(\bar{x}+)\mathbf{S}_{+,p}\mathbf{e}_{0}\mathbf{e}_{1,-}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+h (f(\bar{x}+)^3 \boldsymbol{\Gamma}_{+,p}+ f(\bar{x}-)^3 \boldsymbol{\Gamma}_{-,p}).
\end{align*}
Next note that
\begin{align*}
\boldsymbol{\Gamma}_{-,p} &= \iint_{[-1,0]^2} (u\wedge v)\mathbf{r}_{-,p}(u)\mathbf{r}_{-,p}(v)'K(u)K(v) dudv\\
&= \iint_{[0,1]^2} ((-u)\wedge (-v))\mathbf{r}_{-,p}(-u)\mathbf{r}_{-,p}(-v)'K(u)K(v) dudv\\
&= \iint_{[0,1]^2} (u\wedge v - u - v)\boldsymbol{\Psi}\mathbf{r}_{+,p}(u)\mathbf{r}_{+,p}(v)'\boldsymbol{\Psi}K(u)K(v) dudv\\
&= \boldsymbol{\Psi}\boldsymbol{\Gamma}_{+,p}\boldsymbol{\Psi} - \boldsymbol{\Psi}\mathbf{S}_{+,p}\mathbf{e}_{1,+}\mathbf{e}_{0}'\mathbf{S}_{+,p}\boldsymbol{\Psi}-\boldsymbol{\Psi}\mathbf{S}_{+,p}\mathbf{e}_{0}\mathbf{e}_{1,+}'\mathbf{S}_{+,p}\boldsymbol{\Psi}\\
&= \boldsymbol{\Psi}\boldsymbol{\Gamma}_{+,p}\boldsymbol{\Psi} + \mathbf{S}_{-,p}\mathbf{e}_{1,-}\mathbf{e}_{0}'\mathbf{S}_{-,p}+\mathbf{S}_{-,p}\mathbf{e}_{0}\mathbf{e}_{1,-}'\mathbf{S}_{-,p},
\end{align*}
then
\begin{align*}
\text{(I)} &= \Big(f(\bar{x}+)\mathbf{S}_{+,p} + f(\bar{x}+)\mathbf{S}_{-,p}\Big)\mathbf{e}_{0}\mathbf{e}_{0}'\Big(f(\bar{x}+)\mathbf{S}_{+,p} + f(\bar{x}+)\mathbf{S}_{-,p}\Big)'\\
& \qquad- h f(\bar{x}-) F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{1,-}\mathbf{e}_{0}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+ h \frac{F^{(2)}(\bar{x})}{f(\bar{x}-)}F(\bar{x})(1 - F(\bar{x})) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{1,-}\mathbf{e}_{0}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad- h f(\bar{x}-) F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,-}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+h \frac{F^{(2)}(\bar{x})}{f(\bar{x}-)}F(\bar{x}) (1-F(\bar{x})) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,-}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad- h f(\bar{x}+) F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{1,+}\mathbf{e}_{0}^\prime(f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+h \frac{F^{(2)}(\bar{x})}{f(\bar{x}+)}(1-F(\bar{x}))F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{1,+}\mathbf{e}_{0}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad- h f(\bar{x}+) F(\bar{x}) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,+}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+h \frac{F^{(2)}(\bar{x})}{f(\bar{x}+)}F(\bar{x})(1-F(\bar{x})) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,+}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+hf(\bar{x}-) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{1,-}\mathbf{e}_{0}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+h f(\bar{x}-) (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\mathbf{e}_{0}\mathbf{e}_{1,-}^\prime (f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})\\
& \qquad+h (f(\bar{x}+)^3 \boldsymbol{\Gamma}_{+,p}+ f(\bar{x}-)^3 \boldsymbol{\Psi}\boldsymbol{\Gamma}_{+,p}\boldsymbol{\Psi}).
\end{align*}
Therefore,
\begin{align*}
&\ \mathbb{V}\left[ (\mathbf{e}_{1,+}-\mathbf{e}_{1,-})'\sqrt{\frac{n}{h}}(f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})^{-1}\hat{\mathbf{L}} \right]\\
&= (\mathbf{e}_{1,+}-\mathbf{e}_{1,-})'(f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})^{-1}
(f(\bar{x}+)^3 \boldsymbol{\Gamma}_{+,p}\\
&\qquad+ f(\bar{x}-)^3 \boldsymbol{\Psi}\boldsymbol{\Gamma}_{+,p}\boldsymbol{\Psi})
(f(\bar{x}+)\mathbf{S}_{+,p}+f(\bar{x}-)\mathbf{S}_{-,p})^{-1}(\mathbf{e}_{1,+}-\mathbf{e}_{1,-}) + O(h).
\end{align*}
Proof of Lemma (ref)
proofThis follows from Lemma (ref) by splitting the bias calculation for the two subsamples, below and above the cutoff $\bar{x}$.
Proof of Corollary (ref)
proofThis follows from the previous lemmas and verifying the Lindeberg condition. See also the proof of Lemma (ref), Theorem (ref) and Theorem (ref).
}
\addcontentsline{toc}{section}{Tables}
\newcounter{K}
\newcounter{M}
\newcounter{P}
\forloop{K}{1}{\value{K} < 2} {
\forloop{M}{1}{\value{M} < 5} {
\forloop{P}{1}{\value{P} < 3} {
\ifthenelse{\theK=1}{\ifthenelse{\theK=2}}
\ifthenelse{\theM=1}{\ifthenelse{\theM=2}{\ifthenelse{\theM=3}{\ifthenelse{\theM=4}}}}
\ifthenelse{\theP=1}
table[table omitted — 868 chars of source]
}}}
\forloop{K}{1}{\value{K} < 2} {
\forloop{M}{1}{\value{M} < 4} {
\forloop{P}{1}{\value{P} < 3} {
\ifthenelse{\theK=1}{\ifthenelse{\theK=2}}
\ifthenelse{\theM=1}{\ifthenelse{\theM=2}{\ifthenelse{\theM=3}{\ifthenelse{\theM=4}}}}
\ifthenelse{\theP=1}
table[table omitted — 859 chars of source]
}}}