Boundary Adaptive Local Polynomial Conditional Density Estimators

Supplementary material to “Boundary adaptive local polynomial conditional density estimators”

Setup

Let $\mathbf{x}_i\in\mathbb{R}^d$ and $y_i\in\mathbb{R}$ be continuously distributed random variables supported on $\mathcal{X}=[0,1]^d$ and $\mathcal{Y}=[0,1]$, respectively. We are interested in estimating the conditional distribution function and its derivatives:

where $\mu\in\mathbb{N}$, and $\boldsymbol{\nu}\in\mathbb{N}^d$ representing multi-indices. (In the main paper we only consider the estimation of conditional density and derivatives thereof, that is, we set $\boldsymbol{\nu}=0$ and $\mu = \vartheta + 1\geq 1$.)

To present our estimation strategy, we start from $\theta_{0,\boldsymbol{\nu}}$, the conditional distribution function and its derivatives with respect to the conditioning variable, and apply the local polynomial method:

where and $\mathbf{e}_{\boldsymbol{\nu}}^\intercal$ is a basis vector extracting the corresponding estimate. We can write the solution in closed form as

where

To estimate $\theta_{\mu, \boldsymbol{\nu}}$, we further smooth via local polynomials along the $y$-direction:

We can write the solution in closed-form as

where

While in the above we considered local polynomial regressions along both the $\mathbf{x}$- and $y$-directions, it is also possible to employ a local smoothing technique. To be precise, let $G$ be some function such that the following Lebesgue-Stieltjes integration is well-defined, then an alternative estimator can be constructed as

which has the solution

where

Notation

Limits are taken with respect to the sample size tending to infinity and the bandwidth shrinking to zero (i.e., $n\to\infty$ and $h\to 0$). For two positive sequences, $a_n\precsim b_n$ implies that $\limsup_{n\to\infty}|a_n/b_n|<\infty$. Similarly, $a_n\precsim_\mathbb{P} b_n$ means $|a_n/b_n|$ is asymptotically bounded in probability. We also adopt the small-o and big-O notation: $a_n = O_{\mathbb{P}}(b_n)$ is just $a_n\precsim_\mathbb{P} b_n$, and $a_n = o_{\mathbb{P}}(b_n)$ means $a_n/b_n$ converges to zero in probability. Constants that do not depend on the sample size or the bandwidth will be denoted by $\mathfrak{c}$, $\mathfrak{c}_1$, $\mathfrak{c}_2$, etc.

We introduce another notation, $O_{\mathtt{TC}}$, which not only provides an asymptotic order but also controls the tail probability. To be specific, $a_n = O_{\mathtt{TC}}(b_n)$ if for any $\mathfrak{c}_1 > 0$, there exists some $\mathfrak{c}_2$ such that \[\limsup_{n\to\infty} n^{\mathfrak{c}_1}\mathbb{P}\left[ a_n \geq \mathfrak{c}_2b_n \right] < \infty.\] Here the subscript, $\mathtt{TC}$, stands for “tail control.” Finally, let $\mathbf{X} = (\mathbf{x}_1^\intercal,\cdots,\mathbf{x}_n^\intercal)^\intercal$ and $\mathbf{Y} = (y_1,\cdots,y_n)^\intercal$ be the data matrices.

Overview

In this subsection we provide an overview of the main results. Underlying assumptions and precise statements of the lemmas and theorems will be given in later sections. First consider $\check{\theta}_{\mu,\boldsymbol{\nu}}(y|\mathbf{x})$, with a conditional expectation decomposition:

As we will show in Section (ref), the first term above consists of the centering of the estimator (i.e., the parameter of interest $\theta_{\mu,\boldsymbol{\nu}}(y|\mathbf{x})$) and the smoothing bias. The second term, on the other hand, gives the asymptotic representation of the estimator. To be precise, we have

As a result, we can focus on establishing properties of the the first term, which provides an equivalent kernel expression. Denote its variance by $\mathsf{V}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x})$. Then we show that the standardized process,

is approximately normally distributed both pointwise and uniformly for $y\in\mathcal{Y}$ and $\mathbf{x}\in\mathcal{X}$. To be even more precise, we establish a strong approximation result, meaning that there exists a copy $\bar{\mathbb{S}}^{\prime}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x})$, and a Gaussian process $\mathbb{G}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x})$ with the same covariance structure, such that \[ \sup_{ y \in \mathcal{Y}, \mathbf{x} \in \mathcal{X}}\left|\bar{\mathbb{S}}^{\prime}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x}) - \mathbb{G}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x})\right| = O_{\mathbb{P}}\left(\frac{\log^{d+1} (n)}{nh^{d+1}}\right)^{\frac{1}{2d+2}}. \] Together with a feasible variance-covariance estimator, the strong approximation result not only allows us to construct confidence bands for the target parameter and test shape restrictions, but also provides an explicit characterization of the coverage error probability for those procedures.

Inside the remainder term, $h^{\mathfrak{q}+1-|\boldsymbol{\nu}|} + h^{\mathfrak{p}+1-\mu}$ is the order of the leading smoothing bias, and $\log(n) \sqrt{\mathsf{V}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x})/(nh^d)}$ arises from the linearization step which replaces the random matrix $\hat{\mathbf{S}}_{\mathbf{x}}$ by its large-sample analogue $\mathbf{S}_{\mathbf{x}}$. It is worth mentioning that the order of the remainder term is uniformly valid for $y\in\mathcal{Y}$ and $\mathbf{x}\in\mathcal{X}$, which is why an extra logarithmic factor is present.

Now consider the other estimator, $\hat{\theta}_{\mu,\boldsymbol{\nu}}(y|\mathbf{x})$. While it is not possible to take a conditional expectation, we can still “center” the estimator with the conditional distribution function. That is,

As before, the first term captures the target parameter and the smoothing bias. The analysis of the second term is more involved. Besides the asymptotic linear representation term, it also consists of a leave-in bias term (since the same observation is used twice) and a second order U-statistic. We show that the following expansion holds uniformly for $y\in\mathcal{Y}$ and $\mathbf{x}\in\mathcal{X}$:

Here, the contribution of the U-statistic is represented by the order $\log (n)/\sqrt{n^2h^{d+2\mu+2|\boldsymbol{\nu}|+1}}$ in the remainder term. Interestingly, this term is negligible compared to the standard error, $\sqrt{\mathsf{V}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x})}$, provided that $\log (n)/(nh^2)\to \infty$.

The above demonstrates that important large-sample properties of the local regression based estimator, $\hat{\theta}_{{\mu, \boldsymbol{\nu}}}(y|\mathbf{x})$ --- such as pointwise and uniform normal approximation --- stem from the equivalent kernel representation. Here we note that the representation holds by setting $G = F_y$. In other words, $\hat{\theta}_{{\mu, \boldsymbol{\nu}}}(y|\mathbf{x})$ is first-order asymptotically equivalent to $\check{\theta}_{{\mu, \boldsymbol{\nu}}}(y|\mathbf{x})$ with the (infeasible) local smoothing using the marginal distribution $F_y$.

Assumptions

We make the following assumptions on the joint distribution, the kernel function, and the weighting $ G$.

Pointwise large-sample properties

We first present several uniform convergence results which will be used later to establish pointwise and uniform properties of our estimators.

We now follow the decomposition in Section (ref) and study the leading bias of our estimators.

For future reference, we define the order of the leading bias as

Next we study the leading variance of our estimator, defined as

For future reference, we will define

Now, we propose two estimators for the variance that are valid for all three cases of Lemma (ref), and hence will be useful for establishing a self-normalized distributional approximation later. Define

Note that $\hat{\mathsf{V}}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x})$ is simply the plug-in variance estimator for $\hat{\theta}_{{\mu, \boldsymbol{\nu}}}(y|\mathbf{x})$ and $\check{\mathsf{V}}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x})$ is the plug-in variance estimator for $\check{\theta}_{{\mu, \boldsymbol{\nu}}}(y|\mathbf{x})$. The next lemma provides pointwise convergence results for the two variance estimators.

Next, we study the large-sample distributional properties of the infeasible, standardized statistic

Note that this is equivalent to the scaled asymptotic linear representation of the estimator.

While the theorem focuses on asymptotic normality of the infeasible t-statistic, $\bar{\mathbb{S}}_{{\mu, \boldsymbol{\nu}}}^{\circ} (y, \mathbf{x})$, we show in the following remark that similar conclusions can be made for the t-statistics constructed with the estimators, $\hat{\theta}_{{\mu, \boldsymbol{\nu}}}(y|\mathbf{x})$ and $\check{\theta}_{{\mu, \boldsymbol{\nu}}}(y|\mathbf{x})$.

Uniform large-sample properties

To conduct statistical inference on the entire function $\theta_{{\mu, \boldsymbol{\nu}}}$, such as constructing confidence bands or testing shape restrictions, we need uniform distributional approximations to our estimators. In this section, we will consider large-sample properties of our estimator which hold uniformly on $\mathcal{Y} \times \mathcal{X} = [0,1]^{d+1}$. In the following remark, we demonstrate that the local sample size is uniformly large on the support $\mathcal{Y} \times \mathcal{X}$.

We now establish the uniform convergence rate of our estimator.

In the next lemma, we characterize the uniform convergence rate of the variance estimators introduced in the previous section.

Now, we introduce the Studentized processes for each of the estimators, $\hat{\theta}_{{\mu, \boldsymbol{\nu}}}$ and $\check{\theta}_{{\mu, \boldsymbol{\nu}}}$:

In the following lemma we study the error that arises from the Studentization of our estimators.

Our next goal is to establish a uniform normal approximation to the process $\bar{\mathbb{S}}_{{\mu, \boldsymbol{\nu}}}(y,\mathbf{x})$. See Appendix A.4 of the main paper for important properties of the equivalent kernel.

The Gaussian approximation in the above lemma is not feasible, as its covariance kernel depends on unknowns. To be more precise, the covariance kernel takes the form

where

We consider two estimators of the covariance kernel

and

Applications

Confidence bands

A natural corollary of Theorem (ref) is that one can employ critical values computed from $\check{\mathbb{G}}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x}) $ and $\hat{\mathbb{G}}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x})$ to construct confidence bands. To be very precise, define

Then level $(1-\alpha)$ confidence bands can be constructed as

whose coverage error is given in the following theorem.

Parametric specification testing

In applications, it is not uncommon to estimate conditional densities or higher-order derivatives by specifying a parametric family of distributions. While such parametric restrictions may provide reasonable approximations, it is still worthwhile to conduct specification testing. To be specific, assume the researcher postulates the following class

where $\mathsf{\Gamma}_{\mu,\boldsymbol{\nu}}$ is some compact parameter space. We abstract away from the specifics of the estimation technique, and assume that the researcher also picks some estimator (maximum likelihood, minimum distance, etc.) $\hat{\boldsymbol{\gamma}}$. Under fairly mild conditions, the estimator will converge in probability to some (possibly pseudo-true) parameter $\bar{\boldsymbol{\gamma}}$ in the parameter space $\mathsf{\Gamma}_{\mu,\boldsymbol{\nu}}$. As before, we will denote the true parameter as $\theta_{\mu,\boldsymbol{\nu}}(y|\mathbf{x})$, and consider the following competing hypotheses:

The test statistics we employ takes the following form

Testing shape restrictions

Now consider shape restrictions on the conditional density or its derivatives. Let $c(y,\mathbf{x})$ be a pre-specified function, and we study the following one-sided competing hypotheses.

The statistic we employ takes the form

and we will reject the null hypothesis if the test statistic exceeds a critical value.

Bandwidth selection

We assume throughout this section that $\mu >0$. Using the bias expression derived in Lemma (ref), and the leading variance is as characterized in Lemma (ref), we can derive precise expressions for bandwidth selection.

Pointwise asymptotic MSE minimization

Following from Fan-Gijbels_1996_Book, the pointwise MSE-optimal bandwidth is defined as a minimizer of the following optimization problem

The solution to this equation gives an MSE-optimal bandwidth that depends on (i) the order of the polynomials, (ii) the order of the derivative to be estimated, and (iii) the position of the evaluation point. \

Case 1: $\mathfrak{q}-|\boldsymbol{\nu}| = \mathfrak{p}-\mu$, odd

In this case, both the leading bias constants, $B_{(i),\mathfrak{q}+1}(y,\mathbf{x})$ and $B_{(ii),\mathfrak{p}+1}(y,\mathbf{x})$, are nonzero. Therefore, the MSE-optimal bandwidth is

Case 2: $\mathfrak{q}-|\boldsymbol{\nu}| = \mathfrak{p}-\mu$, even; either $\mathbf{x}$ or $y$ is at or near the boundary

In this case, at least one of the leading bias constants, $B_{(i),\mathfrak{q}+1}(y,\mathbf{x})$ and $B_{(ii),\mathfrak{p}+1}(y,\mathbf{x})$, is nonzero. Therefore, the MSE-optimal bandwidth is the same as in Case 1:

Case 3: $\mathfrak{q}-|\boldsymbol{\nu}| = \mathfrak{p}-\mu \neq 0$, even; both $\mathbf{x}$ and $y$ are interior

In this case, both leading bias constants are zero. Therefore, the MSE-optimal bandwidth will depend on higher-order bias terms:

Case 4: $\mathfrak{q}-|\boldsymbol{\nu}| = \mathfrak{p}-\mu = 0$, even; both $\mathbf{x}$ and $y$ are interior

As in Case 3, both leading bias constants are zero. The difference, however, is that the leading bias will involve an extra term:

Case 5: $\mathfrak{q}-|\boldsymbol{\nu}| < \mathfrak{p}-\mu$, $\mathfrak{q}-|\boldsymbol{\nu}|$ odd

In this case, the leading bias will involve only one term:

Case 6: $\mathfrak{q}-|\boldsymbol{\nu}| = \mathfrak{p}-\mu-1$, $\mathfrak{q}-|\boldsymbol{\nu}|$ even; $\mathbf{x}$ is interior

In this case, the leading bias will involve two terms:

Case 7: $\mathfrak{q}-|\boldsymbol{\nu}| < \mathfrak{p}-\mu-1$, $\mathfrak{q}-|\boldsymbol{\nu}|$ even; $\mathbf{x}$ is interior

In this case, the leading bias will involve only one term:

Case 8: $\mathfrak{q}-|\boldsymbol{\nu}| > \mathfrak{p}-\mu$, $\mathfrak{p}-\mu$ odd

In this case, the leading bias will involve only one term:

Case 9: $\mathfrak{q}-|\boldsymbol{\nu}|-1 = \mathfrak{p}-\mu$, $\mathfrak{q}-|\boldsymbol{\nu}|$ even; $y$ is interior

In this case, the leading bias will involve two terms:

Case 10: $\mathfrak{q}-|\boldsymbol{\nu}|-1 > \mathfrak{p}-\mu$, $\mathfrak{p}-\mu$ even; $y$ is interior

In this case, the leading bias will involve only one term:

Rule-of-thumb bandwidth selection

This section outlines the methodology that the companion R package, lpcde, uses to construct the rule-of-thumb bandwidth selection.

The rule-of-thumb estimation uses the following assumptions in order to compute the optimal bandwidth:

Using these assumptions, each of the terms in the formula given in Case 1 of Section (ref) are computed as follows:

Alternative variance estimators

V-statistic variance estimator

We propose here an alternative variance estimator that is quick to implement in practice. We start by first observing that the estimator $\hat{\theta}_{\mu, \boldsymbol{\nu}}(y|\mathbf{x})$ is a V-statistic. That is,

where,

Note that $a(\cdot)$ and $b(\cdot)$ are scalar functions that are non-zero only for data points that are within $h$ distance of the evaluation point. The second term in (ref) can now be symmetrized and treated as a U-statistic. Applying the Hoeffding decomposition to the symmetrized version of the second term and plugging it back into Equation (ref), we get

where

and

Dependence on polynomial orders is suppressed for notational simplicity. Since each of the terms in (ref) are orthogonal, the variance of the estimator can be expressed as the sum of the variance of each of the terms on the right hand side. Furthermore, we note that the first three terms and $W_{\mu, \boldsymbol{\nu}}(y, \mathbf{x})$ have higher-order variance. Thus, we only need to look at the variance of $L_{\mu, \boldsymbol{\nu}}(y, \mathbf{x})$.

where we know

We can expand and simplify this to get

Note that this expression is identical to the variance expression derived in the proof of Lemma (ref). This leads to a natural alternative jackknife covariance estimator,

where

In particular, note that if the two evaluation points are equivalent, we return the variance estimator,

Asymptotic variance estimator

Another alternative variance estimator is a sample version of the asymptotic variance derived in Lemma (ref). That is, each of the matrices in the formula are replaced with sample analogs. That is, \\ (i) $\mu=0$:

(ii) $\mu > 0$:

The covariance can be estimated using similar idea.

Proofs

Proof of Lemma (ref)

Part (i). See the proof of Lemma 1 in the main paper (i.e., Appendix A.3).

Part (ii). Next consider $\mathbf{e}_{\mu}^{\intercal}\mathbf{S}_y^{-1}\big(\bar{\mathbf{R}}_{y, \mathbf{x}} - \mathbb{E}[\bar{\mathbf{R}}_{y, \mathbf{x}} |\mathbf{X}]\big)$, which takes the form

It is straightforward to see that

for some $C'$ that holds uniformly for $y\in \mathcal{Y}$ and $\mathbf{x}\in\mathcal{X}$. We also have the following bound on the variance

Consider the first case above ($\mu=0$). By a discretization $\{ (y_\ell,\mathbf{x}_\ell): 1\leq \ell\leq M_n \}$ of $\mathcal{Y}\times\mathcal{X}$, we have the probabilistic bound due to Bernstein's inequality

provided that we set $\mathtt{r} = \sqrt{\log (n)/(nh^d)}$.$M_n$ is at most polynomial in $n$, and the error from discretization can be ignored. This concludes the proof for the $\mu = 0$ case.

For $\mu > 0$, we set $\mathtt{r} = \sqrt{\log (n)/(nh^{d-1})}$, and the probabilistic bound takes the form

This concludes the proof for the second case, where $\mu>0$.

Proof of Lemma (ref)

The conditional expectation of $\bar{\mathbf{R}}_{y, \mathbf{x}}$ in $\check{\theta}_{{\mu, \boldsymbol{\nu}}}$ is

To proceed, we employ a Taylor expansion of the conditional distribution function to order $s$:

Then, the conditional expectation can be simplified as

We note that

and

Therefore,

By Lemma (ref), the second term on the right-hand side satisfies

which means we can denote the leading bias as

For the second claim of this lemma, we again consider a Taylor expansion

Then

Proof of Lemma (ref)

Let $\mathbf{c}_1 = \mathbf{S}_y^{-1}\mathbf{e}_{\mu}$ and $\mathbf{c}_2 = \mathbf{S}_{\mathbf{x}}^{-1} \mathbf{e}_{\boldsymbol{\nu}}$.

We make a further expansion:

Note that the remainder term, $O(h^2)$, holds uniformly for $y\in\mathcal{Y}$ and $\mathbf{x}_i\in\mathcal{X}$ since the conditional distribution function is assumed to have bounded second derivative. Therefore, {

} To conclude the proof, we note that two scenarios can arise: $\mu=0$ and $\mu > 0$. In the second case,

The first case is more involved. If $\theta_{0,\mathbf{0}}(y|\mathbf{x})\neq 0,1$, then

If $\theta_{0,\mathbf{0}}(y|\mathbf{x})=0$ or $1$, then a further expansion is needed, which is why an extra $h$ will be present in the leading variance.

Proof of Lemma (ref)

Consistency of $\check{\mathsf{V}}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$. For the purposes of this proof, let $\mathbf{c}_1 = \mathbf{S}_{y}^{-1}\mathbf{e}_{\mu}$, $\hat{\mathbf{c}}_2 = \hat{\mathbf{S}}_{\mathbf{x}}^{-1}\mathbf{e}_{\boldsymbol{\nu}}$, and $\mathbf{c}_2 = \mathbf{S}_{\mathbf{x}}^{-1}\mathbf{e}_{\boldsymbol{\nu}}$. To start, consider

First consider term (III). With the uniform convergence result for the estimated conditional distribution function, it is clear that {

}

Now we study term (I), which is clearly unbiased for $\mathsf{V}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$. Therefore, we compute its variance.

With iterative expectation (by conditioning on $\mathbf{x}_i$), the above further reduces to

In other words,

Finally, we consider (II). Using the Cauchy-Schwartz inequality, we have

As a result,

To conclude the proof for $\check{\mathsf{V}}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$, we note that replacing $\hat{\mathbf{c}}_2$ by $\mathbf{c}_2$ only leads to an additional multiplicative factor $1 + O_{\mathbb{P}}(1/\sqrt{nh^d})$. See Lemma (ref).

Consistency of $\hat{\mathsf{V}}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$. For the purposes of this proof, let $\hat{\mathbf{c}}_1 = \hat{\mathbf{S}}_{y}^{-1}\mathbf{e}_{\mu}$, $\mathbf{c}_1 = \mathbf{S}_{y}^{-1}\mathbf{e}_{\mu}$, $\hat{\mathbf{c}}_2 = \hat{\mathbf{S}}_{\mathbf{x}}^{-1}\mathbf{e}_{\boldsymbol{\nu}}$, and $\mathbf{c}_2 = \mathbf{S}_{\mathbf{x}}^{-1}\mathbf{e}_{\boldsymbol{\nu}}$. We first consider the following decomposition

By the uniform convergence rate of the estimated conditional distribution function, we have

Next we consider (II). Using the Cauchy-Schwartz inequality, we have

Finally, consider term (I), which has the expansion

Then,

Using similar techniques, one can show that

To streamline the remaining derivation, define

Then

We have studied the term (I.1.3) in the proof for $\check{\mathsf{V}}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$. In particular,

Term (I.1.1) is a mean zero third order U-statistic. Consider its variance

The above expectation is non-zero only in three scenarios: $(j=j',k=k',i\neq i')$, $(j=j',k=k',i=i')$ or $(j=i',k=k',i=j')$. Therefore,

Finally consider (I.1.2), which has a mean of zero. Its variance is

In addition, an extra $h$ factor emerges if $\mu > 0$, or if $\theta_{0, \mathbf{0}}= 0$ or $1$. As a result,

To conclude the proof for $\hat{\mathsf{V}}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$, we note that replacing $\hat{\mathbf{c}}_1$ by $\mathbf{c}_1$ and $\hat{\mathbf{c}}_2$ by $\mathbf{c}_2$ only leads to an additional multiplicative factor $1 + O_{\mathbb{P}}(1/\sqrt{nh^d})$. See Lemma (ref).

Proof of Theorem (ref)

We will write

Define $\mathbf{c}_1 = \mathbf{S}_{y}^{-1}\mathbf{e}_{\mu}$ and $\mathbf{c}_2 = \mathbf{S}_{\mathbf{x}}^{-1}\mathbf{e}_{\boldsymbol{\nu}}$.

To apply the Berry-Esseen theorem, we first compute the third moment

The leading term in the above is simply

Note that the above will be exactly zero in cases (ii) and (iii) of Lemma (ref).

Omitted details of Remark (ref)

Approximation and coverage error of $\check{\mathbb{S}}_{{\mu, \boldsymbol{\nu}}}^\circ(y, \mathbf{x})$. To start, {

} By allowing the constant $\mathfrak{c}_1$ to take possibly different values in each term, we have

where the conclusions follow from the uniform rates established in Lemma (ref) and the variance calculations in Lemma (ref). Next, we consider the normal approximation error. Note that

which means

where $\mathtt{r_{BE}}$ is defined in Theorem (ref).

\noindentApproximation and coverage error of $\hat{\mathbb{S}}_{{\mu, \boldsymbol{\nu}}}^\circ(y, \mathbf{x})$. To begin with, we decompose the double sum into

where we set $G = F_y$. Term (I) represents the leave-in bias, and it is straightforward to show that

for some constants $\mathfrak{c}_{1}$, $\mathfrak{c}_{2}$, and $\mathfrak{c}_{3}$. See Lemma (ref) for the proof strategy.

Term (II) is a degenerate U-statistic. Define

where $\mathbf{c}_1$ and $\mathbf{c}_2$ are arbitrary (fixed) vectors of conformable dimensions. Then we apply Equation (3.5) of Gine-Latala-Zinn_2000_Ustat (Lemmas 7 and 8 in the main paper), which gives (the value of $C'$ may change for each line)

As a result,

for some constants $\mathfrak{c}_{1}$, $\mathfrak{c}_{2}$, and $\mathfrak{c}_{3}$.

We now collect the pieces. The difference between $\hat{\mathbb{S}}_{{\mu, \boldsymbol{\nu}}}^\circ(y, \mathbf{x})$ and $\check{\mathbb{S}}_{{\mu, \boldsymbol{\nu}}}^\circ(y, \mathbf{x})$ is

and the conclusion follows from Lemmas (ref) and (ref).

Omitted details of Remark (ref)

To show this result, we first partition the support $\mathcal{Y} \times \mathcal{X}$ into cubes with edge length $\mathfrak{c}_3h$, where the constant $\mathfrak{c}_3$ is chosen so that, for any $(y,\mathbf{x})$ in $\mathcal{Y} \times \mathcal{X}$, at least one of the cubes will be contained in the ball $\{ y': |y'-y|\leq \mathfrak{c}_1h\}\times \{ \mathbf{x}': |\mathbf{x}'-\mathbf{x}|\leq \mathfrak{c}_1h\}$. The number of cubes in this partition is $\lceil1/(\mathfrak{c}_3h)^{d+1}\rceil$. Then the conclusion follows from Lemma (ref).

Proof of Lemma (ref). For simplicity let $c_n = (1-\delta_n) \frac{n}{J_n}$. We first employ the union bound

Note that $z_j\sim \mathrm{Binomial}(n;\frac{1}{J_n})$, and therefore

Then we have

Therefore, the above will vanish faster than any polynomial of $n$ provided that $\delta_n^2\pi_n\to\infty$.

Proof of Lemma (ref)

Part (i) Convergence of $\check{\theta}_{\mu,\boldsymbol{\nu}} - \theta_{\mu,\boldsymbol{\nu}}$. Recall that we have the following decomposition of our estimator

(I) is simply the conditional bias, whose order is given in Lemma (ref). The convergence rate of (II) can be easily deduced from that of $\mathbf{e}_{\mu}^{\intercal}\mathbf{S}_y^{-1}\left(\bar{\mathbf{R}}_{y, \mathbf{x}} - \mathbb{E} \left[\bar{\mathbf{R}}_{y, \mathbf{x}}|\mathbf{X}\right]\right)$ in Lemma (ref). Finally, it should be clear that (III) is negligible relative to (II).

Part (ii) Convergence of $\hat{\theta}_{\mu,\boldsymbol{\nu}} - \theta_{\mu,\boldsymbol{\nu}}$. This part follows from Remark (ref).

Proof of Lemma (ref)

Uniform consistency of $\check{\mathsf{V}}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$. For the purposes of this proof, let $\mathbf{c}_1 = \mathbf{S}_{y}^{-1}\mathbf{e}_{\mu}$, $\hat{\mathbf{c}}_2 = \hat{\mathbf{S}}_{\mathbf{x}}^{-1}\mathbf{e}_{\boldsymbol{\nu}}$, and $\mathbf{c}_2 = \mathbf{S}_{\mathbf{x}}^{-1}\mathbf{e}_{\boldsymbol{\nu}}$. To start, consider

First consider term (I). Clearly this term is unbiased for $\mathsf{V}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$. In the proof of Lemma (ref), we showed that

Also note that

In the above, the constants, $C_1$ and $C_2$, can be chosen to be independent of the evaluation point, the sample size, and the bandwidth. Then by a proper discretization of $\mathcal{Y}\times \mathcal{X}$, and applying the union bound and Bernstein's inequality, one has

for some constants $\mathfrak{c}_1$, $\mathfrak{c}_2$, and $\mathfrak{c}_3$. In addition, $\mathfrak{c}_3$ can be made arbitrarily large by appropriate choices of $\mathfrak{c}_1$. See the proof of Lemma (ref) for an example of this proof strategy.

Next consider term (III). With the uniform convergence result for the estimated conditional distribution function, it is clear that

Finally, we consider (II). Using the Cauchy-Schwartz inequality, we have

As a result,

To conclude the proof for $\check{\mathsf{V}}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$, we note that replacing $\hat{\mathbf{c}}_2$ by $\mathbf{c}_2$ only leads to an additional multiplicative factor $1 + O_{\mathbb{P}}(\sqrt{\log (n)/(nh^d)})$. See Lemma (ref).

\noindentUniform consistency of $\hat{\mathsf{V}}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$. For the purposes of this proof, let $\hat{\mathbf{c}}_1 = \hat{\mathbf{S}}_{y}^{-1}\mathbf{e}_{\mu}$, $\mathbf{c}_1 = \mathbf{S}_{y}^{-1}\mathbf{e}_{\mu}$, $\hat{\mathbf{c}}_2 = \hat{\mathbf{S}}_{\mathbf{x}}^{-1}\mathbf{e}_{\boldsymbol{\nu}}$, and $\mathbf{c}_2 = \mathbf{S}_{\mathbf{x}}^{-1}\mathbf{e}_{\boldsymbol{\nu}}$. We consider the same decomposition used in the proof of Lemma (ref):

By the uniform convergence rate for the estimated conditional distribution function, we have

Employing the Cauchy-Schwartz inequality gives

As a result, a probabilistic order for term (II) follows that of terms (I) and (III).

Finally, consider term (I), which has the expansion

Then,

which means

Using similar techniques, one can show that

To streamline the remaining derivation, define

Then

By employing the same techniques in the proof for $\check{\mathsf{V}}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$, we have that

Term (I.1.2) admits the following decomposition:

Using the same techniques of Lemmas (ref) and (ref), we have

Term (I.1.2.2) is a degenerate second order U-statistic. We adopt Equation (3.5) of Gine-Latala-Zinn_2000_Ustat (Lemmas 7 and 8 in the main paper), which implies (see Remark (ref) and its proof for an example)

To handle term (I.1.1), first consider the quantity $\phi_{j,i} - \phi_i$, which takes the form

Then it is straightforward to show that

As a result,

To conclude the proof for $\hat{\mathsf{V}}_{\mu,\boldsymbol{\nu}}(y,\mathbf{x})$, we note that replacing $\hat{\mathbf{c}}_1$ by $\mathbf{c}_1$ and $\hat{\mathbf{c}}_2$ by $\mathbf{c}_2$ only leads to an additional multiplicative factor $1 + O_{\mathbb{P}}(\sqrt{\log (n)/nh^d})$. See Lemma (ref).

Proof of Lemma (ref)

First consider $\check{\mathbb{T}}_{{\mu, \boldsymbol{\nu}}}^\circ(y, \mathbf{x})$. The difference between $\check{\mathbb{T}}_{{\mu, \boldsymbol{\nu}}}^\circ(y, \mathbf{x})$ and $\check{\mathbb{S}}_{{\mu, \boldsymbol{\nu}}}^\circ(y, \mathbf{x})$ is

From Lemma (ref), we have

To close the proof, it is straightforward to verify that

which follows from the uniform convergence rate in Lemma (ref). The same technique applies to the analysis of $\hat{\mathbb{T}}_{{\mu, \boldsymbol{\nu}}}^\circ(y, \mathbf{x})-\hat{\mathbb{S}}_{{\mu, \boldsymbol{\nu}}}^\circ(y, \mathbf{x})$.

Proof of Theorem (ref)

See the proof of Theorem 2 in the main paper (i.e., Appendix A.6).

Proof of Lemma (ref)

Consider $\check{\rho}_{{\mu, \boldsymbol{\nu}}}(y,\mathbf{x},y',\mathbf{x}')$. Note that we can decompose the difference into

The probabilistic order of the second term is given in Lemma (ref).

Using similar techniques as in the proof of Lemma (ref) or (ref), it is also straightforward to verify that term (I) has the same order. That is,

Proof of Lemma (ref)

Consider an $\varepsilon$ discretization of $\mathcal{Y}\times \mathcal{X}$, which is denoted by $\mathcal{A}_{\varepsilon} = \{(y_\ell,\mathbf{x}_{\ell}^\intercal)^\intercal:\ 1\leq \ell \leq L\}$. Then one can define two Gaussian vectors, $\mathbf{z},\check{\mathbf{z}}\in\mathbb{R}^L$, such that

Then we apply the Gaussian comparison result in Corollary 5.1 of chernozhukov2022central (Lemma 11 in the main paper) and the error rate in Lemma (ref), which lead to

Since $\varepsilon$ only enters the above error bound logarithmically, one can choose $\varepsilon = n^{-c}$ for some $c$ large enough, so that the error that arises from discretization becomes negligible. The same applies to $\hat{\mathbb{G}}_{{\mu, \boldsymbol{\nu}}}(y_{\ell},\mathbf{x}_{\ell})$.

Proof of Theorem (ref)

First consider $\check{\mathbb{T}}_{{\mu, \boldsymbol{\nu}}}^\circ(y, \mathbf{x})$. Since

then with Lemma (ref),

In the above, we also used the fact that the difference $\check{\mathbb{S}}_{{\mu, \boldsymbol{\nu}}}^\circ(y, \mathbf{x})-\bar{\mathbb{S}}_{{\mu, \boldsymbol{\nu}}}(y, \mathbf{x})$ is negligible compared to $\mathtt{r_{SE}}$ (see Remark (ref)).

By applying Lemma (ref),

Finally, we apply the Gaussian comparison result in Lemma (ref), which implies that

Finally, due to Theorem 2.1 of Chernozhukov-Chetverikov-Kato_2014b_AoS (Lemma 12 in the main paper), we have

Proof of Theorem (ref)

Note that $\theta_{\mu,\boldsymbol{\nu}}(y|\mathbf{x})$ falls into the confidence band $\check{\mathcal{C}}_{{\mu, \boldsymbol{\nu}}}(1-\alpha)$ if and only if

A sufficient condition would then be

The conclusion then follows from Theorem (ref) and the bias calculation in Lemma (ref). The same analysis applies to $\hat{\mathcal{C}}_{{\mu, \boldsymbol{\nu}}}(1-\alpha)$.

Proof of Theorem (ref)

To start, we decompose the test statistic into

Then by the leading bias order in Lemma (ref) and the leading variance order in Lemma (ref), we have that

Similarly, under the null hypothesis,

Then we have the following error bound

As a result,

The same strategy can be employed to establish results for $\hat{\mathbb{T}}_{\mathtt{PS}}(y, \mathbf{x})$.

Proof of Theorem (ref)

The conclusion follows directly from Theorem (ref).