The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
83,650 characters
improved inference for nonparametric regression -0.25cmand regression-discontinuity designs
\date{March 3, 2026}
\title{\textsc{improved inference for nonparametric regression\\\vspace*{-0.25cm}and regression-discontinuity designs}\thanks{We are grateful for comments from Matias Cattaneo, Bruce Hansen, Phillip Heiler, Guido Imbens, Michael Jansson, Sid Kankanala, Dennis Kristensen, Xinwei Ma, Mikkel S{\o}lvsten, Guo Yan, seminar participants at Aarhus~U, Bank of Portugal, Harvard-MIT, MacQuarie~U, Queen's~U, Rutgers~U, Tilburg~U, TU~Dortmund, U~Bologna, UC3M, U~Connecticut, U~Copenhagen, U~Exeter, UNSW, U~Sydney, and participants at several conferences and workshops.
Cavaliere and Zanelli acknowledge support from the Italian Ministry of University and Research (PRIN 2020 Grant 2020B2AKFW and PRIN 2022 Grant 2020B2AKFW). Nielsen and Zanelli are grateful for support from the Danish National Research Foundation (DNRF Chair grant DNRF154) and the Aarhus Center for Econometrics (ACE) funded by the Danish National Research Foundation grant DNRF186. Gon\c{c}alves acknowledges support from the Natural Science and Engineering Research Council of Canada (NSERC) grant RGPIN-2021-02663.}}
\author{Giuseppe Cavaliere\\University of Bologna, Italy\\\& Exeter Business School, UK
\and S\'{\i}lvia Gon\c{c}alves\\McGill University, Canada\\\& CIREQ, Canada
\and Morten {\O}rregaard Nielsen\\Aarhus Center for Econometrics\\Aarhus University, Denmark
\and Edoardo Zanelli\\Aarhus Center for Econometrics\\Aarhus University, Denmark}
\maketitle
\begin{abstract}
Nonparametric regression and regression-discontinuity designs suffer from smoothing bias that distorts conventional confidence intervals. Solutions based on robust bias correction (\textsf{RBC}) are now central to the economist's toolbox. In this paper, we establish a novel connection between \textsf{RBC} methods and bootstrap prepivoting. Revisiting RBC through the lens of bootstrapping allows us to develop a novel bias correction procedure which delivers improved nonparametric inference. The resulting confidence intervals are 17\% shorter than the usual intervals employed in curve estimation and regression discontinuity designs, without compromising asymptotic coverage. This holds regardless of evaluation point location, bandwidth choice, or regressor and error distribution.
\bigskip\noindent\textsc{Keywords}: Asymptotic bias; bootstrap; local polynomial estimation; nonparametric regression; regression-discontinuity; robust bias correction.
\bigskip\noindent\textsc{JEL codes}: C14; C21.
\end{abstract}
\newpage
\section{introduction}
Economists routinely rely on nonparametric regression and regression-discontinuity designs (RDD) to estimate causal effects. Inference in this setting is challenging due to the presence of a bias component affecting the nonparametric estimators and statistics usually employed by practitioners, even asymptotically. Several solutions to this inference problem have been proposed in the econometrics and statistics literatures. These include the so-called undersmoothing and bias correction approaches; see, inter alia, Hall (1992), Imbens and Lemieux (2008), Calonico, Cattaneo, and Titiunik (2014), Calonico, Cattaneo and Farrell (2018), and the references therein. A recent leading solution is the robust bias correction (\textsf{RBC}) approach of Calonico et al.\ (2014, 2018), which explicitly resolves this problem by debiasing the reference nonparametric estimator, while adjusting its standard error to account for the additional uncertainty introduced by the bias correction. Bootstrap and resampling methods, though widely used by economists as a bias correction tool, are rarely employed in the present setting due to the fact that they generally fail to correct for smoothing bias (e.g.,\ Hall,~1992).
The main objective of this paper is to propose a new robust procedure for inference in the context of nonparametric regression and~RDD. Our approach is based on a non-standard application of the bootstrap that implicitly corrects for the bias, thus challenging the above fact. Standard bootstrap methods fail in the present context because the bootstrap cannot correctly mimic the asymptotic bias of the nonparametric estimator, leading to confidence intervals with incorrect coverage and bootstrap p-values which are not uniformly distributed (see also Cavaliere and Georgiev,~2020). We show that prepivoting, originally proposed by Beran (1987) as a means to obtain refinements for asymptotically unbiased estimators, can provide a solution to this failure, as motivated by recent work on asymptotically biased estimators in Cavaliere et al.\ (2024). The idea is that if the asymptotic distribution of the boostrap p-value can be consistently estimated, then we can use it to restore validity by transforming (or prepivoting) a non-uniformly distributed bootstrap p-value into one that is uniformly distributed. In our context, we can use its quantiles to adjust the original bootstrap confidence intervals.
In this framework, our first contribution is to show that prepivoting performs an implicit bias correction and at the same time adjusts the standard error of the nonparametric regression estimator to account for the additional uncertainty introduced by debiasing. That is, the prepivoted bootstrap interval is asymptotically equivalent to an \textsf{RBC}-style interval. We also show that, for certain bootstrap schemes, bootstrap inference based on prepivoting can be computationally straightforward, as it does not require any resampling algorithm. This desirable property occurs because both the bootstrap bias correction and the standard error adjustment are known functions of the regressors (and the kernel). When the external random variable used to construct the bootstrap confidence interval is drawn from a normal distribution, the prepivoted bootstrap interval is in fact identical to an \textsf{RBC}-style interval.
As our second contribution, we show that the \textsf{RBC} interval of Calonico et al.\ (2014, 2018) is asymptotically equivalent to prepivoting a confidence interval based on the bootstrap scheme considered by Bartalotti et al.\ (2017) and He and Bartalotti (2020) for RDD. For inference based on a local linear estimator at a given evaluation point (e.g., the cutoff for RDD), this bootstrap method generates bootstrap observations using a local quadratic estimator estimated at that point. For this reason, we call this the `global polynomial' (GP) bootstrap. The GP bootstrap is asymptotically invalid when used to construct conventional bootstrap intervals. Bartalotti et al.\ (2017) and He and Bartalotti (2020) estimate the bias by a first level bootstrap, and then use a second level bootstrap to obtain standard errors. However, we show that prepivoting solves the invalidity problem directly by modifying the bootstrap quantiles to reflect the presence of bias. In particular, we show that the prepivoted GP bootstrap (hereafter, \textsf{PGP}) interval is asymptotically equivalent to the \textsf{RBC} interval of Calonico et al.\ (2014, 2018) based on local quadratic estimates of the second-order derivative that enters the bias component. More generally, by changing the estimator of the conditional mean function used in the bootstrap DGP, we obtain different bias corrections and different studentizations.
By exploiting the asymptotic equivalence between prepivoting and \textsf{RBC}-style intervals, our third and main contribution is to revisit the classical bootstrap in nonparametric statistics (e.g., Härdle and Bowman, 1988, Härdle and Marron, 1991, and Hall and Horowitz, 2013). This bootstrap smooths the conditional mean function at each regressor value using a local polynomial estimator, and we label it the `local polynomial' (LP) bootstrap. The inconsistency of naive, i.e., non-prepivoted, implementations of the LP bootstrap when used with `large' (including MSE-optimal) bandwidths is well documented in the statistics literature. Different solutions have been proposed, and they typically require choosing additional bandwidth or tuning parameters. We propose and analyze new prepivoted LP bootstrap schemes (hereafter, \textsf{PLP}), and we show that they (i) deliver novel confidence intervals of the \textsf{RBC} style, not explored in the existing literature, and, more importantly, (ii) lead to efficiency gains (i.e., shorter confidence intervals) with respect to the classic \textsf{RBC} interval. Specifically, we show that the bias correction mechanism implicitly generated by prepivoting these bootstrap methods is more efficient than the bias estimator employed in the \textsf{RBC} procedure, resulting in confidence intervals of shorter length and correct coverage, asymptotically. Importantly, \textsf{PLP} allows using the same bandwidth for point estimation and in the bootstrap data generating process~(DGP). As far as we are aware, ours is the first approach in the nonparametrics literature that delivers a valid and efficient implementation of the LP bootstrap without the need for additional tuning parameters.
As is usually the case for nonparametric estimators, the asymptotic distribution of \textsf{PLP} statistics depends on whether the evaluation point is interior or boundary. For interior points, prepivoting can be applied directly as discussed above. However, for boundary points, the implicit bias correction mimics the leading bias term of the original statistic only up to a multiplicative factor. The consequence is that the standard prepivoting approach studied by Cavaliere et al.\ (2024) requires modification in this context. Since the multiplicative term depends on the regressors (and the kernel) only, and is therefore known, our solution involves reweighting the bootstrap statistic before applying the prepivoting approach. We show that this `modified' \textsf{PLP} approach (hereafter, \textsf{mPLP}) adapts to the nature of the evaluation point and yields asymptotically valid inference across the entire support of the regressor, including boundary points.
Importantly, the new \textsf{mPLP} method not only has correct asymptotic coverage, but it also delivers shorter interval length compared to the \textsf{RBC} interval. Intuitively, \textsf{mPLP} is based on a convolution of the original observations, and such convolution introduces an additional layer of smoothing. We show that this additional smoothing induces a more efficient bias-correction in the sense that the overall variance of the debiased statistic is smaller. Crucially, this efficiency result is not tied to a specific DGP, as the asymptotic relative lengths of the confidence intervals based on \textsf{mPLP} and \textsf{RBC} depend only on the choice of kernel in the nonparametric regression and on whether the evaluation point is interior or boundary. Table~\ref{tab:lengths} gives the asymptotic relative lengths of the \textsf{mPLP} bootstrap and \textsf{RBC} intervals for five common kernels for the special case of the local linear estimator. For the popular Epanechnikov kernel, \textsf{mPLP} bootstrap intervals are 17\% shorter than \textsf{RBC} intervals for both interior and boundary evaluation points. These efficiency results are supported by our Monte Carlo simulation experiments.
\begin{table}[t]
\caption{Comparison of relative (asymptotic) interval lengths}
\label{tab:lengths}
\vskip -6pt
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lccccc}
\toprule
& Triangular & Uniform & Epanechnikov & Biweight & Triweight \\
\midrule
Interior point & 0.84 & 0.86 & 0.83 & 0.85 & 0.86 \\
Boundary point & 0.84 & 0.86 & 0.83 & 0.84 & 0.85 \\
\bottomrule
\end{tabular*}
\vskip 4pt
{\footnotesize \emph{Notes}: Asymptotic relative length of confidence intervals based on \textsf{RBC} and \textsf{mPLP} for five different kernels for the local linear estimator. Entries less than one imply that \textsf{mPLP} is shorter.}
\end{table}
The remainder of this paper is organized as follows. Section~\ref{section RBC and prepivoting} presents some general results on prepivoting in nonparametric regression and its relation to \textsf{RBC}-style intervals, as well as to the GP bootstrap scheme. In Section~\ref{Section main}, we discuss how improved inference can be based on the LP bootstrap. We show asymptotic validity of the corresponding prepivoted confidence interval and analyze its efficiency, in terms of length, compared with the \textsf{RBC} interval. In Section~\ref{Section boundary} we extend our approach to boundary problems and, in particular, to (sharp)~RDD. In Section~\ref{Section Monte Carlo} we assess the performance of our methods in finite samples via Monte Carlo simulation. Section~\ref{Section guidance} provides guidance on practical implementation for applied researchers. Finally, Section~\ref{Section conclusions} concludes. All notation, technical derivations, and proofs are included in the appendix and in the accompanying supplemental material. \textsf{R} packages that implement the procedures in this paper are available at \url{https://pppackages.github.io}.
\section{robust bias correction and prepivoting}
\label{section RBC and prepivoting}
\subsection{preliminaries}
\label{section preliminaries}
Given a random sample $\mathscr{D}_n := \{(y_i, x_i): i=1,\dots,n\}$, consider the problem of inference on the conditional expectation $g(x):=\operatorname{\mathbb{E}}[y_i|x_i=x]$ at a fixed point $x=\mathsf{x}$, using a local ($p$-th order, $p$ odd) polynomial estimator $\hat{g}_n(x): = \hat{g}_n(x;h,K)$, where $h=h_n>0$ is the bandwidth and $K$ the kernel function; that is, with $r_p(u)=(1,u,\ldots ,u^{p})^{\prime }$,
\begin{equation}
\label{eq LP estimator}
\hat{g}_n(x):=\iota _0^{\prime }\hat{\beta}_{p,n}(x),\quad \hat{\beta}_{p,n}(x):=\operatorname*{\arg\!\min}_{b\in \operatorname{\mathbb{R}} ^{p+1}}\sum_{i=1}^{n}
\big(y_i-r_p(x_i-x)^{\prime }b\big) ^2 K\big( \tfrac{x_i-x}{h}\big) .
\end{equation}
The estimator admits the linear closed-form expression $\hat{g}_n(x) = \sum_{i=1}^n w_i (x) y_i /(nh)$ with weights $w_i(x)$ depending on the data, the polynomial order $p$, and the kernel function~$K$; see Appendix~\ref{app:definition}. Inference based on $\hat{g}_n(\mathsf{x})$ is challenging due to the presence of an asymptotic bias. In particular, we typically have
\begin{equation}
\label{biasedCLT}
T_n:= \sqrt{nh} (\hat{g}_n(\mathsf{x}) - g(\mathsf{x})) \operatorname*{\overset{d}{\to}} N(B,v^2_1),
\end{equation}
where, letting $\mathscr{X}_n:=\{x_i : i=1,\dots,n\}$, the asymptotic bias $B:=\operatorname*{\mathrm{plim}}_{n \to \infty} B_n$, $B_n:=\operatorname{\mathbb{E}}[T_n|\mathscr{X}_n]$, is non-zero when bandwidths with MSE-optimal rate are employed. Sufficient conditions for \hyperref[{biasedCLT}]{\textup{\tagform@{\ref*{biasedCLT}}}} are given by Assumptions~\ref{Ass_g}--\ref{Ass_K,h} below; e.g., Ruppert and Wand (1994) and Fan and Gijbels (1996).
\begin{assumption}\label{Ass_g}
$(y_i,x_i)$ are i.i.d.\ and $x_i$ has bounded support $\mathbb{S}_x$ and continuous density~$f$. Let $\mathsf{N}$ be an open neighborhood of~$\mathsf{x}$. Then, on $\mathsf{N} \cap \mathbb{S}_x$: (i)~$f(x)>0$; (ii)~$\sigma^2 (x):=\operatorname{\mathbb{V}}[y_i|x_i=x]>0$ is continuous and $\sup_x \operatorname{\mathbb{E}}[y_i^4|x_i=x]<\infty$; (iii)~$g^{(p+1)}$ is H\"{o}lder continuous with exponent~$\eta>0$.
\end{assumption}
\begin{assumption}\label{Ass_K,h}
(i)~The kernel $K$ has support $(-1,1)$ on which it is symmetric, positive, and Lipschitz continuous; (ii)~the bandwidth $h=h_n$ is such that $h + (nh)^{-1} \log n \to 0$ and $nh^{2p+3}\to\kappa \in [0,+\infty)$.
\end{assumption}
Assumption~\ref{Ass_K,h} covers the most widely used kernel functions, e.g., the uniform, Epanechnikov, and triangular kernels, as well as `large' ($\kappa >0$) and undersmoothing ($\kappa = 0$) bandwidths. Under Assumptions~\ref{Ass_g}--\ref{Ass_K,h}, the bias term $B_n$ satisfies
\begin{equation} \label{eq:Bn}
B_n = \sqrt{nh^{2p+3}}g^{(p+1)}(\mathsf{x})C_n /(p+1)!+\operatorname*{o_{\mathrm{p}}} (1),
\end{equation}
where $C_n =C_n (\mathsf{x})$ is a random variable satisfying $C_n \operatorname*{\to_{\mathrm{p}}} C$ for a non-zero constant $C$ defined in Appendix~\ref{app:definition}. Note that for undersmoothing bandwidths, $B_n=\operatorname*{o_{\mathrm{p}}} (1)$ and $B=0$.
The presence of bias invalidates confidence intervals that ignore it. For instance, with $\alpha \in (0,1)$, a $100(1-\alpha)$\% nominal level conventional interval that ignores the bias is
\begin{equation*}
\mathrm{CI}_n := \big[ \hat{g}_n(\mathsf{x}) \pm z_{1-\alpha/2}(nh)^{-1/2}\hat{v}_{1n}\big],
\end{equation*}
where $\hat{v}^2_{1n}$ is a consistent estimator of~$v^2_1$. However, $\operatorname{\mathbb{P}} ( g(\mathsf{x}) \in \mathrm{CI}_n )$ does not converge to $1-\alpha$ as $n\to \infty$. To restore asymptotic validity of the confidence interval, the \textsf{RBC} approach of Calonico et al.\ (2014, 2018) recenters $\mathrm{CI}_n$ by an estimate of the leading bias term and adjusts the standard error to account for the variability of the bias estimator. With $\hat{B}_{\textsf{RBC},n}$ and $\hat{v}_{\textsf{RBC},n}$ denoting the bias estimate and the adjusted standard error, respectively, the $100(1-\alpha)$\% nominal level \textsf{RBC} interval is
\begin{equation}
\label{eq RBC int}
\mathrm{CI}_{\mathsf{RBC},n} := \big[ (\hat{g}_n(\mathsf{x}) - (nh)^{-1/2}\hat{B}_{\textsf{RBC},n}) \pm z_{1-\alpha/2}(nh)^{-1/2} \hat{v}_{\textsf{RBC},n} \big].
\end{equation}
For local polynomial estimators, $\hat{B}_{\textsf{RBC},n}$ is based on nonparametric estimators of higher-order derivatives, which may involve the choice of additional bandwidths. It has the form $\hat{B}_{\textsf{RBC},n} = \sqrt{nh^{2p+3}}\hat{g}_n^{(p+1)}(\mathsf{x})C_n/(p+1)!$, where $C_n$ is as introduced previously, and $\hat{g}_n^{(p+1)}(\mathsf{x})$ is an estimator of~$g^{(p+1)}(\mathsf{x})$. The estimator $\hat{v}^2_{\textsf{RBC},n}$ estimates the variance of the debiased statistic $\sqrt{nh}\hat{g}_n(\mathsf{x})-\hat{B}_{\textsf{RBC},n}$. As such, it is equal to $\hat{v}^2_{1n}$ (the variance estimator of $\sqrt{nh}\hat{g}_n(\mathsf{x})$) plus an additional term that estimates the variance of $\hat{B}_{\textsf{RBC},n}$ and the covariance between $\hat{B}_{\textsf{RBC},n}$ and $\sqrt{nh}\hat{g}_n(\mathsf{x})$; see Calonico et al.\ (2018, Section~3.1).
\subsection{an equivalence result}
\label{section equivalence}
As we show next, \textsf{RBC} intervals can be derived by prepivoting specific (invalid) bootstrap schemes. In Section~\ref{Section main} we exploit the link between the two approaches to derive new \textsf{RBC}-type intervals with improved length properties.
Let $\hat{g}^*_n(\mathsf{x}): = \hat{g}^*_n(\mathsf{x};h,K)$ denote a bootstrap analogue of $\hat{g}_n(\mathsf{x})$ based on the same kernel function $K$ and the same bandwidth~$h$. Conventional percentile bootstrap intervals are typically based on the quantiles of the bootstrap distribution of $T^*_n$, a given bootstrap analogue of $T_n$ (for instance, $T^*_n = \sqrt{nh} ( \hat{g}^*_n(\mathsf{x}) - \hat{g}_n(\mathsf{x}))$, see Section~\ref{Section main}). Letting $\hat{L}_n(u):= \operatorname{\mathbb{P}}^* ( T_n^* \leq u )$, a standard $100(1-\alpha)$\% nominal level equal-tailed bootstrap interval for $g(\mathsf{x})$ is
\begin{equation} \label{eq general invalid BS}
\mathrm{CI}_{\mathsf{Boot},n} := \big[ \hat{g}_n(\mathsf{x}) - (nh)^{-1/2} \hat{L}^{-1}_n(1-\alpha/2), \hat{g}_n (\mathsf{x}) - (nh)^{-1/2} \hat{L}^{-1}_n(\alpha/2) \big] ,
\end{equation}
which, however, does not contain $g(\mathsf{x})$ with probability converging to $1-\alpha$. Just as the original nonparametric estimator $\hat{g}_n(\mathsf{x})$ is biased, its bootstrap analogue $\hat{g}^*_n(\mathsf{x})$ is also biased, implying the invalidity of~$\mathrm{CI}_{\mathsf{Boot},n}$. As discussed by Cavaliere et al.\ (2024, Remark~3.2), the root of the problem lies in the fact that the standard bootstrap p-value $\hat{p}_n:=\operatorname{\mathbb{P}}^*(T^*_n\leq T_n)$ is not asymptotically distributed as $U_{[0,1]}$ (standard uniform) in the presence of asymptotic bias.
Prepivoting solves this problem by changing the level of the bootstrap quantiles. Let $\hat{H}_n$ be a (uniformly) consistent estimator of $H$, the asymptotic cumulative distribution function (cdf) of~$\hat{p}_n$. A $100(1-\alpha)$\% nominal level equal-tailed bootstrap interval based on prepivoting is
\begin{equation} \label{BTprepCI}
\mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n} := \big[ \hat{g}_n (\mathsf{x}) - (nh)^{-1/2} \hat{L}_n^{-1} (\hat{H}_n^{-1}(1-\alpha/2)), \hat{g}_n (\mathsf{x}) - (nh)^{-1/2} \hat{L}_n^{-1} (\hat{H}_n^{-1}(\alpha/2 ))\big].
\end{equation}
When $\hat{p}_n\operatorname*{\to_{\mathrm{d}}} U_{[0,1]}$ then $H(u)=u$ and $\hat{H}_n^{-1}(\alpha )\operatorname*{\to_{\mathrm{p}}}\alpha$, in which case the prepivoted interval $\mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n}$ is asymptotically equivalent to $\mathrm{CI}_{\mathsf{Boot},n}$. But when asymptotic uniformity of $\hat{p}_n$ fails, $\mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n}$ remains asymptotically valid while $\mathrm{CI}_{\mathsf{Boot},n}$ becomes invalid. The main reason why prepivoting restores validity is that if $\hat{H}_n$ is a uniformly consistent estimator of $H$ (regardless of whether the latter is the uniform cdf or not), the prepivoted p-value $\tilde{p}_n:=\hat{H}_n(\hat{p}_n)$ is asymptotically distributed as $U_{[0,1]}$ by the probability integral transform. This implies asymptotic validity of $\mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n}$, as shown by Cavaliere et al.\ (2024, Remark~3.4); that is, for any $\alpha \in (0,1)$,
\begin{align}
\operatorname{\mathbb{P}} ( g(\mathsf{x}) \in \mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n} ) &= \operatorname{\mathbb{P}} ( \hat{L}_n^{-1} ( \hat{H}_n^{-1}(\alpha/2 ))\leq T_n \leq \hat{L}_n^{-1} ( \hat{H}_n^{-1}(1 - \alpha/2 )) ) \nonumber \\
&= \operatorname{\mathbb{P}} ( \hat{H}_n^{-1}(\alpha/2 )\leq \hat{p}_n \leq \hat{H}_n^{-1}(1 - \alpha/2 ) )
\to 1- \alpha.
\label{CIp coverage}
\end{align}
Prepivoting avoids an explicit bias correction whenever $\hat{H}_n$ does not require estimating~$B$. Denote by $\hat{B}_n :=\operatorname{\mathbb{E}}^*[T_n^*]$ and $\hat{v}^2_n:=\operatorname{\mathbb{V}}^*[T_n^*]$ the bias and variance, respectively, of the bootstrap statistic $T^*_n$ under the bootstrap probability measure. As we show next, we may regard $\hat{B}_n$ as a bias estimator implicitly induced by the bootstrap. Suppose the following high level conditions hold.
\setcounter{assnLetter}{7}
\begin{assnLetter}
\label{assn 1}
(i) $T^*_n-\hat{B}_n \operatorname*{\overset{d*}{\to_{\mathrm{p}}}} N(0,v^2)$ with $v^2=\operatorname*{\mathrm{plim}} \hat v^2_n>0$; (ii) $(T_n - B_n, \hat{B}_n - B_n)^\prime \operatorname*{\overset{d}{\to}} N(0,V)$, where $V$ is positive definite.
\end{assnLetter}
Assumption~\ref{assn 1}(i) requires the bootstrap statistic $T^*_n$ to be conditionally asymptotically Gaussian after we remove the bootstrap bias~$\hat{B}_n$. Note that the bootstrap does not need to replicate the asymptotic variance of~$T_n$, i.e., $v^2$ may be different from~$v^2_1$. Under Assumption~\ref{assn 1}(i), as in Theorem~3.4 of Cavaliere et al.\ (2024), we can show that $\hat{p}_n=\Phi(\hat{v}_n^{-1}(T_n-\hat{B}_n))+\operatorname*{o_{\mathrm{p}}}(1)$. If in addition Assumption~\ref{assn 1}(ii) holds, then $T_n-\hat{B}_n\operatorname*{\overset{d}{\to}} N(0,v^2_\mathsf{P})$ with $v_{\mathsf{P}}^2>0$, and the asymptotic cdf of $\hat p_n$ is $H (u):= \Phi (m^{-1} \Phi^{-1} (u))$, $m:= {v}_{\mathsf{P}}/{v}$. A consistent (plug-in) estimator of $H$ can be obtained as
\begin{equation}\label{Hhat_def}
\hat{H}_n(u):=\Phi(\hat{m}_n^{-1}\Phi^{-1}(u)), \quad \hat{m}_n:=\hat{v}_{\mathsf{P},n}/\hat{v}_n,
\end{equation}
where $\hat{v}^2_{\mathsf{P},n}$ is a consistent estimator of~$v^2_{\mathsf{P}}$. Different bootstrap methods induce different $\hat{B}_n$ and consequently require different estimators~$\hat{v}^2_{\mathsf{P},n}$.
Although $\hat{H}_n$ of \hyperref[{Hhat_def}]{\textup{\tagform@{\ref*{Hhat_def}}}}, and hence $\mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n}$, does not involve $\hat{B}_n$ explicitly, we can show that $\mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n}$ is asymptotically equivalent to an interval that bears a close resemblance to the \textsf{RBC} interval $\mathrm{CI}_{\mathsf{RBC},n}$ in that it too contains a bias correction and adjusts the standard error for the additional uncertainty appropriately. We describe this result next.
For simplicity, suppose that $T^*_n$ is conditionally Gaussian; that is, conditionally on $\mathscr{D}_n$, $T^*_n \sim N(\hat{B}_n,\hat{v}^2_n)$ and hence $\hat{L}_n(u) = \Phi ((u-\hat{B}_n )/\hat{v}_n )$. In this case, $\mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n}$ in \hyperref[{BTprepCI}]{\textup{\tagform@{\ref*{BTprepCI}}}} is \emph{exactly} equal to
\begin{equation} \label{CIprbc}
\mathrm{CI}_{\mathsf{P} \textsf{-} \mathsf{RBC},n} =\big[ (\hat{g}_n(\mathsf{x}) - (nh)^{-1/2}\hat{B}_n) \pm z_{1-\alpha/2}(nh)^{-1/2} \hat{v}_{\mathsf{P},n}\big].
\end{equation}
In general, $T^*_n$ is Gaussian only asymptotically, and we obtain the following result.
\begin{theorem} \label{th 1}
Let $\mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n}$, $\hat{H}_n$, and $\mathrm{CI}_{\mathsf{P} \textsf{-} \mathsf{RBC},n}$ be defined in \hyperref[{BTprepCI}]{\textup{\tagform@{\ref*{BTprepCI}}}}, \hyperref[{Hhat_def}]{\textup{\tagform@{\ref*{Hhat_def}}}}, and \hyperref[{CIprbc}]{\textup{\tagform@{\ref*{CIprbc}}}}, respectively.
\begin{itemize}
\item[(i)] If Assumption~\ref{assn 1}(i) holds then $\mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n}=\mathrm{CI}_{\mathsf{P} \textsf{-} \mathsf{RBC},n}+\operatorname*{o_{\mathrm{p}}}((nh)^{-1/2})$.
\item[(ii)] If Assumption~\ref{assn 1}(i) holds and $T_n^*$ is Gaussian conditional on~$\mathscr{D}_n$, then $\mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n}=\mathrm{CI}_{\mathsf{P} \textsf{-} \mathsf{RBC},n}$~a.s.
\item[(iii)] If Assumption~\ref{assn 1} holds and $\hat{v}^2_{\mathsf{P},n}\operatorname*{\to_{\mathrm{p}}} v^2_{\mathsf{P}}$, then $ \operatorname{\mathbb{P}} ( g(\mathsf{x}) \in \mathrm{CI}_{\mathsf{P}\hspace{-0.5pt},n} ) \to 1- \alpha$.
\end{itemize}
\end{theorem}
An interesting implication of Theorem~\ref{th 1} is that we can obtain $\mathrm{CI}_{\mathsf{P} \textsf{-} \mathsf{RBC},n}$ without any resampling as it depends only on $\hat{B}_n$ and $\hat{v}^2_{\mathsf{P},n}$ (which we derive analytically for different residual-based bootstrap methods in the next sections), and standard normal critical values~$z_\alpha$. More importantly, this asymptotic equivalence result implies that we can view \textsf{RBC} through the lens of prepivoting (or vice versa), as we show in the next section.
\subsection{standard rbc inference as bootstrap prepivoting}
\label{section GP}
Assume that the bootstrap data $\mathscr{D}_n^*:=\{(y_i^*,x_i^*): i=1,\dots, n\}$ are generated by setting $x_i^* = x_i$ for $i=1,\dots, n$ (`fixed regressor' bootstrap) and
\begin{equation}
y_i^{\ast}=\tilde{g}_{\mathsf{x},n}(x_i)+\varepsilon _i^{\ast}\text{,\quad }
i=1,\ldots ,n , \label{btsDGPLQ}
\end{equation}
where $\varepsilon _i^{\ast}:=\tilde\varepsilon_i e_i^{\ast}$ with $\tilde\varepsilon_i$ being the bias-corrected HC\emph{k} residuals considered in Calonico et al.\ (2018) and $e_i^{\ast}$ being (conditionally on the data)\ i.i.d.\ with $\operatorname{\mathbb{E}}^{\ast}[e_i^{\ast}]=0$ and $\operatorname{\mathbb{E}}^{\ast}[e_i^{\ast 2}]=1$. The bootstrap conditional expectation $\operatorname{\mathbb{E}}^{\ast}[y_i^{\ast}|x^*_i=x_i]$ in \hyperref[{btsDGPLQ}]{\textup{\tagform@{\ref*{btsDGPLQ}}}} is
\begin{equation*}
\tilde{g}_{\mathsf{x},n} (x_i):= r_{p+1}(x_i - \mathsf{x})'\hat{\beta}_{p+1,n}(\mathsf{x})=\sum_{j=0}^{p+1}\iota _{j}^{\prime }\hat{\beta}_{p+1,n}(\mathsf{x})(x_i-\mathsf{x})^{j},
\end{equation*}
where $\hat{\beta}_{p+1,n}(\mathsf{x})$ is obtained as in \hyperref[{eq LP estimator}]{\textup{\tagform@{\ref*{eq LP estimator}}}} using local polynomial regression of order $p+1$ at the \emph{fixed} point~$\mathsf{x}$; see Bartalotti et al.\ (2017) and He and Bartalotti (2020) for the case $p=1$. Since $\iota _{j}^{\prime }\hat{\beta}_{p+1,n}(\mathsf{x})$ is an estimator of $g^{(j)}(\mathsf{x})/j!$, the function $\tilde{g}_{\mathsf{x},n} (x)$ is an estimator of a ($p+1$)-th order Taylor expansion of the conditional expectation $g(x)=\operatorname{\mathbb{E}}[y_i|x_i=x]$ around the fixed point~$\mathsf{x}$. That is, to generate the bootstrap data we use \emph{the same} function $\tilde{g}_{\mathsf{x},n} (x)$ estimated locally at $\mathsf{x}$ but evaluated globally at each $x=x_i$ for $i =1,\dots, n$. Hence, we use `global polynomial' (GP) bootstrap to label the DGP~\hyperref[{btsDGPLQ}]{\textup{\tagform@{\ref*{btsDGPLQ}}}}. Importantly, this method generates bootstrap data using a polynomial order ($p+1$) that is larger than the order $p$ used in estimation. This mismatch is key to inducing an asymptotic bias under the bootstrap probability measure that satisfies Assumption~\ref{assn 1}.
Let $\hat{g}_n^{\ast}(x)$ be the local ($p$-th order) polynomial estimator applied to the bootstrap sample generated as \hyperref[{btsDGPLQ}]{\textup{\tagform@{\ref*{btsDGPLQ}}}}. The GP bootstrap analogue of $T_n$ is $T_{\mathsf{GP},n}^{\ast}=\sqrt{nh}(\hat{g}_n^{\ast}(\mathsf{x})-\tilde{g}_{\mathsf{x},n} (\mathsf{x}))$, and its bootstrap bias and variance are denoted $\hat{B}_{\mathsf{GP},n}:=\operatorname{\mathbb{E}}^*[T_{\mathsf{GP},n}^\ast]$ and $\hat{v}_{\mathsf{GP},n}^2:=\operatorname{\mathbb{V}}^*[T_{\mathsf{GP},n}^\ast]$, respectively. Because the GP bootstrap is a fixed-design scheme (i.e., the regressors are fixed across bootstrap samples), these moments can be computed in closed form. Specifically, we show that
\begin{equation}
\hat{B}_{\mathsf{GP},n} =\sqrt{nh^{2p+3}}\hat{g}_n^{(p+1)}(\mathsf{x})\tfrac{C_n}{(p+1)!},
\label{BhatGP}
\end{equation}
where $\hat{g}_n^{(p+1)}(\mathsf{x}):=(p+1)!\iota_{p+1}' \hat{\beta}_{p+1,n}(\mathsf{x})$ is an estimator of~$g^{(p+1)}(\mathsf{x})$, clearly implying that $\hat{B}_{\mathsf{GP},n}$ is an estimator of the dominant term of $B_n$ in~\hyperref[{eq:Bn}]{\textup{\tagform@{\ref*{eq:Bn}}}}. Interestingly, $\hat{B}_{\mathsf{GP},n}$ is identical to $\hat{B}_{\mathsf{RBC},n}$ as given in Calonico et al.\ (2018, Section~3); see the discussion after~\hyperref[{eq RBC int}]{\textup{\tagform@{\ref*{eq RBC int}}}}.
As explained in the previous section, a prepivoted bootstrap confidence interval \hyperref[{BTprepCI}]{\textup{\tagform@{\ref*{BTprepCI}}}} based on $(T_n,T_{\mathsf{GP},n}^{\ast})$, say $\mathrm{CI}_{\mathsf{PGP}\hspace{-0.5pt},n}$, can be constructed by choosing $\hat{v}^2_{\mathsf{PGP},n}$ as a consistent estimator of $v^2_{\mathsf{PGP}}$, which is the asymptotic variance of $T_n-\hat{B}_{\mathsf{GP},n}$. Because $\hat{B}_{\mathsf{GP},n}=\hat{B}_{\mathsf{RBC},n}$, the natural plug-in estimator for the GP case coincides with the variance estimator in Calonico et al.\ (2018). In Theorem~\ref{Th RBC PGP equivalence} below we show (i)~validity of such confidence interval, (ii)~its asymptotic equivalence to the \textsf{RBC} interval $\mathrm{CI}_{\mathsf{RBC},n}$ in~\hyperref[{eq RBC int}]{\textup{\tagform@{\ref*{eq RBC int}}}}, and (iii)~its almost sure equality to the \textsf{RBC} interval for Gaussian bootstraps.
The proof requires verifying that Assumption~\ref{assn 1} holds for $(T_n,T_{\mathsf{GP},n}^{\ast})$ under Assumptions~\ref{Ass_g} and~\ref{Ass_K,h}. The theorem then follows by Theorem~\ref{th 1}.
\begin{theorem} \label{Th RBC PGP equivalence}
Under Assumptions~\ref{Ass_g} and~\ref{Ass_K,h}, for any $\mathsf{x}\in \mathbb{S}_{x}$,
$\operatorname{\mathbb{P}}(g(\mathsf{x})\in \mathrm{CI}_{\mathsf{PGP}\hspace{-0.5pt},n} )\to 1-\alpha$ and $\mathrm{CI}_{\mathsf{PGP}\hspace{-0.5pt},n} = \mathrm{CI}_{\mathsf{RBC},n} + \operatorname*{o_{\mathrm{p}}} ((nh)^{-1/2})$. If, in addition, $T_{\mathsf{GP},n}^{\ast}$ is Gaussian conditional on~$\mathscr{D}_n$, then $\mathrm{CI}_{\mathsf{PGP}\hspace{-0.5pt},n}=\mathrm{CI}_{\mathsf{RBC},n}$~a.s.
\end{theorem}
\section{improved nonparametric inference}
\label{Section main}
The main feature of the GP bootstrap discussed above is that the conditional mean function under the bootstrap probability measure is based on a \emph{single} local approximation of the curve~$g$, using a polynomial expansion of $g$ around~$\mathsf{x}$, which is then applied globally to all~$x_i$'s. In contrast, in this section we consider the `local polynomial' (LP) bootstrap scheme traditionally employed in the statistics literature, where $g$ is approximated locally for each of the $x_i$'s, implying that the bootstrap conditional expectation function resembles $g$ more closely over the entire support (Härdle and Marron, 1988, Härdle and Bowman, 1991, Hall and Horowitz, 2013). This point is illustrated in Figure~\ref{figure curves}, where we contrast the true regression curve with the GP and LP bootstrap curves under the setup considered in Section~\ref{Section Monte Carlo}.
\begin{figure}[t]
\includegraphics[width=\textwidth]{FLLPplot.pdf}
\caption{Nonparametric curves:\ original (\rule[0.5ex]{1em}{0.8pt}), LP bootstrap (\rule[0.5ex]{0.3em}{0.8pt}\hspace{0.3em}\rule[0.5ex]{0.3em}{0.8pt}\hspace{0.3em}\rule[0.5ex]{0.3em}{0.8pt}), and GP bootstrap (\rule[0.5ex]{0.3em}{0.8pt}\hspace{0.3em}\rule[0.5ex]{0.1em}{0.8pt}\hspace{0.3em}\rule[0.5ex]{0.3em}{0.8pt}\hspace{0.3em}\rule[0.5ex]{0.1em}{0.8pt})}
\vskip 2pt
{\footnotesize \emph{Notes}: The vertical bar corresponds to the evaluation point $\mathsf{x}=-1/3$; the shaded area corresponds to the selected bandwidth; $n=1,000$ observations are summarized into $n/5$ bins (circles); polynomial order $p=1$.}
\label{figure curves}
\end{figure}
In this section, we show how the prepivoted LP bootstrap leads to improved inference. We focus here on the case where $\mathsf{x}$ is an interior point. Boundary points will be considered in Section~\ref{Section boundary}.
\subsection{the local polynomial bootstrap}
\label{Section main-LP}
The LP bootstrap sets the bootstrap conditional expectation to $\operatorname{\mathbb{E}}^{\ast}[y_i^{\ast}|x^*_i=x] = \hat{g}_n(x)$, where $\hat{g}_n(x)$ is as in \hyperref[{eq LP estimator}]{\textup{\tagform@{\ref*{eq LP estimator}}}}. The bootstrap data $\mathscr{D}_n^*:=\{(y_i^*,x_i^*): i=1,\dots, n\}$ satisfy $x_i^* = x_i$, $i=1,\dots, n$, and
\begin{equation}
y_i^{\ast}=\hat{g}_n(x_i)+\varepsilon _i^{\ast}\text{,\quad }i=1,\ldots ,n, \label{btsDGPLP}
\end{equation}
where $\varepsilon_i^* := \tilde\varepsilon_i e_i^*$ with $\tilde\varepsilon_i$ and $e_i^*$ as defined in Section~\ref{section GP} to facilitate comparison with the existing \textsf{RBC} approach. Other choices are possible, e.g., $\tilde\varepsilon_i = y_i - \hat{g}_n(x_i)$, without changing our results. Bootstrap DGPs of this or similar forms have been widely applied to nonparametric regression without prepivoting; see, e.g., Härdle and Marron (1988), Härdle and Bowman (1991), and Hall and Horowitz (2013).
The LP bootstrap uses a $p$-th order local polynomial regression estimate of $g(x)$ calculated at each sample point $x=x_i$ to generate $y^*_i$, rather than extrapolating values for $g(x_i)$ from a single $(p+1)$-th order local polynomial regression estimated at $x=\mathsf{x}$, as for the GP bootstrap. In this section, we exploit this key difference and show the following. First, we demonstrate that prepivoting can also be used to restore asymptotic validity of the LP bootstrap without requiring undersmoothing or the choice of additional bandwidths, a result that stands in contrast to the existing bootstrap literature. Second, we show that the resulting prepivoted LP bootstrap (\textsf{PLP}) generates a novel RBC-type confidence interval which differs from the existing \textsf{RBC} intervals in that it does not rely on a higher-order local polynomial estimate of the derivatives entering the bias term. Third, and most importantly, we show that in large samples, the new \textsf{PLP} interval is shorter than the classical \textsf{RBC} CIs in Calonico et al.\ (2014, 2018), while maintaining correct coverage. Put differently, with respect to the usual bias estimators used in \textsf{RBC} procedures, the bias estimator implicitly generated by the \textsf{PLP} method leads to improved inference on the regression curve.
\subsection{theory}
\label{Section main-theory}
Consider the prepivoted bootstrap confidence interval \hyperref[{BTprepCI}]{\textup{\tagform@{\ref*{BTprepCI}}}} based on $(T_n,T_{\mathsf{LP},n}^{\ast})$, say $\mathrm{CI}_{\mathsf{PLP}\hspace{-0.5pt},n}$, where $T_{\mathsf{LP},n}^{\ast}=\sqrt{nh}(\hat{g}_n^{\ast}(\mathsf{x})-\hat{g}_n(\mathsf{x}))$ and $\hat{g}_n^{\ast}(\mathsf{x})$ is the local ($p$-th order) polynomial estimator applied to the bootstrap sample generated as~\hyperref[{btsDGPLP}]{\textup{\tagform@{\ref*{btsDGPLP}}}}. The bootstrap bias and variance of $T_{\mathsf{LP},n}^{\ast}$ are $\hat{B}_{\mathsf{LP},n}:=\operatorname{\mathbb{E}}^*[T_{\mathsf{LP},n}^\ast]$ and $\hat{v}_{\mathsf{LP},n}^2:=\operatorname{\mathbb{V}}^*[T_{\mathsf{LP},n}^\ast]$, respectively.
Recall from Section~\ref{section equivalence} that in order to construct $\mathrm{CI}_{\mathsf{PLP}\hspace{-0.5pt},n}$, a consistent estimator $\hat{v}_{\mathsf{PLP},n}^2$ of the asymptotic variance $v^2_{\mathsf{PLP}}$ of $T_n - \hat{B}_{\mathsf{LP},n}$ is required.
To find $\hat{v}_{\mathsf{PLP},n}^2$, we first derive the bootstrap bias~$\hat{B}_{\mathsf{LP},n}$. It follows from the closed-form expression for $\hat{g}^{\ast}_n(\mathsf{x})$ that
\begin{align} \label{eq:BhatLP0}
\hat{B}_{\mathsf{LP},n}
= \sqrt{nh} \Big( \frac{1}{nh} \sum_{i=1}^n w_i ( \mathsf{x})\hat{g}_n (x_i) - \hat{g}_n ( \mathsf{x}) \Big).
\end{align}
Since $B_n = \sqrt{nh}\big((nh)^{-1}\sum_{i=1}^n w_i(\mathsf{x}) g(x_i) - g(\mathsf{x})\big)$, \hyperref[{eq:BhatLP0}]{\textup{\tagform@{\ref*{eq:BhatLP0}}}} makes clear that the LP bootstrap bias $\hat{B}_{\mathsf{LP},n}$ differs from $\hat{B}_{\mathsf{RBC},n}$ by targeting the bias directly, rather than estimating the $(p+1)$-th derivative that appears in its asymptotic expansion. We can rewrite $\hat{B}_{\mathsf{LP},n}$ as
\begin{equation}\label{eq Bhat LP}
\hat{B}_{\mathsf{LP},n} = \frac{1}{\sqrt{nh}} \sum_{i=1}^n w_{\mathsf{LP} \text{-}\mathsf{bc} ,i} (\mathsf{x})y_i
\end{equation}
with weights defined by the convolution $w_{\mathsf{LP}\text{-}\mathsf{bc},i} (\mathsf{x}):= (nh)^{-1} \sum_{j=1}^n w_{j} ( \mathsf{x}) w_i ( x_j)- w_i( \mathsf{x})$. Since $T_n = \sqrt{nh}\big( (nh)^{-1} \sum_{i=1}^n w_i(\mathsf{x})y_i - g(\mathsf{x})\big)$, we obtain
\begin{equation*}
T_n - \hat{B}_{\mathsf{LP},n} = \sqrt{nh} \Big( \frac{1}{{nh}} \sum_{i=1}^n w_{\mathsf{PLP},i}(\mathsf{x}) y_i - g(\mathsf{x})\Big), \hspace{0.5em} w_{\mathsf{PLP},i} (\mathsf{x}) := w_i(\mathsf{x}) - w_{\mathsf{LP}\text{-}\mathsf{bc},i} (\mathsf{x}) .
\end{equation*}
Let $v^2_{\mathsf{PLP}}$ denote the asymptotic variance of~$T_n - \hat{B}_{\mathsf{LP},n}$. Because the weights are measurable with respect to $\mathscr{X}_n$ and $\operatorname{\mathbb{V}} [y_i | \mathscr{X}_n] = \sigma^2 (x_i)$, by conditioning on the regressors we have $\operatorname{\mathbb{V}} [T_n - \hat{B}_{\mathsf{LP},n}| \mathscr{X}_n] = (nh)^{-1} \sum_{i=1}^n w_{\mathsf{PLP},i}(\mathsf{x})^2 \sigma^2 (x_i)$. This suggests the estimator
\begin{equation} \label{vPest}
\hat{v}_{\mathsf{PLP},n}^2:= \frac{1}{nh} \sum_{i=1}^n w_{\mathsf{PLP},i}(\mathsf{x})^2 \tilde\varepsilon_i^2 .
\end{equation}
A formal proof of consistency of $\hat{v}_{\mathsf{PLP},n}^2$ under Assumptions~\ref{Ass_g} and~\ref{Ass_K,h} is provided next.
\begin{lemma} \label{Lemma v2plp}
Under Assumptions~\ref{Ass_g} and~\ref{Ass_K,h}, $\hat{v}_{\mathsf{PLP},n}^2 \operatorname*{\to_{\mathrm{p}}} v_{\mathsf{PLP}}^2 >0$.
\end{lemma}
In view of Lemma~\ref{Lemma v2plp}, an application of Theorem~\ref{th 1} implies the asymptotic validity of the \textsf{PLP} interval~$\mathrm{CI}_{\mathsf{PLP}\hspace{-0.5pt},n}$. In addition, $\mathrm{CI}_{\mathsf{PLP}\hspace{-0.5pt},n}$ is asymptotically equivalent to
\begin{equation*}
\mathrm{CI}_{\mathsf{PLP} \textsf{-} \mathsf{RBC},n} = [(\hat{g}_n(\mathsf{x})-(nh)^{-1/2} \hat{B}_{\mathsf{LP},n})\pm z_{1-\alpha /2}(nh)^{-1/2}\hat{v}_{\mathsf{PLP},n}],
\end{equation*}
which is a new $\mathsf{RBC}$-type interval that does not require resampling. The following theorem states these results.
\begin{theorem} \label{Th RBC PLP equivalence}
Let $\mathsf{x}$ be an interior point. Under Assumptions~\ref{Ass_g} and~\ref{Ass_K,h}, $\operatorname{\mathbb{P}}(g(\mathsf{x})\in \mathrm{CI}_{\mathsf{PLP}\hspace{-0.5pt},n})\to 1-\alpha$ and $\mathrm{CI}_{\mathsf{PLP}\hspace{-0.5pt},n}=\mathrm{CI}_{\mathsf{PLP} \textsf{-} \mathsf{RBC},n} + \operatorname*{o_{\mathrm{p}}} ((nh)^{-1/2})$. If, in addition, $T_{\mathsf{LP},n}^{\ast}$ is Gaussian conditional on~$\mathscr{D}_n$, then $\mathrm{CI}_{\mathsf{PLP}\hspace{-0.5pt},n}=\mathrm{CI}_{\mathsf{PLP} \textsf{-} \mathsf{RBC},n}$~a.s.
\end{theorem}
\subsection{improved inference}
\label{Section main-improved}
In this section, we compare the asymptotic efficiency of the new prepivoted \textsf{PLP} intervals and the existing \textsf{RBC}-type intervals proposed by Calonico et al.\ (2014, 2018). We show that the former are asymptotically shorter than the latter. This efficiency result is based on the asymptotic equivalence between prepivoting and the \textsf{RBC} approach in Section~\ref{section equivalence}. In particular, Theorems~\ref{Th RBC PGP equivalence} and~\ref{Th RBC PLP equivalence} imply that the asymptotic relative length of $\mathrm{CI}_{\mathsf{PLP}\hspace{-0.5pt},n}$ compared to $\mathrm{CI}_{\mathsf{RBC},n}$ is equal to the ratio between the asymptotic standard deviations of their respective debiased statistics, i.e., $v_{\mathsf{PLP}}/v_{\mathsf{RBC}}$.
To shed some light on this ratio, let $\mathsf{w}(u), u\in \operatorname{\mathbb{R}}$, be the so-called `equivalent kernel' associated with $K$ and~$p$. Specifically, for $\mathcal{X}=\mathbb{R}$, $\mathsf{w}(u) := \iota_0^{\prime}(\int_{\mathcal{X}} r_p(s) r_p^{\prime}(s)K(s) ds)^{-1} r_p(u) K(u)$, the asymptotic analogue of the sample weights~$w_i(x)$; see Appendix~\ref{app:definition}. Similarly, define the asymptotic analogue of $w_{\mathsf{PLP},i}(x)$ as $\mathsf{w}_{\mathsf{PLP}}(u) := 2 \mathsf{w}(u) - \int_{\mathcal{X}} \mathsf{w}(u-r) \mathsf{w}(r) dr$. This contrasts with that of the \textsf{RBC} method given by $\mathsf{w}_{\mathsf{RBC}}(u) := \mathsf{w}(u) - C\iota_{p+1}^{\prime}(\int_{\mathcal{X}} r_{p+1}(s) r_{p+1}^{\prime}(s)K(s) ds)^{-1} r_{p+1}(u) K(u)$; see Lemma~\ref{lemmaKernelConstants} and Calonico et al.\ (2022, Section~4.2).
\begin{corollary} \label{corollary v^2_P comparisons}
Under Assumptions~\ref{Ass_g} and~\ref{Ass_K,h},
\begin{equation*}
(v^2_{\mathsf{PLP}}, v^2_{\mathsf{RBC}})= \frac{\sigma^2 (\mathsf{x})}{f(\mathsf{x})} (\mathcal{K}_{\mathsf{PLP}}, \mathcal{K}_{\mathsf{RBC}}) ,
\end{equation*}
where
\begin{equation}
\label{eq K PLP}
\mathcal{K}_{\mathsf{PLP}} := \int_{\mathcal{X}} \mathsf{w}_{\mathsf{PLP}}(u)^2 du >0
\quad\text{and}\quad
\mathcal{K}_{\mathsf{RBC}} := \int_{\mathcal{X}} \mathsf{w}_{\mathsf{RBC}}(u)^2 du >0
\end{equation}
are functions only of the kernel $K$ and the polynomial order~$p$.
\end{corollary}
As shown in Section~\ref{Section main-theory}, $v^2_{\mathsf{PLP}}$ is the asymptotic variance of a (scaled) weighted average of~$y_i$ with weights~$w_{\mathsf{PLP},i}(\mathsf{x})$. Using fairly standard arguments in the nonparametric literature (e.g., Fan and Gijbels, 1996, p.~66), Corollary~\ref{corollary v^2_P comparisons} shows that this variance is proportional to $\sigma^2 (\mathsf{x})/f(\mathsf{x})$ where the constant of proportionality is~$\mathcal{K}_{\mathsf{PLP}}$, which is a known function of the kernel and the polynomial order,~$p$. A similar proportionality holds for $v^2_{\mathsf{RBC}}$ with the difference that the constant of proportionality is~$\mathcal{K}_{\mathsf{RBC}}$, which is a different function of the kernel and~$p$. Hence, the relative asymptotic length of the $\mathsf{PLP}$ and $\mathsf{RBC}$ intervals is the square root of $\mathcal{K}_{\mathsf{PLP}}/\mathcal{K}_{\mathsf{RBC}}$.
\begin{table}[t]
\caption{Relative (asymptotic) interval length comparison}
\label{table relative variances}
\vskip -6pt
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lccc}
\toprule
Kernel & $\mathcal{K}_{\mathsf{PLP}}$ & $\mathcal{K}_{\mathsf{RBC}}$ & $(\mathcal{K}_{\mathsf{PLP}}/\mathcal{K}_{\mathsf{RBC}})^{1/2}$ \\
\midrule
Triangular & 0.95 & 1.33 & 0.84\\
Uniform & 0.83 & 1.13 & 0.86\\
Epanechnikov & 0.85 & 1.25 & 0.83\\
Biweight & 1.01 & 1.41 & 0.85\\
Triweight & 1.15 & 1.55 & 0.86 \\
\bottomrule
\end{tabular*}
\vskip 4pt
{\footnotesize \emph{Notes}: Polynomial order $p=1$.}
\end{table}
Table~\ref{table relative variances} reports the values of $\mathcal{K}_{\mathsf{PLP}}$, $\mathcal{K}_{\mathsf{RBC}}$, and the square root of their ratio for five common kernels and $p=1$, which is by far the most commonly used value in practice. For all the kernels considered, the decrease in interval length is substantial, ranging from $14\%$ to~$17\%$. Figure~\ref{figure equiv_ker_int} provides a comparison of the functions $\mathsf{w}_{\mathsf{PLP}}$ and $\mathsf{w}_{\mathsf{RBC}}$ for the five kernels previously considered and $p=1$. While neither function is always closer to zero than the other, $\mathsf{w}_{\mathsf{PLP}}(u)$ clearly shows less variation around zero, resulting in $\mathcal{K}_{\mathsf{PLP}}<\mathcal{K}_{\mathsf{RBC}}$.
\begin{figure}[t]
\vskip -6pt
\includegraphics[trim={2.5cm 0 2cm 0},clip,width=1\textwidth]{equiv_ker_int.jpg}
\vskip -4pt
\caption{The functions $\mathsf{w}_{\mathsf{PLP}}$ (\rule[0.5ex]{1em}{0.8pt}) and $\mathsf{w}_{\mathsf{RBC}}$ (\rule[0.5ex]{0.3em}{0.8pt}\hspace{0.3em}\rule[0.5ex]{0.3em}{0.8pt}\hspace{0.3em}\rule[0.5ex]{0.3em}{0.8pt}) for five kernels}
{\footnotesize \emph{Notes}: Polynomial order $p=1$.}
\label{figure equiv_ker_int}
\end{figure}
\section{inference at the boundary and rdd}
\label{Section boundary}
In this section, we extend the results in Section~\ref{Section main} to address inference on $g(\mathsf{x})$ when $\mathsf{x}$ is a boundary point. This includes standard applications of nonparametric regression to RDD, which is prominent in empirical work. As it turns out, the key challenge when $\mathsf{x}$ is on the boundary is that the LP bootstrap bias is no longer centered around the original bias~$B_n$, not even asymptotically. Therefore, Assumption~\ref{assn 1}(ii) breaks down and Theorem~\ref{Th RBC PLP equivalence} no longer applies; see Section~\ref{sec failure}. We show that a simple modification of the prepivoting method ($\mathsf{mPLP}$) restores the validity results established in Section~\ref{Section main} for interior points. Importantly, confidence intervals based on the proposed $\mathsf{mPLP}$ procedure automatically account for cases where $\mathsf{x}$ is at the boundary, while remaining asymptotically equivalent to the $\mathsf{PLP}$ procedure from Section~\ref{Section main} when $\mathsf{x}$ is an interior point. The proposed modification is presented in Section~\ref{sec adaptive}, with its application to regression discontinuity designs discussed in Section~\ref{sec RDD}.
\subsection{Challenges of standard prepivoting at the boundary}
\label{sec failure}
A crucial assumption underlying the results of Section~\ref{Section main-theory} is Assumption~\ref{assn 1}(ii), which requires that $\hat{B}_{\mathsf{LP},n} - B_n$ converges in distribution to a mean-zero Gaussian random variable. This assumption fails when $\mathsf{x}$ is a boundary point since $\hat{B}_{\mathsf{LP},n} - B_n$ is no longer (asymptotically) centered at zero. The reason is that, in contrast with the interior case, $C_n$ and its LP bootstrap analogue $C_{\mathsf{LP},n}$ do not converge to the same constant $C$ when $\mathsf{x}$ is on the boundary. Instead, Lemma~\ref{lemmaLPbias} shows that
\begin{equation}
\label{eq Bhat An boundary}
\hat{B}_{\mathsf{LP},n} - B_n - A_n = \xi_{2n,\mathsf{LP}} +\operatorname*{o_{\mathrm{p}}} (1), \quad A_n \operatorname*{\to_{\mathrm{p}}} A := \sqrt{\kappa} g^{(p+1)}(\mathsf{x}) (C_{\mathsf{LP}} - C) /(p+1)!,
\end{equation}
where $\xi_{2n,\mathsf{LP}}:= (nh)^{-1/2} \sum_{i=1}^n w_{\mathsf{LP}\text{-}\mathsf{bc},i} (\mathsf{x}) \varepsilon_i$ has zero mean, $A_n:=\sqrt{nh^{2p+3}}{g^{(p+1)}(\mathsf{x})}({C}_{\mathsf{LP},n}-C_n)/(p+1)!$, $C_{\mathsf{LP},n}:=(nh)^{-1}\sum_{i=1}^{n}w_i(\mathsf{x})C_n(x_i)$ is a smoothed version of~$C_n$, and $g^{(p+1)}(\mathsf{x})$ denotes the (directional) ($p+1$)-order derivative at~$\mathsf{x}$. Importantly, $C \neq C_{\mathsf{LP}}$ when $\mathsf{x}$ is on the boundary, but they are uniquely determined by the kernel~$K$ and the polynomial order~$p$; see Lemma~\ref{lemmaKernelConstants}.
Our modified method restores validity of the prepivoted bootstrap confidence interval by simply rescaling the bootstrap statistic $T^*_{\mathsf{LP},n}$ with a known function of the data, say~$Q_n$. The scaling factor $Q_n$ is chosen such that the term $A_n$ is eliminated from the asymptotic distribution of the bootstrap bias. As we shall see, it does not involve any additional tuning parameters or unknown quantities.
\subsection{a boundary-adaptive prepivoting method}
\label{sec adaptive}
Let $T^*_{\mathsf{mLP},n}:=Q_nT^*_{\mathsf{LP},n}$, where $Q_n:=C_n/C_{\mathsf{LP},n}$ depends only on $K$, $h$, and~$\mathscr{X}_n$. By definition, $\hat{B}_{\mathsf{mLP},n} :=\operatorname{\mathbb{E}}^* [T^*_{\mathsf{mLP},n}] = Q_n \hat{B}_{\mathsf{LP},n}$ with $\hat{B}_{\mathsf{LP},n}$ as in Section~\ref{Section main-theory}. In contrast with $T^*_{\mathsf{LP},n}$, Assumption~\ref{assn 1}(ii) holds for $T^*_{\mathsf{mLP},n}$ for any $\mathsf{x}$, including boundary points, because we have defined $Q_n$ precisely such that $Q_n C_{\mathsf{LP},n} = C_n$. Specifically,
\begin{equation}
\label{eq:BmLP-B}
\hat{B}_{\mathsf{mLP},n}-B_n=\sqrt{nh^{2p+3}}\tfrac{g^{(p+1)}(\mathsf{x})}{(p+1)!}(Q_nC_{\mathsf{LP},n}-C_n) + Q_n\xi_{2n,\mathsf{LP}}+\operatorname*{o_{\mathrm{p}}} (1)=Q_n\xi_{2n,\mathsf{LP}}+\operatorname*{o_{\mathrm{p}}} (1).
\end{equation}
Consequently, the modified LP bootstrap bias $\hat{B}_{\mathsf{mLP},n}$ is asymptotically centered at~$B_n$, as required in Assumption~\ref{assn 1}(ii). Crucially, this result holds for any $\mathsf{x} \in \mathbb{S}_x$, including boundary points.
Taken together, these results imply that an asymptotically valid modified confidence interval, denoted $\mathrm{CI}_{\mathsf{mPLP}\hspace{-0.5pt},n}$, can be constructed provided a consistent estimator $\hat{v}_{\mathsf{mPLP},n}^2$ is available for the asymptotic variance $v^2_{\mathsf{mPLP}}$ of $T_n-\hat{B}_{\mathsf{mLP},n}$. By relying on arguments similar to those used for $\hat{v}^2_{\mathsf{PLP},n}$ in \hyperref[{vPest}]{\textup{\tagform@{\ref*{vPest}}}}, we can show that a consistent estimator of $v_{\mathsf{mPLP}}^2$ is given by
\begin{align} \label{eq:hatvmPLP}
\hat{v}_{\mathsf{mPLP},n}^2:= \frac{1}{{nh}} \sum_{i=1}^n w_{\mathsf{mPLP},i}^2 (\mathsf{x}) \tilde\varepsilon_i^2 ,\quad w_{\mathsf{mPLP},i} (\mathsf{x}) := w_i(\mathsf{x}) - w_{\mathsf{mLP}\text{-}\mathsf{bc},i} (\mathsf{x}),
\end{align}
where $w_{\mathsf{mLP}\text{-}\mathsf{bc},i} (\mathsf{x}):=Q_n w_{\mathsf{LP}\text{-}\mathsf{bc},i} (\mathsf{x})$ with $w_{\mathsf{LP}\text{-}\mathsf{bc},i} (\mathsf{x})$ and $\tilde\varepsilon_i$ as defined in Section~\ref{Section main-LP}. The following lemma formalizes this result.
\begin{lemma} \label{Lemma v2mplp}
Under Assumptions~\ref{Ass_g} and~\ref{Ass_K,h}, $\hat{v}_{\mathsf{mPLP},n}^2 \operatorname*{\to_{\mathrm{p}}} v_{\mathsf{mPLP}}^2 >0$.
\end{lemma}
Because $T^*_{\mathsf{mLP},n}$ satisfies Assumption~\ref{assn 1}, we can apply Theorem~\ref{th 1} to conclude that the \textsf{mPLP} interval~$\mathrm{CI}_{\mathsf{mPLP}\hspace{-0.5pt},n}$ is asymptotically valid and equivalent to
\begin{equation*}
\mathrm{CI}_{\mathsf{mPLP} \textsf{-} \mathsf{RBC},n} = [(\hat{g}_n(\mathsf{x})-(nh)^{-1/2} \hat{B}_{\mathsf{mLP},n})\pm z_{1-\alpha /2}(nh)^{-1/2}\hat{v}_{\mathsf{mPLP},n}],
\end{equation*}
which is an RBC-type interval involving no resampling. This interval is based on a new bootstrap bias correction, $\hat{B}_{\mathsf{mLP},n}$, and a new studentization, $\hat{v}_{\mathsf{mPLP},n}$, and it is valid for both interior and boundary points as stated in the next theorem.
\begin{theorem} \label{Th RBC mPLP equivalence}
Under Assumptions~\ref{Ass_g} and~\ref{Ass_K,h}, for any $\mathsf{x} \in \mathbb{S}_{x}$, $\operatorname{\mathbb{P}}(g(\mathsf{x})\in \mathrm{CI}_{\mathsf{mPLP}\hspace{-0.5pt},n})\to 1-\alpha$ and $\mathrm{CI}_{\mathsf{mPLP}\hspace{-0.5pt},n}=\mathrm{CI}_{\mathsf{mPLP} \textsf{-} \mathsf{RBC},n} + \operatorname*{o_{\mathrm{p}}} ((nh)^{-1/2})$. If, in addition, $T_{\mathsf{LP},n}^{\ast}$ is Gaussian conditional on~$\mathscr{D}_n$, then $\mathrm{CI}_{\mathsf{mPLP}\hspace{-0.5pt},n}=\mathrm{CI}_{\mathsf{mPLP} \textsf{-} \mathsf{RBC},n}$~a.s.
\end{theorem}
An important feature of the \textsf{mPLP} approach is that it automatically adapts to the location of the point~$\mathsf{x}$. That is, it coincides with the \textsf{PLP} approach asymptotically when $\mathsf{x}$ lies in the interior of $\mathbb{S}_x$ (where $Q_n\operatorname*{\to_{\mathrm{p}}} Q=1$), and it provides a valid generalization when $\mathsf{x}$ is on the boundary.
Finally, as in Section~\ref{Section main-improved}, we can show that confidence intervals for $g(\mathsf{x})$ are asymptotically shorter when based on \textsf{mPLP} compared to the existing \textsf{RBC} intervals. This result, which holds for both interior and boundary points, follows by comparing the asymptotic variances $v^2_{\mathsf{mPLP}}$ and~$v^2_{\mathsf{RBC}}$.
\begin{corollary} \label{corollary v^2_mP comparisons}
Let Assumptions~\ref{Ass_g} and~\ref{Ass_K,h} hold.
If $\mathsf{x}$ is an interior point then $v^2_{\mathsf{mPLP}} = v^2_{\mathsf{PLP}}$, and the conclusions of Corollary~\ref{corollary v^2_P comparisons} hold.
If $\mathsf{x}$ is a boundary point, then
\begin{equation*}
(v^2_{\mathsf{mPLP}}, v^2_{\mathsf{RBC}})= \frac{\sigma^2 (\mathsf{x})}{f(\mathsf{x})} (\mathcal{K}_{\mathsf{mPLP}}, \mathcal{K}_{\mathsf{RBC}}) ,
\end{equation*}
where $\mathcal{K}_{\mathsf{mPLP}}$ and $\mathcal{K}_{\mathsf{RBC}}$ are functions only of the kernel $K$ and the polynomial order~$p$.
\end{corollary}
In contrast to Corollary~\ref{corollary v^2_P comparisons}, Corollary~\ref{corollary v^2_mP comparisons} applies to any evaluation point~$\mathsf{x}$, including both interior and boundary points. The values of the constants $\mathcal{K}$ depend on the kernel function~$K$, the polynomial order~$p$, as well as on whether $\mathsf{x}$ is in the interior or at the boundary of~$\mathbb{S}_x$. When $\mathsf{x}$ is in the interior of~$\mathbb{S}_x$, these values coincide with those of $\mathcal{K}_{\mathsf{PLP}}$ and $\mathcal{K}_{\mathsf{RBC}}$ in Corollary~\ref{corollary v^2_P comparisons}, implying that Corollary~\ref{corollary v^2_P comparisons} is a special case of Corollary~\ref{corollary v^2_mP comparisons}.
\begin{table}[t]
\caption{Relative (asymptotic) interval length comparison -- boundary case}
\label{table relative variances boundary}
\vskip -6pt
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lccc}
\toprule
Kernel & $\mathcal{K}_{\mathsf{mPLP}}$ & $\mathcal{K}_{\mathsf{RBC}}$ & $(\mathcal{K}_{\mathsf{mPLP}}/\mathcal{K}_{\mathsf{RBC}})^{1/2}$ \\
\midrule
Triangular & 7.17 & 10.29 & 0.84\\
Uniform & 6.62 & 9.00 & 0.86\\
Epanechnikov & 6.78 & 9.82 & 0.83\\
Biweight & 7.67 & 10.87 & 0.84\\
Triweight & 8.54 & 11.87 & 0.85 \\
\bottomrule
\end{tabular*}
\vskip 4pt
{\footnotesize \emph{Notes}: Polynomial order $p=1$.}
\end{table}
When $\mathsf{x}$ is a boundary point, $\mathcal{K}_{\mathsf{RBC}}$ takes different values compared with the interior point case, and $\mathcal{K}_{\mathsf{mPLP}}$ is different from~$\mathcal{K}_{\mathsf{PLP}}$. Table~\ref{table relative variances boundary} is the boundary analogue of Table~\ref{table relative variances} and confirms that, for the five kernels previously considered, the modified prepivoted interval~$\mathrm{CI}_{\mathsf{mPLP}\hspace{-0.5pt},n}$ is asymptotically shorter than $\mathrm{CI}_{\mathsf{RBC},n}$. As in Figure~\ref{figure equiv_ker_int}, we also plot in Figure~\ref{figure equiv_ker_bnd} the equivalent kernels of \textsf{mPLP} and \textsf{RBC} (see Appendix~\ref{app:kernels}). The comparison reveals that the conclusions from Section~\ref{Section main-improved} also hold for \textsf{mPLP}, i.e., $\mathcal{K}_{\mathsf{mPLP}}<\mathcal{K}_{\mathsf{RBC}}$.
\begin{figure}[t]
\vskip -6pt
\includegraphics[trim={2.5cm 0 2cm 0},clip,width=1\textwidth]{equiv_ker_bnd.jpg}
\vskip -4pt
\caption{The functions $\mathsf{w}_{\mathsf{mPLP}}$ (\rule[0.5ex]{1em}{0.8pt}) and $\mathsf{w}_{\mathsf{RBC}}$ (\rule[0.5ex]{0.3em}{0.8pt}\hspace{0.3em}\rule[0.5ex]{0.3em}{0.8pt}\hspace{0.3em}\rule[0.5ex]{0.3em}{0.8pt}) for five kernels -- boundary case}
{\footnotesize \emph{Notes}: Polynomial order $p=1$.}
\label{figure equiv_ker_bnd}
\end{figure}
\subsection{application to regression-discontinuity designs}
\label{sec RDD}
In this section, we apply the modified prepivoting approach from Section~\ref{sec adaptive} to the problem of inference in~RDD. Adopting the standard potential outcome framework, let $y_i(1)$ and $y_i(0)$ denote the potential outcomes for individual $i$ ($i=1, \ldots,n$) with and without the treatment, respectively. For each individual, we observe the associated treatment indicator $d_i$, which equals 1 if unit $i$ is treated (0 otherwise), and the observed outcome, $y_i=y_i(0)+(y_i(1)-y_i(0))d_i$. We also observe the `forcing' variable~$x_i$, a scalar covariate which is not affected by the treatment and determines whether~$d_i=1$. Hence, the observed data are $\mathscr{D}_n := \{(y_i, x_i, d_i): i=1,\dots,n\}$.
We consider (sharp) RDD, where $d_i$ is determined by $x_i$ being above a given known cutoff $\mathsf{x}$; that is, $d_i:=\mathbb{I}_{\{x_i\geq \mathsf{x}\}}$. Interest is in estimating the average treatment effect at the cutoff, namely $\mathsf{ATE}(\mathsf{x}) := \operatorname{\mathbb{E}}[y_i(1)-y_i(0)|x_i= \mathsf{x}] = g_{+}(\mathsf{x})-g_{-}(\mathsf{x}) $, where $g_{+}(x) := \operatorname{\mathbb{E}}[y_i(1)|x_i = x]$ and $g_{-}(x) := \operatorname{\mathbb{E}}[y_i(0)|x_i = x]$ are the regression functions of the potential outcomes. We consider the following modification of Assumption~\ref{Ass_g}, where $\sigma^2_{+}(x) := \operatorname{\mathbb{V}}[y_i(1)|x_i = x]$ and $\sigma^2_{-}(x) := \operatorname{\mathbb{V}}[y_i(0)|x_i = x]$.
\begin{assumption}\label{Ass_g_RDD}
$(y_i,x_i,d_i)$ are i.i.d.\ and $x_i$ has bounded support $\mathbb{S}_x$ and continuous density~$f$. For all $x$ in an open neighborhood of $\mathsf{x}$, it holds that (i)~$f(x)>0$; (ii)~$\sigma_{+}^2$ and $\sigma_{-}^2$ are continuous and $\sup_{x\in\mathbb{S}_x}\operatorname{\mathbb{E}}[y_i^4|x_i=x]<\infty$; (iii)~$g_{+}^{(p+1)}$ and $g_{-}^{(p+1)}$ are H\"{o}lder continuous with exponent~\mbox{$\eta>0$}.
\end{assumption}
This assumption, which will be used for the asymptotic analysis, also identifies $\mathsf{ATE}(\mathsf{x})$ as the jump in the regression function $g$ at the cutoff~$\mathsf{x}$, i.e., as $\tau(\mathsf{x}):= \lim_{x\to \mathsf{x}^{+}}g(x)-\lim_{x\to \mathsf{x}^{-}}g(x)=g_{+}(\mathsf{x})-g_{-}(\mathsf{x}) $; see Hahn et al.\ (2001) and Imbens and Kalyanaraman (2012). Hence, $\tau(\mathsf{x})$ is the parameter of interest, and the average treatment effect can be estimated as the difference of two local polynomial regressions at the cutoff~$\mathsf{x}$. Formally, we can write this estimator~as
\begin{equation} \label{LPrd}
\hat{\tau}_n (\mathsf{x}) := \hat{g}_{+,n} (\mathsf{x}) - \hat{g}_{-,n}(\mathsf{x}),
\end{equation}
where $\hat{g}_{+,n} (x)$ and $\hat{g}_{-,n} (x)$ are defined as in \hyperref[{eq LP estimator}]{\textup{\tagform@{\ref*{eq LP estimator}}}} with $K((x_i-x)/h)$ replaced by $K_{+}((x_i-x)/h):=K((x_i-x)/h)\mathbb{I}_{\{x_i \geq \mathsf{x}\}}$ and $K_{-}((x_i-x)/h):=K((x_i-x)/h)\mathbb{I}_{\{x_i < \mathsf{x}\}}$, respectively.
To construct valid confidence intervals using the modified prepivoting method of Section~\ref{sec adaptive}, let $T_n:=(nh)^{1/2}(\hat{\tau}_n(\mathsf{x})-\tau (\mathsf{x}))$ and its LP bootstrap analogue $T_n^{\ast}:=(nh)^{1/2}(\hat{\tau}_n^{\ast}(\mathsf{x})-\hat{\tau}_n(\mathsf{x}))$, where $\hat{\tau}_n^{\ast}(\mathsf{x})$ is calculated from bootstrap data, $\mathscr{D}_n^*:=\{(y_i^*,x_i^*,d_i^*): i=1,\dots, n\}$. Specifically, $(x_i^*,d_i^*) = (x_i, d_i)$, $i=1,\dots, n$, and the $y^*_i$'s are as in \hyperref[{btsDGPLP}]{\textup{\tagform@{\ref*{btsDGPLP}}}} with $\hat{g}_n(x_i)$ replaced by $\hat{g}_{+,n}(x_i)\mathbb{I}_{\{x_i \geq \mathsf{x}\}}+\hat{g}_{-,n}(x_i)\mathbb{I}_{\{x_i < \mathsf{x}\}}$.
Next, we define the modified bootstrap statistic,~$T_{\mathsf{rd},n}^{\ast}$. Since estimation of $\tau$ requires fitting two distinct LP regressions, it is useful to decompose $T_n=T_{+,n}-T_{-,n}$, where $T_{+,n}:=(nh)^{1/2}(\hat{g}_{+,n}(\mathsf{x})-g_{+}(\mathsf{x}))$ and $T_{-,n}:=(nh)^{1/2}(\hat{g}_{-,n}(\mathsf{x})-g_{-}(\mathsf{x}))$. In the same way, we can decompose $T_n^{\ast}=T_{+,n}^{\ast}-T_{-,n}^{\ast}$ and define the scaling factors $Q_{+,n}$ and $Q_{-,n}$ as in Section~\ref{sec adaptive} with $K$ replaced by $K_{+}$ and $K_{-}$, respectively. Thus, the modified bootstrap statistic is $T_{\mathsf{rd},n}^{\ast}:=Q_{+,n}T_{+,n}^{\ast}-Q_{-,n}T_{-,n}^{\ast}$.
The bootstrap bias is
\begin{equation}
\hat{B}_{\mathsf{rd},n}:=\operatorname{\mathbb{E}}^{\ast}[T_{\mathsf{rd},n}^{\ast}]=Q_{+,n}\hat{B}_{+,n}-Q_{-,n}\hat{B}_{-,n},
\end{equation}
where $\hat{B}_{+,n}:=\operatorname{\mathbb{E}}^{\ast}[T_{+,n}^{\ast}]$ and $\hat{B}_{-,n}:=\operatorname{\mathbb{E}}^{\ast}[T_{-,n}^{\ast}]$ are computed as in \hyperref[{eq Bhat LP}]{\textup{\tagform@{\ref*{eq Bhat LP}}}} with $K$ replaced by $K_{+}$ and $K_{-}$, respectively, in the definition of the weights~$w_i(x)$. Finally, we need a consistent estimator $\hat{v}_{\mathsf{rd},n}^2$ of the asymptotic variance $v_{\mathsf{rd}}^2$ of $T_n-\hat{B}_{\mathsf{rd},n}$. To this end, we note that $\operatorname{\mathbb{V}} [T_n-\hat{B}_{\mathsf{rd},n}|\mathscr{X}_n] = \operatorname{\mathbb{V}} [T_{+,n}-Q_{+,n}\hat{B}_{+,n}|\mathscr{X}_n]+\operatorname{\mathbb{V}} [T_{-,n}-Q_{-,n}\hat{B}_{-,n}|\mathscr{X}_n]$. Consistent estimators of these two variance components can be obtained as in \hyperref[{eq:hatvmPLP}]{\textup{\tagform@{\ref*{eq:hatvmPLP}}}} with $K$ replaced by $K_{+}$ and $K_{-}$, respectively, thus yielding~$\hat{v}_{\mathsf{rd},n}^2$.
We show in the next theorem that $T^{\ast}_{\mathsf{rd},n}$ satisfies Assumption~\ref{assn 1}, so we can apply Theorem~\ref{th 1} to conclude that a modified confidence interval for $\tau (\mathsf{x})$, denoted $\mathrm{CI}_{\mathsf{rd}\hspace{-0.5pt},n}$ and constructed using $\hat{v}_{\mathsf{rd},n}^2$ as in Sections~\ref{section equivalence} and~\ref{Section main-theory}, is asymptotically valid and equivalent to
\begin{equation*}
\mathrm{CI}_{\mathsf{rd} \textsf{-} \mathsf{RBC},n} = [(\hat{\tau}_n(\mathsf{x})-(nh)^{-1/2} \hat{B}_{\mathsf{rd},n})\pm z_{1-\alpha /2}(nh)^{-1/2}\hat{v}_{\mathsf{rd},n}],
\end{equation*}
which is an RBC-type interval involving no resampling.
\begin{theorem} \label{Th RBC mPLP equivalence rdd}
Under Assumptions~\ref{Ass_K,h} and~\ref{Ass_g_RDD} it holds that (i)~$\hat{v}_{\mathsf{rd},n}^2\operatorname*{\to_{\mathrm{p}}} v_{\mathsf{rd}}^2$, (ii)~$\operatorname{\mathbb{P}}( \tau (\mathsf{x})\in \mathrm{CI}_{\mathsf{rd}\hspace{-0.5pt},n})\to 1-\alpha$, (iii)~$\mathrm{CI}_{\mathsf{rd}\hspace{-0.5pt},n}=\mathrm{CI}_{\mathsf{rd} \textsf{-} \mathsf{RBC},n} + \operatorname*{o_{\mathrm{p}}} ((nh)^{-1/2})$, and (iv)~if, in addition, $T_{\mathsf{rd},n}^{\ast}$ is Gaussian conditional on~$\mathscr{D}_n$, then $\mathrm{CI}_{\mathsf{rd}\hspace{-0.5pt},n}=\mathrm{CI}_{\mathsf{rd} \textsf{-} \mathsf{RBC},n}$~a.s.
\end{theorem}
Finally, we note that all the results in this section can be generalized to allow for different kernel and bandwidth choices on each side of the cutoff.
This changes the efficiency results in Sections~\ref{Section main-improved} and~\ref{sec adaptive} only slightly. Specifically, let $\mathcal{K}_{+,\mathsf{rd}}$ and $\mathcal{K}_{-,\mathsf{rd}}$ denote $\mathcal{K}_{\mathsf{mPLP}}$ for the kernel used to the right and to the left of the cutoff, respectively. We similarly define $\mathcal{K}_{+,\mathsf{RBC}}$ and $\mathcal{K}_{-,\mathsf{RBC}}$. By the i.i.d.\ assumption, the following result then follows from Theorem~\ref{Th RBC mPLP equivalence} and Corollary~\ref{corollary v^2_mP comparisons}.
\begin{corollary} \label{cor RDD}
Under Assumptions~\ref{Ass_K,h} (applied to both sides of the cutoff) and~\ref{Ass_g_RDD},
\begin{equation*}
\frac{v^2_{\mathsf{rd}}}{v^2_{\mathsf{RBC}\textsf{-}\mathsf{rd}}}
= \frac{\sigma_{+}^2(\mathsf{x})\mathcal{K}_{+,\mathsf{rd}}+\sigma_{-}^2(\mathsf{x})\mathcal{K}_{-,\mathsf{rd}}}{\sigma_{+}^2(\mathsf{x})\mathcal{K}_{+,\mathsf{RBC}}+\sigma_{-}^2(\mathsf{x})\mathcal{K}_{-,\mathsf{RBC}}},
\end{equation*}
where $v^2_{\mathsf{RBC}\textsf{-}\mathsf{rd}}$ is the asymptotic variance from the existing $\mathsf{RBC}$ interval for~RDD.
\end{corollary}
As previously, the asymptotic relative length of the confidence intervals is the square root of the ratio given in Corollary~\ref{cor RDD}. Because $\mathcal{K}_{+,\mathsf{rd}}<\mathcal{K}_{+,\mathsf{RBC}}$ and $\mathcal{K}_{-,\mathsf{rd}}<\mathcal{K}_{-,\mathsf{RBC}}$ for all the kernels considered (as seen in Table~\ref{table relative variances boundary}), it follows from Corollary~\ref{cor RDD} that our new confidence intervals $\mathrm{CI}_{\mathsf{rd}\hspace{-0.5pt},n}$ and $\mathrm{CI}_{\mathsf{rd} \textsf{-} \mathsf{RBC},n}$ are asymptotically shorter compared to the existing \textsf{RBC} intervals for RDD. Indeed, when the same kernel is applied on both sides of the cutoff, as is common in applied work, the result in Corollary~\ref{cor RDD} simplifies to that obtained in Corollary~\ref{corollary v^2_mP comparisons} and displayed in Table~\ref{table relative variances boundary}. Specifically, the new intervals are shorter than the existing \textsf{RBC} intervals by the same amount as in Table~\ref{table relative variances boundary}.
It is likely that these results can be extended to other types of RDDs, e.g., fuzzy or kinked~RDD. Although relatively straightforward conceptually, such extensions are nontrivial. For example, the fuzzy RDD estimator is a ratio of two sharp RDD estimators. Thus, limit theory would require joint convergence of those two estimators, which in turn requires a generalization of the prepivoting theory of Cavaliere et al.\ (2024) to vector-valued statistics. With such theory in hand, the fuzzy RDD estimator can be analyzed using the arguments above combined with the delta method.
\section{monte carlo}
\label{Section Monte Carlo}
We now discuss the finite sample performance of the proposed CIs and compare them with the CIs based on existing \textsf{RBC} using Monte Carlo simulation. We also include the invalid (i.e., not prepivoted) CIs for comparison. We consider two distinct inference problems: a nonparametric regression curve evaluated at both an interior and a boundary point and a sharp RDD.
We report results for $5,000$ Monte Carlo replications and nominal level~$0.95$. Estimators are based on the Epanechnikov kernel for the nonparametric regression setup and on the triangular kernel for the RDD setup, as those represent popular kernel choices. Two relevant bandwidth choices are considered: the infeasible MSE-optimal bandwidth~($h$), as a theoretical benchmark, and the feasible (plug-in based) coverage-error-optimal bandwidth~($\hat h$) of Calonico et al.\ (2018, 2020, 2022), representing the reference bandwidth if the aim is to minimize the coverage error of \textsf{RBC} intervals. For each bandwidth choice, we report the average bandwidth~($\bar h$) across Monte Carlo replications. We report empirical coverage probabilities and average lengths for four methods: $\mathrm{CI}_{\mathsf{GP}\hspace{-0.5pt},n}$, $\mathrm{CI}_{\mathsf{LP}\hspace{-0.5pt},n}$, $\mathrm{CI}_{\mathsf{RBC},n}$, and~$\mathrm{CI}_{\mathsf{mPLP}\hspace{-0.5pt},n}$. The first two methods are based on the interval $\mathrm{CI}_{\mathsf{Boot},n}$ in~\hyperref[{eq general invalid BS}]{\textup{\tagform@{\ref*{eq general invalid BS}}}}, using the GP and LP bootstrap schemes, respectively. $\mathrm{CI}_{\mathsf{RBC},n}$ is the interval \hyperref[{eq RBC int}]{\textup{\tagform@{\ref*{eq RBC int}}}} and $\mathrm{CI}_{\mathsf{mPLP}\hspace{-0.5pt},n}$ is the interval defined in Section~\ref{sec adaptive}. For all methods, we use HC3 residuals. The bootstrap intervals are all based on a Gaussian wild bootstrap scheme, and hence implemented analytically without resampling; see Section~\ref{section equivalence}. For the same reason, results for the prepivoted GP interval $\mathrm{CI}_{\mathsf{PGP}\hspace{-0.5pt},n}$ defined in Section~\ref{section GP} would be identical to $\mathrm{CI}_{\mathsf{RBC},n}$, and are thus not reported.
\begin{table}[t]
\small
\caption{coverage and length of 95\% confidence intervals - nonparametric regression}
\vskip -6pt
\label{Table nonpar}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lrlrrrrrrrrr}
\toprule
& & & & \multicolumn{4}{c}{Coverage} & \multicolumn{4}{c}{Length} \\
\cline{5-8} \cline{9-12}
\multicolumn{1}{c}{eval.} & \multicolumn{1}{c}{$n$} & \multicolumn{1}{c}{$h$} & \multicolumn{1}{c}{$\bar{h}$} &
\multicolumn{1}{c}{\textsf{GP}} &
\multicolumn{1}{c}{\textsf{LP}} &
\multicolumn{1}{c}{\textsf{RBC}} &
\multicolumn{1}{c}{\textsf{mPLP}} &
\multicolumn{1}{c}{\textsf{GP}} &
\multicolumn{1}{c}{\textsf{LP}} &
\multicolumn{1}{c}{\textsf{RBC}} &
\multicolumn{1}{c}{\textsf{mPLP}} \\
\midrule
int & 250 & $h$ & 0.189 & 81.5 & 88.9 & 93.6 & 94.2 & 0.635 & 0.635 & 0.915 & 0.765 \\
& ~ & $\hat h$ & 0.359 & 80.3 & 79.3 & 93.1 & 86.5 & 0.461 & 0.461 & 0.664 & 0.554\\
& 500 & $h$ & 0.165 & 82.1 & 89.2 & 94.2 & 94.4 & 0.479 & 0.479 & 0.690 & 0.573 \\
& ~ & $\hat h$ & 0.303 & 81.6 & 85.4 & 94.1 & 91.1 & 0.353 & 0.353 & 0.510 & 0.422 \\
& 1000 & $h$ & 0.143 & 81.9 & 90.0 & 94.4 & 94.9 & 0.361 & 0.361 & 0.521 & 0.431 \\
& ~ & $\hat h$ & 0.256 & 81.9 & 87.7 & 94.6 & 93.8 & 0.270 & 0.270 & 0.390 & 0.322 \\
& 2000 & $h$ & 0.125 &83.3&90.7&95.1&95.4&0.273&0.273&0.394&0.326\\
& & $\hat h$ & 0.217 &82.5&89.4&94.9&94.5&0.207&0.207&0.299&0.247\\
\midrule
bnd & 250 & $h$ & 0.353 & 77.7 & 85.6 & 90.9 & 92.2 & 1.277 & 1.277 & 1.901 & 1.596 \\
& ~ & $\hat h$ & 0.373 & 77.8 & 83.2 & 91.0 & 91.2 & 1.271 & 1.271 & 1.895& 1.602 \\
& 500 & $h$ & 0.307 & 79.9 & 87.6 & 93.3 & 93.8 & 0.961 & 0.961 & 1.426 & 1.191 \\
& ~ & $\hat h$ & 0.327 & 80.6 & 85.6 & 93.2 & 92.8 & 0.946 & 0.946 & 1.403 & 1.174 \\
& 1000 & $h$ & 0.267 & 80.6 & 87.8 & 93.8 & 94.0 & 0.725 & 0.725 & 1.071 & 0.894 \\
& ~ & $\hat h$ & 0.285 & 80.7 & 86.7 & 93.9 & 93.3 & 0.711 & 0.711 & 1.051 & 0.878 \\
& 2000 & $h $ & 0.233 & 81.4&89.5&94.7&94.8&0.549&0.549&0.812&0.676\\
& & $\hat h$ & 0.257 & 81.5&88.1&94.5&94.3&0.526&0.526&0.779&0.649\\
\bottomrule
\end{tabular*}
\vskip 4pt
\end{table}
\medskip
\noindent \textsc{nonparametric regression.} We consider i.i.d.\ data generated as $y_i =g(x_i)+\varepsilon_i$ with $x_i\sim U_{[-1,1]}$ and $\varepsilon_i \sim N(0,\sigma^2)$, where $\sigma=1$ and
\begin{equation*}
g(x)= \frac{ \sin(3\pi x/2)}{1+18x^2 (\mathrm{sign}(x)+1)}.
\end{equation*}
This DGP was previously considered in Berry, Carroll, and Ruppert (2001), Hall and Horowitz (2013), and Calonico et al.\ (2018, 2022), among others. We consider inference both at an interior point, $\mathsf{x}=-1/3$, and a boundary point,~$\mathsf{x}=-1$, denoted int and bnd, respectively, in Table~\ref{Table nonpar}.
\begin{table}[t]
\small
\caption{coverage and length of 95\% confidence intervals - rdd}
\vskip -6pt
\label{Table rd}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}crlrrrrrrrrr}
\toprule
& & & & \multicolumn{4}{c}{Coverage} & \multicolumn{4}{c}{Length} \\
\cline{5-8} \cline{9-12}
\multicolumn{1}{c}{DGP} & \multicolumn{1}{c}{$n$} & \multicolumn{1}{c}{$h$} & \multicolumn{1}{c}{$\bar{h}$} &
\multicolumn{1}{c}{\textsf{GP}} &
\multicolumn{1}{c}{\textsf{LP}} &
\multicolumn{1}{c}{\textsf{RBC}} &
\multicolumn{1}{c}{\textsf{mPLP}} &
\multicolumn{1}{c}{\textsf{GP}} &
\multicolumn{1}{c}{\textsf{LP}} &
\multicolumn{1}{c}{\textsf{RBC}} &
\multicolumn{1}{c}{\textsf{mPLP}} \\
\midrule
1 & 500 & $h $ & 0.082 & 80.6 & 86.7 & 93.3 & 93.2 & 0.345 & 0.345 & 0.589 & 0.470 \\
& & $\hat h$ & 0.073 & 79.6 & 86.0 & 93.5 & 94.0 & 0.373 & 0.373 & 0.674 & 0.539 \\
& 1000 & $h $ & 0.072 & 81.1 & 87.8 & 94.3 & 94.2 & 0.249 & 0.249 & 0.388 & 0.318 \\
& & $\hat h$ & 0.060 & 81.0 & 86.8 & 93.8 & 93.9 & 0.276 & 0.276 & 0.441 & 0.358 \\
& 2000 & $h $ & 0.063 & 81.2 & 88.0 & 94.5 & 94.7 & 0.185 & 0.185 & 0.279 & 0.231 \\
& & $\hat h$ & 0.049 & 80.4 & 87.8 & 94.4 & 94.4 & 0.210 & 0.210 & 0.321 & 0.265 \\
& 4000 & $h $ & 0.054 & 81.3 & 88.3 & 94.8 & 95.1 & 0.138 & 0.138 & 0.205 & 0.170 \\
& & $\hat h$ & 0.041 & 81.2 & 88.2 & 94.8 & 94.8 & 0.160 & 0.160 & 0.240 & 0.199 \\
\midrule
2 & 500 & $h $ & 0.260 & 80.7 & 88.0 & 94.4 & 94.2 & 0.183 & 0.183 & 0.274 & 0.228 \\
& & $\hat h$ & 0.126 & 80.2 & 86.9 & 94.1 & 94.0 & 0.271 & 0.271 & 0.434 & 0.354 \\
& 1000 & $h $ & 0.226 & 81.6 & 89.0 & 94.7 & 94.8 & 0.136 & 0.136 & 0.202 & 0.168 \\
& & $\hat h$ & 0.115 & 81.9 & 88.7 & 94.9 & 95.1 & 0.194 & 0.194 & 0.297 & 0.244 \\
& 2000 & $h $ & 0.197 & 82.9 & 88.5 & 94.8 & 94.8 & 0.102 & 0.102 & 0.150 & 0.125 \\
& & $\hat h$ & 0.101 & 82.0 & 88.8 & 94.7 & 94.7 & 0.144 & 0.144 & 0.214 & 0.178 \\
& 4000 & $h $ & 0.172 & 81.0 & 89.1 & 94.7 & 94.7 & 0.077 & 0.077 & 0.113 & 0.094 \\
& & $\hat h$ & 0.087 & 82.4 & 89.2 & 94.8 & 95.1 & 0.108 & 0.108 & 0.160 & 0.134 \\
\bottomrule
\end{tabular*}
\vskip 4pt
\end{table}
Results for samples of size $n \in \{250, 500, 1000, 2000\}$ are presented in Table~\ref{Table nonpar}. The numerical evidence supports the anticipated decrease in interval lengths of~17\%, even for small sample sizes, suggesting that the asymptotic efficiency of \textsf{mPLP} is rapidly achieved as \( n \) increases. We observe that \textsf{RBC} and \textsf{mPLP} perform similarly in terms of coverage under both bandwidth choices, with empirical coverage probabilities remaining close to the nominal level. A deviation from the nominal level is detected for \textsf{mPLP} with $\hat{h}$, when $\mathsf{x}$ is an interior point and $n$ is small, but this deviation rapidly vanishes as the sample size increases. Finally, because they are not robust to the `large' bandwidth choices considered, non-prepivoted methods exhibit much smaller average lengths, resulting in severe undercoverage even for large~$n$.
\medskip
\noindent \textsc{rdd.} In the context of RDD, we generate i.i.d.\ data from $y_i = g (x_i)+\varepsilon_i$ with $x_i\sim 2 \mathcal{B}(2,4) -1$ and $\varepsilon_i \sim N(0,\sigma^2)$, where $\mathcal B$ is the Beta distribution and~$\sigma=0.1295$. We consider two choices for the regression function~$g$ based on the datasets in Ludwig and Miller (2007) and Lee (2008), as in Calonico et al. (2014). Specifically, the two DGPs are
\begin{align*}
\text{DGP1:} \quad g (x) &= \begin{cases}
3.71 + 2.30x + 3.28 x^2 + 1.45 x^3 + 0.23 x^4 + 0.03 x^5 & \hspace{2.69cm} \text{if } x < 0 , \\
0.26 + 18.49x - 54.81 x^2 +74.30 x^3 - 45.02x^4 + 9.83x^5 & \hspace{2.69cm} \text{if } x \geq 0 ,
\end{cases}\\
\text{DGP2:} \quad g(x) &= \begin{cases}
0.48 + 1.27x - 0.5 \cdot 7.18 x^2 + 0.7 \cdot 20.21 x^3 + 1.1 \cdot 21.54 x^4 + 1.5 \cdot 7.33 x^5 &\text{if } x < 0 , \\
0.52 + 0.84x - 0.1 \cdot 3.00 x^2 - 0.3 \cdot 7.99 x^3 - 0.1 \cdot 9.01x^4 + 3.56x^5 &\text{if } x \geq 0 ,
\end{cases}
\end{align*}
where, as in Calonico et al.\ (2014), the fifth-order polynomial from the Lee (2008) data in DGP2 is modified to increase the curvature of $g$ and hence increase bias.
Results for samples of size $n\in\{500,1000, 2000,4000\}$ (total sample sizes are twice those in the nonparametric regression example) are presented in Table~\ref{Table rd}. The results confirm the conclusions drawn for the nonparametric regression example: the \textsf{mPLP} method produces shorter intervals than \textsf{RBC} for all the considered sample sizes and bandwidth choices, with coverage levels of both prepivoted methods close to the nominal value. As expected, non-prepivoted methods fail to deliver valid inference.
\begin{table}[t]
\caption{practical implementation guide for \textsf{mPLP}}
\label{tab:guidance}
\vskip -6pt
\begin{tabularx}{\textwidth}{@{}lX@{}}
\toprule
Choice & Recommendation \\
\midrule
Bandwidth & Any that can be used with \textsf{RBC}; e.g.\ coverage-error-optimal \\
\addlinespace
Kernel & Any standard kernel; efficiency gains range from 14\%--17\% across common choices (Epanechnikov, triangular, uniform, biweight, triweight) \\
\addlinespace
Polynomial order & $p=1$ (local linear) recommended; higher orders supported \\
\addlinespace
Interior vs.\ boundary/RDD \quad & Use \textsf{mPLP} (adapts to boundary automatically) \\
\addlinespace
Computation & Fully analytic; no resampling required \\
\addlinespace
Software & \texttt{R} packages at \url{https://pppackages.github.io} \\
\bottomrule
\end{tabularx}
\vskip 4pt
\end{table}
\section{guidance for applied researchers}
\label{Section guidance}
The \textsf{mPLP} method proposed in this paper can be used as an alternative or complement to the standard \textsf{RBC} interval in any nonparametric regression or RDD application. Its implementation requires no additional tuning parameters beyond those already needed for \textsf{RBC}: the same bandwidth~$h$, kernel~$K$, and polynomial order~$p$ are used throughout. \textsf{R} packages implementing our procedures are available at \url{https://pppackages.github.io}.
Table~\ref{tab:guidance} summarizes the key choices and considerations for practitioners. We offer the following specific recommendations.
\medskip
\noindent\textsc{bandwidth.} Our method is compatible with any bandwidth selection rule used with \textsf{RBC}, including the coverage-error-optimal bandwidth of Calonico et al.\ (2018) implemented in the \texttt{rdrobust} package, the MSE-optimal bandwidth, or cross-validation. No other bandwidth is required.
\medskip
\noindent\textsc{kernel.} Our efficiency results hold for all standard kernels. For practitioners already using \textsf{RBC} with the Epanechnikov or triangular kernel, the two most common choices,\textsf{mPLP} delivers 17\% and 16\% shorter intervals, respectively (Tables~\ref{table relative variances} and~\ref{table relative variances boundary}). There is no reason to switch kernels when moving from \textsf{RBC} to~\textsf{mPLP}.
\medskip
\noindent\textsc{polynomial order.} As with \textsf{RBC}, local linear estimation ($p=1$) is the standard choice and the one that requires least smoothness of the conditional mean function. Higher-order polynomials are supported by the theory and can be implemented if one can assume additional smoothness.
\medskip
\noindent\textsc{interior vs.\ boundary/rdd.} We recommend using the modified \textsf{mPLP} interval for both interior evaluation points, boundary points, and RDD (where the cutoff is a boundary point) because it adapts automatically to the boundary and requires no additional input from the user.
\medskip
\noindent\textsc{computation.} Because the bootstrap moments (mean and variance) entering the \textsf{mPLP} interval are available in closed form as functions of the kernel weights and residuals, implementation is fully analytic and requires no resampling.
\section{conclusions}
\label{Section conclusions}
This paper proposes novel procedures for inference in nonparametric regression and regression-discontinuity designs based on non-standard implementations of the bootstrap. New confidence intervals based on the concept of prepivoting (Beran 1987, 1988; Cavaliere et al., 2024) are shown to deliver asymptotically correct coverage under general conditions that allow for the presence of bias and thus do not require undersmoothing. We show that prepivoting different choices of bootstrap DGPs yield different robust bias correction (RBC)-type confidence intervals. This connection between the prepivoting and robust bias correction approaches allows us to identify the specific bootstrap DGP underlying the RBC approach of Calonico et al.\ (2014, 2018). More importantly, it enables us to propose a new alternative approach based on a local polynomial bootstrap algorithm that delivers improved inference. Specifically, the prepivoted local polynomial method yields 14--17\% shorter intervals than existing RBC intervals, depending only on the choice of kernel function.
It is worth noting that, because the conditional moments (expectation and variance) of our reference bootstrap statistics can be computed analytically, the practical implementation of our improved bias correction procedure does not require simulating any bootstrap samples. Although we could use resampling instead of the analytical formulas to make the implementation more automatic, this would be much more computationally costly.
The application of prepivoting to inference in the presence of nonnegligible bias has potential applications beyond those explored here. A natural extension is to sieve regression estimators, which have a long tradition in the nonparametrics literature (see, e.g., Andrews, 1991; Huang, 2003; Chen, 2007; Belloni et al., 2015; Chen and Christensen, 2015). Although bias is often assumed away through undersmoothing conditions, Cattaneo et al.\ (2020) have recently proposed RBC inference methods for the particular case of partitioning-based nonparametric estimators. It would be interesting to explore whether prepivoting can improve inference in this setting.
Other potential applications include two-step semiparametric estimators where the first step involves nonparametric regression, potentially affecting the asymptotic distribution of the second step estimator (e.g., Andrews, 1994; Newey, 1994; Chen, Linton, and van Keilegom, 2003). While this literature has primarily focused on how to adjust the variance of the second step estimator, Cattaneo and Jansson (2018) show that asymptotic bias emerges under `small bandwidth' asymptotics when the first step uses kernel regression. Although they show that a particular bootstrap automatically corrects for this bias, their results assume away smoothing bias. Exploring whether prepivoting can improve inference in this setting would be an interesting contribution. Finally, extensions to time series and high-dimensional frameworks (Gupta and Seo, 2023) or spatial data (Hallin et al., 2004) are also of interest.