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.
89,194 characters · 15 sections · 41 citation commands
Optimal estimation for regression discontinuity design with binary outcomes
A large-sample approximation is the basis for the leading estimators for regression discontinuity (RD) designs Imbens.Kalyanaraman2012,Calonico.Cattaneo.Titiunik2014. RD designs involve the estimation of conditional expectation functions at a cutoff point on the support of a running variable. Hence, effective observations are limited to the neighborhood of the cutoff, and the number of these observations can be small, even if the total sample size is large Cattaneo.Frandsen.Titiunik2015,Canay.Kamat2017. For example, an effective sample can be small for designs with multiple cutoffs, with a cutoff at the tail of the distribution, or with subgroup analyses. In small samples, the large-sample asymptotics may not provide good approximations of the behaviors of the existing estimators; hence, their desirable properties may be lost.
Some studies have considered finite-sample minimax estimators for RD designs.\footnote{Throughout the manuscript, we compare our estimator with existing finite-sample minimax estimators. Another notable approach is finite-sample valid estimation and inference based on the local randomization of the RD design Cattaneo.Frandsen.Titiunik2015, cattaneoInferenceRegressionDiscontinuity2016, cattaneoComparingInferenceApproaches2017. The local randomization approach is based on the assumption that the running variable is randomly assigned with a constant regression function within a given small window around the threshold Cattaneo_Idrobo_Titiunik_2024, whereas we consider a smooth but nonconstant regression function within the window.} For example, Armstrong.Kolesar2018 and Imbens.Wager2019 propose finite-sample minimax linear estimators under the smoothness of the regression function. However, these minimax estimators require knowledge of the conditional variance function, which is generally unavailable in practice. Although the variance can be estimated, we cannot guarantee the theoretical validity of the plug-in estimators with the estimated variance in finite samples. Furthermore, the construction of finite-sample valid confidence intervals based on these estimators additionally requires the normality of the regression errors.
In this study, we propose finite-sample estimation and inference methods for RD designs with binary outcomes. For a binary dependent variable, all features of its conditional distribution, including its conditional variance, are known functions of its conditional mean function. We establish the finite-sample validity of our methods under a smoothness restriction on the conditional mean function by considering the implicit restrictions imposed on the entire conditional distribution. Therefore, our procedure is both feasible and theoretically valid without knowledge or estimation of the conditional variance or, more generally, any features of the conditional distribution, except for the smoothness of the conditional mean.
Specifically, we consider a minimax optimal estimator among a class of {\it linear shrinkage} estimators for the regression function at a boundary point, under the assumption that the regression function satisfies Lipschitz continuity. The class of linear shrinkage estimators has the form $\sum_{i=1}^nw_i(Y_i-1/2)+1/2$ with $\sum_{i=1}^nw_i\le 1$ and $w_i\ge 0$, where $Y_1,...,Y_n$ are the observed outcomes on either side of the boundary. The shrinkage toward $1/2$ is motivated by the fact that the regression function is bounded and takes values in $[0,1]$, leading to a scope of efficiency gain by shrinkage. Given this class of linear shrinkage estimators, we derive a linear shrinkage estimator that minimizes the maximum mean squared error (MSE) under Lipschitz continuity with a known Lipschitz constant. In other words, we assume the researcher's a priori knowledge of the bound of how much the function value can change if the running variable is changed by one unit. We emphasize that the Lipschitz constant is the only tuning parameter used. Furthermore, we show that the minimax estimator is a solution to a convex optimization problem that is computationally feasible. Thus, we provide a practical estimator that achieves finite-sample optimality in the binary-outcome setting.
Our estimator is applicable to many practical RD designs. Binary outcomes are among the most common types of outcomes in empirical applications. For example, the following outcome variables are all binary: an indicator for winning the next election in the famous U.S. House election study by Lee2008; a corruption indicator in Brollo.Nannicini.Perotti.Tabellini2013; a mortality indicator in card_does_2009; and indicators for student's enrollment and dropout in Melguizo2016 and Cattaneo2021multi. Furthermore, the first stage in fuzzy RD designs often involves treatment status as a binary dependent outcome. Moreover, the minimax optimality of our estimator for binary outcomes immediately extends to that for bounded outcomes because the variance of any linear estimator is maximized when the outcomes are Bernoulli, given the conditional mean function. Hence, our estimator can be applied not only to binary outcomes but also to bounded outcomes, which are also frequently used in the RD design.
Our method also complements existing minimax estimators. We compare our estimator to a version of the existing minimax estimators Armstrong.Kolesar2018,Imbens.Wager2019 and demonstrate that our method has better finite-sample performance than the existing approach, although their asymptotic behaviors are similar. Specifically, we consider a minimax linear estimator obtained under a misspecified model in which the conditional mean and variance are unrelated, the variance is known, and the regression function lies in a Lipschitz class with no bounds on the function values. This estimator is not directly feasible in our binary-outcome setting, in which the variance is unknown. As a feasible version of this estimator, we consider one obtained under the assumption of a constant variance of $1/4$, which is the maximum possible variance of a binary variable. For binary outcomes, we theoretically show that the efficiency gain from our estimator, relative to the alternative estimator above, tends to vanish as the sample size increases. Nevertheless, for small samples, we numerically demonstrate that the alternative method can result in a $5\%$ to $20\%$ increase in the worst-case root MSE owing to model misspecification. Hence, our method supplements the existing minimax estimators with better finite-sample performance and similar asymptotic behaviors in a binary-outcome setting.
We also propose confidence intervals with correct coverage in finite samples uniformly over the Lipschitz class. We construct confidence intervals by inverting one- or two-sided uniformly valid tests that use a linear estimator as the test statistic. To construct a uniformly valid test, we propose a simulation-based approximation to the distribution of the test statistic by drawing samples from a multivariate Bernoulli distribution that satisfies the null restriction. We then numerically optimize the critical value so that the worst-case rejection probability is equal to or smaller than the significance level. A computational challenge with this approach is the calculation of the worst-case rejection probability, which involves the optimization of an $(n+1)$-dimensional parameter. We overcome this challenge by deriving a simple characterization of the worst-case rejection probability under Lipschitz continuity, which significantly reduces the computational burden. We also emphasize that our confidence intervals are valid in finite samples for binary outcomes. This contrasts with existing inference methods, which are based on either a large-sample approximation or a restrictive assumption of Gaussian errors with a known variance.
The same inference approach does not apply to bounded outcomes because the simple characterization of the worst-case rejection probability relies on the fact that the outcome is binary. For bounded outcomes, we provide an alternative finite-sample inference procedure based on a uniform bound on the rejection probability obtained using Hoeffding's inequality. The resulting confidence intervals have correct coverage in finite samples but can be conservative, similar to Hoeffding's inequality-based confidence intervals in other contexts.
We demonstrate the performance of our methods through simulations and an empirical application. In the simulations, our estimator achieves substantially smaller MSEs relative to the leading large-sample estimators when the sample size is small. Furthermore, our estimator behaves similarly to the large-sample estimators when the sample size is large; the differences in MSE decrease as the number of observations increases. Our proposed inference method also achieves guaranteed coverage rates with shorter confidence intervals when the sample size is small. Hence, our methods, while theoretically valid, are also useful in practice.
We illustrate our methods by revisiting Brollo.Nannicini.Perotti.Tabellini2013, who estimate the impact of additional government revenue on corruption. They exploit the regional fiscal rule in Brazil, where federal transfers to municipal governments change exogenously at given population thresholds. This is a multi-cutoff RD design with a small sample size near each cutoff point. We demonstrate that our estimates are similar to the conventional estimates for large-sample pooling of multiple cutoffs. Nevertheless, our inference method provides much shorter confidence intervals than conventional methods when we focus on a small sample near each cutoff value. Consequently, our estimates provide more informative results than conventional methods.
Both the simulation and application results indicate that the small-sample estimations are generally challenging, whereas our estimator has the potential to provide informative estimates. Hence, our estimator is a practical last resort for an empirical researcher facing a research question with a small effective sample size for an RD design.
In addition to contributing to estimation in RD design, we contribute to the vast literature on minimax estimation. donoho1994 considers minimax affine estimation and inference on linear functionals in nonparametric regression models with Gaussian errors. Recently, this framework has been applied to the estimation and inference of treatment effects in various settings, including RD designs Armstrong.Kolesar2018,Armstrong2021ATE,Gao2018,Imbens.Wager2019,Kwon.Kwon2020,Chaisemartin2021,rambachan2023parallel. We complement these studies by examining nonparametric regression models with Bernoulli dependent variables, which are not covered by their frameworks. To the best of our knowledge, no general minimax estimator under the squared error loss has been established for the problem of estimating linear functionals in this setting.\footnote{DeRouen1974 derive a $\Gamma$-minimax estimator for a linear combination of the success probabilities of multiple independent binomial variables when the class of prior distributions consists of distributions with the same, known means.} No solution is known, even for the estimation of the difference in the success probability between two independent binomial variables with unequal numbers of trials Lehmann1998.\footnote{For the estimation of the success probability of a single binomial variable, a linear shrinkage (toward $1/2$) estimator is minimax among all estimators Lehmann1998. Marchand2000 consider this problem with a restricted parameter space. They show that, when the success probability is known to lie in a symmetric interval around $1/2$, a linear shrinkage estimator is minimax among all linear estimators.} We contribute to this underexplored literature by developing a minimax estimator for a regression function at a point within the class of linear shrinkage estimators under the Lipschitz continuity of the regression function.
RD designs exploit a discontinuous change in treatment status when a running variable exceeds a cutoff point. For example, Brollo.Nannicini.Perotti.Tabellini2013 exploit discontinuous increases in the amount of central government subsidies for a local government when its population equals or exceeds a threshold level. The target parameter of the RD design is the average treatment effect at the cutoff point, which is identified as the difference between the conditional expectation functions evaluated at the cutoff point. Hence, its estimation involves the nonparametric estimation of the conditional mean functions at their boundary points.
Suppose that we have a random sample $\{Y_i,D_i,R_i\}_{i=1}^N$, where $R_i \in \mathbb{R}^{d_r}$ is a $d_r (\geq 1)$-dimensional vector of running variables, $Y_i$ is a binary outcome, $D_i$ is a binary treatment assigned as $D_i = 1\{R_i \in \mathcal{T}\}$, and $\mathcal{T} \subset \mathbb{R}^{d_r}$ is a known treated region. The leading case is that in which $R_i$ is univariate ($d_r = 1$) and $\mathcal{T} = [c,\infty)$ for some known cutoff $c$; however, the following arguments also apply to a multidimensional case (i.e., $d_r > 1$). Suppose
for some unknown function $f:\{0,1\}\times\mathbb{R}^{d_r}\rightarrow[0,1]$. Let $R_0$ be a fixed boundary point in the treatment region $\mathcal{T}$. When $f(d,r)$ represents the conditional expectation function of the underlying potential outcome $Y_{i}(d)$ conditional on $R_i = r$ for each $d \in \{0,1\}$, $f(1,R_0) - f(0,R_0)$ is interpreted as the average treatment effect at the boundary point $R_0$ hahnIdentificationEstimationTreatment2001. The data $\{Y_i,D_i,R_i\}_{i=1}^N$ can be divided into $\{ Y_{i,+}, R_{i,+} \}_{i=1}^{n_{+}}$ and $\{ Y_{i,-}, R_{i,-} \}_{i=1}^{n_{-}}$, where the former is the data from the treatment group and the latter is the data from the control group. We use the two samples separately to estimate $f(1,R_0)$ and $f(0,R_0)$, respectively.
Without loss of generality, we consider the estimation of $f(1,R_0)$ throughout this section, except in Remark (ref) at the end of this section, where we discuss the estimation of $f(1,R_0) - f(0,R_0)$. To simplify the notation, we use $\{Y_i,R_i\}_{i=1}^n$ to denote $\{ Y_{i,+}, R_{i,+} \}_{i=1}^{n_{+}}$, so that $R_i\in\mathcal{T}$ for all $i=1,...n$. Furthermore, we use $f(\cdot)$ to denote $f(1,\cdot)$. Additionally, our analysis conditions on the realization of $\{R_i\}_{i=1}^n$, and we treat $\{R_i\}_{i=1}^n$ as deterministic, so that $P(Y_i=1)=f(R_i)$ for all $i=1,\ldots,n$. Let $p_i \equiv f(R_i)$ for $i=0, 1, \ldots, n$ and $\bm{p} \equiv (p_0,p_1, \ldots, p_n)' \in [0,1]^{n+1}$. Without loss of generality, we assume that $R_0=0$ and $\|R_0\| \leq \|R_1\| \leq \cdots \leq \|R_n\|$, where $\|\cdot\|$ is a norm on $\mathbb{R}^{d_r}$. The following theoretical result holds for any norm. We focus on the Euclidean norm in our numerical exercises, simulations, and empirical application.
For the parameter of interest $p_0 = f(R_0) = f(0)$, we consider the following linear shrinkage estimator:
where $\mathcal{W} \equiv \left\{ \bm{w} \in \mathbb{R}^n : \sum_{i=1}^n w_i \leq 1 \ \text{and} \ w_i \geq 0 \ \text{for all $i$} \right\}$. When $\sum_{i=1}^n w_i =1$, $\hat{p}_0(\bm{w})=\sum_{i=1}^n w_iY_i$, and no shrinkage occurs. When $\sum_{i=1}^n w_i < 1$, $\hat{p}_0(\bm{w})$ is an estimator that shrinks toward $1/2$.
We assume that $f$ belongs to a Lipschitz class:
where $C$ denotes the Lipschitz constant and is known. This assumption implies that $\bm{p} \in [0,1]^{n+1}$ satisfies $|p_i-p_j| \leq C \|R_i - R_j\|$ for all $i$ and $j$. Conversely, if $|p_i-p_j| \leq C \|R_i - R_j\|$ for all $i$ and $j$, we can find a function $f\in\mathcal{F}_{\text{Lip}}(C)$ such that $f(R_i)=p_i$ for all $i$ Beliakov2006. Hence, the parameter space of $\bm{p}$ can be expressed as follows:
Since $Y_1,...,Y_n$ are independent binary variables, the mean squared error (MSE) of $\hat{p}_0(\bm{w})$ is given by
We consider the linear shrinkage estimator whose corresponding weight vector solves the following problem:
To simplify the expression in (ref), we redefine $p_i$ as $\theta_i \equiv p_i-1/2$ for $i=0, 1, \ldots, n$ and let $\bm{\theta} \equiv (\theta_0, \theta_1, \ldots , \theta_n)'$; thus, the problem is
where $\Theta \equiv \{\bm{\theta} \in [-1/2,1/2]^{n+1} : |\theta_i - \theta_j| \leq C \| R_i-R_j \| \ \text{for all $i$ and $j$} \}$ and
Hence, we aim to obtain a weight vector that minimizes the maximum MSE by solving (ref).
Our goal is to obtain a weight vector $\bm{w}$ that minimizes the maximum MSE. First, we consider the maximization part of ((ref)) for a given weight vector $\bm{w}\in \mathcal{W}$. Note that the objective function $\text{MSE}(\bm{w},\bm{\theta})$ is generally nonconcave in $\bm{\theta} = (\theta_0, \ldots, \theta_n)'$, as the squared bias $\left( \sum_{i=1}^n w_i \theta_i - \theta_0 \right)^2$ is convex in $\bm{\theta}$, whereas the variance $\sum_{i=1}^n w_i^2 \left( \frac{1}{4} - \theta_i^2 \right)$ is concave in $\bm{\theta}$. Nevertheless, we show that this nonconvex ($n+1$)-dimensional optimization problem can be reduced to an optimization problem over a single parameter $\theta_0$, and is therefore computationally tractable.
Note first that $\Theta$ is centrosymmetric (i.e., $\bm{\theta}\in \Theta$ implies $-\bm{\theta}\in\Theta$) and that $\text{MSE}(\bm{w},\bm{\theta})=\text{MSE}(\bm{w},-\bm{\theta})$ for all $\bm{\theta}\in \Theta$. Therefore, it suffices to consider maximizing the MSE over $\bm{\theta}\in\Theta$ such that $\theta_0\leq 0$. In addition, the following lemma implies that it suffices to consider $\bm{\theta} = (\theta_0, \ldots, \theta_n)'$ satisfying $\theta_i \geq \theta_0$ for all $i$.
The proofs of all theoretical results in the main text are provided in Appendix (ref). In the proof of Lemma (ref), we show that $\tilde{\bm{\theta}} = (\theta_0, \theta_1 + 2\cdot \max\{0, \theta_0 - \theta_1\}, \ldots, \theta_n + 2\cdot \max\{0, \theta_0 - \theta_n\})'$ satisfies $\text{MSE}(\bm{w},\bm{\theta}) \leq \text{MSE}(\bm{w}, \tilde{\bm{\theta}})$. We construct $\tilde{\bm{\theta}}$ by increasing $\theta_i$ to $\theta_0+\theta_0-\theta_i$ for each $i$ if $\theta_i$ is less than $\theta_0$. The new value is larger than $\theta_0$ by $\theta_0-\theta_i$. The change from $\bm{\theta}$ to $\tilde{\bm{\theta}}$ increases the variance while maintaining the Lipschitz constraint. Furthermore, we can show that this change results in a positive bias whose absolute value is larger than that of the bias at the original $\bm{\theta}$.
In view of Lemma (ref), we may consider the maximization of the MSE over $\bm{\theta}\in\Theta$ subject to the following restriction:
By calculating the derivatives of the MSE, we can show that $\text{MSE}(\bm{w},\bm{\theta})$ is nondecreasing in $\theta_j$ under ((ref)). To see this, note that
Because we have $\sum_{i \neq j} w_i \theta_i - \theta_0 \geq \left( \sum_{i \neq j} w_i - 1 \right) \theta_0 \geq 0$ for all $\bm{w} \in \mathcal{W}$ under ((ref)), it follows from ((ref)) that $\text{MSE}(\bm{w},\bm{\theta})$ is nondecreasing in $\theta_j$ under ((ref)). This monotonicity of the MSE implies that $\text{MSE}(\bm{w},(\theta_0,\theta_1,\ldots,\theta_n)')$ is maximized by setting $\theta_1,\ldots,\theta_n$ to the largest possible values that satisfy the Lipschitz constraint for each fixed value of $\theta_0$.
Formally, we define the largest possible values of $\theta_0, \theta_1,\ldots, \theta_n$ given $\theta_0 = t$ as \[ \tilde{\bm{\theta}}(t) \equiv \left( \tilde{\theta}_0(t), \tilde{\theta}_1(t), \ldots, \tilde{\theta}_n(t) \right)' \ \text{and} \ \tilde{\theta}_i(t) \equiv \min \{t + C \|R_i\|, 1/2 \} \ \text{for $i = 0,1, \ldots, n$,} \] as illustrated in Figure (ref). For any $\bm{\theta}=(\theta_0,\theta_1,\ldots,\theta_n)'\in\Theta$, we have $\theta_0=\tilde\theta_0(\theta_0)$ and $\theta_i \leq \tilde{\theta}_i(\theta_0)$ for $i=1,\ldots,n$. From ((ref)), if $\bm{\theta}\in\Theta$ satisfies ((ref)), we can increase the MSE by increasing $\theta_i$ to $\tilde{\theta}_i(\theta_0)$:
Also, $\tilde{\bm{\theta}}(t) \in \Theta$ for any $t \in [-1/2,1/2]$ because $\tilde{\bm{\theta}}(t)$ satisfies $\tilde{\bm{\theta}}(t) \in [-1/2, 1/2]^{n+1}$ and
where the second inequality follows from the reverse triangle inequality. Hence, we can reduce the ($n+1$)-dimensional maximization problem in ((ref)) to a one-dimensional problem with a single parameter $\theta_0$, as in the following theorem.
Next, we derive a weight vector that minimizes the maximum MSE. The following two lemmas show that the optimal weight vector is nonincreasing and that the $i$th element of the optimal weight vector is zero if $R_i$ is sufficiently far away from $R_0$.
Lemma (ref) shows that the optimal weight vector must be nonincreasing. In the proof of Lemma (ref), we show that if $\bm{w} \in \mathcal{W}$ satisfies $w_j < w_{j+1}$, then the maximum MSE can be reduced by swapping the positions of $w_j$ and $w_{j+1}$. By repeating this procedure until the weight vector becomes monotone, we can obtain $\tilde{\bm{w}} \in \mathcal{W}_0$ such that $\max_{\bm{\theta} \in \Theta} \text{MSE}(\tilde{\bm{w}},\bm{\theta}) \leq \max_{\bm{\theta} \in \Theta} \text{MSE}(\bm{w},\bm{\theta})$. Lemma (ref) shows that the $i$th element of the optimal weight vector is zero if $C\|R_i\| \geq 1/2$. By calculating the derivative of $\text{MSE}(\bm{w}, \tilde{\bm{\theta}}(\theta_0))$ with respect to $w_{i}$, we can show that $\text{MSE}(\bm{w}, \tilde{\bm{\theta}}(\theta_0))$ is nondecreasing in $w_i$ when $C\|R_i\| \geq 1/2$, and hence, setting $w_i = 0$ is optimal.
These two lemmas allow us to restrict our search space for the optimal $\bm{w}$ to nonincreasing vectors that place no weight on the observations with $C\|R_i\| \geq 1/2$. For notational simplicity, without loss of generality, we assume that our sample includes observations with $C \|R_i\| < 1/2$ only, so that $\mathcal{W}_0=\mathcal{W}_1$. Theorem (ref) and Lemma (ref) then imply that the minimax problem is reduced to
where
We now present a method for numerically solving the minimax problem ((ref)). We define $g(\bm{w};\theta_0) \equiv \text{MSE}(\bm{w},\tilde{\bm{\theta}}(\theta_0))$ and $\overline{g}(\bm{w}) \equiv \max_{\theta_0 \in [-1/2,0]} g(\bm{w};\theta_0)$. Because both $\bm{w} \mapsto \left( \sum_{i=1}^n w_i \theta_i - \theta_0 \right)^2$ and $\bm{w} \mapsto \sum_{i=1}^n w_i^2 \left( \frac{1}{4} - \theta_i^2 \right)$ are convex for any $\bm{\theta} \in \Theta$, $g(\bm{w};\theta_0)$ is also convex with respect to $\bm{w}$ for any $\theta_0 \in [-1/2, 0]$. As the maximum of convex functions is also convex, $\overline{g}(\bm{w})$ is a convex function. Therefore, the minimax problem ((ref)) becomes the following convex optimization problem with linear constraints:
We use nonlinear optimization via the augmented Lagrange method package_rsolnp,Ye_1987 in the implementation in simulations and applications.
Many existing studies have considered minimax estimation problems for unbounded outcomes with known variances, primarily motivated by the Gaussian model. We compare our proposed estimator with a Gaussian-motivated minimax linear estimator when the underlying data-generating process is a binary-outcome model.
As a Gaussian-motivated estimator, we consider the minimax linear estimator for an unbounded space of mean vectors with known variances under a smoothness restriction, following the existing minimax analysis in RD designs Armstrong.Kolesar2018,Imbens.Wager2019. Note that if the outcome $Y_i$ is normally distributed, that is, $Y_i \sim N(p_i, \sigma_i^2)$, the MSE of a linear estimator $\hat{p}_0(\bm{w}) = \frac{1}{2} + \sum_{i=1}^n w_i \left( Y_i-\frac{1}{2} \right)$ with $\bm{w}\in \mathbb{R}^{n}$ is given by \[ E\left[ (\hat{p}_0(\bm{w}) - p_0)^2 \right] \ = \ \left\{ \frac{1}{2} + \sum_{i=1}^n w_i \left( p_i - \frac{1}{2} \right) - p_0 \right\}^2 + \sum_{i=1}^n w_i^2 \sigma_i^2. \] Letting $\theta_i = p_i - 1/2$, the MSE can be written as follows: \[ \left( \sum_{i=1}^n w_i \theta_i - \theta_0 \right)^2 + \sum_{i=1}^n w_i^2 \sigma_i^2. \] As a smoothness restriction, we impose the Lipschitz constraint, leading to the following parameter space: \[ \Theta_g \ \equiv \ \left\{ \bm{\theta} \in \mathbb{R}^{n+1} : |\theta_i - \theta_j| \leq C \| R_i-R_j \| \ \text{for all $i$ and $j$} \right\}. \] The minimax linear estimator is a solution of the following problem:
We refer to the linear estimator that solves ((ref)) as the Gaussian estimator.\footnote{Note that this estimator is a minimax linear estimator without normality of $Y_i$ as long as the variance is known and the parameter space is $\Theta_g$. The normality of $Y_i$ is exploited for finite-sample valid inference based on a linear estimator.} The above minimax problem ((ref)) differs from the original binary-outcome problem ((ref)) in three respects. First, the minimum in ((ref)) is considered among all linear estimators, including those with negative weights. Second, the parameter space in ((ref)) is unbounded. Finally, but most importantly, the variance in ((ref)) does not depend on the parameter $\bm{\theta}$; hence, the maximum MSE is attained at the parameter values that maximize the squared bias.
In Appendix (ref), we derive the form of the optimal weights that solve the minimax problem ((ref)) by applying the results of donoho1994 to our Gaussian setting. We show that the optimal weights satisfy $\sum_{i=1}^n w_i = 1$ and $w_i \geq 0$ for all $i$. Hence, the minimax problem ((ref)) can be solved by minimizing the maximum MSE over $\mathcal{W}$. Specifically, the Gaussian estimator is obtained by solving the following quadratic program: \[ \min_{\bm{w}} \left\{ C^2 \left( \sum_{i=1}^n w_i \|R_i\| \right)^2 + \sum_{i=1}^n w_i^2 \sigma_i^2 \right\} \ \ \text{s.t.} \ \ \sum_{i=1}^n w_i = 1 \ \text{and} \ w_i \geq 0 \ \text{for all $i$}, \] where $C^2\left( \sum_{i=1}^n w_i \|R_i\| \right)^2$ is the maximum squared bias of the estimator $\hat{p}_0(\bm{w})$ with $\sum_{i=1}^n w_i = 1$ over $\Theta_g$.
We compare the maximum MSE of the proposed estimator with that of the Gaussian estimator in a setting in which the true model is the binary-outcome model described in Section (ref). The Gaussian estimator requires the specification of variance. In the following, we focus on the Gaussian estimator with $\sigma_1^2 = \cdots = \sigma_n^2 = 1/4$ because the variance of a binary variable is less than or equal to $1/4$. Define \[ \hat{\bm{w}} \in \mathrm{arg} \min_{\bm{w} \in \mathcal{W}} \max_{\bm{\theta} \in \Theta} \mathrm{MSE}(\bm{w},\bm{\theta}) \ \ \text{and} \ \ \tilde{\bm{w}} \in \mathrm{arg} \min_{\bm{w} \in \mathcal{W}} \max_{\bm{\theta} \in \Theta_g} \mathrm{MSE}_g(\bm{w},\bm{\theta}), \] where
Then, $\hat{p}_0(\hat{\bm{w}})$ is the minimax linear shrinkage estimator when $Y_i$ is binary, and $\hat{p}_0(\tilde{\bm{w}})$ is the minimax linear estimator when $Y_i \sim N(p_i,1/4)$. The following lemma compares the maximum MSEs of $\hat{p}_0(\hat{\bm{w}})$ and $\hat{p}_0(\tilde{\bm{w}})$ when $Y_i$ is binary and the parameter space is bounded.
Lemma (ref) provides lower and upper bounds of the ratio of the maximum MSEs. As $\hat{\bm{w}}$ minimizes $\max_{\bm{\theta} \in \Theta} \mathrm{MSE}(\bm{w},\bm{\theta})$ over $\mathcal{W}$, the lower bound is trivial. In the proof of Lemma (ref), we derive the upper bound based on an upper bound on the numerator and a lower bound on the denominator.
Although the finite-sample bounds in Lemma (ref) may be loose, sharper bounds can be obtained if we consider an asymptotic setting in which the sample size increases. In the following, we consider a triangular array $\{(R_{n,1},\ldots,R_{n,n})\}_{n\in \mathbb{N}}$, where $(R_{n,1},\ldots,R_{n,n})$ is a deterministic vector that collects the values of the running variable when the sample size is $n$. We fix the value of the Lipschitz constant $C$ as $n$ varies. In this asymptotic regime, we show that under mild conditions, the convergence rate of $\hat{p}_0(\hat{\bm{w}})$ is $O_p(n^{-1/3})$ and the ratio of the maximum MSEs of $\hat{p}_0(\hat{\bm{w}})$ and $\hat{p}_0(\tilde{\bm{w}})$ converges to one as $n \to \infty$. For the brevity of the notation, we suppress the first index $n$ of $(R_{n,1},\ldots,R_{n,n})$ below.
To establish the asymptotic result, we consider a univariate running variable $R_i$ and assume that it is bounded and that the empirical distribution of $\|R_i\|$ is bounded above and below by linear functions.\footnote{The convergence holds under a weaker condition, which may be plausible for a multivariate running variable. See Remark (ref) for a discussion of the general case.}
Figure (ref) illustrates Assumption (ref) (ii). For example, when $R_i = i/n$ for all $i=1, \ldots, n$, this assumption is satisfied for $0 < c_0 < 1 < c_1$. More generally, Assumption (ref) (ii) requires the empirical distribution $F_n(x)$ to be bounded by a pair of linear functions.
Theorem (ref) shows that the convergence rate of $\hat{p}_0(\hat{\bm{w}})$ is $O_p(n^{-1/3})$. This convergence rate is the same as that of standard nonparametric estimators under the Lipschitz constraint in univariate RD designs. Theorem (ref) also shows that the maximum MSE of $\hat{p}_0(\tilde{\bm{w}})$ is asymptotically identical to that of $\hat{p}_0(\hat{\bm{w}})$. The Gaussian estimator $\hat{p}_0(\tilde{\bm{w}})$ minimizes the maximum MSE when $Y_i \sim N(p_i, 1/4)$ and the parameter space is unbounded. This result implies that the Gaussian estimator is asymptotically optimal in terms of the maximum MSE for a particular sequence of distributions of the running variable, even when outcomes are binary.
Although the efficiency gain from our estimator relative to the Gaussian estimator can be small in large samples, their behaviors are quite different in finite samples. We demonstrate the finite-sample comparisons of our estimator with the Gaussian estimator through numerical analyses. Figures (ref) and (ref) plot the weights $w_1,\ldots,w_n$ for samples of observations whose running-variable values are equally spaced between $0$ and $1$. Figure (ref) plots the weights of our estimator (rdbinary) and the Gaussian estimator (gauss) for a sample size of $50$ and four Lipschitz constant values. Figure (ref) shows the corresponding plots for a sample size of $500$. The weights of the Gaussian estimator are computed under the assumption that the variance is homoskedastic and $1/4$ for all units, as in Section (ref). For a small sample size of $50$, our estimator exhibits moderate size of shrinkage, whereas the Gaussian estimator exhibits no shrinkage. For $C > 0$, the weights of the Gaussian estimator are triangular, while our estimator's weights display mild nonlinearity. Also, the Gaussian weights have thicker tails than ours. These differences in shape arise because the Gaussian estimator is constructed under homoskedasticity and the maximum possible variance of $1/4$, whereas our estimator optimizes the weights under potential heteroskedasticity.
By contrast, the two estimators appear almost equivalent for a sufficiently large sample size of $500$. The shape of our estimator's weights is still sharper than that of the Gaussian weights for $C = 1$; however, the differences between the two weights are negligible compared to the case with a small sample size of $50$.
More distinct differences appear in the maximum root MSEs in small samples. Figure (ref) reports the ratio of the maximum root MSE of the Gaussian estimator to that of our estimator, calculated in the binary-outcome model. For a small sample size of $50$, the Gaussian estimator has $5\%$ to $20\%$ larger root MSEs than our estimator. Hence, our estimator achieves substantial improvement over the Gaussian estimator in small samples.
Nevertheless, these ratios shrink as the sample size increases and the gaps shrink below $5\%$ when $N = 500$. This property is consistent with the theoretical result that the ratio of the worst-case MSEs converges to one as the sample size increases. In summary, our estimator is substantially different from and superior to the Gaussian estimator in finite samples, whereas the two estimators behave similarly in large samples.
In this section, we return to the original setup introduced in Section (ref), where we observe both the treated sample $\{ Y_{i,+}, R_{i,+} \}_{i=1}^{n_{+}}$ and the untreated sample $\{ Y_{i,-}, R_{i,-} \}_{i=1}^{n_{-}}$. We propose an inference procedure with respect to $\tau \equiv f(1,R_0)-f(0,R_0)$ based on a given linear shrinkage estimator. Let $p_{i,+} \equiv f(1,R_{i,+})$, $p_{i,-} \equiv f(0,R_{i,-})$, and $R_{0,+}=R_{0,-}=0$, so that $Y_{i,+}$ and $Y_{i,-}$ follow Bernoulli distribution with parameters $p_{i,+}$ and $p_{i,-}$, respectively. Similar to the previous sections, we assume that $p_{i,+}$ and $p_{i,-}$ satisfy $\bm{p}_{+} \equiv (p_{0,+}, p_{1,+}, \ldots, p_{n_{+},+})' \in \mathcal{P}_{+}$ and $\bm{p}_{-} \equiv (p_{0,-}, p_{1,-}, \ldots, p_{n_{-},-})' \in \mathcal{P}_{-}$, where
We propose an inference procedure of $\tau = p_{0,+}-p_{0,-}$ based on the estimator $\hat{\tau} \equiv \hat{p}_{0,+}(\bm{w}_{+}) - \hat{p}_{0,-}(\bm{w}_{-})$, where
Our inference procedure is valid for any linear estimator with nonnegative weights (even if $\sum_{i=1}^{n_{+}}w_{i,+}>1$ or $\sum_{i=1}^{n_{-}}w_{i,-}>1$) when the outcome is binary. Hence, we can conduct an inference using the linear shrinkage estimator proposed in the previous sections. The following argument does not extend to general bounded outcomes. In Appendix (ref), we consider an inference procedure for general bounded outcomes based on Hoeffding's inequality.
We provide finite-sample valid confidence intervals by inverting tests that are uniformly valid over the Lipschitz class. We begin our analysis with a one-sided test. Using a uniformly valid one-sided test, we then construct a uniformly valid two-sided test and confidence interval.
Specifically, we consider a one-sided test for the following null and alternative hypotheses:
We propose the following testing procedure based on the linear estimator $\hat{\tau}$:
where $\gamma$ is a critical value. The critical value $\gamma$ must satisfy $P_{\bm{p}}(\hat{\tau} - \tau_0 > \gamma) \leq \alpha$ for any parameter $\bm{p} \equiv (\bm{p}_{+}', \bm{p}_{-}')' \in \mathcal{P}_{\ast} \equiv \mathcal{P}_{+} \times \mathcal{P}_{-}$ satisfying $H_0$. Hence, we must choose the critical value $\gamma^{\ast}(\tau_0)$ such that
where $\mathcal{P}(\tau_0) \equiv \left\{ \bm{p} \in \mathcal{P}_{\ast} : p_{0,+} - p_{0,-} = \tau_0 \right\}$. The critical value $\gamma^*(\tau_0)$ provides a uniformly valid one-sided test for finite samples.
To obtain an appropriate critical value, we must calculate $\max_{\bm{p} \in \mathcal{P}(\tau_0)} P(\hat{\tau} > \gamma)$. The following theorem shows that we can calculate $\max_{\bm{p} \in \mathcal{P}(\tau_0)} P(\hat{\tau} > \gamma)$ by optimizing a single parameter.
Theorem (ref) is obtained using first-order stochastic dominance. Suppose that $(Y_1, \ldots, Y_n)' \in \{0,1\}^n$ and $(\tilde{Y}_1, \ldots, \tilde{Y}_n)' \in \{0,1\}^n$ follow $n$-dimensional independent Bernoulli distributions with parameters $\bm{p} \in \mathbb{R}^n$ and $\tilde{\bm{p}} \in \mathbb{R}^n$, respectively, and each element of $\bm{p}$ is smaller than or equal to that of $\tilde{\bm{p}}$. Then, if $w_i$ is nonnegative for all $i$, $\sum_{i=1}^n w_i \tilde{Y}_i$ has first-order stochastic dominance over $\sum_{i=1}^n w_i Y_i$. Hence, if we fix $p_{0,+}$ and $p_{0,-}$, then $P_{\bm{p}}(\hat{\tau} > \gamma)$ is maximized at $\bm{p} = \left( \tilde{\bm{p}}_{+}(p_{0,+})', \tilde{\bm{p}}_{-}(p_{0,-})' \right)'$, implying ((ref)) holds.
From Theorem (ref), we can obtain the critical value $\gamma^{\ast}(\tau_0)$ satisfying ((ref)) using the following algorithm:
Next, we construct a uniformly valid two-sided test and confidence interval using the one-sided test proposed in Section (ref). We consider the following null and alternative hypotheses:
Similar to the one-sided test, we propose the following testing procedure based on the linear estimator $\hat{\tau}$:
where the critical values $\gamma_l$ and $\gamma_r$ must satisfy $P_{\bm{p}}(\hat{\tau} \not\in [\gamma_l,\gamma_r] ) \leq \alpha$ under $H_0$. Hence, we must choose the critical values $\gamma_l^{\ast}(\tau_0)$ and $\gamma_r^{\ast}(\tau_0)$ such that
However, it is challenging to derive a simple expression for the maximum of the probability $P_{\bm{p}}\left( \hat{\tau}\not\in [\gamma_l,\gamma_r] \right)$, unlike for the one-sided testing. Therefore, we instead calculate an upper bound on the maximum of $P_{\bm{p}}\left( \hat{\tau} \not\in [\gamma_l,\gamma_r] \right)$:
where $\pi_r(\gamma_r) \equiv \max_{\bm{p} \in \mathcal{P}(\tau_0)} P_{\bm{p}}( \hat{\tau} > \gamma_r)$ and $\pi_l(\gamma_l) \equiv \max_{\bm{p} \in \mathcal{P}(\tau_0)} P_{\bm{p}}( \hat{\tau} < \gamma_l)$. We can calculate $\pi_r(\gamma_r)$ as in Section (ref) and $\pi_l(\gamma_l)$ similarly. We then propose the following critical values $\gamma_r^{\ast}(\tau_0)$ and $\gamma_l^{\ast}(\tau_0)$:
By construction, these critical values $\gamma_r^{\ast}(\tau_0)$ and $\gamma_l^{\ast}(\tau_0)$ satisfy ((ref)).
We obtain the confidence region of $\tau$ by inverting the testing procedure. We define $\widehat{CR}_{1-\alpha}$ as the set of hypothesized values that are not rejected by the proposed two-sided test, that is, \[ \widehat{CR}_{1-\alpha} \ \equiv \ \left\{ \tau_0 \in [0,1] : \gamma_l^{\ast}(\tau_0) \leq \hat{\tau} \leq \gamma_r^{\ast}(\tau_0) \right\}. \] By construction, $\widehat{CR}_{1-\alpha}$ satisfies \[ \min_{\bm{p} \in \mathcal{P}_{\ast}} P_{\bm{p}} \left( \tau \in \widehat{CR}_{1-\alpha} \right) \ \geq \ 1-\alpha. \] In other words, this confidence region is uniformly valid over the Lipschitz class.
This confidence region is an interval. As discussed in Remark (ref), $\gamma_r^{\ast}(\tau_0)$ is increasing in $\tau_0$. Similarly, $\gamma_l^{\ast}(\tau_0)$ is increasing in $\tau_0$. Suppose that $t_1 < t_2$ and $t_1, t_2 \in \widehat{CR}_{1-\alpha}$. Then, for any $t \in [t_1,t_2]$, we obtain \[ \gamma_l^{\ast}(t) \leq \gamma_l^{\ast}(t_2) \leq \hat{\tau} \ \ \text{and} \ \ \hat{\tau} \leq \gamma_r^{\ast}(t_1) \leq \gamma_r^{\ast}(t). \] Hence, any $t$ within the interval $[t_1,t_2]$ must be contained in the confidence region $\widehat{CR}_{1-\alpha}$, which means that $\widehat{CR}_{1-\alpha}$ is an interval. Consequently, searching for the boundary points of $\widehat{CR}_{1-\alpha}$ is sufficient for constructing the confidence interval.
We demonstrate the performance of our estimator relative to existing estimators in Monte Carlo simulations. We compare our estimator (rdbinary) with three different estimators: (1) the Gaussian estimator (gauss) with homoskedastic variance $\sigma_i^2 = 1/4$ as in Section (ref); (2) the XU20171's estimator (rd.mnl), which is specific to multinomial outcomes, including binary outcomes as a special case; and (3) the Calonico.Cattaneo.Titiunik2014's estimator (rdrobust).\footnote{We obtained the rd.mnl package files from the author's personal website (\url{https://sites.google.com/view/kelixu/home}). For rd.mnl and \textit{rdrobust}, we use their default specifications with bias-corrected robust estimation and inference.}
We compare their performance for three sample sizes $(N \in \{50, 100, 500\})$ of observations whose running-variable values are equally spaced between $-1$ and $1$. We consider the following three different models of the conditional mean of a binary dependent variable: (1) the Lee2008 model, which is a polynomial approximation of the conditional mean for Lee2008's data and is frequently used in simulation studies for RD designs; (2) the “worst-case” model, which is the parameter value $\bm{p}$ maximizing the MSE of any linear shrinkage estimator among parameter values such that $p_{0,+}=p_{0,-}=1/2$;\footnote{Note that the worst-case MSE of a linear shrinkage estimator is not necessarily attained at the parameter values of this model, since $p_{0,+}$ and $p_{0,-}$ are fixed at $1/2$.} and (3) a flat model in which the conditional probability is constant at $0.5$. The three designs are illustrated in Figures (ref)--(ref). For each model, the dependent variable takes $1$ with the probability specified as mean and otherwise $0$.
We consider the estimation and inference of $\tau=p_{0,+}-p_{0,-}$. We use the true value of the Lipschitz constant $C$ for each design to implement our proposed method and the Gaussian method. Our proposed estimator for $\tau$ is given by $\hat{\tau} = \hat{p}_{0,+}(\hat{\bm{w}}_{+}) - \hat{p}_{0,-}(\hat{\bm{w}}_{-})$, where $\hat{\bm{w}}_{+}$ and $\hat{\bm{w}}_{-}$ are chosen to minimize the worst-case MSE for the estimation of $p_{0,+}$ and $p_{0,-}$, respectively, as in Section (ref). We then use $\hat{\tau}$ to construct a two-sided confidence interval for $\tau$ following the procedure in Section (ref).\footnote{We computed the pair of critical values $\gamma_{r}^*(\tau_0)$ and $\gamma_{l}^*(\tau_0)$ by computing $\pi_{r}(\gamma_r)$ and $\pi_l(\gamma_l)$ with $3000$ draws of an $n$-dimensional Bernoulli random vector for each. The confidence intervals were constructed by inverting tests evaluated at $300$ grid points of $\tau_0$.} An alternative, the Gaussian estimator, is $\tilde{\tau} = \hat{p}_{0,+}(\tilde{\bm{w}}_{+}) - \hat{p}_{0,-}(\tilde{\bm{w}}_{-})$, where $\tilde{\bm{w}}_{+}$ and $\tilde{\bm{w}}_{-}$ minimize the worst-case MSE for the estimation of $p_{0,+}$ and $p_{0,-}$, respectively, under the misspecified model where $Y_i \sim N(p_i,1/4)$, as in Section (ref). Following kolesar2018discrete and Armstrong2021ATE, we construct a two-sided fixed-length confidence interval centered at $\tilde{\tau}$ with finite-sample validity under the Gaussian model. Specifically, the $100\cdot (1-\alpha)\%$ confidence interval is given by $\left(\tilde{\tau}\pm \rm{cv}_\alpha\left(\rm{maxbias}(\tilde{\tau})/\rm{sd}(\tilde{\tau})\right)\cdot \rm{sd}(\tilde{\tau})\right)$, where $\rm{maxbias}(\tilde{\tau})$ denotes the maximum bias of $\tilde{\tau}$ under the Lipschitz class and $\rm{cv}_\alpha(b)$ denotes the $1-\alpha$ quantile of $|N(b,1)|$, the folded normal distribution with location and scale parameters $(b,1)$.
First, we demonstrate the point-estimation properties of our proposed estimator. Tables (ref) and (ref) compare the root MSE and bias for the estimation of the ATE at the cutoff, computed from $3000$ replication draws. Table (ref) compares three different sample sizes under the Lee model. For all sample sizes, our estimator has substantially smaller MSEs than the other estimators. Furthermore, the differences decrease as the sample size increases, and the MSEs are relatively similar for $N = 500$. We note that the same pattern is confirmed for different designs with different Lipschitz constants $C$. Hence, our estimator is superior to the existing estimators in small samples, whereas their behaviors resemble in larger samples.
Second, we demonstrate the inference properties of our proposed method. Tables (ref) and (ref) compare the average length and coverage probability of the four confidence intervals, computed from $5000$ replication draws. Table (ref) shows that our confidence interval has shorter lengths with guaranteed coverage than rd.mnl and rdrobust for different sample sizes. Unlike the point estimation results, the differences in lengths remain similar as the sample size increases. Note that the Gaussian confidence interval has slightly shorter lengths while achieving the $95\%$ coverage for the Lee design. Nevertheless, the Gaussian confidence interval does not guarantee uniform coverage across designs; for instance, it falls below $95\%$ for the flat design when $N=100$. This behavior is consistent with the fact that the Gaussian confidence interval is constructed under a missspecified model, where the outcomes, and hence the linear estimators, are assumed to follow a normal distribution. In contrast, our proposed confidence interval provides guaranteed coverage, a feature that makes it preferable. We also note that the rdrobust confidence interval is based on large-sample asymptotics and not specifically designed for binary outcomes, which results in unsatisfactory coverage properties for all designs, especially with small samples.
We apply our estimator to a small-sample RD study by Brollo.Nannicini.Perotti.Tabellini2013 brolloReplicationDataPolitical2019. Brollo.Nannicini.Perotti.Tabellini2013 exploit a regional fiscal rule in Brazil to study the impact of an additional government fiscal transfer on the frequency of corruption in local politics. In Brazil, 40 percent of municipal revenue comes from the Fundo de Participação dos Municipios (FPM), which is allocated based on the population size of municipalities. Specifically, each municipality is allocated to one of the nine brackets based on its population. The bracketing fiscal rule induces population thresholds that discontinuously alter the amount of FPM transfers. Following Brollo.Nannicini.Perotti.Tabellini2013, we focus on the first seven thresholds because of sample size limitations.
This study is particularly suited for our method for two reasons. First, their primary dependent variables are binary indicators. Specifically, they study the impact of the fiscal rule on two measures of corruption indicators:
Second, the sample size is relatively small. In particular, within each cutoff neighborhood, the sample size is limited to less than $400$ and is mostly around $100$ to $200$. In these small samples, our estimator is expected to be superior to other estimators that are based on asymptotic approximations.
The following tables present our rdbinary estimates and rdrobust estimates.\footnote{The original study runs global polynomial estimations for each cutoff neighborhood, as well as for the whole sample by pooling across cutoff neighborhoods. Their primary estimation is the fuzzy design, but we focus on the reduced-form sharp design estimates.} Tables (ref) and (ref) report the pooled estimates over multiple cutoffs for the broad and narrow corruption indicators. $Crot$ denotes the rule-of-thumb value for the Lipschitz constant $C$, which is the largest (in absolute value) slope estimate from the binscatter estimation using the binsreg package CattaneoCrumpFarrellFeng_2024. In all tables, we report the point estimates and confidence intervals for three different values of the constant $C$: $Crot$; one-half of $Crot$; and $1.5$ times $Crot$.
For both indicators, our rdbinary estimates appear to be similar to the rdrobust estimates, which are valid for large samples. The sample size is $1,202$ for the entire pooling sample and, hence, is sufficiently large for the rdrobust estimator.\footnote{We do not report rd.mnl estimates because rd.mnl sometimes fails to select a bandwidth in this dataset, particularly for small samples.} For both methods and outcomes, the $95\%$ confidence intervals include $0$. This finding differs from that of the original study, which reports significant positive effects on the frequency of corruption. This difference highlights the importance of using local nonparametric estimations in RD designs.
By pooling samples across multiple cutoffs, we obtain a sufficiently large sample. Nevertheless, heterogeneity across different cutoffs may be of interest, as the original study explores cutoff-specific estimates. However, only a few hundred observations are available around each individual cutoff, and the asymptotic approximation may not perform well in such small samples.
Tables (ref), (ref), and (ref) present our rdbinary and rdrobust estimates of the impact on the broad corruption indicator for seven different subsamples around each individual cutoff. See the Online Appendix for qualitatively similar results for the narrow corruption indicator. For all specifications, the confidence intervals for each subsample are much wider than those for the pooled sample. Nevertheless, our rdbinary method tends to yield much shorter confidence intervals than rdrobust. For example, Cutoff 3 has a sample size of $225$, which is too small for rdrobust to offer any meaningful implications from its confidence interval. By contrast, our rdbinary confidence intervals offer reasonable lower bounds for the impact on the broad corruption measure, which are not very negative compared with the lower bound of the confidence interval from \textit{rdrobust}.
Empirical studies often use RD designs with small samples, where estimation and inference are particularly challenging. Methods that rely on large-sample properties may not perform well in such settings. Several finite-sample minimax estimators have been proposed. However, these estimators typically require knowledge of the variance, which is generally unavailable. Furthermore, the finite-sample validity of inference procedures based on these estimators requires the normality of the outcome variable.
In this study, we provide a minimax optimal estimator for RD designs with a binary outcome variable, together with its inference procedure. The only tuning parameter of our proposed estimator is the Lipschitz constant; our estimator does not require specification of the conditional variance function, unlike existing minimax estimators in RD designs. Our estimator is also applicable to any bounded outcome variable. Thus, we offer a practical method that serves as a last resort for RD studies with relatively small effective sample sizes.
We demonstrate that the proposed estimator is superior to existing estimators in finite samples through numerical and simulation exercises. In a numerical exercise, we show that our estimator is $5$ to $20$% more efficient in the worst-case root MSEs than a feasible version of the existing minimax optimal estimators for extremely small samples. Through simulation studies, we show that our estimator has much smaller MSEs than existing methods for sufficiently small sample sizes. Furthermore, we demonstrate that our inference procedure generates shorter confidence intervals with guaranteed coverage rates compared with existing large-sample methods. In the empirical application to a small-sample RD study, we document that our method generates similar results to the standard large-sample procedure for sufficiently large samples but provides much more informative results for small samples.
Our contribution provides a critical baseline for developing estimation procedures for binary or limited outcome variables in RD designs. Recent studies, such as noack_bias_aware_2024, have considered bias-aware inference for fuzzy RD designs. As mentioned in the Introduction, the binary treatment status is a primary dependent variable in the first stage of the fuzzy design. Applying our results is not necessarily straightforward in that context because the first-stage estimand appears in the denominator of the target parameter. We reserve further extensions and generalizations of our results for future research.