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.
77,210 characters · 11 sections · 40 citation commands
Coverage Error Optimal Confidence Intervals for Local Polynomial Regression
\addtocontents{toc}{\setcounter{tocdepth}{2}}
\addtocontents{toc}{\setcounter{tocdepth}{0}}
We study local polynomial inference in the general heteroskedastic nonparametric regression model:
where $(Y,X)$ is a pair of random variables with distribution $F$. The parameter of interest is the level or derivative of the regression function at $X = \mathsf{x}$:
where the evaluation point $\mathsf{x}$ may be in the interior or on the boundary of the support of $X$. We drop the evaluation point from the notation when possible, and employ the usual convention $\mu=\mu^{(0)}$. Derivatives at boundary points are defined as one-sided derivatives from the interior. Given a random sample $(Y_1,X_1),\dots,(Y_n,X_n)$ of size $n$ from $F$, we investigate the quality of statistical inference for $\mu^{(\nu)}$ when using kernel-based local polynomial regression methods Fan-Gijbels1996_book,Fan-Yao2005_book, focusing in particular on higher-order distributional properties of $t$ statistics as well as on coverage error and length of Wald-type confidence interval estimators. We also employ our results to compare and optimize inference procedures for empirical practice and to shed light on the sometimes underappreciated gap between point estimation and inference.
Our main technical contributions are novel Edgeworth expansions for local polynomial based Wald-type $t$ statistics of the form
for different choices of point estimator $\hat{\theta}$ and standard error estimator $\hat{\vartheta}$ detailed in Section (ref). We study the accuracy of the Gaussian approximation to the distribution of such $T$ and the error in coverage probability of their dual confidence interval estimators. Our expansions capture the dependence on implementation choices including the polynomial order, the kernel function, and the bandwidth sequence.
Edgeworth expansions are a long-standing tool for more detailed (higher-order) analyses of asymptotic distributional approximations, named after the author of a series of papers on the idea, beginning with Edgeworth1883 and treated more extensively in Edgeworth1906_JRSS. See Hall1992_book for a textbook review. Informally, an Edgeworth expansion characterizes the leading terms of the difference between the distribution of $T$ and the Gaussian distribution, denoted $\Phi(z)$, for $z \in \mathbb{R}$. That is, an Edgeworth expansion gives the leading terms $E_{T,F}(z)$ and the rate $r_{T,F}$ (which depend on the distribution generating the data, the specific $t$ statistic at issue, along with $n$, $\mathsf{x}$, and other particulars) such that
where $\mathbb{P}_F$ is the probability law when $F$ is the true data generating process.
We improve on prior work on valid Edgeworth expansions for nonparametric kernel-based regression in three ways: (i) the expansions hold uniformly over a class of data-generating processes (instead of only for one $F$), (ii) the uniform kernel is allowed (instead of only for kernel functions with sufficient variation), and (iii) the expansions hold for any derivative $\nu \geq 0$ (instead of only for the level $\nu = 0$). As discussed below, these improvements offer new theoretical and practical conclusions.
Edgeworth expansions are almost always established pointwise in the underlying distribution, that is, for a single, fixed $F$, as in (ref). Indeed, standard references on the subject Bhattacharya-Rao1976_book,Hall1992_book do not even mention uniformity. However, $F$ is unknown, and a researcher would like some assurances that their inference is equally accurate regardless of the specific underlying data generating process. This motivates expansions that are valid uniformly over a class of plausible distributions for the data, denoted $\mathscr{F}_S$, encoding the researcher's statistical model, accompanying assumptions, and the empirical regularities of the application of interest. Thus, instead of (ref), in Section (ref), we prove
We also characterize the worst-case rate $r_{T} = \inf_{F \in \mathscr{F}_S} r_{T,F}$ of distributional approximation over $\mathscr{F}_S$. The specific class $\mathscr{F}_S$ we consider, defined precisely in Section (ref), matches standard empirical settings employing kernel-based nonparametric inference for $\mu^{(\nu)}$, and therefore our theoretical and methodological results speak directly to common practice. Uniformly valid expansions have some precedence in the literature when studying notions of optimality, perhaps originating with Beran1982_AoS, but these results are rare and confined to parametric models. Our corresponding uniform results for nonparametric kernel-smoothing do not appear to have a direct antecedent in the literature.
Second, the uniform kernel is ruled out in all prior work on Edgeworth expansions for kernel-based nonparametrics, both for density estimation Hall1991_Statistics,Hall1992_AoS_density,Hall1992_book and regression Chen-Qin2000_Bmka,Chen-Qin2002_SJS,Calonico-Cattaneo-Farrell2018_JASA, due to a technical limitation in the proofs that we overcome. Other work on nonparametric regression has assumed away the issue by studying non-random designs Hall1992_AoS_regression,Neumann1997_Statistics. In fact, Hall1991_Statistics conjectured that valid Edgeworth expansions would require techniques for lattice-valued random variables if the uniform kernel was used. On the contrary, we show that such techniques are not needed. Allowing for the uniform kernel is important for empirical work because it is the optimal kernel shape in terms of minimizing interval length (as discussed in Section (ref)) and because unweighted local least squares regression is a popular choice in some applications.
Finally, inference on derivatives of the regression function, again ignored in prior work, is a common task in empirical work and therefore it is valuable to have valid Edgeworth expansions and implementation guidance specifically for this case, including inference-optimal bandwidth selection. Moreover, considering derivatives yields several interesting theoretical conclusions, highlighting the difference between point estimation and inference: we find not only that the rate of the inference-optimal bandwidth does not depend on the specific derivative order $\nu$ being considered, analogous to the well-known result for mean squared error (MSE) optimal bandwidth, but also that the rate for inference itself does not depend on $\nu$, in sharp contrast to the MSE of the point estimator.
The main result, a generic Edgeworth expansion encompassing all three of these contributions, is Theorem (ref) in Section (ref). We then use this general result to examine the error in coverage probability of confidence interval estimators dual to each $t$ statistic, and the roles of smoothing bias and Studentization in both the distributional approximation and the coverage error of the confidence intervals.
The role of bias is concretized in Section (ref). Given a level of smoothness of the unknown function $\mu$ and polynomial order of the local polynomial procedure $\hat{\theta}$, the nonparametric bias must be removed for valid inference. The method of robust bias correction (RBC) addresses this issue by incorporating explicit bias estimation into the centering $\hat{\theta}$ and then also adjusting the scale $\hat{\vartheta}$ to account for the additional variability introduced by the bias estimation Calonico-Cattaneo-Titiunik_2014_ECMA,Calonico-Cattaneo-Farrell2018_JASA. An alternative principled inference method relies on removing the bias by shrinking the bandwidth used when conducting inference, often called undersmoothing. Other ad-hoc inference approaches rely on either upper bounding the bias, inflating the scale of the $t$ statistic, or simply ignoring the bias altogether. Using our higher-order expansions, we show that RBC leads to demonstrable higher-order superior inference for $\mu^{(\nu)}$ relative to the other approaches in the literature.
Our results also show that the choice of Studentization $\hat{\vartheta}$ is crucial for good higher-order properties. This is in contrast to first order approximations, where only consistency of the standard errors is required. An important finding here is that using asymptotic approximations to the variance of $\hat{\theta}$ will increase the leading remainder terms $E_{T,F}(z)$ and hence also coverage error. Using fixed-$n$ Studentization, where $\hat{\vartheta}$ directly estimates the variability of $\hat{\theta}$, completely removes these errors. This result was first proved in Calonico-Cattaneo-Farrell2018_JASA, but only pointwise in $F$ and excluding the uniform kernel and derivatives of $\mu$. Section (ref) shows this in full generality. Failure to account for the effect of using asymptotic variance approximations has lead to some confusion in the prior literature: for example, Chen-Qin2002_SJS found inflated coverage error at boundary points and Hall1992_AoS_density found that undersmoothing provides more accurate coverage than bias correction, but both conclusions are due to improper Studentization.
A key practical consequence of the foregoing is that RBC with fixed-$n$ Studentization has leading remainder terms $E_{T,F}(z)$ and rate $r_{T,F}$, of the expansion (ref), that vanish at least as fast as, and often strictly faster than, undersmoothing-based approaches, both at interior and boundary evaluation points and for any derivative $\nu$. Intuitively, this holds because RBC exploits all available smoothness to remove bias, but is not punished (in rates) if no additional smoothness is available to remove bias. Section (ref) discusses novel implementation of RBC intervals, giving inference-optimal bandwidth and kernel choices that further improve the coverage properties and length of RBC intervals.
More broadly, our results speak to the sometimes neglected gap between point estimation and inference. Implementations focused on optimizing point estimation may not deliver optimal, or even valid, inference. In particular, they need not proceed at the same rate, and perhaps more surprisingly, the inference rate can be faster: the rate $r_{T,F}$ at which the distribution of $\hat{\theta}$ collapses to its asymptotic value (namely $\Phi(\cdot)$) can be faster than the rate at which $\hat{\theta}$ itself collapses to its asymptotic value ($\mu^{(\nu)}$). Indeed, there are cases where a bandwidth choice yields the fastest possible inference rate $r_{T,F}$ but yields invalid point estimation. This is the reverse of the better-known fact that using the estimation-optimal bandwidth (minimizing mean squared error) yields invalid inference. Rate optimality is not as well studied for inference as it is for estimation, but Section (ref) follows Hall-Jing1995_AoS to develop minimax optimal rates in the sense of achieving the fastest (minimal) rate at which the worst-case (maximal) coverage error vanishes and finds that RBC attains this rate.
The paper closes with simulation evidence supporting our theoretical and methodological work reported in Section (ref), and a brief conclusion in Section (ref). An appendix contains formulas omitted to improve the exposition, while an online supplement gives all proofs, detailed simulation results, and other methodological results. Software implementing our main results is provided in R and Stata Calonico-Cattaneo-Farrell2019_JSS. Last but not least, some of the ideas in this paper have been applied to causal inference and treatment effect estimation in the context of regression discontinuity designs in Calonico-Cattaneo-Farrell2020_EJ.
We define the class $\mathscr{F}_S$ of distributions for the pair $(Y,X)$ and make precise the local polynomial point estimator $\hat{\theta}$ and scale estimator $\hat{\vartheta}$ of the $t$ statistic (ref). The class $\mathscr{F}_S$ is determined through the following assumption. (Recall that derivatives at the boundary of the support of $X$ correspond to one-sided derivatives from the interior of the support.)
These conditions are not materially stronger than usual in kernel-based nonparametric settings. The restrictions on densities and moments are imposed to achieve uniform validity of Edgeworth expansions. The smoothness condition on $\mu$ plays a key role: the assumed smoothness, captured by $S$ and $s$, and its relationship to the smoothness utilized in estimation, will be important for coverage error.
We consider several options for the elements of the $t$ statistic $T = (\hat{\theta} - \mu^{(\nu)})/\hat{\vartheta}$ given in (ref). The starting point is the standard local polynomial regression point estimate of $\mu^{(\nu)}$. See Fan-Gijbels1996_book for an introduction. We index the classical local polynomial estimate by $p\in\mathbb{Z}_+$, the order of the polynomial used, assumed to be at least $\nu$. Suppressing the dependence on $\mathsf{x}$ to simplify notation, we therefore set
where $K$ is a kernel or weighting function, $h =h(n) \to 0$ is a bandwidth sequence, $X_{h,i} = (X_i - \mathsf{x})/h$, $\bm{r}_p(u) = (1, u, u^2, \ldots, u^p)'$, \[\bm{\bm{\Gamma}} = \frac{1}{nh}\sum_{i=1}^n K(X_{h,i})\bm{r}_p(X_{h,i})\bm{r}_p( X_{h,i} )', \quad \bm{\bm{\Omega}} = \frac{1}{h} \left[ K( X_{h,1})\bm{r}_p( X_{h,1} ), \ldots, K( X_{h,n} )\bm{r}_p( X_{h,n} )\right], \] $\bm{e}_\nu$ is the $(p+1)$-vector with a one in the $(\nu+1)^{\text{th}}$ position and zeros in the rest, and $\bm{Y} = (Y_1, \ldots, Y_n)'$.
The point estimator $\hat{\theta}$ in $T$ is then finalized depending on how the smoothing bias is to be accounted for. The traditional approach takes $\hat{\theta}=\hat{\mu}_p^{(\nu)}$, and then for inference to be valid undersmoothing is required. Explicit bias correction incorporates into $\hat{\theta}$ an estimate of the leading bias term of $\hat{\mu}_p^{(\nu)}$. Both approaches are motivated by the fact that the conditional bias of $\hat{\mu}_p^{(\nu)}$ is of order $h^{p+1 - \nu}$ and given by
with $\bm{\bm{\Lambda}} = \bm{\bm{\Omega}} [ X_{h,1}^{p+1}, \cdots, X_{h,n}^{p+1}]'/n$, provided $p-\nu$ is odd and $p\leqS-1$, the standard setting in the literature. Section (ref) details other cases for $p$ and $S$. Throughout, asymptotic orders and their in-probability versions always hold uniformly in $\mathscr{F}_S$, as required by our framework: for example, $A_n = o_\mathbb{P}(a_n)$ means $\sup_{F \in \mathscr{F}_S} \mathbb{P}_F [\vert A_n/a_n \vert > \epsilon] \to 0$ for every $\epsilon > 0$. Limits are taken as $n\to\infty$ unless stated otherwise.
Undersmoothing leaves the center of the interval at $\hat{\theta}=\hat{\mu}_p^{(\nu)}$ unchanged and assumes that the bandwidth $h$ vanishes rapidly enough to render the leading term of (ref) negligible relative to the standard error of the point estimator. The term undersmoothing refers to using less nonparametric smoothing than would be optimal from a mean squared error (MSE) point estimation point of view Fan-Gijbels1996_book. The MSE-optimal bandwidth choice is the most common by far, and indeed, the default in most software. With $p\leqS-1$, the MSE-optimal bandwidth for $\hat{\mu}_p^{(\nu)}$ is well-defined whenever $\mu^{(p+1)}(\mathsf{x}) \neq 0$. However, the MSE-optimal bandwidth is too “large” for standard Gaussian inference: the bias remains first-order important when scaled by the standard deviation of the point estimator, and so valid inference requires a bandwidth that vanishes faster.
Explicit bias correction, on the other hand, subtracts an estimate of the leading term of (ref), of which only $\mu^{(p+1)}$ is unknown. Thus we have:
where $\bm{\bm{\Omega}}_{\mathtt{rbc}} = \bm{\bm{\Omega}} - \rho^{p+1} \bm{\bm{\Lambda}} \bm{e}_{p+1}' \bm{\bar{\bm{\Gamma}}}^{-1} \bm{\bar{\bm{\Omega}}}$ and $\bm{\hat{\beta}}_{p+1}$, $\bm{\bar{\bm{\Gamma}}}$, and $\bm{\bar{\bm{\Omega}}}$ are defined akin to $\bm{\hat{\beta}}_p$, $\bm{\bm{\Gamma}}$, and $\bm{\bm{\Omega}}$ of (ref), but with $p+1$ in place of $p$ and a bandwidth $b := \rho^{-1}h$ instead of $h$. The parameter $\rho$ will play a key role in the Edgeworth and coverage error expansions and we will derive optimal choices below.
With the point estimator $\hat{\theta}$ defined, we now define the choice of standard errors $\hat{\vartheta}$. We will focus primarily on “fixed-$n$” Studentization, also called “preasymptotic” by Fan-Yao2005_book, which means choosing the Studentization to directly estimate $\mathbb{V}[ \hat{\theta} | X_1, \ldots, X_n ]$, a population quantity but not an asymptotic one. Such choices have superior coverage, as shown below, particularly compared to employing an estimator of an asymptotic representation of $\mathbb{V}[ \hat{\theta} | X_1, \ldots, X_n ]$. Importantly, when $\hat{\theta} = \hat{\theta}_\mathtt{rbc}$, a fixed-$n$ approach makes bias correction robust, because the Studentization accounts for the variability of bias estimation.
These fixed-$n$ variances are easy to compute based on standard least squares logic. Referring to (ref), for $\hat{\theta}=\hat{\mu}_p^{(\nu)}$,
where $\bm{\Sigma}$ is the $n$-diagonal matrix of conditional variances $v(X_i)$. This formula applies to $\hat{\theta}_\mathtt{rbc}$ as well, upon replacing $\bm{\bm{\Omega}}$ with $\bm{\bm{\Omega}}_{\mathtt{rbc}}$, because the two estimators share the same structure, as shown by comparing the second form in (ref) to (ref). The fixed-$n$ Studentization is obtained by replacing $\bm{\Sigma}$ with an appropriate plug-in estimator, and we then obtain the final $\hat{\vartheta}$ as follows:
where $\bm{\hat{\Sigma}}_p$ and $\bm{\hat{\Sigma}}_\mathtt{rbc}$ are the $n$-diagonal matrices of the squared residuals $\hat{v}(X_i) = (Y_i - \bm{r}_p(X_i)'\bm{\hat{\beta}}_p)^2$ and $\hat{v}(X_i) = (Y_i - \bm{r}_{p+1}(X_i)'\bm{\hat{\beta}}_{p+1})^2$, respectively. The above variance estimators separate explicitly the “constant” portions, denoted $\hat{\sigma}_p^2$ and $\hat{\sigma}_\mathtt{rbc}^2$, which will be used in Section (ref) for interval length optimization. More precisely, $\hat{\sigma}_p^2$ and $\hat{\sigma}_\mathtt{rbc}^2$ will both be bounded and bounded away from zero in probability under our assumptions.
To complete the set of $t$ statistics under consideration, we impose the following standard conditions on the kernel function. This assumption allows for standard choices such as not only the triangular and Epanechnikov kernels, but also the uniform kernel.
We now give the main technical result of this paper: a uniformly valid, generic Edgeworth expansion as in (ref), for the $t$-statistic $T$ in (ref) when using local polynomial regression methods as described in the previous section. To state the result we need some notation. Here we give only what is needed conceptually, leaving cumbersome formulas to the appendix. The terms of the Edgeworth expansion are defined as
where $z$ is the point of evaluation of the distribution, $\Psi_{T,F}$ denotes the generic non-random (fixed-$n$) bias of the $\sqrt{nh^{1+2\nu}}$-scaled numerator of $T$, $\lambda_{T,F}$ denotes the mismatch between the variance of the numerator of the $t$-statistic and the population standardization used, and the six terms $\omega_{k,T,F}(z)$, $k=1,2, \ldots, 6$, are non-random functions bounded uniformly in $\mathscr{F}_S$, and bounded away from zero for at least one $F\in\mathscr{F}_S$. Section (ref) provides further details on $\Psi_{T,F}$ and Section (ref) discusses $\lambda_{T,F}$. The quantities $\omega_{k,T,F}(z)$, $k=1,2, \ldots, 6$ are relatively less important, beyond their parity, because they cannot be altered by implementation choices.
We then have the following result (Theorem (ref)), establishing (ref). This result is general, covering interior and boundary points, $p-\nu$ even and odd, any derivative $\nu \geq 0$, and any combination of $p$ and $S$. Different cases for each of these primarily affect the expansion, and the final rates, through the bias $\Psi_{T,F}$, as explored in the next section. The conditions imposed are strengthened relative to typical pointwise first-order analyses only by $\log(nh)$ factors on the bandwidth(s) and the other uniformity requirements of Assumption (ref). (Recall that asymptotic orders and their in-probability versions are always required to hold uniformly in $\mathscr{F}_S$ throughout.)
A crucial piece in the proof of Theorem (ref) is establishing that the appropriate Cram\'er's condition holds under Assumption (ref), and in particular the linear independence condition. Such linear independence fails when $K$ is uniform and $u$ runs over the support of $K(u)$, and this failure has prevented the uniform kernel from being covered by past work. Our key insight is that previous approaches ignored the region outside the support of $K(\cdot)$ but inside the neighborhood of Assumption (ref). Loosely speaking, $(1, K(u), uK(u), \ldots)'$ may be linearly dependent on $u \in [-1,1]$ (when $K$ is uniform), but $(1, K(\frac{x - \mathsf{x}}{h}), (\frac{x - \mathsf{x}}{h}) K(\frac{x - \mathsf{x}}{h}), \ldots)'$ is linearly independent on $x$ in a fixed neighborhood of $\mathsf{x}$. This allows us to verify Cram\'er's condition. See the supplement for details.
In practice, the error in coverage probability of two-sided interval estimators may be more directly relevant than the distributional approximation of the Edgeworth expansion. We therefore turn to interval estimators dual to each $t$ statistic, given by
where $z_l$ and $z_u$ denote chosen quantiles. Our starting point is a generic coverage error expansion for confidence intervals $I$, dual to a given $T$, which follows immediately from Theorem (ref) by evaluating the Edgeworth expansion at the interval quantiles (see the supplement).
This result is as general as Theorem (ref). The uniform-in-$\mathscr{F}_S$ rate $r_I$ is the slowest vanishing of the rates of each term in the Edgeworth expansion (ref), which without specifying any elements further, can only be known to vanish at least as fast as $r_{T} = \sup_{F \in \mathscr{F}_S} r_{T,F}$ from Theorem (ref). However, even at this level of generality, several conclusions are already evident due to the parity of the functions $\omega_k$ making up $E_{T,F}(z)$ and hence $C_{I,F}(z_l,z_u)$. First, regarding the choice of quantiles, we recover the classical finding that symmetric intervals, where $z_l = -z_u$, have superior coverage properties, because $\omega_1$ and $\omega_2$ are even functions of $z$. Asymmetric choices that still have $\Phi(z_u) - \Phi(z_l) = 1-\alpha$ can yield correct coverage, but the error will vanish more slowly, whereas other choices will not yield uniformly correct coverage. Bootstrap-based quantiles will, in general, not improve coverage error rates in nonparametric contexts beyond the symmetric case Hall1992_AoS_density, and can in fact be detrimental for coverage error Hall-Kang2001_AoS. Second, the remaining $w_k$ functions are odd, and therefore to obtain better coverage properties we should focus on intervals with small (rapidly vanishing) $\Psi_{T,F}$ and $\lambda_{T,F}$. The upcoming subsections discuss each of these pieces in turn.
Our expansions highlight the conceptual gap between point estimation and inference. The rate at which the distribution of $\hat{\theta}$ collapses to its asymptotic value ($\Phi(\cdot)$) can be faster than the rate at which the point estimator $\hat{\theta}$ itself collapses to its asymptotic value ($\mu^{(\nu)}$). Moreover, it is possible that coverage error may vanish even if mean squared error does not, and vice versa. One direction of this phenomenon captures the well-known result that the coverage error of a confidence interval centered at the MSE-optimal point estimator does not vanish. That is, $\hat{\mu}_p^{(\nu)}$ in (ref) using the MSE-optimal bandwidth $h_\mathtt{mse} = H_\mathtt{mse} n^{-1/(2p+3)}$, for some constant $H_\mathtt{mse}$, is optimal for point estimation given a fixed $p$, but \[ \sup_{F \in \mathscr{F}_S} \left| \mathbb{P}_F \left[ \mu^{(\nu)} \in \left\{\hat{\mu}_p^{(\nu)} \pm z_{\alpha/2} \hat{\sigma}_p H_\mathtt{mse}^{-1/2} n^{-1/2 + (1 + 2\nu)/(4p+6)} \right\} \right] - (1-\alpha) \right| \; \asymp \; 1,\] where $a \asymp b$ denotes that $a \leq C_1 b$ and $b \leq C_1 a$ for some constants $C_1$ and $C_2$.
The other direction may be more surprising and novel: we find that the variance of $\hat{\theta}$ can be too large for mean-square consistency, but nonetheless be captured well enough by $\hat{\vartheta}$ for valid inference. For example, consider inference on $\mu^{(1)}(\mathsf{x})$ using $I_p$ with local linear regression ($p=1$). Choosing $h \asymp n^{-1/3}$ yields $r_{I_p} \asymp n^{-2/3}$, which is the fastest attainable rate for $I_p$ in this case, but also gives $\mathbb{V}[\hat{\mu}_p^{(\nu)} | X_1, \ldots, X_n ] \asymp_\mathbb{P} (nh^{1 + 2v})^{-1} \asymp 1$, and thus $\hat{\mu}_p^{(1)}$ is not consistent in mean square. Therefore, we found a confidence interval that is optimal for coverage of $\mu^{(1)}$, but implicitly relies on a point estimator that is not even consistent in mean square.
We now give details for the bias term, $\Psi_{T,F}$, highlighting three main points. First, the rate at which $\Psi_{T,F}$ vanishes does not depend on the derivative $\nu$. Second, we establish that performing bias correction never slows the rate at which $\Psi_{T,F}$ vanishes. The third goal is then practical: we spell out several cases of the rates and constants for the bias of $\hat{\theta}_\mathtt{rbc}$ so that we may use these for bandwidth and kernel selection later.
To describe $\Psi_{T_p,F}$, the bias term for $T_p$, let $\bm{\beta}_{p}$ be the $p+1$ vector with $(j+1)$ element equal to $\mu^{(j)}(\mathsf{x})/j!$ for $j = 0, 1, \ldots, p$ as long as $j \leq S$, and zero otherwise, and $\bm{B}_{p}$ as the $n$-vector with $i^{\text{th}}$ entry $[\mu(X_i) - \bm{r}_{p}(X_i - \mathsf{x})'\bm{\beta}_{p}]$. Then,
Turning to bias correction, define $\bm{\beta}_{p+1}$ and $\bm{B}_{p+1}$ as above, but with $p+1$ in place of $p$ in all cases. Then, using the definition of $\bm{\bm{\Omega}}_{\mathtt{rbc}}$ in (ref),
These bias terms are non-random but otherwise non-asymptotic: all expectations are fixed-$n$ and we have not done the typical Taylor expansion. The derivative $\nu$ only appears in the constant term $\nu! \bm{e}_\nu$, and therefore the rate at which $\Psi_{T_p,F}$ vanishes does not depend on the derivative being estimated. Intuitively, this can be seen from the second form for $\hat{\mu}_p^{(\nu)}$ in (ref), $n^{-1} h^{-\nu} \nu! \bm{e}_\nu'\bm{\bm{\Gamma}}^{-1} \bm{\bm{\Omega}} \bm{Y}$, coupled with rate $\sqrt{nh^{1+2\nu}}$ of the Studentizations of (ref): together, these account for the derivative, and distributional properties of $\bm{\bm{\Gamma}}^{-1} \bm{\bm{\Omega}} \bm{Y}$ are left independent of $\nu$; the first conclusion of this subsection.
The rate of convergence (to zero) of $\Psi_{T_p,F}$ or $\Psi_{T_\mathtt{rbc},F}$ can be deduced by first expanding $\mu(X_i)$ entering $\bm{B}_p$ and $\bm{B}_{p+1}$ around $\mathsf{x}$, and then specializing to a given $p$ and $S$. For any $p$, we have \[ \mu(X_i) - \bm{r}_p(X_i - \mathsf{x})'\bm{\beta}_p = \sum_{k= S \wedge p + 1}^S \frac{1}{k!} (X_i - \mathsf{x})^k \mu^{(k)}(\mathsf{x}) + \frac{1}{S!} (X_i - \mathsf{x})^S \left( \mu^{(S)}(\bar{x}) - \mu^{(S)}(\mathsf{x}) \right), \] where the summation is taken to be zero if $p\geqS$. To obtain the final rate, this expansion is substituted into $\bm{B}_p$ and the leading terms are identified by stabilizing the expectation of the terms involving $(X_i - \mathsf{x})^k$ by writing $h^k (X_{h,i})^k$, thus isolating the rate. The rate will depend on the smoothness, location of $\mathsf{x}$, parity of $p-\nu$, and the bandwidth $h$. For $\bm{B}_{p+1}$, replace $p$ with $p+1$ everywhere and use $b$ in place of $h$ in the second term. The supplement gives complete details.
Our second point is that $\Psi_{T_\mathtt{rbc},F} = O(\Psi_{T_p,F})$, which follows from the expansion above and taking $\rho$ bounded and bounded away from zero. First, observe from the Taylor expansion applied to (ref) that $\Psi_{T_\mathtt{rbc},F}$ depends on higher order derivatives than $\Psi_{T_p,F}$, which follows from applying the Taylor expansion to (ref), and therefore stabilizing leads to higher powers of $h$. Intuitively, the bias of $\hat{\mu}_p^{(\nu)}$ is the product of the rate $h^{p+1}$ and the constant targeted by bias correction. Therefore, the bias of $\hat{\theta}_\mathtt{rbc}$ is at most $h^{p+1}$ times the bias of the bias correction plus the higher order term of (ref). For a fixed sequence $h$, neither of these can be greater than $h^{p+1}$. Second, the rate for $\Psi_{T_\mathtt{rbc},F}$ cannot be improved by letting $\rho = h/b$ vanish or diverge: $\rho$ vanishing decreases the second term, but the first term is unchanged, while letting $\rho$ diverge can only inflate the second term. Further, diverging $\rho$ renders the effective sample size $nb$, which is smaller than $nh$, which would only inflate the Edgeworth expansion terms without reducing bias (hence the restriction in Theorem (ref) to bounded $\rho$).
Therefore, in optimizing inference later on, we will focus on $\hat{\theta}_\mathtt{rbc}$ and take $\rho$ bounded and bounded away from zero. We need the leading bias constants for this case, which follow from carrying on the Taylor expansion completely in (ref). The bias is always of the form \[\Psi_{T_\mathtt{rbc},F} = O(\sqrt{nh} h^\zeta)\] for an exponent $\zeta$ that depends on the location of $\mathsf{x}$, the parity of $p-\nu$, and the smoothness $S$. A complete list of $\zeta$ is shown in Table (ref). From there, we see that if $p$ is large enough relative to $S$ (how large depends on the specific case), then $\zeta = S + s$, implying $\Psi_{T_\mathtt{rbc},F} = O(\sqrt{nh} h^{S + s})$.
The more empirically relevant case is to treat $p$ as fixed and smaller than $S$, specifically $p \leq S - 3$ for interior $\mathsf{x}$ with $p-\nu$ odd and $p \leq S-2$ otherwise (i.e. for boundary points or if $\mathsf{x}$ is an interior point with $p-\nu$ even). In these cases, we can use the Taylor expansion above to characterize the leading term, and write \[\Psi_{T_\mathtt{rbc},F} = \sqrt{nh} h^{\zeta} \psi_{T_\mathtt{rbc},F} [1 + o(1)],\] where $\zeta = p+3$ for interior $\mathsf{x}$ with $p-\nu$ odd and $p+2$ otherwise. The term $\psi_{T_\mathtt{rbc},F}$ will be referred to as the constant term for simplicity, though technically it is a non-random sequence with known form, uniformly bounded in $\mathscr{F}_S$, and nonzero for some $F \in \mathscr{F}_S$. Referring to Table (ref) for the different cases, $\psi_{T_\mathtt{rbc},F}$ can be
where $\bm{\bm{\Lambda}}_k = \bm{\bm{\Omega}} [ X_{h,1}^{p+k}, \cdots, X_{h,n}^{p+k}]'/n$ and $\bm{\bar{\bm{\Lambda}}}_k = \bm{\bar{\bm{\Omega}}} [ X_{b,1}^{p+1+k}, \ldots, X_{b,n}^{p+1+k}]'/n$, and hence in particular $\bm{\bm{\Lambda}}_1 \equiv \bm{\bm{\Lambda}}$ as defined in Section (ref).
In contrast to first order distributional analysis, where only consistency is required, the choice of scaling, or Studentization, is crucial for higher order properties. Our detailed expansions show that, in general, there are two types of higher-order terms that arise due to Studentization. One is the unavoidable estimation error incurred when replacing any population quantity with a feasible counterpart. The second error arises from the difference between the population variability of the centering $\hat{\theta}$ and the population standardization chosen as the target. This second type of error is what is captured by $\lambda_{T,F}$, and the most important conclusion is that the fixed-$n$ standard errors in (ref) achieve $\lambda_{T,F} \equiv 0$, and are therefore demonstrably superior choices for inference. That is, there should not be a “mismatch” between the population variability of the $t$ statistic numerator and the population standardization.
Using an asymptotic approximation to $\mathbb{V}[ \hat{\theta} | X_1, \ldots, X_n ]$ may yield nonzero $\lambda_{T,F}$, and thus the distributional approximation (and coverage) will suffer. There are too many options to treat comprehensively, but several points warrant discussion. In general, if the chosen standard errors are consistent, $\lambda_{T,F}$ has the form $\lambda_{T,F} = l_n L$, for a rate $l_n \to 0$ and a sequence $L$ that is bounded and bounded away from zero, a “constant”, capturing the difference between the variance of the numerator of the $t$-statistic and the population standardization chosen.
At boundary points the use of asymptotic approximations can be particularly deleterious for coverage, and this has lead to some confusion in the literature. A headline finding of Chen-Qin2000_Bmka is that an empirical likelihood confidence interval estimator has coverage error of the same order at interior and boundary points, which is claimed (in the abstract) to be a “significant improvement over confidence intervals based directly on the asymptotic normal distribution”. This claim is based on work by the same authors Chen-Qin2002_SJS who study, in our notation, the interval with centering $\hat{\theta} = \hat{\mu}_1^{(0)}$ and scaling $\hat{\vartheta} = (nh)^{-1/2} \hat{v}(\mathsf{x}) \hat{f}(\mathsf{x})^{-1} \mathcal{V}$, for $\hat{v}(\mathsf{x})$, $\hat{f}(\mathsf{x})$, and $\mathcal{V}$ given therein, where $v(\mathsf{x}) f(\mathsf{x})^{-1} \mathcal{V}$ is the probability limit of $\mathbb{V}[ (nh)^{1/2} \hat{\mu}_1^{(0)} \mid X_1, \ldots, X_n ]$. They find that $\lambda_{T,F} = l_n L$ holds with $l_n = h$ at boundary points, meaning greatly increased coverage error. Concerned that this result is due to estimation error, they confirm that $l_n = h$ holds with the infeasible standardization $\vartheta = (nh)^{-1/2} v(\mathsf{x}) f(\mathsf{x})^{-1} \mathcal{V}$. This neglects the fact that $\lambda_{T,F}$ captures only the estimation error, not the “mismatch” error, and their result is entirely due to using an asymptotic standardization as opposed to a fixed-$n$ one, and thus empirical likelihood, in particular, does not offer higher-order improvements over normality-based intervals.
Explicit bias correction was claimed by Hall1992_AoS_density to be inferior to undersmoothing for inference; a finding also based entirely on using an asymptotic standardization. In this case, nonrobust bias correction was studied, which pairs $\hat{\theta}_\mathtt{rbc}$ with $\hat{\sigma}_p$. This is valid to first order if $\rho = o(1)$, because then $\mathbb{V}[ \hat{\theta}_\mathtt{rbc} | X_1, \ldots, X_n ] = \mathbb{V}[ \hat{\mu}_p^{(\nu)} \mid X_1, \ldots, X_n ] = o_\mathbb{P}(n^{-1}h^{-1-2\nu})$. However, higher order expansions find $\lambda_{T,F} = \rho ^{p+2}(L_1 + \rho^{p+2} L_2)$, where $L_1$ captures the (scaled) covariance between $\hat{\mu}^{(\nu)}$ and $\hat{\mu}^{(p+1)}$ and $L_2$ the variance of $\hat{\mu}^{(p+1)}$. These terms lead Hall1992_AoS_density to conclude that bias correction is inferior to undersmoothing, which Calonico-Cattaneo-Farrell2018_JASA later showed is not true for robust bias correction. Our results extend this conclusion to hold for derivatives, boundary points, all smoothness cases, and uniformly in $\mathscr{F}_S$, while also allowing for the uniform kernel.
We turn to optimizing inference in practice, using the conclusions from the previous sections. Collectively, the previous sections imply that the best coverage will be from using symmetric RBC intervals, i.e. those with $z_l = -z_u = z_{\alpha/2}=\Phi^{-1}(\alpha/2)$, $\hat{\theta}_\mathtt{rbc}$ as in (ref), $\hat{\vartheta}_\mathtt{rbc}$ as in (ref), and $\rho$ bounded and bounded away from zero (implying $h = \rho b$). With an eye toward empirical work, we assume in this section that $p$ is fixed and small compared to $S$. The other cases detailed in Section (ref) are of relatively little practical value. In practice researchers first choose $p$ and then conduct inference based on that choice (witness the ubiquity of local linear regression and cubic splines).
Letting \[I_\mathtt{rbc}(h) = \Big[ \hat{\theta}_\mathtt{rbc} + z_{\alpha/2} \; \hat{\vartheta}_\mathtt{rbc} \ , \ \hat{\theta}_\mathtt{rbc} - z_{\alpha/2} \; \hat{\vartheta}_\mathtt{rbc} \Big] \] denote the recommended RBC confidence interval, now with its dependence on the bandwidth $h$ explicit to enhance the exposition, we readily deduce from Corollary (ref) that
where the coverage error rate is $r_{\mathtt{rbc}} = \max\{(nh)^{-1}, nh^{1+2\zeta}, h^\zeta \}$, with $\zeta = p+3$ if $p-\nu$ is odd and $\mathsf{x}$ is a boundary point, or $\zeta =p+2$ otherwise. Furthermore, its interval length is
Notice that the rate of contraction of length does depend on $\nu$, while the coverage error rate does not.
In the next two subsection we use the above two displays, (ref) and (ref), to choose the bandwidth parameters $h$ and $\rho = h/b$, and the kernel shape. Before any choices can be made, the researcher must decide on the usual size versus power trade off. In our context, this translates to the relative value they place on coverage error, the discrepancy from nominal level, versus interval length. Because we give the first characterizations of coverage error in many cases, and the first uniformly valid ones, this issue can now be studied in detail: our theoretical ideas can inform this trade off, providing new insights to consider, as well as guiding implementation given a preference for coverage error and length.
At one extreme is the approach that requires only that the interval is not anti-conservative, and then minimizes (expected) length. In this case, a shorter interval that uniformly over-covers is preferred to an interval that is longer but has correct coverage asymptotically. Our results lead one to consider the other extreme: minimize the coverage error directly, and only after optimize length. That is, seek for the confidence interval $I$ such that, in the notation of Corollary (ref), $r_I$ vanishes as fast as possible. In applications, an interval with a faster decaying coverage error may approximate its nominal level more closely in finite samples. Such approach focuses on the accuracy of the Gaussian approximation for coverage error, and thus for inference. However, both of these extremes may be unappealing in practice because neither may be optimal from a coverage-length (or, perhaps, size-power for the dual hypothesis test) perspective. Therefore, we will also consider compromises, trading off between coverage error and interval length. One option is to minimize length among consistent interval estimators: seek the shortest interval such that $r_I = o(1)$. In the context of kernel-based nonparametrics, interval estimators with good control of worst-case coverage are able to use larger bandwidths in general, and are thus shorter in large samples; an analogue to the adage that “similar tests have higher power”. In general, we will let the user determine a trade off between the two and thus find a bandwidth choice to implement their preference.
We now focus on choosing the bandwidth $h$ optimally, leaving $\rho$ and $K$ to the next section. With pragmatism in mind, we restrict attention to bandwidth sequences that are polynomial in $n$, that is, of the form $h = H n^{-\eta}$ for some constants $H>0$ and $\eta >0$. For implementation purposes, we optimize $C_{I_\mathtt{rbc}(h),F}(z_{\alpha/2}, -z_{\alpha/2})$ pointwise in $F$. The optimal bandwidths will be functions of $F$ and their implementations are functions of the data, which are draws from $F$; neither depend explicitly upon $\mathscr{F}_S$. The resulting coverage error rates still hold uniformly, because the bandwidths are of the form $h = H n^{-\eta}$, where $\eta$ does not depend on $F$ and $H$ is well-behaved uniformly in $\mathscr{F}_S$. We will focus on cases where coverage is consistent, leveraging our new higher-order results in this paper.
An obvious candidate for $h$ in applications is the classical MSE-optimal choice, denoted $h_\mathtt{mse}$, for the point estimator $\hat{\mu}_p^{(v)}(x)$ used as part of the centering of the confidence interval $I_\mathtt{rbc}(h)$. This bandwidth choice is popular and readily available in most statistical software. Although designed to optimize point estimation, our theoretical results show that it yields valid robust bias corrected inference, that is, $\sup_{F \in \mathscr{F}_S} \; | \mathbb{P}_F [ \mu^{(\nu)}(\mathsf{x}) \in I_\mathtt{rbc}(h_\mathtt{mse}) ] - (1-\alpha) | \to 0$, in contrast to the traditional interval $I_p(h_\mathtt{mse})$, which undercovers. This gives a principled endorsement for using $h_\mathtt{mse}$ coupled with robust bias correction in applications, if a researcher wishes to optimize point estimation instead of inference when choosing the bandwidth $h$. To be more precise, our results give formal justification (and demonstrate higher-order coverage improvements) for reporting $\hat{\mu}_p^{(\nu)}(\mathsf{x})$ along with $I_\mathtt{rbc}(h_\mathtt{mse})$, both implemented using the same bandwidth $h_\mathtt{mse}$, that is, pairing an MSE-optimal point estimator with a valid measure of uncertainty that uses the same samples. In fact, an interesting consequence of our results is that for interior points and local linear regression ($p=1$), $I_\mathtt{rbc}(h_\mathtt{mse})$ has coverage error that vanishes as fast as possible: for this special case, both the mean squared error and coverage error are optimal in rates upon setting $h = H n^{-1/(2p+3)}$ for a constant $H > 0$. In other cases, coverage of confidence intervals implemented using $h_\mathtt{mse}$ remains consistent but the coverage rate is suboptimal.
To see this, we now turn to inference-optimal bandwidths. We start with the point of view that minimizing coverage error alone is the goal and therefore we choose $h$ by minimizing the terms of (ref). This means setting $h_\mathtt{rbc} = H n^{-\eta_\mathtt{rbc}}$ for $\eta_\mathtt{rbc} = 1/(p+4)$ for interior $\mathsf{x}$ with $p-\nu$ odd and $\eta_\mathtt{rbc} = 1/(p+3)$ otherwise: Corollary (ref) holds for $I_\mathtt{rbc}(h_\mathtt{rbc})$ with rates $r_\mathtt{rbc} = n^{-(p+3)/(p+4)}$ and $r_\mathtt{rbc} = n^{-(p+2)/(p+3)}$, respectively. In terms of rates, $h_\mathtt{rbc}$ balances the variance and bias of the point estimator, instead of the squared bias as in MSE optimality.
A natural way of choosing the constant $H$ in practice is to minimize the constant portion of the coverage error of (ref). Plugging in $h_\mathtt{rbc} = H n^{-\eta_\mathtt{rbc}}$ and factoring out the rate we get \[ H_\mathtt{rbc} = \operatorname*{arg\,min}_{H > 0} \left\vert H^{-1} \big\{ 2 \omega_{4,\mathtt{rbc},F} \big\} + H^{1+2\zeta} \big\{ 2 \psi_{T_\mathtt{rbc},F}^2 \omega_{5,\mathtt{rbc},F} \big\} + H^\zeta \big\{ 2 \psi_{T_\mathtt{rbc},F} \omega_{6,\mathtt{rbc},F} \big\} \right\vert.\] It is straightforward to give a data-driven version of $H_\mathtt{rbc}$, and therefore of $h_\mathtt{rbc}$, because all quantities involved can be estimated. We defer the details to the supplement to conserve space. In a nutshell, plug-in estimators can be constructed, denoted by $\hat{\omega}_{4,\mathtt{rbc},F}$, $\hat{\omega}_{5,\mathtt{rbc},F}$ and $\hat{\omega}_{6,\mathtt{rbc},F}$, as well as an estimate of the bias constant, $\hat{\psi}_{\mathtt{rbc},F}$. We then numerically solve \[ \hat{H}_\mathtt{rbc} = \operatorname*{arg\,min}_{H > 0} \left\vert H^{-1} \big\{ 2 \hat{\omega}_{4,\mathtt{rbc},F} \big\} + H^{1+2\zeta} \big\{ 2 \hat{\psi}_{\mathtt{rbc},F}^2 \hat{\omega}_{5,\mathtt{rbc},F} \big\} + H^\zeta \big\{ 2 \hat{\psi}_{\mathtt{rbc},F} \hat{\omega}_{6,\mathtt{rbc},F} \big\} \right\vert.\] Because this bandwidth depends on the specific data-generating process $F$, we view it as a rule-of-thumb implementation.
As discussed above, we can also seek for a shorter interval (more power) by sacrificing coverage error (size control). Interval length (ref) is reduced for larger bandwidths, meaning smaller exponents $\eta$. Corollary (ref), or Equation (ref) specifically, shows that the smallest $\eta$ (i.e., the slowest vanishing bandwidth) such that the coverage of $I_\mathtt{rbc}(n^{-\eta})$ to be (uniformly) asymptotically correct is $\eta > (1/(1 + 2 \zeta)$, where recall that $\zeta = p+3$ for interior points with $p-\nu$ odd and $\zeta = p+2$ otherwise. Therefore, taking $h = H n^{-\eta}$ for any $\eta > (1/(1 + 2 \zeta)$ and $H>0$ results in the ideal interval given these preferences over coverage error and length.
This same idea can be extended to accomplish a trade-off between coverage error and length. Researchers may want to have an interval that is closer to nominal level, and therefore may be concerned that in finite samples an interval with coverage error only known to obey $r_\mathtt{rbc} = o(1)$ will not be satisfactory. We can therefore take $h_\mathtt{to} = H_\mathtt{to} n^{-\eta_\mathtt{to}}$ for some $\eta_\mathtt{to} \in (1/(1 + 2 \zeta) , \eta_\mathtt{rbc}]$. Note that if $\eta > \eta_\mathtt{rbc}$ (i.e.\ $h = o(h_\mathtt{rbc})$), both the rate of coverage error decay and interval length contraction can be improved. There is no well-defined optimal choice in this range of asymptotically valid options, as the choice must reflect each researcher's preference for length vs.\ coverage error. This range does not depend on $\nu$, even though the resulting length will, see (ref). This may affect how the researcher wishes to trade off the two quantities. The endpoints of the range for $\eta_\mathtt{to}$ represent preferences for only optimizing coverage error or only length.
To select the constant for this trade off, $H_\mathtt{to}$, note first that for $\eta < \eta_\mathtt{rbc}$ the middle term of the coverage error (ref) is dominant. This term, $n ^{ 1 - \eta_\mathtt{to}(1+2\zeta)} \big\{ 2 \psi_{T,F}^2 \omega_{5,T,F} \big\}$, shares the rate of the scaled, squared bias. Therefore, it is natural to balance this against the square of interval length, to match the trade off that $h_\mathtt{to}$ represents. The feasible choice of this constant, $\hat{H}_\mathtt{to}$, will also be a direct plug-in rule that uses the estimators above and a pilot version of $\hat{\sigma}_\mathtt{rbc}^2$, as well a researcher's choice of weight $\mathcal{H} \in (0,1)$ capturing their trade off between the two. Put altogether, we can then set
The resulting data-driven bandwidth choice is $\hat{h}_\mathtt{to} = \hat{H}_\mathtt{to} n^{-\eta_\mathtt{to}}$, for a choice $\eta_\mathtt{to} \in(1/(1 + 2 \zeta), \eta_\mathtt{rbc}]$, and weight $\mathcal{H} \in (0,1)$. The supplement contains details and some additional results.
To complete the implementation of $I_\mathtt{rbc}(h)$ we need to select the bias-correction bandwidth $b$, which we do in the form of $\rho = h/b$, and the kernel function $K(\cdot)$. We choose these to optimize the length (ref). With $\rho$ bounded and bounded away from zero, this choice affects only the constant portions of the coverage error expansion of $I_\mathtt{rbc}(h_\mathtt{rbc})$, in particular changing the shape of the equivalent kernel of $\hat{\theta}_\mathtt{rbc}$. For more details on equivalent kernels, see Fan-Gijbels1996_book. To find this equivalent kernel, begin by writing $\hat{\theta}_\mathtt{rbc} = \nu! \bm{e}_\nu'\bm{\bm{\Gamma}}^{-1} \bm{\bm{\Omega}}_{\mathtt{rbc}} \bm{Y} / n h^\nu$ as a weighted average of the $Y_i$. Recall that $X_{h,i} = (X_i - \mathsf{x})/h$ and similarly for $X_{b,i}$. Then,
The weights here depend on the sample, as $\bm{\bm{\Gamma}}$, $\bm{\bm{\Lambda}}$, and $\bm{\bar{\bm{\Gamma}}}$ are sample quantities. The equivalent kernel replaces these with their limiting versions (not, as elsewhere, their fixed-$n$ expectations), which we shall denote $\bm{\mathsf{G}} = f(\mathsf{x}) \int K(u)\bm{r}_p(u)\bm{r}_p(u)'du$, $\bm{\mathsf{L}} = f(\mathsf{x}) \int K(u)\bm{r}_p(u) u^{p+1} du$, and $\bm{\bar{\mathsf{G}}} = f(\mathsf{x}) \int K(u)\bm{r}_{p+1}(u)\bm{r}_{p+1}(u)'du$, respectively. The integrals are over $[-1,1]$ if $\mathsf{x}$ is an interior point and appropriately truncated when $\mathsf{x}$ is a boundary point. Under our assumptions, convergence to these limits is fast enough that, for the equivalent kernel $\mathcal{K}_\mathtt{rbc}(u; K, \rho, \nu)$ defined as \[ \mathcal{K}_\mathtt{rbc}(u; K, \rho, \nu) = \nu! \bm{e}_\nu'\bm{\mathsf{G}}^{-1} \left( K(u)\bm{r}_p( u) - \rho^{p+2} \bm{\mathsf{L}} \bm{e}_{p+1}' \bm{\bar{\mathsf{G}}}^{-1} K(u \rho)\bm{r}_{p+1}( u \rho) \right), \] and we have the representation
It follows that the (constant portion of the) asymptotic length of $I_\mathtt{rbc}(h)$ depends on $K(\cdot)$ and $\rho$ only through the specific functional $\int\big(\mathcal{K}_\mathtt{rbc}(u; K, \rho, \nu)\big)^2du$, which corresponds to the asymptotic variance.
The asymptotic variance of a local polynomial point estimator at a boundary or interior point is minimized by employing the uniform kernel Cheng-Fan-Marron1997_AoS. Therefore, to minimize the constant term of interval length we choose $\rho$, depending on $K$, to make $\mathcal{K}_\mathtt{rbc}(u; K, \rho, \nu)$ as close as possible to the optimal equivalent kernel, i.e. the $\mathcal{K}^*_{p}(u)$ induced by the uniform kernel for a given $p$. If the uniform kernel is used initially, then $\rho^*=1$ is optimal: that is, $\mathcal{K}_\mathtt{rbc}(\cdot; \mathbbm{1}\{|u|<1\}/2, 1, \nu) \equiv \mathcal{K}^*_{p+1}(\cdot)$. This highlights the importance of being able to accommodate the uniform kernel in our higher-order expansions. If a kernel other than uniform is used, we look for the optimal choice of $\rho$ by minimizing the $L_2$ distance between the induced equivalent kernel and the optimal variance-minimizing equivalent kernel, solving \[\rho^* = \operatorname*{arg\,min}_{\rho>0} \int \left| \mathcal{K}_\mathtt{rbc}\big(u; K, \rho, \nu \big) - \mathcal{K}^*_{p+1}(u) \right| ^2 du.\] This is not a sample-dependent problem, only computational. For $p-\nu$ odd, the standard case in practice, Table (ref) shows the optimal $\rho^*$, for boundary and interior points, respectively, the triangular kernel ($K(u)=(1-|u|)\mathbbm{1}(|u|\leq 1)$) and the Epanechnikov kernel ($K(u)=0.75(1-u^2)\mathbbm{1}(|u|\leq 1)$). These two are popular choices and are MSE-optimal at boundary and interior points, respectively. The shapes of the resulting equivalent kernel, $\mathcal{K}_\mathtt{rbc}(u; K, \rho^*, \nu)$, are shown in Figure (ref) for $\nu=\{0,1\}$. Note that although $\rho^*$ itself does not vary with $\nu$, the equivalent kernel shape does. Additional choices of $p$ are illustrated in the supplement.
In this section we build on Hall-Jing1995_AoS and look for a minimax result: characterizing the fastest (minimal) rate at which the worst-case (maximal) coverage error vanishes. The “optimal” interval estimator is one for which this maximal error is minimized. At an intuitive level, this corresponds to the desire for similarity in testing: the confidence interval should have “similar” coverage over the set of plausible distributions. Hall-Jing1995_AoS proposed this inference-specific notion of minimax optimality and studied it in the case of one-sided confidence intervals in the i.i.d.\ parametric location model. This problem is different from the more typical minimaxity considered for point estimation, though the latter is established for robust bias correction by Tuvaandorj2020_JoE and is discussed more broadly for local polynomials by Cheng-Fan-Marron1997_AoS and Fan-etal1997_AISM.
To state the problem more formally, let $\mathscr{I}_p$ denote a class of confidence interval estimators. We then define the minimax coverage error as
where the dependence on the fixed quantities, such as the classes $\mathscr{I}_p$ and $\mathscr{F}_S$ or the level $\alpha$, are suppressed. Our goal is to characterize the minimax optimal coverage error decay rate bound, which is the fastest vanishing sequence $r_\star=r_\star(n)$, $n\in\mathbb{N}$, such that for constants $c_1$ and $c_2$,
We have already characterized the worst-case coverage error in Corollary (ref) for the class of distributions defined in Section (ref). The key point here is that if we take $\mathscr{I}_p$ to be the class of intervals for which we studied worst-case coverage error in Corollary (ref), then we can characterize the minimax rate $r_\star$ as well as intervals which attain it. Specifically, we take $\mathscr{I}_p$ to be the Wald-type intervals of the form (ref), based on a local polynomial of degree $p$, with any choice of centering, scaling, bandwidth(s), kernel shape, and quantiles, discussed in Section (ref). This includes all those intervals dual to $t$ statistics covered by Theorem (ref), but also includes other choices which are not asymptotically level $1-\alpha$. Examples include trivial cases such improper choices of quantiles or inconsistent variance estimators, but also choices such as $I_p(h_\mathtt{mse})$, i.e., using the MSE-optimal bandwidth sequence for with centering $\hat{\mu}_p^{(\nu)}$ and scaling $\hat{\sigma}_p^2/(nh^{1+2\nu})$. We could also include other procedures, such as bootstrap based quantiles, empirically chosen bandwidths, or empirical likelihood methods, as these will not improve on the worst-case coverage error Hall1992_AoS_density,Hall-Kang2001_AoS,Chen-Qin2000_Bmka.
Crucial to proving that such an interval is minimax optimal is that the bias vanishes at the best possible rate, given the smoothness assumed ($S$) and utilized ($p$), and this in turn depends on whether $\mathsf{x}$ is an interior or boundary point. Collecting all the smoothness cases studied in Section (ref), we immediately obtain the following result (see the supplement for omitted details).
For the classes $\mathscr{F}_S$ and $\mathscr{I}_p$ considered herein, this result establishes the minimax rate bounds. The interplay between the two classes is crucial: they should be neither too “large” nor too “small” in order to obtain useful and interesting results. The larger is $\mathscr{F}_S$, the more plausible a given data set is generated by some $F \in \mathscr{F}_S$, but well known results dating back at least to Bahadur-Savage1956_AoMS show that if $\mathscr{F}_S$ is too large it is impossible to construct an “effective confidence interval” that controls the worst-case coverage. Our particular $\mathscr{F}_S$ captures common restrictions in the setting of nonparametric regression, and therefore matches empirical practice. The class $\mathscr{I}_p$ is restricted to contain Wald-type interval estimators commonly employed in practice using nonparametric kernel-based regression methods (but can be trivially extended to cover alternatives mentioned above). Recall that our goal is to identify if RBC confidence intervals improve over other options in a uniform sense, and this result is tailored to that goal.
The main message of Corollary (ref) is that $I_\mathtt{rbc}(h_\mathtt{rbc})$ is minimax optimal in all cases. This strengthens the pointwise improvement offered by robust bias correction to optimality within the class $\mathscr{I}_p$ considered here. Intuitively, this is because robust bias correction successfully exploits additional smoothness if it exists, but is not punished (in rates) if there is no such smoothness due to the change in Studentization. This can be compared to $I_p$, the classical interval that requires undersmoothing. This interval is optimal only in the case when $S$ is known so that $p$ can be chosen large enough; for a fixed $p$ that is small relative to $S$ this interval is dominated in the minimax sense.
This section presents results from a simulation study to examine the finite-sample performance of our methods. Additional results and implementation details can be found in the supplement. We focus on the performance of confidence intervals for $\mu(\mathsf{x})$ and $\mu^{(1)}(\mathsf{x})$ based on robust bias correction and traditional undersmoothing. Data is generated from model (ref), with $X_i$ uniformly distributed on $[-1,1]$, $\varepsilon$ standard normal, and \[ \mu(x) = \frac{ \sin(3\pi x/2 ) }{ 1+18x^2(\operatorname*{sgn}(x) +1) }, \] where $\operatorname*{sgn}(x)=-1$, $0$, or $-1$ according to $x>0$, $x=0$ or $x<0$, respectively. This function, which was also analyzed in Calonico-Cattaneo-Farrell2018_JASA, is displayed in Figure (ref) together with $\mu^{(1)}(x)$. By looking at different evaluation points, we will be able to capture the performance of the methods under different levels of complexity.
We show results for sample sizes $n\in\{100,250,500,750,1000,2000\}$, always with $5,000$ replications. We study inference at three evaluation points: $\mathsf{x}=-1$ (boundary point), $\mathsf{x}=-0.6$ (low curvature), and $\mathsf{x}=-0.2$ (high curvature). The supplement shows results for $\mathsf{x}\in\{0.2,0.6,1\}$. For implementation, we use $p=1$ (for $\nu=0$) and $p=2$ (for $\nu=1$) with the Epanechnikov kernel (the supplement gives results for the uniform kernel). Finally, we evaluate the performance of the confidence intervals using several bandwidth choices. First, following the results from Section (ref), we use $\hat{h}_\mathtt{rbc}$, a data-driven version of the inference-optimal bandwidth $h_\mathtt{rbc}$. We also consider the analogous version for undersmoothed confidence intervals, denoted $\hat{h}_\mathtt{us}$ (detailed in the supplement), and the standard choice in practice, $\hat{h}_\mathtt{mse}$. Robust bias correction is implemented using $\rho=\rho^*$ according to Table (ref). All implementation details are available for R and Stata Calonico-Cattaneo-Farrell2019_JSS.
Figures (ref) and (ref) present empirical coverage probabilities for $\nu=0$ and $\nu=1$, respectively, for each evaluation point and choice of bandwidth, as a function of the sample size. Overall, we can see that robust bias correction yields close to accurate coverage, improving over undersmoothing in almost every case. Performance is highly superior at points where the functions present high curvature and also at the boundary. Performance is never worse even when the function is quite linear and optimal bandwidths are (close to) ill-defined.
We also compare confidence interval performance in terms of length in Figure (ref). We take coverage into account by looking at RBC and US confidence intervals implemented using their corresponding coverage error optimal bandwidth choices ($\hat{h}_\mathtt{rbc}$ and $\hat{h}_\mathtt{us}$, respectively), which is when they perform best in terms of coverage. We also include other valid, but non optimal choices $I_\mathtt{rbc}(\hat{h}_{\mathtt{mse}})$, $I_\mathtt{rbc}(\hat{h}_{\mathtt{us}})$. We find that RBC confidence intervals are, on average, not larger than US, and sometimes even shorter. Lastly, Figure (ref) shows the average estimated bandwidths at each point for each sample size, which behave as expected following our theory.
This paper derived higher order expansions for inference in nonparametric local polynomial regression. We provided new Edgeworth expansions and associated error in coverage probability expansions for standard and robust bias corrected methods, showing that the latter have superior coverage properties. Our results hold uniformly in the data generating process, cover derivative estimation, and allow for the uniform kernel. Using our results we developed novel bandwidth selections that target inference directly, achieving lower coverage error and/or shorter length.
Our main results measured coverage error symmetrically, but it is worth mentioning that the absolute loss function may be replaced by the “check” loss function, and thus studying the maximal coverage error $\sup_{F \in \mathscr{F}_S} \mathcal{L}( \mathbb{P}_F [ \theta_F \in I ] - (1-\alpha) )$, with $\mathcal{L}(e) =\mathcal{L}_\tau(e)=e(\tau - \mathbbm{1}\{e<0\})$, and where $\tau \in (0,1)$ encodes the researcher's weight for over- and under-coverage. Setting $\tau = 1/2$ recovers the above, symmetric measure of coverage error. Guarding more against undercoverage (a preference for conservative intervals) requires choosing a $\tau < 1/2$. For example, setting $\tau = 1/3$ encodes the belief that undercoverage is twice as bad as the same amount of overcoverage. All our results can be established for this loss function.
Finally, this paper studied the properties of confidence intervals at a fixed evaluation point $\mathsf{x}$, but it would be of theoretical and practical interest to extent our results to the case of confidence band construction. Robust bias correction has recently been used to construct valid confidence bands for local polynomial estimation Cheng-Chen2019_EJS and linear sieve estimation Cattaneo-Farrell-Feng2020_AoS. Because the underlying distributional approximations for confidence band constructions are substantially more complex, obtaining results similar to those presented herein will require substantial extension of our technical work.