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.
60,597 characters
The moment is here: a generalized class of estimators for fuzzy regression discontinuity designs
\title{\Large{The moment is here: a generalized class of estimators for fuzzy regression discontinuity designs}\thanks{I give special thanks to Sami Stouli for his extensive comments, suggestions, conversations, and guidance on this paper. I also thank Vincent Han, Arthur Lewbel, Adam McCloskey, Claudia Noack, Senay Sokullu, Frank Windmeijer, and participants in seminars at University College London for helpful comments, suggestions, and conversations. Code for estimation and inference in R and Python is available at \href{https://github.com/stuart-lane/LambdaFRD}{https://github.com/stuart-lane/LambdaFRD}.}}
\author{\large{Stuart Lane}\thanks{School of Economics, University of Bristol, UK, [email removed]}\\}
\date{}
\maketitle
\begin{abstract}
The standard fuzzy regression discontinuity (FRD) estimator is a ratio of differences of local polynomial estimators. I show that this estimator does not possess any finite integer moments, regardless of local polynomial degree, kernel function, or bandwidth. The estimator is heavy-tailed in small samples or when the treatment probability discontinuity at the cutoff is small. I present a generalized class of FRD estimators which preserves all finite moments from the data, indexed by a single tuning parameter, and nesting both standard FRD and sharp (SRD) estimators. Simple deterministic values of the tuning parameter lead to substantial improvements in median bias, median absolute deviation, and root mean squared error. Confidence intervals typically give reliable small-sample coverage in simulations. Estimator stability and performance are demonstrated using data on class size effects on educational attainment. \\
\\
\noindent \textbf{Keywords:} Fuzzy regression discontinuity, treatment effect estimation, instrumental variables, moments, finite sample, causal inference \\
\end{abstract}
\thispagestyle{empty}
\newpage
\section{Introduction}
Regression discontinuity (RD) designs are important tools for treatment effect estimation and causal inference when the probability of treatment assignment depends on an observed running variable, and is discontinuous at some known cutoff \citep{thistlethwaite1960regression, hahn2001identification}. If treatment assignment is deterministic conditional on the running variable, then the design is sharp (SRD), and if treatment is random conditional on the running variable, then the design is fuzzy (FRD). These methods have been applied to a wide range of topics, including the economics of education and educational outcomes \citep{angrist1999using}, industrial organization \citep{busse20061}, political economy \citep{lee2001electoral}, environmental economics \citep{salman2022paris}, the effects of social/financial aid programs \citep{ludwig2007does}, and medical applications \citep{cattaneo2023guide}.
The treatment effect of interest in the fuzzy design is identified by the ratio of the differences of conditional expectation functions \citep{hahn2001identification}. The standard estimator uses four local polynomial regressions, applied separately to the outcome and treatment conditional expectation functions, on either side of the cutoff. Here, the numerator and denominator are treated as separate SRD objects; this is justified asymptotically via the delta method, but ratio estimators potentially lack finite moments in finite samples if the denominator can be close to 0 with sufficient probability (e.g., the just-identified linear instrumental variables (IV) estimator under joint error normality \citep{chao2013expository}). Despite a significant theoretical literature discussing optimal choices regarding the degree of polynomial regression \citep{gelman2019high}, choice of kernel \citep{cheng1997automatic, imbens2008regression} and bandwidth selection \citep{imbens2012optimal, calonico2020optimal}, this particular issue of finite-sample instability has received limited attention; though \citet{noack2024bias} develop a bias-aware test statistic that bypasses the delta method and is robust to identification strength, but they do not consider estimation.
In this paper, I show that the standard FRD estimator does not have finite integer moments in finite samples, regardless of the degree of local polynomial regression, the kernel function, or the bandwidth used. This follows as the density of the denominator is bounded away from zero at the origin. Consequently, the sampling distribution has heavy tails, which can create significant problems for estimation and inference in small samples or with a small discontinuity in the treatment assignment probability at the cutoff. Heavy tails lead to substantial bias and large variability, and cause the finite-sample distribution to poorly approximate the asymptotic normal distribution, breaking down the delta method justification. This result provides an analogue to the lack of finite integer moments in finite samples for the ratio-IV estimator, as the local linear estimator with uniform kernel is well-known to be numerically equivalent to a ratio-IV estimator \citep{hahn2001identification, imbens2008regression}.
To fix the lack of finite moments for the FRD estimator, I introduce a class of estimators dependent on a single tuning parameter $\lambda$, nesting the standard estimator. I show that a continuum of estimators within this class preserve all finite moments within the data in finite samples; if the errors have finite moments of all orders (e.g., errors are normal, Laplacian etc.), then so will the new estimators. These estimators are asymptotically equivalent to the standard FRD estimator. Further, they are computationally simple, and are akin to a local nonparametric version of the $k$-class estimators from linear IV models \citep{nagar1959bias}. Simple choices for $\lambda$ dependent only on known model parameters (the effective sample size within the bandwidth and the degree of local polynomial regression) give strong performance in simulations; therefore, the researcher need not use selection algorithms such as cross-validation for $\lambda$.\footnote{While the degree of local polynomial regression $p$ is itself in principle a tuning parameter, it is typically treated as fixed at $p=1$ (with higher degrees sometimes employed as a robustness check) rather than as a data-driven tuning parameter (see e.g., \citet{gelman2019high}).} These estimators are similar to Fuller estimators in linear IV \citep{fuller1977some}. This class further demonstrates that the aforementioned equivalence of the local linear estimator with uniform kernel and the ratio-IV estimator is generalizable to any degree of local polynomial regression and choice of kernel function, thereby deepening the already well-established links between the linear IV and the FRD literature.
In simulations, the new estimators are shown to substantially improve on the standard FRD estimator in terms of median bias, median absolute deviation and root mean squared error, particularly with small sample sizes or a small jump in the treatment assignment probability (interpretable as a weaker instrument). The new estimators strictly dominate on root mean squared error, and have lower median bias and median absolute deviation for most parameter configurations. Confidence intervals exploiting the IV structure of the class typically have good small-sample coverage in simulations. While coverage is generally good in simulations, these confidence intervals are only pointwise asymptotically valid with undersmoothing; they do not have the MSE-optimal bandwidth pointwise asymptotic validity of bias-corrected confidence intervals \citep{calonico2014robust}, or the uniform asymptotic validity of bias-aware Anderson-Rubin confidence intervals \citep{noack2024bias}.
This improved finite-sample stability of the new class is demonstrated in an empirical application, looking at class size effects using the \citet{angrist1999using} dataset; point estimates and confidence intervals based on existing estimators show larger variation across different bandwidths, with some wide confidence intervals symptomatic of small sample sizes or potentially indicating weak identification. The point estimates and confidence intervals from the new class are stable across bandwidths, with both the smallest variation in point estimates and also the tightest confidence intervals in general. This stability across bandwidths and sample sizes is in line with the extensive simulation results of this paper. Results are also qualitatively the same for both verbal and mathematics test scores.
In Section \ref{sec model}, I introduce the model and assumptions, and show that the standard FRD estimator does not have finite integer moments in finite samples. In Section 3, I introduce a generalized class of FRD estimators which preserves all finite moments from the data in finite samples. Section \ref{sec lambda} provides details on selecting the tuning parameter, and shows the asymptotic equivalence of the standard and new estimators under minimal conditions. Section \ref{sec monte_carlo} gives a Monte Carlo study, and Section \ref{sec empirical} gives an empirical application with the \citet{angrist1999using} dataset looking at class size effects on educational performance. Proofs of theorems and supporting lemmas are presented in the Appendix. Further extensive simulation results are presented in the Online Appendix.
\section{Model and the moment problem}\label{sec model}
The baseline model is
\begin{align}
Y_i = m(X_i) + \tau D_i + U_i, \label{eq main model}
\end{align}
where $Y_i\in\mathbb{R}$ is the outcome of interest, $X_i\in \mathcal{X}\subseteq \mathbb{R}$ is the running variable, $m(\cdot): \mathcal{X} \to \mathbb{R}$ is some unknown continuous function, $\tau \in \mathcal{T}$ is the treatment effect of interest, where $\mathcal{T}\subset \mathbb{R}$ is compact, $D_i\in\{0,1\}$ is the treatment assignment, and $U_i\in\mathbb{R}$ is an unobserved error.
\begin{assum} \label{as fx}
(i) $\{(Y_i,X_i,D_i)\}_{i=1}^n$ is an i.i.d. sample. (ii) $X_i$ has density $f_X(x)$ for each $x\in \mathcal{X}$. (iii) $\{U_i\}_{i=1}^n$ are independent conditional on $(X_1,\hdots,X_n)$, and for each $i$, $U_i|X_i \sim f_{U|X}(\cdot|X_i)$ such that: (a) $\sup_{x \in \mathcal{X}}\mathbb{E}[|U_i|^r|X_i = x] < \infty$ for all $r < r^*$, for some $r^* \in (0,\infty]$, and $\sup_{x \in \mathcal{X}}\mathbb{E}[|U_i|^r|X_i = x] = \infty$ for all $r \geq r^*$; (b) for any compact $\mathcal{K} \subset \mathbb{R}$, then $\inf_{(u,x)\in \mathcal{K} \times \mathcal{X}}f_{U|X}(u|x) > 0.$
\end{assum}
Parts (i)-(ii) of Assumption \ref{as fx} are standard, assuming an i.i.d. sample with a continuous running variable. Errors have unbounded support with strictly positive density over compact subsets in $\mathbb{R}$; this is satisfied for common distributions such as the normal, $t$, and Laplace distributions. They may also possess finite moments of all orders e.g., normal, or no finite integer moments at all e.g., Cauchy. Define $\pi: \mathcal{X}\to [0,1]$ such that $\pi(x) := \mathbb{P}\{D_i = 1 | X_i = x\} = \mathbb{E}[D_i | X_i = x]$ for each $x\in\mathcal{X}$. The function $\pi(x)$ is continuous everywhere except for a discontinuity at a single point $x_0 \in \mathcal{X}$, referred to as the cutoff. Then, define
\begin{align*}
\mu_+(x_0) = \lim_{\nu\to 0}\mathbb{E}[Y_i|X_i = x_0 + \nu],&\,\,\, \mu_-(x_0) = \lim_{\nu\to 0}\mathbb{E}[Y_i|X_i = x_0 - \nu], \\
\pi_+(x_0) = \lim_{\nu\to 0}\mathbb{E}[D_i|X_i = x_0 + \nu],&\,\,\, \pi_-(x_0) = \lim_{\nu\to 0}\mathbb{E}[D_i|X_i = x_0 - \nu].
\end{align*}
\begin{assum}\label{as pi}
(i) $\mu_+(x_0)$, $\mu_-(x_0)$,
$\pi_+(x_0)$, and $\pi_-(x_0)$ exist. (ii) $\pi_+(x_0) \neq$ $ \pi_-(x_0)$.
\end{assum}
Assumption \ref{as pi} ensures that the conditional probability of assignment is discontinuous at $x_0$, and therefore the cutoff is well-defined. Given Assumptions \ref{as fx} and \ref{as pi}, \citet{hahn2001identification} show that $\tau$ is identified by the ratio
\begin{align}\label{eq tau}
\tau = \frac{\tau^Y}{\tau^D} = \frac{\mu_+(x_0) - \mu_-(x_0)}{\pi_+(x_0) - \pi_-(x_0)}.
\end{align}
\subsection{The standard estimator and no finite integer moments}\label{subsec standard est}
The standard approach to estimating \eqref{eq tau} uses local polynomial regressions to separately estimate each of $\mu_+(x_0)$, $\mu_-(x_0)$, $\pi_+(x_0)$, and $\pi_-(x_0)$. The four individual local polynomial estimators of degree $p\geq 0$ are
\begin{equation} \label{eq main estimators}
\begin{aligned}
\hat{\mu}_{p,+}(x_0) &= \frac{\displaystyle\sum_{i\in \mathcal{N}_+}Y_i K_{h}(X_i-x_0)\omega_{p,+,i}}{\displaystyle\sum_{i\in \mathcal{N}_+} K_{h}(X_i-x_0)\omega_{p,+,i}}, \,\,\,\, \hat{\mu}_{p,-}(x_0) = \frac{\displaystyle\sum_{i\in \mathcal{N}_-}Y_i K_{h}(X_i-x_0)\omega_{p,-,i}}{\displaystyle\sum_{i\in \mathcal{N}_-} K_{h}(X_i-x_0)\omega_{p,-,i}}, \\
\hat{\pi}_{p,+}(x_0) &= \frac{\displaystyle\sum_{i\in \mathcal{N}_+}D_i K_{h}(X_i-x_0)\omega_{p,+,i}}{\displaystyle\sum_{i\in \mathcal{N}_+} K_{h}(X_i-x_0)\omega_{p,+,i}}, \,\,\,\, \hat{\pi}_{p,-}(x_0) = \frac{\displaystyle\sum_{i\in \mathcal{N}_-}D_i K_{h}(X_i-x_0)\omega_{p,-,i}}{\displaystyle\sum_{i\in \mathcal{N}_-} K_{h}(X_i - x_0)\omega_{p,-,i}},
\end{aligned}
\end{equation}
where $K(\cdot)$ is some kernel function, $K_h(\upsilon) = K(\upsilon/h)/h$ for bandwidth $h$, $\mathcal{N}_+ = \{i: X_i \in \mathcal{X}_{h,+}\}$, $\mathcal{N}_- = \{i: X_i \in \mathcal{X}_{h,-} \}$, $\mathcal{X}_{h,-} = [x_0-h,x_0)$, $\mathcal{X}_{h,+} = [x_0,x_0+h]$, and where $\mathcal{X}_h = \mathcal{X}_{h,-} \cup \mathcal{X}_{h,+}$. Bandwidth $h > 0$ is taken as fixed.\footnote{Taking $h$ as fixed is natural for finite-sample analysis. The lack of moments is not a consequence of bandwidth choice, and cannot be remedied via appropriate bandwidth selection. In simulations, I provide clear evidence of the moment problem for $\hat{\tau}_1$ estimated with multiple data-driven bandwidths.} The generalized kernel weights in \eqref{eq main estimators} for $i \in \mathcal{N}_+$ are constructed as $\omega_{p,+,i} = e_1' \mathcal{S}^{-1}_{p,+} H_{p,i}$, where $e_1' = (1,0,\hdots,0)$ is a $(p+1)$-vector with first element 1 and 0 otherwise, $H_{p,i}$ is a $(p+1)$-vector with $j^{th}$ element $(X_i - x_0)^{j-1}$ and $\mathcal{S}_{p,+} = \sum_{i \in \mathcal{N}_+} K_h(X_i - x_0) H_{p,i} H_{p,i}'$. The analogous definitions for $\mathcal{S}_{p,-}$ and $\omega_{p,-,i}$ are readily adapted by exchanging $i \in \mathcal{N}_+$ for $i \in \mathcal{N}_-$. These four local polynomial estimators are then combined to form the estimator
\begin{align}\label{eq tau_hat}
\hat{\tau}_p = \frac{\hat{\tau}_p^Y}{\hat{\tau}_p^D} = \frac{\hat{\mu}_{p,+}(x_0) - \hat{\mu}_{p,-}(x_0)}{\hat{\pi}_{p,+}(x_0) - \hat{\pi}_{p,-}(x_0)}.
\end{align}
Define the effective sample to be the set of individuals $i \in \mathcal{N}_h$ for $\mathcal{N}_h = \mathcal{N}_+ \cup \mathcal{N}_-$, denote $n_+ = |\mathcal{N}_+|$, $n_- = |\mathcal{N}_-|$, and $n_h = n_+ + n_-$, and consider $Y_i$, $D_i$ and $H_{p,i}'$ for $i \in \mathcal{N}_h$. These individuals can be stacked into the $n_h\times 1$ vectors $Y$ and $D$, and $n_h\times(p+1)$ matrix $H_p$, and let $K := \textup{diag}(K_h(X_i - x_0))$ for $i \in \mathcal{N}_h$.\footnote{For notational clarity, $K(\cdot)$ denotes the scalar kernel function, $K_h(\cdot)$ denotes the scaled version, and $K:=\textup{diag}(K_h(X_i - x_0))$ denotes the associated diagonal matrix.} Then, \eqref{eq tau_hat} can also be expressed as
\begin{equation}\label{eq tau p matrix}
\hat{\tau}_p = \frac{e_1'(\mathcal{S}_{p,+}^{-1}H_{p,+}'K_+Y_+ - \mathcal{S}_{p,-}^{-1}H_{p,-}'K_-Y_-)}{e_1'(\mathcal{S}_{p,+}^{-1}H_{p,+}'K_+D_+ - \mathcal{S}_{p,-}^{-1}H_{p,-}'K_-D_-)},
\end{equation}
where for vector or matrix $A$ with $i^{th}$ row $A_i'$, $A_+$ denotes just the rows of $A$ relating to individuals $i \in \mathcal{N}_+$, and $A_-$ denotes just the rows relating to individuals $i \in \mathcal{N}_-$ \citep{fan1996local}. I treat $p \geq 0$ as fixed, as is common in practice.\footnote{Researchers typically fix $p = 1$ following \citet{gelman2019high} (sometimes using higher degrees as a robustness check), rather than using e.g., cross-validation to select $p$.}
\begin{assum} \label{as kernel}
(i) $K(\cdot)$ is a symmetric second-order $L_K$-Lipschitz kernel with support $[-1,1]$, strictly positive on $(-1,1)$. (ii) $K(\cdot) \in C_b^{\varsigma}((-1,1)\setminus \{0\})$ for some $\varsigma \geq 1$, for $C_b^s(\mathcal{A})$ the space of $s$-times continuously differentiable functions with bounded derivatives on $\mathcal{A}$. (iii) $K(\upsilon) \leq C_K$ for every $\upsilon\in[-1,1]$ for some $C_K < \infty$.
\end{assum}
All common kernels used for FRD estimation (the triangular, uniform and Epanechnikov kernels) satisfy Assumption \ref{as kernel}. Part (ii) requires $\varsigma$-times continuous differentiability with bounded derivatives for $\varsigma \geq 1$; while this is not usually explicitly assumed, it is satisfied for $\varsigma=\infty$ for all kernels used empirically, and so is a technical assumption rather than a meaningful restriction in practice. Differentiability of the kernel gives differentiability of $\hat{\tau}_p^D$ (see Lemma \ref{lem differentiability}), and is required as I show that configurations of $(X,D)$ exist that give exactly $\hat{\tau}_p^D = 0$, and then use geometric arguments to show that the density of $\hat{\tau}_p^D$ is bounded away from zero in an open neighborhood of $\hat{\tau}_p^D = 0$. This leads to the expectation integral diverging. The punctured interval in part (ii) allows for the triangular kernel $K(\upsilon) = (1 - |\upsilon|)1\{|\upsilon| < 1\}$, which is the most popular choice as it has optimal properties for boundary estimation \citep{cheng1997automatic, imbens2012optimal}, but is not differentiable at $\upsilon = 0$.
\begin{assum}\label{as full rank}
(i) For each $x \in \mathcal{X}_h$, $f_X(x)$ is bounded away from both 0 and $\infty$, and $\mathbb{E}[X_i^j]<\infty$ for $j\in\{1,\hdots,2p\}$.
(ii) For each $x \in \mathcal{X}_h$, $0 <\underline{\pi}\leq \pi(x) \leq \overline{\pi} < 1$.
(iii) There exist treated and untreated observations both above and below the cutoff within the effective sample.
(iv) The eigenvalues of $\mathcal{S}_{p,+}$ and $\mathcal{S}_{p,-}$
are bounded away from both 0 and $\infty$.
\end{assum}
Parts (i) and (ii) mildly strengthen standard regularity conditions to the bandwidth window $\mathcal{X}_h$ instead of an arbitrary open neighborhood of $x_0$ usually imposed for asymptotic results. This is necessary for finite-sample results. Part (iii) rules out perfect compliance within the observed effective sample, and excludes the pathological case of homogeneous treatment status; e.g., suppose $D_i = 0$ for every $i \in \mathcal{N}_h$ (which is theoretically compatible with the data generating process assumption in part (ii)). This trivially means $\hat{\tau}_p^D = 0$, hence $\hat{\tau}_p$ is undefined. Importantly, the divergence of moments for $\hat{\tau}_p$ persists even after imposing this mild structure on the observed effective sample, reinforcing the practical relevance of the results that follow. Part (iv) assumes that the sample kernel moment matrices are well-conditioned, implying $n_+,n_- \geq p+1$. Also note that $\sum_{i\in \mathcal{N}_+} K_{h}(X_i-x_0)\omega_{p,+,i} = \sum_{i\in \mathcal{N}_-} K_{h}(X_i-x_0)\omega_{p,-,i} = 1$ under Assumptions \ref{as kernel} and \ref{as full rank}, so the denominators can be dropped from the estimators in \eqref{eq main estimators}.
\begin{assum}\label{as gradient}
For fixed $D$, define $g : \mathcal{X}_h^{n_h} \to \mathbb{R}$ such that $g(X) = \hat{\tau}_p^D$. Then, $\norm{\nabla g(X)}_2 > 0$ for each $X\in (\mathcal{X}_h^\circ \setminus \{x_0\})^{n_h}$, where $\mathcal{A}^\circ$ denotes the interior of $\mathcal{A}$, $\nabla$ denotes the gradient of $g(\cdot)$, and $\norm{\cdot}_2$ is the Euclidean norm.
\end{assum}
Assumption \ref{as gradient} ensures that the gradient of $g(X)$ is strictly positive for $X\in (\mathcal{X}_h^\circ \setminus \{x_0\})^{n_h}$, necessary for the aforementioned geometric arguments to obtain finite-sample results. By the continuous distribution of $X$ and the polynomial nature of $\hat{\tau}_p^D$, either of the following conditions is individually sufficient for Assumption \ref{as gradient}: (i) $p\geq 1$, (ii) $|K'(\upsilon)|>0$ for each $\upsilon \in (-1,1)\setminus \{0\}$.\footnote{For the $p=0$ local constant estimator with uniform kernel, $\hat{\tau}_0^D = n_+^{-1}\sum_{i \in \mathcal{N}_+}D_i - n_-^{-1}\sum_{i \in \mathcal{N}_-}D_i$ is unconditionally discrete, and $\mathbb{P}\left\{n_+ = 2, n_- = 2, \sum_{i \in \mathcal{N}_+} D_i = 1, \sum_{i \in \mathcal{N}_-} D_i= 1\right\} > 0$ for $n \geq 4$ under the maintained assumptions. Therefore $\mathbb{P}\{\hat{\tau}_0^D = 0\} > 0$. Clearly $\mathbb{P}\{\hat{\tau}_0^Y \neq 0|X,D\} = 1$ since $Y_i$ has a continuous conditional distribution, then $\mathbb{P}\{|\hat{\tau}_0| = \infty\} > 0$, and so $\mathbb{E}[|\hat{\tau}_0|^r] = \infty$ for all $r > 0$.} The punctured interior $\mathcal{X}_h^\circ \setminus \{x_0\}$ is again used as Assumption \ref{as kernel} does not require differentiability of $K(\upsilon)$ at $\upsilon \in \{-1,0,1\}$. These assumptions lead to the first theorem; that $\hat{\tau}_p$ has no finite integer moments.
\begin{theorem} \label{thm no moments}
Let Assumptions \ref{as fx}-\ref{as gradient} hold. Then, $\mathbb{E}[|\hat{\tau}_{p}|^r] < \infty$ if and only if $r < \min\{r^*,1\}$, $\mathbb{E}[|\hat{\tau}_{p}|^r] = \infty$ otherwise.
\end{theorem}
Theorem \ref{thm no moments} holds regardless of the degree of local polynomial regression, kernel function, or bandwidth used. Even with normal errors with finite moments of all orders, $\hat{\tau}_p$ does not have a finite mean, because the density of the denominator $\hat{\tau}_p^D$ of $\hat{\tau}_p$ is bounded away from zero over some open neighborhood containing 0. This causes $\mathbb{E}[|\hat{\tau}_p^Y/\hat{\tau}_p^D|^r]$ to behave like $C^r\cdot\int_0^a x^{-r}dx$ for some $C \neq 0$, which will diverge if $r \geq 1$ for $a>0$. Due to this, the estimator may be highly dispersed in small samples, and potentially lead to inaccurate inferences. Simulations in Section \ref{sec monte_carlo} provide clear evidence of the moment problem for the local linear estimator $\hat{\tau}_{1}$ in small samples or with a small discontinuity in the treatment assignment probability.
\section{A generalized class of FRD estimators}\label{sec class}
It is well known \citep[see e.g.,][]{imbens2008regression} that when (\ref{eq tau_hat}) is estimated using $p = 1$ and the uniform kernel, then $\hat{\tau}_1$ is numerically equivalent to the IV estimator $\tilde{\tau}_{IV}$ with instrument $Z_i = 1\{X_i \geq x_0\}$ in the model
\begin{equation}\label{eq simple_reg}
Y_i = V_i'\delta + \tau D_i + U_i,
\end{equation}
where $V_i' = ( 1, \, Z_i(X_i - x_0), \, (1-Z_i)(X_i - x_0) )$ are included exogenous regressors and $\delta = (\delta_0, \, \delta_1, \, \delta_2)'$, and \eqref{eq simple_reg} is restricted to just the effective sample $i\in \mathcal{N}_h$ \citep[e.g.,][]{hahn2001identification}. This gives the IV estimator
\begin{align}\label{eq tau iv imbens}
\tilde{\tau}_{IV} = (Z'M_{V}D)^{-1}Z'M_{V}Y,
\end{align}
where $Z$ is the $n_h\times 1$ vector with $i^{th}$ element $Z_i$, $V$ is the $n_h\times 3$ matrix with $i^{th}$ row $V_i'$, and for some general $m \times q$ matrix $A$, then $M_A = I_m - P_A$ for $P_A = A(A'A)^{-1}A'$ and $I_m$ the $m \times m$ identity matrix. $\tilde{\tau}_{IV}$ is computationally simple and typically numerically similar to $\hat{\tau}_1$ with triangular kernel, and only requires one pass through the data instead of computing four separate local polynomial regressions. For this reason, $\tilde{\tau}_{IV}$ is also used in practice.
To my knowledge, the fact that the estimator in \eqref{eq tau iv imbens} can be generalized has not been explicitly noted. Numerous papers such as \citet{hahn2001identification}, \citet{imbens2008regression}, and \citet{noack2024bias} state that the equivalence holds for the local linear estimator with uniform kernel. However, the equivalence in fact holds for any kernel function and any degree of local polynomial regression. For general $p$, denote the weighted data by $\tilde{Y} = K^{1/2}Y$, $\tilde{D} = K^{1/2}D$, $\tilde{Z} = K^{1/2}Z$, $\tilde{U} = K^{1/2}U$, and $\tilde{V}_p = K^{1/2}V_p$, where $V_p$ has $i^{th}$ row $V_{p,i}' = (1, \, Z_i\tilde{H}_{p,i}', \, (1-Z_i)\tilde{H}_{p,i}')$ for $\tilde{H}_{p,i}$ the $p$-vector with $j^{th}$ element $(X_i - x_0)^j$. Therefore, $\tilde{Y}$, $\tilde{D}$, $\tilde{Z}$ and $\tilde{U}$ are $n_h$-vectors, and $\tilde{V}_p$ is an $n_h\times (2p+1)$ matrix. The weighted model then becomes
\begin{equation}\label{eq simple_reg p_degree}
\tilde{Y} = \tilde{V}_{p}\delta + \tau\tilde{D}+ \tilde{U},
\end{equation}
with $\delta\in \mathbb{R}^{2p+1}$, and the IV estimator of $\tau$ in \eqref{eq simple_reg p_degree} with instrument vector $\tilde{Z}$ is
\begin{align}\label{eq tau iv general kernel}
\hat{\tau}_{IV,p} = \left(\tilde{Z}'M_{\tilde{V}_p}\tilde{D}\right)^{-1} \tilde{Z}M_{\tilde{V}_p}\tilde{Y}.
\end{align}
\begin{prop} \label{prop equivalence iv}
Let Assumptions \ref{as fx}-\ref{as full rank} hold. Then $\hat{\tau}_{IV,p} \equiv \hat{\tau}_p$.
\end{prop}
This proposition shows that the well-known numerical equivalence between \eqref{eq tau_hat} and \eqref{eq tau iv imbens} when using a uniform kernel and local linear regression is generalizable to any kernel function and degree of local polynomial regression. Hence, standard FRD estimators can be implemented as a single weighted local IV estimator rather than the ratio of differences of four separate local polynomial regressions. Therefore, the general form in \eqref{eq tau iv general kernel} offers a computationally simple way of computing FRD estimators, and deepens the links between the FRD and linear IV literature. However, numerical equivalence implies that $\hat{\tau}_{IV,p}$ will also lack finite integer moments in finite samples.
\begin{corol}\label{cor tau_iv}
Let Assumptions \ref{as fx}-\ref{as gradient} hold. Then, $\mathbb{E}[|\hat{\tau}_{IV,p}|^r] < \infty$ if and only if $r < \min\{r^*,1\}$, $\mathbb{E}[|\hat{\tau}_{IV,p}|^r] = \infty$ otherwise.
\end{corol}
Corollary \ref{cor tau_iv} follows immediately by combining Theorem \ref{thm no moments} and Proposition \ref{prop equivalence iv}, and provides an analogue to the lack of finite moments of the just-identified ratio-IV estimator in linear IV models with joint normal errors (see e.g., \citet{chao2013expository}).
\subsection{The $\lambda$-class of generalized FRD estimators}
Consider a generalized class of estimators with parameter $\lambda \in [0,1]$ of the form
\begin{align}\label{eq lambda class estimators}
\hat{\tau}_{\lambda,p} = \left(\tilde{D}'M_{\tilde{V}_p}(I_{n_h} - \lambda M_{M_{\tilde{V}_p}\tilde{Z}})M_{\tilde{V}_p}\tilde{D}\right)^{-1} \tilde{D}'M_{\tilde{V}_p}(I_{n_h} - \lambda M_{M_{\tilde{V}_p}\tilde{Z}})M_{\tilde{V}_p}\tilde{Y}.
\end{align}
\eqref{eq lambda class estimators} has a structure akin to the $k$-class estimators from linear IV models \citep[e.g.,][]{nagar1959bias}, and the class nests $\hat{\tau}_p$ when setting $\lambda = 1$.
\begin{prop} \label{prop equivalence 2sls}
Let Assumptions \ref{as fx}-\ref{as full rank} hold. Then $\hat{\tau}_{1,p} \equiv \hat{\tau}_{p}$.
\end{prop}
The proof follows from least squares algebra. $\hat{\tau}_{\lambda,p}$ can be expressed in terms of the numerator $\hat{\tau}_p^Y$ and denominator $\hat{\tau}_p^D$ of $\hat{\tau}_p$ as
\begin{align}\label{eq tau lambda expanded}
\hat{\tau}_{\lambda,p} = \left(\lambda \tilde{\Gamma}_p(\hat{\tau}_p^D)^2 + (1-\lambda)\tilde{D}'M_{\tilde{V}_p}\tilde{D}\right)^{-1}\left(\lambda \tilde{\Gamma}_p\hat{\tau}_p^Y\hat{\tau}_p^D + (1-\lambda)\tilde{D}'M_{\tilde{V}_p}\tilde{Y}\right),
\end{align}
where $\tilde{\Gamma}_p$ is a scalar function of $X$, $K(\cdot)$, and $p$.\footnote{$\tilde{\Gamma}_p = \Gamma_{p,+}\Gamma_{p,-}/(\Gamma_{p,+}+\Gamma_{p,-})$, where $\Gamma_{p,+}$ and $\Gamma_{p,-}$ are the Schur complements corresponding to the lower $p\times p$ blocks in the $(p+1)\times (p+1)$ sample kernel moment matrices $\mathcal{S}_{p,+}$ and $\mathcal{S}_{p,-}$, respectively. See the proof of Proposition \ref{prop equivalence iv} in the Appendix for more details.} In \eqref{eq tau lambda expanded}, $\lambda$ is a mixing parameter, where $\lambda = 1$ yields the standard FRD estimator $\hat{\tau}_p = \hat{\tau}_p^Y/\hat{\tau}_p^D$ in \eqref{eq tau_hat}, and $\lambda = 0$ gives the ordinary least squares estimator of $\tau$ in \eqref{eq simple_reg p_degree}, which is the standard SRD estimator for \eqref{eq main model} when treatment assignment is deterministic conditional on $X_i$. The class $\hat{\tau}_{\lambda,p}$ for $\lambda\in[0,1]$ therefore defines a continuum of estimators mixing between the standard FRD and SRD estimators.
\begin{assum}\label{as deltasample}
There exists a set $\mathcal{N}_+(\delta,\kappa) \subseteq \mathcal{N}_{+}$ with $|\mathcal{N}_+(\delta,\kappa)| = 2p+1$, such that:
(i) $|X_i - X_j| \geq \delta$ for each $i\neq j$, $i,j \in \mathcal{N}_+(\delta,\kappa)$, for some $\delta > 0$;
(ii) $K_h(X_i - x_0) \geq \kappa$ for each $i \in \mathcal{N}_+(\delta, \kappa)$, for some $\kappa > 0$;
(iii) $D_i = 1$ and $D_j = 0$ for some $i, j \in \mathcal{N}_+(\delta, \kappa)$.
There exists some set $\mathcal{N}_-(\delta,\kappa) \subseteq \mathcal{N}_{-}$ with $|\mathcal{N}_-(\delta,\kappa)| = 2p+1$ with the equivalent properties.
\end{assum}
Assumption \ref{as deltasample} requires mild structure on the observed effective sample, stating that on each side of the cutoff, there are at least $2p+1$ individuals with distinct values $X_i$ with non-negligible kernel weights (i.e., not boundary points if the kernel is boundary-vanishing, such as the triangular or Epanechnikov kernels), and treatment variation among these individuals. Part (iii) mildly strengthens Assumption \ref{as full rank}(iii), as two individuals with heterogeneous treatment statuses must have distinct $X_i$ values. This minimal structure will hold in scenarios of practical relevance (for $p=1$, this requires just three individuals on each side of the cutoff with distinct $X_i$, non-negligible kernel weights, and varying treatments between two of the three), and is sufficient to preserve all finite moments in the data when $\lambda < 1$.
\begin{theorem} \label{thm lambda moments}
Let Assumptions \ref{as fx}-\ref{as gradient} and \ref{as deltasample} hold, and $ \lambda \in [0,1)$. Then,
$\mathbb{E}[|\hat{\tau}_{\lambda,p}|^r] < \infty$ for all $r < r^*$, $\mathbb{E}[|\hat{\tau}_{\lambda,p}|^r] = \infty$ otherwise.
\end{theorem}
The estimator $\hat{\tau}_{\lambda,p}$ preserves all finite moments in the data, and therefore will have finite moments of all orders with e.g., normal errors. Setting $\lambda < 1$ acts as a form of regularization that bounds the denominator away from zero. This is because the term $\tilde{D}'M_{\tilde{V}_p}\tilde{D}$ is bounded away from zero (see Lemma \ref{lem positivity}), and is added to the non-negative term $\lambda \tilde{\Gamma}_p (\hat{\tau}_p^D)^2$. This is similar to how Tikhonov regularization shifts the spectrum of a linear operator to bound eigenvalues away from zero, although unlike the Tikhonov parameter, $\lambda$ appears in both the numerator and denominator of \eqref{eq lambda class estimators}. With appropriate choices of $\lambda$, $\hat{\tau}_{\lambda,p}$ has excellent properties in the large-scale Monte Carlo study of Section \ref{sec monte_carlo} and the Online Appendix.
\section{Asymptotics and the selection of $\lambda$}\label{sec lambda}
The estimator $\hat{\tau}_{p}$ has well-studied asymptotic properties \citep[e.g.,][]{fan1996local,hahn2001identification,imbens2012optimal}. Only mild assumptions are needed to establish asymptotic equivalence between $\hat{\tau}_{\lambda,p}$ and $\hat{\tau}_{p}$.
\begin{assum}\label{as asymptotics}
(i) $h \to 0$, $nh \to \infty$, and $f_X(x_0)$ is continuous at $x_0$. (ii) $\lambda - 1 = o_p(1)$. (iii) $\lambda-1 = o_p(1/\sqrt{n_h})$, $nh^{2p+3} = O(1)$, and $m(\cdot),\pi(\cdot) \in C_b^{p+1}(\mathcal{B}_{\varepsilon}(x_0)\setminus\{x_0\})$ for some $\varepsilon > 0$. (iv) Assumption \ref{as fx}(iii)(a) is satisfied with $r^* \in (2,\infty]$.
\end{assum}
\begin{theorem} \label{thm lambda consistent}
Let Assumptions \ref{as fx}-\ref{as full rank} and \ref{as asymptotics}(i)-(ii) hold. Then, $\hat{\tau}_{\lambda,p} - \hat{\tau}_{p} \overset{p}{\to} 0$. If Assumptions \ref{as asymptotics}(iii)-(iv) also hold, then $\sqrt{n_h}(\hat{\tau}_{\lambda,p}-\tau) = \sqrt{2f_{X}(x_0) nh} \ (\hat{\tau}_p - \tau) + o_p(1)$.
\end{theorem}
Theorem \ref{thm lambda consistent} states that $\hat{\tau}_{\lambda,p}$ has the same probability limit and asymptotic distribution as $\hat{\tau}_{p}$ (up to a consistently estimable scaling factor, reflecting normalization by the observed effective sample $n_h$) under standard bandwidth and regularity conditions, given mild assumptions on $\lambda$. If $nh^{2p+3} = o(1)$, a simple test statistic $T(\hat{\tau}_{\lambda,p})$ for testing $H_0: \tau = \tau_0$ vs. $H_1: \tau \neq \tau_0$ is then
\begin{align*}
T(\hat{\tau}_{\lambda,p}) = \frac{\sqrt{n_h}(\hat{\tau}_{\lambda,p}-\tau_0)}{\sqrt{\hat{\mathbb{V}}(\hat{\tau}_{\lambda,p})}}, \,\,\,
\hat{\mathbb{V}}(\hat{\tau}_{\lambda,p}) = \frac{\tilde{D}'P_{M_{\tilde{V}_p}\tilde{Z}}\hat{\Omega}_{\tilde{U}} P_{M_{\tilde{V}_p}\tilde{Z}}\tilde{D}}{(\tilde{D}'M_{\tilde{V}_p}(I_{n_h} - \lambda M_{M_{\tilde{V}_p}\tilde{Z}})M_{\tilde{V}_p}\tilde{D})^2},
\end{align*}
where $\hat{\Omega}_{\tilde{U}}$ is some estimator of $\mathbb{E}[\tilde{U}_i^2(M_{\tilde{V}_p}\tilde{Z})_i(M_{\tilde{V}_p}\tilde{Z})_i'|X_i]$, and the variance estimator $\hat{\mathbb{V}}(\hat{\tau}_{\lambda,p})$ exploits the simple IV structure of \eqref{eq lambda class estimators}. Then, pointwise asymptotically valid confidence intervals with coverage probability $1-\alpha$ can be constructed as
\begin{align}\label{eq confidence interval}
\mathcal{C}(\hat{\tau}_{\lambda,p}) = \left[\hat{\tau}_{\lambda,p} - z_{1-\alpha/2}\sqrt{\frac{\hat{\mathbb{V}}(\hat{\tau}_{\lambda,p})}{n_h}}\,,\, \hat{\tau}_{\lambda,p} + z_{1-\alpha/2}\sqrt{\frac{\hat{\mathbb{V}}(\hat{\tau}_{\lambda,p})}{n_h}}\,\right],
\end{align}
where $n_{h,eff} := n_h - 2(p+1)$ is the residual degrees of freedom, and $z_{1-\alpha/2}$ denotes the $1-\alpha/2$ critical value of the $N(0,1)$ distribution. The simple confidence interval in \eqref{eq confidence interval} has in general good finite-sample coverage in simulations, competitive with alternatives such as the bias-corrected confidence intervals of \citet{calonico2014robust} and the bias-aware confidence intervals of \cite{noack2024bias}.\footnote{It is important to note however that the confidence intervals in \eqref{eq confidence interval} require undersmoothing, unlike the \citet{calonico2014robust} confidence intervals, and are not uniformly asymptotically valid, unlike the \cite{noack2024bias} confidence intervals.} Without undersmoothing, $\hat{\tau}_{\lambda,p}$ and $\hat{\tau}_{1,p}$ have the same asymptotic bias (up to scaling); constructing bias-corrected estimators and confidence intervals is an important topic for future study, but is beyond the scope of this paper.
\subsection{Choosing $\lambda$}
Selection of $\lambda$ is important for estimator performance. $\lambda = 1$ gives an estimator without finite integer moments in finite samples, and $\lambda = 0$ gives the SRD estimator consistent only in sharp designs. Define the function
\begin{align} \label{eq lambda function}
\Lambda(\psi) = 1 - \psi/n_{h,eff}
\end{align}
for $\psi \in [0, n_{h,eff}]$.\footnote{$\Lambda(\psi)$ is also a function of $n_h$ and $p$, but this is suppressed for notational convenience.} $\Lambda(\psi)$ is akin to the Fuller correction in linear IV models \citep{fuller1977some, hahn2004estimation}. For the linear-IV Fuller analogue, popular choices are $\psi = 1$ and $\psi = 4$. I provide strong simulation evidence that the estimators based on setting either $\psi = 1$ and $\psi = 4$ in \eqref{eq lambda function} perform very well; $\psi = 4$ gives the best performance across a wide variety of models in terms of median bias, median absolute deviation, and root mean squared error. Although $\psi = 4$ provides the best performance, $\psi = 1$ still provides large improvements relative to $\hat{\tau}_{1,p}$. I leave the study of $\lambda$-optimality (in the sense of asymptotic MSE minimization or some other desired property) to future research, with clear evidence in simulations that the simple choices detailed above are sufficient for immediate, substantial improvements in estimator performance.
\section{Monte Carlo simulations}\label{sec monte_carlo}
I consider the model $Y_i = m(X_i) + \tau D_i + U_i$ with
\begin{equation*}
m(X_i) =
\begin{cases}
0.48 + 1.27X_i + 7.18X_i^2 + 20.21X_i^3 + 21.54X_i^4 + 7.33X_i^5 & \textup{ if} \,\, X_i < 0, \\
0.48 + 0.84X_i - 3.00X_i^2 + 7.99X_i^3 - 9.01X_i^4 + 3.56X_i^5 & \textup{ if} \,\, X_i \geq 0,
\end{cases}
\end{equation*}
which is calibrated to \citet{lee2008randomized}. I set $X_i \sim N(0,1)$, $U_i\sim N(0,0.09)$, $x_0 = 0$ and $\tau = 0.04$. Three treatment assignment functions are used, given by
\begin{align*}
\pi_1(x) &= \pi_-1\{x < 0 \} + \pi_+1\{x \geq 0 \}, \\
\pi_2(x) &= (\pi_-x + \pi_-)1\{-1 \leq x <0 \} + (\pi_-x + \pi_+)1\{0 \leq x < 1 \} + 1\{x \geq 1 \}, \\
\pi_3(x) &= \pi_- \exp(0.2x) 1\{x < 0 \} + (\pi_+ + \pi_-(1-\exp(-0.2x)))1\{x \geq 0 \},
\end{align*}
where $\pi_- = 1-\pi_+$. I set $\pi_+\in\{0.6,0.7,0.8,0.9\}$, which gives treatment probability discontinuities of $(\pi_+-\pi_-)\in\{0.2,0.4,0.6,0.8\}$, and denote $\pi_+-\pi_-$ as $\pi_0$ for convenience. Identification becomes stronger as $\pi_0$ increases, and $\pi_0 = 1$ gives a sharp design. Each experiment uses $n\in \{300, 600\}$, and 10,000 repetitions. Further simulations with $X_i\sim 2Beta(2,4)-1$, $U_i \sim t(2.5)$ (which has finite variance, but heavy tails and no finite third moment), and $m(\cdot)$ calibrated to \citet{ludwig2007does} are reported in the Online Appendix. In total, there are 192 parameter configurations across the full $(n, \pi(\cdot), \pi_0, f_X(\cdot), f_U(\cdot), m(\cdot))$ grid.
\subsection{Estimator results}
\begin{table}[h!]
\centering
\addtolength{\tabcolsep}{-2pt}
\small
\begin{threeparttable}
\caption{Estimator results}
\vspace{-0.5em}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} c c c c c c c c c c c c c c}
\toprule
\multicolumn{2}{c}{$n=300$} & \multicolumn{4}{c}{Median bias} & \multicolumn{4}{c}{MAD} & \multicolumn{4}{c}{RMSE} \\
\cmidrule(r){1-2} \cmidrule(r){3-6} \cmidrule(l){7-10} \cmidrule(l){11-14}
$\pi_j$ & $\pi_0$ & $\hat{\tau}^{CCF}_{1}$ & $\hat{\tau}^{IK}_{1}$ & $\hat{\tau}^{CCF}_{\Lambda(4)}$ & $\hat{\tau}^{IK}_{\Lambda(4)}$ &
$\hat{\tau}^{CCF}_{1}$ & $\hat{\tau}^{IK}_{1}$ & $\hat{\tau}^{CCF}_{\Lambda(4)}$ & $\hat{\tau}^{IK}_{\Lambda(4)}$ &
$\hat{\tau}^{CCF}_{1}$ & $\hat{\tau}^{IK}_{1}$ & $\hat{\tau}^{CCF}_{\Lambda(4)}$ & $\hat{\tau}^{IK}_{\Lambda(4)}$ \\
\midrule
\multirow{4}{*}{1} & 0.2 & 0.09 & 0.12 & -0.01 & -0.00 & 0.50 & 0.49 & 0.09 & 0.09 & 89.22 & 126 & 0.14 & 0.14 \\
& 0.4 & 0.09 & 0.12 & 0.02 & 0.03 & 0.32 & 0.29 & 0.10 & 0.11 & 16.52 & 43.83 & 0.16 & 0.16 \\
& 0.6 & 0.08 & 0.09 & 0.04 & 0.05 & 0.22 & 0.19 & 0.12 & 0.11 & 12.74 & 36.02 & 0.18 & 0.17 \\
& 0.8 & 0.04 & 0.05 & 0.04 & 0.05 & 0.16 & 0.14 & 0.12 & 0.11 & 60.11 & 0.75 & 0.19 & 0.17 \\
\midrule
\multirow{4}{*}{2} & 0.2 & 0.09 & 0.14 & -0.01 & 0.00 & 0.50 & 0.48 & 0.09 & 0.09 & 143 & 41.39 & 0.15 & 0.15 \\
& 0.4 & 0.10 & 0.13 & 0.02 & 0.04 & 0.32 & 0.28 & 0.11 & 0.11 & 35.95 & 10.72 & 0.17 & 0.17 \\
& 0.6 & 0.07 & 0.08 & 0.04 & 0.05 & 0.22 & 0.19 & 0.12 & 0.12 & 79.85 & 6.23 & 0.19 & 0.18 \\
& 0.8 & 0.04 & 0.05 & 0.04 & 0.05 & 0.16 & 0.14 & 0.12 & 0.11 & 4.02 & 1.28 & 0.19 & 0.17 \\
\midrule
\multirow{4}{*}{3} & 0.2 & 0.07 & 0.12 & -0.01 & -0.00 & 0.52 & 0.50 & 0.09 & 0.09 & 837 & 76.48 & 0.14 & 0.14 \\
& 0.4 & 0.10 & 0.12 & 0.02 & 0.04 & 0.32 & 0.28 & 0.10 & 0.10 & 18.56 & 375 & 0.16 & 0.16 \\
& 0.6 & 0.07 & 0.09 & 0.04 & 0.06 & 0.22 & 0.19 & 0.12 & 0.12 & 17.79 & 3.02 & 0.18 & 0.18 \\
& 0.8 & 0.05 & 0.05 & 0.04 & 0.05 & 0.16 & 0.14 & 0.12 & 0.11 & 3.60 & 0.47 & 0.19 & 0.17 \\
\midrule
\midrule
\multicolumn{2}{c}{$n=600$} & \multicolumn{4}{c}{Median bias} & \multicolumn{4}{c}{MAD} & \multicolumn{4}{c}{RMSE} \\
\cmidrule(r){1-2} \cmidrule(r){3-6} \cmidrule(l){7-10} \cmidrule(l){11-14}
$\pi_j$ & $\pi_0$ & $\hat{\tau}^{CCF}_{1}$ & $\hat{\tau}^{IK}_{1}$ & $\hat{\tau}^{CCF}_{\Lambda(4)}$ & $\hat{\tau}^{IK}_{\Lambda(4)}$ &
$\hat{\tau}^{CCF}_{1}$ & $\hat{\tau}^{IK}_{1}$ & $\hat{\tau}^{CCF}_{\Lambda(4)}$ & $\hat{\tau}^{IK}_{\Lambda(4)}$ &
$\hat{\tau}^{CCF}_{1}$ & $\hat{\tau}^{IK}_{1}$ & $\hat{\tau}^{CCF}_{\Lambda(4)}$ & $\hat{\tau}^{IK}_{\Lambda(4)}$ \\
\midrule
\multirow{4}{*}{1} & 0.2 & 0.14 & 0.20 & 0.01 & 0.03 & 0.47 & 0.45 & 0.09 & 0.09 & 33.85 & 10.44 & 0.14 & 0.15 \\
& 0.4 & 0.11 & 0.14 & 0.04 & 0.07 & 0.25 & 0.22 & 0.10 & 0.11 & 52.17 & 14.69 & 0.16 & 0.17 \\
& 0.6 & 0.07 & 0.09 & 0.05 & 0.07 & 0.16 & 0.14 & 0.11 & 0.11 & 1.35 & 0.32 & 0.17 & 0.16 \\
& 0.8 & 0.04 & 0.05 & 0.04 & 0.06 & 0.12 & 0.10 & 0.10 & 0.09 & 2.25 & 0.16 & 0.15 & 0.14 \\
\midrule
\multirow{4}{*}{2} & 0.2 & 0.15 & 0.21 & 0.01 & 0.04 & 0.46 & 0.44 & 0.09 & 0.10 & 56.41 & 3954 & 0.14 & 0.16 \\
& 0.4 & 0.12 & 0.14 & 0.05 & 0.08 & 0.25 & 0.22 & 0.11 & 0.12 & 19.39 & 5.45 & 0.17 & 0.18 \\
& 0.6 & 0.07 & 0.08 & 0.06 & 0.07 & 0.16 & 0.14 & 0.11 & 0.11 & 0.70 & 0.45 & 0.17 & 0.17 \\
& 0.8 & 0.05 & 0.05 & 0.05 & 0.06 & 0.11 & 0.10 & 0.10 & 0.09 & 0.20 & 0.16 & 0.15 & 0.14 \\
\midrule
\multirow{4}{*}{3} & 0.2 & 0.13 & 0.19 & 0.01 & 0.03 & 0.46 & 0.43 & 0.09 & 0.09 & 44.08 & 70.92 & 0.14 & 0.15 \\
& 0.4 & 0.11 & 0.14 & 0.04 & 0.07 & 0.25 & 0.22 & 0.11 & 0.11 & 23.93 & 6.08 & 0.16 & 0.17 \\
& 0.6 & 0.07 & 0.09 & 0.05 & 0.07 & 0.16 & 0.14 & 0.11 & 0.11 & 2.31 & 0.68 & 0.17 & 0.16 \\
& 0.8 & 0.04 & 0.05 & 0.04 & 0.06 & 0.11 & 0.10 & 0.10 & 0.09 & 0.42 & 0.16 & 0.15 & 0.14 \\
\bottomrule
\end{tabular*}
\caption*{\citet{lee2008randomized} design, $X_i \sim N(0,1)$, $U_i \sim N(0,0.09)$. Each experiment is repeated 10,000 times. Any values greater than 100 are rounded to integers for space.}
\label{table estimation lee normal concise}
\end{threeparttable}
\end{table}
In Table \ref{table estimation lee normal concise}, I report the median bias, median absolute deviation (MAD), and root mean squared error (RMSE) of both $\hat{\tau}_{1,1}$ and $\hat{\tau}_{\lambda,1}$. Median bias and MAD are chosen as they are robust to heavy-tailed distributions. $\hat{\tau}_{1,1}$ is estimated using the triangular kernel, and $\hat{\tau}_{\lambda,1}$ is estimated using the uniform kernel, which appears to give superior performance for this class in simulations.\footnote{$\hat{\tau}_{\lambda,1}$ with triangular kernel still significantly improves on $\hat{\tau}_{1,1}$ in simulations.} As all estimators use local polynomial degree $p=1$, I suppress this subscript. For $\lambda$, I use the function $\Lambda(\psi)$ defined in \eqref{eq lambda function}, with $\psi = 1$ and $\psi = 4$, and therefore the three estimators under consideration are $\hat{\tau}_1$, $\hat{\tau}_{\Lambda(1)}$ and $\hat{\tau}_{\Lambda(4)}$. I use the MSE-optimal bandwidth of \citet{imbens2012optimal} and coverage-optimal bandwidth of \citet{calonico2020optimal}, denoted by $IK$ and $CCF$ respectively and indicated via superscripts. Each estimator is computed with both bandwidths. As $\hat{\tau}_{\Lambda(4)}$ gives better performance than $\hat{\tau}_{\Lambda(1)}$ (which itself still significantly improves on $\hat{\tau}_1$), I report only $\hat{\tau}_{\Lambda(4)}$ and $\hat{\tau}_1$ for conciseness in the main text, and present the full six-estimator results in the Online Appendix.
The top panel presents results for $n = 300$. With the largest discontinuity at $\pi_0 = 0.8$, median bias is essentially determined by bandwidth instead of $\lambda$. As the jump size decreases, the median bias for $\hat{\tau}_{\Lambda(4)}$ remains stable, but increases somewhat for $\hat{\tau}_1$, particularly with $CCF$ bandwidth. The new estimator $\hat{\tau}_{\Lambda(4)}$ also performs well for median absolute deviation, and remains essentially unchanged across the different $\pi(\cdot)$ functions and jump discontinuities. Even with location and scale measures robust to heavy tails, $\hat{\tau}_1$ exhibits an increasingly large median bias and MAD as $\pi_0$ decreases. $\hat{\tau}_{\Lambda(4)}$ controls location and scale much better than $\hat{\tau}_1$ across the range of sample sizes and probability discontinuity jumps considered.
The difference in performance is even clearer when considering RMSE; the new $\hat{\tau}_{\Lambda(4)}$ gives significant improvements in RMSE, even with a large discontinuity. With a small discontinuity, performance of $\hat{\tau}_{\Lambda(4)}$ remains stable, and even in small samples, the estimator retains a small RMSE. The lack of finite moments is clear here for $\hat{\tau}_1$, which displays some extremely large values. Even when $\pi_0 = 0.8$, $\hat{\tau}_1$ often attains a significantly larger RMSE than $\hat{\tau}_{\Lambda(4)}$. Overall, the new estimator demonstrates excellent performance in both sample sizes, and performance is good even for $n=300$. While $\hat{\tau}_1$ unsurprisingly does improve with sample size, the improvements are modest, and the moment problem is still evident at $n = 600$.
Table \ref{table estimator summary} summarizes the frequency for which each estimator (including $\hat{\tau}_{\Lambda(1)}$ with results reported in the Online Appendix) gives the best performance for each metric over the entire parameter grid (192 total combinations), and provides a strong justification for using $\hat{\tau}_{\Lambda(4)}$; $\hat{\tau}_{\Lambda(4)}$ strictly dominates in terms of RMSE, and gives the best median bias and median absolute deviation performance in 87.0\% and 96.4\% of parameter configurations respectively. Furthermore, $\hat{\tau}_{\Lambda(4)}$ performs well with both bandwidths; this is useful in practice, where researchers often consider multiple bandwidths to check the robustness of results. $\hat{\tau}_{\Lambda(1)}$ still improves on $\hat{\tau}_1$ for all three metrics in most data generating processes, but $\hat{\tau}_{\Lambda(4)}$ clearly provides the largest and most effective improvements.
\begin{table}[h]
\centering
\caption{Best performing estimator over 192 parameter configurations}
\begin{tabularx}{\textwidth}{c*{6}{>{\centering\arraybackslash}X}}
\hline
\toprule
Metric & $\hat{\tau}^{CCF}_{1}$ & $\hat{\tau}^{IK}_{1}$ & $\hat{\tau}^{CCF}_{\Lambda(1)}$ & $\hat{\tau}^{IK}_{\Lambda(1)}$ & $\hat{\tau}^{CCF}_{\Lambda(4)}$ & $\hat{\tau}^{IK}_{\Lambda(4)}$ \\
\midrule
Med. Bias & 12 & 6 & 7 & 0 & 119 & 48 \\
MAD & 0 & 7 & 0 & 0 & 82 & 103 \\
RMSE & 0 & 0 & 0 & 0 & 64 & 128 \\
\bottomrule
\label{table estimator summary}
\end{tabularx}
\vspace{-15pt}
\caption*{Summary of estimator performance by metric. For a given parameter configuration and metric, the best performing estimator is the one with the lowest absolute value of that metric. The new estimators are columns 3-6. There are 192 parameter configurations considered in total.}
\end{table}
\subsection{Inference}
Table \ref{table inference lee normal concise} reports coverage probabilities for a range of confidence intervals in the same models as Table \ref{table estimation lee normal concise}. $BC_1$ and $BC_2$ denote the bias-corrected confidence intervals from \citet{calonico2014robust}, using the $CCF$ and the $IK$ bandwidth selection algorithms, respectively. $\mathcal{C}_{\hat{\tau}_{\Lambda(4)}}^{CCF}$ and $\mathcal{C}_{\hat{\tau}_{\Lambda(4)}}^{IK}$ denote confidence intervals constructed using \eqref{eq confidence interval} formed with $\hat{\tau}^{CCF}_{\Lambda(4)}$ and $\hat{\tau}^{IK}_{\Lambda(4)}$, respectively. 95\% critical values are taken from the $t(n_{h,eff})$ distribution for constructing these intervals as a heuristic finite-sample correction. $AR_2$ denotes the bias-aware Anderson-Rubin confidence intervals of \citet{noack2024bias} using $ROT_2$ to estimate the smoothness bounds.\footnote{Results for the $AR$ test using $ROT_1$ are also presented in the Online Appendix.}
\begin{table}[h!]
\small
\centering
\begin{threeparttable}
\caption{Coverage rates of confidence intervals}
\vspace{-0.5em}
\begin{tabularx}{\textwidth}{cc *{5}{Y} @{\hspace{1.5em}} *{5}{Y}}
\toprule
& & \multicolumn{5}{c}{$n = 300$} & \multicolumn{5}{c}{$n = 600$} \\
\cmidrule(r){3-7} \cmidrule(l){8-12}
$\pi_j$ & $\pi_0$ & $BC_1$ & $BC_2$ & $\mathcal{C}_{\Lambda(4)}^{CCF}$ & $\mathcal{C}_{\Lambda(4)}^{IK}$ & $AR_2$ &
$BC_1$ & $BC_2$ & $\mathcal{C}_{\Lambda(4)}^{CCF}$ & $\mathcal{C}_{\Lambda(4)}^{IK}$ & $AR_2$ \\
\midrule
\multirow{4}{*}{1}
& 0.2 & 98.2 & 98.4 & 93.8 & 94.8 & 96.2 & 98.4 & 98.6 & 94.4 & 94.6 & 95.9 \\
& 0.4 & 97.0 & 97.0 & 95.3 & 95.7 & 96.3 & 96.6 & 96.6 & 96.2 & 95.5 & 95.6 \\
& 0.6 & 94.9 & 94.9 & 96.6 & 96.2 & 96.5 & 94.2 & 94.7 & 95.8 & 94.6 & 95.6 \\
& 0.8 & 92.8 & 92.7 & 96.3 & 95.8 & 96.3 & 93.2 & 93.4 & 95.8 & 94.1 & 96.0 \\
\midrule
\multirow{4}{*}{2}
& 0.2 & 98.2 & 98.2 & 93.4 & 93.8 & 96.1 & 98.4 & 98.8 & 94.6 & 94.3 & 95.4 \\
& 0.4 & 97.0 & 97.0 & 95.0 & 95.5 & 96.3 & 96.5 & 96.6 & 96.0 & 95.1 & 95.6 \\
& 0.6 & 94.6 & 94.5 & 96.3 & 96.2 & 95.9 & 94.7 & 94.8 & 95.9 & 94.3 & 95.6 \\
& 0.8 & 92.5 & 92.5 & 96.6 & 95.6 & 96.2 & 93.2 & 93.4 & 95.3 & 93.7 & 95.7 \\
\midrule
\multirow{4}{*}{3}
& 0.2 & 98.3 & 98.3 & 93.1 & 93.4 & 95.7 & 98.3 & 98.6 & 95.3 & 94.7 & 95.6 \\
& 0.4 & 96.8 & 96.6 & 95.4 & 95.7 & 96.2 & 97.0 & 96.9 & 96.2 & 94.9 & 95.6 \\
& 0.6 & 94.9 & 94.7 & 96.6 & 96.5 & 96.4 & 94.7 & 94.4 & 96.2 & 94.6 & 96.0 \\
& 0.8 & 92.3 & 92.4 & 96.4 & 96.1 & 95.9 & 93.1 & 93.3 & 95.9 & 94.2 & 96.0 \\
\bottomrule
\end{tabularx}
\caption*{\citet{lee2008randomized} design, $X_i \sim N(0,1)$, $U_i \sim N(0,0.09)$. Each experiment is repeated 10,000 times.}
\label{table inference lee normal concise}
\end{threeparttable}
\end{table}
The new confidence intervals based on $\hat{\tau}_{\Lambda(4)}$ appear to perform well, with coverage probabilities close to the nominal 95\%. $\mathcal{C}_{\hat{\tau}_{\Lambda(4)}}^{IK}$ is typically closer to $95\%$ coverage than $\mathcal{C}_{\hat{\tau}_{\Lambda(4)}}^{\textsc{CCF}}$, but both are mostly within simulation error of correct coverage. The $BC_1$ and $BC_2$ confidence intervals have good coverage as expected, although they are slightly conservative when $\pi_0 \in \{0.2, 0.4\}$. $AR_2$ has good coverage, especially for $n = 600$. Coverage is stable across discontinuity jumps, which is unsurprising, given that the $AR_2$ test is robust to the strength of identification. The additional tables in the Online Appendix give generally qualitatively similar results, although both the $\mathcal{C}$ and $AR$ confidence intervals can have large size distortions in the high-curvature \citet{ludwig2007does} design. It must be noted that the new confidence intervals are only pointwise asymptotically valid with undersmoothing, whereas the \citet{calonico2014robust} confidence intervals are pointwise valid with MSE-optimal bandwidths, and the \citet{noack2024bias} confidence intervals are uniformly valid; however, the $\lambda$-class estimators are designed to fix a fundamentally small-sample problem, and simulations suggest that $\mathcal{C}_{\hat{\tau}_{\Lambda(4)}}^{CCF}$ often gives good coverage in small samples. The $\mathcal{C}$ confidence intervals should therefore be seen as complementary to existing confidence intervals, and applied researchers may benefit from reporting $BC$ or $AR$ confidence intervals in conjunction with $\mathcal{C}$ when using $\lambda$-class estimators.
\section{Empirical application}\label{sec empirical}
I revisit the data of \citet{angrist1999using}, who estimate the effect of class sizes on test scores in Israel.\footnote{Data available at: \href{https://economics.mit.edu/people/faculty/josh-angrist/angrist-data-archive}{https://economics.mit.edu/people/faculty/josh-angrist/angrist-data-archive}.} Class sizes are determined using the Maimonides' rule, where class size should increase mechanically up to a maximum size of 40 students per class. When a $41^{st}$ student is enrolled, the class should split into two classes with average size 20.5. Class sizes should then mechanically increase until 80 students are enrolled. This process continues with a new cutoff at every multiple of 40. However, because compliance with the rule is not perfect, the setting is fuzzy.
\begin{sidewaystable}
\centering
\addtolength{\tabcolsep}{-2pt}
\small
\begin{threeparttable}
\caption{Class size effects on test scores}
\vspace{-0.5em}
\centering
\begin{tabular}{c c c c c c c c c c}
\hline
\multicolumn{10}{c}{(a) Verbal test scores} \\
\hline
\hline
$h$ & $n_h$ & $\hat{\tau}_{1}$ & $\hat{\tau}_{\Lambda(1)}$ & $\hat{\tau}_{\Lambda(4)}$ & $STD$ & $BC$ & $\mathcal{C}_{\Lambda(1)}$ & $\mathcal{C}_{\Lambda(4)}$ & $AR_2$ \rule{0pt}{2.3ex} \\
\hline
6 & 149 & -0.12 & -0.10 & -0.07 & [-0.39, 0.14] & [-0.79, 0.32] & [-0.23, 0.03] & [-0.15, 0.01] & $(-\infty, \infty)$ \\
8 & 229 & -0.10 & -0.09 & -0.08 & [-0.23, 0.02] & [-0.38, 0.09] & [-0.16, -0.01] & [-0.14, -0.02] & [-0.78, 0.10] \\
10 & 295 & -0.08 & -0.06 & -0.06 & [-0.16, -0.00] & [-0.27, 0.00] & [-0.11, -0.01] & [-0.10, -0.01] & [-0.26, -0.02] \\
12 & 379 & -0.07 & -0.05 & -0.05 & [-0.12, -0.01] & [-0.21, -0.02] & [-0.08, -0.01] & [-0.08, -0.01] & [-0.26, 0.06] \\
14 & 445 & -0.06 & -0.05 & -0.05 & [-0.10, -0.02] & [-0.17, -0.03] & [-0.08, -0.02] & [-0.08, -0.02] & [-0.14, 0.02] \\
16 & 527 & -0.05 & -0.03 & -0.03 & [-0.08, -0.01] & [-0.14, -0.03] & [-0.05, -0.01] & [-0.05, -0.01] & [-0.22, 0.10] \\
18 & 609 & -0.04 & -0.03 & -0.03 & [-0.07, -0.01] & [-0.12, -0.03] & [-0.05, -0.01] & [-0.05, -0.01] & [-0.10, 0.02] \\
\hline
\multicolumn{10}{c}{(b) Mathematics test scores} \\
\hline
\hline
$h$ & $n_h$ & $\hat{\tau}_{1}$ & $\hat{\tau}_{\Lambda(1)}$ & $\hat{\tau}_{\Lambda(4)}$ & $STD$ & $BC$ & $\mathcal{C}_{\Lambda(1)}$ & $\mathcal{C}_{\Lambda(4)}$ & $AR_2$ \rule{0pt}{2.3ex} \\
\hline
6 & 149 & -0.10 & -0.08 & -0.05 & [-0.32, 0.12] & [-0.62, 0.34] & [-0.20, 0.04] & [-0.12, 0.02] & $(-\infty, \infty)$ \\
8 & 229 & -0.09 & -0.07 & -0.06 & [-0.20, 0.02] & [-0.31, 0.10] & [-0.15, 0.00] & [-0.13, -0.00] & [-0.74, 0.14] \\
10 & 295 & -0.07 & -0.05 & -0.05 & [-0.14, 0.01] & [-0.22, 0.01] & [-0.10, 0.00] & [-0.09, 0.00] & [-0.22, -0.02] \\
12 & 379 & -0.05 & -0.03 & -0.03 & [-0.11, 0.00] & [-0.18, -0.02] & [-0.07, 0.01] & [-0.07, 0.01] & [-0.22, 0.06] \\
14 & 445 & -0.04 & -0.03 & -0.03 & [-0.09, 0.00] & [-0.15, -0.01] & [-0.07, 0.00] & [-0.07, 0.00] & [-0.14, 0.02] \\
16 & 527 & -0.03 & -0.02 & -0.02 & [-0.07, 0.00] & [-0.13, -0.02] & [-0.05, 0.01] & [-0.05, 0.01] & [-0.14, 0.06] \\
18 & 609 & -0.03 & -0.02 & -0.02 & [-0.06, 0.00] & [-0.11, -0.01] & [-0.04, 0.01] & [-0.04, 0.01] & [-0.10, 0.02] \\
\hline
\hline
\end{tabular}
\caption*{The MSE- and coverage optimal bandwidths are $h_{IK} = 9.75$ and $h_{CCF} = 6.66$ for verbal test scores, and $h_{IK} = 10.83$ and $h_{CCF} = 7.39$ for mathematics test scores}
\label{table empirical verbal}
\end{threeparttable}
\end{sidewaystable}
The outcome of interest is the average score for verbal and mathematics tests in $4^{th}$ grade classes, and the treatment of interest is the class size. The running variable is the total $4^{th}$ grade cohort enrollment for each school, and the number of children considered disadvantaged in each class is included as a control. I use the cutoff at 40 and consider $h\in\{6,8,10,12,14,16,18\}$. For reference, the MSE- and coverage optimal bandwidths of \cite{imbens2012optimal} and \cite{calonico2020optimal} are reported in the table notes. Results for the verbal and mathematics test scores are presented in Panels (a) and (b) of Table \ref{table empirical verbal} respectively.\footnote{While the theoretical framework developed in this paper assumes a binary treatment, the $\lambda$-class remains valid with non-binary treatments, as in this application. The target parameter of interest is identified by the ratio of left- and right-limit discontinuities of conditional mean functions $\mathbb{E}[Y|X]$ and $\mathbb{E}[D|X]$, and does not require $D$ to be binary. The support of class size here is also rich enough to approximate a continuous random variable. This dataset is widely used for empirical applications in papers studying FRD designs \citep[e.g.,][]{otsu2015empirical, feir2016weak, arai2022testing}.}
For estimation, I consider three estimators: the standard FRD estimator $\hat{\tau}_{1}$, and the new $\lambda$-class estimators $\hat{\tau}_{\Lambda(1)}$ and $\hat{\tau}_{\Lambda(4)}$, where the subscript for $p = 1$ is again dropped. I further report the standard and robust bias-corrected confidence intervals of \citet{calonico2014robust} for $\hat{\tau_1}$, denoted $STD$ and $BC$ respectively, the confidence intervals in \eqref{eq confidence interval} estimated with $\psi = 1$ and $\psi = 4$ with robust errors, and the $AR_2$ test of \citet{noack2024bias}. I also give the bandwidth $h$ and overall effective sample size used for each specification.
The point estimates from the $\lambda$-class estimators in Table \ref{table empirical verbal} are highly stable, particularly for $\hat{\tau}_{\Lambda(4)}$, with all estimates in [-0.08, -0.02]. These estimates also exhibit low variation across different bandwidths, which is of empirical relevance since researchers often report estimates with a range of bandwidths as a sensitivity check. The other estimators have greater variation in the point estimates, particularly the bias-corrected estimator. The confidence intervals $\mathcal{C}_{\Lambda(1)}$ and $\mathcal{C}_{\Lambda(4)}$ are also stable across specifications, and in general provide similar confidence regions. Potential weak identification is suggested for $h = 6$, as $AR_2$ is unbounded, and the $BC$ confidence intervals are wide. The new estimators give point estimates that remain stable across bandwidths relative to the other estimators, strengthening the robustness of the results and interpretations, and also produce tight confidence intervals; these findings align with both the theory and simulations of this paper.
\section{Conclusion}
In this paper, I show that the standard FRD estimator does not have finite integer moments in finite samples, leading to poor finite-sample performance. I present a new $\lambda$-class of estimators which preserves all moments in the data. This class is computationally simple and delivers significant improvements over existing estimators, particularly with small samples or a small treatment probability discontinuity. The choice $\lambda = \Lambda(4)$ provides especially good performance. I also show that simple confidence intervals typically have good coverage in small samples.
This work leads to many potential areas for future research. Of particular interest is the development of bias-corrected confidence intervals using the $\lambda$-class estimators. Bias-corrected confidence intervals use the square of the standard estimator denominator; as the $\lambda$-class and standard FRD estimators have the same asymptotic bias up to scaling, it may be possible to develop confidence intervals to use the more stable $\lambda$-class estimators, and in particular the denominator, for bias estimation. Another important avenue for future research is to consider theoretical $\lambda$-optimality results, or to develop finite-sample bandwidth selection algorithms that can specifically exploit the small-sample properties of the $\lambda$-class.
\bibliographystyle{apalike}
\bibliography{references}