Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
60,607 characters · 12 sections · 55 citation commands
The moment is here: a generalized class of estimators for fuzzy regression discontinuity designs
\thispagestyle{empty}
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 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 angrist1999using, industrial organization busse20061, political economy lee2001electoral, environmental economics salman2022paris, the effects of social/financial aid programs ludwig2007does, and medical applications cattaneo2023guide.
The treatment effect of interest in the fuzzy design is identified by the ratio of the differences of conditional expectation functions 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 chao2013expository). Despite a significant theoretical literature discussing optimal choices regarding the degree of polynomial regression gelman2019high, choice of kernel cheng1997automatic, imbens2008regression and bandwidth selection imbens2012optimal, calonico2020optimal, this particular issue of finite-sample instability has received limited attention; though 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 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 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., gelman2019high).} These estimators are similar to Fuller estimators in linear IV 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 calonico2014robust, or the uniform asymptotic validity of bias-aware Anderson-Rubin confidence intervals noack2024bias.
This improved finite-sample stability of the new class is demonstrated in an empirical application, looking at class size effects using the 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), 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) provides details on selecting the tuning parameter, and shows the asymptotic equivalence of the standard and new estimators under minimal conditions. Section (ref) gives a Monte Carlo study, and Section (ref) gives an empirical application with the 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.
The baseline model is
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.
Parts (i)-(ii) of Assumption (ref) 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
Assumption (ref) ensures that the conditional probability of assignment is discontinuous at $x_0$, and therefore the cutoff is well-defined. Given Assumptions (ref) and (ref), hahn2001identification show that $\tau$ is identified by the ratio
The standard approach to estimating (ref) 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
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 (ref) 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
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, (ref) can also be expressed as
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}_-$ fan1996local. I treat $p \geq 0$ as fixed, as is common in practice.\footnote{Researchers typically fix $p = 1$ following gelman2019high (sometimes using higher degrees as a robustness check), rather than using e.g., cross-validation to select $p$.}
All common kernels used for FRD estimation (the triangular, uniform and Epanechnikov kernels) satisfy Assumption (ref). 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)), 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 cheng1997automatic, imbens2012optimal, but is not differentiable at $\upsilon = 0$.
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) and (ref), so the denominators can be dropped from the estimators in (ref).
Assumption (ref) 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): (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) 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.
Theorem (ref) 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) 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.
It is well known imbens2008regression that when ((ref)) 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
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 (ref) is restricted to just the effective sample $i\in \mathcal{N}_h$ hahn2001identification. This gives the IV estimator
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 (ref) can be generalized has not been explicitly noted. Numerous papers such as hahn2001identification, imbens2008regression, and 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
with $\delta\in \mathbb{R}^{2p+1}$, and the IV estimator of $\tau$ in (ref) with instrument vector $\tilde{Z}$ is
This proposition shows that the well-known numerical equivalence between (ref) and (ref) 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 (ref) 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.
Corollary (ref) follows immediately by combining Theorem (ref) and Proposition (ref), 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., chao2013expository).
Consider a generalized class of estimators with parameter $\lambda \in [0,1]$ of the form
(ref) has a structure akin to the $k$-class estimators from linear IV models nagar1959bias, and the class nests $\hat{\tau}_p$ when setting $\lambda = 1$.
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
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) in the Appendix for more details.} In (ref), $\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 (ref), and $\lambda = 0$ gives the ordinary least squares estimator of $\tau$ in (ref), which is the standard SRD estimator for (ref) 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.
Assumption (ref) 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)(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$.
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)), 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 (ref). With appropriate choices of $\lambda$, $\hat{\tau}_{\lambda,p}$ has excellent properties in the large-scale Monte Carlo study of Section (ref) and the Online Appendix.
The estimator $\hat{\tau}_{p}$ has well-studied asymptotic properties fan1996local,hahn2001identification,imbens2012optimal. Only mild assumptions are needed to establish asymptotic equivalence between $\hat{\tau}_{\lambda,p}$ and $\hat{\tau}_{p}$.
Theorem (ref) 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
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 (ref). Then, pointwise asymptotically valid confidence intervals with coverage probability $1-\alpha$ can be constructed as
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 (ref) has in general good finite-sample coverage in simulations, competitive with alternatives such as the bias-corrected confidence intervals of calonico2014robust and the bias-aware confidence intervals of noack2024bias.\footnote{It is important to note however that the confidence intervals in (ref) require undersmoothing, unlike the calonico2014robust confidence intervals, and are not uniformly asymptotically valid, unlike the 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.
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
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 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 (ref) 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.
I consider the model $Y_i = m(X_i) + \tau D_i + U_i$ with
which is calibrated to 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
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 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.
In Table (ref), 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 (ref), 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 imbens2012optimal and coverage-optimal bandwidth of 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) 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.
Table (ref) reports coverage probabilities for a range of confidence intervals in the same models as Table (ref). $BC_1$ and $BC_2$ denote the bias-corrected confidence intervals from 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 (ref) 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 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.}
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 ludwig2007does design. It must be noted that the new confidence intervals are only pointwise asymptotically valid with undersmoothing, whereas the calonico2014robust confidence intervals are pointwise valid with MSE-optimal bandwidths, and the 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.
I revisit the data of 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.
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 imbens2012optimal and 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) 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 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 calonico2014robust for $\hat{\tau_1}$, denoted $STD$ and $BC$ respectively, the confidence intervals in (ref) estimated with $\psi = 1$ and $\psi = 4$ with robust errors, and the $AR_2$ test of 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) 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.
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.