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.
333,999 characters · 50 sections · 112 citation commands
On the Effect of Bias Estimation on Coverage Accuracy in Nonparametric Inference Sebastian Calonico is Assistant Professor of Economics, Department of Economics, University of Miami, Coral Gables, FL 33124 (email: [email removed]). Matias D. Cattaneo is Professor of Economics and Statistics, Department of Economics and Department of Statistics, University of Michigan, Ann Arbor, MI 48109 (email: [email removed]). Max H. Farrell is Assistant Professor of Econometrics and Statistics, Booth School of Business, University of Chicago, Chicago, IL 60637 (email: [email removed]). The second author gratefully acknowledges financial support from the National Science Foundation (SES 1357561 and SES 1459931). We thank Ivan Canay, Xu Cheng, Joachim Freyberger, Bruce Hansen, Joel Horowitz, Michael Jansson, Francesca Molinari, Ulrich M\"uller, and Andres Santos for thoughtful comments and suggestions, as well as seminar participants at Cornell, Cowles Foundation, CREST Statistics, London School of Economics, Northwestern, Ohio State University, Princeton, Toulouse School of Economics, University of Bristol, and University College London. The Associate Editor and three reviewers also provided very insightful comments that improved this manuscript.
\setcounter{page}{0}
\thispagestyle{empty}
\setcounter{page}{0} \thispagestyle{empty}
Keywords: Edgeworth expansion, coverage error, kernel methods, local polynomial regression.
\doublespacing
Nonparametric methods are widely employed in empirical work, as they provide point estimates and inference procedures that are robust to parametric misspecification bias. Kernel-based methods are commonly used to estimate densities, conditional expectations, and related functions nonparametrically in a wide variety of settings. However, these methods require specifying a bandwidth and their performance in applications crucially depends on how this tuning parameter is chosen. In particular, valid inference requires the delicate balancing act of selecting a bandwidth small enough to remove smoothing bias, yet large enough to ensure adequate precision. Tipping the scale in either direction can greatly skew results. This paper studies kernel density and local polynomial regression estimation and inference based on the popular Wald-type statistics and demonstrates (via higher-order expansions) that by coupling explicit bias correction with a novel, yet simple, Studentization, inference can be made substantially more robust to bandwidth choice, greatly easing implementability.
Perhaps the most common bandwidth selection approach is to minimize the asymptotic mean-square error (MSE) of the point estimator, and then use this bandwidth choice even when the goal is inference. So difficult is bandwidth selection perceived to be, that despite the fact that the MSE-optimal bandwidth leads to invalid confidence intervals, even asymptotically, this method is still advocated, and is the default in most popular software. Indeed, Hall-Kang2001_AoS write: “there is a growing belief that the most appropriate approach to constructing confidence regions is to estimate [the density] in a way that is optimal for pointwise accuracy\ldots. [I]t has been argued that such an approach has advantages of clarity, simplicity and easy interpretation.”
The underlying issue, as formalized below, is that bias must be removed for valid inference, and the MSE-optimal bandwidth (in particular) is “too large”, leaving a bias that is still first order. Two main methods have been proposed to address this: undersmoothing and explicit bias correction. We seek to compare these two, and offer concrete ways to better implement the latter. Undersmoothing amounts to choosing a bandwidth smaller than would be optimal for point estimation, then arguing that the bias is smaller than the variability of the estimator asymptotically, leading to valid distributional approximations and confidence intervals. In practice this method often involves simply shrinking the MSE-optimal bandwidth by an ad-hoc amount. The second approach is to bias correct the estimator with the explicit goal of removing the bias that caused the invalidity of the inference procedure in the first place.
It has long been believed that undersmoothing is preferable for two reasons. First, theoretical studies showed inferior asymptotic coverage properties of bias-corrected confidence intervals. The pivotal work was done by Hall1992_AoS_density, and has been relied upon since. Second, implementation of bias correction is perceived as more complex because a second (usually different) bandwidth is required, deterring practitioners. However, we show theoretically that bias correction is always as good as undersmoothing, and better in many practically relevant cases, if the new standard errors that we derive are used. Further, our findings have important implications for empirical work because the resulting confidence intervals are more robust to bandwidth choice, including to the bandwidth used for bias estimation. Indeed, the two bandwidths may be set equal, a simple and automatic choice that performs well in practice and is optimal in certain objective senses.
Our proposed robust bias correction method delivers valid confidence intervals (and related inference procedures) even when using the MSE-optimal bandwidth for the original point estimator, the most popular approach in practice. Moreover, we show that at interior points, when using second-order kernels or local linear regressions, the coverage error of such intervals vanishes at the best possible rate. (Throughout, the notion of “optimal” or “best” rate is defined as the fastest achievable coverage error decay for a fixed kernel order or polynomial degree; and is also different from optimizing point estimation.) When higher-order kernels are used, or boundary points are considered, we find that the corresponding MSE-optimal bandwidth leads to asymptotically valid intervals, but with suboptimal coverage error decay rates, and must be shrunk (sometimes considerably) for better inference.
Heuristically, employing the MSE-optimal bandwidth for the original point estimator, prior to bias correction, is like undersmoothing the bias-corrected point estimator, though the latter estimator employs a possibly random, $n$-varying kernel, and requires a different Studentization scheme. It follows that the conventional MSE-optimal bandwidth commonly used in practice need not be optimal, even after robust bias correction, when the goal is inference. Thus, we present new coverage error optimal bandwidths and a fully data-driven direct plug-in implementation thereof, for use in applications. In addition, we study the important related issue of asymptotic length of the new confidence intervals.
Our comparisons of undersmoothing and bias correction are based on Edgeworth expansions for density estimation and local polynomial regression, allowing for different levels of smoothness of the unknown functions. We prove that explicit bias correction, coupled with our proposed standard errors, yields confidence intervals with coverage that is as accurate, or better, than undersmoothing (or, equivalently, yields dual hypothesis tests with lower error in rejection probability). Loosely speaking, this improvement is possible because explicit bias correction can remove more bias than undersmoothing, while our proposed standard errors capture not only the variability of the original estimator but also the additional variability from bias correction. To be more specific, our robust bias correction approach yields higher-order refinements whenever additional smoothness is available, and is asymptotically equivalent to the best undersmoothing procedure when no additional smoothness is available.
Our findings contrast with well established recommendations: Hall1992_AoS_density used Edgeworth expansions to show that undersmoothing produces more accurate intervals than explicit bias correction in the density case and Neumann1997_Statistics repeated this finding for kernel regression. The key distinction is that their expansions, while imposing the same levels of smoothness as we do, crucially relied on the assumption that the bias correction was first-order negligible, essentially forcing bias correction to remove less bias than undersmoothing. In contrast, we allow the bias estimator to potentially have a first order impact, an alternative asymptotic experiment designed to more closely mimic the finite-sample behavior of bias correction. Therefore, our results formally show that whenever additional smoothness is available to characterize leading bias terms, as is usually the case in practice where MSE-optimal bandwidth are employed, our robust bias correction approach yields higher-order improvements relative to standard undersmoothing.
Our standard error formulas are based on fixed-$n$ calculations, as opposed to asymptotics, which also turns out to be important. We show that using asymptotic variance formulas can introduce further errors in coverage probability, with particularly negative consequences at boundary points. This turns out to be at the heart of the “quite unexpected” conclusion found by Chen-Qin2002_SJS that local polynomial based confidence intervals are not boundary-adaptive in coverage error: we prove that this is not the case with proper Studentization. Thus, as a by-product of our main theoretical work, we establish higher-order boundary carpentry of local polynomial based confidence intervals that use a fixed-$n$ standard error formula, a result that is of independent (but related) interest.
This paper is connected to the well-established literature on nonparametric smoothing, see Wand-Jones1995_book, Fan-Gijbels1996_book, Horowitz2009_book, and Ruppert-Wand-Carroll2009_book for reviews. For more recent work on bias and related issues in nonparametric inference, see Hall-Horowitz2013_AoS, Calonico-Cattaneo-Titiunik2014_Ecma, Armstrong-Kolesar2015_minimax, Schennach2015_bias, and references therein. We also contribute to the literature on Edgeworth expansions, which have been used both in parametric and, less frequently, nonparametric contexts; see, e.g., Bhattacharya-Rao1976_book and Hall1992_book. Fixed-$n$ versus asymptotic-based Studentization has also captured some recent interest in other contexts, e.g., Mykland-Zhang2015_WP. Finally, see Calonico-Cattaneo-Farrell2016_EE-RD for uniformly valid Edgeworth expansions and optimal inference.
The paper proceeds as follows. Section (ref) studies density estimation at interior points and states the main results on error in coverage probability and its relationship to bias reduction and underlying smoothness, as well as discussing bandwidth choice and interval length. Section (ref) then studies local polynomial estimation at interior and boundary points. Practical guidance is explicitly discussed in Sections (ref) and (ref), respectively; all methods are available in {\sf R} and {\tt STATA} via the {\tt nprobust} package, see Calonico-Cattaneo-Farrell2017_nprobust. Section (ref) summarizes the results of a Monte Carlo study, and Section (ref) concludes. Some technical details, all proofs, and additional simulation evidence are collected in a lengthy online supplement.
We first present our main ideas and conclusions for inference on the density at an interior point, as this requires relatively little notation. The data are assumed to obey the following.
The parameter of interest is $f(x)$ for a fixed scalar point $x$ in the interior of the support. (In the supplemental appendix we discuss how our results extend naturally to multivariate $X_i$ and derivative estimation.) The classical kernel-based estimator of $f(x)$ is
for a kernel function $K$ that integrates to $1$ and positive bandwidth $h \to 0$ as $n \to \infty$. The choice of $h$ can be delicate, and our work is motivated in part by the standard empirical practice of employing the MSE-optimal bandwidth choice for $\hat{f}(x)$ when conducting inference.
In this vein, let us suppose for the moment that $K$ is a kernel of order $\mathcal{k}$, where $\mathcal{k} \leq S$ so that the MSE-optimal bandwidth can be characterized. The bias is then given by
where $f^{(\mathcal{k})}(x):=\partial^\mathcal{k} f(x)/\partial x^\mathcal{k}$ and $\mu_{K,\mathcal{k}} = \int u^\mathcal{k} K(u) du / \mathcal{k}!$. Computing the variance gives
which is non-asymptotic: $n$ and $h$ are fixed in this calculation. Using other, first-order valid approximations, e.g.\ $(nh) \mathbb{V}[\hat{f}(x)] \approx f(x) \int K(u)^2 du$, will have finite sample consequences that manifest as additional terms in the Edgeworth expansions. In fact, Section (ref) shows that using an asymptotic variance for local polynomial regression removes automatic coverage-error boundary adaptivity.
Together, the prior two displays are used to characterize the MSE-optimal bandwidth, $h^*_\mathtt{mse}\propto n^{-1/(1 + 2 \mathcal{k})}$. However, using this bandwidth leaves a bias that is too large, relative to the variance, to conduct valid inference for $f(x)$. To address this important practical problem, researchers must either undersmooth the point estimator (i.e., construct $\hat{f}(x)$ with a bandwidth smaller than $h^*_\mathtt{mse}$) or bias-correct the point estimator (i.e., subtract an estimate of the leading bias). Thus, the question we seek to answer is this: if the bias is given by (ref), is one better off estimating the leading bias (explicit bias correction) or choosing $h$ small enough to render the bias negligible (undersmoothing) when forming nonparametric confidence intervals?
To answer this question, and to motivate our new robust approach, we first detail the bias correction and variance estimators. Explicit bias correction estimates the leading term of Eqn.\ (ref), denoted by $B_f$, using a kernel estimator of $f^{(\mathcal{k})}(x)$, defined as:
for a kernel $L(\cdot)$ of order $\ell$ and a bandwidth $b \to 0$ as $n \to \infty$. Importantly, $\hat{B}_f$ takes this form for any $\mathcal{k}$ and $S$, even if (ref) fails; see Sections (ref) and (ref) for discussion. Conventional Studentized statistics based on undersmoothing and explicit bias correction are, respectively, \[T_\mathtt{us}(x) = \frac{\sqrt{nh}\bigl(\hat{f}(x) - f(x)\bigr)}{\hat{\sigma}_\mathtt{us}} \qquad \text{ and } \qquad T_\mathtt{bc}(x) = \frac{\sqrt{nh}\bigl(\hat{f}(x) - \hat{B}_f - f(x)\bigr)}{\hat{\sigma}_\mathtt{us}},\] where $\hat{\sigma}_\mathtt{us}^2 := \hat{\mathbb{V}}[ \hat{f}(x)]$ is the natural estimator of the variance of $\hat{f}(x)$ which only replaces the two expectations in (ref) with sample averages, thus maintaining the nonasymptotic spirit. These are the two statistics compared in the influential paper of Hall1992_AoS_density, under the same assumption imposed herein.
From the form of these statistics, two points are already clear. First, the numerator of $T_\mathtt{us}$ relies on choosing $h$ vanishing fast enough so that the bias is asymptotically negligible after scaling, whereas $T_\mathtt{bc}$ allows for slower decay by virtue of the manual estimation of the leading bias. Second, $T_\mathtt{bc}$ requires that the variance of $h^\mathcal{k} \hat{f}^{(\mathcal{k})}(x) \mu_{K,\mathcal{k}}$ be first-order asymptotically negligible: $\hat{\sigma}_\mathtt{us}$ in the denominator only accounts for the variance of the main estimate, but $\hat{f}^{(\mathcal{k})}(x)$, being a kernel-based estimator, naturally has a variance controlled by its bandwidth. That is, even though $\hat{\sigma}_\mathtt{us}^2$ is based on a fixed-$n$ calculation, the variance of the numerator of $T_\mathtt{bc}$ only coincides with the denominator asymptotically. Under this regime, Hall1992_AoS_density showed that the bias reduction achieved in $T_\mathtt{bc}$ is too expensive in terms of noise and that undersmoothing dominates explicit bias correction for coverage error.
We argue that there need not be such a “mismatch” between the numerator of the bias-corrected statistic and the Studentization, and thus consider a third option corresponding to the idea of capturing the finite sample variability of $\hat{f}^{(\mathcal{k})}(x)$ directly. To do so, note that we may write, after setting $\rho=h/b$,
We then define the collective variance of the density estimate and the bias correction as $\sigma_\mathtt{rbc}^2 = (nh) \mathbb{V}[\hat{f}(x) - \hat{B}_f]$, exactly as in Eqn.\ (ref), but with $M(\cdot)$ in place of $K(\cdot)$, and its estimator $\hat{\sigma}_\mathtt{rbc}^2$ exactly as $\hat{\sigma}_\mathtt{us}^2$. Therefore, our proposed robust bias corrected inference approach is based on \[T_\mathtt{rbc} = \frac{\sqrt{nh}\bigl(\hat{f}(x) - h^\mathcal{k} \hat{f}^{(\mathcal{k})}(x) \mu_{K,\mathcal{k}} - f(x)\bigr)}{ \hat{\sigma}_\mathtt{rbc} }.\] That is, our proposed standard errors are based on a fixed-$n$ calculation that captures the variability of both $\hat{f}(x)$ and $\hat{f}^{(\mathcal{k})}(x)$, and their covariance. As shown in Section (ref), the case of local polynomial regression is analogous, but notationally more complicated.
The quantity $\rho=h/b$ is key. If $\rho \to 0$, then the second term of $M$ is dominated by the first, i.e.\ the bias correction is first-order negligible. In this case, $\sigma_\mathtt{us}^2$ and $\sigma_\mathtt{rbc}^2$ (and their estimators) will be first-order, but not higher-order, equivalent. This is exactly the sense in which traditional bias correction relies on an asymptotic variance, instead of a fixed-$n$ one, and pays the price in coverage error. To more accurately capture finite sample behavior of bias correction we allow $\rho$ to converge to any (nonnegative) finite limit, allowing (but not requiring) the bias correction to be first-order important, unlike prior work. We show that doing so yields more accurate confidence intervals (i.e., higher-order corrections).
We first present generic Edgeworth expansions for all three procedures (undersmoothing, traditional bias correction, and robust bias correction), which are agnostic regarding the level of available smoothness (controlled by $S$ in Assumption (ref)). To be specific, we give higher-order expansions of the error in coverage probability of the following $(1-\alpha)\%$ confidence intervals based on Normal approximations for the statistics $T_\mathtt{us}$, $T_\mathtt{bc}$, and $T_\mathtt{rbc}$:
where $z_{\alpha}$ is the upper $\alpha$-percentile of the Gaussian distribution. Here and in the sequel we omit the point of evaluation $x$ for simplicity. Equivalently, our results can characterize the error in rejection probability of the corresponding hypothesis tests. In subsequent sections, we give specific results under different smoothness assumptions and make direct comparisons of the methods.
We require the following standard conditions on the kernels $K$ and $L$.
The boundary conditions are needed for the derivative estimation inherent in bias correction, even if $x$ is an interior point, and are satisfied if the support of $f$ is the whole real line. Higher order results also require a standard $n$-varying Cram\'er's condition, given in the supplement to conserve space (see Section S.I.3). Altogether, our assumptions are identical to those of Hall1991_Statistics,Hall1992_AoS_density.
To state the results some notation is required. First, let the (scaled) biases of the density estimator and the bias-corrected estimator be $\eta_\mathtt{us} = \sqrt{nh}(\mathbb{E}[\hat{f}] - f)$ and $\eta_\mathtt{bc} = \sqrt{nh}(\mathbb{E}[\hat{f} - \hat{B}_f] - f)$. Next, let $\phi(z)$ be the standard Normal density, and for any kernel $K$ define
where $\vartheta_{K,k} = \int K(u)^k du$. All that is conceptually important is that these functions are known, odd polynomials in $z$ with coefficients that depend only on the kernel, and not on the sample or data generating process. Our main theoretical result for density estimation is the following.
This result leaves the scaled biases $\eta_\mathtt{us}$ and $\eta_\mathtt{bc}$ generic, which is useful when considering different levels of smoothness $S$, the choices of $\mathcal{k}$ and $\ell$, and in comparing to local polynomial results. In the next subsection, we make these quantities more precise and compare them, paying particular attention to the role of the underlying smoothness assumed.
At present, the most visually obvious feature of this result is that all the error terms are of the same form, except for the notable presence of $\rho^{1+\mathcal{k}}(\Omega_1 + \rho^\mathcal{k} \Omega_2)$ in part (b). These are the leading terms of $\sigma_\mathtt{rbc}^2 / \sigma_\mathtt{us}^2 - 1$, consisting of the covariance of $\hat{f}$ and $\hat{B}_f$ (denoted by $\Omega_1$) and the variance of $\hat{B}_f$ (denoted by $\Omega_2$), and are entirely due to the “mismatch” in the Studentization of $T_\mathtt{bc}$. Hall1992_AoS_density showed how these terms prevent bias correction from performing as well as undersmoothing in terms of coverage. In essence, the potential for improved bias properties do not translate into improved inference because the variance is not well-controlled: in any finite sample, $\hat{B}_f$ would inject variability (i.e., $\rho=h/b>0$ for each $n$) and thus $\rho \to 0$ may not be a good approximation. Our new Studentization does not simply remove these leading $\rho$ terms; the entire sequence is absent. As explained below, allowing for $\bar{\rho}=\infty$ can not reduce bias, but will inflate variance; hence restricting to $\bar{\rho}<\infty$ capitalizes fully on the improvements from bias correction.
Theorem (ref) makes no explicit assumption about smoothness beyond the requirement that the scaled biases vanish asymptotically. The fact that the error terms in parts (a) and (c) of Theorem (ref) take the same form implies that comparing coverage error amounts to comparing bias, for which the smoothness $S$ and the kernel orders $\mathcal{k}$ and $\ell$ are crucial. We now make the biases $\eta_\mathtt{us}$ and $\eta_\mathtt{bc}$ concrete and show how coverage is affected.
For $I_\mathtt{us}$, two cases emerge: (a) enough derivatives exist to allow characterization of the MSE-optimal bandwidth ($\mathcal{k} \leq S$); and (b) no such smoothness is available ($\mathcal{k} > S$), in which case the leading term of Eqn.\ (ref) is exactly zero and the bias depends on the unknown H\"older constant. These two cases lead to the following results.
The first result is most directly comparable to Hall1992_AoS_density, and many other past papers, which typically take as a starting point that the MSE-optimal bandwidth can be characterized. This shows that $T_\mathtt{us}$ must be undersmoothed, in the sense that the MSE-optimal bandwidth is “too large” for valid inference. In fact, we know that $I_\mathtt{us}(h^*_\mathtt{mse})$ will asymptotically undercover because $T_\mathtt{us}(h^*_\mathtt{mse}) \to_d \mathscr{N}((2\mathcal{k})^{-1/2},1)$ (see the supplement). Instead, the optimal $h$ for coverage error, which can be characterized and estimated, is equivalent in rates to balancing variance against bias, not squared bias as in MSE. Part (b) shows that a faster rate of coverage error decay can be obtained by taking a sufficiently high order kernel, relative to the level of smoothness $S$, at the expense of feasible bandwidth selection.
Turning to robust bias correction, characterization of $\eta_\mathtt{bc}$ is more complex as it has two pieces: the second-order bias of the original point estimator, and the bias of the bias estimator itself. The former is the $o(h^\mathcal{k})$ term of Eqn.\ (ref) and is not the target of explicit bias correction; it depends either on higher derivatives, if they are available, or on the H\"older condition otherwise. To be precise, if $\mathcal{k} \leq S-2$, this term is $[h^{\mathcal{k} + 2} + o(1)] f^{(\mathcal{k} + 2)} \mu_{K_\mathtt{bc}, \mathcal{k} + 2}$, while otherwise is known only to be $O(h^{S + \varsigma})$. Importantly, the bandwidth $b$ and order $\ell$ do not matter here, and bias reduction beyond $O(\min\{h^{\mathcal{k} + 2},h^{S + \varsigma}\})$ is not possible; there is thus little or no loss in fixing $\ell=2$, which we assume from now on to simplify notation.
The bias of the bias estimator also depends on the smoothness available: if enough smoothness is available the corresponding bias term can be characterized, otherwise only its order will be known. To be specific, when smoothness is not binding ($\mathcal{k} \leq S-2$), arguably the most practically-relevant case, the leading term of $\mathbb{E}[\hat{B}_f] - B_f$ will be $h^\mathcal{k} b^2 f^{(\mathcal{k} + 2)} \mu_{K,\mathcal{k}} \mu_{L,2}$. Smoothness can be exhausted in two ways, either by the point estimate itself ($\mathcal{k} > S$) or by the bias estimation ($S-1 \leq \mathcal{k} \leq S$), and these two cases yield $O(h^\mathcal{k} b^{S - \mathcal{k}})$ and $O(h^\mathcal{k} b^{S + \varsigma - \mathcal{k}})$, respectively, which are slightly different in how they depend on the total H\"older smoothness assumed. (Complete details are in the supplement.) Note that regardless of the value of $\mathcal{k}$, we set $\hat{B}_f = h^\mathcal{k} \hat{f}^{(\mathcal{k})} \mu_{K,\mathcal{k}}$, even if $\mathcal{k} > S$ and $B_f \equiv 0$.
With these calculations for $\eta_\mathtt{bc}$, we have the following result.
Part (a) is the most empirically-relevant setting, which reflects the idea that researchers first select a kernel order, then conduct inference based on that choice, taking the unknown smoothness to be nonbinding. The most notable feature of this result, beyond the formalization of the coverage improvement, is that the coverage error terms share the same structure as those of Corollary (ref), with $\mathcal{k}$ replaced by $\mathcal{k}+2$, and represent the same conceptual objects. By virtue of our new Studentization, the leading variance remains order $(n h)^{-1}$ and the problematic correlation terms are absent. We explicitly discuss the advantages of robust bias correction relative to undersmoothing in the following section.
Part (a) also argues for a bounded, positive $\rho$. First, because bias reduction beyond $O(h^{\mathcal{k}+2})$ is not possible, $\rho \to \infty$ will only inflate the variance. On the other hand, $\bar{\rho} = 0$ requires a delicate choice of $b$ and $\ell > 2$, else the second bias term dominates $\eta_\mathtt{bc}$, and the full power of the variance correction is not exploited; that is, more bias may be removed without inflating the variance rate. Hall1992_AoS_density remarked that if $\mathbb{E}[\hat{f}] - f - B_f$ is (part of) the leading bias term, then “explicit bias correction [\ldots] is even less attractive relative to undersmoothing.” We show, on the contrary, that with our proposed Studentization, it is optimal that $\mathbb{E}[\hat{f}] - f - B_f$ is part of the dominant bias term.
Finally, in both Corollaries above the best possible coverage error decay rate (for a given $S$) is attained by exhausting all available smoothness. This would also yield point estimators attaining the bound of Stone1982_AoS; robust bias correction can not evade such bounds, of course. In both Corollaries, coverage is improved relative to part (a), but the constants and optimal bandwidths can not be quantified. For robust bias correction, Corollary (ref) shows that to obtain the best rate in part (b) the unknown $f^{(\mathcal{k})}$ must be consistently estimated and $\rho$ must be bounded and positive, while in part (c), bias estimation merely adds noise, but this noise is fully accounted for by our new Studentization, as long as $\rho \to 0$ ($b \not\to 0$ is allowed).
We now employ Corollaries (ref) and (ref) to directly compare nonparametric inference based on undersmoothing and robust bias correction. To simplify the discussion we focus on three concrete cases, which illustrate how the comparisons depend on the available smoothness and kernel order; the messages generalize to any $S$ and/or $\mathcal{k}$. For this discussion we let $\mathcal{k}_\mathtt{us}$ and $\mathcal{k}_\mathtt{bc}$ be the kernel orders used for point estimation in $I_\mathtt{us}$ and $I_\mathtt{rbc}$, respectively, and restrict attention to sequences $h\to0$ where both confidence intervals are first-order valid, even though robust bias correction allows for a broader bandwidth range. Finally, we set $\ell=2$ and $\bar{\rho} \in (0,\infty)$ based on the above discussion.
For the first case, assume that $f$ is twice continuously differentiable ($S=2$) and both methods use second order kernels ($\mathcal{k}_\mathtt{us} = \mathcal{k}_\mathtt{bc} = \ell = 2$). In this case, both methods target the same bias. The coverage errors for $I_\mathtt{us}$ and $I_\mathtt{rbc}$ then follow directly from Corollaries (ref)(a) and (ref)(b) upon plugging in these kernel orders, yielding \[ \bigl| \mathbb{P}[f \in I_\mathtt{us}] - (1 - \alpha) \bigr| \asymp \frac{1}{n h} + nh^5 + h^2 \quad \text{and} \quad \bigl| \mathbb{P}[f \in I_\mathtt{rbc}] - (1 - \alpha) \bigr| \asymp \frac{1}{n h} + n h^{5+2\varsigma} + h^{2 +\varsigma} .\] Because $h \to 0$ and $\bar{\rho} \in (0,\infty)$, the coverage error of $I_\mathtt{rbc}$ vanishes more rapidly by virtue of the bias correction. A higher order kernel ($\mathcal{k}_\mathtt{us} > 2$) would yield this rate for $I_\mathtt{us}$.
Second, suppose that the density is four-times continuously differentiable ($S=4$) but second order kernels are maintained. The relevant results are now Corollaries (ref)(a) and (ref)(a). Both methods continue to target the same leading bias, but now the additional smoothness available allows precise characterization of the improvement shown above, and we have \[ \bigl| \mathbb{P}[f \in I_\mathtt{us}] - (1 - \alpha) \bigr| \asymp \frac{1}{n h} + nh^5 + h^2 \quad \text{ and } \quad \bigl| \mathbb{P}[f \in I_\mathtt{rbc}] - (1 - \alpha) \bigr| \asymp \frac{1}{n h} + n h^9 + h^4 .\] This case is perhaps the most empirically relevant one, where researchers first choose the order of the kernel (here, second order) and then conduct/optimize inference based on that choice. Indeed, for this case optimal bandwidth choices can be derived (Section (ref)).
Finally, maintain $S=4$ but suppose that undersmoothing is based on a fourth-order kernel while bias correction continues to use two second-order kernels ($\mathcal{k}_\mathtt{us} = 4$, $\mathcal{k}_\mathtt{bc} = \ell = 2$). This is the exact example given by Hall1992_AoS_density. Now the two methods target different biases, but utilize the same amount of smoothness. In this case, the relevant results are again Corollaries (ref)(a) and (ref)(a), now with $\mathcal{k}=4$ and $\mathcal{k}=2$, respectively. The two methods have the same coverage error decay rate: \[ \bigl| \mathbb{P}[f \in I_\mathtt{us}] - (1 - \alpha) \bigr| \asymp \bigl| \mathbb{P}[f \in I_\mathtt{rbc}] - (1 - \alpha) \bigr| \asymp \frac{1}{n h} + n h^9 + h^4 .\] Indeed, more can be said: with the notation of Eqn.\ (ref), the difference between $T_\mathtt{us}$ and $T_\mathtt{rbc}$ is the change in “kernel” from $K$ to $M$, and since $\mathcal{k}_\mathtt{bc} + \ell = \mathcal{k}_\mathtt{us}$, the two kernels are the same order. ($M$ acts as a $n$-varying, higher-order kernel for bias, but may not strictly fit the definition, as explored in the supplement.) This tight link between undersmoothing and robust bias correction does not carry over straightforwardly to local polynomial regression, as we discuss in more detail in Section (ref).
In the context of this final example, it is worth revisiting traditional bias correction. The fact that undersmoothing targets a different, and asymptotically smaller, bias than does explicit bias correction, coupled with the requirement that $\rho \to 0$, implicitly constrains bias correction to remove less bias than undersmoothing. This is necessary for traditional bias correction, but on the contrary, robust bias correction attains the same coverage error decay rate as undersmoothing under the same assumptions.
In sum, these examples show that under identical assumptions, bias correction is not inferior to undersmoothing and if any additional smoothness is available, can yield improved coverage error. These results are confirmed in our simulations.
The prior sections established that robust bias correction can equal, or outperform, undersmoothing for inference. We now show how the method can be implemented to deliver these results in applications. We mimic typical empirical practice where researchers first choose the order of the kernel, then conduct/optimize inference based on that choice. Therefore, we assume the smoothness is unknown but taken to be large and work within Corollary (ref)(a), that is, viewing $\mathcal{k} \leq S-2$ and $\ell=2$ as fixed and $\rho$ bounded and positive. This setup allows characterization of the coverage error optimal bandwidth for robust bias correction.
We can use this result to give concrete methodological recommendations. At the end of this section we discuss the important issue of interval length. Construction of the interval $I_\mathtt{rbc}$ from Eqn.\ (ref) requires bandwidths $h$ and $b$ and kernels $K$ and $L$. Given these choices, the point estimate, bias correction, and variance estimators are then readily computable from data using the formulas above. For the kernels $K$ and $L$, we recommend either second order minimum variance (to minimize interval length) or MSE-optimal kernels Gasser-Muller-Mammitzsch1985_JRSSB.
The bandwidth selections are more important in applications. For the bandwidth $h$, Corollary (ref)(a) shows that the MSE-optimal choice $h^*_\mathtt{mse}$ will deliver valid inference, but will be suboptimal in general (Corollary (ref)). From a practical point of view, the robust bias corrected interval $I_\mathtt{rbc}(h)$ is attractive because it allows for the MSE-optimal bandwidth and kernel, and hence is based on the MSE-optimal point estimate, while using the same effective sample for both point estimation and inference. Interestingly, although $I_\mathtt{rbc}(h^*_\mathtt{mse})$ is always valid, its coverage error decays as $n^{-\min\{4,\mathcal{k}+2\}/(1+2\mathcal{k})}$ and is thus rate optimal only for second order kernels ($\mathcal{k}=2$), while otherwise being suboptimal, with a rate that is slower the larger is the order $\mathcal{k}$.
Corollary (ref) gives the coverage error optimal bandwidth, $h^*_\mathtt{rbc}$, which can be implemented using a simple direct plug-in (DPI) rule: $\hat{h}_{\mathtt{dpi}} = \hat{H}_{\mathtt{dpi}} \; n^{-1/(\mathcal{k}+3)}$, where $\hat{H}_{\mathtt{dpi}}$ is a plug-in estimate of $H^*_\mathtt{rbc}$ formed by replacing the unknown $f^{(\mathcal{k}+2)}$ with a pilot estimate (e.g., a consistent nonparametric estimator based on the appropriate MSE-optimal bandwidth). In the supplement we give precise implementation details, as well as an alternative rule-of-thumb bandwidth selector based on rescaling already available data-driven MSE-optimal choices.
For the bandwidth $b$, a simple choice is $b = h$, or, equivalently, $\rho = 1$. We show in the supplement that setting $\rho=1$ has good theoretical properties, minimizing interval length of $I_\mathtt{rbc}$ or the MSE of $\hat{f} - \hat{B}_f$, depending on the conditions imposed. In our numerical work, we found that $\rho=1$ performed well. As a result, from the practitioner's point of view, the choice of $b$ (or $\rho$) is completely automatic, leaving only one bandwidth to select.
An extensive simulation study, reported in the supplement, illustrates our findings and explores the numerical performance of these choices. We find that coverage of $I_\mathtt{rbc}$ is robust to both $h$ and $\rho$ and that our data-driven bandwidth selectors work well in practice, but we note that estimating bandwidths may have higher-order implications Hall-Kang2001_AoS.
Finally, an important issue in applications is whether the good coverage properties of $I_\mathtt{rbc}$ come at the expense of increased interval length. When coverage is asymptotically correct, Corollaries (ref) and (ref) show that $I_\mathtt{rbc}$ can accommodate (and will optimally employ) a larger bandwidth (i.e.\ $h \to 0$ more slowly), and hence $I_\mathtt{rbc}$ will have shorter average length in large samples than $I_\mathtt{us}$. Our simulation study (see below and the supplement) gives the same conclusion.
We study a plug-in bias correction method, but there are alternatives. In particular, as pointed out by a reviewer, a leading alternative is the generalized jackknife method of Schucany-Sommers1977_JASA (see Cattaneo-Crump-Jansson_2013_JASA for an application to kernel-based semiparametric inference and for related references). We will briefly summarize this approach and show a tight connection to our results, restricting to second-order kernels and $S\geq 2$ only for simplicity.
The generalized jackknife estimator is $\hat{f}_{\mathtt{GJ},R} := ( \hat{f}_1 - R \hat{f}_2 ) / (1 - R)$, where $\hat{f}_1$ and $\hat{f}_2$ are two initial kernel density estimators, with possibly different bandwidths ($h_1,h_2$) and second-order kernels ($K_1,K_2$). From Eqn.\ (ref), the bias of $\hat{f}_{\mathtt{GJ},R}$ is $(1-R)^{-1} f^{(2)} \left( h_1^2 \mu_{K_1,2} - R h_2^2 \mu_{K_2,2} \right) + o(h_1^2 + h_2^2)$, whence choosing $R = (h_1^2 \mu_{K_1,2} )/( h_2^2 \mu_{K_2,2})$ renders the leading bias term exactly zero. Further, if $S \geq 4$, $\hat{f}_{\mathtt{GJ},R}$ has bias $O(h_1^4 + h_2^4)$; behaving as a point estimator with $\mathcal{k}=4$. To connect this approach to ours, observe that with this choice of $R$ and $\tilde{\rho} = h_1 / h_2$, \[\hat{f}_{\mathtt{GJ},R} = \frac{1}{n h_1} \sum_{i=1}^n \tilde{M}\left(\frac{X_i - x}{h_1} \right) , \quad \tilde{M}(u) = K_1(u) - \tilde{\rho}^{1+2} \left\{ \frac{K_2(\tilde{\rho}u) - \tilde{\rho}^{-1} K_1(u) }{\mu_{K_2,2}(1-R)} \right\} \mu_{K_1, 2},\] exactly matching Eqn.\ (ref); alternatively, write $\hat{f}_{\mathtt{GJ},R} = \hat{f}_1 - h_1^2 \tilde{f}^{(2)} \mu_{K_1,2}$, where \[\tilde{f}^{(2)} = \frac{1}{nh_2^{1+2}}\sum_{i=1}^n \tilde{L}\left( \frac{X_i - x}{h_2} \right), \qquad \tilde{L}(u) = \frac{ K_2(u) - \tilde{\rho}^{-1} K_1( \tilde{\rho}^{-1} u) }{\mu_{K_2,2}(1-R)}, \] is a derivative estimator. Therefore, we can view $\hat{f}_{\mathtt{GJ},R}$ as a specific kernel $M$ or a specific derivative estimator, and all our results directly apply to $\hat{f}_{\mathtt{GJ},R}$; hence our paper offers a new way of conducting inference (new Studentization) for this case as well. Though we omit the details to conserve space, this is equally true for local polynomial regression (Section (ref)).
More generally, our main ideas and generic results apply to many other bias correction methods. For a second example, Singh1977_AoS also proposed a plug-in bias estimator, but without using the derivative of a kernel. Our results cover this approach as well; see the supplement for further details and references. The key, common message in all cases is that to improve inference one must account for the additional variability introduced by any bias correction method (i.e., to avoid the mismatch present in $T_\mathtt{bc}$).
This section studies local polynomial regression Ruppert-Wand1994_AoS,Fan-Gijbels1996_book, and has two principal aims. First, we show that the conclusions from the density case, and their implications for practice, carry over to odd-degree local polynomials. Second, we show that with proper fixed-$n$ Studentization, coverage error adapts to boundary points. We focus on what is novel relative to the density, chiefly variance estimation and boundary points. For interior points, the implications for coverage error, bandwidth selection, and interval length are all analogous to the density case, and we will not retread those conclusions.
To be specific, throughout this section we focus on the case where the smoothness is large relative to the local polynomial degree $p$, which is arguably the most relevant case in practice. The results and discussion in Sections (ref) and (ref) carry over, essentially upon changing $\mathcal{k}$ to $p+1$ and $\ell$ to $q-p$ (or $q-p+1$ for interior points with $q$ even). Similarly, but with increased notational burden, the conclusions of Section (ref) also remain true. The present results also extend to multivariate data and derivative estimation.
To begin, we define the regression estimator, its bias, and the bias correction. Given a random sample $\{(Y_i,X_i):1\leq i\leq n\}$, the local polynomial estimator of $m(x)=\mathbb{E}[Y_i|X_i=x]$, temporarily making explicit the evaluation point, is \[\hat{m}(x) = \mathbf{e}_0' \boldsymbol{\hat{\beta}}_p, \qquad\quad \boldsymbol{\hat{\beta}}_p = \operatorname*{arg\,min}_{\boldsymbol{b} \in \mathbb{R}^{p+1}} \sum_{i=1}^n ( Y_i - \mathbf{r}_p(X_i - x)'\boldsymbol{b})^2 K \left( \frac{X_i - x}{h}\right), \] where, for an integer $p \geq 1$, $\mathbf{e}_0$ is the $(p+1)$-vector with a one in the first position and zeros in the rest, and $\mathbf{r}_p(u) = (1, u, u^2, \ldots, u^p)'$. We restrict attention to $p$ odd, as is standard, though the qualifier may be omitted. We define $\mathbf{Y} = (Y_1, \cdots, Y_n)'$, $\mathbf{R}_p = [ \mathbf{r}_p( (X_1 - x) / h), \cdots, \mathbf{r}_p( (X_n - x) / h)]'$, $\mathbf{W}_p = \operatorname*{diag}(h^{-1} K((X_i - x)/h): i = 1, \ldots, n)$, and $\boldsymbol{\Gamma}_p = \mathbf{R}_p' \mathbf{W}_p \mathbf{R}_p/n$ (here $\operatorname*{diag}(a_i:i = 1, \ldots, n)$ denotes the $n\times n$ diagonal matrix constructed using $a_1, a_2, \cdots, a_n$). Then, reverting back to omitting the argument $x$, the local polynomial estimator is $\hat{m} = \mathbf{e}_0'\boldsymbol{\Gamma}_p^{-1} \mathbf{R}_p' \mathbf{W}_p \mathbf{Y} / n$.
Under regularity conditions below, the conditional bias satisfies
where $\boldsymbol{\Lambda}_p = \mathbf{R}_p' \mathbf{W}_p [ ((X_1 - x)/h)^{p+1}, \cdots, ((X_n - x)/h)^{p+1}]'/n.$ Here, the quantity $\mathbf{e}_0' \boldsymbol{\Gamma}_p^{-1} \boldsymbol{\Lambda}_p / (p+1)!$ is random, unlike in the density case (c.f. (ref)), but it is known and bounded in probability. Following Fan-Gijbels1996_book, we will estimate $m^{(p+1)}$ in ((ref)) using a second local polynomial regression, of degree $q > p$ (even or odd), based on a kernel $L$ and bandwidth $b$. Thus, $\mathbf{r}_q(u)$, $\mathbf{R}_q$, $\mathbf{W}_q$, and $\boldsymbol{\Gamma}_q$ are defined as above, but substituting $q$, $L$, and $b$ in place of $p$, $K$, and $h$, respectively. Denote by $\mathbf{e}_{p+1}$ the $(q+1)$-vector with one in the $p+2$ position, and zeros in the rest. Then we estimate the bias with \[\hat{B}_m = h^{p+1} \hat{m}^{(p+1)} \frac{1}{(p+1)!} \mathbf{e}_0' \boldsymbol{\Gamma}_p^{-1} \boldsymbol{\Lambda}_p, \qquad \hat{m}^{(p+1)} = b^{-p-1} (p+1)!\mathbf{e}_{p+1}'\boldsymbol{\Gamma}_q^{-1} \mathbf{R}_q' \mathbf{W}_q \mathbf{Y} / n.\] Exactly as in the density case, $\hat{B}_m$ introduces variance that is controlled by $\rho$ and will be captured by robust bias correction.
The Studentizations in the density case were based on fixed-$n$ expectations, and we will show that retaining this is crucial for local polynomials. The fixed-$n$ versus asymptotic distinction is separate from, and more fundamental than, whether we employ feasible versus infeasible quantities. The advantage of fixed-$n$ Studentization also goes beyond bias correction.
To begin, we condition on the covariates so that $\boldsymbol{\Gamma}_p^{-1}$ is fixed. Define $v(\cdot) = \mathbb{V}[Y \vert X = \cdot]$ and $\boldsymbol{\Sigma} = \operatorname*{diag}( v(X_i): i = 1,\ldots, n)$. Straightforward calculation gives
One can then show that $\sigma_\mathtt{us}^2 \to_P v(x) f(x)^{-1} \mathscr{V}(K,p)$, with $\mathscr{V}(K,p)$ a known, constant function of the kernel and polynomial degree. Importantly, both the nonasymptotic form and the convergence hold in the interior or on the boundary, though $\mathscr{V}(K,p)$ changes.
To first order, one could use $\sigma_\mathtt{us}^2$ or the leading asymptotic term; all that remains is to make each feasible, requiring estimators of the variance function, and for the asymptotic form, also the density. These may be difficult to estimate when $x$ is a boundary point. Concerned by this, Chen-Qin2002_SJS consider feasible and infeasible versions but conclude that “an increased coverage error near the boundary is still the case even when we know the values of $f(x)$ and $v(x)$.” Our results show that this is not true in general: using fixed-$n$ Studentization, feasible or infeasible, leads to confidence intervals with the same coverage error decay rates at interior and boundary points, thereby retaining the celebrated boundary carpentry property.
For robust bias correction, $\sigma_\mathtt{rbc}^2 = \ (n h) V[\hat{m} - \hat{B}_m \vert X_1, \ldots, X_n]$ captures the variances of $\hat{m}$ and $\hat{m}^{(p+1)}$ as well as their covariance. A fixed-$n$ calculation gives
To make the fixed-$n$ scalings feasible, $\hat{\sigma}_\mathtt{us}^2$ and $\hat{\sigma}_\mathtt{rbc}^2$ take the forms (ref) and (ref) and replace $\boldsymbol{\Sigma}$ with an appropriate estimator. First, we form $\hat{v}(X_i) = ( Y_i - \mathbf{r}_p(X_i - x)'\boldsymbol{\hat{\beta}}_p )^2$ for $\hat{\sigma}_\mathtt{us}^2$ or $\hat{v}(X_i) = ( Y_i - \mathbf{r}_q(X_i - x)'\boldsymbol{\hat{\beta}}_q )^2$ for $\hat{\sigma}_\mathtt{rbc}^2$. The latter is bias-reduced because $\mathbf{r}_p(X_i - x)'\boldsymbol{\beta}_p$ is a $p$-term Taylor expansion of $m(X_i)$ around $x$, and $\boldsymbol{\hat{\beta}}_p$ estimates $\boldsymbol{\beta}_p$ (similarly with $q$ in place of $p$), and we have $q>p$. Next, motivated by the fact that least-squares residuals are on average too small, we appeal to the HC$k$ class of estimators (see MacKinnon2013_BookChap for a review), which are defined as follows. First, $\hat{\sigma}_\mathtt{us}^2$-HC0 uses $\boldsymbol{\hat{\Sigma}}_\mathtt{us} = \operatorname*{diag}(\hat{v}(X_i): i = 1,\ldots, n)$. Then, $\hat{\sigma}_\mathtt{us}^2$-HC$k$, $k=1, 2, 3$, is obtained by dividing $\hat{v}(X_i)$ by, respectively, $(n-2\operatorname*{tr}(\mathbf{Q}_p)+\operatorname*{tr}(\mathbf{Q}_p'\mathbf{Q}_p))/n$, $(1-\mathbf{Q}_{p,ii})$, or $(1-\mathbf{Q}_{p,ii})^2$, where $\mathbf{Q}_p:=\mathbf{R}_p '\boldsymbol{\Gamma}_p^{-1} \mathbf{R}_p' \mathbf{W}_p/n$ is the projection matrix and $\mathbf{Q}_{p,ii}$ its $i$-th diagonal element. The corresponding estimators $\hat{\sigma}_\mathtt{rbc}^2$-HC$k$ are the same, but with $q$ in place of $p$. For theoretical results, we use HC0 for concreteness and simplicity, though inspection of the proof shows that simple modifications allow for the other HC$k$ estimators and rates do not change. These estimators may perform better for small sample sizes. Another option is to use a nearest-neighbor-based variance estimators with a fixed number of neighbors, following the ideas of Muller-Stadtmuller1987_AoS and Abadie-Imbens2008_AdES. Note that none of these estimators assume local or global homoskedasticity nor rely on new tuning parameters. Details and simulation results for all these estimators are given in the supplement, see \S S.II.2.3 and Table S.II.9.
Recycling notation to emphasize the parallel, we study the following three statistics: \[T_\mathtt{us} = \frac{\sqrt{nh}( \hat{m} - m)}{\hat{\sigma}_\mathtt{us}}, \qquad T_\mathtt{bc} = \frac{\sqrt{nh}( \hat{m} - \hat{B}_m - m)}{\hat{\sigma}_\mathtt{us}}, \qquad T_\mathtt{rbc} = \frac{\sqrt{nh}( \hat{m} - \hat{B}_m - m)}{\hat{\sigma}_\mathtt{rbc}},\] and their associated confidence intervals $I_\mathtt{us}$, $I_\mathtt{bc}$, and $I_\mathtt{rbc}$, exactly as in Eqn.\ (ref). Importantly, all present definitions and results are valid for an evaluation point in the interior and at the boundary of the support of $X_i$. The following standard conditions will suffice, augmented with the appropriate Cram\'er's condition given in the supplement to conserve space.
We now give our main, generic result for local polynomials, analogous to Theorem (ref). For notation, the polynomials $q_1$, $q_2$, and $q_3$ and the biases $\eta_\mathtt{us}$ and $\eta_\mathtt{bc}$, are cumbersome and exact forms are deferred to the supplement. All that matters is that the polynomials are known, odd, bounded, and bounded away from zero and that the biases have the usual convergence rates, as detailed below.
This theorem, which covers both interior and boundary points, establishes that the conclusions found in the density case carry over to odd-degree local polynomial regression. (Although we focus on $p$ odd, part (a) is valid in general and (b) and (c) are valid at the boundary for $p$ even.) In particular, this shows that robust bias correction is as good as, or better than, undersmoothing in terms of coverage error. Traditional bias correction is again inferior due to the variance and covariance terms $\rho^{p+2} (\Omega_{1,\mathtt{bc}} + \rho^{p+1} \Omega_{2,\mathtt{bc}})$. Coverage error optimal bandwidths can be derived as well, and similar conclusions are found. Best possible rates are defined for fixed $p$ here, the analogue of $\mathcal{k}$ above; see Section (ref) for further discussion on smoothness.
Before discussing bias correction, one aspect of the undersmoothing result is worth mentioning. The fact that Theorem (ref) covers both interior and boundary points, without requiring additional assumptions, is in some sense, expected: one of the strengths of local polynomial estimation is its adaptability to boundary points. In particular, from Eqn.\ (ref) and $p$ odd it follows that $\eta_\mathtt{us} \asymp \sqrt{nh} h^{p+1}$ at the interior and the boundary. Therefore, part (a) shows that the decay rate in coverage error does not change at the boundary for the standard confidence interval (but the leading constants will change). This finding contrasts with the result of Chen-Qin2002_SJS who studied the special case $p=1$ without bias correction (part (a) of Theorem (ref)), and is due entirely to our fixed-$n$ Studentization.
Turning to robust bias correction, we will, in contrast, find rate differences between the interior and the boundary, no matter the parity of $q$. As before, $\eta_\mathtt{bc}$ has two terms, representing the higher-order bias of the point estimator and the bias of the bias estimator. The former can be viewed as the bias if $m^{(p+1)}$ were zero, and since $p+1$ is even, we find that it is of order $\sqrt{nh} h^{p+3}$ in the interior but $\sqrt{nh} h^{p+2}$ at the boundary. The bias of the bias correction depends on both bandwidths $h$ and $b$, as well as $p$ and $q$, in exact analogy to the density case. For $q$ odd, it is of order $h^{p+1} b^{q-p}$ at all points, whereas for $q$ even this rate is attained at the boundary, but in the interior the order increases to $h^{p+1} b^{q+1-p}$. Collecting these facts: in the interior, $\eta_\mathtt{bc} \asymp \sqrt{nh} h^{p+3} ( 1 + \rho^{-2} b^{q-p-2} ) $ for odd $q$ or with $b^{q-p-1}$ for $q$ even; at the boundary, $\eta_\mathtt{bc} \asymp \sqrt{nh} h^{p+2} ( 1 + \rho^{-1} b^{q-p-1})$. Further details are in the supplement.
In light of these rates, the same logic of Section (ref) leads us to restrict attention to bounded, positive $\rho$ and $q=p+1$, and thus even. Calonico-Cattaneo-Titiunik2014_Ecma point out that in the special case of $q=p+1$, $K=L$, and $\rho=1$, $\hat{m} - \hat{B}_m$ is identical to a local polynomial estimator of order $q$; this is the closest analogue to $M$ being a higher-order kernel. If the point of interest is in the interior, then $q = p+2$ yields the same rates.
For notational ease, let $\tilde{\eta}_\mathtt{bc}^{\tt int}$ and $\tilde{\eta}_\mathtt{bc}^{\tt bnd}$ be the leading constants for the interior and boundary, respectively, so that e.g. $\eta_\mathtt{bc} = \sqrt{n h} h^{p+3} [\tilde{\eta}_\mathtt{bc}^{\tt int} + o(1)]$ in the interior (exact expressions are in the supplement). We then have the following, precise result; the analogue of Corollary (ref)(a).
There are differences in both the rates and constants between parts (a) and (b) of this result, though most of the changes to constants are “hidden” notationally by the definitions of $\tilde{\eta}_\mathtt{bc}^{\tt bnd}$ and the polynomials $q_{k,\mathtt{rbc}}$. Part (a) most closely resembles Corollary (ref) due to the symmetry yielding the corresponding rate improvement (recall that $\mathcal{k}$ in the density case is replaced with $p+1$ here), and hence all the corresponding conclusions hold qualitatively for local polynomials.
As we did for the density, we now derive bandwidth choices, and data-driven implementations, to optimize coverage error in applications.
To implement these results, we first set $\rho=1$ and the kernels $K$ and $L$ equal to any desired second order kernel, typical choices being triangular, Epanechnikov, and uniform. The variance estimator $\hat{\sigma}_\mathtt{rbc}^2$ is defined in Section (ref), and is fully implementable, and thus so is $I_\mathtt{rbc}$, once the bandwidth $h$ is chosen.
For selecting $h$ at an interior point, the same conclusions from density estimation apply: (i) coverage of $I_\mathtt{rbc}$ is quite robust with respect to $h$ and $\rho$, (ii) feasible choices for $h$ are easy to construct, and (iii) an MSE-optimal bandwidth only delivers the best coverage error for $p=1$ (that is, $\mathcal{k}=2$ in the density case). On the other hand, for a boundary point, an interesting consequence of Corollary (ref) is that an MSE-optimal bandwidth never delivers optimal coverage error decay rates, even for local linear regression: $h^*_\mathtt{mse} \propto n^{-1/(2p+3)} \gg h^*_\mathtt{rbc} \propto n^{-1/(p + 3)}$.
Keeping this in mind, we give a fully data-driven direct plug-in (DPI) bandwidth selector for both interior and boundary points: $\hat{h}^{\tt int}_\mathtt{dpi} = \hat{H}^{\tt int}_{\mathtt{dpi}} \; n^{-1/(p+4)}$ and $\hat{h}^{\tt bnd}_\mathtt{dpi} = \hat{H}^{\tt bnd}_{\mathtt{dpi}} \; n^{-1/(p+3)}$, where $\hat{H}^{\tt int}_{\mathtt{dpi}}$ and $\hat{H}^{\tt bnd}_{\mathtt{dpi}}$ are estimates of (the appropriate) $H^*_\mathtt{rbc}$ of Corollary (ref), obtained by estimating unknowns by pilot estimators employing a readily-available pilot bandwidth. The complete steps to form $\hat{H}^{\tt int}_{\mathtt{dpi}}$ and $\hat{H}^{\tt bnd}_{\mathtt{dpi}}$ are in the supplement, as is a second data-driven bandwidth choice, based on rescaling already-available MSE-optimal bandwidths. All our methods are available in the {\tt nprobust} package: see \url{http://sites.google.com/site/nppackages/nprobust}.
We now report a representative sample of results from a simulation study to illustrate our findings. We drew 5,000 replicated data sets, each being $n=500$ i.i.d.\ draws from the model $Y_i = m(X_i) + \varepsilon_i$, with $m(x) = \sin(3\pi x/2) (1+18x^2[\operatorname*{sgn}(x)+1])^{-1}$, $X_i \sim \mathcal{U}[0,1]$, and $\varepsilon_i \sim \mathscr{N}(0,1)$. We consider inference at the five points $x \in \{-2/3, -1/3, 0, 1/3, 2/3\}$. The function $m(x)$ and the five evaluation points are plotted in Figure (ref); this function was previously used by Berry-Carroll-Ruppert2002_JASA and Hall-Horowitz2013_AoS. The supplement gives results for other models, bandwidth selectors and their simulation distributions, alternative variance estimators, and more detailed studies of coverage and length.
We compared robust bias correction to undersmoothing, traditional bias correction, the off-the-shelf {\sf R} package {\tt locfit} locfit, and the procedure of Hall-Horowitz2013_AoS. In all cases the point estimator is based on local linear regression with the data-driven bandwidth $\hat{h}^{\tt int}_\mathtt{dpi}$, which shares the rate of $\hat{h}_\mathtt{mse}$ in this case, and $\rho=1$. The {\tt locfit} package has a bandwidth selector, but it was ill-behaved and often gave zero empirical coverage. Hall-Horowitz2013_AoS do not give an explicit optimal bandwidth, but do advocate a feasible $\hat{h}_\mathtt{mse}$, following Ruppert-Sheather-Wand1995_JASA. To implement their method, we used 500 bootstrap replications and we set $1-\xi=0.9$ over a sequence $\{x_1,...,x_N\}=\{-0.9,-0.8,\ldots,0,\ldots,0.8,0.9\}$ to obtain the final quantile $\hat{\alpha}_\xi(\alpha_0)$, and used their proposed standard errors $\hat{\sigma}_{\tt HH}^2=\kappa\hat{\sigma}^2 / \hat{f}_X$, where $\hat{\sigma}^2=\sum_{i=1}^n\hat{\varepsilon}_i^2/n$ for $\hat{\varepsilon}_i=\tilde{\varepsilon}_i-\bar{\varepsilon}$, with $\tilde{\varepsilon}_i=Y_i-\hat{m}(X_i)$ and $\bar{\varepsilon}=\sum_{i=1}^n \tilde{\varepsilon}_i / n$.
Table (ref) shows empirical coverage and average length at all five points for all five methods. Robust bias correction yields accurate coverage throughout the support; performance of the other methods varies. For $x=-2/3$, the regression function is nearly linear, leaving almost no bias, and the other methods work quite well. In contrast, at $x=-1/3$ and $x=0$, all methods except robust bias correction suffer from coverage distortions due to bias. Indeed, Hall-Horowitz2013_AoS report that “[t]he `exceptional' 100$\xi$% of points that are not covered are typically close to the locations of peaks and troughs, [which] cause difficulties because of bias.” Finally, bias is still present, though less of a problem, for $x=1/3$ and $x=2/3$, and coverage of the competing procedures improves somewhat. Motivated by the fact that the data-driven bandwidth selectors may be “too large” for proper undersmoothing, we studied the common practice of ad-hoc undersmoothing of the MSE-optimal bandwidth choice $\hat{h}_\mathtt{mse}$: the results in Table S.II.8 of the supplement show this to be no panacea.
To illustrate our findings further, Figures (ref)(a) and (ref)(b) compare coverage and length of different inference methods over a range of bandwidths. Robust bias correction delivers accurate coverage for a wide range of bandwidths, including larger choices, and thus can yield shorter intervals. For undersmoothing, coverage accuracy requires a delicate choice of bandwidth, and for correct coverage, a longer interval. Figure (ref)(c), in color online, reinforces this point by showing the “average position” of $I_\mathtt{us}(h)$ and $I_\mathtt{rbc}(h)$ for a range of bandwidths: each bar is centered at the average bias and is of average length, and then color-coded by coverage (green indicates good coverage, fading to red as coverage deteriorates). These results show that when $I_\mathtt{us}$ is short, bias is large and coverage is poor. In contrast, $I_\mathtt{rbc}$ has good coverage at larger bandwidths and thus shorter length.
This paper has made three distinct, but related points regarding nonparametric inference. First, we showed that bias correction, when coupled with a new standard error formula, performs as well or better than undersmoothing for confidence interval coverage and length. Further, such intervals are more robust to bandwidth choice in applications. Second, we showed theoretically when the popular empirical practice of using MSE-optimal bandwidths is justified, and more importantly, when it is not, and we gave concrete implementation recommendations for applications. Third, we proved that confidence intervals based on local polynomials do have automatic boundary carpentry, provided proper Studentization is used. These results are tied together through the themes of higher order expansions and the importance of finite sample variance calculations and the key, common message that inference procedures must account for additional variability introduced by bias correction.
\singlespacing
\begingroup \endgroup
\numberwithin{section}{part} \numberwithin{table}{part} \numberwithin{figure}{part}
\setenumerate[1]{label=\bf(\alph*)} \setenumerate[2]{label=\bf(\roman*)}
{2pt} {2em}
{3.5em} {3.75em} {4em}
{2em} {3em} {4em}
\partfont
\setcounter{page}{0}
\thispagestyle{empty}
\singlespacing
This supplement contains technical and notational details omitted from the main text, proofs of all results, further technical details and derivations, and additional simulations results and numerical analyses. The main results are Edgeworth expansions of the distribution functions of the $t$-statistics $T_\mathtt{us}$, $T_\mathtt{bc}$, and $T_\mathtt{rbc}$, for density estimation and local polynomial regression. Stating and proving these results is the central purpose of this supplement. The higher-order expansions of confidence interval coverage probabilities in the main paper follow immediately by evaluating the Edgeworth expansions at the interval endpoints.
Part (ref) contains all material for density estimation at interior points, while Part (ref) treats local polynomial regression at both interior and boundary points, as in the main text. Roughly, these have the same generic outline:
All our methods are implemented in {\sf R} and {\tt STATA} via the {\tt nprobust} package, available from \url{http://sites.google.com/site/nppackages/nprobust} (see also \url{http://cran.r-project.org/package=nprobust}). See Calonico-Cattaneo-Farrell2017_nprobust for a complete description.
\setcounter{tocdepth}{2} \singlespacing
\onehalfspacing
\part{Kernel Density Estimation and Inference}
Here we collect notation to be used throughout this section, even if it is restated later. Throughout this supplement, let $X_{h,i} = (x - X_i)/h$ and similarly for $X_{b,i}$. The evaluation point is implicit here. In the course of proofs we will frequently write $s=\sqrt{nh}$.
To begin, recall that the original and bias-corrected density estimators are \[\hat{f}(x) = \frac{1}{n h} \sum_{i=1}^n K\left( X_{h,i} \right)\] and
for symmetric kernel functions $K(\cdot)$ and $L(\cdot)$ that integrate to one on their compact support, $h$ and $b$ are bandwidth sequences that vanish as $n \to \infty$, and where \[\rho = h /b, \qquad \qquad \hat{B}_f = h^\mathcal{k} \hat{f}^{(\mathcal{k})}(x) \mu_{K,\mathcal{k}}, \qquad \qquad \hat{f}^{(\mathcal{k})}(x) = \frac{1}{n b^{1 + \mathcal{k}} } \sum_{i=1}^n L^{(\mathcal{k})}\left( X_{b,i} \right),\] and integrals of the kernel are denoted \[\mu_{K,k} = \frac{(-1)^k}{k!}\int u^k K(u) du, \quad\qquad\text{ and }\quad\qquad \vartheta_{K,k} = \int K(u)^k du.\]
The three statistics $T_\mathtt{us}$, $T_\mathtt{bc}$, and $T_\mathtt{rbc}$ share a common structure that is exploited to give a unified theorem statement and proof. For $v \in \{1,2\}$, define \[\hat{f}_v = \frac{1}{n h} \sum_{i=1}^n N_v \left( X_{h,i} \right), \quad\qquad \text{where} \quad\qquad N_1(u) = K(u) \text{ and } N_2(u) = M(u),\] and $M$ is given in Eqn.\ (ref). Thus, $\hat{f}_1 = \hat{f}$ and $\hat{f}_2 = \hat{f} - \hat{B}_f$. In exactly the same way, define \[\sigma^2_v := nh \mathbb{V}[\hat{f}_v] = \frac{1}{h} \left\{ \mathbb{E} \left[ N_v \left( X_{h,i} \right)^2 \right] - \mathbb{E} \left[ N_v \left( X_{h,i} \right) \right]^2 \right\}\] and the estimator \[\hat{\sigma}^2_v = \frac{1}{h} \left\{ \frac{1}{n}\sum_{i=1}^n \left[ N_v \left( X_{h,i} \right)^2 \right] - \left[ \frac{1}{n}\sum_{i=1}^n N_v \left( X_{h,i} \right) \right]^2 \right\}.\]
The statistic of interest for the generic Edgeworth expansion is, for $1 \leq w \leq v\leq 2$, \[T_{v,w} := \frac{ \sqrt{n h} (\hat{f}_v - f) }{ \hat{\sigma}_w }.\] In this notation, \[T_\mathtt{us} = T_{1,1}, \qquad T_\mathtt{bc} = T_{2,1}, \qquad \text{ and } \qquad T_\mathtt{rbc} = T_{2,2}.\]
The scaled bias is $\eta_v = \sqrt{n h} (\mathbb{E}[\hat{f}_v] - f)$. The Standard Normal distribution and density functions are $\Phi(z)$ and $\phi(z)$, respectively.
The Edgeworth expansion for the distribution of $T_{v,w}$ will consist of polynomials with coefficients that depend on moments of the kernel(s). To this end, continuing with the generic notation, for nonnegative integers $j, k, p$, define \[\gamma_{v,p} = h^{-1} \mathbb{E}\left[N_v\left(X_{h,i}\right)^p\right], \qquad \quad \qquad \Delta_{v,j} = \frac{1}{s} \sum_{i=1}^n \left\{ N_v\left(X_{h,i}\right)^j - \mathbb{E}\left[N_v\left(X_{h,i}\right)^j\right]\right\},\] and \[\nu_{v,w}(j,k,p) = \frac{1}{h}\mathbb{E}\left[ \left(N_v\left(X_{h,i}\right) - \mathbb{E}\left[N_v\left(X_{h,i}\right)\right] \right)^j \left(N_w\left(X_{h,i}\right)^p - \mathbb{E}\left[N_w\left(X_{h,i}\right)^p\right] \right)^k \right].\] We abbreviate $\nu_{v,w}(j,0,p) = \nu_v(j)$.
To expand the distribution function, additional polynomials are needed beyond those used in the main text for coverage error. These are
Next, recall from the main text the polynomials used in coverage error expansions, here with an explicit argument for a generic quantile $z$ rather than the specific $z_{\alpha/2}$:
The corresponding polynomials for expansions of the distribution function are \[q_{v,w}^{(k)}(z) = \frac{1}{2} \frac{\phi(z)}{f} q_k (z; N_w), \qquad k = 1,2,3.\]
Finally, the precise forms of $\Omega_1$ and $\Omega_2$ are: \[\Omega_1 = - 2 \frac{\mu_{K,\mathcal{k}}}{\nu_1(2)} \left\{ \int f(x - uh )K(u) L^{(\mathcal{k})}(u \rho) du - b \int f(x - uh) K(u) du \int f(x - ub) L^{(\mathcal{k})}(u) du \right\}\] and $\Omega_2 = \mu_{K,\mathcal{k}}^2 \vartheta_{K,2}^{-2} \vartheta_{L^{(\mathcal{k})},2}$. These only appear for $T_\mathtt{bc}$, and so are not indexed by $\{v,w\}$.
All these are discussed in Section (ref).
We maintain $\ell=2$ and recommend $\mathcal{k}=2$. For the kernels $K$ and $L$, we recommend either the second order minimum variance (to minimize interval length) or the MSE-optimal kernels; see Sections (ref) and (ref). In the next two subsections we discuss choice of $h$ and $\rho$.
As argued below in Section (ref), we shall maintain $\rho=1$. In the main text we give a direct plug-in (DPI) rule to implement the coverage-error optimal bandwidth. Here we we give complete details for this procedure as well as document a second practical choice, based on a rule-of-thumb (ROT) strategy. Both choices yield the optimal coverage error decay rate of $n^{-(\mathcal{k}+2)/(1+(\mathcal{k}+2))}$.
All our methods are implemented in {\sf R} and {\tt STATA} via the {\tt nprobust} package, available from \url{http://sites.google.com/site/nppackages/nprobust} (see also \url{http://cran.r-project.org/package=nprobust}). See Calonico-Cattaneo-Farrell2017_nprobust for a complete description.
Motivated by the fact that estimating $\hat{H}_{\mathtt{dpi}}$ might be difficult in practice, while data-driven MSE-optimal bandwidth selectors are readily-available, the ROT bandwidth choice is to simply rescale any feasible MSE-optimal bandwidth $\hat{h}_\mathtt{mse}$ to yield optimal coverage error decay rates (but sub-optimal constants): \[\hat{h}_\mathtt{rot} = \hat{h}_\mathtt{mse} \; n^{-(\mathcal{k}-2)/((1+2\mathcal{k})(\mathcal{k}+3))}.\] When $\mathcal{k}=2$, $\hat{h}_\mathtt{rot}=\hat{h}_\mathtt{mse}$, which is optimal (in rates) as discussed previously.
To detail the direct plug-in (DPI) rule from the main text, it is useful to first simplify the problem. Recall from the main text that the optimal choice is $h^*_\mathtt{rbc} = H^*_\mathtt{rbc} (\rho) n^{-1/(\mathcal{k} + 3)}$, where
With $\ell=2$ and $\rho = 1$, and using the definitions of $q_k(M_1)$, $k=1,2,3$, from the main text or Section (ref), this simplifies to:
where $z = z_{\alpha/2}$ the appropriate upper quantile of the Normal distribution. However, $H^*_\mathtt{rbc} (\rho)$ still depends on the unknown density through $f^{(\mathcal{k}+2)}$.
Our recommendation is a DPI rule of order one, which uses a pilot bandwidth to estimate $f^{(\mathcal{k}+2)}$ consistently. A simple and easy to implement choice is the MSE-optimal bandwidth appropriate to estimating $f^{(\mathcal{k}+2)}$, say $h^*_{\mathcal{k}+2,\mathtt{mse}}$, which is different from $h^*_\mathtt{mse}$ for the level of the function; see e.g., Wand-Jones1995_book. Let us denote a feasible MSE-optimal pilot bandwidth by $\hat{h}_{\mathcal{k}+2,\mathtt{mse}}$. Then we have:
This is now easily solved numerically (see note below). Further, if $\mathcal{k}=2$, the most common case in practice, and $K$ and $L$ are either the respective second order minimum variance or MSE-optimal kernels (Sections (ref) and (ref)), then the above may be simplified to:
Continuing with $\mathcal{k}=2$, a second option is a DPI rule of order zero, which uses a reference model to build the rule of thumb, more akin to Silverman1986_book. Using the Normal distribution, so that $f(x) = \phi(x)$ and derivatives have known form, we obtain:
where $\tilde{x} = (x - \hat{\mu}) / \hat{\sigma}_X$ is the point of interest centered and scaled.
First, we expand on the argument that $\rho$ should be bounded and positive. Intuitively, the standard errors $\hat{\sigma}_\mathtt{rbc}^2$ control variance up to order $(n h)^{-1}$, while letting $b \to 0$ faster removes more bias. If $b$ vanishes too fast, the variance is no longer controlled. Setting $\bar{\rho} \in (0,\infty)$ balances these two. Let us simplify the discussion by taking $\ell=2$, reflecting the widespread use of symmetric kernels. This does not affect the conclusions in any conceptual way, but considerably simplifies the notation. With this choice, Eqn.\ (ref) yields the tidy expression \[ \eta_\mathtt{bc} = \sqrt{nh} h^{\mathcal{k}+2} f^{(\mathcal{k}+2)} \left( \mu_{K,\mathcal{k}+2} - \rho^{-2} \mu_{K,\mathcal{k}} \mu_{L,2} \right) \; \{ 1 + o(1)\}. \] Choice of $\ell$ and $b$ (or $\rho$) cannot reduce the first term, which represents $\mathbb{E}[\hat{f}] - f - B_f$, and further, if $\bar{\rho} = \infty$, the bias rate is not improved, but the variance is inflated beyond order $(n h)^{-1}$. On the other hand, if $\bar{\rho} = 0$, then not only is a delicate choice of $b$ needed, but $\ell > 2$ is required, else the second term above dominates $\eta_\mathtt{bc}$, and the full power of the variance correction is not exploited; that is, more bias may be removed without inflating the variance rate. Hall1992_AoS_density remarked that if $\mathbb{E}[\hat{f}] - f - B_f$ is (part of) the leading bias term, then “explicit bias correction [\ldots] is even less attractive relative to undersmoothing.” We show that, on the contrary, when using our proposed Studentization, it is optimal that $\mathbb{E}[\hat{f}] - f - B_f$ is (part of) the dominant bias term. This reasoning is not an artifact of choosing $\mathcal{k}$ even and $\ell=2$, but in other cases $\rho \to 0$ can be optimal if the convergence is sufficiently slow to equalize the two bias terms.
The following result which makes the above intuition precise.
By virtue of our new studentization, the leading variance remains order $(n h)^{-1}$ and the problematic correlation terms are absent, however by forcing $\rho \to 0$, the $\rho^{-2}$ terms of $\eta_\mathtt{bc}$ are dominant (the bias of $\hat{B}_f$), and in light of our results, unnecessarily inflated. This verifies that $\bar{\rho} = 0$ or $\infty$ will be suboptimal.
We thus restrict to bounded and positive, $\rho$. Therefore, $\rho$ impacts only the shape of the “kernel” $M_\rho(u) = K(u) - \rho^{1 + \mathcal{k}} L^{(\mathcal{k})}(\rho u) \mu_{K,\mathcal{k}}$, and hence the choice of $\rho$ depends on what properties the user desires for the kernel. It happens that $\rho=1$ has good theoretical properties and performs very well numerically (see Section (ref)). As a result, from the practitioner's point of view, choice of $\rho$ (or $b$) is completely automatic.
To see the optimality of $\rho=1$, consider two cogent and well-studied possibilities: finding the kernel shape to minimize (i) interval length and (ii) MSE. The following optimal shapes are derived by Gasser-Muller-Mammitzsch1985_JRSSB and references therein. Given the above results, we set $\mathcal{k}=2$. Indeed, the optimality properties here do not extend to higher order kernels.
Minimizing interval length is (asymptotically) equivalent to finding the minimum variance fourth-order kernel, as $\sigma_\mathtt{rbc}^2 \to f \vartheta_{M,2}$. Perhaps surprisingly, choosing $K$ and $L^{(2)}$ to be the second-order minimum variance kernels for estimating $f$ and $f^{(2)}$ respectively, yields an $M_1(u)$ that is exactly the minimum variance kernel. The fourth order minimum variance kernel for estimating $f$ is $K_{\mathtt{mv}}(u) = (3/8)(-5u^2 + 3)$, which is identical to $M_1(u)$ when $K$ is the uniform kernel and $L^{(2)} = (15/4) (3u^2 -1)$, the minimum variance kernels for $f$ and $f^{(2)}$ respectively.
The result is similar for minimizing MSE: choosing $K$ and $L^{(2)}$ to be the MSE-optimal kernels for their respective point estimation problems yields an MSE-optimal $M_1(u)$. The optimal fourth order kernel is $K_{\mathtt{mse}}(u) = (15/32)(7u^4 - 10u^2 + 3)$, and the respective second-order MSE optimal kernels are $K(u) = (3/4)(1-u^2)$ and $L^{(2)}(u) = (105/16)(6u^2 - 5u^4 - 1)$. A practitioner might use the MSE-optimal kernels (along with $h^*_\mathtt{mse}$) to obtain the best possible point estimate. Our results then give an accompanying measure of uncertainty that both has correct coverage and the attractive feature of using the same effective sample.
In Section (ref) we numerically compare several kernel shapes, focusing on: (i) interval length, measured by $\vartheta_{M,2}$, (ii) bias, given by $\tilde{\mu}_{M,4}$, and (iii) the associated MSE, given by $(\vartheta_{M,2}^8 \tilde{\mu}_{M,4}^2 )^{1/9}$. These results, and the discussion above, give the foundations for our recommendation of $\rho=1$, which delivers an easy-to-implement, fully automatic choice for implementing robust bias-correction that performs well numerically, as in Section (ref).
The following assumptions are sufficient for our results. The first two are copied directly from the main text (see discussion there) and the third is the appropriate Cram\'er's condition.
It will cause no confusion (as the notations never occur in the same place), but in the course of proofs we will frequently write $s=\sqrt{nh}$.
This section accomplishes three things. First, we first carefully derive the bias of the initial estimator and the bias correction. Second, we explicate the properties of the induced kernel $M_\rho$ in terms of bias reduction and how exactly this kernel is “higher-order”. Finally, we examine two other methods of bias reduction: (i) estimating the derivatives without using derivatives of kernels Singh1977_AoS, and (ii) the generalized jackknife approach Schucany-Sommers1977_JASA. Further methods are discussed and compared by Jones-Signorini1997_JASA. The message from both alternative methods echoes our main message: it is important to account for any bias correction when doing inference, i.e., to avoid the mismatch present in $T_\mathtt{bc}$.
Recall that the biases of the two estimators are as follows:
and
The following Lemma gives a rigorous proof of these statements.
As made precise below, $M_\rho$ is a higher-order kernel. The choices of $K$, $L$, and $\rho$ determine the shape of $M_\rho$, which in turn effects the variance and bias constants. In standard kernel analyses, these constants are used to determine optimal kernel shapes for certain problems (see Gasser-Muller-Mammitzsch1985_JRSSB and references therein). For several choices of $K$, $L$, and $\rho$, Table (ref) shows numerical results for the various constants of the induced kernel $M_\rho$. The table includes (i) the variance, given by $\vartheta_{M,2}$ and relevant for interval length, (ii) a measure of bias given by $\tilde{\mu}_{M,4}$, and finally (iii) the resulting mean square error constant, $ [ \vartheta_{M,2}^8 \tilde{\mu}_{M,4}^2 ]^{1/9}$ ($\tilde{\mu}_{M,4}=(k!)(-1)^{k}\mu_{M,4}$). These specific constants are due to $M_\rho$ being a fourth order kernel, as discussed next, and would otherwise remain conceptually the same but rely on different moments. A more general, but more cumbersome procedure would be to choose $\rho$ numerically to minimize some notation of distance (e.g., $L_2$) between the resulting kernel $M_\rho$ and the optimal kernel shape already available in the literature. However, using $\rho=1$ as a simple rule-of-thumb exhibits very little lost performance, as shown in the Table and discussed in the paper.
It is worthwhile to make precise the sense in which the $n$-varying “kernel” $M_\rho(\cdot)$ of Eqn.\ (ref) is a higher-order kernel. Comparing Equations (ref) and (ref) shows exactly what is meant by this statement: the bias rate attained agrees with a standard estimate using a kernel of order $\mathcal{k}+2$ (if $\bar{\rho} > 0$), as $\ell \geq 2$. For example, if $\mathcal{k} = \ell = 2$ and $\bar{\rho} > 0$, then $M_{\bar{\rho}}(\cdot)$ behaves as a fourth-order kernel in terms of bias reduction.
However, it is not true in general that $M(\cdot)$ is a higher-order kernel in the sense that its moments below $\mathcal{k} + 2$ are zero. That is, for any $k < \mathcal{k}$, by the change of variables $w = \rho u$,
Now, $ L(u) = L(-u)$ implies that $L^{(k)}(u) = (-1)^k L^{(k)}(-u) $. Since $\mathcal{k}$ is even, $L^{(\mathcal{k})}(w)$ is symmetric, therefore if $k$ is odd $0=\int_{-\rho}^\rho w^k L^{(\mathcal{k})}(w) du$ for any $\rho$. But this fails for $k$ even, even for $\rho=1$, and hence $\int_{-1}^1 u^k M(u) du \neq 0$. For example, in the leading case of $\mathcal{k}=\ell=2$, $\int_{-1}^1 u^2 M(u) du \neq 0$ in general, and so $M(\cdot)$ is not a fourth-order kernel in the traditional sense.
Instead, the bias reduction is achieved differently. The proof of Lemma (ref) makes explicit use of the structure imposed by estimating $f^{(\mathcal{k})}$ using the derivative of the kernel $L(\cdot)$. From a technical standpoint, an integration by parts argument shows how the properties of the kernel $L(\cdot)$ (not the function $L^{(\mathcal{k})}(\cdot)$) are used to reduce bias. This argument precedes the Taylor expansion of $f$, and thus moments of $M$ are never encountered and there is no requirement that they be zero. This approach is simple, intuitive, and leads to natural restrictions on the kernel $L$, and for this reason it is commonly employed in the literature and in practice Hall1992_AoS_density.
We now examine two other methods of bias reduction: (i) estimating the derivatives without using derivatives of kernels Singh1977_AoS, and (ii) the generalized jackknife approach Schucany-Sommers1977_JASA. Further methods are discussed and compared by Jones-Signorini1997_JASA. Both methods are shown to be tightly connected to our results. Further, a more general message is that it is important to account for any bias correction when doing inference, i.e., to avoid the mismatch present in $T_\mathtt{bc}$.
The first method, which dates at least to Singh1977_AoS, is to introduce a class of kernel functions directly for derivative estimation, more closely following the standard notion of a higher-order kernel rather than using the derivative of a kernel to estimate the density derivative and proving bias reduction via integration by parts. Jones1994_CSTM expands on this method and gives further references. This class of kernels is used in the derivation of optimal kernel shapes (for derivative estimation) by Gasser-Muller-Mammitzsch1985_JRSSB. It is worthwhile to show how this class of kernel achieves bias correction and how this approach fits into our Edgeworth expansions.
Consider estimating $f^{(\mathcal{k})}$ with \[\tilde{f}^{(\mathcal{k})}(x) = \frac{1}{n b^{1 + \mathcal{k}} } \sum_{i=1}^n J \left( X_{b,i} \right),\] for some kernel function $J(\cdot)$. Note well that $J$ is generic, it need not itself be a derivative, but this is the only difference here. A direct Taylor expansion (i.e. without first integrating by parts) then gives \[\mathbb{E}[\tilde{f}^{(\mathcal{k})}] = b^{-\mathcal{k}} \sum_{k=0}^S b^k \mu_{J,k} f^{(k)} + O(b^{S + \varsigma}).\] Thus, if $J$ satisfies $\mu_{J,k} = 0$ for $k=0, 1, \ldots, \mathcal{k}-1, \mathcal{k}+1, \mathcal{k}+2, \ldots, \mathcal{k}+(\ell-1)$, $\mu_{J,\mathcal{k}} = 1$, and $\mu_{J,\mathcal{k} + \ell} \neq 0$, and $S$ is large enough then \[\mathbb{E}[\tilde{f}^{(\mathcal{k})}] = f^{(\mathcal{k})} + b^\ell f^{(\mathcal{k} + \ell)} \mu_{J,\mathcal{k} + \ell} + o(b^\ell),\] just as achieved by $\hat{f}^{(\mathcal{k})}$ and exactly matching Eqn.\ (ref). Note that $\mu_{J,0} = 0$, that is, the kernel $J$ does not integrate to one. In the language of Gasser-Muller-Mammitzsch1985_JRSSB, $J$ is a kernel of order $(\mathcal{k},\mathcal{k} + \ell)$.
Given this result, bias correction can of course be performed using $\tilde{f}^{(\mathcal{k})}(x)$ (based on $J$) rather than $\hat{f}^{(\mathcal{k})}$ (based on $L^{(\mathcal{k})}$). Much will be the same: the structure of Eqn.\ (ref) will hold with $J$ in place of $L^{(\mathcal{k})}$ and the results in Eqn.\ (ref) are achieved with modifications to the constants (e.g., in the first line, $\mu_{J,\mathcal{k} + \ell}$ appears in place of $\mu_{L,\ell}$). In either case, the same bias rates are attained. Our Edgeworth expansions will hold for this class under the obvious modifications to the notation and assumptions, and all the same conclusions are obtained.
When studying optimal kernel shapes, Gasser-Muller-Mammitzsch1985_JRSSB actually further restrict the class, by placing a limit on the number of sign changes over the support of the kernel, which ensures that the MSE and variance minimization problems have well-defined solutions. Collectively, these differences in the kernel classes explain why it is possible to demonstrate “super-optimal” MSE and variance performance for certain choices of $K$, $L^{(\mathcal{k})}$, and $\rho$, as in Table (ref).
A second alternative is the generalized jackknife method of Schucany-Sommers1977_JASA, and expanded upon by Jones-Foster1993_JNPS. To simplify the notation and ease exposition, we describe this approach for second order kernels ($\mathcal{k}=2$), but the method, and all the conclusions below, generalize fully. We thank an anonymous reviewer for encouraging us to include these details.
Begin with two estimators $\hat{f}_1$ and $\hat{f}_2$, with (possibly different) bandwidths and second-order kernels $h_j$ and $K_j$, $j=1, 2$; thus Eqn.\ (ref) gives \[\mathbb{E}[\hat{f}_j] - f(x) = h_j^2 f^{(2)} \mu_{K_j,2} + o(h_j^2), \qquad \qquad j=1, 2.\] Schucany-Sommers1977_JASA propose to estimate $f$ with $\hat{f}_{\mathtt{GJ},R} := ( \hat{f}_1 - R \hat{f}_2 ) / (1 - R)$, the bias of which is \[\mathbb{E}[ \hat{f}_{\mathtt{GJ},R} - f ] = \frac{f^{(2)}}{1-R} \left( h_1^2 \mu_{K_1,2} - R h_2^2 \mu_{K_2,2} \right) + o(h_1^2 + h_2^2).\] Hence, setting $R = (h_1^2 \mu_{K_1,2} )/( h_2^2 \mu_{K_2,2})$ renders the leading bias exactly zero. Moreover, if $S \geq 4$, $\hat{f}_{\mathtt{GJ},R}$ has bias $O(h_1^4 + h_2^4)$; behaving as a single estimator with $\mathcal{k}=4$. To put this in context of our results, observe that with this choice of $R$, if we let $\tilde{\rho} = h_1 / h_2$, then \[\hat{f}_{\mathtt{GJ},R} = \frac{1}{n h_1} \sum_{i=1}^n \tilde{M}\left(\frac{X_i - x}{h_1} \right) , \quad M(u) = K_1(u) - \tilde{\rho}^{1+2} \left\{ \frac{K_2(\tilde{\rho}u) - \tilde{\rho}^{-1} K_1(u) }{\mu_{K_2,2}(1-R)} \right\} \mu_{K_1, 2},\] exactly matching Eqn.\ (ref). Or equivalently, $\hat{f}_{\mathtt{GJ},R} = \hat{f}_1 - h_1^2 \tilde{f}^{(2)} \mu_{K_1,2}$, for the derivative estimator \[\tilde{f}^{(2)} = \frac{1}{nh_2^{1+2}}\sum_{i=1}^n \tilde{L}\left( \frac{X_i - x}{h_2} \right), \quad \tilde{L}(u) = \frac{ K_2(u) - \tilde{\rho}^{-1} K_1( \tilde{\rho}^{-1} u) }{\mu_{K_2,2}(1-R)}. \] Therefore, we can view $\hat{f}_{\mathtt{GJ},R}$ as a change in the kernel $M(\cdot)$ or an explicit bias estimation described directly above with a specific choice of $J(\cdot)$ (depending on $\tilde{\rho}$ in either case). Again, Eqn.\ (ref) holds exactly. Thus, our results cover the generalized jackknife method as well, and the same lessons apply.
Finally, we note that these bias correction methods can be applied to nonparametric regression as well, and local polynomial regression in particular, and that the same conclusions are found. We will not repeat this discussion however.
Here we briefly state the first-order properties of $T_\mathtt{us}$, $T_\mathtt{bc}$, and $T_\mathtt{rbc}$, using the common notation $T_{v,w}$ defined in Section (ref). Recall that $\eta_v = \sqrt{n h} (\mathbb{E}[\hat{f}_v] - f)$ is the scaled bias in either case. With this notation, we have the following result.
The conditions on $h$ and $b$ behind the generic assumption that the scaled bias vanishes can be read off of (ref) and (ref): $T_\mathtt{us}$ requires $\sqrt{nh} h^\mathcal{k} \to 0$ whereas $T_\mathtt{bc}$ and $T_\mathtt{rbc}$ require only $\sqrt{n h} h^\mathcal{k} (h^2 \vee b^\ell) \to 0$, and thus accommodate $\sqrt{n h} h^\mathcal{k} \not\to 0$ or $b \not\to 0$ (but not both). However, bias correction requires a choice of $\rho=h/b$. One easily finds that $\mathbb{V}[\sqrt{nh} \hat{B}_f] = O(\rho^{1 + 2\mathcal{k}})$, whence $\rho \to 0$ is required for $T_\mathtt{bc}$. But $T_\mathtt{rbc}$ does not suffer from this requirement because of our proposed, new Studentization. From a first-order point of view, traditional bias correction allows for a larger class of sequences $h$, but requires a delicate choice of $\rho$ (or $b$), and Hall1992_AoS_density shows that this constraint prevents $T_\mathtt{bc}$ from improving inference. Our novel standard errors remove these constraints, allowing for improvements in bias to carry over to improvements in inference. The fact that a wider range of bandwidths is allowed hints at the robustness to tuning parameter choice discussed above and formalized by our Edgeworth expansions.
Recall the generic notation: \[T_{v,w} := \frac{\sqrt{nh} (\hat{f}_v - f) }{ \hat{\sigma}_w },\] for $1 \leq w \leq v\leq 2$. The Edgeworth expansion for the distribution of $T_{v,w}$ will consist of polynomials with coefficients that depend on moments of the kernel(s). Additional polynomials are needed beyond those used in the main text for coverage error. These are:
The polynomials $p_{v,w}^{(k)}$ are even, and hence cancel out of coverage probability expansions, but are used in the expansion of the distribution function itself (or equivalently, the coverage of a one-sided confidence interval).
Next, recall from the main text the polynomials used in coverage error expansions:
The corresponding polynomials for expansions of the distribution function are \[q_{v,w}^{(k)}(z) = \frac{1}{2} \frac{\phi(z)}{f} q_k (z; N_w), \qquad k = 1,2,3.\] As before, the $q_{v,w}^{(k)}$ are odd and hence do not cancel when computing coverage: the $q_k (z; N_w)$ in the main text are doubled for just this reason.
Note that, despite the notation, $q_{v,w}^{(k)}(z)$ depends only on the “denominator” kernel $N_w$. The notation comes from the fact that when first computed, the terms which enter into the $q_{v,w}^{(k)}(z)$ depend on both kernels, but the simplifications in Eqn.\ (ref) reduce the dependence to $N_w$. This is because for undersmoothing and robust bias correction, $v=w$, and for traditional bias correction $N_2 = M = K + o(1) = N_1 + o(1)$, as $\rho \to 0$ is assumed. Thus, when computing $\vartheta_{M,q}$ the terms with the lowest powers of $\rho$ will be retained. These can be found by expanding \[\vartheta_{M,q} = \int \left(K(u) - \rho^{1 + \mathcal{k}}\mu_{K,\mathcal{k}} L^{(\mathcal{k})}(u)\right)^q du = \sum_{j=0}^q {q \choose j}\left(-\mu_{K,\mathcal{k}}\rho^{1+\mathcal{k}}\right)^{q-j}\int K(u)^j L^{(\mathcal{k})}(\rho u)^{q-j}du,\] and hence we can write $\vartheta_{M,q} = \vartheta_{K,q} - \rho^{1+\mathcal{k}} q \mu_{K,\mathcal{k}} L^{(\mathcal{k})}(0) \vartheta_{K,q-1} + O(h + \rho^{2 + \mathcal{k}})$. We can thus write $q_j(z ; M) = q_j(z ; K) + o(1)$ in this case. If the expansions were carried out beyond terms of order $(nh)^{-1} + (nh)^{-1/2}\eta_v + \eta_v^2 + \mathbbm{1}\{v \!\neq\! w\} \rho^{1+2\mathcal{k}}$ this would not be the case.
Finally, for traditional bias correction, there are additional terms in the expansion (see discussion in the main text) representing the covariance of $\hat{f}$ and $\hat{B}_f$ (denoted by $\Omega_1$) and the variance of $\hat{B}_f$ ($\Omega_2$). We now state their precise forms. These arise from the mismatch between the variance of the numerator of $T_\mathtt{bc}$ and the standardization used, $\sigma_\mathtt{us}^2$, that is $\sigma_\mathtt{rbc}^2/\sigma_\mathtt{us}^2$ is given by \[ \frac{ nh \mathbb{V}[\hat{f} - \hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } = \frac{ nh \mathbb{V}[\hat{f}] - 2 nh \mathbb{C}[\hat{f},\hat{B}_f] + nh \mathbb{V}[\hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } = 1 - 2 \frac{ nh \mathbb{C}[\hat{f},\hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } + \frac{ nh \mathbb{V}[\hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } .\] This makes clear that $\Omega_1$ and $\Omega_2$ are the constant portions of the last two terms. We have \[ - 2 \frac{ nh \mathbb{C}[\hat{f},\hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } = \rho^{1+\mathcal{k}} \Omega_1, \] where \[\Omega_1 = - 2 \frac{\mu_{K,\mathcal{k}}}{\nu_1(2)} \left\{ \int f(x - uh )K(u) L^{(\mathcal{k})}(u \rho) du - b \int f(x - uh) K(u) du \int f(x - ub) L^{(\mathcal{k})}(u) du \right\}.\] Note $\nu_1(2) = \sigma_\mathtt{us}^2$. Turning to $\Omega_2$, using the calculations in Section (ref) (recall $\tilde{\mathcal{k}} = \mathcal{k} \vee S$), we find that \[\frac{ nh \mathbb{V}[\hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } = \rho^{1 + 2\mathcal{k}} \Omega_2 \quad \text{ where }\quad \Omega_2 = \frac{ \mu_{K,\mathcal{k}}^2}{ \nu_1(2) } \left\{ \int f(x - ub) L^{(\mathcal{k})}(u)^2 du - b^{1 + 2\tilde{\mathcal{k}}} \left( \int L^{(\mathcal{k}-\tilde{\mathcal{k}})}(u) f^{(\tilde{\mathcal{k}})}(x - u b) du \right)^2 \right\}.\] Fully simplifying would yield \[\Omega_2 = \mu_{K,\mathcal{k}}^2 \vartheta_{K,2}^{-2} \vartheta_{L^{(\mathcal{k})},2},\] which can be used in Theorem (ref).
As a last piece of notation, define the scaled bias as $\eta_v = \sqrt{n h}(\mathbb{E}[\hat{f}_v] - f)$.
We can now state our generic Edgeworth expansion, from whence the coverage probability expansion results follow immediately.
To use this result to find the expansion of the error in coverage probability of the Normal-based confidence interval, the function $F_{v,w}(z)$ is simply evaluated at the two endpoints of the interval. (Note: if the confidence interval were instead constructed with the bootstrap, a few additional steps are needed, but these do not alter any conclusions or results outside of constant terms.)
In general, we have assumed that the level of smoothness was large enough to be inconsequential in the analysis, and in particular this allowed for characterization of optimal bandwidth choices. In this section, in contrast, we take the level of smoothness to be binding, so that we can fully utilize the $S$ derivatives and the H\"older condition to obtain the best possible rates of decay in coverage error for both undersmoothing and robust bias correction, but at the price of implementability: the leading bias constants can not be characterized, and hence feasible “optimal” bandwidths are not available.
For undersmoothing, the lowest bias is attained by setting $\mathcal{k}>S$ (see Eqn.\ (ref)), in which case the bias is only known to satisfy $\mathbb{E}[\hat{f}] - f = O(h^{S+\varsigma})$ (i.e., $B_f$ is identically zero) and bandwidth selection is not feasible. Note that this approach allows for $\sqrt{n h} h^S \not\to 0$, as $\eta_\mathtt{us} = O(\sqrt{n h} h^{S+\varsigma})$.
Robust bias correction has several interesting features here. If $\mathcal{k} \leq S-2$ (the top two cases in Eqn.\ (ref)), then the bias from approximating $\mathbb{E}[\hat{f}] - f$ by $B_f$, that is not targeted by bias correction, dominates $\eta_\mathtt{bc}$ and prevents robust bias correction from performing as well as the best possible infeasible (i.e., oracle) undersmoothing approach. That is, even bias correction requires a sufficiently large choice of $\mathcal{k}$ in order to ensure the fastest possible rate of decay in coverage error: if $\mathcal{k} \geq S-1$, robust bias correction can attain error decay rate as the best undersmoothing approach, and allow $\sqrt{n h} h^S \not\to 0$.
Within $\mathcal{k} \geq S-1$, two cases emerge. On the one hand, if $\mathcal{k} =S-1$ or $S$, then $B_f$ is nonzero and $f^{(\mathcal{k})}$ must be consistently estimated to attain the best rate. Indeed, more is required. From Eqn.\ (ref), we will need a bounded, positive $\rho$ to equalize the bias terms. This (again) highlights the advantage of robust bias correction, as the classical procedure would enforce $\rho \to 0$, and thus underperform. On the other hand, $\rho \to 0$ will be required if $\mathcal{k}>S$ because (from the final case of (ref)) we require $\rho^{\mathcal{k}-S} = O(h^\varsigma)$ to attain the same rate as undersmoothing. Note that we can accommodate $b \not\to 0$ (but bounded). Interestingly, $B_f$ is identically zero and $\hat{B}_f$ merely adds noise to the problem, but this noise is fully accounted for by the robust standard errors, and hence does not affect the rates of coverage error (though the constants of course change). The $\hat{f}^{(\mathcal{k})}$ in $\hat{B}_f$ is inconsistent ($f^{(\mathcal{k})}$ does not exist), but the nonvanishing bias of $\hat{f}^{(\mathcal{k})}$ is dominated by $h^\mathcal{k}$.
This discussion is summarized by the following result:
We now briefly present state analogues of our results, both for distributional convergence and Edgeworth expansions, that cover multivariate data and derivative estimation. The conceptual discussion and implications are similar to those in the main text, once adjusted notationally to the present setting, and are hence omitted.
For a nonnegative integral $d$-vector $q$ we adopt the notation that: (i) $[q] = q_1 + \cdots + q_d$, (ii) $g^{(q)}(x) = \partial^{[q]} g(x)/(\partial^{q_1} x_1 \cdots \partial^{q_d} x_d)$, (iii) $k! = q_1!\cdots q_d!$, and (iv) $\sum_{[q] = Q}$ for some integer $Q \geq 0$ denotes the sum over all indexes in the set $\{q : [q] = Q\}$.
The parameter of interest is $f^{(q)}(x)$, for $x \in \mathbb{R}^d$ and $[q] \leq S$. The estimator is \[\hat{f}^{(q)}(x) = \frac{1}{n h^{d + [q]} } \sum_{i=1}^n K^{(q)}\left( X_{h,i} \right).\] Note that here, and below for bias correction, we use a constant, diagonal bandwidth matrix, e.g. $h \times I_d$. This is for simplicity and comparability, and could be relaxed at notational expense.
The bias, for a given kernel of order $\mathcal{k} \leq S - [q]$ (we restrict attention to the case where $S$ is large enough), is \[h^\mathcal{k} \sum_{k: [k + q] = \mathcal{k}} \mu_{K,k} f^{(q+k)}(x) + o(h^\mathcal{k} ),\] exactly mirroring Eqn.\ (ref), where now $\mu_{K,k}$ represents a $d$-dimensional integral. Bias estimation is straightforward, relying on estimates $\hat{f}^{(q+k)}(x)$, for all $[k] = \mathcal{k} - [q]$. The form of $\hat{f}^{(q)}_2(x) = \hat{f}^{(q)}(x) - \hat{B}_{f^{(q)}}(x)$ is now given by \[ \hat{f}^{(q)}_2(x) = \frac{1}{n h^{d + [q]}} \sum_{i=1}^n M_{(q)}\left( X_{h,i}\right) \quad \text{where} \quad M_{(q)}(u)= K^{(q)}(u) - \left(\rho \right)^{d + [q] + \mathcal{k}} \sum_{[k] = \mathcal{k} } \mu_{K,k} L^{(q+k)}(u), \] exactly analogous to Eqn.\ (ref).
With these changes in notation out of the way, we can (re-)define the generic framework for both estimators exactly as above. Dropping the point of evaluation $x$, for $v \in \{1,2\}$, define the estimator as \[\hat{f}^{(q)}_v = \frac{1}{n h^{d + [q]} } \sum_{i=1}^n N_v \left( X_{h,i} \right), \quad\qquad \text{where} \quad\qquad N_1(u) = K^{(q)}(u) \text{ and } N_2(u) = M_{(q)}(u);\] the variance \[\sigma^2_v := n h^{d + [q]} \mathbb{V}[\hat{f}^{(q)}_v] = \frac{1}{h^d} \left\{ \mathbb{E} \left[ N_v \left( X_{h,i} \right)^2 \right] - \mathbb{E} \left[ N_v \left( X_{h,i} \right) \right]^2 \right\}\] and its estimator as \[\hat{\sigma}^2_v = \frac{1}{h^d} \left\{ \frac{1}{n}\sum_{i=1}^n \left[ N_v \left( X_{h,i} \right)^2 \right] - \left[ \frac{1}{n}\sum_{i=1}^n N_v \left( X_{h,i} \right) \right]^2 \right\};\] and the $t$-statistics, for $1 \leq w \leq v\leq 2$, as, \[T_{v,w} := \frac{ \sqrt{n h^{d + 2[q]} } \left(\hat{f}^{(q)}_v - f^{(q)} \right) }{ \hat{\sigma}_w }.\] As before, $T_\mathtt{us} = T_{1,1}$, $T_\mathtt{bc} = T_{2,1}$, and $T_\mathtt{rbc} = T_{2,2}$.
The scaled bias $\eta_v$ has the same general definition as well: the bias of the numerator of the $T_{v,w}$. In this case, given by \[\eta_v = \sqrt{n h^{d + 2[q]} } \left( \mathbb{E} \left[\hat{f}^{(q)}_v\right] - f^{(q)}(x) \right). \] The asymptotic order of $\eta_v$ for different settings can be obtained straightforwardly via the obvious multivariate extensions of Equation (ref) and the corresponding conclusion of Lemma (ref).
First-order convergence is now given by the following result. the proof of which is standard.
For the Edgeworth expansion, redefine \[\nu_{v,w}(j,k,p) = \frac{1}{h^{d + [q]\mathbbm{1}\{j + pk=1\}} } \mathbb{E}\left[ \left(N_v(u_i) - \mathbb{E}[N_v(u_i)] \right)^j \left(N_w(u_i)^p - \mathbb{E}[N_v(u_i)^p] \right)^k \right],\] where $u_i = (x - X_i ) / h$. The polynomials $p_{v,w}^{(k)}(z)$ and $q_{v,w}^{(k)}(z)$ are as given above, but using multivariate moments. The analogue of Theorem (ref) is given by the following result, which can be proven following the same steps as in Section (ref).
The same conclusions reached in the main text continue to hold for multivariate and/or derivative estimation, both in terms of comparing undersmoothing, bias correction, and robust bias correction, as well as for inference-optimal bandwidth choices. In particular, it is straightforward that the MSE optimal bandwidth in general has the rate $n^{- 1 / (d + 2 \mathcal{k} + 2[q])}$, whereas the coverage error optimal choice is of order $n^{- 1 / (d + \mathcal{k} + [q])}$. Note that these two fit the same patter as in the univariate, level case, with $\mathcal{k} + [q]$ in place of $\mathcal{k}$ and $d$ in place of one. One intuitive reason for the similarity is that the number of derivatives in question does not impact that variance or higher order moment terms of the expansion, once the scaling is accounted for. That is, for all averages beyond the first, for example of the kernel squared, $\sqrt{nh^d}$ can be thought of as the effective sample size, since that is the multiplier which stabilizes averages.
Throughout $C$ shall be a generic constant that may take different values in different uses. If more than one constant is needed, $C_1$, $C_2$, \ldots, will be used. It will cause no confusion (as the notations never occur in the same place), but in the course of proofs we will frequently write $s=\sqrt{nh}$, which overlaps with the order of the kernel $L$.
The first step is to write $T_{v,w}$ as a smooth function of sums of i.i.d.\! random variables plus a remainder term that is shown to be of higher order. In addition to the notation above, define \[\gamma_{v,p} = h^{-1} \mathbb{E}\left[N_v\left(X_{h,i}\right)^p\right] \qquad \text{ and } \qquad \Delta_{v,j} = \frac{1}{s} \sum_{i=1}^n \left\{ N_v\left(X_{h,i}\right)^j - \mathbb{E}\left[N_v\left(X_{h,i}\right)^j\right]\right\}.\] With this notation $\hat{f}_v - \mathbb{E}[\hat{f}_v] = s^{-1} \Delta_{v,1}$, $\sigma^2_w = \mathbb{E}[\Delta_{w,1}^2] = \gamma_{w,2} - h \gamma_{w,1}^2$ and
By a change of variables \[\gamma_{v,p} = h^{-1} \int N_v\left(X_{h,i}\right)^p f(X_i) d X_i = \int N_v(u)^{p} f(x - u h) d u = O(1).\] Further, by construction $\mathbb{E}[\Delta_{w,j}] = 0$ and
Returning to Eqn.\ (ref) and applying Markov's inequality, we find that $h s^{-2} \Delta_{w,1}^2 = n^{-1} \Delta_{w,1}^2 = O_p(n^{-1})$ and $\hat{\sigma}^2_w - \sigma^2_w = s^{-1} O_p(1) - h O(1) s^{-1} O_p(1) - h s^{-2} O_p(1) = O_p(s^{-1})$, whence $\left| \hat{\sigma}^2_w - \sigma^2_w \right|^2 = O_p(s^{-2})$. Using these results preceded by a Taylor expansion, we have
Combining this result with the fact that \[T_{v,w} = \frac{\Delta_{v,1} + \eta_v}{\hat{\sigma}_w} = \frac{\Delta_{v,1}}{\hat{\sigma}_w} + \frac{\eta_v}{\sigma_w} \left( \frac{\hat{\sigma}^2_w}{\sigma^2_w} \right)^{-1/2},\] we have
where \[\tilde{T}_{v,w} = \frac{\Delta_{v,1}}{\hat{\sigma}_w} - \frac{\eta_v}{2 \sigma^3_w} \left( s^{-1} \Delta_{w,2} - h 2 \gamma_{w,1} s^{-1} \Delta_{w,1} \right) \] and is a smooth function of sums of i.i.d.\! random variables and the remainder term is \[R_{v,w} = \frac{\eta_v}{\sigma_w} \left( h s^{-2} \frac{\Delta_{w,1}^2}{2 \sigma^2_w} + \frac{3}{8} \frac{(\hat{\sigma}^2_w - \sigma^2_w)^2}{\sigma^4_w} + o_p((\hat{\sigma}^2_w - \sigma^2_w)^2) \right). \]
Next we apply the delta method, see Hall1992_book or Andrews2002_Ecma. It will be true that
if it can be shown that $s^2 \mathbb{P}[|R_{v,w}| > \varepsilon^2 s^{-2} \log(s)^{-1} ] = o(1)$.\footnote{Here, $s^{-2} \log(s)^{-1}$ may be replaced with any sequence that is $o(s^{-2} + \eta_v^2 + s^{-1} \eta_v)$.} This can be demonstrated by applying Bernstein's inequality to each piece of $R_{v,w}$, as the kernels $K$ and $L$, and their derivatives, are bounded.
To apply this inequality to the first term of $R_{v,w}$, note that $| N_w ((x - X_i)/h) | \leq C_1$ and that $\mathbb{V}[N_w ((x - X_i)/h)] \leq C_2 h$, for different constants, and so for $\varepsilon > 0$ we have
which tends to zero because $\eta_v \to 0$ as $n \to \infty$ is assumed. To see why, note first that the second term of the denominator automatically vanishes, as $\eta_v \to 0$ and $\log(s)^3/n \to 0$. Second, suppose $\eta_v^2 \asymp n h^\omega$ (for example, if $\eta_\mathtt{us} \asymp s h^\mathcal{k}$, then $\omega = 1 + 2\mathcal{k}$) and the first term diverges, it must be that $h$ is at least as large (in order) as \[\left( \frac{1}{n \log(s)^4} \right)^{1/(2 + \omega)},\] which makes the requirement that $\eta_v \to 0$ equivalent to \[\eta_v^2 \asymp n h^\omega = n^{1 - \omega/(2 + \omega)} \log(s)^{ - 4 \omega/(2 + \omega)} \to 0,\] which is impossible. The remaining terms of $R_{v,w}$, characterized using Eqn.\ (ref), are handled in exactly the same way. This establishes Eqn.\ (ref).
Next, the proofs of Hall1992_book show that $\tilde{T}_{v,w}$ has an Edgeworth expansion valid through $o(s^{-2} + s^{-1}\eta_v + \eta_v^2)$. Thus, for a smooth function $G(z)$ we can write $\mathbb{P}[\tilde{T}_{v,w} < z ] = G(z) + o(s^{-2} + s^{-1}\eta_v + \eta_v^2)$. Therefore
The final result now follows by combining Equations (ref), (ref), and (ref) with the terms of the expansion computed below.\qed
Identifying the terms of the expansion is a matter of straightforward, if tedious, calculation. The first four cumulants of $T_{v,w}$ must be calculated, which are functions of the first four moments. In what follows, we give a short summary. Note well that we always discard higher-order terms for brevity, and to save notation we will write $\stackrel{o}{=}$ to stand in for “equal up to $o((nh)^{-1} + (nh)^{-1/2}\eta_v + \eta_v^2 + \mathbbm{1}\{v \!\neq\! w\} \rho^{1+2\mathcal{k}} )$”.
Referring to the Taylor expansion above, for the purpose of computing moments and cumulants, we can use \[T_{v,w} \approx \left( \frac{\Delta_{v,1}}{\sigma_w} + \frac{\eta_v}{\sigma_w} \right) \left(1 - \frac{ s^{-1} \Delta_{w,2}}{2 \sigma_w} + \frac{h \gamma_{w,1} s^{-1} \Delta_{w,1}}{\sigma_w} + \frac{3}{8} \frac{s^{-2} \Delta_{w,2}^2}{\sigma^2_w} \right).\] Moments of the two sides agree up to the requisite order. Straightforward moment calculations then give
and,
The expansion now follows, formally, from the following steps. First, combining the above moments into cumulants. Second, these cumulants may be simplified using that \[\frac{\sigma_v^2}{\sigma_w^2} = 1 + \mathbbm{1}(w\!\neq\! v) \left( \rho^{1+ \mathcal{k}} \Omega_1 + \rho^{1 + 2\mathcal{k}} \Omega_2\right) \] and in all cases present
The second relation is readily proven for $v=w$, as $\nu_{v,v}(i,j,p) = \mathbb{E}[N_v(X_{h,i})^{i + jp}] + O(h)$, where the remainder represents products of expectations. In the case for $v \neq w$, we find $\nu_{2,1}(i,j,p) = f \vartheta_{N_1,i + jp} + O( \rho^{1+\mathcal{k}} + h)$, and in this case $\rho \to 0$ is assumed. For any term of a cumulant with a rate of $(nh)^{-1}$, $(nh)^{-1/2}\eta_v$, $\eta_v^2$, or $\rho^{1+2\mathcal{k}}$ (i.e., the extent of the expansion), these simplifications may be inserted as the remainder will be negligible. Note that this is exactly why the polynomials $p_{v,w}^{(k)}$ do not simplify, while the $q_{v,w}^{(k)}$ do. Third, with the cumulants in hand, the terms of the expansion are determined as described by e.g., Hall1992_book.
Finally, for traditional bias correction, there are additional terms in the expansion (see discussion in the main text) representing the covariance of $\hat{f}$ and $\hat{B}_f$ (denoted by $\Omega_1$) and the variance of $\hat{B}_f$ ($\Omega_2$). We now state their precise forms. These arise from the mismatch between the variance of the numerator of $T_\mathtt{bc}$ and the standardization used, $\sigma_\mathtt{us}^2$, that is $\sigma_\mathtt{rbc}^2/\sigma_\mathtt{us}^2$ is given by \[ \frac{ nh \mathbb{V}[\hat{f} - \hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } = \frac{ nh \mathbb{V}[\hat{f}] - 2 nh \mathbb{C}[\hat{f},\hat{B}_f] + nh \mathbb{V}[\hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } = 1 - 2 \frac{ nh \mathbb{C}[\hat{f},\hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } + \frac{ nh \mathbb{V}[\hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } .\] This makes clear that $\Omega_1$ and $\Omega_2$ are the constant portions of the last two terms. First, for $\Omega_1$,
Therefore \[ - 2 \frac{ nh \mathbb{C}[\hat{f},\hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } = \rho^{1+\mathcal{k}} \Omega_1, \] where \[\Omega_1 = - 2 \frac{\mu_{K,\mathcal{k}}}{\nu_1(2)} \left\{ \int f(x - uh )K(u) L^{(\mathcal{k})}(u \rho) du - b \int f(x - uh) K(u) du \int f(x - ub) L^{(\mathcal{k})}(u) du \right\}.\] Note $\nu_1(2) = \sigma_\mathtt{us}^2$. If we did not include $\Omega_2$ in the Edgeworth expansion, i.e. we stopped at order $\rho^{1+\mathcal{k}}$, then we could capture only the leading terms of $\Omega_1$, as follows, using that kernel integrates to 1 and $\rho \to 0$,
Note that this matches the term Hall1992_AoS_density calls $w_2$. We do not do this, for completeness. There are no other terms of up to order $\rho^{1+2\mathcal{k}}$, so capturing the full contribution of $\sigma_2^2 / \sigma_1^2 - 1 = \sigma_\mathtt{rbc}^2 / \sigma_\mathtt{us}^2 - 1$ is natural and informative.
Turning to $\Omega_2$, using the calculations in Section (ref) (recall $\tilde{\mathcal{k}} = \mathcal{k} \vee S$), we find that
and hence \[\frac{ nh \mathbb{V}[\hat{B}_f] }{ nh \mathbb{V}[\hat{f}] } = \rho^{1 + 2\mathcal{k}} \Omega_2 \quad \text{ where }\quad \Omega_2 = \frac{ \mu_{K,\mathcal{k}}^2}{ \nu_1(2) } \left\{ \int f(x - ub) L^{(\mathcal{k})}(u)^2 du - b^{1 + 2\tilde{\mathcal{k}}} \left( \int L^{(\mathcal{k}-\tilde{\mathcal{k}})}(u) f^{(\tilde{\mathcal{k}})}(x - u b) du \right)^2 \right\}.\] The final piece will be $b^{1 + 2S} f^{(\mathcal{k})}(x)^2 [ 1+o(1)]$ if $\mathcal{k} \leq S$. Substituting this is permitted because $\rho^{1+2\mathcal{k}}$ is the limit of the expansion, though it is not necessary to do, because this term is always higher order. Fully simplifying would yield \[\Omega_2 = \mu_{K,\mathcal{k}}^2 \vartheta_{K,2}^{-2} \vartheta_{L^{(\mathcal{k})},2},\] which can be used in Theorem (ref).
To illustrate the gains from robust bias correction we conduct a Monte Carlo study to compare undersmoothing, traditional bias correction, and robust bias correction in terms coverage accuracy and interval length using several data-driven procedures to select the bandwidth. We generate $n=500$ observations from a density $f$ given by:
We evaluate the density at $x=\{-2,-1,0,1,2\}$. These models were previously analyzed in Marron-Wand1992_AoS and they are plotted in Figure (ref). In this simulation study we compare the performance of the confidence intervals defined by $T_\mathtt{us}$, $T_\mathtt{bc}$, and $T_\mathtt{rbc}$. For $T_\mathtt{us}$, we take $K$ to be the Epanechnikov kernel, while bias correction uses the Epanechnikov and MSE-optimal kernels for $K$ and $L^{(2)}$, respectively. The bandwidth $h$ is chosen in three different ways:
Empirical coverage and length are reported in Tables (ref)--(ref) (Panel A) using our two proposed data-driven bandwidth selectors, as well as the infeasible $h_\mathtt{mse}$. The most obvious finding is that robust bias correction has accurate coverage for all bandwidth choices in all models. The intervals are generally longer than for undersmoothing, but neither undersmoothing nor traditional bias correction yield correct coverage outside of a few special cases (e.g., undersmoothing at the infeasible MSE-optimal bandwidth in Model 4). The DPI bandwidth selector generally results in slightly smaller bandwidths (on average). Summary statistics for the two fully data-driven bandwidths are shown in Panel B. The fact that the DPI bandwidth is slightly smaller is born out. It is also, in general, more variable.
To illustrate the robustness to tuning parameter selection, Figures (ref)--(ref) show coverage and length for all four models. The dotted vertical line shows the population MSE-optimal bandwidth for reference. These figures demonstrate the delicate balance required for undersmoothing to provide correct coverage, whereas for a wide range of bandwidths robust bias correction provides correct coverage. Further, interval length is not unduly inflated for bandwidths that provide correct coverage. Recall that robust bias correction can accommodate, and will optimally employ, a larger bandwidth, yielding higher precision. Further emphasizing the point of robustness, we depart from $\rho = 1$ in Figures (ref) and (ref) to show coverage and length over a grid of $h$ and $\rho$.
The simulation results for local polynomial regression reported in Section (ref) below bear out these same conclusions and study these issues in more detail, in particular interval length.
All our methods are implemented in {\sf R} and {\tt STATA} via the {\tt nprobust} package, available from \url{http://sites.google.com/site/nppackages/nprobust} (see also \url{http://cran.r-project.org/package=nprobust}). See Calonico-Cattaneo-Farrell2017_nprobust for a complete description.
\newcounter{model} \forloop{model}{1}{\value{model} < 5}{
}
\part{Local Polynomial Estimation and Inference}
Local polynomial regression is notationally demanding, and the Edgeworth expansions will be substantially more so. For ease of reference, we collect all notation here regardless of where it is introduced and used. Much of the notation is fully restated later, when needed. As such, this subsection is designed more for reference, and is not easily readable.
Throughout, a subscript $p$ will generally refer to a quantity used to estimate $m(x)=\mathbb{E}[Y_i|X_i=x]$, while a subscript $q$ will refer to the bias correction portion (the vectors $e_0$ and $e_{p+1}$ below are notable exceptions to this rule). Recall that $p \geq 1$ is odd and $q>p$ may be even or odd.
Throughout this section let $X_{h,i} = (X_i - x)/h$ and similarly for $X_{b,i}$. The evaluation point is implicit here.
To save notation, products of functions will be written together, with only one argument. For example \[ (K r_p r_p')(X_{h,i}) := K(X_{h,i}) r_p(X_{h,i}) r_p(X_{h,i})' = K\left(\frac{X_i - x}{h}\right) r_p \left(\frac{X_i - x}{h}\right) r_p \left(\frac{X_i - x}{h}\right)' , \] and similarly for $(K r_p)(X_{h,i})$, $(L r_q) (X_{b,i})$, etc.
All expectations are fixed-$n$ calculations. To give concrete examples of this notation ($\boldsymbol{\Lambda}_{p,k}$, $R_p$, and $W_p$ are redefined below): \[\boldsymbol{\Lambda}_{p,k} = R_p' W_p [ ((X_1 - x)/h)^{p+k}, \cdots, ((X_n - x)/h)^{p+k}]'/n = \frac{1}{nh} \sum_{i=1}^n (K r_p)(X_{h,i}) X_{h,i}^{p+k}\] and \[\tilde{\boldsymbol{\Lambda}}_{p,k} = \mathbb{E}[\boldsymbol{\Lambda}_{p,k}] = h^{-1} \mathbb{E}[(K r_p)(X_{h,j}) X_{h,i}^{p+k}] = h^{-1} \int_{\operatorname*{supp}\{X\}} K\left(\frac{X_i - x}{h}\right) r_p \left(\frac{X_i - x}{h}\right) \left(\frac{X_i - x}{h}\right)^{p+k} f(X_i) dX_i.\] Here the range of integration is explicit, but in general it will not be. This is important for boundary issues, where the notation is generally unchanged, and it is to be understood that moments and moments of the kernel be replaced by the appropriate truncated version. Continuing this example, if $\operatorname*{supp}\{X\} = [0,\infty)$ and $x = 0$, then by a change of variables \[\tilde{\boldsymbol{\Lambda}}_{p,k} = h^{-1} \int_{\operatorname*{supp}\{X\}}(K r_p)(X_{h,j}) X_{h,i}^{p+k} f(X_i) dX_i = \int_0^\infty (K r_p)(u) u^{p+k} f(-uh)du,\] whereas if $\operatorname*{supp}\{X\} = (-\infty,0]$ and $x = 0$, then \[\tilde{\boldsymbol{\Lambda}}_{p,k} = \int_{-\infty}^0 (K r_p)(u) u^{p+k} f(-uh)du.\] For the remainder of this section, the notation is left generic.
For the proofs (Section (ref)) we will frequently abbreviate $s = \sqrt{n h}$.
To define the estimator $\hat{m}$ of $m$ and the bias correction, begin by defining:
where $\operatorname*{diag}(a_i:i = 1, \ldots, n)$ denote the $n\times n$ diagonal matrix constructed using the elements $a_1, a_2, \cdots, a_n$. Note that in the main text $\boldsymbol{\Lambda}_{p,1}$ is denoted by $\boldsymbol{\Lambda}_p$.
Similarly, define
These are identical, but substituting $q$, $L$, and $b$ in place of $p$, $K$, and $h$, respectively. Note that some dimensions change but other do not: for example, $W_p$ and $W_q$ are both $n \times n$, but $\boldsymbol{\Gamma}_p$ is $(p+1)$ square whereas $\boldsymbol{\Gamma}_q$ is $(q+1)$.
Denote by $e_0$ the $(p+1)$-vector with a one in the first position and zeros in the remaining and $Y = (Y_1, \cdots, Y_n)'$. The local polynomial estimator of $m(x)=\mathbb{E}[Y_i|X_i=x]$ is \[\hat{m} = e_0' \boldsymbol{\hat{\beta}}_p = e_0' H_p \boldsymbol{\Gamma}_p^{-1} R_p' W_p Y / n,\] where \[\boldsymbol{\hat{\beta}}_p = \operatorname*{arg\,min}_{b \in \mathbb{R}^{p+1}} \frac{1}{n h} \sum_{i=1}^n ( Y_i - r_p(X_i - x)'b)^2 K \left( X_{h,i} \right) = H_p \boldsymbol{\Gamma}_p^{-1} R_p' W_p Y / n. \] If we define $\check{R} = \left[ r_p( X_1 - x), \cdots, r_p( X_n - x ) \right]'$ and $M = [m(X_1), \ldots, m(X_n)]'$, then we can split $\hat{m} - m$ into the variance and bias terms \[\hat{m} - m = e_0' \boldsymbol{\Gamma}_p^{-1} R_p' W_p (Y-M) / n + e_0' \boldsymbol{\Gamma}_p^{-1} R_p' W_p (M - \check{R}\beta_p) / n.\] This will be useful in the course of the proofs.
The conditional bias is given by
(Recall that in the main paper, $\boldsymbol{\Lambda}_{p,1}$ is denoted $\boldsymbol{\Lambda}_p$.) This result is valid for $p$ odd, our main focus, but also for $p$ even at boundary points.
Denote by $e_{p+1}$ the $(q+1)$-vector with one in the $p+2$ position, and zeros in the rest. Then we estimate the bias as \[\hat{B}_m = h^{p+1} \hat{m}^{(p+1)} \frac{1}{(p+1)!} e_0' \boldsymbol{\Gamma}_p^{-1} \boldsymbol{\Lambda}_{p,1}, \qquad \text{ where } \qquad \hat{m}^{(p+1)} = [(p+1)!] e_{p+1}'H_q \boldsymbol{\Gamma}_q^{-1} R_q' W_q Y / n.\] The bias corrected estimator can then be written
using the fact that $e_{p+1}'H_q = b^{p+1} e_{p+1}'$.
The fixed-$n$ variances are
and
where \[\Sigma = \operatorname*{diag}(v(X_i): i = 1,\ldots, n), \qquad \text{ with } \qquad v(x) = \mathbb{V}[Y \vert X = x].\]
These are the closest analogue to the density case, but are still random due to the conditioning on the covariates. Their respective estimators are
and
The conditional variance matrixes are estimated as \[\boldsymbol{\hat{\Sigma}}_p = \operatorname*{diag}(\hat{v}(X_i): i = 1,\ldots, n), \qquad \text{ with } \qquad \hat{v}(X_i) = ( Y_i - r_p(X_i - x)'\boldsymbol{\hat{\beta}}_p )^2,\] and \[\boldsymbol{\hat{\Sigma}}_q = \operatorname*{diag}(\hat{v}(X_i): i = 1,\ldots, n), \qquad \text{ with } \qquad \hat{v}(X_i) = ( Y_i - r_q(X_i - x)'\boldsymbol{\hat{\beta}}_q )^2.\]
The Studentized statistics of interest are then: \[T_\mathtt{us} = \frac{\sqrt{nh}( \hat{m} - m)}{\hat{\sigma}_\mathtt{us}}, \qquad T_\mathtt{bc} = \frac{\sqrt{nh}( \hat{m} - \hat{B}_m - m)}{\hat{\sigma}_\mathtt{us}}, \qquad T_\mathtt{rbc} = \frac{\sqrt{nh}( \hat{m} - \hat{B}_m - m)}{\hat{\sigma}_\mathtt{rbc}}.\] The main result of this section is an Edgeworth expansion of the distribution function of these statistics.
The terms of the Edgeworth expansion require further notation and discussion. The expressions are not nearly as compact as in the density case (cf. Section (ref)).
Define the expectations of $\boldsymbol{\Gamma}_p$, $\boldsymbol{\Gamma}_q$, $\boldsymbol{\Lambda}_{p,k}$, and $\boldsymbol{\Lambda}_{q,k}$ as $\tilde{\boldsymbol{\Gamma}}_p$, $\tilde{\boldsymbol{\Gamma}}_q$, $\tilde{\boldsymbol{\Lambda}}_{p,k}$, and $\tilde{\boldsymbol{\Lambda}}_{q,k}$, such as \[\tilde{\boldsymbol{\Gamma}}_p = \mathbb{E} \left[ \boldsymbol{\Gamma}_p \right] = \mathbb{E} \left[ h^{-1} (K r_p r_p')(X_{h,i}) \right]. \] These will be used to define nonrandom biases and variances that appear in the expansions.
The biases are defined in Eqn.\ (ref), and are given by
Further discussion and leading terms are found in Section (ref).
The fixed-$n$ variances are computed conditionally, and we must replace them with their nonrandom analogues (just as $\eta_\mathtt{us}$ and $\eta_\mathtt{bc}$ must be nonrandom). Recalling Equations (ref) and (ref), define
where \[ \boldsymbol{\tilde{\Psi}}_p = \mathbb{E} \left[ \boldsymbol{\check{\Psi}}_p \right] \qquad \text{ and } \qquad \boldsymbol{\check{\Psi}}_p := h R_p' W_p \Sigma W_p R_p / n, \] and
where \[ \boldsymbol{\tilde{\Psi}}_q = \mathbb{E} \left[ \boldsymbol{\check{\Psi}}_q \right] \quad \text{and} \quad \boldsymbol{\check{\Psi}}_q := h \left( R_p' W_p - \rho^{p+1} \tilde{\boldsymbol{\Lambda}}_{p,1} \tilde{\boldsymbol{\Gamma}}_q^{-1} R_q' W_q\right) \Sigma \left( R_p' W_p / n - \rho^{p+1} \tilde{\boldsymbol{\Lambda}}_{p,1} \tilde{\boldsymbol{\Gamma}}_q^{-1} R_q' W_q / n \right)' . \] In the course of the proofs, we will also use $\boldsymbol{\hat{\Psi}}_p = h R_p' W_p \boldsymbol{\hat{\Sigma}}_p W_p R_p / n$ and the analogously-defined $\boldsymbol{\hat{\Psi}}_q$.
We now give the precise forms of the polynomials in the Edgeworth expansion. As with the density, there will be both even and odd polynomials. These are not as compact or simple as the density case. Further, we will not attempt to simplify these functions by making use of limiting versions of moments. For example, we will not replace $\tilde{\boldsymbol{\Lambda}}_{p,1}$ by $f(x) \int (K r_p)(u) u^{p+1}du$, and similarly for other pieces. The only simplification made will be the use of $q_{k,\mathtt{us}}(z)$ in the expansion for $T_\mathtt{bc}$, which otherwise would require further notation than what is below (along the lines of $p_{1,\mathtt{us}}(z)$ below).
First, define the following functions, which depend on $n$, $p$, $q$, $h$, $b$, $K$ and $L$, but this is generally suppressed:
With this notation, we can write
and
We will define the Edgeworth expansion polynomials first for the undersmoothing case. The standard Normal density is $\phi(z)$. First, the even polynomials are \[p_{1,\mathtt{us}}(z) = \phi(z) \tilde{\sigma}_\mathtt{us}^{-3} \mathbb{E} \left[ h^{-1} \ell^0_\mathtt{us}(X_i)^3 \varepsilon_i^3 \right] \left\{ (2z^2 - 1)/6 \right\} \] and \[p_{3,\mathtt{us}}(z) = - \phi(z) \tilde{\sigma}_\mathtt{us}^{-1}.\] The absence of $p^{(2)}(z)$ is noteworthy: there is no version of this term for local polynomial estimation, because $\varepsilon_i$ is conditionally mean zero.
Next, the odd polynomials for undersmoothing are defined as follows:
\[q_{2,\mathtt{us}}(z) = - \phi(z) \tilde{\sigma}_\mathtt{us}^{-2} z / 2 ; \] \[q_{3,\mathtt{us}}(z) = \phi(z) \tilde{\sigma}_\mathtt{us}^{-4} \mathbb{E} [ h^{-1} \ell^0_\mathtt{us}(X_i)^3 \varepsilon_i^3 ] ( z^3 / 3 ). \] For robust bias correction, both the even polynomials, $p_{1,\mathtt{rbc}}(z)$ and $p_{3,\mathtt{rbc}}(z)$, and the odd polynomials, $q_{1,\mathtt{rbc}}(z)$, $q_{2,\mathtt{rbc}}(z)$, and $q_{3,\mathtt{rbc}}(z)$ are defined in the exact same way, but changing the $\tilde{\sigma}_\mathtt{us}$ to $\tilde{\sigma}_\mathtt{rbc}$, $\ell^k_\mathtt{us}(\cdot)$ to $\ell^k_\mathtt{bc}(\cdot)$, $K$ to $L$, and $p$ to $q$, and so forth. For $q_{1,\mathtt{us}}(z)$ and $q_{1,\mathtt{rbc}}(z)$, the seventh term can be rewritten by rearranging the terms and factoring the expectation, as follows:
The polynomials defined here are for distribution function expansions, and are different from those used for coverage error. The polynomials $q_{1,\mathtt{us}}$, $q_{2,\mathtt{us}}$, and $q_{3,\mathtt{us}}$ and $q_{1,\mathtt{rbc}}$, $q_{2,\mathtt{rbc}}$, and $q_{3,\mathtt{rbc}}$, which do not have an argument, used for coverage error in the main text and in Corollary (ref) below, are defined in terms of those given above, which do have an argument. Specifically, the polynomials above should be doubled, divided by the standard Normal density, and evaluated at the Normal quantile $z_{\alpha/2}$, that is, \[ q_{k,\boldsymbol{\bullet}} := \left. \frac{2}{\phi(z)} q_{k,\boldsymbol{\bullet}}(z) \right|_{z = z_{\alpha/2} }, \qquad\qquad k=1,2,3, \quad \boldsymbol{\bullet} = \mathtt{us}, \mathtt{rbc} \]
For traditional bias correction, $q_{1,\mathtt{us}}(z)$, $q_{2,\mathtt{us}}(z)$, and $q_{3,\mathtt{us}}(z)$ are used, but such simplification can not be done for $p_{1,\mathtt{bc}}(z)$ and $p_{3,\mathtt{bc}}(z)$, which must be defined as
and \[p_{3,\mathtt{bc}}(z) = - \phi(z) \tilde{\sigma}_\mathtt{us}^{-1}.\]
Lastly, traditional bias correction also exhibits additional terms in the expansion (see discussion in the main text) representing the covariance of $\hat{m}$ and $\hat{B}_m$ (denoted by $\Omega_{1,\mathtt{bc}}$) and the variance of $\hat{B}_m$ ($\Omega_{2,\mathtt{bc}}$). We now state their precise forms. These arise from the mismatch between the variance of the numerator of $T_\mathtt{bc}$ and the standardization used, $\sigma_\mathtt{us}^2$, but these are random, and so $\Omega_{1\mathtt{bc}}$ and $\Omega_{2,\mathtt{bc}}$ must be derived from the nonrandom versions, $\tilde{\sigma}_\mathtt{rbc}^2$ and $\tilde{\sigma}_\mathtt{us}^2$ (cf. Section (ref); for the same reason $\eta_\mathtt{us}$ and $\eta_\mathtt{bc}$ must be nonrandom). Recalling the definitions above,
Therefore \[\Omega_{1,\mathtt{bc}} = -2 \tilde{\sigma}_\mathtt{us}^{-2} \mathbb{E}[ h^{-1} \{\rho^{-p-2} \ell^0_\mathtt{us}(X) (\ell^0_\mathtt{bc}(X) - \ell^0_\mathtt{us}(X))\} v(X) ] \] and \[\Omega_{2,\mathtt{bc}} = \tilde{\sigma}_\mathtt{us}^{-2} \mathbb{E}[ b^{-1} \{\rho^{-p-2}(\ell^0_\mathtt{bc}(X) - \ell^0_\mathtt{us}(X))\}^2 v(X) ] . \]
In the main text we give a direct plug-in (DPI) rule to implement the coverage-error optimal bandwidth. Here we we give complete details for this procedure as well as document a second practical choice, based on a rule-of-thumb (ROT) strategy. Both choices yield the optimal coverage error decay rate at interior and boundary points.
All our methods are implemented in {\sf R} and {\tt STATA} via the {\tt nprobust} package, available from \url{http://sites.google.com/site/nppackages/nprobust} (see also \url{http://cran.r-project.org/package=nprobust}). See Calonico-Cattaneo-Farrell2017_nprobust for a complete description.
As in the density case, the MSE-optimal bandwidth undercovers when used in the undersmoothing confidence interval; that is, Remark (ref) applies directly. See also Hall-Horowitz2013_AoS.
As with the density case, a simple rule-of-thumb based on rescaling the MSE-optimal bandwidth is: \[\hat{h}^{\tt int}_\mathtt{rot} = \hat{h}^{\tt int}_\mathtt{mse} \; n^{-(p-1)/((2p+3)(p+4))} \qquad \text{ and } \qquad \hat{h}^{\tt bnd}_\mathtt{rot} = \hat{h}^{\tt bnd}_\mathtt{mse} \; n^{-p/((2p+3)(p+3))}.\] where $\hat{h}^{\tt int}_\mathtt{mse}$ and $\hat{h}^{\tt bnd}_\mathtt{mse}$ denote readily-available implementations of the MSE-optimal bandwidth for interior and boundary points, respectively. See, e.g., Fan-Gijbels1996_book. Again, when $p=1$ in the interior, no scaling is needed ($\hat{h}^{\tt int}_\mathtt{rot} = \hat{h}^{\tt int}_\mathtt{mse}$), but for $p>1$ any data-driven MSE-optimal bandwidth should always be shrunk to improve inference at the boundary (i.e., reduce coverage errors of the robust bias-corrected confidence intervals).
The ROT selector may be especially attractive for simplicity, if estimating the constants described below in the DPI case is prohibitive.
Remark (ref) applies to this case as well, though less transparently and without consequences that are as dramatic.
We now detail the required steps to implement the plug-in bandwidth $\hat{h}_{\mathtt{dpi}}$ for interior and boundary points. We always set $K=L$, $\rho=1$, and $q=p+1$. The steps are:
As argued in the main text, using variance forms other than (ref) and (ref) can be detrimental to coverage. Within these forms however, two alternative estimates of $\Sigma$ are natural. First, motivated by the fact that the least-squares residuals are on average too small, the well-known HC$k$ class of heteroskedasticity consistent estimators can be used; see MacKinnon2013_BookChap for details and a recent review. In our notation, these are defined as follows. First, $\hat{\sigma}_\mathtt{us}^2$-HC0 is the estimator above. Then, for $k=1, 2, 3$, the $\hat{\sigma}_\mathtt{us}^2$-HC$k$ estimator is obtained by dividing $\hat{\varepsilon}_i^2$ by, respectively, $(n-2\operatorname*{tr}(Q_p)+\operatorname*{tr}(Q_p'Q_p))/n$, $(1-Q_{p,ii})$, and $(1-Q_{p,ii})^2$, where $Q_{p,ii}$ is the $i$-th diagonal element of the projection matrix $Q_p:=R_p '\boldsymbol{\Gamma}_p^{-1} R_p' W_p/n$. The corresponding estimators $\hat{\sigma}_\mathtt{rbc}^2$-HC$k$ are the same way, with $q$ in place of $p$. As is well-known in the literature, these estimators perform better for small sample sizes, a fact we confirm in our simulation study below.
A second option is to use a nearest-neighbor-based variance estimators with a fixed number of neighbors, following the ideas of Muller-Stadtmuller1987_AoS,Abadie-Imbens2008_AdES. To define these, let $J$ be a fixed number and $j(i)$ be the $j$-th closest observation to $X_i$, $j=1, \ldots, J$, and set $ \hat{v}(X_i) = \frac{J}{J+1} ( Y_i - \sum_{j=1}^J Y_{j(i)} / J )^2$. This “estimate” is unbiased (but inconsistent) for $v(X_i)$.
Both types of residual estimators could be handled in our results. The constants will change, but the rates will not. This is because, in all cases, the errors in estimating $v(X_i)$ are no greater than in the original $\hat{m}(x)$. Inspection of the proof shows that simple modifications allow for the HC$k$ estimators: only the terms of Eqn.\ (ref) will change, and indeed, we conjecture that the HC$k$ estimators will result in fewer terms and a reduced coverage error. This is consistent with the improved finite-sample behavior of these estimators and the fact that they are asymptotically equivalent. Accommodating the nearest-neighbor estimates require slightly more work and a modified version of Assumption (ref).
One crucial property of our method, in the context of Edgeworth expansions, is that the bias in estimation of $\Sigma$ is of the same order as the original $\hat{m}(x)$. Using other methods may result in additional terms, with possibly distinct rates, appearing in the Edgeworth expansions. Some examples that may have this issue are (i) using $\hat{v}(X_i) = ( Y_i - \hat{m}(x) )^2$; (ii) using local or assuming global heteroskedasticity; (iii) using other nonparametric estimators for $v(X_i)$, relying on new tuning parameters.
The following assumptions are sufficient for our results. The first two are copied directly from the main text (see discussion there) and the third is the appropriate Cram\'er's condition.
The random variables of Assumption (ref) are defined follows. For two kernels $K_1$ and $K_2$, two polynomial orders (i.e. positive integers) $p_1$ and $p_2$, a bandwidth $b$, and a scalar $\rho$, let
and
The subscripts are intended to make clear that $Z_m(\cdot)$ collects quantities from the numerator of the Studentized statistic, while $Z_\sigma(\cdot)$ gathers additional variables required for the variance estimation. With this notation, we define \[ Z_\mathtt{us}(u) = \bigl( Z_m(u; K, p, p, h, 1)' , \ Z_\sigma(u; K, K, p, p, h, 1)' \bigr)', \] \[ Z_\mathtt{bc}(u) = \bigl( Z_m(u; K, p, p+1, h, 1)' , \ Z_m(u; L, q, q, b, \rho)' , \ \operatorname*{vech}(K(u) r_p(u) u^{p+1})' , \ Z_\sigma(u; K, K, p, p, h, 1)' \bigr)', \] and
{\bf Discussion.} This notation is quite compact, and while it emphasizes the simplicity of Cram\'er's condition and the fact that it puts mild restrictions on the kernels, it does obscure the full notational breadth, particularly for $Z_\mathtt{rbc}$. I is also mostly repetitive: what holds for the kernel $K$ and order $p$ fit must also hold for $L$ and $q$, and for their squares and cross products. To make this clear, we can expand all the $Z_m$ and $Z_\sigma$, to write out the full random variables as
and
Finally, the precise random variables $Z_\mathtt{us}(u)$, $Z_\mathtt{bc}(u)$, and $Z_\mathtt{rbc}(u)$ used can be replaced with slightly different constructions without altering the conclusions of Theorem (ref): there are other potential functions $\tilde{T}$ that satisfy Eqn.\ (ref) in the proof. Such changes necessarily involve asymptotically negligible terms, and do not materially alter the severity of the restrictions imposed.
We will not present a detailed discussion of bias issues, along the lines of Section (ref), for brevity; we focus only on the case of nonbinding smoothness.
The biases $\eta_\mathtt{us}$ and $\eta_\mathtt{bc}$ are not as conceptually simple as in the density case. The closest parallel to the density case would be (for example) $\eta_\mathtt{us} = \sqrt{nh} (\mathbb{E}[\hat{m}] - m)$, but this can not be used due to the presence of $\boldsymbol{\Gamma}_p^{-1}$ inside the expectation, and the next natural choice, the conditional bias $\sqrt{nh} (\mathbb{E}[\hat{m} \vert X_1, \ldots X_n] - m)$, is still random. Instead, $\eta_\mathtt{us}$ and $\eta_\mathtt{bc}$ are biases computed after replacing $\boldsymbol{\Gamma}_p$, $\boldsymbol{\Gamma}_q$, and $\boldsymbol{\Lambda}_{p,1}$ with their expectations, denoted $\tilde{\boldsymbol{\Gamma}}_p$, $\tilde{\boldsymbol{\Gamma}}_q$, and $\tilde{\boldsymbol{\Lambda}}_{p,1}$. We thus define
For the generic results of coverage error or the generic Edgeworth expansions of Theorem (ref) below, the above definitions of $\eta_\mathtt{us}$ and $\eta_\mathtt{bc}$ are suitable. For the Corollaries detailing specific cases, and to understand the behavior at different points, it is useful to make the leading terms precise, that is, analogues of Equations (ref) and (ref). We must consider interior and boundary point estimation, and even and odd $q$. We depart slightly from other terms of the expansion in that we do retain only the leading term for some pieces. This is done in order to capture the rate of convergence explicitly and to give practicable results. These results are derived by Fan-Gijbels1996_book and similar calculations (though our expressions differ slightly as fixed-$n$ expectations are retained as much as possible).
Since $p$ is odd, both at boundary and interior points we have \[\eta_\mathtt{us} = \sqrt{nh} h^{p+1} \frac{ m^{(p+1)} } { (p+1)! } e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1} \tilde{\boldsymbol{\Lambda}}_{p,1} \left[1 + o(1) \right].\]
Moving to $\eta_\mathtt{bc}$, consider the first term, which in the present notation is: $ \sqrt{n h} \mathbb{E}[ h^{-1} \ell^0_\mathtt{us}(X) (m(X) - r_{p+1}(X - x)' \beta_{p+1})]$. With $p+1$ even, we find that in the interior the leading terms are \[ \sqrt{n h} h^{p+3} e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1}\left( \frac{ m^{(p+2)} } { (p+2)! } \tilde{\boldsymbol{\Lambda}}_{p,2} h^{-1} + \frac{ m^{(p+3)} } { (p+3)! } \tilde{\boldsymbol{\Lambda}}_{p,3} \right) \left[1 + o(1) \right],\] due to the well-known symmetry properties of local polynomials that result in the cancellation of the leading terms of $\tilde{\boldsymbol{\Gamma}}_p^{-1}$ and $\tilde{\boldsymbol{\Lambda}}_{p,2}$. The rate of $h^{p+3}$ accounts for this. At the boundary, no such cancellation occurs and we have only \[ \sqrt{n h} h^{p+2} \frac{ m^{(p+2)} } { (p+2)! } e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1} \tilde{\boldsymbol{\Lambda}}_{p,2} \left[1 + o(1) \right].\] Next, turn to the bias of the bias estimate: \[\sqrt{nh} \rho^{p+1} e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1} \tilde{\boldsymbol{\Lambda}}_{p,1} e_{p+1}' \tilde{\boldsymbol{\Gamma}}_q^{-1} \int L(u) r_q(u) \left( m(x - u b) - r_q(ub)'\beta_q \right) f(x - ub)du.\] If $q$ is odd (so that $q-(p+1)$ is also odd), then at the interior or boundary the leading term will be \[\sqrt{nh} b^{q+1} \rho^{p+1} \frac{ m^{(q+1)} } { (q+1)! } e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1} \tilde{\boldsymbol{\Lambda}}_{p,1} e_{p+1}' \tilde{\boldsymbol{\Gamma}}_q^{-1} \tilde{\boldsymbol{\Lambda}}_{q,1} \left[1 + o(1) \right] \asymp \sqrt{nh}h^{p+1} b^{q-p}.\] The same expression applies at the boundary for $q$ even. However, for the interior, if $q$ is even, which it is in the leading case of $q=p+1$, then we again have cancellation of certain leading terms, resulting in the bias of the bias estimate being \[\sqrt{nh} b^{q+2} \rho^{p+1} e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1} \tilde{\boldsymbol{\Lambda}}_{p,1} e_{p+1}' \tilde{\boldsymbol{\Gamma}}_q^{-1} \left( \frac{ m^{(q+1)} } { (q+1)! } \tilde{\boldsymbol{\Lambda}}_{q,1} b^{-1} + \frac{ m^{(q+2)} } { (q+2)! } \tilde{\boldsymbol{\Lambda}}_{q,2} \right) \left[1 + o(1) \right] \asymp \sqrt{nh}h^{p+1} b^{q+1-p}.\] Combining all these results, we find the following. For an interior point $\eta_\mathtt{bc}^{\tt int} = \sqrt{n h} h^{p+3} \left[ \tilde{\eta}_\mathtt{bc}^{\tt int} + o(1)\right]$, where, if $q$ is even
while if $q$ is odd, \[\tilde{\eta}_\mathtt{bc}^{\tt int} = e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1}\biggl( \frac{ m^{(p+2)} } { (p+2)! } \tilde{\boldsymbol{\Lambda}}_{p,2} h^{-1} + \frac{ m^{(p+3)} } { (p+3)! } \tilde{\boldsymbol{\Lambda}}_{p,3} \biggr) - \rho^{-2} b^{q-(p+2)} \frac{ m^{(q+1)} } { (q+1)! } e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1} \tilde{\boldsymbol{\Lambda}}_{p,1} e_{p+1}' \tilde{\boldsymbol{\Gamma}}_q^{-1} \tilde{\boldsymbol{\Lambda}}_{q,1} .\] At the boundary, for any $q$, $\eta_\mathtt{bc}^{\tt bnd} = \sqrt{n h} h^{p+2} \left[ \tilde{\eta}_\mathtt{bc}^{\tt bnd} + o(1)\right]$, with \[ \tilde{\eta}_\mathtt{bc}^{\tt bnd} = \frac{ m^{(p+2)} } { (p+2)! } e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1} \tilde{\boldsymbol{\Lambda}}_{p,2} - \rho^{-1}b^{q-(p+1)} \frac{ m^{(q+1)} } { (q+1)! } e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1} \tilde{\boldsymbol{\Lambda}}_{p,1} e_{p+1}' \tilde{\boldsymbol{\Gamma}}_q^{-1} \tilde{\boldsymbol{\Lambda}}_{q,1} .\]
We now state our generic Edgeworth expansion, from whence the coverage probability expansion results follow immediately. We have opted to state separate results for undersmoothing, bias correction, and robust bias correction, rather than the unified statement of Theorem (ref), for clarity. The unified structure is still present, and will be used in the proof of the result below, but is too cumbersome to use here. The Standard Normal distribution and density functions are $\Phi(z)$ and $\phi(z)$, respectively.
For undersmoothing estimators, we have the following result, which is valid for both interior and boundary points, with moments appropriately truncated if necessary. This result is the analogue of the robust bias correction corollary in the main text, and follows directly from the generic theorem there or Theorem (ref) above. Exponents such as $1 + 2(p+1)$ are intentionally not simplified to ease comparison to other results, particularly the density case.
The polynomials $q_{1,\mathtt{us}}$, $q_{2,\mathtt{us}}$, and $q_{3,\mathtt{us}}$, which do not have an argument, are defined in terms of those given in Section (ref) and used in Theorem (ref), which do have an argument. Specifically, the polynomials in Section (ref) and Theorem (ref) should be doubled, divided by the standard Normal density, and evaluated at the Normal quantile $z_{\alpha/2}$, that is, \[ q_{k,\mathtt{us}} := \left. \frac{2}{\phi(z)} q_{k,\mathtt{us}}(z) \right|_{z = z_{\alpha/2} }, \qquad\qquad k=1,2,3. \]
We will first prove Theorem (ref)(a), as it is notationally simplest. From a technical and conceptual point of view, proving the remainder of Theorem (ref) is identical, simply more involved notationally due to the additional complexity of the bias correction. Outlines of these proofs are found below.
Let $s = \sqrt{n h}$.
Throughout this proof, we will generally omit the subscripts $\mathtt{us}$ and $p$ when this causes no confusion. This entire proof focuses on the undersmoothing statistic, $T_\mathtt{us} = \hat{\sigma}_\mathtt{us}^{-1} s(\hat{m} - m)$, and since bias correction is not involved at all, the associated constructions such as $\boldsymbol{\Gamma}_q$, $W_q$, etc, do not appear, and hence there is no need to carry the additional notation to distinguish $W_p$ from $W_q$, or $\hat{\sigma}_\mathtt{us}$ from $\hat{\sigma}_\mathtt{rbc}$, for example, and we will simply write $\boldsymbol{\Gamma}$ for $\boldsymbol{\Gamma}_p$, $W$ for $W_p$, $\hat{\sigma}$ for $\hat{\sigma}_\mathtt{us}$, etc.
Our goal is to expand $\mathbb{P}[T_\mathtt{us}<z]$, where $T_\mathtt{us} = \hat{\sigma}^{-1} s(\hat{m} - m)$. The proof proceeds by identifying a smooth function $\tilde{T} = \tilde{T}(z)$ such that, for the random variable $Z_\mathtt{us} := Z_\mathtt{us}(u)$ that obeys Cram\'er's condition (Assumption (ref)), $\tilde{T}(\mathbb{E}[Z_\mathtt{us}]) = 0$ and
where $\bar{Z} = \sum_{i=1}^n Z_i/n$ and $\tilde{z}$ is a known, nonrandom quantity that depends on the original quantile $z$ and the remainder $T_\mathtt{us} - \tilde{T}$. An Edgeworth expansion for $\tilde{T}$ holds under Assumption (ref), and a Taylor expansion of this function around $\tilde{z}$ yields the final result. As in the density case, $\tilde{z}$ will capture the bias terms of $T_\mathtt{us}$: in that case $\tilde{z} = z - \eta/\tilde{\sigma}$, but here bias is present in both the numerator and the Studentization.
To begin, define the notation $\check{R} = \left[ r_p( X_1 - x), \cdots, r_p( X_n - x ) \right]'$ and $M = [m(X_1), \ldots, m(X_n)]'$, and use this to split $T$ into variance and bias terms, as follows: \[T = \hat{\sigma}^{-1} s e_0' \boldsymbol{\Gamma}^{-1} R'W(Y-M)/n + \hat{\sigma}^{-1} s e_0' \boldsymbol{\Gamma}^{-1} R'W(M - \check{R}\beta)/n.\] We use this decomposition to rewrite $\mathbb{P}[ T_\mathtt{us} < z ]$ as
The first three lines in the last equality obey the desired properties of $\tilde{T}$ by the orthogonality of $\varepsilon_i$, the definition of $\eta_\mathtt{us}$ in Eqn.\ (ref) as $\mathbb{E}\left[ s e_0' \tilde{\boldsymbol{\Gamma}}^{-1} R'W(M - \check{R}\beta)/n \right]$, and the fact that $\boldsymbol{\Gamma}^{-1} - \tilde{\boldsymbol{\Gamma}}^{-1} = \tilde{\boldsymbol{\Gamma}}^{-1} \left( \tilde{\boldsymbol{\Gamma}} - \boldsymbol{\Gamma} \right) \boldsymbol{\Gamma}^{-1}$. For the final two (which are $T_\mathtt{us} - \tilde{\sigma}^{-1} s (\hat{m} - m) = \hat{\sigma}^{-1} - \tilde{\sigma}^{-1} s (\hat{m} - m) $), we must expand the difference $\hat{\sigma}^{-1} - \tilde{\sigma}^{-1}$. Accounting for the resulting terms will constitute the bulk of the remainder of the proof, as well as complete the construction of $\tilde{z}$ and the remainder terms of Eqn.\ (ref).\footnote{Technically, to obtain a $\tilde{T}$ with the desired properties, one need not expand $\hat{\sigma}^{-1} - \tilde{\sigma}^{-1}$ for the variance term: that is, in Eqn.\ (ref), $\tilde{\sigma}^{-1} s e_0' \boldsymbol{\Gamma}^{-1} R'W(Y-M)/n$ and $\left( \hat{\sigma}^{-1} - \tilde{\sigma}^{-1} \right) s e_0' \boldsymbol{\Gamma}^{-1} R'W(Y-M)/n$ may be collapsed. This requires strengthening Cram\'er's condition (see Section (ref)), and since $\hat{\sigma}^{-1} - \tilde{\sigma}^{-1}$ must be accounted for in the final bias term, $\left( \hat{\sigma}^{-1} - \tilde{\sigma}^{-1} \right) s e_0' \boldsymbol{\Gamma}^{-1} R'W(M - \check{R}\beta)/n$, there is little reason not to do both terms.}
To begin, with $\tilde{\sigma}^2 = e_0 ' \tilde{\boldsymbol{\Gamma}}^{-1} \boldsymbol{\tilde{\Psi}} \tilde{\boldsymbol{\Gamma}}^{-1} e_0$ defined in Section (ref),
and hence a Taylor expansion gives \[\frac{1}{\hat{\sigma}} = \frac{1}{\tilde{\sigma}} \left[ 1 - \frac{1}{2} \frac{ \hat{\sigma}^2 - \tilde{\sigma}^2}{\tilde{\sigma}^2} + \frac{3}{8} \left( \frac{ \hat{\sigma}^2 - \tilde{\sigma}^2}{\tilde{\sigma}^2} \right)^2 - \frac{1}{3!} \frac{15}{8} \left( \frac{ \hat{\sigma}^2 - \tilde{\sigma}^2}{\tilde{\sigma}^2} \right)^3 \frac{\tilde{\sigma}^7}{\bar{\sigma}^7} \right],\] for a point $\bar{\sigma}^2 \in [\tilde{\sigma}^2, \hat{\sigma}^2]$, and so
We thus focus on $\hat{\sigma}^2 - \tilde{\sigma}^2$. Recall the definition of $\boldsymbol{\check{\Psi}} = h R' W \Sigma W R/n$. Then define the two terms $A_1$ and $A_2$ through the following:
For $A_1$, recall that $\hat{\varepsilon}_i = y_i - r_p(X_i - x)'\boldsymbol{\hat{\beta}}_p$ and so
where
is due to the approximation of the (average over the) conditional variance by the squared residuals (i.e.\ $A_{1,1}$ is the sole remainder that would arise if the true residuals were known and used in place of $\hat{\varepsilon}_i^2$), and, using $r_p(X_i - x)' \boldsymbol{\hat{\beta}} = r_p(X_i - x)' H_p \boldsymbol{\Gamma}^{-1} R'W Y/n = r_p(X_{h,i})' \boldsymbol{\Gamma}^{-1} R'W Y/n$, the terms $A_{1,k}$, $k=2, 3, \ldots, 8$ are:
With this notation, we can write $A_1 = e_0' \boldsymbol{\Gamma}^{-1} \left( \boldsymbol{\hat{\Psi}} - \boldsymbol{\check{\Psi}} \right) \boldsymbol{\Gamma}^{-1} e_0 = e_0' \boldsymbol{\Gamma}^{-1} \left( \sum_{k=1}^8 A_{1,k} \right) \boldsymbol{\Gamma}^{-1} e_0$. The terms $A_{1,1}$ to $A_{1,5}$ will be incorporated into $\tilde{T}$: notice that these terms obey $A_{1,k} = A_{1,k}(\bar{Z}_\mathtt{us})$ and $A_{1,k}(\mathbb{E}[Z_\mathtt{us}]) = 0$, and hence these properties will be inherited in the final two lines of Eqn.\ (ref). However, $A_{1,6}$, $A_{1,7}$, and $A_{1,8}$ do not have these properties, and will thus be incorporated into $\tilde{z}$ and the remainder. Details are below.
Turning to $A_2$ in Eqn.\ (ref), using the identity $\boldsymbol{\Gamma}^{-1} - \tilde{\boldsymbol{\Gamma}}^{-1} = \tilde{\boldsymbol{\Gamma}}^{-1} \left( \tilde{\boldsymbol{\Gamma}} - \boldsymbol{\Gamma} \right) \boldsymbol{\Gamma}^{-1}$ and that $\boldsymbol{\Gamma}$ and $\Psi$ are symmetric, we find that
All of these terms obey the required properties of $\tilde{T}$.
We now collect the terms from expanding $\hat{\sigma}^{-1} - \tilde{\sigma}^{-1}$ and return to Eqn.\ (ref). Plugging the terms $A_{1,1}$--$A_{1,8}$ and $A_2$ into the Taylor expansion in Eqn.\ (ref), by way of Eqn.\ (ref), and collecting terms appropriately (i.e. those that belong in $\tilde{T}$ as described above), we have the following, which picks up from Eqn.\ (ref) and is a precursor to Eqn.\ (ref):
In this statement, we have made the following constructions:
and
In $U$ and $\tilde{z}$, each $\tilde{A}_{1,k} $ is $A_{1,k}$ where all elements have been replaced by their respective fixed-$n$ expected values, that is,
and \[\tilde{A}_{1,8} = \mathbb{E}\left[ h^{-1} (K^2 r_p r_p ' )(X_{h,i}) \mathbb{E}\left[ \left. h^{-1} r_p(X_{h,i}) ' \tilde{\boldsymbol{\Gamma}}^{-1} (K r_p)(X_{h,j}) \left[m(X_j) - r_p(X_j - x)'\beta_p \right] \right\vert X_i \right]^2 \right].\]
The next step in the proof is to show that, for $r_* = \max\{s^{-2}, \eta^2, h^{p+1} \}$ (i.e., the slowest decaying), it holds that
This result is established by Lemma (ref) in Section (ref) below. This, together with Eqn.\ (ref), implies Eqn.\ (ref).
Under Assumption (ref), an Edgeworth expansion holds for $\tilde{T}$ up to $o(s^{-2} + s^{-1}\eta + \eta^2)$. Thus, for a smooth function $G(z)$, we have $\mathbb{P}[\tilde{T} < z] = G(z) + o(s^{-2} + s^{-1}\eta + \eta^2)$. Therefore, a Taylor expansion gives \[ \mathbb{P}[\tilde{T} < \tilde{z}] = G(z) - G^{(1)}(z) \left\{\tilde{\sigma}^{-1} - \frac{1}{2 \tilde{\sigma}^3} e_0' \boldsymbol{\Gamma}^{-1} \left( \tilde{A}_{1,6} + \tilde{A}_{1,7} + \tilde{A}_{1,8} \right) \boldsymbol{\Gamma}^{-1} e_0 \right\} + o(s^{-2} + s^{-1}\eta + \eta^2),\] which together with Eqn.\ (ref) establishes the validity of the Edgeworth expansion. The terms of the expansion are computed in Section (ref) below.\qed
To prove parts (b) and (c) of Theorem (ref) the same steps are required, and so we will not pursue all the details here. Indeed, the same expansions are performed and the same bounds computed on objects which are conceptually similar, only taking into account the bias correction (in the numerator for (b), and also in the denominator for (c)). The bias correction will result in essentially two changes: first, many more terms like $\boldsymbol{\Gamma} - \tilde{\boldsymbol{\Gamma}}$ appear, and second, the bias expressions and rates change. To illustrate, we will list several key points where these changes manifest. This list is not exhaustive, but it will show that the same methods used above still apply.
First, for the numerator of $T_\mathtt{bc}$ and $T_\mathtt{rbc}$, recall that the estimator $\hat{m}$ is \[\hat{m} = \Bigl\{ e_0' \boldsymbol{\Gamma}_p^{-1} R_p' W_p \Bigr\}Y / n, \] while the bias corrected estimator is \[\hat{m} - \hat{B}_m = \Bigl\{ e_0' \boldsymbol{\Gamma}_p^{-1} \left( R_p' W_p - \rho^{p+1} \boldsymbol{\Lambda}_{p,1} e_{p+1}' \boldsymbol{\Gamma}_q^{-1} R_q' W_q \right) \Bigr\}Y / n. \] Comparing these two expressions, it can be seen that the terms in the proof above that involve $\boldsymbol{\Gamma}_p - \tilde{\boldsymbol{\Gamma}}_p$ will now additionally involve $\boldsymbol{\Gamma}_q - \tilde{\boldsymbol{\Gamma}}_q$ and $\boldsymbol{\Lambda}_{p,1} - \tilde{\boldsymbol{\Lambda}}_{p,1}$, whereas those that with $e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1} R_p'W_p$ will now have $e_0' \tilde{\boldsymbol{\Gamma}}_p^{-1} \left( R_p' W_p - \rho^{p+1} \tilde{\boldsymbol{\Lambda}}_{p,1} e_{p+1}' \tilde{\boldsymbol{\Gamma}}_q^{-1} R_q' W_q \right)$ instead. To give a concrete example, consider the third line of Eqn.\ (ref), \[\tilde{\sigma}_\mathtt{us}^{-1} s e_0' \left( \boldsymbol{\Gamma}_p^{-1} - \tilde{\boldsymbol{\Gamma}}_p^{-1} \right) R_p'W_p(M - \check{R}_p\beta_p)/n,\] which becomes a piece of the function $\tilde{T}$. For part (b) Theorem (ref), treating $T_\mathtt{bc}$, this will become
and part (c) will have the same but with $\tilde{\sigma}_\mathtt{rbc}^{-1}$. Then, since
this term is handled identically, since the appropriate Cram\'er's condition is assumed.
Consider now the denominator of the Studentized statistics. For part (b), there is no change as $\hat{\sigma}_\mathtt{us}^2$ is still used, and so the terms involving $A_{1,k}$ and $A_2$ will be identical. However, for $T_\mathtt{rbc}$, we must account for changes of the above form, but also that the residuals are estimated with the degree $q$ fit: $\hat{\varepsilon}_ i = y_i - r_q(X_i - x)'\boldsymbol{\hat{\beta}}_q$ instead of degree $p$. With these changes in mind, the analogue of Eqn.\ (ref) will be
The second term will proceed as above, though $\boldsymbol{\check{\Psi}}_p - \boldsymbol{\tilde{\Psi}}_p$ will be replaced by \[\boldsymbol{\check{\Psi}}_q - \boldsymbol{\tilde{\Psi}}_q = \frac{1}{n h} \sum_{i=1}^n \left\{\tilde{\ell}^0_\mathtt{bc}(X_i) \tilde{\ell}^0_\mathtt{bc}(X_i)' v(X_i) - \mathbb{E} \left[\tilde{\ell}^0_\mathtt{bc}(X_i) \tilde{\ell}^0_\mathtt{bc}(X_i)' v(X_i)\right]\right\}, \] where $\tilde{\ell}^0_\mathtt{bc}(X_i) = (K r_p)(X_{h,i}) - \rho^{p+1} \tilde{\boldsymbol{\Lambda}}_{p,1} \tilde{\boldsymbol{\Gamma}}_q^{-1} (L r_p)(\rho X_{h,i})$ (cf. Section (ref), the function $\ell^0_\mathtt{bc}$ therein is $\ell^0_\mathtt{bc}(X_i) = e_0'\tilde{\boldsymbol{\Gamma}}_p^{-1} \tilde{\ell}^0_\mathtt{bc}(X_i)$). To use similar notation, \[\boldsymbol{\check{\Psi}}_p - \boldsymbol{\tilde{\Psi}}_p = \frac{1}{n h} \sum_{i=1}^n \left\{\tilde{\ell}^0_\mathtt{us}(X_i) \tilde{\ell}^0_\mathtt{us}(X_i)' v(X_i) - \mathbb{E} \left[\tilde{\ell}^0_\mathtt{us}(X_i) \tilde{\ell}^0_\mathtt{us}(X_i)' v(X_i)\right]\right\}.\] Then, expanding $\tilde{\ell}^0_\mathtt{bc}(X_i)$ shows that $\boldsymbol{\check{\Psi}}_q - \boldsymbol{\tilde{\Psi}}_q$ is equal to
and since all these terms still obey the appropriate Cram\'er's condition, the same steps apply. (The extra factor of $\rho$ in $\rho^{2(p+1)+1}$ and $\rho^{(p+1) + 1}$ accounts for the fact that $\hat{\sigma}_\mathtt{rbc}^2$ is scaled by $(nh)$ instead of $(nb)$, but the $W_q$ matrixes contribute a $b^{-1}$.)
The first term of Eqn.\ (ref) will also follow by the same method as in the prior proof, but more care must be taken as many more terms will be present because $\boldsymbol{\hat{\Psi}}_q - \boldsymbol{\check{\Psi}}_q$ consists of the following three terms, representing the variance of $\hat{m}$, the variance of $\hat{B}_m$, and their covariance, respectively:
The first of these three is as in the prior proof, and yields the same $A_{1,1}$--$A_{1,8}$, only with the bias of a $q$-degree fit: $m(X_i) - r_q(X_i - x)'\beta_q$. If we define \[\check{\boldsymbol{\check{\Psi}}}_q := \frac{1}{nb} \sum_{i=1}^n (L^2 r_q r_q')(X_{b,i}) v(X_i)\] then the second term of $\boldsymbol{\hat{\Psi}}_q - \boldsymbol{\check{\Psi}}_q$ is equal to
The first of these terms will also give rise to versions of $A_{1,1}$--$A_{1,8}$, only with the bias of a $q$-degree fit and changing $K$ to $L$, $p$ to $q$, $h$ to $b$, etc, and will thus be treated exactly as above. The rest of these are incorporated into $\tilde{T}_\mathtt{rbc}$, similar to how $A_2$ is treated, because Cram\'er's condition is satisfied. The third and final piece of $\boldsymbol{\hat{\Psi}}_q - \boldsymbol{\check{\Psi}}_q$ is equal to
and thus is entirely analogous, with yet another version of $A_{1,1}$--$A_{1,8}$ defined for the remainder in the first line, and the second two easily incorporated into $\tilde{T}_\mathtt{rbc}$.
From these arguments, it is clear that the analogue of Lemma (ref) will hold for these cases as well: the same fundamental pieces are involved, and thus the same arguments will apply, just as above.
Our proof of Theorem (ref) relies on the following lemmas. The first gives generic results used to derive rate bounds on the probability of deviations of the necessary terms. Some such results are collected in Lemma (ref). Lemma (ref) shows how to use the previous results to establish negligibility of the remainder terms required for Eqn.\ (ref).
As above, we will generally omit the details required for Theorem (ref) parts (b) and (c), to save space. These are entirely analogous, as can be seen from the steps in Lemma (ref). Indeed, the first results are stated in terms of the kernel $K$ and bandwidth $h$, but continue to hold for $L$ and $b$ under the obvious substitutions and appropriate assumptions.
Throughout proofs $C$ shall be a generic conformable constant that may take different values in different places. If more than one constant is needed, $C_1$, $C_2$, \ldots, will be used.
To illustrate how the above Lemma is used for the objects under study, we present the following collection of results. This is not meant to be an exhaustive list of all such results needed to prove all parts of Theorem (ref), but any and all omitted terms follow by identical reasoning.
We next state, without proof, the following fact about the rates appearing in all these Lemmas, which follows from elementary inequalities.
The next Lemma proves Eqn.\ (ref), a crucial step in the proof of Theorem (ref)(a). Because this result only involves undersmoothing, we will omit the subscript $p$ as above.
Identifying the terms of the expansion is a matter of straightforward, if tedious, calculation. The first four cumulants of the Studentized statistics must be calculated (due to James-Mayne1962_Sankhya), which are functions of the first four moments. In what follows, we give a short summary. Note well that we always discard higher-order terms for brevity, and to save notation we will write $\stackrel{o}{=}$ to stand in for “equal up to $o((nh)^{-1} + (nh)^{-1/2}\eta + \eta^2)$”, and including $o(\rho^{1+2(p+1)})$ for $T_\mathtt{bc}$.
The computations will be aided by putting all three estimators into a common structure. In close parallel to the density case, let us define $\hat{m}_1 := \hat{m}$ and $\hat{m}_2 = \hat{m} - \hat{m}_m$, $\sigma_1^2 := \sigma_\mathtt{us}^2$, and $\sigma_2^2 := \sigma_\mathtt{rbc}^2$, so that subscripts $1$ and $2$ generically stand in for undersmoothing and bias correction, respectively. With this in mind, we write \[T_\mathtt{us} = T_{1,1}, \qquad T_\mathtt{bc} = T_{2,1}, \qquad \text{ and } \qquad T_\mathtt{rbc} = T_{2,2},\] again paralleling the density case, so that the first subscript refers to the numerator and the second to the denominator. In the same vein, with some abuse of notation, we will also use\footnote{Throughout Section (ref), we use only generic polynomial orders $p$ and $q$, and so this notation will not conflict with the local linear or local quadratic fits, which would also be denoted $r_1(u)$ and $r_2(u)$, respectively.} $r_1(u) = r_p(u)$, $r_2(u) = r_q(u)$, $K_1(u) = K(u)$, $K_2(u) = L(u)$, $h_1 = h$, and $h_2 = b$, as well as
For the purpose of computing the expansion terms (i.e. moments of the two sides agree up to the requisite order), recalling the Taylor series expansion above, we will use
where we define, for $v \in \{1,2\}$,
where the final line defines $\ell^2_\mathtt{us}(X_i, X_j, X_k)$ in the obvious way following $\ell^1_\mathtt{us}$. To concretize the notation, for undersmoothing we are defining
In a similar way,
and specifically for undersmoothing and bias correction, let \[B_{1,1} = s \frac{1}{nh} \sum_{i=1}^n \ell_1^0(X_i) [m(X_i) - r_p(X_i - x)'\beta_p]\] and
Note that $\eta_\mathtt{us} = \mathbb{E}[B_{1,1}]$ and $\eta_\mathtt{bc} = \mathbb{E}[B_{2,1}]$.
Straightforward moment calculations yield
and
Computing each term in turn, we have
The expansion now follows, formally, from the following steps. First, combining the above moments into cumulants. Second, these cumulants may be simplified using that \[\frac{\sigma_v^2}{\sigma_w^2} = 1 + \mathbbm{1}(w\!\neq\! v) \left( \rho^{1+ (p+1)} \Omega_{1,\mathtt{bc}} + \rho^{1 + 2(p+1)} \Omega_{2,\mathtt{bc}} \right) \] and that in all cases present products such as $\ell_w^0(X_i)^{k_1} \ell_v^0(X_i)^{k_2}$ and $\ell_w^1(X_i, X_j)^{k_1} \ell_v^1(X_i, X_j)^{k_2}$ may be replaced with $\ell_v^0(X_i)^{k_1 + k_2}$ and $\ell_v^1(X_i, X_j)^{k_1 + k_2}$, respectively, provided the arguments match. This is immediate for $v=w$, and for $v \neq w$, follows because $\rho \to 0$ is assumed. This is the analogous step to Eqn.\ (ref) in the density case. For any term of a cumulant with a rate of $(nh)^{-1}$, $(nh)^{-1/2}\eta_v$, $\eta_v^2$, or $\rho^{1+2(p+1)}$ (i.e., the extent of the expansion), these simplifications may be inserted as the remainder will be negligible. Third, with the cumulants in hand, the terms of the expansion are determined as described by e.g., Hall1992_book.
In this section we present the results of a simulation study addressing the finite-sample performance of the methods described in the main paper. As with the density estimator, we report empirical coverage probabilities and average interval length of nominal 95% confidence interval for different estimators of a regression functions $m(x)$ evaluated at values $x=\{-2/3,-1/3,0,1/3,2/3\}$. For each replication, the data is generated as i.i.d. draws, $i=1,2,...,n$, $n=500$ as follows:
Models 1 to 3 were used by Fan-Gijbels1996_book and Cattaneo-Farrell2013_JoE, while Models 4 to 6 are from Hall-Horowitz2013_AoS, with some originally studied by Berry-Carroll-Ruppert2002_JASA. The regression functions are plotted in Figure (ref) together with the evaluation points used.
We compute confidence intervals for $m(x)$ using five alternative approaches:
In all cases the Epanechnikov kernel is used. The bandwidth $h$ is chosen in three different ways:
For the construction of the variance estimators $\hat{\sigma}_\mathtt{us}^2$ and $\hat{\sigma}_\mathtt{rbc}^2$ we consider HC$3$ plug-in residuals when forming the $\Sigma$ matrix. In Table (ref) we report empirical coverage and average interval length of RBC 95% Confidence Intervals (only for Model 5) using $\hat{h}_\mathtt{mse}$ for different variance estimators. The results reflect the robustness of the findings to this choice.
The results are presented in detail in the tables and figures below to give a complete picture of the performance of robust bias correction. First, Tables (ref)-(ref) show, for each regression model, respectively, the performance of the five methods above, in terms of empirical coverage and interval length, for all evaluation points and bandwidth choices (recall that $I_\mathtt{us}$ and $I_\mathtt{bc}$ have the same length). Panel A of each shows the coverage and length, while Panel B gives summary statistics for the two fully data-driven bandwidths. Note that in some cases, the population MSE-optimal bandwidth is not defined or is not computable numerically; usually because the bias is too small or other values are too extreme.
The broad conclusion from these tables is that robust bias correction provides excellent coverage and that the data-driven bandwidths perform well and are numerically stable. In almost all cases robust bias correction provides correct coverage, whereas the other methods often, but not always, fail to do so. In cases where there is little to no bias all the methods give good coverage. This can be seen in results for Models 2 and 4, at $|x|=2/3$, far enough away from the “hump” in the center of each, where the true regression function is (nearly) linear. But despite the encouraging results away from the center, only robust bias correction yields good coverage closer to the center ($|x| = 1/3$), when there is more bias. Going further, considering $x=0$, the center of the sharp peak in these models, we see that even robust bias correction fails to provide accurate coverage for $\hat{h}_\mathtt{rot}$, although $\hat{h}_\mathtt{dpi}$ performs slightly better. At this point, for these models, the bias is too extreme even for robust bias correction to overcome. The results for the other models yield similar lessons.
It is somewhat more difficult to compare interval length using these tables. The comparison is invited for a fixed bandwidth, in which case, by construction, undersmoothing will have a shorter length. However, this ignores the fact that robust bias correction can accommodate a larger range of bandwidths, and in particular will optimally use a larger bandwidth. For example, robust bias correction has excellent coverage in many cases for $\hat{h}_\mathtt{rot}$, which is in this case a data-driven MSE-optimal choice (i.e. they coincide). This bandwidth is generally larger than $\hat{h}_\mathtt{dpi}$, and hence undersmoothing generally covers better with the latter. However, if you compare the length of $I_\mathtt{us}(\hat{h}_\mathtt{rot})$ to the length of $I_\mathtt{us}(\hat{h}_\mathtt{dpi})$, we see that robust bias correction compares favorably in terms of length.
Both to better make this point and to illustrate the robustness of $I_\mathtt{rbc}$ to tuning parameter selection, Figures (ref)--(ref) show empirical coverage and length for all six models, and all evaluation points, across a range of bandwidths. The dotted vertical line shows the population MSE-optimal bandwidth (whenever available) for reference. The coverage figures highlight the delicate balance required for undersmoothing to provide correct coverage, and the generally poor performance of traditional bias correction, but show that for a wide range of bandwidths robust bias correction provides correct coverage. Further, interval length is not unduly inflated for bandwidths that provide correct coverage. Again, by construction, undersmoothing will yield shorter intervals for a fixed bandwidth, and this is clear from Figures (ref)--(ref), but it is also clear that robust bias correction can use much larger bandwidths while still maintaining correct coverage.
To further illustrate this idea, in Tables (ref)--(ref) we compare average interval length of US and RBC 95% confidence intervals but at different bandwidths. First, in Table (ref) we compute average interval length at the largest bandwidth that provides close to correct coverage for each method separately. Note that in all cases these bandwidths are not feasible: these are ex-post findings. Next, in Table (ref) we evaluate the performance of US and RBC confidence intervals at certain alternative bandwidths likely to be chosen in practice. First, we evaluate the performance of US confidence intervals at $h=\lambda\hat{h}_\mathtt{mse}$ for $\lambda=\{0.5;0.7\}$. We then compare the performance with RBC confidence intervals computed using the optimal, fully data-driven choices $\hat{h}_\mathtt{rot}$ and $\hat{h}_\mathtt{dpi}$. Both tables reflect that, once we control for coverage, intervals lengths do not differ systematically between both approaches.
Figures (ref)-(ref) make this same point in a different way. For a range of bandwidths, as in the previous figures, we show the “average position” of $I_\mathtt{us}$ and $I_\mathtt{rbc}$, where the center of the bar is placed at the average bias and the length of each bar is the average interval length across the simulations. The bars are then color-coded by coverage (green bars having good coverage, fading to red showing undercoverage). These make visually clear that although undersmoothing provides shorter intervals in general, that this comes at the expense of coverage, while robust bias correction provides good coverage for a range of bandwidths, many of which are “large” enough to yield narrow intervals.
All our methods are implemented in {\sf R} and {\tt STATA} via the {\tt nprobust} package, available from \url{http://sites.google.com/site/nppackages/nprobust} (see also \url{http://cran.r-project.org/package=nprobust}). See Calonico-Cattaneo-Farrell2017_nprobust for a complete description.
\forloop{model}{1}{\value{model} < 7}{
}