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.
101,601 characters
Testing for Heterogeneous Treatment Effects in Regression Discontinuity Designs
\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1}
\if00
{
\title{\bf Testing for Heterogeneous Treatment Effects in Regression Discontinuity Designs}
\author{Xiaojun Song\thanks{Corresponding author. Email: \texttt{[email removed]}. This work was supported by the National Natural Science Foundation of China [Grant Numbers 72373007, 72333001, and 72621002]. The author also gratefully acknowledges the research support from the Center for Statistical Science of Peking University and the Key Laboratory of Mathematical Economics and Quantitative Finance (Peking University) of the Ministry of Education, China.}
\hspace{.2cm} \\
Department of Business Statistics and Econometrics \\
Guanghua School of Management, Peking University \\ \\
Haojiao Zhao\thanks{Email: \texttt{[email removed]}.} \\
Department of Business Statistics and Econometrics \\
Guanghua School of Management, Peking University}
\maketitle
} \fi
\if10
{
\bigskip
\bigskip
\bigskip
\begin{center}
{\LARGE\bf Testing for Heterogeneous Treatment Effects in Regression Discontinuity Designs}
\end{center}
\medskip
} \fi
\bigskip
\begin{abstract}
We propose a nonparametric test for unobserved treatment effect heterogeneity in regression discontinuity designs.
Under the null of no unobserved heterogeneity, a transformed outcome that imputes treated potential outcomes for untreated units must have a continuous conditional distribution at the cutoff.
We convert this implication into an integrated conditional-moment restriction using characteristic functions, thereby allowing the conditional local average treatment effect to be an unrestricted function of covariates.
We derive the asymptotic distribution of the test statistics via a $U$-process and establish the validity of a multiplier bootstrap procedure for calculating critical values.
Monte Carlo experiments show well-controlled size and increasing power. Two empirical applications illustrate how the test distinguishes between heterogeneity explained by observables and that explained by unobserved factors.
\end{abstract}
\noindent
{\it Keywords:} Heterogeneous treatment effect, Regression discontinuity, $U$-process, Multiplier bootstrap.
\vfill
\newpage
\spacingset{1.45}
\begin{bibunit}
\section{Introduction} \label{sec:introduction}
Empirical studies increasingly document heterogeneous treatment effects, motivating targeted policy interventions, yet a fundamental question remains unanswered: does the observed variation reflect heterogeneity explained by observable covariates, or a residual component driven by unobserved factors?
Treatment effects may vary with observed covariates or with latent characteristics that persist even after conditioning on those covariates.
Existing inference methods primarily address the first source. The second source, whether treatment effects still vary within subgroups defined by observed covariates, has received far less formal attention.
This gap matters for policy: if unobserved heterogeneity exists, interventions designed solely on the basis of observable covariates will be incomplete. A formal diagnostic is therefore needed.
Current approaches to exploring treatment effect heterogeneity suffer from two limitations.
First, applied studies typically rely on ad hoc analyses, which are sensitive to model misspecification and multiple-testing concerns.
Second, existing methods are mainly designed for instrumental-variable settings and cannot readily be extended to regression discontinuity (RD) designs, where identification is inherently local at the cutoff.
Moreover, existing methods focus on detecting heterogeneity attributable to observables, leaving formal tests of unobserved treatment effect heterogeneity scarce.
Consequently, applied researchers currently lack a principled diagnostic for assessing unobserved treatment effect heterogeneity, especially in RD designs.
This paper fills this gap by developing a nonparametric test for unobserved treatment effect heterogeneity in RD designs.
The key insight is that if all treatment effect variation is explained by the covariates, then imputing treated potential outcomes for untreated units should yield a transformed outcome with a continuous conditional distribution at the cutoff.
We convert this continuity implication into an integrated conditional moment (ICM) restriction, thereby avoiding parametric restrictions on the functional form of the conditional local average treatment effect (CLATE), denoted by $\tau(X)$ with $X$ the observed covariates.
Unlike existing methods that focus on whether $\tau(X)$ is constant across subgroups, our framework directly targets the presence of unobserved heterogeneity beyond the covariates and is valid regardless of the functional form of $\tau(X)$.
In doing so, the proposed test provides a principled way to assess whether observed covariates can plausibly account for treatment effect variation without imposing restrictions on the functional form.
The proposed test has two theoretical advantages.
First, the proposed process converges at the rate $1/\sqrt{nh}$ (with $n$ the sample size and $h$ the testing bandwidth), the natural rate associated with RD designs.
The ICM approach integrates over the covariates, reducing the testing problem to an effectively one-dimensional smoothing exercise and thereby alleviating the curse of dimensionality.
Second, it exhibits favorable power against local alternatives.
The global-smoothing structure of the ICM statistic can detect Pitman-type alternatives shrinking at the nonparametric rate $1/\sqrt{nh}$.
From a practical standpoint, the test also serves as a diagnostic for covariate sufficiency.
A failure to reject means that the data do not reveal residual heterogeneity within the class of alternatives detectable under the maintained assumptions.
A rejection indicates a violation of the homogeneity restriction, which can be interpreted as residual treatment effect heterogeneity under the maintained assumptions.
Moreover, the test adapts to both sharp and fuzzy RD designs and accommodates both local constant and local linear estimators, enhancing the applicability of the proposed testing framework.
The work most closely related to ours is \cite{Hsu2019}, which tests whether the CLATE is constant across covariate-defined subgroups.
In contrast, our test examines whether the heterogeneity can be plausibly explained by observed covariates.
The two procedures are therefore complementary. One can first use \cite{Hsu2019} to assess whether the CLATE varies with the observed covariates, and then apply our procedure to determine whether those covariates exhaust the heterogeneity.
Together, the two tests provide a systematic diagnostic: the first assesses whether $\tau(X)$ varies with $X$, and the second assesses whether residual treatment effect heterogeneity persists.
In this respect, our testing framework explicitly investigates the existence of unobserved treatment effect heterogeneity, which is difficult to examine and has received limited attention.
The rest of the paper is organized as follows.
Section \ref{sec:literature} reviews the related literature.
Section \ref{sec:framework} introduces the testing framework and the $U$-process formulated under the null hypothesis.
Section \ref{sec:asymptotic} establishes the asymptotic properties of the test statistics.
Section \ref{sec:bootstrap} develops a multiplier bootstrap procedure to compute critical values and proves its validity.
Sections \ref{sec:simulation} and \ref{sec:empirical} examine the finite-sample performance through numerical simulations and empirical applications, respectively.
Section \ref{sec:conclusion} concludes.
Proofs of the theoretical results and supplementary simulation results are provided in the Appendix.
\section{Related Literature} \label{sec:literature}
A growing number of empirical studies have examined heterogeneous treatment effects across a wide range of contexts, such as institutional reforms and college enrollment.
Substantial variation in treatment effects has been reported along observable dimensions including demographic characteristics, employment histories, and family status.
A common feature of these studies is their reliance on ad hoc strategies---splitting the sample into subgroups \citep{LopesdeFonseca2020, Brock2023, Giupponi2023, Chetty2026, Mountjoy2026} or adding interaction terms to regression models \citep{AbouDaher2025, Sorrenti2025}---to probe heterogeneity.
While such approaches are informative, none provides statistical evidence on whether the observable covariates exhaust the sources of treatment effect variation or whether significant unobserved heterogeneity persists.
This methodological gap motivates the formal testing procedure developed in the present paper.
A substantial literature studies treatment effect heterogeneity across various identification frameworks.
Early contributions focus on identification in selection models \citep{Heckman1997, Heckman2001, Heckman2005} and non-separable models \citep{Chesher2003, Imbens2009, Torgovitsky2015}.
A growing econometric literature investigates treatment effect heterogeneity under unconfoundedness \citep{Crump2008, Hsu2017, SantAnna2021, Cai2024} and IV identification \citep{Chang2015, Hsu2023}.
Most existing tests assess heterogeneity attributable to observable covariates, whereas \cite{Hsu2023} directly studies unobserved heterogeneity in an IV model.
Nevertheless, none formally addresses unobserved heterogeneity in fuzzy RD designs.
This paper fills this gap by developing the corresponding testing framework for fuzzy RD designs.
Since the seminal work of \cite{Thistlethwaite1960} and \cite{Angrist1999}, RD designs have been widely used to identify causal treatment effects in observational studies.
Within the RD literature, two strands are most relevant to our work.
The first concerns estimation of (heterogeneous) treatment effects in RD designs, including RD estimation methods with covariates \citep{Calonico2019, Frolich2019, kreiss2023, Calonico2025, Noack2025} and continuous treatments \citep{Dong2023, Xie2024}.
The second strand develops specification tests for the identifying assumptions in RD designs, such as continuity of the running variable density \citep{McCrary2008, Otsu2013, Bugni2021}, monotonicity \citep{Hsu2021, Arai2022, Hsu2024}, and other identification conditions \citep{Dong2018, Bertanha2020a}.
Neither strand, however, offers a formal test for unobserved treatment effect heterogeneity, the question addressed in this paper.
Methodologically, our test contributes to the literature on ICM tests \citep{Bierens1982, Bierens1990, Bierens1997, Delgado2001} by developing an ICM test to detect unobserved treatment effect heterogeneity in fuzzy RD designs.
We adapt and extend the ICM framework to fuzzy RD designs, in which both discontinuity-based identification and compliance behavior pose additional challenges.
The ICM framework allows us to test for unobserved heterogeneity without imposing parametric restrictions, making the test attractive in RD designs.
\section{Testing Framework} \label{sec:framework}
\subsection{Model} \label{subsec:model}
Consider the following nonparametric and nonseparable RD model:
$$Y=g(D, R, X, \epsilon),$$
where $Y\in\mathbb{R}$ is the observed outcome, $D\in\{0,1\}$ is the binary treatment, $R\in\mathbb{R}$ is the running variable, $X\in\mathcal{X} \subset \mathbb{R}^d$ collects observed covariates, and $\epsilon\in\mathbb{R}^{d_\epsilon}$ is a general disturbance vector.
The structural function $g$ is a measurable function that determines $Y$ given $D$, $R$, $X$, and $\epsilon$.
Under the potential outcome framework, the observed outcome $Y$ can be written as
$$ Y = DY_1 + (1-D) Y_0,$$
where the potential outcomes $Y_1 = g(1,R,X,\epsilon)$ and $Y_0 = g(0,R,X,\epsilon)$ cannot be observed simultaneously for the same individual.
In the function $g$, $R$ may affect $Y$ directly, but RD identification exploits only the discontinuous jump in $D$ at the cutoff.
Without loss of generality, define $R=0$ as the cutoff.
Let the propensity score be $p(x,r) = \mathbb{P}(D=1|X=x, R=r) = \mathbb{E}[D|X=x, R=r]$, i.e., the conditional probability of treatment, and let $\mu(x,r) = \mathbb{E}[Y|X=x, R=r]$ denote the conditional expectation function of $Y$.
For any such function $\psi(x,r)$, denote its one-sided limits by $\psi(x,0^-) = \lim_{r\uparrow 0} \psi(x,r)$ and $\psi(x,0^+) = \lim_{r\downarrow 0} \psi(x,r)$.
Under Assumption A3 below, the propensity score $p(x,r)$ is continuous in $r$ for $r\ne 0$, exhibiting a jump at the cutoff, i.e., $p(x,0^-) \ne p(x,0^+)$ for all $x\in\mathcal{X}$.
In fuzzy RD designs, the propensity score exhibits a discontinuous but incomplete jump at the cutoff.
Our identification naturally focuses on compliers, whose treatment status switches from untreated to treated as the running variable crosses the cutoff.\footnote{In sharp RD designs, the treatment $D$ is determined by the running variable $R$: $D=1$ whenever $R\ge0$. The propensity score jumps at the cutoff with $p(x,0^+) - p(x,0^-) = 1$ exactly for all $x\in\mathcal{X}$, indicating a deterministic treatment assignment at the cutoff. The conditional average treatment effect at the cutoff can be identified as $\mu(x,0^+)-\mu(x,0^-)$. In this sense, sharp RD is a degenerate special case in which all units near the cutoff are compliers. The theory of this paper is established for the more general fuzzy RD designs, but the testing procedure still accommodates sharp RD with proper estimators.}
The following assumptions guarantee identification of the CLATE in fuzzy RD designs.
\textbf{Assumption A1 (Continuity)} The conditional density $f_{R|X}(r|x)$ is continuous at $r=0$ uniformly in $x\in\mathcal X$ and strictly positive in a neighborhood of $r=0$.
\textbf{Assumption A2 (Probability)} The probability $\mathbb{P}(R\ge0 \mid X=x)$ is bounded away from zero and one for all $x\in \mathcal{X}$.
\textbf{Assumption A3 (Monotonicity)} The treatment is determined by a threshold function $D = \mathbf{1}(\theta(X,R)\ge\eta)$ with an unknown function $\theta$ that is continuous in $r$ on each side of the cutoff and satisfies $\theta(x,0^+)>\theta(x,0^-)$ for all $x\in\mathcal{X}$.
The latent variable satisfies $(\eta,\epsilon) \perp R \mid X$ and the conditional CDF $F_{\eta|X}(\cdot|x)$ is strictly increasing for each $x\in\mathcal X$.
Assumptions A1--A3 combine standard RD overlap and continuity requirements with the stronger latent-index restriction in Assumption A3 for the structural characterization.
For Assumption A1, the continuity of the conditional density $f_{R|X}(r|x)$ at $r=0$ ensures that units on either side of the cutoff are locally comparable for every $x\in\mathcal{X}$.
Assumption A2 ensures that comparable groups exist on both sides of the cutoff for every $x\in\mathcal{X}$.
Assumption A3 rules out defiers, units whose treatment status switches from treated to untreated as $R$ crosses the cutoff.
Under the monotone selection mechanism in Assumption A3, units can be partitioned into always-takers, never-takers, and compliers via $\eta$.
The three groups are defined as $\mathcal{A}_x = \left\{ \eta\in\mathbb{R}: \eta \le \theta(x,0^-) \right\}$, $\mathcal{N}_x = \left\{ \eta\in\mathbb{R}: \eta > \theta(x,0^+) \right\}$, and $\mathcal{C}_x = \left\{ \eta\in\mathbb{R}: \theta(x,0^-) < \eta \le \theta(x,0^+) \right\}$.
This mechanism exhausts the support of $\eta$ and inherently rules out defiers.
Consequently, the propensity score is $p(x,r) = F_{\eta|X}(\theta(x,r)|x)$, and the strict increase of $F_{\eta|X}(\cdot|x)$ together with $\theta(x,0^+)>\theta(x,0^-)$ implies $p(x,0^+)>p(x,0^-)$ for every $x\in\mathcal X$.
The discontinuity in the propensity score at the cutoff confirms the identifiability of the CLATE $\tau(x)$ with a nonzero denominator.
Under Assumptions A1--A3, the CLATE is identified at the cutoff:
$$\tau(x) = \frac{\mu(x,0^+)-\mu(x,0^-)}{p(x,0^+)-p(x,0^-)} = \mathbb{E}[Y_1 - Y_0 | X=x, R=0, \eta\in\mathcal{C}_X]. $$
As is standard, the denominator is the proportion of compliers, while the numerator corresponds to the intent-to-treat effect.
The CLATE $\tau(X)$ is the average treatment gain for the compliers conditional on the covariate $X$, excluding the always-takers and never-takers from the identification.
The above identification formula is inherently local: it is valid only at the cutoff, where the discontinuity in $p(x,r)$ provides exogenous variation.
It loses validity outside the cutoff, where the discontinuity disappears, and the denominator becomes zero.
The remainder of this section develops a formal test for unobserved heterogeneity.
\subsection{Hypothesis of interest} \label{subsec:hypothesis}
Building on the CLATE identification above, we now formulate the hypothesis of interest.
The CLATE summarizes the average treatment effect within each covariate cell but is silent about the dispersion of individual treatment effects within the cell.
Detecting unobserved heterogeneity, therefore, requires examining the distribution of treatment effects, not merely their mean.
In this paper, we test whether the heterogeneity in treatment effects is driven solely by observable covariates $X$, or, alternatively, whether the unobservable factor $\epsilon$ also plays a role.
Formally, the null hypothesis of no unobserved heterogeneity is
$$\mathbb{H}_0: \mathbb{P}[ Y_{1i} - Y_{0i} = \tau(X_i) | X_i, R_i = 0] = 1 \quad a.s., $$
i.e., conditional on $X_i$, the difference in potential outcomes at the cutoff is deterministic.
Conversely, the alternative hypothesis posits that unobserved disturbances also contribute to treatment effect heterogeneity.
Within the structural function $g$, the null hypothesis $\mathbb{H}_0$ is equivalent to an additively separable structure at the cutoff: the unobserved term does not interact with $D$ to produce heterogeneity.
The following proposition provides a structural characterization of the null hypothesis, thereby formalizing a testable implication.
\begin{proposition} \label{prop1}
The null hypothesis of no unobserved heterogeneity in treatment effects, i.e., for some measurable function $\tau(\cdot): \mathcal{X} \mapsto \mathbb{R}$,
$$ \mathbb{H}_0: g(1,0,X,\cdot) - g(0,0,X,\cdot) = \tau(X)$$
holds if and only if the function $g$ at $R=0$ is additive in $D$ and $\epsilon$, i.e.,
$$ g(D,0,X,\epsilon) = m(D,X) + \nu(X,\epsilon), $$
where $m: \mathcal{S}_{DX} \mapsto \mathbb{R}$ and $\nu: \mathcal{S}_{X\epsilon} \mapsto \mathbb{R}$ are measurable functions with $\mathcal{S}$ denoting the support of the relevant variables.
\end{proposition}
Proposition \ref{prop1} adapts the result of \cite{Lu2014} to RD designs: the absence of unobserved heterogeneity is structurally equivalent to additive separability in the function $g$ at the cutoff.
Crucially, the separability property is required only at the cutoff, which is a weaker restriction than the global structural conditions widely used in the literature.
At the cutoff ($R=0$), $Y$ is determined by $D$ and $\epsilon$ given the covariate $X$. If the structural function $g$ satisfies the separability condition in Proposition \ref{prop1}, the treatment effect at the cutoff is $g(1,0,X,\epsilon) - g(0,0,X,\epsilon) = m(1,X) - m(0,X) = \tau(X)$.
Proposition~\ref{prop1} establishes a structural equivalence, but additive separability of $g$ is not directly testable: the structural function $g$ is latent and cannot be identified from data alone.
Assumptions A4 and A5 below translate this structural condition into a testable restriction that can be verified from the observable data.
\textbf{Assumption A4 (Single-index error)}
The structural function $g(d,r,x,\epsilon)$ is continuous in $r$ in a neighborhood of $r=0$ for all $d,x,\epsilon$.
There exists a measurable function $\tilde{g}: \mathcal{S}_{DX} \times \mathbb{R} \mapsto \mathbb{R}$ and a scalar-valued function $\nu: \mathcal{S}_{X\epsilon} \mapsto \mathbb{R}$ such that $g(D,0,X,\epsilon) = \tilde{g}(D,X,\nu(X,\epsilon))$, where $\tilde{g}$ strictly increases in $\nu$.
\textbf{Assumption A5 (Support invariance)} For each $d\in\{0,1\}$ and $x\in\mathcal{X}$, the support of $g(d,0,x,\epsilon)$ conditional on compliers $\mathcal{C}_x$ equals the support in the overall population conditional on the covariate $X$, i.e., $\mathcal{S}_{g(d,0,x,\epsilon) | X=x, \eta\in\mathcal{C}_x} = \mathcal{S}_{g(d,0,x,\epsilon) | X=x }$.
Assumptions A4 and A5 jointly transform the structural property of Proposition \ref{prop1} into an observable restriction at the cutoff.
Assumption A4 is a direct primitive condition on $g$, which complements the continuity of the density in Assumption A1.
When $d_\epsilon>1$, Assumption A4 reduces the multidimensional unobserved heterogeneity to a scalar index, which is the economically relevant source of heterogeneity for our test.
Note that Assumption A4 is naturally satisfied under $\mathbb{H}_0$ since $g(D,0,X,\epsilon) = m(D,X) + \nu(X,\epsilon) $ provides an explicit single-index representation that is monotone in $\nu$.
Monotonicity is essential for identifying the latent scalar index $\nu$ from the observable distribution, and it is a standard regularity condition in nonseparable models \citep{Chesher2003, Matzkin2003}.
Under the alternative, however, the single-index and monotonicity restrictions jointly determine the class of detectable departures from the null.
The proposed test has power against alternatives where the treatment effect depends on a scalar unobservable in a monotone fashion; for instance, when the structural function takes the form $g(d,0,x,\epsilon) = m(d,x) + h(d,x) \cdot \nu(x,\epsilon)$ with $h(0,x)$ and $h(1,x)$ nonzero, of the same sign, and different from each other, so that treatment effects vary systematically with $\nu$ in one direction.
Conversely, the test may lack power when unobserved heterogeneity arises from separate latent factors entering the treated and untreated potential outcomes, as in a sharp RD design with $g(d,0,x,\epsilon) = m(d,x) + d \cdot \epsilon_1 + (1-d) \cdot \epsilon_2$, where $\epsilon_1$ and $\epsilon_2$ are independent and identically distributed.
In this case, the treatment effect $Y_1 - Y_0$ depends on $(\epsilon_1, \epsilon_2)$ and cannot be compressed into a single scalar index while preserving monotonicity.\footnote{Because $\epsilon_1$ and $\epsilon_2$ share the same distribution, however, the conditional distribution of the transformed outcome $W=Y+(1-D)\tau(X)$ defined below remains continuous at the cutoff, so the test may have limited power against this form of heterogeneity.}
Thus, Assumption A4 narrows the alternatives, translating the structural null hypothesis into a testable implication and distinguishing the alternatives from the null.
An important consequence is that the test controls size under $\mathbb{H}_0$ without requiring Assumption A4 to hold under the alternative: the additively separable structure implied by the null automatically satisfies the single-index condition.
This creates a useful asymmetry: size control is robust to violations of A4, while power is concentrated against alternatives that respect the scalar-index structure.
Assumption A5 is a support assumption first introduced in \cite{Vuong2017} in the IV setting.
It ensures that the range of potential outcomes realized by compliers covers that of the whole population.
It also implies that $\mathcal{S}_{g(d,0,x,\epsilon) | X=x, \eta\in\mathcal{C}_x} = \mathcal{S}_{g(d,0,x,\epsilon) | D=d, X=x }$ since $ \mathcal{S}_{g(d,0,x,\epsilon) | X=x, \eta\in\mathcal{C}_x} \subseteq \mathcal{S}_{g(d,0,x,\epsilon) | D=d, X=x} \subseteq \mathcal{S}_{g(d,0,x,\epsilon) | X=x}. $
This condition guarantees that testing the continuity condition for the compliers is equivalent to testing the continuity condition for the whole sample, thereby enabling a global diagnostic of unobserved heterogeneity.
The distribution of $g(d,0,x,\epsilon)$ among compliers can be identified as
\begin{align*}
& \mathbb{P}\left[ g(d,0,x,\epsilon)\le y | X=x, \eta\in\mathcal{C}_x \right] \\
=& \frac{\lim_{r\downarrow0} \mathbb{P}\left[ Y\le y, D=d | X=x, R=r \right] - \lim_{r\uparrow0} \mathbb{P}\left[ Y\le y, D=d | X=x, R=r \right]}{\lim_{r\downarrow0} \mathbb{P}\left[ D=d | X=x, R=r \right] - \lim_{r\uparrow0} \mathbb{P}\left[ D=d | X=x, R=r \right]} \,\,\text{ for any } y.
\end{align*}
The support $\mathcal{S}_{g(d,0,x,\epsilon) | X=x, \eta\in\mathcal{C}_x}$ can be identified from this distribution. Some primitive conditions may be sufficient to imply Assumption A5, e.g.~$\mathcal{S}_{\epsilon|X=x,\eta\in\mathcal{C}_x} = \mathcal{S}_{\epsilon|X=x}$, which holds if the selection margin $\eta$ affects treatment assignment but does not truncate the support of the unobserved determinants relative to the full population.
This condition is plausible when the complier subpopulation is sufficiently representative of the overall population within each covariate cell, for instance, when the instrument induces a shift in treatment take-up without systematically excluding individuals from the extremes of the unobserved distribution.
In practice, the plausibility of Assumption A5 should be assessed on a case-by-case basis.
If compliers systematically differ from the full population in ways that affect the range of potential outcomes (e.g., only individuals with moderate unobserved gains self-select into treatment), then the support equality may be violated.
When A5 fails, the test may produce spurious rejections: a discontinuity in the conditional distribution of $W$ could emerge not from genuine unobserved treatment effect heterogeneity, but from the truncation of the complier support relative to the full population.
Conversely, in some cases the test could become conservative if support differences mask genuine heterogeneity.
A useful diagnostic is to compare the empirical support of the complier distribution (identified via the RD ratio above) with that of the full sample, and assess whether discrepancies are substantively large.
Developing a formal sensitivity analysis for violations of Assumption A5 remains an important avenue for future research.
Together, Assumptions A4 and A5 yield a valid, observable implication regarding the absence of unobserved treatment effect heterogeneity.
In particular, they are required only at the cutoff, rather than across the entire support of $R$.
To derive a testable implication from Proposition \ref{prop1}, we define the transformed outcome $W = Y + (1-D)\tau(X)$.
This construction has a natural interpretation: adding $(1-D)\tau(X)$ to $Y$ imputes the treated potential outcome $Y_1$ for untreated units.
For treated units ($D=1$), $W = Y = Y_1$ directly; for untreated units ($D=0$), $W = Y_0 + \tau(X) = Y_1$ under $\mathbb{H}_0$ at the cutoff.
Intuitively, if $W$ were discontinuous at the cutoff, its density function would exhibit a jump.
Based on Proposition \ref{prop1} and Assumptions A1--A5, the following proposition further connects the null hypothesis and the continuity condition in RD designs.
\begin{proposition} \label{prop2}
Suppose Assumptions A1--A5 hold. Then $\mathbb{H}_0$ holds if and only if the conditional distribution of $W$ is continuous at the cutoff for all $(w,x)$, i.e., $$f_{W|XR}(w|x,0^+) = f_{W|XR}(w|x,0^-) \text{ for all } (w,x).$$
\end{proposition}
Proposition~\ref{prop2} transforms the additive structure into a testable condition in RD designs.
The ``only if'' direction (from $\mathbb{H}_0$ to continuity) follows from Proposition~\ref{prop1} together with Assumptions A1--A3 and the continuity of $g$ in $r$ imposed in the first part of Assumption A4; the latter rules out a jump in the distribution of $W$ that is unrelated to unobserved heterogeneity.
The ``if'' direction (from continuity back to $\mathbb{H}_0$) additionally requires the single-index monotonicity in the second part of Assumption A4 and the support invariance of Assumption A5, as detailed in the proof in the Appendix.
Under $\mathbb{H}_0$, $W$ equals $Y_1$ at the cutoff.
The continuity of the conditional density function is equivalent to the continuity of $Y_1$'s conditional density, which is stricter than the continuity of the conditional expectation commonly used in the literature since the unobserved heterogeneity cannot be captured by the expectation.
Proposition \ref{prop2} exploits such distributional continuity as a characterization of the null implied by the absence of unobserved heterogeneity.
Furthermore, denote the conditional characteristic function of $W$ by $\varphi_{W|XR}(w|x,r)$, which is the Fourier transform of the conditional density function $f_{W|XR}(w|x,r)$:
$$\varphi_{W|XR}(w|x,r) = \mathbb{E}[e^{\mathrm{i} wW}|X=x, R=r] = \int e^{\mathrm{i}w\bar{w}} f_{W|XR}(\bar{w}|x,r) d\bar{w}, $$
with $\mathrm{i}$ denoting the imaginary unit in the characteristic function.
Because characteristic functions uniquely determine distributions, the continuity of the conditional density function is equivalent to the continuity of the conditional characteristic function.
Therefore, by Proposition \ref{prop2}, $\mathbb{H}_0$ is equivalent to the continuity of the conditional characteristic functions:
$$\varphi_{W|XR}(w|x,0^+) = \varphi_{W|XR}(w|x,0^-) \text{ for all } (w,x).$$
A natural approach would be to estimate and compare the conditional functions on either side of the cutoff.
However, the required nonparametric estimation introduces the random-denominator problem and suffers from the curse of dimensionality when $X$ is multivariate.
Rather than using local smoothing, we aggregate the discontinuity into a global moment by integrating it against an exponential weight over the distribution of the covariates:
$$U(w,x) = \int e^{\mathrm{i} x'\bar{x}} \left( \varphi_{W|XR}(w|\bar{x}, 0^+) - \varphi_{W|XR}(w|\bar{x}, 0^-) \right) f_{XR}(\bar{x}, 0) d\bar{x}=0 \text{ for all } (w,x).$$
In the spirit of \cite{Bierens1982}, $U(w,x)$ is a global distance that can be recognized as a weighted aggregation of the discontinuity in integral form with an exponential weighting function.\footnote{Readers are referred to \cite{stinchcombe1998} and \cite{Escanciano2006}, among others, for alternative weighting functions frequently used in the ICM literature.}
Under the null, $\varphi_{W|XR}(w|x,r)$ is continuous at $r=0$ for all $(w,x)$ and the discontinuity at the cutoff equals zero everywhere, leading to $U(w,x)=0$ for all $(w,x)$.
Under the alternatives, the discontinuity of $\varphi_{W|XR}(w|x,r)$ is aggregated in the integral, making $U(w,x)\neq0$ for some $(w,x)$.
Under the ICM approach, $U(w,x)$ captures any discontinuity and avoids the random-denominator problem in estimating the conditional density, reducing the testing problem to an effectively one-dimensional smoothing exercise.
This property preserves all the distributional information required for testing while mitigating the curse of dimensionality.
Note that, by applying the Fourier representation of the conditional density and interchanging the order of integration via Fubini's theorem (the integrand is bounded by the integrable density $f_{WX|R}f_R$), $U(w,x)$ admits the following equivalent representations:
\begin{align*}
U(w,x) =& \int e^{\mathrm{i} x'\bar{x}} \left( \varphi_{W|XR}(w|\bar{x}, 0^+) - \varphi_{W|XR}(w|\bar{x}, 0^-) \right) f_{XR}(\bar{x}, 0) d\bar{x} \\
=& \int e^{\mathrm{i} x'\bar{x}} e^{\mathrm{i} w\bar{w}} \left( f_{W|XR}(\bar{w}|\bar{x}, 0^+) - f_{W|XR}(\bar{w}|\bar{x}, 0^-) \right) f_{XR}(\bar{x},0) d\bar{w}d\bar{x} \\
=& \int e^{\mathrm{i} x'\bar{x}} e^{\mathrm{i} w\bar{w}} \left( f_{WX|R}(\bar{w}, \bar{x}| 0^+) f_R(0^+) - f_{WX|R}(\bar{w}, \bar{x}| 0^-) f_R(0^-) \right) d\bar{w}d\bar{x} \\
=& f_R(0) \left[ \lim_{r\downarrow0} \mathbb{E} \left[ e^{\mathrm{i} \left( x'X + wW \right)} | R=r \right] - \lim_{r\uparrow0} \mathbb{E} \left[ e^{\mathrm{i} \left( x'X + wW \right)} | R=r \right] \right].
\end{align*}
The passage from the third to the fourth line uses the law of iterated expectation, the identity $f_{WX|R}(w,x|r)f_R(r)=f_{WXR}(w,x,r)$, and the fact that $f_R(0^+)=f_R(0^-)=f_R(0)$ under Assumption A1.
This equality of the marginal density of $R$ at the cutoff follows from integrating the continuity of the conditional density $f_{R|X}$ over the compact support $\mathcal{X}$: $f_R(r) = \int f_{R|X}(r|x) f_X(x) dx$, which is continuous at $r=0$ whenever $f_{R|X}(\cdot|x)$ is continuous at $0$ uniformly in $x$.
The final expression in fact shows that $U(w,x)$ equals the jump in the conditional characteristic function of $(X,W)$ at $R=0$, multiplied by $f_R(0)$.
Thus, $U(w,x)=0$ for all $(w,x)$ if and only if the discontinuity is zero almost everywhere, which is the completeness property behind the ICM transformation.
Up to the one-sided kernel normalization, a natural sample analog of $U(w,x)$ is
$$ U_n(w,x) = \frac{1}{n} \sum_{i=1}^n e^{\mathrm{i} \left( x'X_i + wW_i \right)} K_h(R_i) \delta_i,$$
where $\delta_i = \mathbf{1}(R_i\ge 0) - \mathbf{1}(R_i<0)$ encodes the signed one-sided difference and $K_h(\cdot) = h^{-1} K(\cdot/h)$ is the kernel function with bandwidth $h=h_n\to 0$ as $n\to\infty$.
Note that $U_n(w,x)$ is infeasible because $W_i$ cannot be observed for all individuals.
Specifically, for untreated units ($D_i=0$), $W_i = Y_i + \tau(X_i)$ requires the CLATE $\tau(X_i)$, which is unknown and must be estimated.
We therefore construct a feasible counterpart $\hat{U}_n(w,x)$ by plugging in $\hat{W}_i = Y_i + (1-D_i) \hat\tau(X_i)$ with a consistent estimator
$$ \hat\tau(X_i) = \frac{\hat\mu(X_i, 0^+) - \hat\mu(X_i, 0^-)}{\hat p(X_i, 0^+) - \hat p(X_i, 0^-)},$$
where $\hat\mu(X_i,0^\pm)$ and $\hat{p}(X_i,0^\pm)$ are leave-one-out estimators from local constant or local linear regression for the functions $\mu(x,0^\pm)$ and $p(x,0^\pm)$ at the sample points using bandwidth parameters $h_x$ and $h_r$.
Using the estimated values $\hat{W}_i$, we construct the test statistics from the following feasible $U$-process:
$$ \hat{U}_n(w,x) = \frac{1}{n} \sum_{i=1}^n e^{\mathrm{i} \left( x'X_i + w \hat{W}_i \right)} K_h(R_i) \delta_i, $$
whose asymptotic behavior will be investigated in Section \ref{sec:asymptotic}.
\begin{remark}
When the covariates are absent, the testing problem reduces to a test of pure homogeneity of the treatment effect.
The null hypothesis reduces to $$\widetilde{\mathbb{H}}_0: \mathbb{P}[ Y_{1i} - Y_{0i} = \tau | R_i = 0] = 1 \quad a.s., $$
i.e., the treatment effect is homogeneous among individuals at the cutoff.
Let $Y=g(D,R,\epsilon)$ denote the structural function of $Y$, where $\epsilon$ may include some unobservable factors.
The null hypothesis $\widetilde{\mathbb{H}}_0$ can be stated as $g(1,0,\cdot) - g(0,0,\cdot) = \tau$ for some constant $\tau\in\mathbb{R}$.
It holds if and only if the function $g$ is additively separable in $D$ and $\epsilon$ at the cutoff, i.e.,
$$ g(D,0,\epsilon) = m(D) + \nu(\epsilon), $$
where the mean function $m: \{0,1\} \mapsto \mathbb{R}$ and the single-index error function $\nu: \mathcal{S}_{\epsilon} \mapsto \mathbb{R}$ are measurable.
If the analogues of Assumptions A1--A5 hold in the unconditional sense, $\widetilde{\mathbb{H}}_0$ is equivalent to the continuity of the density function $f_{W|R}(w|r)$ at $r=0$ for all $w$ using the imputed treated outcome $W=Y+(1-D)\tau$, which can be assessed by the test statistics constructed from the $U$-process $\hat{U}_n(w) = n^{-1}\sum_{i=1}^n e^{\mathrm{i}w\hat{W}_i} K_h(R_i) \delta_i$.
\end{remark}
\subsection{Test Statistics}\label{subsec:statistics}
Based on the feasible $U$-process $\hat{U}_n(w,x)$, one can construct various types of test statistics to assess how close the process is to 0.
In this paper, we consider two popular choices: Kolmogorov--Smirnov (KS) statistics based on the coordinatewise supremum and Cram\'{e}r--von Mises (CvM) statistics based on the $L_2$ norm.
Both statistics aggregate deviations over the appropriate index set: the KS statistic via a supremum, and the CvM statistic via integration.
They therefore retain sensitivity to a broad class of alternatives while avoiding pointwise nonparametric estimation of the conditional functions.
Neither statistic imposes a parametric model on $\tau(\cdot)$ or on the conditional distribution of $W$.
In practice, we evaluate the $U$-process $\hat{U}_n(w,x)$ on a finite grid whose points are typically the sample points $(\hat{W}_i, X_i)_{i=1}^n$ or a fixed grid $(w_l,x_l)_{l=1}^L$.
Based on this specification, the KS statistic is $\sqrt{nh}$ times the supremum of the maximum of the absolute real and imaginary parts of the process $\hat{U}_n(w,x)$ over a given grid.
For the CvM statistic, a closed-form expression can be derived with a weighting function $g(w,x) =(2\pi)^{-\frac{1+d}{2}} e^{-\frac{w^2 + \Vert x \Vert^2}{2}}$.
Specifically, the formulas for the two types of test statistics are
\begin{align*}
\mathrm{KS}_n = \sqrt{nh} \sup_{(w,x)\in\Omega} \max\left\{ \left\vert \mathrm{Re} (\hat{U}_n(w,x)) \right\vert, \left\vert \mathrm{Im} (\hat{U}_n(w,x)) \right\vert \right\},
\end{align*}
where, here and henceforth, $\Omega$ denotes any compact subset of $\mathbb R\times\mathcal X$, and
\begin{align*}
\mathrm{CvM}_n = \int \left| \sqrt{nh} \hat{U}_n(w,x) \right|^2 g(w,x)dwdx = \frac{h}{n} \sum_{i=1}^n \sum_{j=1}^n e^{-\frac{(\hat{W}_i - \hat{W}_j)^2 + \lVert X_i-X_j \rVert^2}{2}} K_h(R_i)K_h(R_j) \delta_i\delta_j.
\end{align*}
These specifications can simplify computation and improve efficiency.
Alternative norms or weighting functions could also be employed to construct other forms of test statistics. We note a practical consideration: the double-sum expression for $\mathrm{CvM}_n$ above is an exact evaluation of the integrated squared process for any bandwidth configuration, because it is obtained by interchanging the integral with the double sum over the sample; in particular, it does not rely on $h=o(h_r)$.
By contrast, when the estimation effect is non-negligible (i.e., $h=O(h_r)$), the limiting null distribution of $\mathrm{CvM}_n$ depends on the estimation effect, and the corresponding bootstrap statistic $\mathrm{CvM}_n^\ast$ involves quadruple summations (see Section~\ref{sec:bootstrap}), which become computationally expensive for moderate to large sample sizes.
For this reason, we recommend the KS statistic as the primary test when $h=O(h_r)$ and report only KS-based results in the numerical studies.
The CvM statistic remains available as a complementary option whenever $h=o(h_r)$ or when computational resources permit.
In the next section, we give the asymptotic properties of the proposed test statistics.
\section{Asymptotic Theory}\label{sec:asymptotic}
In this section, we study the asymptotic behavior of the proposed test statistics.
First, we investigate the estimation effect of the $U$-process $\hat{U}_n(w,x)$ and obtain its linear representation.
We then analyze the limiting behavior of the $U$-process under the null and the alternatives in Theorems \ref{thm1}--\ref{thm3}, and derive the limiting distributions of the proposed test statistics.
The asymptotic results are organized into three cases reflecting the different regimes of the bandwidth ratio $h/h_r$: cases (i) and (ii) cover the local constant and local linear estimators under $h=O(h_r)$, respectively, where the estimation effect is non-negligible and enters the limiting distribution; case (iii) covers either estimator under $h=o(h_r)$, where the estimation effect is asymptotically dominated by the main term and drops out, yielding a simpler limiting process.
Before presenting these asymptotic results, several technical assumptions are introduced as follows:
\textbf{Assumption B1 (Bandwidth)} As $n\to\infty$, $h, h_r, h_x\to0$, $h/h_r\to c\in[0, \infty)$, $nh\to\infty$, $nh^3\to0$, $nhh_r^2\to0$, $nhh_x^{2l}\to0$, $nh_x^{2d}h_r^2\to\infty$.
\textbf{Assumption B1' (Bandwidth)} As $n\to\infty$, $h, h_r, h_x\to0$, $h/h_r\to c\in[0, \infty)$, $nh\to\infty$, $nh^3\to0$, $nhh_r^4\to0$, $nhh_x^{2l}\to0$, $nh_x^{2d}h_r^2\to\infty$.
\textbf{Assumption B2 (Kernel)} The kernel function $K$ is defined on a compact support. It is nonnegative, symmetric, and of order $l$ (that is $\int v^\alpha K(v)dv =0$ for $\alpha = 1, \cdots, l-1$, and $0 < \int v^l K(v)dv < \infty$). Its roughness is finite: $R_K = \int K^2(v)dv < \infty$.
The bandwidth conditions in Assumptions B1 and B1' control bias and variance for the local constant and local linear CLATE estimators, respectively.
Recall that $h_x$ and $h_r$ are bandwidth parameters for estimating the CLATE $\tau(X)$, which can be flexibly chosen in the regression and can be different from the testing bandwidth $h$ in $U_n(w,x)$.
In particular, under-smoothing is used so that the smoothing bias is asymptotically negligible.
The bandwidth condition $nh_x^{2d}h_r^2\to\infty$ for variance replaces the standard rate $nh_x^d h_r\to\infty$ from conditional nonparametric estimation because the density-weighted kernel estimator in $\hat\tau(X_i)$ incorporates an additional $h_x^d h_r$ factor from the density weighting term.
The condition $h/h_r\to c\in[0,\infty)$ allows the testing bandwidth $h$ to be asymptotically proportional to or smaller than the estimating bandwidth $h_r$, preventing the estimation effect from dominating the main term $U_n(w,x)$.
Assumption B1 is more restrictive than Assumption B1' because the local constant estimator has a larger boundary bias at the cutoff. This is the cost of the simple local constant estimator.
For transparency, write $h\asymp n^{-a}$, $h_r\asymp n^{-b_r}$, and $h_x\asymp n^{-b_x}$.
A sufficient exponent formulation of Assumption B1 is
$$ 0<a<1, \quad 3a>1, \quad a+2b_r>1, \quad a+2lb_x>1, \quad 2db_x+2b_r<1, \quad a\ge b_r,$$
with $a+4b_r>1$ replacing $a+2b_r>1$ under Assumption B1'.
These inequalities are sufficient, and the feasible set they define depends on the covariate dimension $d$ and the kernel order $l$ in Assumption B2.
When no exponent triple $(a,b_r,b_x)$ satisfies the system with $a=b_r$ ($c\in(0,\infty)$), the bandwidths may instead be selected from a feasible region with $a>b_r$, i.e., $h=o(h_r)$.
That case is covered by Theorem~\ref{thm1}(iii), under which the estimation effect is asymptotically negligible.
In practice, the bandwidth can be selected following \cite{Imbens2012a}, \cite{Calonico2014}, and \cite{Calonico2019} when Assumption B1 or B1' is satisfied. Assumption B2 is standard and is satisfied by commonly used kernels.
\textbf{Assumption B3 (Bound)} The support $\mathcal{X}$ is compact and there exists a positive constant $C$ such that $\inf_{x\in\mathcal{X}} |f_{XR}(x,0)(p(x,0^+) - p(x,0^-))| \ge C^{-1}$.
\textbf{Assumption B4 (Smoothness)} Functions $\varphi_{WD|XR}(w,d|x,r)$ (defined in Lemma \ref{Lemma1}), $p(x,r)$ and $\mu(x,r)$ are continuously differentiable to order $l$ in $x$ and up to second order in $r$ on both sides of $r=0$, uniformly in $x\in\mathcal{X}$.
\textbf{Assumption B5 (Moment)} The outcome has a finite second moment: $\mathbb{E}[Y^2]<\infty$.
Assumption B3 ensures that $f_{XR}(x,0)>0$ almost everywhere in $\mathcal{X}$ and that the denominator of $\tau(x)$ is uniformly bounded away from 0, making the treatment effect $\tau(x)$ identifiable.
The compact support restriction can be relaxed at the cost of stronger conditions on the distribution of $X$.
Assumption B4 is a smoothness condition that guarantees the validity of a Taylor expansion around the cutoff.
Assumption B5 imposes a mild moment condition to ensure that the second moments required in the analysis of the $U$-process are finite.
Although our test is motivated by the infeasible process $U_n(w,x)$, we can only obtain $\hat{U}_n(w,x)$.
Therefore, it is necessary to consider the difference between them when investigating the asymptotic behavior of the test statistics.
The (nonparametric) estimation effect is given by
$$\hat{U}_n(w,x) - U_n(w,x)= \frac{1}{n} \sum_{i=1}^n e^{\mathrm{i} x'X_i} \left( e^{\mathrm{i} w\hat{W}_i} - e^{\mathrm{i} wW_i} \right) K_h(R_i) \delta_i. $$
Its asymptotic properties are given in the following lemma.
\begin{lemma} \label{Lemma1}
(i) Suppose Assumptions A1--A5 and B1--B5 hold. For the local constant estimator $\hat\tau(X_i)$, the estimation effect admits the expression
$$ \hat{U}_n(w,x) - U_n(w,x) = \mathrm{i}w \frac{1}{n} \sum_{j=1}^n e^{\mathrm{i} x'X_j} \kappa(w,X_j) (Y_j-D_j\tau(X_j)) K_{h_r}(R_j) \delta_j + o_p\left(\frac{1}{\sqrt{nh}}\right) $$
uniformly in $(w,x)\in\Omega$, where
$$ \kappa(w,x) = \frac{\varphi_{WD|XR}(w, 0|x, 0^+)-\varphi_{WD|XR}(w, 0|x, 0^-)}{p(x,0^+)-p(x,0^-)}$$
with $\varphi_{WD|XR}(w, d|x, r) = \mathbb{E}[e^{\mathrm{i}wW} \mathbf{1}(D=d) | X=x, R=r]$.
(ii) For the local linear estimator $\hat\tau(X_i)$, under Assumption B1' (instead of Assumption B1), the estimation effect admits the expression
\begin{align*}
\hat{U}_n(w,x) - U_n(w,x) =& \mathrm{i}w \frac{1}{2n} \sum_{j=1}^n e^{\mathrm{i}x'X_j} \kappa(w,X_j) f_{XR}(X_j,0) (u_{Yj} - u_{Dj}\tau(X_j)) K_{h_r}(R_j) \\
& e_0' \left( \mathbf{1}(R_j\ge0)B_+^{-1}(X_j) - \mathbf{1}(R_j<0) B_-^{-1}(X_j) \right) \mathbf{r}(R_j) + o_p\left( \frac{1}{\sqrt{nh}}\right),
\end{align*}
uniformly in $(w,x)\in\Omega$, where $u_{Yj}=Y_j-\mu(X_j,R_j)$ and $u_{Dj}=D_j-p(X_j,R_j)$ are the nonparametric residuals of $Y_j$ and $D_j$. The matrices $B_+(x)$ and $B_-(x)$ are defined as $B_+(x)=f_X(x)\mathbb{E}[\mathbf{r}(R)\mathbf{r}(R)'K_{h_r}(R)\mathbf{1}(R\ge0)|X=x]$ and $B_-(x)=f_X(x)\mathbb{E}[\mathbf{r}(R)\mathbf{r}(R)'K_{h_r}(R)\mathbf{1}(R<0)|X=x]$ with $\mathbf{r}(R) = (1, R/h_r)'$.
\end{lemma}
Lemma \ref{Lemma1} provides an explicit representation of the estimation effect $\hat{U}_n(w,x) - U_n(w,x)$ with either a local constant or a local linear estimator for the CLATE $\tau(\cdot)$.
The order of the estimation effect $(nh_r)^{-1/2}$ is slightly different from the order of the infeasible $U$-process $U_n(w,x)$, namely $(nh)^{-1/2}$. The estimation effect may affect the asymptotic distribution, depending on the relative rates of $h$ and $h_r$.
When $h=h_r$, the estimation effect and $U_n(w,x)$ are of the same order, and thus the estimation effect contributes to the asymptotic distribution of $\hat{U}_n(w,x)$. If $h/h_r \to 0$ (i.e., $h=o(h_r)$), the estimation effect is asymptotically negligible relative to $U_n(w,x)$.
In this sense, Lemma \ref{Lemma1} bridges the feasible $U$-process $\hat{U}_n(w,x)$ and the infeasible $U$-process $U_n(w,x)$, guiding the investigation of the limiting distribution under the null hypothesis $\mathbb{H}_0$, the fixed alternative $\mathbb{H}_1$, and a sequence of local alternatives $\mathbb{H}_{1n}$.
\subsection{Asymptotic Consistency}\label{subsec:null}
\begin{theorem} \label{thm1}
Suppose Assumptions A1--A5 and B2--B5 are satisfied. Under the null hypothesis $\mathbb{H}_0: \varphi_{W|XR}(w|x,0^+) = \varphi_{W|XR}(w|x,0^-), \,\,\, \forall (w,x),$
(i) when Assumption B1 holds for the local constant estimator $\hat\tau(X_i)$, $\hat{U}_n(w,x)$ converges weakly on $\Omega$:
$$ \sqrt{nh} \hat{U}_n(w,x) \Longrightarrow U_{\infty,1}(w,x),$$
where, here and henceforth, $\Longrightarrow$ denotes weak convergence;
(ii) when Assumption B1' holds for the local linear estimator $\hat\tau(X_i)$, $\hat{U}_n(w,x)$ converges weakly on $\Omega$: $$ \sqrt{nh} \hat{U}_n(w,x) \Longrightarrow U_{\infty,2}(w,x);$$
(iii) when Assumption B1 holds for the local constant estimator $\hat\tau(X_i)$ or Assumption B1' holds for the local linear estimator $\hat\tau(X_i)$, if $h=o(h_r)$, $\hat{U}_n(w,x)$ converges weakly on $\Omega$:
$$ \sqrt{nh} \hat{U}_n(w,x) \Longrightarrow U_{\infty,0}(w,x).$$
Here, $U_{\infty,1}(w,x)$ is a zero-mean Gaussian process with covariance structure $\mathcal{K}_1(w_1, x_1; w_2, x_2) = \mathbb{E} [U_{\infty,1}(w_1,x_1) \bar{U}_{\infty,1}(w_2,x_2)]$, with $\bar{U}_{\infty,1}(w,x)$ denoting the conjugate process of $U_{\infty,1}(w,x)$.
In addition, $U_{\infty,2}(w,x)$ and $U_{\infty,0}(w,x)$ are zero-mean Gaussian processes with covariance structures $\mathcal{K}_2(w_1, x_1; w_2, x_2) = \mathbb{E} [U_{\infty,2}(w_1,x_1) \bar{U}_{\infty,2}(w_2,x_2)]$ and $\mathcal{K}_0(w_1, x_1; w_2, x_2) = \mathbb{E} [U_{\infty,0}(w_1,x_1) \bar{U}_{\infty,0}(w_2,x_2)]$.
In particular, when $h=o(h_r)$ so that the estimation effect is asymptotically negligible, the covariance structure $\mathcal{K}_0(w_1, x_1; w_2, x_2)$ simplifies to
\begin{align*}
\mathcal{K}_0(w_1, x_1; w_2, x_2) = R_K \int e^{\mathrm{i} (x_1-x_2)'x} \varphi_{W|XR}(w_1-w_2|x, 0) f_{XR}(x, 0) dx.
\end{align*}
The expressions for $\mathcal{K}_1(w_1, x_1; w_2, x_2)$ and $\mathcal{K}_2(w_1, x_1; w_2, x_2)$ additionally incorporate terms arising from the estimation effect and are provided in the Appendix.
\end{theorem}
Under $\mathbb{H}_0$, the conditional characteristic function $\varphi_{W|XR}$ is continuous at $r=0$ for all $(w,x)\in\Omega$ and the one-sided limits coincide.
Consequently, $\sqrt{nh} \hat{U}_n(w,x)$ converges weakly to the centered Gaussian process $U_{\infty, j}$ for $j=0,1,2$ under the corresponding under-smoothing conditions.
The covariance structure of the limiting Gaussian process $\mathcal{K}_j$ for $j=0,1,2$ is complicated.
In particular, the covariance structures $\mathcal{K}_1$ and $\mathcal{K}_2$ additionally incorporate the first-order estimation effect and the intrinsic stochastic variation in the RD model.
The convergence rate in Theorem \ref{thm1} is $1/\sqrt{nh}$, instead of the parametric rate $1/\sqrt{n}$, reflecting the local identification of RD designs.
This rate is analogous to that of a one-dimensional boundary kernel estimator, in which the effective sample size is $nh$ rather than $n$.
We note, however, that the first-stage estimation of $\tau(X)$ still requires $d$-dimensional smoothing over the covariates, so the overall procedure shares the curse-of-dimensionality limitations common to all RD methods that involve covariates in the estimation step.
In practice, this restricts the feasible dimension of $X$ to $d=2$ or $3$ with realistic sample sizes, as discussed in Assumptions B1 and B1'.
Thus, procedures for accommodating higher-dimensional covariate vectors, such as semiparametric index restrictions on $\tau(X)$, constitute a useful extension.
Theorem \ref{thm1} and the continuous mapping theorem yield the asymptotic null distributions of continuous functionals of $\sqrt{nh} \hat{U}_n(w,x)$, including our test statistics $\mathrm{KS}_n$ and $\mathrm{CvM}_n$ in the following corollary.
\begin{corollary} \label{corollary}
Under the assumptions of Theorem \ref{thm1} and $\mathbb{H}_0$, for any continuous functional $\mathcal{T}(\cdot)$ (with respect to the supremum norm),
$$\mathcal{T}(\sqrt{nh} \hat{U}_n) \overset{d}{\rightarrow} \mathcal{T}(U_{\infty,j}), $$
for $j=0,1,2$, where $\overset{d}{\rightarrow}$ denotes convergence in distribution. In particular, for the KS and CvM statistics,
\begin{align*}
\mathrm{KS}_n &\overset{d}{\rightarrow} \sup_{(w,x)\in\Omega} \max \left\{ \left| \mathrm{Re} (U_{\infty,j}(w,x))\right|, \left| \mathrm{Im} (U_{\infty,j}(w,x)) \right| \right\}, \\
\mathrm{CvM}_n &\overset{d}{\rightarrow} \int \left| U_{\infty,j}(w,x) \right|^2 g(w,x)dwdx.
\end{align*}
\end{corollary}
The test statistics $\mathrm{KS}_n$ and $\mathrm{CvM}_n$ converge to the images of $U_{\infty,j}(w,x)$ under continuous mappings for $j=0,1,2$, depending on the bandwidth choice and the estimator. Their limits are obtained via the continuous mapping theorem (see, e.g., Theorem 1.3.6 in \cite{vdv1996}), and a detailed proof is provided in the Appendix.
Importantly, the limiting distributions of the test statistics $\mathrm{KS}_n$ and $\mathrm{CvM}_n$ under the null depend in a complex way on the underlying data-generating process through the covariance structure $\mathcal{K}_j(w_1, x_1; w_2, x_2)$ of the centered Gaussian process $U_{\infty,j}$.
Therefore, critical values cannot be tabulated, which motivates the multiplier bootstrap procedure developed in Section \ref{sec:bootstrap}.
\subsection{Asymptotic Power} \label{subsec:power}
\subsubsection{Against the fixed alternative} \label{subsubsec:H1}
The following theorem delivers the asymptotic power against the fixed alternative
$$\mathbb{H}_1: \varphi_{W|XR}(w|x,0^+) - \varphi_{W|XR}(w|x,0^-) =\gamma(w,x), $$
where $\gamma(w,x)\ne0$ for some $(w,x) \in \Omega$.
\begin{theorem} \label{thm2}
Suppose Assumptions A1--A5 and B2--B5 hold. Under $\mathbb{H}_1$, when Assumption B1 holds for the local constant estimator or Assumption B1' holds for the local linear estimator,
$$\sup_{(w,x) \in \Omega} \left| \hat{U}_n(w,x) - \Gamma(w,x) \right| = o_p(1), $$
where $\Gamma(w,x) = \int e^{\mathrm{i} x'\bar{x}} \gamma(w,\bar{x}) f_{XR}(\bar{x},0) d\bar{x}$ is not equal to 0 for some $(w,x)\in\Omega$.
\end{theorem}
Theorem \ref{thm2} indicates that $\sqrt{nh} \hat{U}_n(w,x)$ diverges when $\Gamma(w,x)\ne 0$ for some $(w,x)\in\Omega$.
Under $\mathbb{H}_1$, the uniform law of large numbers guarantees that $\hat{U}_n(w,x)$ converges to the deterministic function $\Gamma(w,x)$ in probability.
Since $\Gamma(w,x) \ne 0$ whenever the characteristic function of $W$ is discontinuous at the cutoff for some $(w,x)$, $\sqrt{nh} \hat{U}_n(w,x)$ diverges at some points $(w,x)\in\Omega$ with $\Gamma(w,x)\ne0$.
Consequently, the test is consistent against fixed alternatives that generate a nonzero ICM signal under the maintained assumptions, with rejection probability approaching 1 as $n$ increases.
\subsubsection{Against the sequence of local alternatives} \label{subsubsec:H1n}
When the alternatives are close to the null hypothesis, the performance of the test may be sharply different from that under the fixed alternative $\mathbb{H}_1$.
We therefore study the asymptotic power under a sequence of local alternatives that approach the null, thereby exploring the theoretical limit of the test.
Consider a Pitman-type sequence of alternatives:
$$\mathbb{H}_{1n}: \varphi_{W|XR}(w|x,0^+) - \varphi_{W|XR}(w|x,0^-) = \frac{\lambda(w,x) }{\sqrt{nh}}, $$
where $\lambda(w,x)\ne 0$ for some $(w,x)\in\Omega$. This sequence of alternatives converges to the null at the rate $1/\sqrt{nh}$, which is the detection boundary of our procedure in the RD context.
The following theorem reveals the asymptotic behavior of $\sqrt{nh}\hat{U}_n(w,x)$ under the sequence $\mathbb{H}_{1n}$ and guarantees non-trivial power against local alternatives of this type.
\begin{theorem} \label{thm3}
Suppose Assumptions A1--A5 and B2--B5 hold. Under the sequence of local alternatives $\mathbb{H}_{1n}$,
(i) when Assumption B1 holds for the local constant estimator $\hat\tau(X_i)$, $$ \sqrt{nh} \hat{U}_n(w,x) \Longrightarrow U_{\infty,1}(w,x) + \Lambda(w,x);$$
(ii) when Assumption B1' holds for the local linear estimator $\hat\tau(X_i)$, $$ \sqrt{nh} \hat{U}_n(w,x) \Longrightarrow U_{\infty,2}(w,x) + \Lambda(w,x);$$
(iii) when Assumption B1 holds for the local constant estimator $\hat\tau(X_i)$ or Assumption B1' holds for the local linear estimator $\hat\tau(X_i)$, if $h=o(h_r)$,
$$ \sqrt{nh} \hat{U}_n(w,x) \Longrightarrow U_{\infty,0}(w,x) + \Lambda(w,x),$$
where the drift term $\Lambda(w,x) = \int e^{\mathrm{i} x'\bar{x}} \lambda(w,\bar{x}) f_{XR}(\bar{x},0) d\bar{x} \ne 0$ for some $(w,x)\in\Omega$.
\end{theorem}
Under $\mathbb{H}_{1n}$, the limiting distribution of $\sqrt{nh} \hat{U}_n(w,x)$ differs from the null limit by the deterministic drift $\Lambda(w,x)\neq0$ for some $(w,x)\in\Omega$.
While the stochastic fluctuation remains asymptotically identical to that under $\mathbb{H}_0$, the additional drift term $\Lambda(w,x)$ introduces a systematic deviation and shifts the limiting distribution away from the centered Gaussian process defined in Theorem \ref{thm1}.
Consequently, the test statistics exhibit non-negligible deviations from $\mathbb{H}_0$ under $\mathbb{H}_{1n}$.
This behavior is characteristic of global-smoothing tests, whose integrated structure allows weak but systematic discrepancies from the null hypothesis to accumulate asymptotically.
In contrast, local-smoothing tests evaluate departures pointwise and require the alternative to converge to the null at a slower rate than the nonparametric smoothing rate.
This difference in mechanism explains the power advantage of global-smoothing tests against local alternatives. In light of the continuous mapping theorem, under $\mathbb{H}_{1n}$, the test statistics satisfy
$$\operatorname{KS}_{n} \overset{d}{\rightarrow} \sup_{(w,x)\in\Omega} \max \left\{ \left\vert \mathrm{Re} \left( U_{\infty,j}(w,x) + \Lambda(w,x) \right) \right\vert, \left\vert \mathrm{Im} \left( U_{\infty,j}(w,x) + \Lambda(w,x) \right) \right\vert \right\}$$
and $$\operatorname{CvM}_{n} \overset{d}{\rightarrow} \int\left| U_{\infty,j}(w,x) + \Lambda(w,x) \right|^2 g(w,x)dwdx $$ for $j=0,1,2$.
The test statistics converge to a non‑centered Gaussian limit, guaranteeing non‑trivial local power at the $1/\sqrt{nh}$ detection boundary.
\section{Critical Values} \label{sec:bootstrap}
In this section, we introduce a multiplier bootstrap procedure to approximate the limiting distributions of the proposed test statistics and calculate their critical values.
In Section \ref{sec:asymptotic}, Lemma \ref{Lemma1} provides linear representations of the $U$-process $\hat{U}_n(w,x)$, and Theorems \ref{thm1}--\ref{thm3} establish its limiting behavior under different hypotheses.
However, the limiting distributions of the test statistics are difficult to evaluate analytically because they depend on the underlying data-generating process in a complex way.
This motivates the use of a multiplier bootstrap procedure, which is both computationally efficient and theoretically valid.
The resulting bootstrap sample allows us to compute the critical values of the test statistics.
In practice, we provide the following multiplier bootstrap procedure for the proposed tests.
\begin{itemize}
\item \textbf{Step 1} \quad
Compute the feasible $U$-process $\hat U_n(w,x)$ using $\hat\tau(X_i)$ and obtain the test statistics $\mathrm{KS}_n$ and $\mathrm{CvM}_n$ as described in Section \ref{sec:framework}.
\item \textbf{Step 2} \quad Draw multipliers $\{V_i\}_{i=1}^n$ and compute the bootstrap statistics $\mathrm{KS}_n^\ast$ and $\mathrm{CvM}_n^\ast$ based on the bootstrap process $\hat{U}_n^\ast(w,x)$.
\item \textbf{Step 3} \quad Repeat Step 2 $B$ times to generate the bootstrap statistics $\{\mathrm{KS}^\ast_{n,b}\}_{b=1}^B$ and $\{\mathrm{CvM}^\ast_{n,b}\}_{b=1}^B$.
\item \textbf{Step 4} \quad Reject the null hypothesis at significance level $\alpha$ if the test statistic exceeds the $(1-\alpha)$ quantile of the bootstrap statistics.
\end{itemize}
In Step 2, the multipliers $\{V_i\}_{i=1}^n$ are i.i.d. random variables with zero mean, unit variance, and bounded support, independent of the original sample path $\{(Y_i, D_i, X'_i, R_i)'\}_{i=1}^n$.
Following \cite{Mammen1993}, one can adopt the i.i.d. two-point Bernoulli random variables $\{V_i\}_{i=1}^n$ with $\mathbb{P}(V_i = 1- \iota) = \iota/\sqrt{5}$ and $\mathbb{P}(V_i = \iota) = 1-\iota/\sqrt{5}$, where $\iota = (\sqrt{5} + 1)/2$. This sequence $\{V_i\}_{i=1}^n$ ensures that $\mathbb{E}[V_i]=0$, $\mathbb{E}[V_i^2]=1$, and $\mathbb{E}[V_i^3]=1$.
As discussed in \cite{Escanciano2014} and \cite{santanna2019}, the multiplier bootstrap avoids repeated estimation in each bootstrap repetition.
This feature makes the multiplier bootstrap straightforward to implement and computationally efficient, compared with alternative resampling methods.
With the multipliers $\{V_i\}_{i=1}^n$, we define the bootstrap process to construct the bootstrap statistics.
Note that the asymptotic linear representation varies across different estimators and bandwidth conditions.
We introduce the bootstrap process and the test statistics here for different scenarios. The bootstrap process based on the local constant estimator is
$$\hat{U}_n^\ast(w,x) = \frac{1}{n} \sum_{i=1}^n V_i e^{\mathrm{i} x'X_i} \left( e^{\mathrm{i} w\hat{W}_i } K_h(R_i) \delta_i + \mathrm{i} w \hat\kappa(w,X_i) (Y_i-D_i\hat\tau(X_i)) K_{h_r}(R_i) \delta_i \right), $$
with the leave-one-out estimator of $\kappa(w,X_i)$ given by
$$ \hat\kappa(w,X_i) = \frac{\frac{1}{n-1} \sum_{j\ne i}^n e^{\mathrm{i}w\hat{W}_j} (1-D_j) K_{h_x}(X_j-X_i) K_{h_r}(R_j) \delta_j}{\frac{1}{n-1} \sum_{j\ne i}^n D_j K_{h_x}(X_j-X_i) K_{h_r}(R_j) \delta_j}.$$
As in $\hat\tau(X_i)$, the product kernel $K_{h_x}(x) = h_x^{-d} \prod_{k=1}^d K(x_k/h_x)$ is employed for the $d$-dimensional vector $x$ in conditional local smoothing.
When $h=o(h_r)$, the estimation effect is asymptotically dominated by the main term, and the bootstrap process based on either the local constant or the local linear estimator reduces to
$$\hat{U}_n^\ast(w,x) = \frac{1}{n} \sum_{i=1}^n V_i e^{\mathrm{i} \left( x'X_i + w\hat{W}_i \right) } K_h(R_i) \delta_i.$$
Based on the bootstrap process $\hat{U}_n^\ast(w,x)$, the bootstrap statistics, $\mathrm{KS}_n^\ast$ and $\mathrm{CvM}_n^\ast$, can be computed as follows:
\begin{align*}
\mathrm{KS}_n^\ast =& \sqrt{nh} \sup_{(w,x)\in\Omega} \max \left\{ \left| \mathrm{Re} (\hat{U}_n^\ast(w,x)) \right|, \left| \mathrm{Im} (\hat{U}_n^\ast(w,x)) \right| \right\}, \\
\mathrm{CvM}_n^\ast =& \frac{h}{n} \sum_{i=1}^n \sum_{j=1}^n V_iV_j e^{-\frac{(\hat{W}_i - \hat{W}_j)^2 + \lVert X_i-X_j \rVert^2}{2}} K_h(R_i)K_h(R_j) \delta_i\delta_j.
\end{align*}
Note that the above expression for $\mathrm{CvM}_n^\ast$ applies when $h=o(h_r)$ and the estimation effect is asymptotically negligible.
Otherwise, the corresponding bootstrap expression of the CvM statistic involves quadruple summations, which are computationally expensive.
The bootstrap samples closely mimic the theoretical distribution of the test statistics under the null, facilitating robust inference in finite samples.
To establish the validity of the multiplier bootstrap, we first investigate the estimation effect of the bootstrap process in the following lemma.
\begin{lemma} \label{Lemma2}
Suppose Assumptions A1--A5 and B2--B5 hold.
(i) Suppose Assumption B1 holds. Then the infeasible bootstrap process with the estimation effect of the local constant estimator is
$$U_n^\ast(w,x) = \frac{1}{n} \sum_{i=1}^n V_i e^{\mathrm{i} x'X_i} \left( e^{\mathrm{i} wW_i } K_h(R_i) \delta_i + \mathrm{i} w \kappa(w,X_i) (Y_i-D_i\tau(X_i)) K_{h_r}(R_i) \delta_i \right), $$
where $W, \kappa$ and $\tau$ in the formula are the true values.
(ii) Suppose Assumption B1 holds for the local constant estimator or Assumption B1' holds for the local linear estimator. When $h=o(h_r)$, the infeasible bootstrap process is
$$U_n^\ast(w,x) = \frac{1}{n} \sum_{i=1}^n V_i e^{\mathrm{i} \left( x'X_i + wW_i \right) } K_h(R_i) \delta_i. $$
In either case (i) or (ii), the difference between the bootstrap process $\hat{U}_n^\ast(w,x)$ and its infeasible version $U_n^\ast(w,x)$ is uniformly bounded:
$$ \sqrt{nh} \sup_{(w,x)\in\Omega} \left| \hat{U}_n^\ast(w,x) - U_n^\ast(w,x) \right| = o_p(1). $$
\end{lemma}
\textbf{Remark} \quad
Lemma \ref{Lemma2} establishes the asymptotic equivalence between the feasible and infeasible bootstrap processes.
The equivalence holds for the local constant estimator under Assumption B1 and, when $h=o(h_r)$, for both the local constant and local linear estimators.
The corresponding result for the local linear estimator when $h=O(h_r)$ can be established, but is not demonstrated here to avoid excessive technical complexity.
The next theorem establishes the asymptotic validity of the multiplier bootstrap procedure above.
\begin{theorem} \label{thm4}
Suppose Assumptions A1--A5 and B2--B5 hold.
(i) Suppose Assumption B1 holds. Then the bootstrap process $\hat{U}_n^\ast(w,x)$ with the local constant estimator converges weakly on $\Omega$ under $\mathbb{H}_0$, $\mathbb{H}_1$, or $\mathbb{H}_{1n}$:
$$\sqrt{nh} \hat{U}_n^\ast(w,x) \Longrightarrow_\ast U_{\infty,1}(w,x), $$
where $\Longrightarrow_\ast$ denotes weak convergence under the bootstrap law.
(ii) Suppose Assumption B1 holds for the local constant estimator or Assumption B1' holds for the local linear estimator. When $h=o(h_r)$, the bootstrap process $\hat{U}_n^\ast(w,x)$ converges weakly on $\Omega$ under $\mathbb{H}_0$, $\mathbb{H}_1$, or $\mathbb{H}_{1n}$:
$$\sqrt{nh} \hat{U}_n^\ast(w,x) \Longrightarrow_\ast U_{\infty,0}(w,x). $$
\end{theorem}
The bootstrap validity established above is guaranteed by the conditional multiplier central limit theorem for $U$-processes.
An interesting feature of Theorem~\ref{thm4} is that the bootstrap process converges weakly to the same centered Gaussian process under $\mathbb{H}_0$, $\mathbb{H}_1$, and $\mathbb{H}_{1n}$.
This unification is uncommon in conventional settings, where the bootstrap process under the alternative typically converges to a different limit, although still a centered Gaussian process.
The distinctive behavior here follows from the nonparametric convergence rate $1/\sqrt{nh}$.
Under $\mathbb{H}_1$, the feasible $U$-process $\hat{U}_n(w,x)$ contains a nonzero deterministic component $\Gamma(w,x)$, which diverges when scaled by $\sqrt{nh}$.
However, the key insight is that the multiplier bootstrap $U$-process perturbs each summand by a mean-zero multiplier $V_i$; the signal $\Gamma(w,x)$ resides in the empirical mean of the summands, and because $\mathbb{E}[V_i]=0$ it does not enter the bootstrap mean.
Formally, conditional on the original sample, the expectation of the bootstrap process $\hat{U}_n^\ast(w,x)$ is approximately zero regardless of whether $\mathbb{H}_0$ or $\mathbb{H}_1$ holds, because $V_i$ is independent of the sample path with zero mean.
Therefore, the signal $\Gamma(w,x)$ does not enter its mean. The component of the bootstrap process that carries the signal is of order $\sqrt{h}$ (and thus asymptotically negligible).
The detailed calculation in the Appendix shows that the bootstrap limit is governed by the same functional form as the limiting null process, with the covariance kernel $\mathcal{K}_j$ of Theorem~\ref{thm1}.
Consequently, $\sqrt{nh}\,\hat{U}_n^\ast(w,x)$ converges to the same centered Gaussian limit under $\mathbb{H}_1$ as under $\mathbb{H}_0$.
A practical consequence is that the bootstrap critical values remain stochastically bounded regardless of which hypothesis holds, while the test statistics diverge at the rate $\sqrt{nh}$ under $\mathbb{H}_1$, ensuring that the rejection probability of the proposed test approaches 1 against any fixed alternative as $n$ increases.
This phenomenon highlights a distinct feature of nonparametric inference in RD designs: the local identification at the cutoff forces a slow convergence rate, which in turn produces a multiplier bootstrap procedure that is uniformly valid across the null, the alternative, and the local alternatives.
\section{Numerical Study} \label{sec:simulation}
In this section, Monte Carlo experiments are conducted to evaluate the finite-sample performance of the proposed test.
We consider a set of data-generating processes (DGPs) adapted from \cite{Hsu2019, Hsu2021}.
In all DGPs, the running variable $R$, covariate $X$, and error term $\epsilon$ are independently generated from the following distributions:
$$ R \sim 2\text{Beta}(2,2)-1, \quad\quad X \sim \text{Unif}[0,1], \quad\quad \epsilon \sim N(0,1). $$
For each DGP, the experiments are repeated 1000 times.
We follow the finite-sample guidance on bandwidth selection in \cite{Calonico2014} and \cite{Calonico2019}, which provide MSE-optimal bandwidths for the RD estimator.
For the estimating bandwidths $h_r$ and $h_x$, we start with the MSE-optimal bandwidth $\hat{h}_{\text{MSE}}$ obtained from the procedure of \cite{Calonico2014} and apply under-smoothing: we set $h_r = h_x = \hat{h}_{\text{MSE}} \times n^{1/5-1/k_1}$ with $k_1 \in [4, 4.5]$, where larger $k_1$ produces a larger bandwidth.
For the testing bandwidth $h$, we set $h = k_2 h_r$ with $k_2 \in [0.9, 1.1]$. When $k_2 = 1$, the estimation effect and the main $U$-process contribute equally to the limiting distribution; when $k_2 < 1$, the estimation effect is partially attenuated relative to the main term.
For the configurations used in our simulations, the effective sample size $nh$ ranges from approximately 100 to 700.
The triangular kernel function $K(u) = (1-|u|) \mathbf{1}(|u|\le 1)$ is used for both the testing and estimating procedures, following the recommendations in \cite{Imbens2012a}, \cite{Calonico2014}, \cite{Calonico2019} and \cite{Hsu2019}.
A total of 1000 bootstrap simulations are conducted within each repetition to compute the critical values.
The test statistics are calculated using the $U$-process, which incorporates the estimation effect of the local constant estimator.
The KS statistic is calculated over a grid of sample points.
Additional results using the local linear estimator are reported in the Appendix.
The outcome $Y$ and binary treatment $D$ are generated according to DGP-specific mechanisms listed below:
DGP1: \textit{Homogeneous Zero Treatment Effects}.
\begin{align*}
Y &= -0.555 - 0.553X + 0.581R + 0.060XR - 0.058R^2 + 1.074X^2 + 0.1\epsilon. \\
D &= \mathbf{1}(R\ge0) \mathbf{1}(0.596 - 2.103X + 0.128R + 0.352XR + 0.013R^2 + 2.454X^2 + \epsilon > 0).
\end{align*}
The treatment effect under DGP1 is $\tau(x)=0$.
DGP2: \textit{Homogeneous Nonzero Constant Treatment Effects}.
\begin{align*}
Y &= \begin{cases}
-0.373 + 0.545R - 0.056R^2 + 0.1\epsilon, & R\ge0, \\
-0.531 + 0.556R - 0.192R^2 + 0.1\epsilon, & R<0. \end{cases} \\
D &= \mathbf{1}(R\ge0) \mathbf{1}(0.331 + 0.277R + 0.049R^2 + \epsilon > 0).
\end{align*}
The nonzero constant treatment effect in DGP2 is
$$\tau(x) = \frac{\mu(x,0^+)-\mu(x,0^-)}{p(x,0^+)-p(x,0^-)}= \frac{0.158}{1-\Phi(-0.331)} \approx 0.25,$$
where $\Phi(\cdot)$ denotes the CDF of the standard normal distribution.
DGP3: \textit{Homogeneous Functional Treatment Effects}.
\begin{align*}
Y &= \begin{cases}
-0.921 - 4X + 0.584R - 0.054R^2 + 5X^2 + 0.1\epsilon, & R\ge0, \\
-0.705 + 0.264X + 0.580R + 0.191R^2 + 0.1\epsilon, & R<0. \end{cases} \\
D &= \mathbf{1}(R\ge0) \mathbf{1}(0.331 + 0.277R + 0.049R^2 + \epsilon > 0).
\end{align*}
The functional treatment effect in DGP3 is
$$\tau(x)=\frac{5x^2-4.264x-0.216}{1-\Phi(-0.331)}.$$
DGP4: \textit{Heterogeneous Treatment Effects}: a 1:1 mixture of DGP1 (zero treatment effect) and DGP3 (functional treatment effect).
DGP5: \textit{Heterogeneous Treatment Effects}: a 1:1 mixture of DGP2 (constant treatment effect) and DGP3 (functional treatment effect).
DGP6: \textit{Heterogeneous Treatment Effects}: a 1:1:1 mixture of DGP1 (zero treatment effect), DGP2 (nonzero constant treatment effect) and DGP3 (functional treatment effect).
DGPs 1--3, adapted from \cite{Hsu2019, Hsu2021}, represent scenarios with zero, nonzero constant, and functional treatment effects, respectively, all in the absence of unobserved heterogeneity.
These DGPs are used to evaluate the size of the proposed test under the null hypothesis of no unobserved heterogeneity.
DGPs 4--6 are constructed as mixtures of DGPs 1--3 to introduce unobserved heterogeneity and assess the test's power under complex scenarios.
Specifically, for each individual $i$, a latent group membership variable $G_i$ is drawn independently of $(R_i, X_i, \epsilon_i)$, with the appropriate probabilities to match the stated mixture proportions (e.g., $G_i \in \{1,3\}$ each with probability $1/2$ for DGP4).
For a given value of $G_i$, the pair $(Y_i, D_i)$ is generated from the corresponding component DGP using a common draw of the base variables $(R_i, X_i, \epsilon_i)$. In all mixture DGPs, the treatment effect at the cutoff is $Y_{1i} - Y_{0i} = \tau_{G_i}(X_i)$, where $\tau_g(X_i)$ denotes the treatment effect under component DGP $g$. Because $G_i$ is independent of the observables, the conditional local average treatment effect $\tau(X_i) = \mathbb{E}[Y_{1i} - Y_{0i} \mid X_i, R_i=0]$ is the weighted average of the component-specific conditional treatment effects.
Nevertheless, conditional on $X_i$, the individual-level treatment effect $Y_{1i} - Y_{0i}$ still varies with the unobserved $G_i$, thus violating $\mathbb{H}_0$. Note that each component $g$ specifies $g(d,0,x,\epsilon) = m_g(d,x) + 0.1\epsilon$, which is additively separable in $\epsilon$ with a monotone index $\nu(x,\epsilon) = \epsilon$.
Across the mixture, however, the group indicator $G$ acts as an additional latent factor that shifts both the outcome level $m_G(0,x)$ and the treatment effect $\tau_G(x)$, so the mixture as a whole need not satisfy the single-index restriction in Assumption A4.
As discussed in Section~\ref{subsec:hypothesis}, Assumption A4 is used to characterize the null hypothesis and to delineate the class of alternatives against which the test has power; size control under $\mathbb{H}_0$ does not require it to hold under the alternative, and the mixture DGPs are therefore used purely as alternatives that violate the continuity restriction on the imputed outcome $W$.
DGP4, which mixes zero and functional treatment effects, creates the sharpest contrast between components and is expected to yield the highest power; DGP5, which mixes constant and functional effects, represents a subtler form of heterogeneity; DGP6 combines all three components, further diluting the signal and posing the greatest challenge to the test.
These mixture scenarios are challenging but can provide insights for evaluating the power of the proposed test statistics.
The numerical results are presented in Tables \ref{tab:DGP_null_lc} and \ref{tab:DGP_alt_lc}.
\begin{table}[H]
\centering
\caption{Rejection rates under the null hypothesis (DGPs 1--3)}
\label{tab:DGP_null_lc}
\setlength{\extrarowheight}{-1pt}
\begin{adjustbox}{max width=0.99\textwidth, max height=\textheight}
\begin{tabular}{@{}cccccccccccccccccc@{}}
\toprule
\multirow{2}{*}{DGP} & \multirow{2}{*}{$k_1$} & $k_2$ & \multicolumn{3}{c}{0.9} & \multicolumn{3}{c}{0.95} & \multicolumn{3}{c}{1} & \multicolumn{3}{c}{1.05} & \multicolumn{3}{c}{1.1} \\
\cmidrule{4-18}
& & $n$ & 1\% & 5\% & 10\% & 1\% & 5\% & 10\% & 1\% & 5\% & 10\% & 1\% & 5\% & 10\% & 1\% & 5\% & 10\% \\
\midrule
\multirow{12}{*}{1} & \multirow{4}{*}{4} & 500 & 0.004 & 0.031 & 0.057 & 0.004 & 0.032 & 0.056 & 0.005 & 0.032 & 0.059 & 0.005 & 0.032 & 0.062 & 0.004 & 0.036 & 0.065 \\
& & 1000 & 0.006 & 0.031 & 0.065 & 0.006 & 0.030 & 0.068 & 0.008 & 0.031 & 0.068 & 0.008 & 0.036 & 0.070 & 0.011 & 0.038 & 0.072 \\
& & 2000 & 0.011 & 0.048 & 0.089 & 0.011 & 0.051 & 0.091 & 0.013 & 0.055 & 0.093 & 0.013 & 0.053 & 0.095 & 0.012 & 0.055 & 0.097 \\
& & 4000 & 0.015 & 0.057 & 0.092 & 0.017 & 0.056 & 0.092 & 0.016 & 0.057 & 0.095 & 0.018 & 0.061 & 0.098 & 0.018 & 0.065 & 0.100 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.25} & 500 & 0.005 & 0.029 & 0.061 & 0.004 & 0.031 & 0.064 & 0.005 & 0.033 & 0.066 & 0.004 & 0.036 & 0.068 & 0.006 & 0.037 & 0.068 \\
& & 1000 & 0.009 & 0.029 & 0.076 & 0.010 & 0.031 & 0.076 & 0.010 & 0.035 & 0.074 & 0.010 & 0.038 & 0.074 & 0.012 & 0.041 & 0.077 \\
& & 2000 & 0.011 & 0.055 & 0.097 & 0.013 & 0.054 & 0.097 & 0.012 & 0.053 & 0.095 & 0.012 & 0.052 & 0.097 & 0.012 & 0.054 & 0.104 \\
& & 4000 & 0.014 & 0.057 & 0.097 & 0.015 & 0.058 & 0.098 & 0.016 & 0.061 & 0.100 & 0.017 & 0.062 & 0.100 & 0.019 & 0.061 & 0.106 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.5} & 500 & 0.004 & 0.032 & 0.064 & 0.005 & 0.031 & 0.064 & 0.006 & 0.036 & 0.067 & 0.006 & 0.037 & 0.069 & 0.006 & 0.039 & 0.070 \\
& & 1000 & 0.010 & 0.031 & 0.074 & 0.009 & 0.035 & 0.076 & 0.010 & 0.039 & 0.078 & 0.011 & 0.041 & 0.079 & 0.012 & 0.044 & 0.079 \\
& & 2000 & 0.010 & 0.055 & 0.096 & 0.011 & 0.057 & 0.099 & 0.012 & 0.056 & 0.101 & 0.013 & 0.054 & 0.104 & 0.015 & 0.055 & 0.108 \\
& & 4000 & 0.013 & 0.056 & 0.100 & 0.016 & 0.055 & 0.104 & 0.017 & 0.059 & 0.102 & 0.018 & 0.062 & 0.107 & 0.018 & 0.064 & 0.107 \\
\midrule
\multirow{12}{*}{2} & \multirow{4}{*}{4} & 500 & 0.004 & 0.023 & 0.046 & 0.005 & 0.025 & 0.047 & 0.005 & 0.023 & 0.047 & 0.007 & 0.024 & 0.051 & 0.007 & 0.025 & 0.051 \\
& & 1000 & 0.005 & 0.021 & 0.051 & 0.005 & 0.021 & 0.052 & 0.005 & 0.023 & 0.051 & 0.007 & 0.024 & 0.050 & 0.008 & 0.025 & 0.053 \\
& & 2000 & 0.013 & 0.042 & 0.096 & 0.012 & 0.042 & 0.095 & 0.011 & 0.043 & 0.093 & 0.012 & 0.044 & 0.088 & 0.011 & 0.045 & 0.089 \\
& & 4000 & 0.017 & 0.045 & 0.089 & 0.017 & 0.046 & 0.086 & 0.013 & 0.047 & 0.090 & 0.013 & 0.049 & 0.092 & 0.014 & 0.051 & 0.093 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.25} & 500 & 0.005 & 0.023 & 0.050 & 0.005 & 0.023 & 0.052 & 0.007 & 0.025 & 0.054 & 0.006 & 0.024 & 0.056 & 0.007 & 0.025 & 0.056 \\
& & 1000 & 0.006 & 0.025 & 0.053 & 0.007 & 0.026 & 0.058 & 0.009 & 0.026 & 0.058 & 0.009 & 0.027 & 0.057 & 0.010 & 0.028 & 0.057 \\
& & 2000 & 0.010 & 0.046 & 0.097 & 0.011 & 0.044 & 0.092 & 0.011 & 0.048 & 0.090 & 0.011 & 0.048 & 0.089 & 0.011 & 0.049 & 0.089 \\
& & 4000 & 0.013 & 0.048 & 0.092 & 0.013 & 0.050 & 0.090 & 0.012 & 0.052 & 0.094 & 0.014 & 0.054 & 0.096 & 0.014 & 0.057 & 0.100 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.5} & 500 & 0.006 & 0.026 & 0.054 & 0.006 & 0.028 & 0.054 & 0.006 & 0.027 & 0.054 & 0.006 & 0.027 & 0.057 & 0.006 & 0.025 & 0.058 \\
& & 1000 & 0.008 & 0.026 & 0.059 & 0.009 & 0.026 & 0.058 & 0.008 & 0.029 & 0.057 & 0.009 & 0.029 & 0.060 & 0.010 & 0.029 & 0.060 \\
& & 2000 & 0.009 & 0.046 & 0.095 & 0.010 & 0.050 & 0.096 & 0.010 & 0.051 & 0.097 & 0.010 & 0.048 & 0.096 & 0.010 & 0.048 & 0.097 \\
& & 4000 & 0.015 & 0.050 & 0.093 & 0.015 & 0.053 & 0.090 & 0.014 & 0.054 & 0.092 & 0.013 & 0.055 & 0.098 & 0.014 & 0.058 & 0.099 \\
\midrule
\multirow{12}{*}{3} & \multirow{4}{*}{4} & 500 & 0.003 & 0.019 & 0.046 & 0.003 & 0.021 & 0.046 & 0.003 & 0.023 & 0.049 & 0.004 & 0.024 & 0.052 & 0.005 & 0.023 & 0.050 \\
& & 1000 & 0.007 & 0.027 & 0.052 & 0.007 & 0.026 & 0.054 & 0.009 & 0.026 & 0.055 & 0.009 & 0.028 & 0.056 & 0.008 & 0.029 & 0.057 \\
& & 2000 & 0.009 & 0.038 & 0.082 & 0.009 & 0.038 & 0.085 & 0.009 & 0.041 & 0.085 & 0.008 & 0.041 & 0.085 & 0.008 & 0.041 & 0.085 \\
& & 4000 & 0.015 & 0.045 & 0.087 & 0.014 & 0.044 & 0.083 & 0.016 & 0.046 & 0.084 & 0.016 & 0.047 & 0.086 & 0.017 & 0.047 & 0.089 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.25} & 500 & 0.005 & 0.022 & 0.050 & 0.004 & 0.024 & 0.053 & 0.004 & 0.025 & 0.053 & 0.005 & 0.027 & 0.052 & 0.005 & 0.029 & 0.054 \\
& & 1000 & 0.008 & 0.028 & 0.056 & 0.008 & 0.029 & 0.059 & 0.008 & 0.029 & 0.057 & 0.009 & 0.029 & 0.060 & 0.009 & 0.031 & 0.060 \\
& & 2000 & 0.008 & 0.047 & 0.087 & 0.008 & 0.047 & 0.089 & 0.009 & 0.049 & 0.088 & 0.009 & 0.047 & 0.086 & 0.010 & 0.048 & 0.086 \\
& & 4000 & 0.016 & 0.050 & 0.086 & 0.015 & 0.048 & 0.091 & 0.015 & 0.047 & 0.094 & 0.015 & 0.049 & 0.096 & 0.015 & 0.052 & 0.097 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.5} & 500 & 0.004 & 0.026 & 0.050 & 0.004 & 0.029 & 0.053 & 0.004 & 0.032 & 0.054 & 0.004 & 0.030 & 0.057 & 0.005 & 0.031 & 0.056 \\
& & 1000 & 0.008 & 0.026 & 0.063 & 0.008 & 0.028 & 0.060 & 0.008 & 0.031 & 0.064 & 0.009 & 0.029 & 0.062 & 0.010 & 0.029 & 0.067 \\
& & 2000 & 0.008 & 0.052 & 0.097 & 0.009 & 0.053 & 0.096 & 0.009 & 0.050 & 0.097 & 0.011 & 0.049 & 0.098 & 0.013 & 0.049 & 0.097 \\
& & 4000 & 0.014 & 0.048 & 0.098 & 0.012 & 0.049 & 0.096 & 0.013 & 0.048 & 0.096 & 0.013 & 0.053 & 0.097 & 0.015 & 0.057 & 0.098 \\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}
Table \ref{tab:DGP_null_lc} reports the rejection rates for DGPs 1--3 under the null hypothesis $\mathbb{H}_0$.
The results indicate that the proposed test maintains an appropriate size across varying values of the tuning parameters $k_1$ for the estimating bandwidth and $k_2$ for the testing bandwidth.
The rejection rates approach the nominal levels as the sample size grows under all DGPs, demonstrating size control and robustness of the proposed test.
As is typical for nonparametric tests based on kernel methods, the performance is affected by the choice of bandwidth.
Table \ref{tab:DGP_alt_lc} reports the rejection rates for DGPs 4--6 under the alternative hypothesis $\mathbb{H}_1$.
In all cases, the rejection rates increase with sample size. The rejection rate is relatively insensitive to the choice of $k_2$, but depends on $k_1$.
As $k_1$ increases, the rejection rates increase.
As established in Lemma \ref{Lemma1}, $k_1$ determines the magnitude of the estimation effect, thereby affecting the performance of the test.
A comparison of empirical power across DGPs 4--6 shows that the rejection rates under DGP4 are consistently the highest given the same bandwidth setting, because it creates the sharpest contrast between components.
DGP5 (constant versus functional) and DGP6 (three-way mixture) exhibit lower power, consistent with their more diluted heterogeneity signals, as anticipated in the DGP designs.
Overall, the results in Tables \ref{tab:DGP_null_lc}--\ref{tab:DGP_alt_lc} demonstrate that the proposed test maintains appropriate size under $\mathbb{H}_0$ and achieves increasing power under $\mathbb{H}_1$, which aligns closely with the theoretical properties in Section \ref{sec:asymptotic}.
\begin{table}[H]
\centering
\caption{Rejection rates under the alternative hypothesis (DGPs 4--6)}
\label{tab:DGP_alt_lc}
\setlength{\extrarowheight}{-1pt}
\begin{adjustbox}{max width=0.99\textwidth, max height=\textheight}
\begin{tabular}{@{}cccccccccccccccccc@{}}
\toprule
\multirow{2}{*}{DGP} & \multirow{2}{*}{$k_1$} & $k_2$ & \multicolumn{3}{c}{0.9} & \multicolumn{3}{c}{0.95} & \multicolumn{3}{c}{1} & \multicolumn{3}{c}{1.05} & \multicolumn{3}{c}{1.1} \\
\cmidrule{4-18}
& & $n$ & 1\% & 5\% & 10\% & 1\% & 5\% & 10\% & 1\% & 5\% & 10\% & 1\% & 5\% & 10\% & 1\% & 5\% & 10\% \\
\midrule
\multirow{12}{*}{4} & \multirow{4}{*}{4} & 500 & 0.002 & 0.025 & 0.067 & 0.001 & 0.028 & 0.069 & 0.002 & 0.033 & 0.072 & 0.002 & 0.039 & 0.077 & 0.002 & 0.041 & 0.085 \\
& & 1000 & 0.034 & 0.235 & 0.406 & 0.034 & 0.252 & 0.422 & 0.043 & 0.274 & 0.435 & 0.055 & 0.287 & 0.443 & 0.057 & 0.296 & 0.449 \\
& & 2000 & 0.398 & 0.748 & 0.830 & 0.429 & 0.762 & 0.833 & 0.445 & 0.767 & 0.838 & 0.457 & 0.769 & 0.842 & 0.471 & 0.767 & 0.845 \\
& & 4000 & 0.925 & 0.982 & 0.990 & 0.927 & 0.982 & 0.989 & 0.930 & 0.983 & 0.988 & 0.932 & 0.984 & 0.987 & 0.931 & 0.985 & 0.987 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.25} & 500 & 0.002 & 0.034 & 0.087 & 0.002 & 0.038 & 0.095 & 0.002 & 0.043 & 0.103 & 0.003 & 0.048 & 0.112 & 0.004 & 0.057 & 0.124 \\
& & 1000 & 0.053 & 0.323 & 0.505 & 0.060 & 0.345 & 0.517 & 0.076 & 0.366 & 0.528 & 0.087 & 0.378 & 0.532 & 0.101 & 0.393 & 0.537 \\
& & 2000 & 0.525 & 0.821 & 0.888 & 0.562 & 0.820 & 0.899 & 0.571 & 0.823 & 0.899 & 0.591 & 0.828 & 0.906 & 0.607 & 0.833 & 0.906 \\
& & 4000 & 0.960 & 0.991 & 0.996 & 0.962 & 0.991 & 0.996 & 0.960 & 0.992 & 0.996 & 0.964 & 0.991 & 0.995 & 0.963 & 0.991 & 0.995 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.5} & 500 & 0.002 & 0.040 & 0.120 & 0.003 & 0.046 & 0.131 & 0.003 & 0.056 & 0.147 & 0.004 & 0.066 & 0.155 & 0.005 & 0.071 & 0.173 \\
& & 1000 & 0.079 & 0.397 & 0.600 & 0.093 & 0.417 & 0.604 & 0.108 & 0.444 & 0.613 & 0.122 & 0.457 & 0.616 & 0.149 & 0.472 & 0.621 \\
& & 2000 & 0.646 & 0.868 & 0.929 & 0.662 & 0.875 & 0.933 & 0.675 & 0.879 & 0.936 & 0.688 & 0.887 & 0.939 & 0.698 & 0.889 & 0.939 \\
& & 4000 & 0.979 & 0.996 & 0.999 & 0.983 & 0.996 & 0.999 & 0.982 & 0.996 & 0.999 & 0.981 & 0.996 & 0.998 & 0.983 & 0.996 & 0.998 \\
\midrule
\multirow{12}{*}{5} & \multirow{4}{*}{4} & 500 & 0.001 & 0.020 & 0.050 & 0.001 & 0.020 & 0.055 & 0.001 & 0.023 & 0.059 & 0.003 & 0.027 & 0.063 & 0.003 & 0.027 & 0.068 \\
& & 1000 & 0.015 & 0.129 & 0.273 & 0.016 & 0.142 & 0.285 & 0.021 & 0.150 & 0.299 & 0.027 & 0.154 & 0.313 & 0.030 & 0.168 & 0.317 \\
& & 2000 & 0.228 & 0.585 & 0.727 & 0.244 & 0.597 & 0.735 & 0.263 & 0.608 & 0.735 & 0.280 & 0.615 & 0.740 & 0.294 & 0.630 & 0.740 \\
& & 4000 & 0.836 & 0.951 & 0.972 & 0.841 & 0.954 & 0.972 & 0.843 & 0.956 & 0.974 & 0.846 & 0.956 & 0.978 & 0.852 & 0.956 & 0.978 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.25} & 500 & 0.001 & 0.025 & 0.063 & 0.002 & 0.028 & 0.073 & 0.003 & 0.030 & 0.079 & 0.004 & 0.031 & 0.086 & 0.004 & 0.034 & 0.094 \\
& & 1000 & 0.024 & 0.185 & 0.354 & 0.028 & 0.200 & 0.372 & 0.035 & 0.211 & 0.382 & 0.044 & 0.222 & 0.391 & 0.048 & 0.243 & 0.402 \\
& & 2000 & 0.347 & 0.685 & 0.815 & 0.361 & 0.701 & 0.825 & 0.382 & 0.709 & 0.827 & 0.397 & 0.712 & 0.826 & 0.414 & 0.714 & 0.833 \\
& & 4000 & 0.900 & 0.971 & 0.989 & 0.903 & 0.972 & 0.990 & 0.909 & 0.974 & 0.991 & 0.910 & 0.976 & 0.989 & 0.911 & 0.976 & 0.989 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.5} & 500 & 0.003 & 0.029 & 0.083 & 0.003 & 0.036 & 0.095 & 0.003 & 0.040 & 0.104 & 0.003 & 0.041 & 0.114 & 0.004 & 0.050 & 0.118 \\
& & 1000 & 0.034 & 0.250 & 0.439 & 0.039 & 0.260 & 0.468 & 0.047 & 0.275 & 0.472 & 0.061 & 0.297 & 0.478 & 0.067 & 0.316 & 0.490 \\
& & 2000 & 0.456 & 0.771 & 0.866 & 0.473 & 0.779 & 0.874 & 0.486 & 0.791 & 0.881 & 0.503 & 0.800 & 0.884 & 0.512 & 0.808 & 0.885 \\
& & 4000 & 0.939 & 0.987 & 0.994 & 0.941 & 0.990 & 0.994 & 0.945 & 0.991 & 0.995 & 0.946 & 0.990 & 0.996 & 0.947 & 0.989 & 0.996 \\
\midrule
\multirow{12}{*}{6} & \multirow{4}{*}{4} & 500 & 0.001 & 0.016 & 0.053 & 0.001 & 0.019 & 0.054 & 0.001 & 0.022 & 0.060 & 0.001 & 0.024 & 0.069 & 0.002 & 0.029 & 0.072 \\
& & 1000 & 0.016 & 0.105 & 0.184 & 0.018 & 0.110 & 0.197 & 0.020 & 0.117 & 0.214 & 0.030 & 0.118 & 0.224 & 0.036 & 0.132 & 0.237 \\
& & 2000 & 0.126 & 0.447 & 0.641 & 0.137 & 0.461 & 0.658 & 0.153 & 0.464 & 0.676 & 0.169 & 0.478 & 0.692 & 0.185 & 0.494 & 0.703 \\
& & 4000 & 0.693 & 0.918 & 0.958 & 0.715 & 0.925 & 0.962 & 0.725 & 0.932 & 0.968 & 0.737 & 0.938 & 0.969 & 0.744 & 0.940 & 0.971 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.25} & 500 & 0.001 & 0.021 & 0.069 & 0.002 & 0.024 & 0.076 & 0.002 & 0.026 & 0.081 & 0.003 & 0.032 & 0.087 & 0.003 & 0.036 & 0.093 \\
& & 1000 & 0.017 & 0.128 & 0.250 & 0.026 & 0.136 & 0.268 & 0.031 & 0.141 & 0.278 & 0.040 & 0.155 & 0.284 & 0.046 & 0.167 & 0.297 \\
& & 2000 & 0.170 & 0.554 & 0.738 & 0.193 & 0.566 & 0.754 & 0.214 & 0.578 & 0.767 & 0.228 & 0.597 & 0.775 & 0.237 & 0.612 & 0.785 \\
& & 4000 & 0.806 & 0.952 & 0.974 & 0.825 & 0.958 & 0.979 & 0.832 & 0.963 & 0.979 & 0.842 & 0.966 & 0.980 & 0.848 & 0.969 & 0.982 \\
\cmidrule{2-18}
& \multirow{4}{*}{4.5} & 500 & 0.002 & 0.028 & 0.082 & 0.003 & 0.033 & 0.084 & 0.003 & 0.037 & 0.092 & 0.003 & 0.043 & 0.099 & 0.005 & 0.053 & 0.113 \\
& & 1000 & 0.023 & 0.150 & 0.319 & 0.034 & 0.160 & 0.329 & 0.045 & 0.178 & 0.337 & 0.051 & 0.194 & 0.347 & 0.059 & 0.206 & 0.364 \\
& & 2000 & 0.240 & 0.651 & 0.825 & 0.261 & 0.669 & 0.836 & 0.268 & 0.682 & 0.848 & 0.290 & 0.696 & 0.855 & 0.297 & 0.710 & 0.857 \\
& & 4000 & 0.877 & 0.968 & 0.991 & 0.880 & 0.971 & 0.992 & 0.890 & 0.974 & 0.993 & 0.900 & 0.976 & 0.995 & 0.903 & 0.977 & 0.996 \\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}
\section{Empirical Applications} \label{sec:empirical}
In addition to the numerical evidence in Section \ref{sec:simulation}, we illustrate the performance of the proposed test in empirical applications.
The two empirical applications considered here, Islamic election and high-school attendance, represent canonical sharp and fuzzy RD designs in political and education economics.
Applying our testing procedure to these well-studied datasets allows us to assess whether treatment effect heterogeneity is driven by unobserved characteristics. In each application, we implement the test in two steps.
We first examine unconditional homogeneity at the cutoff and then test the conditional restriction after introducing the key covariates from the original study.
Comparing the two results reveals whether the observed heterogeneity is attributable to covariates or to unobserved factors.
Throughout this section, rejection and non-rejection are interpreted in terms of the identifying and regularity conditions that the test maintains.
\subsection{The treatment effect of the election of an Islamic party} \label{subsec:Meyersson}
In this subsection, we revisit \cite{meyersson2014}, which studies the consequences of the election of an Islamic party on high-school completion via an RD design.
\cite{meyersson2014} finds that the election of the Islamic party in Turkey causally increased female secular high school completion and reveals an underlying mechanism.
The paper further demonstrates treatment effect heterogeneity across subpopulations defined by proxy variables for poverty and Islamic conservatism and shows that the treatment effect is larger in more religiously conservative communities.
In this application, the running variable $R$ is the Islamic win margin---the vote share difference between the largest Islamic party and the largest secular party in the 1994 municipal elections---with a cutoff at zero; the binary treatment $D$ denotes assignment to an Islamic mayor.
The RD design in this case is sharp since $D=1$ when $R\ge0$. The primary outcome $Y$ is the share of women aged 15--20 who had completed a secular high school education in 2000.
A vector of covariates $X = (X_1, X_2, X_3)'$ includes the Islamic vote share in the 1994 election ($X_1$), the illiteracy rate in 2000 ($X_2$), and the proportion of religious buildings in 1990--2000 ($X_3$).
$X_1$ is a proxy for Islamic conservatism of the municipality; $X_2$ is a proxy variable for poverty; and $X_3$ measures the religious piety of the municipality.
We use a dataset containing election and census information for Turkish municipalities from \cite{meyersson2014}, comprising 2630 observations.
The proposed testing procedure is implemented with both local constant and local linear estimators.
The estimating bandwidth $h_r$ is selected based on \cite{Calonico2014} with an undersmoothing factor $n^{1/5-1/k_1}$, while the testing bandwidth $h = h_r \times n^{1/5-1/k_2}$ is undersmoothed.
The KS statistic is computed as the supremum over a grid of $L=1000$ or $2000$ points drawn from a normal distribution.
The number of bootstrap replications and the kernel function are the same as those in Section \ref{sec:simulation}.
The test results are presented in Table \ref{tab:Meyersson2014}.
\begin{table}[H]
\centering
\caption{The bootstrap $p$-values of the proposed test}
\label{tab:Meyersson2014}
\setlength{\extrarowheight}{-1pt}
\begin{adjustbox}{max width=\textwidth, max height=\textheight}
\begin{tabular}{ccccccccccccc}
\toprule & \multicolumn{6}{c}{Local Constant Estimator} & \multicolumn{6}{c}{Local Linear Estimator} \\
\cmidrule(lr){2-7} \cmidrule(lr){8-13}
& \multicolumn{3}{c}{$L=1000$} & \multicolumn{3}{c}{$L=2000$} & \multicolumn{3}{c}{$L=1000$} & \multicolumn{3}{c}{$L=2000$} \\
\cmidrule(lr){2-4} \cmidrule(lr){5-7} \cmidrule(lr){8-10} \cmidrule(lr){11-13}
$k_1 \backslash k_2$ & 4.25 & 4.50 & 4.75 & 4.25 & 4.50 & 4.75 & 4.25 & 4.50 & 4.75 & 4.25 & 4.50 & 4.75 \\
\midrule & \multicolumn{12}{l}{\textit{Panel A: Testing without conditioning on covariates}} \\
4.00 & 0.001 & 0.000 & 0.000 & 0.001 & 0.000 & 0.000 & 0.001 & 0.000 & 0.000 & 0.001 & 0.000 & 0.000 \\
4.25 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 \\
4.50 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 \\
\midrule & \multicolumn{12}{l}{\textit{Panel B: Testing conditional on $X_1$}} \\
4.00 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 \\
4.25 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 \\
4.50 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 \\
\midrule & \multicolumn{12}{l}{\textit{Panel C: Testing conditional on the covariate vector $X$}} \\
4.00 & 0.012 & 0.003 & 0.002 & 0.012 & 0.003 & 0.002 & 0.012 & 0.003 & 0.002 & 0.012 & 0.003 & 0.002 \\
4.25 & 0.002 & 0.001 & 0.001 & 0.002 & 0.001 & 0.001 & 0.002 & 0.001 & 0.001 & 0.002 & 0.001 & 0.001 \\
4.50 & 0.001 & 0.001 & 0.000 & 0.001 & 0.001 & 0.001 & 0.001 & 0.001 & 0.001 & 0.001 & 0.001 & 0.001 \\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}
Table \ref{tab:Meyersson2014} summarizes the test results.
At the 1\% significance level, Panel A rejects unconditional homogeneity, while Panel B rejects the conditional null given the Islamic vote share.
Thus, detectable treatment effect heterogeneity remains after conditioning on $X_1$.
This conclusion is stable across the reported tuning-parameter choices and complements the heterogeneity analysis in \cite{meyersson2014}.
Notably, the null is rejected even in Panel C after controlling for the full covariate vector, indicating that those covariates alone do not fully account for the heterogeneous response to the election of an Islamic party.
Because the conditioning sets differ in dimension across panels, the corresponding statistics are not on a common scale and cannot be used to rank how much of the heterogeneity each covariate set accounts for.
The results establish that unobserved heterogeneity remains under every conditioning set considered, not that the covariates contribute nothing.
Other characteristics appear to play a role in determining female secular high school completion rates, suggesting a more complex underlying mechanism.
We also examine specifications conditioning on each covariate individually (results available upon request): controlling for $X_2$ (illiteracy rate) or $X_3$ (religious buildings) alone yields qualitative conclusions similar to those in Panel B, with $p$-values near zero, confirming that no single observable dimension plausibly exhausts the treatment effect heterogeneity in this setting.
\subsection{The treatment effect of attending a better high school} \label{subsec:popeleches}
In this subsection, we revisit \cite{Pop-Eleches2013}, which studies Romanian administrative and survey data.
In this study, students can be admitted to high school if their transition score exceeds the admission cutoff.
\cite{Pop-Eleches2013} concludes that attending a more selective high school improves average academic performance on the Baccalaureate exam.
Our focus is on whether the treatment effect of attending a better school varies across students and whether it remains heterogeneous after accounting for peer quality.
In this case, the running variable $R$ is the distance between a student's score and the high school admission cutoff, and the binary treatment $D$ is attendance at the more selective high school with a higher cutoff in a given town.
The covariate $X$ is peer quality (average score of the student's class), and the outcome $Y$ is the student's Baccalaureate exam grade upon graduation.
This setting constitutes a fuzzy RD design because not all students whose scores exceed the admission cutoff choose to attend the more selective school.
Following \cite{Pop-Eleches2013}, we restrict the sample to towns with only two schools to avoid the complexity of multiple cutoffs, resulting in 23507 observations.
The test settings remain the same as those in Section \ref{subsec:Meyersson}, and results are reported in Table \ref{tab:popeleches2013}.
\begin{table}[H]
\centering
\caption{The bootstrap $p$-values of the proposed test}
\label{tab:popeleches2013}
\setlength{\extrarowheight}{-1pt}
\begin{adjustbox}{max width=\textwidth, max height=\textheight}
\begin{tabular}{ccccccccccccc}
\toprule & \multicolumn{6}{c}{Local Constant Estimator} & \multicolumn{6}{c}{Local Linear Estimator} \\
\cmidrule(lr){2-7} \cmidrule(lr){8-13}
& \multicolumn{3}{c}{$L=1000$} & \multicolumn{3}{c}{$L=2000$} & \multicolumn{3}{c}{$L=1000$} & \multicolumn{3}{c}{$L=2000$} \\
\cmidrule(lr){2-4} \cmidrule(lr){5-7} \cmidrule(lr){8-10} \cmidrule(lr){11-13}
$k_1 \backslash k_2$ & 4.25 & 4.50 & 4.75 & 4.25 & 4.50 & 4.75 & 4.25 & 4.50 & 4.75 & 4.25 & 4.50 & 4.75 \\
\midrule & \multicolumn{12}{l}{\textit{Panel A: Testing without conditioning on $X$}} \\
4.00 & 0.146 & 0.056 & 0.023 & 0.147 & 0.056 & 0.023 & 0.146 & 0.055 & 0.023 & 0.146 & 0.056 & 0.023 \\
4.25 & 0.052 & 0.016 & 0.003 & 0.052 & 0.016 & 0.003 & 0.052 & 0.016 & 0.003 & 0.052 & 0.016 & 0.003 \\
4.50 & 0.016 & 0.003 & 0.001 & 0.016 & 0.003 & 0.001 & 0.016 & 0.003 & 0.001 & 0.016 & 0.003 & 0.001 \\
\midrule & \multicolumn{12}{l}{\textit{Panel B: Testing conditional on $X$}} \\
4.00 & 0.186 & 0.239 & 0.220 & 0.282 & 0.347 & 0.318 & 0.943 & 0.932 & 0.960 & 0.944 & 0.934 & 0.956 \\
4.25 & 0.672 & 0.616 & 0.588 & 0.513 & 0.398 & 0.360 & 0.870 & 0.915 & 0.946 & 0.926 & 0.945 & 0.973 \\
4.50 & 0.497 & 0.416 & 0.421 & 0.481 & 0.378 & 0.344 & 0.984 & 0.990 & 0.992 & 0.976 & 0.985 & 0.988 \\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}
First, we test whether the treatment effect is homogeneous among all individuals in Panel A of Table \ref{tab:popeleches2013}.
The proposed test rejects the null of homogeneity at the 10\% significance level in most bandwidth configurations, indicating heterogeneity in individual treatment effects.
We then test whether such heterogeneity can be fully explained by the covariate $X$ in Panel B.
The local constant estimator yields $p$-values ranging from 0.186 to 0.672, so the conditional null is not rejected at the 10\% level, while the local linear estimator yields $p$-values from 0.870 to 0.992, a clear failure to reject at any conventional level.
This discrepancy between the two estimators warrants discussion. The local constant estimator is known to have a larger boundary bias at the cutoff in RD settings, which enters the first-stage estimation of $\tau(X)$ and propagates into the test statistic through the estimation effect characterized in Lemma~\ref{Lemma1}.
The local linear estimator, by contrast, achieves a smaller boundary bias, yielding a more accurate estimate of the CLATE.
The substantially larger $p$-values obtained under the local linear estimator are therefore more reliable for inference in this application.
Under either estimator, the heterogeneity detected in the unconditional analysis (Panel A) is no longer detected after conditioning on $X$ (Panel B).
Peer quality, therefore, accounts for the treatment effect heterogeneity that the unconditional test detects.
This finding indicates that peer quality is a relevant summary of detectable heterogeneity in this application and is informative for targeting.
Taken together, the two applications illustrate two qualitatively distinct conclusions that the proposed test can support.
In the Islamic-election application, significant unobserved heterogeneity persists even after conditioning on covariates, suggesting that latent characteristics are additional drivers of heterogeneous treatment response.
In the school-quality application, conditioning on peer quality accounts for the detectable heterogeneity in the unconditional test, consistent with peer quality serving as a summary of that heterogeneity.
These conclusions demonstrate the practical value of the test: it distinguishes settings in which observables leave residual heterogeneity from those in which they account for the heterogeneity.
We emphasize that the empirical results are interpreted as diagnostics under the maintained assumptions rather than as proof that a covariate set is sufficient.
\section{Conclusion}\label{sec:conclusion}
This paper proposes a novel testing method to detect unobserved heterogeneity in RD designs.
The proposed test transforms the null of no unobserved heterogeneity into the continuity of an imputed treated outcome at the cutoff.
An ICM restriction based on characteristic functions is then used to construct the test statistics.
We establish the limiting behavior of the proposed test under the null, the fixed alternative, and a sequence of local alternatives approaching the null at the rate $1/\sqrt{nh}$.
A straightforward multiplier bootstrap procedure provides feasible critical values that are theoretically valid and computationally efficient.
Numerical simulations and empirical applications show that the test controls size under the null and has increasing power against the alternatives considered.
\let\oldthebibliography\thebibliography
\renewenvironment{thebibliography}[1]{
\oldthebibliography{#1}
\linespread{1.1}\selectfont
}{
\endlist
}
\setlength{\bibsep}{3pt plus 1pt minus 1pt}
\putbib
\end{bibunit}
\newpage
\begin{bibunit}