Isotonic Regression Discontinuity Designs

Isotonic Regression Discontinuity Designs

Introduction

Regression discontinuity designs, see thistlethwaite1960regression, are widely recognized as one of the most credible quasi-experimental strategies for identifying and estimating causal effects. In a nutshell, such designs exploit the discontinuity in the treatment assignment probability around the cutoff value of some running variable. The discontinuous treatment assignment probability frequently occurs due to laws and regulations governing economic and political life. A comprehensive list of empirical applications using regression discontinuity designs can be found in lee2010regression; see also imbens2008regression, and cattaneo2019practical for the methodological review and imbens2009recent, abadie2018econometric for their place in the program evaluation literature. On the methodological side, in the seminal paper hahn2001identification translate regression discontinuity designs into the potential outcomes framework and introduce the local polynomial nonparametric estimator to estimate causal effects in sharp and fuzzy designs.

Regression discontinuity designs encountered in the empirical practice are frequently monotone. Indeed, development and educational programs are often prescribed based on poverty or achievement scores that are monotonically related to average outcomes. For instance, when evaluating the effect of subsidies for fertilizers on yields, the yield per acre is expected to be non-decreasing in the size of the farm due to the increasing returns to scale. Alternatively, when evaluating the effectiveness of the cash transfers program on the households food expenditures, we expect that more affluent households spend, on average, more on food, since food is a normal good; see Section (ref) in the Supplementary Material for a sample of other examples drawn from the empirical research.

Despite this prevalence in the empirical practice, little is known about how monotonicity can be incorporated in the estimation procedure.\footnote{Shape restrictions in regression discontinuity designs is a relatively unexplored area. To the best of our knowledge, armstrong2015adaptive is the only existing work focusing on the k-nearest neighbors estimator.} Estimation and inference in regression discontinuity designs require estimating the conditional mean functions at the boundary of the support. In this paper, we focus on the isotonic regression estimator and establish its formal statistical properties for the boundary point.\footnote{This estimator is relatively unknown in econometrics and to the best of our knowledge has not been previously considered in the RDD setting; see an excellent review paper chetverikov2018econometrics for a comprehensive discussion of this estimator and further references on the econometrics of shape restrictions.} The isotonic regression estimator originates from the work of ayer1955empirical, and brunk1956inequality. brunk1970estimation derives its asymptotic distribution at the interior point under restrictive assumptions that the regressor is deterministic and regression errors are homoskedastic. His treatment builds upon the ideas of rao1969estimation, who derived the asymptotic distribution of the monotone density estimator, also known as the Grenander estimator; see grenander1956theory. wright1981asymptotic provides the final characterization of the large sample distribution for the interior point when the regressor is random and regression errors are heteroskedastic.

However, to the best of our knowledge, little is known about the behavior of the isotonic regression at the boundary point, which is a building block of our isotonic regression discontinuity (iRDD) estimators. This situation contrasts strikingly with the local polynomial estimator, whose boundary behavior is well-understood; see fan1992variable. Most of the existing results for isotonic estimators at the boundary are available only for the Grenander estimator; see woodroofe1993penalized, kulikov2006behavior, and balabdaoui2011grenander. More precisely, we know that the Grenander estimator is inconsistent at the boundary of the support and that the consistent estimator can be obtained with additional boundary correction or penalization. At the same time, some isotonic estimators, e.g., in the current status model, are consistent at the boundary without corrections; see durot2018limit. anevski2002monotone discuss the inconsistency of the isotonic regression at the discontinuity point with deterministic equally spaced covariate and homoskedasticity. However, anevski2002monotone do not discuss whether the isotonic regression with random covariate is consistent at the boundary of its support and do not provide a consistent estimator even in the restrictive equally spaced fixed design case.

In this paper, we aim to understand the behavior of the isotonic regression with a random regressor at the boundary of its support. We show formally that the isotonic regression estimator is inconsistent and tends to underestimate the value of the regression at the left boundary. The inconsistency is related to the extreme-value behavior of the closest to the boundary observation. We introduce boundary-corrected estimators and derive large sample approximations to corresponding distributions. The major technical difficulty when deriving asymptotic distributions in this setting is to establish the tightness of the maximizer of a certain empirical process. This condition is typically needed to apply the argmax continuous mapping theorem of kim1990cube. The difficulty stems from the fact that conventional tightness results of kim1990cube and van2000weak do not always apply. For the Grenander estimator, kulikov2006behavior suggest a solution to this problem based on the Koml\'{o}s-Major-Tusn\'{a}dy strong approximation. In our setting, this approach entails the strong approximation to the general empirical process, see koltchinskii1994komlos and chernozhukov2015constrained, which is more problematic to apply due to slower convergence rates. Consequently, we provide the alternative generic proof which does not rely on the strong approximation and which might be applied to other boundary-corrected shape-constrained estimators.

Since the asymptotic distribution is not pivotal, we introduce a novel trimmed wild bootstrap procedure and establish its consistency. The procedure consists of trimming values of the estimated regression function that are very close to the boundary when simulating wild bootstrap samples. Somewhat unexpectedly, we discover that the trimming and the appropriate boundary correction restores the consistency of the wild bootstrap without additional nonparametric smoothing or subsampling, which is typically needed in such settings.\footnote{In contrast, the bootstrap typically fails at the interior point; see kosorok2008bootstrapping, sen2010inconsistency, guntuboyina2018nonparametric, and patra2018consistent for the discussion of this problem and various case-specific remedies, and cattaneo2017bootstrap for generic solutions that apply to a class of cube-root consistent estimators.}

The paper is organized as follows. In Section (ref), we develop the asymptotic theory for the isotonic regression estimator at the boundary point. The reader interested more in the regression discontinuity designs can skip this section and go directly to Section (ref), where we discuss isotonic sharp and fuzzy designs. In Section (ref), we study the finite sample performance of the iRDD estimator with Monte Carlo experiments. Section (ref) estimates the effect of incumbency using the sharp iRDD on the data of lee2008randomized. Section (ref) concludes. Additional results, proofs, and Monte Carlo experiments appear in the Appendix and the Supplementary Material.

Isotonic regression at the boundary

In this section, we study the isotonic regression estimator at the boundary point. Readers interested in regression discontinuity designs only may skip this section and go directly to Section (ref).

Estimator

We focus on the generic nonparametric regression

where the conditional mean function $m(x)=\ensuremath{\mathbb{E}}[Y|X=x]$ is assumed to be montone. For simplicity of presentation, we normalize the support to $[0,1]$ and assume that $m$ belongs to the set of non-decreasing functions on $[0,1]$

For a sample $(Y_i,X_i)_{i=1}^n$, the isotonic regression estimator, denoted $\hat m$, solves the constrained least-squares problem

see brunk1970estimation. Note that the estimator is uniquely determined at data points and is conventionally interpolated as a piecewise constant function elsewhere. Let $X_{(1)}<X_{(2)}<\dots <X_{(n)}$ be the ordered values of the covariate and let $(Y_{(1)},Y_{(2)}\dots,Y_{(n)})$ be the corresponding observations of the outcome variable. It is easy to see that the isotonic regression is alternatively characterized as a solution to

While, the isotonic regression is a solution to the convex optimization problem and can be obtained with standard constrained convex optimization packages, an efficient way to compute the estimator is via the pool adjacent violators algorithm; see ayer1955empirical. Although the isotonic regression features the number of estimated parameters of the same magnitude as the sample size, the pool adjacent violators algorithm is typically computationally cheaper than nonparametric kernel estimators, and its computational complexity is closer to that of the OLS estimator.

Large sample distribution

We are interested in estimating the value of the regression function at the boundaries of its support $[0,1]$. More precisely, we focus on the regression at the left boundary, denoted $m(0)=\lim_{x\downarrow0}m(x)$. The natural estimator of $m(0)$ is $\hat m(X_{(1)})$. Unfortunately, it follows from Theorem (ref) in the appendix that this estimator is inconsistent. The inconsistency occurs because $X_{(1)}$ converges to zero too fast according to laws of the extreme value theory. More precisely, the estimator $\hat m(X_{(1)})$ underestimates $m(0)$ and $\hat m(X_{(n)})$ overestimates $m(1)$ in small samples; see Section (ref) in the Appendix for more details.\footnote{We are grateful to Matt Masten for this observation. More generally, in the isotonic regression discontinuity design, the isotonic regression estimator underestimates the causal effect at the cutoff, which can be used as a tuning-free lower bound on the estimated causal effect.}

To estimate $m(0)$ consistently, we the consider boundary-corrected estimators $\hat m(cn^{-a})$ for some $c,a>0$ and study the large sample behavior of this estimator for a range of values $c,a>0$. The following assumption imposes a mild restriction on the distribution of the data.

For a stochastic process $\{Z_t: t\in A\subset\ensuremath{\mathbf{R}}\}$, let $D^L_{A}(Z_t)(s)$ denote the left derivative of the greatest convex minorant at $s\in A$.\footnote{The greatest convex minorant of a function $g:A\to\ensuremath{\mathbf{R}}$ is defined as the maximal convex function $h$ such that $h(x)\leq g(x),\forall x\in A$.} We say that $m$ is $\gamma$-H\"{o}lder continuous in the neighborhood of zero if there exists a constant $C<\infty$ such that for all $x>0$ in this neighborhood

Our first result describes the large sample behavior of boundary corrected isotonic regression estimators for a range of boundary corrections.

The proof of this result can be found in the appendix. The most challenging part of the proof is establishing tightness when $a\in(1/3,1)$. The difficulty comes from the fact that the quadratic term vanishes in this asymptotic regime and the standard tightness results for isotonic estimators do not apply.\footnote{See kim1990cube and van2000weak, Theorem 3.2.5.} For the Grenander estimator, kulikov2006behavior suggest a solution to this problem based on the Koml\'{o}s-Major-Tusn\'{a}dy strong approximation. In our case, we would need to rely on the strong approximation to the generic empirical process, see koltchinskii1994komlos, which leads to suboptimal results due to slower convergence rates and uniform boundedness restrictions. Our proof does not rely on the strong approximation and is based on a more standard partitioning argument.

Theorem (ref) shows that for "slow" boundary corrections with $a\in(0,1/3)$, the asymptotic distribution coincides with the one at the interior point, cf. wright1981asymptotic. However, since the estimator is centered at $m(cn^{-a})$, it is inconsistent for $m(0)$ because the "bias" $n^{1/3}(m(cn^{-a}) - m(0))$ does not vanish for all $a\in(0,1/3)$.

The choice $a=1/3$ leads to the cube-root consistent estimator similar to the tuning-free isotonic regression estimator at the interior point

For the interior point $x\in(0,1)$, wright1981asymptotic shows that

This suggests that the boundary correction with $a=1/3$ and $c=1$ should behave similarly to the tuning-free isotonic regression estimator at the interior point. Alternatively, one can try to improve the point estimates with a data-driven MSE-optimal choice which we investigate in Section (ref).

For the "fast" boundary corrections with $a\in(1/3,1)$, the convergence rate is slower than $O_P(n^{-1/3})$, which comes with a benefit of relaxing the smoothness requirement for $m$. For instance, for $a=1/2$, we only need the $\gamma$-H\"{o}lder smoothness with $\gamma>1/2$ and do not require that $m'$ exists, in which case, the convergence rate is $O_P(n^{-1/4})$. Due to the slower than the cube-root convergence rate, we do not recommend the "fast" boundary corrections for the point estimation, unless one expects that $m$ is not differentiable in the neighborhood of zero. However, the "fast" corrections are useful for inference, because as we shall show in the following section, the bootstrap works without additional subsampling or nonparametric smoothing.

Trimmed wild bootstrap

It is well-known that the bootstrap fails for various isotonic estimators at the interior point.\footnote{See kosorok2008bootstrapping and sen2010inconsistency for a formal statement and delgado2001subsampling, leger2006bootstrap, an abrevaya2005bootstrap for earlier evidences.} Several resampling schemes are available in the literature, including the smoothed nonparametric or $m$-out-of-$n$ bootstrap, see sen2010inconsistency and patra2018consistent; reshaping the objective function, see cattaneo2017bootstrap; and smoothed residual bootstrap, see guntuboyina2018nonparametric. Interestingly, as we shall show below, for the boundary point, the "fast" boundary corrections restore the consistency of the bootstrap. Consequently, we focus on more conventional bootstrap inferences.

The wild bootstrap, see wu1986jackknife and liu1988bootstrap, is arguably the most natural resampling scheme for the nonparametric regression. Unlike the naive nonparametric bootstrap, the wild bootstrap imposes the structure of the nonparametric regression model in the bootstrap world, so we may expect it to work better in finite samples than resampling methods that do not incorporate such information. At the same time, unlike the residual bootstrap, it does not rule out the heteroskedasticity.

The bootstrap procedure is as follows. First, we obtain the isotonic regression estimator $\hat m$, construct the trimmed estimator\footnote{Note that the estimator $\tilde m$ depends on $c$ and $a$ and that for the sake of simplicity of presentation we suppress this dependence.}

and compute residuals $\tilde\varepsilon_i = Y_i - \tilde m(X_i)$ for $1\leq i\leq n$. Second, we construct the wild bootstrap samples as follows:

where $(\eta_i^*)_{i=1}^n$ are i.i.d.\ random variables, independent of $(Y_i,X_i)_{i=1}^n$, and such that $\ensuremath{\mathbb{E}}\eta_i^* = 0$, $\ensuremath{\mathrm{Var}}(\eta_i^*) = 1$, and $\ensuremath{\mathbb{E}}|\eta_i^*|^{2+\delta}<\infty$.

Let $\Pr^*(.) = \Pr(.|(Y_i,X_i)_{i=1}^\infty)$ denote the bootstrap probability conditionally on the data, and let $\hat m^*$ be the isotonic regression estimator computed from the bootstrapped sample $(Y_i^*,X_i)_{i=1}^n$.

Note that it follows from groeneboom1983, Corollary 2.2 that the limiting distribution in Theorem (ref) is absolutely continuous with respect to the Lebesgue measure for every $a\in(1/3,1)$. The following result establishes the consistency of the trimmed wild bootstrap procedure.

Theorem (ref) shows that the boundary correction alone restores the consistency of the bootstrap. The reason is that for "fast" corrections the quadratic term $m'(0)$ disappears from the asymptotic distribution and this is precisely the term that is not consistently estimated by the bootstrap. Note that the bootstrap is consistent for every $a\in(1/3,1)$. However, there is a trade-off between the consistency rate of the bootstrapped estimator $O_P(n^{(a-1)/2})$ and the rate $O_P(n^{(1-3a)/2})$ at which the quadratic term disappears; see the proof of Theorem (ref). This trade-off should translate into the trade-off between the rates at which the confidence intervals shrink and the rate at which coverage errors vanish. Therefore, it seems reasonable to take $a=1/2$ as a rule-of-thumb, which leads to the same $O_P(n^{-1/4})$ rate balancing the two rates.\footnote{The inference-optimal choice of $c$ is more difficult. One can still use the rule-of-thumb $c=1$. We find in Monte Carlo experiments that there is some scope for improvements over this default choice.} Note that the slower $O_P(n^{-1/4})$ convergence rate is the price to pay for making the bootstrap work. However, the slower rate also comes with a benefit since this result does not require the existence of $m'$ and relies only on the $\gamma$-H\"{o}lder continuity with smoothness index $\gamma>1/2$.

Isotonic regression discontinuity designs

In this section, we apply our results on the isotonic regression at the boundary to monotone regression discontinuity designs.

Setup

Following hahn2001identification, we focus on the potential outcomes framework

where $D\in\{0,1\}$ is a binary treatment indicator (1 if treated and 0 otherwise), $Y_1,Y_0\in\ensuremath{\mathbf{R}}$ are unobserved potential outcomes of treated and untreated units, and $Y$ is the actual observed outcome.

The causal parameter of interest is the average treatment effect at the cutoff $x_0\in\ensuremath{\mathbf{R}}$ of some running variable $X\in\ensuremath{\mathbf{R}}$, denoted

Without further assumptions, $\theta$ is not identified in the sense that it depends on the distribution of unobserved potential outcomes $(Y_0,Y_1)$; see holland1986statistics.\footnote{While regression discontinuity designs identify local effects, there exist several approaches that aim to extrapolate away from the cutoff; see bertanha2020regression, angrist2015wanna, and other references therein. Note that if the causal effect changes monotonically in the running variable, one can obtain bounds on the causal effect between the two points.}

hahn2001identification shows formally that the causal effect is identified as

while for fuzzy designs, the causal effect is identified as\footnote{Identification in regression discontinuity designs might be testable; see mccrary2008manipulation, cattaneo2019simple, bugni2018testing, canay2017approximate.}

see Proposition (ref) in the Supplementary Material for a statement of these results under weak conditions imposed on conditional means of potential outcomes.

Estimators

We are interested in estimating the average causal effect $\theta$ of a binary treatment $D$ on some outcome $Y$. For a unit $i$, with $1\leq i\leq n$, we observe $(Y_i,D_i,X_i)$. Assuming that $\theta$ is identified from the distribution of $(Y,D,X)$ according to Eq. ((ref)), the estimator of $\theta$ is obtained by plugging-in the corresponding boundary corrected isotonic regression estimators before and after the cutoff.

Isotonic regression discontinuity designs (iRDD) exploit the monotonicity of the expected outcome $m(x) = \ensuremath{\mathbb{E}}[Y|X=x]$ and the treatment assignment probability $p(x) = \Pr(D=1|X=x)$. For concreteness, assume that both functions are non-decreasing. We also assume that $X$ has compact support $[-1,1]$ and normalize the cutoff to $x_0=0$. This restriction is without loss of generality up to the monotone transformation. Let $\hat m_-$ and $\hat m_+$ be solutions to

respectively, where $I_-$ and $I_+$ are sets of indices corresponding to negative and positive values of observations of the running variable. Therefore, $\hat m_-$ and $\hat m_+$ are the isotonic regression estimators before and after the cutoff. Note that the isotonic regression is a solution to the constrained least-squares problem over the closed convex cone of monotone functions and should have the projection interpretation when the monotonicity is violated. We might expect that the deviations from monotonicity far from the cutoff should be less harmful as opposed to the deviations at the cutoff.

\paragraph{Choice of $c$ and $a$.} Following the discussion in Section (ref), we recommend setting $c=1$\footnote{We also consider the MSE-optimal choice of $c$ in Monte Carlo experiments and find that it tends to increase the MSE in small samples.} and $a=1/3$ for the point estimation and $a=1/2$ for inference. For the point estimation, this choice delivers the one-sided counterpart to the asymptotic distribution at the interior obtained with a tuning-free isotonic regression, while for inference, this choice leads to the consistent bootstrap estimates of the asymptotic distribution, balancing the accuracy of the asymptotic approximation and the convergence rate of the estimator. The choice $a=1/3$ for point estimation can also be justified under the quadratic loss function as the one that balances the rates of the bias and the variance. On the other hand, we expect that the choice $a=1/2$ for inference should balance the rate of convergence of confidence intervals and the rate of the accuracy of bootstrap approximations. The inference-optimal choice of $c$ is a more difficult problem that we leave for future research. One possible solution is to develop the loss-function based approach that pays attention to the coverage errors and confidence interval length; see calonico2020optimal or lazarus2018har for this type of analysis.\footnote{Recently, han2019berry develop Berry-Esseen inequalities that might be leveraged on to characterize the coverage error.}

Therefore, for point estimates in sharp designs we focus on the following boundary-corrected iRDD estimator

For fuzzy designs, we also need to estimate treatment assignment probabilities before and after the cutoff, denoted $\hat p_-$ and $\hat p_+$, solving

The fuzzy iRDD estimator is computed as

Large sample distribution

Put $p(x)=\Pr(D=1|X=x)$ and for a function $g:[-1,1]\to\ensuremath{\mathbf{R}}$, define $g_+=\lim_{x\downarrow 0}g(x)$ and $g_-=\lim_{x\uparrow0}g(x)$ with some abuse of notation. The following assumption imposes several mild restrictions on the distribution of the data.

Assumption (ref) is comparable to assumptions typically used in the RDD literature, e.g., see hahn2001identification, Theorem 4. Note that we are agnostic about the smoothness of the marginal density of $X$, and only assume the existence of one-sided derivatives of conditional means. The following results describes the large sample approximation to the distribution of the sharp iRDD estimator.

The proof of this result can be found in the appendix. Theorem (ref) shows that the large sample distribution of the sharp iRDD estimator can be approximated by the difference of slopes of the greatest convex minorants of two scaled independent Brownian motions plus the parabola originating from zero and running in opposite directions. To describe the large sample distribution for the fuzzy iRDD, define additionally the conditional covariance function $\rho(x)=\ensuremath{\mathbb{E}}[\varepsilon(D-p(X))|X=x]$.

Both results follow from the CLT for the boundary-corrected isotonic regression estimator obtained in Theorem (ref); see appendix for the formal proof. A consequence of Theorems (ref) and (ref) is that the boundary-corrected iRDD estimators $\hat \theta$ and $\hat\theta^F$ are consistent in probability for the causal effect parameter $\theta$ and provide valid point estimates.

Trimmed wild bootstrap

In this section, we study the trimmed wild bootstrap. The resampling procedure is as follows. First, we construct the trimmed estimator

Second, we simulate the wild bootstrap samples as follows

where $(\eta_i^*)_{i=1}^n$ are i.i.d. multipliers, independent of $(Y_i,D_i,X_i)_{i=1}^n$, and $\tilde\varepsilon_i = Y_i - \tilde m(X_i)$. We call this procedure trimming since it trims the estimator close to boundaries when we generate bootstrap samples. Trimming is needed in addition to the boundary correction of the iRDD estimator

and its bootstrapped counterpart

where $\hat m^{*}_+$ and $\hat m^{*}_-$ are isotonic estimators computed from the bootstrapped sample $(Y_i^*,D_i,X_i)_{i=1}^n$ similarly to $\hat m_-$ and $\hat m_+$.

The following results establishes the bootstrap consistency for the sharp iRDD.

Monte Carlo experiments

In this section, we study the finite-sample performance of our iRDD estimator. We simulate $5,000$ samples of size $n\in\{200,500,1000\}$ as follows:

where $\varepsilon\sim N(0,1)$ and $\varepsilon\perp \! \! \! \perp X$. In the baseline DGP, we set $m(x) = x^3+0.25x$, $\theta = 1$, $X\sim 2\times \mathrm{Beta}(2,2)-1$, and $\sigma(x)=1$ (homoskedasticity). We compute the boundary-corrected estimator using the pool adjacent violators algorithm.

Figures (ref) illustrates the finite sample distribution of the boundary-corrected iRDD estimator for samples of different sizes. The exact finite-sample distribution is centered around the population value of the parameter and concentrates around the population parameter as the sample size increases.

Table (ref) report results of more comprehensive Monte Carlo experiments for several data-generating processes and shows the exact finite-sample bias, variance, and MSE of our iRDD estimator. We consider the following variations of the baseline DGP with two different functional forms and different amount of density near the cutoff:

To see the benefits of shape restrictions in small samples, we benchmark our iRDD estimators with the unrestricted local polynomial\footnote{We rely on the implementation available in the {\vtext} R package; see calonico2014robust, calonico2019regression, and calonico2020optimal for details.} and the k-nearest neighbors estimators.\footnote{Note that the iRDD estimator with the boundary correction $\lfloor n^{1/3}\rfloor$ might resemble the k-nearest neighbors estimator with $k_n = \lfloor n^{1/3}\rfloor$ neighbors. However, such comparison is deceptive because the k-nearest neighbors estimator with $\lfloor n^{1/3}\rfloor$ neighbors is restricted to take a constant value in the entire neighborhood of the boundary, while the isotonic regression estimator is unrestricted.} Results of experiments are consistent with the asymptotic theory. The bias, the variance, and the MSE decrease with the sample size. As expected, the MSE is higher when the density near the cutoff is lower. The shape-constrained iRDD estimator outperforms the unrestricted estimators. In many cases, we achieve more than 5-10 fold reduction of the MSE. In the Supplementary Material, Section (ref), we also consider DGPs, where the regression function violates the monotonicity constraint and regressions that are steep at the cutoff. The results are still favorable towards the iRDD estimator.

We shall acknowledge that the local linear and the k-NN estimators are applicable in more general settings as they do not rely on the monotonicity assumption and we do not recommend using the iRDD estimator for non-monotone designs.\footnote{The monotonicity of regression is a testable; see chetverikov2019testing and references therein.} However, under the monotonicity, our estimator adapts well to the boundary and seems to outperforms the existing methods. We shall also note that the case of the steep regression function is difficult for all existing approaches and that it becomes difficult to distinguish between the discontinuous and the steep continuous regression functions; see bertanha2020regression for related impossibility results.

We asses the validity of the trimmed wild bootstrap approximations, In Figure (ref), we plot the exact distribution $n^{1/4}(\hat\theta - \theta)$ and the bootstrap distribution $n^{1/4}(\hat\theta^* - \hat\theta)$ for samples of size $n\in\{200,1000\}$. Both estimators are computed using $a=1/2$ and $c=1$. As we can see from panels (b) and (e), the naive wild bootstrap without trimming and boundary correction does not work. On the other hand, the trimmed wild bootstrap mimics the finite-sample distribution. In our simulations, we use Rademacher multipliers for the bootstrap in all our experiments, i.e., $\eta_i\in\{-1,1\}$ with equal probabilities.

Lastly, we provide results of more detailed Monte Carlo experiments in the Supplementary Material. These experiments feature DGPs with heteroskedasticity, steep regression function near the cutoff, and regressions that are non-monotone away from the cutoff. The results show 2-10 fold reduction in the MSE across specifications.

Empirical illustration

Do incumbents have any electoral advantage? An extensive literature, going back at least to cummings1966, aims to answer this question. Estimating the causal effect is elusive because incumbents, by definition, are more successful politicians. Using the regression discontinuity design, lee2008randomized, documents that for the U.S. Congressional elections during 1946-1998, the incumbency advantage is 7.7% of the votes share on the next election. The design is plausibly monotone, since we do not expect that candidates with a larger margin should have a smaller vote share on the next election, at least on average. The unconstrained regression estimates presented by lee2008randomized also support the monotonicity empirically.

We revise the main finding of lee2008randomized with our sharp iRDD estimator. The dataset is publicly available as a companion for the book by angrist2008mostly. We use the rule-of-thumb choices of tuning parameters, $c=1$, $a=1/3$ (point estimation), and $a=1/2$ (inference). Figure (ref) presents the isotonic regression estimates\footnote{We use the piecewise-constant interpolation, but a higher-order polynomial interpolation is another alternative that would produce visually more appealing estimates.} of the average vote share for the Democratic party at next elections as a function of the vote share margin at the previous election (left panel). There is a pronounced jump in average outcomes for Democrats who barely win the election, compared to the results for the penultimate election (right panel). We find the point estimate of 13.8% with the 95% confidence interval $[6.6\%,26.5\%]$. While we reject the hypothesis that the incumbency advantage did not exist, our confidence intervals give a wider range of estimates. Our confidence interval may be conservative if the underlying regression is two times differentiable, however, it is robust to the failure of this assumption as well as to the inference after the model selection issues.

Of course, different approaches work differently in finite samples and it is hard to have a definite comparison. lee2008randomized estimates the causal effect by fitting parametric regressions with the global fourth-degree polynomial, which might be unstable at the boundary; see gelman2019high. We, on the other hand, rely on the nonparametric boundary-corrected isotonic regression. Interestingly, the local linear estimator with the data-driven bandwidth parameter as implemented in the STATA package {\vtext} estimates the causal effect of 6.6%.

We also compute iRDD estimates using isotonic regressions without the boundary correction, and evaluating the regression function at the most extreme to the boundary observations. With this approach, we obtain point estimates of 6.6%. As a consequence of the monotonicity, the isotonic regression estimator is biased upwards to the left of the cutoff and downwards to the right of the cutoff; see the Appendix, Section (ref). Therefore, our iRDD estimators without the boundary correction provide a lower bound on the estimated causal effect, which might be of significant interest given that it is obtained completely tuning-free.

Conclusion

This paper studies the isotonic regression estimator at the boundary of its support, an object that is particularly interesting and required in the analysis of monotone regression discontinuity designs. To the best of our knowledge, our paper is the first to develop the asymptotic theory for the isotonic regression at the boundary point. Building on these results, we offer a new perspective on monotone regression discontinuity designs and contribute more broadly to the growing literature on nonparametric identification, estimation, and inference under shape restrictions.

While the present paper focuses on the isotonic regression discontinuity designs in the classical sharp and fuzzy setting, it also opens several directions for future research. First, since various functional relations in economics are expected to be monotone, the isotonic regression might be of interest in other applications beyond the regression discontinuity designs. Second, it could be interesting to investigate how the monotonicity and other shape restrictions can be incorporated in other regression discontinuity designs; see cattaneo2017 for a review. Designing the coverage-optimal choices of the boundary correction is another important point that should be addressed in the future research. Finally, in the large-sample approximations that we use for the inference, we do not assume the existence of derivatives and rely instead on the weaker one-sided H\"{o}lder continuity condition. This indicates that our results might be honest to the relevant H\"{o}lder class, but additional study is needed to confirm this conjecture.

Acknowledgements

We thank the Editor, the Associate Editor, and anonymous referees for comments that helped us to improve significantly the paper. We are also grateful to Alex Belloni, Federico Bugni, Matias Cattaneo, Xiaohong Chen, Peter Hansen, Jonathan Hill, Jia Li, Matt Masten, Andrew Patton, Andres Santos, Valentin Verdier, Jon Wellner, and other participants of various seminars and conferences for helpful comments and conversations. A special thanks goes to Kenichi Nagasawa for his insightful comments on the earlier draft. All remaining errors are ours.

\setcounter{page}{1} \setcounter{section}{0} \setcounter{equation}{0} \setcounter{table}{0} \setcounter{figure}{0}

Inconsistency at the boundary

Put

By barlow1972statistical, Theorem 1.1, $\hat m(x)$ is the left derivative of the greatest convex minorant of the cumulative sum diagram

at $t=x$; see Figure (ref). The natural estimator of $m(0)$ is $\hat m(X_{(1)})$, which corresponds to the slope of the first-segment

where $X_{(j)}$ is the $j^{th}$ order statistics and $Y_{(j)}$ is the corresponding observation of the dependent variable.

The following result shows the isotonic regression estimator is can be inconsistent at the boundary.\footnote{We are grateful to the anonymous referee for suggesting how to improve the statement of this result.} We believe that the exact asymptotic behavior of $\hat m(X_{(1)})$ is an interesting problem to study on its own; see balabdaoui2011grenander for such analysis in the case of the Grenander estimator.

The exact asymptotic behavior of $\hat m(X_{(1)})$ is not clear, but we know that it will underestimate $m(0)$ in general. To see this, note that

where the first inequality follows from $X_{(1)}\leq cn^{-1/3}$ for some $c>0$ and monotonicity and the second by the consistency of $\hat m(cn^{-1/3})$; see Theorem (ref) (ii). Therefore,

and for large $n$, $\hat m(X_{(1)})$ will be smaller than $m(0)$ with high probability. Similarly, one can show that $\hat m(X_{(n)})$ overestimates $1$.

Proofs of main results

\setcounter{page}{1} \setcounter{section}{0} \setcounter{equation}{0} \setcounter{table}{0} \setcounter{figure}{0}

Identification

In this section, we revise the known identification results for regression discontinuity designs and show that in sharp designs, the identification can be achieved under weaker one-sided continuity conditions. We also show that under the additional assumption our identifying conditions are both necessary and sufficient.

The regression discontinuity design postulates that the probability of receiving the treatment changes discontinuously at the cutoff. In the iRDD, we also assume that the expected outcome and the probability of receiving the treatment are both monotone in the running variable. We introduce several assumptions below.

Assumption (M1) ensures that all limits exist at the discontinuity point. For sharp discontinuity design, all individuals with values of the running variable exceeding the cutoff $x_0$ receive the treatment, while all individuals below the cutoff do not. In other words, $D=\ensuremath{\mathds{1}}\{X\geq x_0\}$, and, whence $x\mapsto \Pr(D=1|X=x)$ trivially satisfies (M1), (ii). Assumption (RD) is also trivially satisfied for the sharp design. (M2) requires the local responsiveness to the treatment at the cutoff. It is not necessary for the identification, but as we shall see below, it allows us to characterize in some sense both necessary and sufficient conditions. (OC) is weaker than the full continuity at the cutoff. For more general fuzzy designs, we need additionally the conditional independence assumption.

Proposition (ref) shows that for sharp designs, $\theta$ is identified for a slightly larger class of distributions than are typically discussed in the literature. It shows that the continuity at the cutoff of both conditional mean functions is not needed and that under monotonicity conditions (M1) and (M2), the one-sided continuity turns out to be both necessary and sufficient for the identification. We illustrate this point in Figure (ref). Panel (a) shows that the causal effect $\theta$ can be identified without full continuity. Panel (b) shows that it may happen that the expression in Eq. ((ref)) coincides with $\theta$, yet the two conditional mean functions do not satisfy (OC). Such counterexamples are ruled out by (M2). Inspection of the proof of the Proposition (ref) reveals that monotonicity can be easily relaxed if we assume instead that all limits in Eq. ((ref)) and Eq. ((ref)) exist, in which case we recover the result of hahn2001identification under weaker (OC) condition for the sharp design.

It is also worth mentioning that for the fuzzy design, the local monotonicity of the treatment in the running variable allows to identify the causal effect for local compliers; see hahn2001identification.

Additional technical results

The proof of Theorem (ref) is based on the argmax continuous mapping theorem, kim1990cube, one of the conditions of which is that the argmax is a uniformly tight sequence of random variables. In our setting it is sufficient to show that the argmax of

is $O_P(n^{-1/3})$ for $a\in(0,1/3)$, where $g>0$ is arbitrary small, and that the argmax of

is $O_P(n^{-a})$ for $a\in[1/3,1)$. The following lemma serves this purpose.

The following lemma shows tightness of the argmax of the bootstrapped process:

The following lemma is a conditional argmax continuous mapping theorem for bootstrapped processes.

The following result is is probabilistic statement of the fact that for monotone functions converging pointwise to a continuous limit we also have the uniform convergence.

Additional Monte Carlo experiments

In this section we report additional results of Monte Carlo experiments for a larger class of data generating processes. First, in Table (ref), we present additional results for the DGPs 1-4 under heteroskedasticity, setting $\sigma(x) = \sqrt{x+1}$.

The results are largely similar to the homoskedastic designs in Table (ref). Next, we consider whether the data-driven choice of $c$ can improve the rule-of-thumb $c=1$.

According to Proposition (ref), the MSE-optimal estimator is $\hat m^*(0) = \hat m (c^*\hat An^{-1/3}),$ where $\hat A = \left(\frac{2}{\hat m'(0)}\sqrt{\frac{\hat \sigma^2(0)}{\hat f(0)}}\right)^{2/3}$, $\hat m',\hat \sigma^2,\hat f$ are consistent estimators of $m',\sigma^2,f$, and $c^*$ minimizes $\ensuremath{\mathbb{E}}|D^L_{[0,\infty)}(W_t - t^2)(c)|^2$. kulikov2006behavior find in simulations that $c^*\approx 0.345$. To have an idea of how much the MSE-optimal estimator performs in small samples, we consider the oracle choice of $m',\sigma^2$, and $f$. Results of these Monte Carlo experiments are presented in Table (ref). We find that the MSE-optimal estimator is inferior in small samples compared to the rule-of-thumb choice $c=1$, cf. Tables (ref) and (ref). The increase in the MSE might come from the fact that asymptotic MSE might be a poor approximation to the finite-sample MSE; cf. kulikov2006behavior.

Lastly, we augment the baseline DGPs 1-4 with four additional DGPs featuring 1) non-monotone functions; 2) functions that are steep at the cutoff; 3) functions that change the slope at the cutoff. The simulations designs are as follows:

Note that in DGPs 5-6, the regression function is increasing before the cutoff and steeply decreasing after the cutoff. The DGPs 7-8 feature the regression function that is steeply increasing after the cutoff for all $x\in[0,0.5)$ and decreasing subsequently on $[0.5,1]$. These DGPs are the most difficult because they violate the monotonicity constraint and at the same time feature steep regression functions.\footnote{As explained in Section (ref), since the isotonic regression estimator is the constrained least-squares estimator, we may expect the projection interpretation under the misspecification.} We compare the performance of our iRDD estimator to the local polynomial estimator in Table (ref). The results show 2-10 fold reduction in the MSE across specifications with the most improvement achieved when the sample size is small.

Finally, we report the empirical coverage and the length of wild bootstrap confidence intervals in Table (ref). Since the inference-optimal choice of $c$ is less obvious, we report results for a range of value of $c$. We find that for $c=1$, the coverage is approximately $90\%$ for a nominal coverage $95\%$ and that the coverage gets closer to the nominal $95\%$ coverage for larger values of $c$ as a sample size increases. Designing inference-optimal data-driven methods to select $c$ is an important point left for the future research; see e.g., calonico2020optimal and lazarus2018har.

Examples of monotone discontinuity designs

In Table (ref), we collect a list of empirical papers with monotone regression discontinuity designs. We focus only on papers where the global monotonicity is economically plausible and is empirically supported. It is worth stressing that monotonicity restricts only how the average outcome changes with the running variable and that in some references monotonicity appears due to the restricted set of values of the running variable, e.g., elderly people. However, we do not include papers where we might have global piecewise monotonicity with known change points, so the scope of the empirical applicability is probably larger.