Isotonic Regression Discontinuity Designs

TL;DR

This paper investigates isotonic regression at boundary points within regression discontinuity designs, revealing inconsistency issues and proposing boundary-corrected estimators with practical bootstrap methods, demonstrated through Monte Carlo simulations and an application to U.S. House elections.

Contribution

It introduces boundary-corrected estimators for isotonic regression at boundary points and shows their advantages over traditional methods in regression discontinuity analysis.

Findings

Boundary correction improves estimator consistency.

Bootstrap methods work without subsampling for boundary estimators.

Shape restrictions enhance finite-sample performance.

Abstract

This paper studies the estimation and inference for the isotonic regression at the boundary point, an object that is particularly interesting and required in the analysis of monotone regression discontinuity designs. We show that the isotonic regression is inconsistent in this setting and derive the asymptotic distributions of boundary corrected estimators. Interestingly, the boundary corrected estimators can be bootstrapped without subsampling or additional nonparametric smoothing which is not the case for the interior point. The Monte Carlo experiments indicate that shape restrictions can improve dramatically the finite-sample performance of unrestricted estimators. Lastly, we apply the isotonic regression discontinuity designs to estimate the causal effect of incumbency in the U.S. House elections.

Tables

Table 1. Table 1 : Results of Monte Carlo experiments

		iRDD			LP			k-NN
	$n$	Bias	Var	MSE	Bias	Var	MSE	Bias	Var	MSE
DGP 1	200	0.001	0.043	0.043	0.008	0.319	0.319	0.010	0.500	0.500
	500	-0.009	0.022	0.022	-0.007	0.110	0.110	0.010	0.280	0.280
	1000	-0.008	0.013	0.013	-0.002	0.056	0.056	0.000	0.220	0.220
DGP 2	200	-0.111	0.080	0.092	0.003	0.747	0.747	0.040	0.500	0.500
	500	-0.079	0.042	0.049	-0.002	0.225	0.225	0.010	0.290	0.290
	1000	-0.063	0.027	0.031	0.001	0.102	0.102	0.000	0.220	0.220
DGP 3	200	-0.131	0.041	0.057	-0.002	0.322	0.322	0.030	0.510	0.510
	500	-0.114	0.019	0.031	-0.005	0.109	0.109	0.010	0.290	0.290
	1000	-0.098	0.011	0.020	0.001	0.056	0.056	0.010	0.220	0.220
DGP 4	200	-0.337	0.075	0.158	-0.019	0.773	0.774	0.050	0.490	0.490
	500	-0.261	0.037	0.087	-0.021	0.213	0.213	0.010	0.290	0.290
	1000	-0.213	0.021	0.055	-0.007	0.104	0.104	0.010	0.220	0.220

Table 2. Table SM.1 : MC results: heteroskedasticity

		iRDD			LP			k-NN
	$n$	Bias	Var	MSE	Bias	Var	MSE	Bias	Var	MSE
DGP 1	200	0.005	0.043	0.043	0.009	0.319	0.319	0.010	0.490	0.490
	500	-0.007	0.022	0.022	-0.008	0.111	0.111	0.010	0.280	0.280
	1000	-0.007	0.013	0.013	-0.001	0.056	0.056	0.000	0.220	0.220
DGP 2	200	-0.101	0.083	0.093	0.006	0.791	0.791	0.090	0.200	0.210
	500	-0.075	0.043	0.049	0.000	0.224	0.224	0.040	0.090	0.090
	1000	-0.061	0.027	0.031	-0.000	0.103	0.103	0.030	0.070	0.070
DGP 3	200	-0.127	0.041	0.057	-0.004	0.323	0.323	0.020	0.510	0.510
	500	-0.111	0.019	0.031	-0.005	0.109	0.109	0.020	0.280	0.280
	1000	-0.096	0.011	0.020	0.000	0.057	0.057	0.000	0.230	0.230
DGP 4	200	-0.286	0.073	0.155	-0.019	0.753	0.754	0.020	0.480	0.490
	500	-0.223	0.035	0.084	-0.019	0.213	0.213	0.000	0.280	0.280
	1000	-0.181	0.019	0.052	-0.009	0.104	0.104	0.010	0.210	0.210

Table 3. Table SM.2 : MC experiments for the iRDD estimator with the MSE-optimal boundary correction.

		Homoskedasticity			Heteroskedasticity
	$n$	Bias	Var	MSE	Bias	Var	MSE
DGP 1	200	0.163	0.070	0.096	0.172	0.070	0.100
	500	0.060	0.018	0.021	0.063	0.018	0.022
	1000	0.043	0.011	0.013	0.045	0.011	0.013
DGP 2	200	0.024	0.097	0.097	0.045	0.100	0.102
	500	0.016	0.031	0.031	0.022	0.032	0.032
	1000	0.013	0.020	0.020	0.016	0.020	0.020
DGP 3	200	0.280	0.085	0.163	0.286	0.085	0.167
	500	0.104	0.021	0.032	0.106	0.021	0.032
	1000	0.070	0.012	0.017	0.071	0.012	0.017
DGP 4	200	0.241	0.149	0.207	0.249	0.149	0.211
	500	0.082	0.041	0.048	0.085	0.042	0.049
	1000	0.050	0.025	0.027	0.052	0.025	0.028

Table 4. Table SM.3 : MC results: steep and non-monotone regressions

		iRDD			LP
	$n$	Bias	Var	MSE	Bias	Var	MSE
Homoskedasticity
DGP 5	200	-0.1340	0.0440	0.0620	0.0193	0.3366	0.3369
	500	-0.1190	0.0200	0.0340	0.0446	0.1097	0.1117
	1000	-0.1020	0.0110	0.0220	0.0427	0.0559	0.0577
DGP 6	200	-0.3080	0.0790	0.1740	0.0272	0.8498	0.8503
	500	-0.2380	0.0380	0.0950	0.0559	0.2156	0.2186
	1000	-0.1910	0.0220	0.0590	0.0678	0.1041	0.1086
DGP 7	200	-0.1319	0.0409	0.0583	-0.0170	0.3194	0.3197
	500	-0.1131	0.0185	0.0313	-0.0114	0.1116	0.1117
	1000	-0.0986	0.0104	0.0202	-0.0026	0.0535	0.0535
DGP 8	200	-0.2891	0.0701	0.1537	-0.0061	0.7805	0.7804
	500	-0.2243	0.0344	0.0847	-0.0247	0.2284	0.2290
	1000	-0.1816	0.0197	0.0527	-0.0170	0.1051	0.1054
Heteroskedasticity
DGP 5	200	-0.1300	0.0450	0.0620	0.0202	0.3369	0.3373
	500	-0.1170	0.0200	0.0330	0.0431	0.1104	0.1123
	1000	-0.1000	0.0110	0.0210	0.0422	0.0567	0.0585
DGP 6	200	-0.3020	0.0800	0.1710	0.0273	0.8629	0.8635
	500	-0.2340	0.0390	0.0930	0.0624	0.2149	0.2187
	1000	-0.1890	0.0220	0.0580	0.0682	0.1025	0.1072
DGP 7	200	-0.1280	0.0410	0.0580	-0.0176	0.3194	0.3196
	500	-0.1110	0.0190	0.0310	-0.0111	0.1092	0.1093
	1000	-0.0970	0.0110	0.0200	-0.0021	0.0534	0.0534
DGP 8	200	-0.2810	0.0720	0.1510	-0.0038	0.7862	0.7861
	500	-0.2190	0.0350	0.0830	-0.0233	0.2269	0.2274
	1000	-0.1790	0.0200	0.0520	-0.0168	0.1065	0.1068

Table 5. Table SM.4 : Coverage and length. 5000 experiments.

		$n=200$		$n=500$		$n=1000$
	$c$	Coverage	Length	Coverage	Length	Coverage	Length
DGP 1	1	0.892	1.108	0.910	0.897	0.909	0.763
	2	0.920	0.915	0.925	0.722	0.926	0.603
	3	0.928	0.851	0.938	0.647	0.929	0.536
	4	0.937	0.822	0.946	0.608	0.946	0.497
	5	0.938	0.814	0.930	0.585	0.933	0.472
DGP 2	1	0.852	1.460	0.889	1.242	0.882	1.078
	2	0.888	1.186	0.904	0.980	0.907	0.843
	3	0.897	1.049	0.903	0.854	0.912	0.736
	4	0.891	0.963	0.903	0.782	0.921	0.669
	5	0.886	0.914	0.899	0.731	0.912	0.625
DGP 3	1	0.892	1.139	0.903	0.929	0.913	0.786
	2	0.911	0.965	0.917	0.760	0.924	0.633
	3	0.915	0.897	0.925	0.693	0.932	0.570
	4	0.907	0.875	0.937	0.657	0.940	0.536
	5	0.869	0.872	0.922	0.642	0.937	0.515
DGP 4	1	0.853	1.550	0.886	1.310	0.895	1.134
	2	0.885	1.303	0.915	1.071	0.913	0.906
	3	0.899	1.186	0.913	0.962	0.926	0.808
	4	0.878	1.121	0.920	0.898	0.939	0.754
	5	0.846	1.080	0.917	0.863	0.935	0.718

Table 6. Table SM.5 : Monotone regression discontinuity designs

Study	Outcome(s)	Treatment(s)	Running variable
Lee (2008)	Votes share in next election	Incumbency	Initial votes share
Duflo et al. (2011)	Endline scores	Higher-achieving peers	Intitial attainment
Abdulkadiroǧlu et al. (2014)	Standardized test scores	Attending elite school	Admission scores
Lucas and Mbiti (2014)	Probability of graduation	Attending elite secondary school	Admission scores
Hoekstra (2009)	Earnings	Attending flagship state university	SAT score
Clark (2010)	Test scores, university enrollment	Attending selective high school	Assignment test
Kaniel and Parham (2017)	Net capital flow	Appearance in the WSJ ranking	Returns
Schmieder et al. (2012)	Unemployment duration	Unemployment benefits	Age
Card et al. (2008)	Health care utilization	Coverage under Medicare	Age
Shigeoka (2014)	Outpatient visits	Cost-sharing policy	Age
Carpenter and Dobkin (2009)	Alcohol-related mortality	Ability to drink legally	Age
Jacob and Lefgren (2004)	Academic achievements	Summer school, grade retention	Test scores
Baum-Snow and Marion (2009)	Income, property value	Tax credit program	Fraction of eligible
Buettner (2006)	Business tax rate	Fiscal equalization transfers	Tax base
Card et al. (2007)	Job finding hazard	Severance pay	Job tenure
Chiang (2009)	Medium run test scores	Sanctions threat	School performance
Ferreira (2010)	Probability to move to a new house	Ability to transfer tax benefits	Age
Lalive (2007)	Unemployment duration	Unemployment benefits	Age
Litschig and Morrison (2013)	Education, literacy, poverty	Government transfers	Size of municipality
Ludwig and Miller (2007)	Mortality, educational attainment	Head Start funding	County poverty rate
Matsudaira (2008)	Test scores	Summer school	Test scores
Chay and Greenstone (2005)	Housing prices	Regulatory status	Pollution levels
Greenstone and Gallagher (2008)	Housing prices	Superfund clean-up status	Ranking of hazard

Equations

$$Y=m(X)+\epsilon ,\mathbb{E}[\epsilon |X]=0,$$

$$\mathcal{M}[0,1]=\left\{\varphi :[0,1]\to \mathbf{R}:\varphi (x_1)\le \varphi (x_2),\forall x_1\le x_2\right\}.$$

$$\underset{\varphi \in \mathcal{M}[0,1]}{\min }\sum _{i=1}^n(Y_i-\varphi (X_i){)}^2;$$

$$\underset{{\varphi }_1\le {\varphi }_2\le \cdots \le {\varphi }_n}{\min }\sum _{i=1}^n(Y_{(i)}-{\varphi }_i{)}^2.$$

$$|m(x)-m(0)|\le C|x{|}^{\gamma }.$$

$$n^{1/3}\left(\hat{m}\left(c{n}^{-a}\right)-m(c{n}^{-a})\right)\stackrel{\mathrm{d}}{\to }{\left|\frac{4m^{\prime }(0){\sigma }^2(0)}{f(0)}\right|}^{1/3}\underset{t\in \mathbf{R}}{\mathrm{argmax}}\{W_t-t^2\}.$$

$$n^{(1-a)/2}\left(\hat{m}\left(c{n}^{-a}\right)-m(0)\right)\stackrel{\mathrm{d}}{\to }{D}_{[0,\infty )}^L\left(\sqrt{\frac{{\sigma }^2(0)}{cf(0)}}{W}_t+\frac{t^{2}c}{2}{m}^{\prime }(0)\text{mathds}{1}_{a=1/3}\right)(1),$$

$$n^{1/3}(\hat{m}(c{n}^{-1/3})-m(0))\stackrel{\mathrm{d}}{\to }{D}_{[0,\infty )}^L\left(\sqrt{\frac{{\sigma }^2(0)}{cf(0)}}{W}_t+\frac{t^{2}c}{2}{m}^{\prime }(0)\right)(1).$$

$$n^{1/3}(\hat{m}(x)-m(x))\stackrel{\mathrm{d}}{\to }{D}_{(-\infty ,\infty )}^L\left(\sqrt{\frac{{\sigma }^2(x)}{f(x)}}{W}_t+\frac{t^2}{2}{m}^{\prime }(x)\right)(1).$$

$$n^{(1-a)/2}\left(\hat{m}\left(-c{n}^{-a}\right)-m(0)\right)\stackrel{\mathrm{d}}{\to }{D}_{(-\infty ,0]}^L\left(\sqrt{\frac{{\sigma }^2(0)}{cf(0)}}{W}_t+\frac{t^{2}c}{2}{m}^{\prime }(0)\text{mathds}{1}_{a=1/3}\right)(-1),$$

$$\tilde{m}(x)=\{\begin{array}{ll}\hat{m}(x),& x\in (c{n}^{-a},1)\\ \hat{m}(c{n}^{-a}),& x\in [0,c{n}^{-a}],\end{array}$$

$$Y_i^{\ast }=\tilde{m}(X_i)+{\eta }_i^{\ast }{\tilde{\epsilon }}_i,1\le i\le n,$$

$$\left|{\mathrm{Pr}}^{\ast }\left(n^{\frac{1-a}{2}}\left({\hat{m}}^{\ast }(c{n}^{-a})-\hat{m}(c{n}^{-a})\right)\le u\right)-{\mathrm{Pr}}^{\ast }\left(n^{\frac{1-a}{2}}\left(\hat{m}(c{n}^{-a})-m(0)\right)\le u\right)\right|\stackrel{P}{\to }0.$$

$$Y=Y_{1}D+Y_0(1-D),$$

$$\theta =\mathbb{E}[Y_1-Y_0|X=x_0].$$

$$\theta =\underset{x\downarrow x_0}{\lim }\mathbb{E}[Y|X=x]-\underset{x\uparrow x_0}{\lim }\mathbb{E}[Y|X=x],$$

$$\theta =\frac{{\lim }_{x\downarrow x_0}\mathbb{E}[Y|X=x]-{\lim }_{x\uparrow x_0}\mathbb{E}[Y|X=x]}{{\lim }_{x\downarrow x_0}\mathrm{Pr}(D=1|X=x)-{\lim }_{x\uparrow x_0}\mathrm{Pr}(D=1|X=x)};$$

$$\underset{\varphi \in \mathcal{M}[-1,0)}{\min }\sum _{i\in I_{-}}(Y_i-\varphi (X_i){)}^2\text{and}\underset{\varphi \in \mathcal{M}[0,1]}{\min }\sum _{i\in I_{+}}(Y_i-\varphi (X_i){)}^2$$

$$\hat{\theta }={\hat{m}}_{+}\left(c{n}^{-1/3}\right)-{\hat{m}}_{-}\left(-c{n}^{-1/3}\right).$$

$$\underset{\pi \in \mathcal{M}[-1,0)}{\min }\sum _{i\in I_{-}}(D_i-\pi (X_i){)}^2\text{and}\underset{\pi \in \mathcal{M}[0,1]}{\min }\sum _{i\in I_{+}}(D_i-\pi (X_i){)}^2.$$

$${\hat{\theta }}^F=\frac{{\hat{m}}_{+}\left(c{n}^{-1/3}\right)-{\hat{m}}_{-}\left(-c{n}^{-1/3}\right)}{{\hat{p}}_{+}\left(c{n}^{-1/3}\right)-{\hat{p}}_{-}\left(-c{n}^{-1/3}\right)}.$$

$$n^{1/3}(\hat{\theta }-\theta )\stackrel{\mathrm{d}}{\to }{D}_{[0,\infty )}^L\left(\sqrt{\frac{{\sigma }_{+}^2}{c{f}_{+}}}{W}_t^{+}+\frac{t^{2}c}{2}{m}_{+}^{\prime }\right)(1)-D_{(-\infty ,0]}^L\left(\sqrt{\frac{{\sigma }_{-}^2}{c{f}_{-}}}{W}_t^{-}+\frac{t^{2}c}{2}{m}_{-}^{\prime }\right)(-1),$$

$$n^{1/3}({\hat{\theta }}^F-\theta )\stackrel{\mathrm{d}}{\to }\frac{1}{p_{+}-p_{-}}{\xi }_1-\frac{m_{+}-m_{-}}{(p_{+}-p_{-}{)}^2}{\xi }_2$$

$${\xi }_1$$

$${\xi }_2$$

$$\mathrm{Cov}(W_t^{+},B_s^{+})=\frac{{\rho }_{+}(t\wedge s)}{\sqrt{{\sigma }_{+}^2{p}_{+}(1-p_{+})}},\mathrm{Cov}(W_t^{-},B_s^{-})=\frac{{\rho }_{-}(t\wedge s)}{\sqrt{{\sigma }_{-}^2{p}_{-}(1-p_{-})}}.$$

$$\tilde{m}(x)=\{\begin{array}{ll}{\hat{m}}_{-}(x),& x\in (-1,-c{n}^{-1/2})\\ {\hat{m}}_{-}(-c{n}^{-1/2}),& x\in [-c{n}^{-1/2},0]\\ {\hat{m}}_{+}(c{n}^{-1/2}),& x\in (0,c{n}^{-1/2}]\\ {\hat{m}}_{+}(x),& x\in (c{n}^{-1/2},1).\end{array}$$

$$Y_i^{\ast }=\tilde{m}(X_i)+{\eta }_i^{\ast }{\tilde{\epsilon }}_i,i=1,\dots ,n,$$

$$\check{\theta }={\hat{m}}_{+}\left(c{n}^{-1/2}\right)-{\hat{m}}_{-}\left(-c{n}^{-1/2}\right)$$

$${\check{\theta }}^{\ast }={\hat{m}}_{+}^{\ast }(c{n}^{-1/2})-{\hat{m}}_{-}^{\ast }(-c{n}^{-1/2}),$$

$$\left|{\mathrm{Pr}}^{\ast }\left(n^{1/4}({\check{\theta }}^{\ast }-\check{\theta })\le u\right)-\mathrm{Pr}\left(n^{1/4}(\check{\theta }-\theta )\le u\right)\right|\stackrel{P}{\to }0,$$

Full text

Keywords: regression discontinuity designs, shape constraints, monotonicity, isotonic regression, boundary point, wild bootstrap.

JEL classification: C14, C31.

1 Introduction

Regression discontinuity designs, see Thistlethwaite and Campbell (1960), 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 Lee and Lemieux (2010); see also Imbens and Lemieux (2008), and Cattaneo et al. (2020c) for the methodological review and Imbens and Wooldridge (2009), Abadie and Cattaneo (2018) for their place in the program evaluation literature. On the methodological side, in the seminal paper Hahn et al. (2001) 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 SM.4 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.333Shape restrictions in regression discontinuity designs is a relatively unexplored area. To the best of our knowledge, Armstrong (2015) 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.444This 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 Chetverikov et al. (2018) 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 Ayer et al. (1955), and Brunk (1956). Brunk (1970) 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 Rao (1969), who derived the asymptotic distribution of the monotone density estimator, also known as the Grenander estimator; see Grenander (1956). Wright (1981) 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 Fan and Gijbels (1992). Most of the existing results for isotonic estimators at the boundary are available only for the Grenander estimator; see Woodroofe and Sun (1993), Kulikov and Lopuhaä (2006), and Balabdaoui et al. (2011). 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 Durot and Lopuhaä (2018). Anevski and Hössjer (2002) discuss the inconsistency of the isotonic regression at the discontinuity point with deterministic equally spaced covariate and homoskedasticity. However, Anevski and Hössjer (2002) 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 Kim and Pollard (1990). The difficulty stems from the fact that conventional tightness results of Kim and Pollard (1990) and van der Vaart and Wellner (1996) do not always apply. For the Grenander estimator, Kulikov and Lopuhaä (2006) suggest a solution to this problem based on the Komlós-Major-Tusnády strong approximation. In our setting, this approach entails the strong approximation to the general empirical process, see Koltchinskii (1994) and Chernozhukov et al. (2015), 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.555In contrast, the bootstrap typically fails at the interior point; see Kosorok (2008a), Sen et al. (2010), Guntuboyina and Sen (2018), and Patra et al. (2018) for the discussion of this problem and various case-specific remedies, and Cattaneo et al. (2020b) for generic solutions that apply to a class of cube-root consistent estimators.

The paper is organized as follows. In Section 2, 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 3, where we discuss isotonic sharp and fuzzy designs. In Section 4, we study the finite sample performance of the iRDD estimator with Monte Carlo experiments. Section 5 estimates the effect of incumbency using the sharp iRDD on the data of Lee (2008). Section 6 concludes. Additional results, proofs, and Monte Carlo experiments appear in the Appendix and the Supplementary Material.

2 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 3.

2.1 Estimator

We focus on the generic nonparametric regression

[TABLE]

where the conditional mean function $m(x)=\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]$

[TABLE]

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

[TABLE]

see Brunk (1970). 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

[TABLE]

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 Ayer et al. (1955). 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.

2.2 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\downarrow 0}m(x)$. The natural estimator of $m(0)$ is $\hat{m}(X_{(1)})$. Unfortunately, it follows from Theorem A.1 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 A.1 in the Appendix for more details.666We 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.

Assumption 2.1**.**

$(Y_{i},X_{i})_{i=1}^{n}$* is an i.i.d. sample of $(Y,X)$ such that (i) $\mathbb{E}[|\varepsilon|^{2+\delta}|X]\leq C<\infty$ for some $\delta>0$ and $m$ is uniformly bounded; (ii) the distribution of $X$ has a Lebesgue density $f$, uniformly bounded away from zero and infinity, and $f(0)=\lim_{x\downarrow 0}f(x)$ exists; (iii) the conditional variance $\sigma^{2}(x)=\mathrm{Var}(Y|X=x)$ is uniformly bounded and $\sigma^{2}(0)=\lim_{x\downarrow 0}\sigma^{2}(x)$ exists; (iv) $m$ is continuously differentiable in the neighborhood of zero with $m^{\prime}(0)=\lim_{x\downarrow 0}m^{\prime}(x)>0$; (v) $m\in\mathcal{M}[0,1]$.*

For a stochastic process $\{Z_{t}:t\in A\subset\mathbf{R}\}$, let $D^{L}_{A}(Z_{t})(s)$ denote the left derivative of the greatest convex minorant at $s\in A$.777The greatest convex minorant of a function $g:A\to\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ölder continuous in the neighborhood of zero if there exists a constant $C<\infty$ such that for all $x>0$ in this neighborhood

[TABLE]

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

Theorem 2.1**.**

Suppose that Assumption 2.1 is satisfied and that $c>0$. Then

1. (i)

for $a\in(0,1/3)$

[TABLE] 2. (ii)

for $a\in[1/3,1)$

[TABLE]

where $(W_{t})_{t\in\mathbf{R}}$ is a two-sided Brownian motion,888The two-sided Brownian motion is defined as two independent Brownian motions originating from zero and moving in the opposite directions. and for $a\in(1/3,1)$, we can replace Assumption 2.1 (iv) by the $\gamma$-Höder continuity with $\gamma>(1-a)/2a$.

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.999See Kim and Pollard (1990) and van der Vaart and Wellner (1996), Theorem 3.2.5. For the Grenander estimator, Kulikov and Lopuhaä (2006) suggest a solution to this problem based on the Komlós-Major-Tusnády strong approximation. In our case, we would need to rely on the strong approximation to the generic empirical process, see Koltchinskii (1994), 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 2.1 shows that for ”slow” boundary corrections with $a\in(0,1/3)$, the asymptotic distribution coincides with the one at the interior point, cf. Wright (1981). 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

[TABLE]

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

[TABLE]

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 4.

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ölder smoothness with $\gamma>1/2$ and do not require that $m^{\prime}$ 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.

Remark 2.1**.**

One can show that for a non-decreasing function $m:[-1,0]\to\mathbf{R}$ and $a\in(0,1/3)$ the asymptotic distribution of $\hat{m}(-cn^{-a})$ is the same as in Theorem 3.1, while for $a\in[1/3,1)$

[TABLE]

where $f(0),m(0)$, and $\sigma^{2}(0)$ are defined as the left limits.

2.3 Trimmed wild bootstrap

It is well-known that the bootstrap fails for various isotonic estimators at the interior point.101010See Kosorok (2008a) and Sen et al. (2010) for a formal statement and Delgado et al. (2001), Léger and Macgibbon (2006), an Abrevaya and Huang (2005) for earlier evidences. Several resampling schemes are available in the literature, including the smoothed nonparametric or $m$-out-of-$n$ bootstrap, see Sen et al. (2010) and Patra et al. (2018); reshaping the objective function, see Cattaneo et al. (2020b); and smoothed residual bootstrap, see Guntuboyina and Sen (2018). 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 Wu (1986) and Liu (1988), 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 estimator111111Note that the estimator $\tilde{m}$ depends on $c$ and $a$ and that for the sake of simplicity of presentation we suppress this dependence.

[TABLE]

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:

[TABLE]

where $(\eta_{i}^{*})_{i=1}^{n}$ are i.i.d. random variables, independent of $(Y_{i},X_{i})_{i=1}^{n}$, and such that $\mathbb{E}\eta_{i}^{*}=0$, $\mathrm{Var}(\eta_{i}^{*})=1$, and $\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 Groeneboom (1983), Corollary 2.2 that the limiting distribution in Theorem 2.1 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 2.2**.**

Suppose that Assumption 2.1 is satisfied with (i) strengthened to (i’) $\mathbb{E}[\varepsilon^{4}|X]\leq C$ for some $C<\infty$; and (iv) replaced by (iv’) $m$ is $\gamma$-Hölder continuous in the neighborhood of zero with $\gamma>(1-a)/2a$. Then for every $u\in\mathbf{R}$ and $a\in(1/3,1)$

[TABLE]

Theorem 2.2 shows that the boundary correction alone restores the consistency of the bootstrap. The reason is that for ”fast” corrections the quadratic term $m^{\prime}(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 2.1. 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.121212The 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^{\prime}$ and relies only on the $\gamma$-Hölder continuity with smoothness index $\gamma>1/2$.

Remark 2.2**.**

It is also possible to show that for the non-decreasing function $m:[-1,0]\to\mathbf{R}$ and $a\in(1/3,1)$, the bootstrap is consistent in probability.

3 Isotonic regression discontinuity designs

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

3.1 Setup

Following Hahn et al. (2001), we focus on the potential outcomes framework

[TABLE]

where $D\in\{0,1\}$ is a binary treatment indicator (1 if treated and 0 otherwise), $Y_{1},Y_{0}\in\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\mathbf{R}$ of some running variable $X\in\mathbf{R}$, denoted

[TABLE]

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 Holland (1986).131313While regression discontinuity designs identify local effects, there exist several approaches that aim to extrapolate away from the cutoff; see Bertanha (2020), Angrist and Rokkanen (2015), 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.

Hahn et al. (2001) shows formally that the causal effect is identified as

[TABLE]

while for fuzzy designs, the causal effect is identified as141414Identification in regression discontinuity designs might be testable; see McCrary (2008), Cattaneo et al. (2020a), Bugni and Canay (2020), Canay and Kamat (2018).

[TABLE]

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

3.2 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. (1), 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)=\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

[TABLE]

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.

Choice of $c$ and $a$.

Following the discussion in Section 2, we recommend setting $c=1$151515We 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 Calonico et al. (2020) or Lazarus et al. (2018) for this type of analysis.161616Recently, Han and Kato (2019) 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

[TABLE]

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

[TABLE]

The fuzzy iRDD estimator is computed as

[TABLE]

3.3 Large sample distribution

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

Assumption 3.1**.**

$(Y_{i},D_{i},X_{i})_{i=1}^{n}$* is an i.i.d. sample of $(Y,D,X)$ such that (i) $\mathbb{E}[|\varepsilon|^{2+\delta}|X]\leq C<\infty$ for some $\delta>0$ and $m$ is uniformly bounded; (ii) the distribution of $X$ has Lebesgue density $f$, uniformly bounded away from zero and infinity on the support of $X$, and such that $f_{-}$ and $f_{+}$ exist; (iii) $\sigma^{2}$ is uniformly bounded on $[-1,1]$ and $\sigma^{2}_{+}$ and $\sigma^{2}_{-}$ exist; (iv) $m$ is continuously differentiable in the right and left neighborhoods of zero with $m^{\prime}_{-},m^{\prime}_{+}>0$; (v) $m\in\mathcal{M}[-1,1]$.*

Assumption 3.1 is comparable to assumptions typically used in the RDD literature, e.g., see Hahn et al. (2001), 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.

Theorem 3.1**.**

Suppose that Assumption 3.1 is satisfied. Then

[TABLE]

where $W_{t}^{+}$ and $W_{t}^{-}$ are two independent standard Brownian motions originating from zero and running in opposite directions.

The proof of this result can be found in the appendix. Theorem 3.1 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)=\mathbb{E}[\varepsilon(D-p(X))|X=x]$.

Theorem 3.2**.**

Suppose that Assumption 3.1 is satisfied. Suppose additionally that $p\in\mathcal{M}[-1,1]$ is continuously differentiable in the right and left neighborhoods of zero with $p^{\prime}_{-},p^{\prime}_{+}>0$, and $p_{+},p_{-}\in(0,1)$. Then

[TABLE]

with

[TABLE]

where $W_{t}^{+},W_{t}^{-},B^{+}_{t}$, and $B^{-}_{t}$ are standard Brownian motions such that any two processes with different signs are independent, and

[TABLE]

Both results follow from the CLT for the boundary-corrected isotonic regression estimator obtained in Theorem 2.1; see appendix for the formal proof. A consequence of Theorems 3.1 and 3.2 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.

3.4 Trimmed wild bootstrap

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

[TABLE]

Second, we simulate the wild bootstrap samples as follows

[TABLE]

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

[TABLE]

and its bootstrapped counterpart

[TABLE]

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.

Theorem 3.3**.**

Suppose that Assumptions 3.1 are satisfied with (i) strengthened to (i’) $\mathbb{E}[\varepsilon^{4}|X]\leq C$ for some $C<\infty$; and (iv) replaced by (iv’) $m$ is $\gamma$-Hölder continuous in the left and the right neighborhoods of zero with $\gamma>1/2$ and continuous on $[-1,0)$ and $(0,1]$. If multipliers $(\eta_{i}^{*})_{i=1}^{n}$ are such that $\mathbb{E}\eta_{i}^{*}=0$, $\mathrm{Var}(\eta_{i}^{*})=1$, and $\mathbb{E}|\eta_{i}^{*}|^{2+\delta}<\infty$ for some $\delta>0$, then for every $u\in\mathbf{R}$

[TABLE]

where $\Pr^{*}(.)=\Pr(.|(X_{i},Y_{i})_{i=1}^{\infty})$.

4 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:

[TABLE]

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 1 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 1 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:

DGP 1 sets $X\sim 2\times\mathrm{Beta}(2,2)-1$ and $m(x)=\exp(0.25x)$; 2. 2.

DGP 2 sets $X\sim 2\times\mathrm{Beta}(0.5,0.5)-1$ (low density) and $m(x)=\exp(0.25x)$; 3. 3.

DGP 3 sets $X\sim 2\times\mathrm{Beta}(2,2)-1$ and $m(x)=x^{3}+0.25x$; 4. 4.

DGP 4 sets $X\sim 2\times\mathrm{Beta}(0.5,0.5)-1$ (low density) and $m(x)=x^{3}+0.25x$;

To see the benefits of shape restrictions in small samples, we benchmark our iRDD estimators with the unrestricted local polynomial171717We rely on the implementation available in the \vtext R package; see Calonico et al. (2014), Calonico et al. (2019), and Calonico et al. (2020) for details. and the k-nearest neighbors estimators.181818Note 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 SM.3, 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.191919The monotonicity of regression is a testable; see Chetverikov (2019) 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 Bertanha (2020) for related impossibility results.

We asses the validity of the trimmed wild bootstrap approximations, In Figure 2, 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.

5 Empirical illustration

Do incumbents have any electoral advantage? An extensive literature, going back at least to Cummings Jr. (1967), aims to answer this question. Estimating the causal effect is elusive because incumbents, by definition, are more successful politicians. Using the regression discontinuity design, Lee (2008), 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 Lee (2008) also support the monotonicity empirically.

We revise the main finding of Lee (2008) with our sharp iRDD estimator. The dataset is publicly available as a companion for the book by Angrist and Pischke (2008). We use the rule-of-thumb choices of tuning parameters, $c=1$, $a=1/3$ (point estimation), and $a=1/2$ (inference). Figure 3 presents the isotonic regression estimates202020We 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. Lee (2008) estimates the causal effect by fitting parametric regressions with the global fourth-degree polynomial, which might be unstable at the boundary; see Gelman and Imbens (2019). 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 A.1. 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.

6 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 Escanciano and Cattaneo (2017) 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ölder continuity condition. This indicates that our results might be honest to the relevant Hö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.

A.1 Inconsistency at the boundary

Put

[TABLE]

By Barlow et al. (1972), Theorem 1.1, $\hat{m}(x)$ is the left derivative of the greatest convex minorant of the cumulative sum diagram

[TABLE]

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

[TABLE]

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.212121We 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 Balabdaoui et al. (2011) for such analysis in the case of the Grenander estimator.

Theorem A.1**.**

Suppose that the conditional distribution $\varepsilon|X=0$ is not degenerate and that $x\mapsto\Pr(Y\leq y|X=x)$ is continuous for every $y$. Then there exists $\epsilon>0$ such that

[TABLE]

Proof.

Let $F_{\varepsilon|X}$ be the conditional CDF of $\varepsilon|X$. Note that $\int_{-\infty}^{0}e\mathrm{d}F_{\varepsilon|X}(e|0)<0$ since $\mathbb{E}[\varepsilon|X]=0$ and the distribution $\varepsilon|X=0$ is not degenerate. Then, by continuity $\exists\epsilon>0$ such that $\int_{-\infty}^{-\epsilon}e\mathrm{d}F_{\varepsilon|X}(e|0)<0$. This shows that $F_{\varepsilon|X}(-\epsilon|0)>0$. Next, for any $\epsilon>0$

[TABLE]

where we use the fact that $X_{(1)}\xrightarrow{d}0$. Therefore, there exists $\epsilon>0$ such that

[TABLE]

∎

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

[TABLE]

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 2.1 (ii). Therefore,

[TABLE]

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$.

A.2 Proofs of main results

Proof of Theorem 2.1.

By Barlow et al. (1972), Theorem 1.1, $\hat{m}(x)$ is the left derivative of the greatest convex minorant of the cumulative sum diagram

[TABLE]

at $t=x$. This corresponds to the piecewise-constant left-continuous interpolation. Put

[TABLE]

Then for any222222For a continuous function $\Phi:[0,1]\to\mathbf{R}$, we define $\operatorname*{arg\,max}_{t\in[0,1]}\Phi(t)=\max\{s\in[0,1]:\;\Phi(s)=\max_{t\in[0,1]}\Phi(t)\}$ to accomodate non-unique maximizers. Recall that continuous function attains its maximum on compact intervals and its argmax is a closed set with a well-defined maximal element. $x\in(0,1)$ and $a\in\mathbf{R}$

[TABLE]

as can be seen from Figure A.1, see also Groeneboom and Jongbloed (2014), Lemma 3.2.

Case (i): $a\in(0,1/3)$.

For every $u\in\mathbf{R}$

[TABLE]

where the second equality follows by the switching relation in Eq. (A.1) and the last by the change of variables $s\mapsto tn^{-1/3}+cn^{-a}$.

The location of the argmax is the same as the location of the argmax of the following process

[TABLE]

due to scale and shift invariance with

[TABLE]

where

[TABLE]

We will show that the process $Z_{n1}$ converges weakly to a non-degenerate Gaussian process in $\ell^{\infty}[-K,K]$ for every $K<\infty$.

Under Assumption 2.1 (ii)-(iii) the covariance structure of the process $I_{n1}$ converges pointwise to the one of the two-sided scaled Brownian motion (two independent Brownian motions starting from zero and running in the opposite directions). Indeed, when $s,t\geq 0$

[TABLE]

where we use the mean-value theorem for some $\xi_{n}$ between $cn^{-a}$ and $cn^{-a}+n^{-1/3}(t\wedge s)$. Similarly, it can be shown that

[TABLE]

The class $\mathcal{G}_{n1}$ is VC subgraph with VC index 2 and the envelop

[TABLE]

which is square integrable

[TABLE]

Next, we verify the Lindeberg’s condition under Assumption 2.1 (i)

[TABLE]

Lastly, under Assumption 2.1 (iii), for every $\delta_{n}\to 0$

[TABLE]

Therefore, by van der Vaart and Wellner (1996), Theorem 2.11.22

[TABLE]

Under Assumption 2.1 (ii) and (iv) by Taylor’s theorem

[TABLE]

uniformly over $[-K,K]$. Lastly, by the uniform law of the iterated logarithm

[TABLE]

uniformly over $t\in[-K,K]$. Therefore, for every $K<\infty$

[TABLE]

Next, we verify conditions of the argmax continuous mapping theorem Kim and Pollard (1990), Theorem 2.7. First, note that since

[TABLE]

by Kim and Pollard (1990), Lemma 2.6, the process $t\mapsto Z_{1}(t)$ achieves its maximum a.s. at a unique point. Second, by law of iterated logarithm for the Brownian motion

[TABLE]

which shows that the quadratic term dominates asymptotically, i.e., $Z_{1}(t)\to-\infty$ as $|t|\to\infty$. It follows that the maximizer of $t\mapsto Z_{1}(t)$ is tight. Lastly, by Lemma SM.2.1 in the Supplementary Material, the argmax of $t\mapsto Z_{n1}(t)$ is uniformly tight. Therefore, by the argmax continuous mapping theorem, see Kim and Pollard (1990), Theorem 2.7

[TABLE]

By the change of variables $t\mapsto\left(\frac{a}{b}\right)^{2/3}s+\frac{c}{2b}$, scale invariance of the Brownian motion $W_{\sigma^{2}t-\mu}=\sigma W_{t}-W_{\mu}$, and scale and shift invariance of the argmax

[TABLE]

This allows us to simplify the limiting distribution as

[TABLE]

where we the use symmetry of the objective function.

Case (ii): $a\in[1/3,1)$.

For every $u\in\mathbf{R}$

[TABLE]

where the second equality follows by the switching relation in Eq. (A.1) and the last by the change of variables $s\mapsto cn^{-a}t$.

The location of the argmax is the same as the location of the argmax of the following process

[TABLE]

with

[TABLE]

where

[TABLE]

We will show that the process $Z_{n2}$ converges weakly to a non-degenerate Gaussian process in $\ell^{\infty}[0,K]$ for every $K<\infty$.

Under Assumption 2.1 (ii)-(iii) the covariance structure of the process $I_{n2}$ converges pointwise to the one of the scaled Brownian motion

[TABLE]

The class $\mathcal{G}_{n2}$ is VC subgraph with VC index 2 and envelop

[TABLE]

which is square integrable

[TABLE]

Next, we verify the Lindeberg’s condition under Assumption 2.1 (i)

[TABLE]

Lastly, under Assumption 2.1 (iii), for every $\delta_{n}\to 0$

[TABLE]

Therefore, by van der Vaart and Wellner (1996), Theorem 2.11.22

[TABLE]

Next,

[TABLE]

For $a=1/3$, under Assumption 2.1 (iv), by Taylor’s theorem uniformly over $t$ on compact sets

[TABLE]

while for $a\in(1/3,1)$ under $\gamma$-Hölder continuity of $m$ since $\gamma>(1-a)/2a$

[TABLE]

Next, by the maximal inequality Kim and Pollard (1990), p.199,

[TABLE]

where $J(1)<\infty$ is the uniform entropy integral and

[TABLE]

Since $a<1$

[TABLE]

uniformly over $[0,K]$.

Lastly,

[TABLE]

uniformly over $t\in[0,K]$. Therefore, for every $K<\infty$

[TABLE]

Next, we extend processes $Z_{n2}$ and $Z_{2}$ to the entire real line as follows

[TABLE]

We verify conditions of the argmax continuous mapping theorem Kim and Pollard (1990), Theorem 2.7. First, by the similar argument as before, the argmax of $t\mapsto\tilde{Z}_{2}(t)$ is unique and tight. Second, by Lemma SM.2.1 in the Supplementary Material, the argmax of $\tilde{Z}_{n2}$ is uniformly tight for every $u\in\mathbf{R}$ when $a=1/3$ and for every $u<0$ when $a\in(1/3,1)$. Therefore, by the argmax continuous mapping theorem, see Kim and Pollard (1990), Theorem 2.7,

[TABLE]

where the last line follows by the switching relation similar to the one in Eq. (A.1).

To conclude, it remains to show that when $a\in(1/3,1)$

[TABLE]

This follows since for every $\epsilon>0$

[TABLE]

which tends to zero as $\epsilon\downarrow 0$ as can be seen from

[TABLE]

∎

Proof of Theorem 2.2.

Put

[TABLE]

and note that now $\hat{m}^{*}(x)$ is the left derivative of the greatest convex minorant of the cumulative sum diagram

[TABLE]

at $t=x$. By the argument similar to the one used in the proof of Theorem 2.1 for every $u<0$

[TABLE]

The location of the argmax is the same as the location of the argmax of the following process

[TABLE]

with

[TABLE]

The process $I_{n}^{*}$ is the multiplier empirical process indexed by the class of functions

[TABLE]

This class is of VC subgraph type with VC index 2 and envelop $G_{n}(\epsilon,x)=n^{a/2}|\epsilon|\mathds{1}_{[0,cn^{-a}K]}(x)$, which is square-integrable

[TABLE]

This envelop satisfies Lindeberg’s condition for every $\eta>0$

[TABLE]

and for every $g_{n,t},g_{n,s}\in\mathcal{G}_{n}$ and $\delta_{n}\to 0$

[TABLE]

Next, we show that the covariance structure is

[TABLE]

with

[TABLE]

Since $\mathbb{E}[\varepsilon^{4}|X]\leq C$, the variance of $R_{n}$ tends to zero

[TABLE]

whence by Chebyshev’s inequality $R_{n}(t,s)=o_{P}(1)$. Therefore, the covariance structure converges pointwise to the one of the scaled Brownian motion

[TABLE]

where $\mathbb{E}^{*}[.]=\mathbb{E}[.|(Y_{i},X_{i})_{i=1}^{\infty}]$. By van der Vaart and Wellner (1996), Theorem 2.11.22, the class $\mathcal{G}_{n}$ is Donsker, whence by the multiplier central limit theorem, van der Vaart and Wellner (1996), Theorem 2.9.6

[TABLE]

Next, $II_{n}^{*}$ is a multiplier empirical process indexed by the degenerate class of functions

[TABLE]

Since this class is of VC subgraph type with VC index 2, by the maximal inequality

[TABLE]

where we apply Proposition SM.2.1 in the Supplementary Material.

Next, changing variables $x\mapsto cn^{-a}y$ and using the fact that $\hat{m}$ is non-decreasing

[TABLE]

where the fourth line follows by Proposition SM.2.1 in the Supplementary Material and Theorem 2.1 (ii), and the fifth since the variance of the term inside is $O(1)$.

Next

[TABLE]

in the same way we treat $III_{n2}+IV_{n2}$ in the proof of Theorem 2.1.

Combining all estimates obtained above together

[TABLE]

By Lemma SM.2.2 in the Supplementary Material, the argmax of $Z_{n}^{*}(t)$ is uniformly tight, so by Lemma SM.2.3 in the Supplementary Material

[TABLE]

To conclude, it remains to show that when $a\in(1/3,1)$

[TABLE]

This follows since for every $\epsilon>0$

[TABLE]

which tends to zero as $\epsilon\downarrow 0$; see the proof of Theorem 2.1. The result follows from Theorem 2.1 (ii) with $a>1/3$. ∎

Proof of Theorem 3.1.

Since

[TABLE]

the proof is similar to the proof of Theorem 2.1 and Remark 2.1 with $a=1/3$. Strictly speaking, the proof of Theorem 2.1 and Remark 2.1 change a little. Now $F(0)\neq 0$ and we will have $\tilde{F}(x)=F(x)-F(0)$ and $\tilde{F}_{n}(x)=F_{n}(x)-F_{n}(0)$ instead of $F(x)$ and $F_{n}(x)$ everywhere in the proof of Theorem 2.1 and $\tilde{F}(x)=F(0)-F(x)$ and $\tilde{F}_{n}(x)=F_{n}(0)-F_{n}(x)$ in the proof of Remark 2.1. Recall that if $X_{n}\perp\!\!\!\perp Y_{n},\forall n\geq 1$ and $X\perp\!\!\!\perp Y$ are such that for $X_{n}\xrightarrow{\mathrm{d}}X$ and $Y_{n}\xrightarrow{\mathrm{d}}Y$, then $(X_{n},Y_{n})\xrightarrow{\mathrm{d}}(X,Y)$. In our case, the estimators $\hat{m}_{+}$ and $\hat{m}_{-}$ are independent by the independence of the two samples. This also implies that the processes $I_{n2}^{+}$ and $I_{n2}^{-}$ (i.e. $I_{n2}$ corresponding to $\hat{m}_{+}$ and for $\hat{m}_{-}$), are asymptotically uncorrelated, which implies that $W^{+}_{t}$ and $W^{-}_{t}$ are independent due to the fact that zero-mean Gaussian processes are completely characterized by their covariance function. ∎

Proof of Theorem 3.2.

Given the proof of Theorem 3.1, it is enough to show the weak convergence of

[TABLE]

For the sake of brevity, we sketch the proof for the estimators after the cutoff. To that end, by the argument similar to the one in the proof of Theorem 2.1 for all $(u_{1},u_{2})\in\mathbf{R}^{2}$

[TABLE]

where

[TABLE]

Then

[TABLE]

where

[TABLE]

and the covariance between $\sqrt{\sigma^{2}_{+}f_{+}c}W_{t}^{+}$ and $\sqrt{p_{+}(1-p_{+})f_{+}c}B_{s}^{+}$ is

[TABLE]

Therefore, by the joint argmax continuous mapping theorem, see, e.g., Abrevaya and Huang (2005), Theorem 3,

[TABLE]

Finally, put $\hat{g}=\hat{m}_{+}\left(cn^{-1/3}\right)-\hat{m}_{-}\left(-cn^{-1/3}\right)$, $\hat{h}=\hat{p}_{+}\left(cn^{-1/3}\right)-\hat{p}_{-}\left(-cn^{-1/3}\right)$, $g=m_{+}-m_{-}$, and $h=p_{+}-p_{-}$, and note that

[TABLE]

∎

Proof of Remark 2.1.

We sketch only the most important differences below:

[TABLE]

∎

SUPPLEMENTARY MATERIAL

SM.1 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).**

The following functions are monotone in some neighborhood of $x_{0}$ (i) $x\mapsto\mathbb{E}[Y_{1}|X=x]$ and $x\mapsto\mathbb{E}[Y_{0}|X=x]$; (ii) $x\mapsto\Pr(D=1|X=x)$.

Assumption** (M2).**

$\mathbb{E}[Y_{1}|X=x_{0}]\geq\mathbb{E}[Y_{0}|X=x_{0}]$* in the non-decreasing case or $\mathbb{E}[Y_{1}|X=x_{0}]\leq\mathbb{E}[Y_{0}|X=x_{0}]$ in the non-increasing case.*

Assumption** (RD).**

Suppose that

[TABLE]

Assumption** (OC).**

Under Assumption (M1), suppose that $x\mapsto\mathbb{E}[Y_{1}|X=x]$ is right-continuous and $x\mapsto\mathbb{E}[Y_{0}|X=x]$ is left-continuous at $x_{0}$.

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=\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.

Assumption** (CI).**

Suppose that $D\perp\!\!\!\perp(Y_{1},Y_{0})|X=x$ for all $x$ in some neighborhood of $x_{0}$.

Proposition SM.1.1**.**

Under Assumptions (M1) and (OC), in the sharp design

[TABLE]

exists and equals to $\theta$. Moreover, under (M1) and (M2), if $\theta$ equals to the expression in Eq. (1), then (OC) is satisfied.

If, additionally, Assumptions (RD) and (CI) are satisfied, and instead of (OC), we assume that $x\mapsto\mathbb{E}[Y_{0}|X=x]$ and $x\mapsto\mathbb{E}[Y_{1}-Y_{0}|X=x]$ are continuous at $x_{0}$, then

[TABLE]

exists and equals to $\theta$.

Proof.

Since treatment of non-decreasing and non-increasing cases is similar, we focus only on the former. Under (M1), Rudin (1976), Theorem 4.29, ensures that all limits in Eq. (1) and Eq. (2) exist. In the sharp design

[TABLE]

where the second line follows under Assumption (OC), and the third since for any $x>0$

[TABLE]

which itself is a consequence of $Y=DY_{1}+(1-D)Y_{0}$ and $D=\mathds{1}\{X\geq x_{0}\}$. Now suppose that (M1) and (M2) are satisfied and that $\theta$ coincides with the expression in Eq. (1). Then under (M1)

[TABLE]

and

[TABLE]

Combining the two inequalities under (M2)

[TABLE]

Finally, $\theta$ is defined as the difference of inner quantities and also equals to the difference of outer quantities by assumption, which is possible only if (OC) holds, i.e,

[TABLE]

The proof for the fuzzy design is similar to the proof of Hahn et al. (2001), Theorem 2, with the only difference that monotonicity ensures existence of limits, so their Assumption (RD), (i) can be dropped. ∎

Proposition SM.1.1 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 SM.1. 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. (1) 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 SM.1.1 reveals that monotonicity can be easily relaxed if we assume instead that all limits in Eq. (1) and Eq. (2) exist, in which case we recover the result of Hahn et al. (2001) 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 Hahn et al. (2001).

SM.2 Additional technical results

The proof of Theorem 2.1 is based on the argmax continuous mapping theorem, Kim and Pollard (1990), 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

[TABLE]

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

[TABLE]

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

Lemma SM.2.1**.**

Suppose that Assumption 2.1 is satisfied. Then

• (i)

For $a\in(0,1/3)$ and $u\in\mathbf{R}$ and every $g>0$, $\operatorname*{arg\,max}_{s\in[-g,1-g]}\mathbb{M}_{n1}(s)=O_{P}(n^{-1/3})$.

• (ii)

For $a=1/3$ and $u\in\mathbf{R}$, $\operatorname*{arg\,max}_{s\in[0,1]}\mathbb{M}_{n2}(s)=O_{P}(n^{-1/3})$.

• (iii)

For $a\in(1/3,1]$ and $u<0$, $\operatorname*{arg\,max}_{s\in[0,1]}\mathbb{M}_{n2}(s)=O_{P}(n^{-a})$.

Proof.

Case (i): $a\in(0,1/3)$. Put $\mathbb{M}_{1}(s)\triangleq m(g)[F(s+g)-F(g)]-[M(s+g)-M(g)]$ with $M(s)=\int_{0}^{F(s)}m(F^{-1}(u))\mathrm{d}u$. For $s_{0}=0$, $\mathbb{M}_{1}(s_{0})=\mathbb{M}_{n1}(s_{0})=0$. By Taylor’s theorem, there exists $\xi_{1s}\in(F(s),F(g+s))$ such that

[TABLE]

where the second line follows by the mean-value theorem for some $\xi_{2s}\in(0,s)$. Then for every $s$ in the neighborhood of $s_{0}$

[TABLE]

since under Assumption 2.1 $f$ is bounded away from zero and infinity and $m^{\prime}$ is finite in the neighborhood of zero. Next we will bound the modulus of continuity for some $\delta>0$

[TABLE]

where $h_{s}\in\mathcal{H}_{\delta}=\{h_{s}(y,x)=(y-m(0))[\mathds{1}_{[0,s+g]}(x)-\mathds{1}_{[0,g]}(x)]:\;s\in[0,\delta]\}$ and $R_{n}(s)=n^{-1/3}u(F_{n}(s+g)-F_{n}(g)-F(s+g)+F(g))$. By the maximal inequality, Kim and Pollard (1990), p.199, the first term in the upper bound in Eq. (SM.1) is $\mathbb{E}\sup_{|s|\leq\delta}\left|(P_{n}-P)h_{s}\right|\leq n^{-1/2}J(1)\sqrt{PH_{\delta}^{2}}$, where $J(1)$ is the uniform entropy integral, which is finite since $\mathcal{H}_{\delta}$ is a VC-subgraph class of functions with VC index 2, $H_{\delta}(y,x)=|y-m(0)|\mathds{1}_{[g,g+\delta]}(x)$ is the envelop of $\mathcal{H}_{\delta}$, and

[TABLE]

Next, by the maximal inequality

[TABLE]

where $J(1)<\infty$ is the uniform entropy integral and $H_{\delta}(x)=\mathds{1}_{[g,g+\delta]}(x)$ is the envelop. Next, setting $\delta=1$, we get $\sup_{s\in[-g,1-g]}\left|\mathbb{M}_{n1}(s)-\mathbb{M}_{1}(s)\right|=o_{P}(1).$ Since $m(0)<m(x)$ and $f(x)>0$ for all $x\in(0,1]$, the function $s\mapsto\mathbb{M}_{1}(s)$ is strictly decreasing with a maximum achieved at $-g$, whence by van der Vaart and Wellner (1996), Corollary 3.2.3 (i), $\operatorname*{arg\,max}_{t\in[-g,1-g]}\mathbb{M}_{n1}(t)=o_{P}(1).$ Then $\phi_{n}(\delta)=\delta^{1/2}+n^{1/6}\delta$ is a good modulus of continuity function for $a=3/2$ and $r_{n}=n^{1/3}$. Indeed, for this choice $\delta\mapsto\phi_{n}(\delta)/\delta^{a}$ is decreasing and $r_{n}^{2}\phi_{n}(r_{n}^{-1})=O(n^{1/2})$. Therefore, the result follows by van der Vaart and Wellner (1996), Theorem 3.2.5.

Case (ii): $a=1/3$.

Put $\mathbb{M}_{2}(s)\triangleq m(0)F(s)-M(s)$ with $M(s)=\int_{0}^{F(s)}m(F^{-1}(u))\mathrm{d}u$. For $s_{0}=0$, $\mathbb{M}_{2}(s_{0})=\mathbb{M}_{n2}(s_{0})=0$. By Taylor’s theorem, there exists $\xi_{1s}\in(0,F(s))$ such that

[TABLE]

where the second line follows by the mean-value theorem for some $\xi_{2s}\in(0,s)$. Then

[TABLE]

Next we will bound the modulus of continuity for some $\delta>0$

[TABLE]

where $g_{s}\in\mathcal{G}_{\delta}=\{g_{s}(y,x)=(y-m(0))\mathds{1}_{[0,s]}(x):\;s\in[0,\delta]\}$ and $R_{n}(s)=n^{(a-1)/2}u(F_{n}(s)-F(s).$ By the maximal inequality Kim and Pollard (1990), p.199, the first term in the upper bound in Eq. (SM.2) is $\mathbb{E}\sup_{|s|\leq\delta}\left|(P_{n}-P)g_{s}\right|\leq n^{-1/2}J(1)\sqrt{PG_{\delta}^{2}},$ where $J(1)$ is the uniform entropy integral, which is finite since $\mathcal{G}_{\delta}$ is a VC-subgraph class of functions with VC-index 2, $G_{\delta}(y,x)=|y-m(0)|\mathds{1}_{[0,\delta]}(x)$ is the envelop of $\mathcal{G}_{\delta}$, and

[TABLE]

Next, by the maximal inequality

[TABLE]

where $J(1)<\infty$ is the uniform entropy integral and $H_{\delta}(x)=\mathds{1}_{[0,\delta]}(x)$ is the envelop. Next, setting $\delta=1$, we get $\sup_{s\in[0,1]}\left|\mathbb{M}_{n2}(s)-\mathbb{M}_{2}(s)\right|=o_{P}(1).$ Since $m(0)<m(x)$ and $f(x)>0$ for all $x\in(0,1]$, the function $s\mapsto\mathbb{M}_{2}(s)$ is strictly decreasing with maximum achieved at [math], whence by van der Vaart and Wellner (1996), Corollary 3.2.3 (i), $\operatorname*{arg\,max}_{t\in[0,1]}\mathbb{M}_{n2}(t)=o_{P}(1).$ Then the modulus of continuity is $\phi_{n}(\delta)=\delta^{1/2}+n^{a/2}\delta$. This is a good modulus of the continuity function for $\alpha=3/2$ and $r_{n}=n^{1/3}$. For this choice $\delta\mapsto\phi_{n}(\delta)/\delta^{\alpha}$ is decreasing and $r_{n}^{2}\phi_{n}(r_{n}^{-1})=O(n^{1/2})$. Therefore, the result follows by van der Vaart and Wellner (1996), Theorem 3.2.5.

Case (iii): $a\in(1/3,1]$

Here van der Vaart and Wellner (1996), Theorem 3.2.5 gives the order $O_{P}(n^{-1/3})$ only, so we will show directly using the ”peeling device” that after the change of variables

[TABLE]

Denote the process inside of the argmax as

[TABLE]

Decompose $Z_{n2}=I_{n2}+II_{n2}+III_{n2}+IV_{n2}$ similarly as in the proof of Theorem 2.1 (with $c=1$). Next, partition the set $[0,\infty)$ into intervals $S_{j}=\left\{s:\;2^{j-1}<s\leq 2^{j}\right\}$ with $j$ ranging over integers. Then if the argmax exceeds $2^{K}$, it will be located in one of the intervals $S_{j}$ with $j\geq K$ and $2^{j-1}\leq n^{a}$. Therefore, using the fact that $u<0,II_{n2}\leq 0$, and $Z_{n2}(0)=0$

[TABLE]

where the third line follows by Markov’s inequality and the fourth by the maximal inequality; cf. Theorem 2.1. The last expression can be made arbitrarily small by the choice of $K$. ∎

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

[TABLE]

Lemma SM.2.2**.**

Suppose that assumptions of Theorem 2.2 are satisfied. Then for every $a\in(1/3,1]$ and $u<0$

[TABLE]

Proof.

Decompose $\mathbb{M}_{n}^{*}=I_{n}^{*}+II_{n}^{*}+III_{n}^{*}+IV_{n}^{*}$ similarly to the proof of Theorem 2.2 (with $c=1$). Next, partition the set $[0,\infty)$ into intervals $S_{j}=\left\{s:\;2^{j-1}<s\leq 2^{j}\right\}$ with $j$ ranging over integers. Let $\|.\|_{S_{j}}$ be the supremum norm over $S_{j}$. Then if the argmax exceeds $2^{K}$, it will be located in one of the intervals $S_{j}$ with $j\geq K$ and $2^{j-1}\leq n^{a}$. Therefore, using the fact that $u<0,II_{n2}\leq 0$, and $\mathbb{M}_{n}^{*}(0)=0$

[TABLE]

where the third line follows by Markov’s inequality and the fourth by computations below. The last expression is $o_{P}(1)$ for every $K=K_{n}\to\infty$. To see that all terms above are controlled as was claimed, first note that the process $I_{n}^{*}$ is a multiplier empirical process indexed by the class of functions $\mathcal{G}_{n}=\left\{(\epsilon,x)\mapsto-n^{a/2}\epsilon\mathds{1}_{[0,n^{-a}t]}(x):\quad t\in S_{j}\right\}.$ This class is of VC subgraph type with VC index 2, whence by the maximal inequality

[TABLE]

where the second line follows since $\mathbb{E}[\varepsilon^{4}|X]\leq C$. Next, it follows from the proof of Theorem 2.2 (replacing $K$ by $2^{j}$) that $\mathbb{E}^{*}\left[\sup_{s\in S_{j}}|II_{n}^{*}(t)|\right]=o_{P}(2^{j/2})$ and that $\|III_{n}^{*}\|_{S_{j}}=o_{P}(1)$. Lastly, by the maximal inequality $\sup_{s\in S_{j}}|IV_{n}^{*}(s)-n^{a}uF(n^{-a}s)|=O_{P}(n^{(a-1)/2}2^{j/2}).$ ∎

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

Lemma SM.2.3**.**

Suppose that for every $K<\infty$

1. (i)

$\sup_{h\in BL_{1}(\ell^{\infty}[0,K])}|\mathbb{E}^{*}h(Z_{n}^{*})-\mathbb{E}h(Z)|\xrightarrow{P}0$; 2. (ii)

$\limsup_{n\to\infty}\mathrm{Pr}^{*}\left(\operatorname*{arg\,max}_{t\in[0,n^{a}]}Z_{n}^{*}(t)>K\right)=o_{P}(1),\qquad K\to\infty$; 3. (iii)

$t\mapsto Z(t)$* has unique maximizer on $[0,\infty)$, which is a tight random variable. Then*

[TABLE]

Proof.

For every $K$

[TABLE]

where by (ii)

[TABLE]

by (i) and (iii)

[TABLE]

More precisely, we used the continuous mapping theorem for the bootstrapped process Kosorok (2008b), Proposition 10.7:

[TABLE]

where the convergence is actually uniform over $z$ in arbitrary closed subset of the set of continuity points of $z\mapsto\mathrm{Pr}\left(\operatorname*{arg\,max}_{t\in[0,K]}Z(t)\geq z\right)$; see Kosorok (2008b), Lemma 10.11. ∎

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.

Proposition SM.2.1**.**

Suppose that assumptions of Theorem 2.1 are satisfied. If $m$ is continuous on $[0,1]$, then

[TABLE]

Proof.

For every $y\in[0,1]$, by Theorem 2.1, $|\tilde{m}(cn^{-a}y)-m(0)|\xrightarrow{P}0.$. Since $m$ is uniformly continuous, one can find $0\leq y_{1}\leq\dots\leq y_{p}\leq 1$ such that $|m(cn^{-a}y_{j})-m(cn^{-a}y_{j-1})|<\epsilon/2$ for all $j=2,\dots,p$. Then on the event $\{|\tilde{m}(cn^{-a}y_{j})-m(cn^{-a}y_{j})|\leq\epsilon/2,\;\forall j=1,\dots,p\}$ by the monotonicity of $\tilde{m}$, for every $x$, there exists $j=2,\dots,p$ such that

[TABLE]

whence

[TABLE]

Since $p$ is fixed, the sum of probabilities tends to zero by the pointwise consistency of $\tilde{m}$, which gives the result as $\epsilon>0$ is arbitrary. ∎

Proposition SM.2.2**.**

Suppose that Assumptions of Theorem 2.1 are satisfied. Then for arbitrary $c>0$

[TABLE]

where $A=\left(\frac{2}{m^{\prime}(0)}\sqrt{\frac{\sigma^{2}(0)}{f(0)}}\right)^{2/3}$.

Proof.

By Theorem 2.1 (ii), for $B>0$

[TABLE]

where the first equality follows by the switching relation, the second by the change of variables Brownian scaling and invariance of argmax to the scaling, the third by plugging-in corresponding value of $A$ and $B=\left|\frac{m^{\prime}(0)\sigma^{2}(0)}{2f(0)}\right|^{1/3}$, and the last by another application of the switching relation. ∎

SM.3 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 SM.1, 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 1. Next, we consider whether the data-driven choice of $c$ can improve the rule-of-thumb $c=1$.

According to Proposition SM.2.2, the MSE-optimal estimator is $\hat{m}^{*}(0)=\hat{m}(c^{*}\hat{A}n^{-1/3}),$ where $\hat{A}=\left(\frac{2}{\hat{m}^{\prime}(0)}\sqrt{\frac{\hat{\sigma}^{2}(0)}{\hat{f}(0)}}\right)^{2/3}$, $\hat{m}^{\prime},\hat{\sigma}^{2},\hat{f}$ are consistent estimators of $m^{\prime},\sigma^{2},f$, and $c^{*}$ minimizes $\mathbb{E}|D^{L}_{[0,\infty)}(W_{t}-t^{2})(c)|^{2}$. Kulikov and Lopuhaä (2006) 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^{\prime},\sigma^{2}$, and $f$. Results of these Monte Carlo experiments are presented in Table SM.2. We find that the MSE-optimal estimator is inferior in small samples compared to the rule-of-thumb choice $c=1$, cf. Tables 1 and SM.1. The increase in the MSE might come from the fact that asymptotic MSE might be a poor approximation to the finite-sample MSE; cf. Kulikov and Lopuhaä (2006).

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:

• •

DGP 5 sets $X\sim 2\times\mathrm{Beta}(2,2)-1$ and $m(x)=(1-e^{4x})\mathds{1}_{x\geq 0}+0.2x\mathds{1}_{x<0}$;

• •

DGP 6 sets $X\sim 2\times\mathrm{Beta}(0.5,0.5)-1$ and $m(x)$ is the same as in the DGP 5;

• •

DGP 7 sets $X\sim 2\times\mathrm{Beta}(2,2)-1$ and $m(x)=(-e^{-1.25}+e^{-(x-0.5)^{2}})\mathds{1}_{x\geq 0.5}-(e^{-1.25}-e^{-5(x-.5)^{2}})\mathds{1}_{0\leq x<0.5}+0.2x\mathds{1}_{x<0}$;

• •

DGP 8 sets $X\sim 2\times\mathrm{Beta}(0.5,0.5)-1$ and $m(x)$ is the same as in DGP7.

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.232323As explained in Section 3, 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 SM.3. 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 SM.4. 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., Calonico et al. (2020) and Lazarus et al. (2018).

SM.4 Examples of monotone discontinuity designs

In Table SM.5, 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.

[FIGURE:]