Optimized Regression Discontinuity Designs

TL;DR

This paper introduces a data-driven, optimization-based approach for regression discontinuity designs that directly estimates the treatment effect with guaranteed finite-sample properties, improving inference accuracy.

Contribution

It develops a convex optimization method to obtain minimax linear estimators for RDDs, addressing calibration issues of traditional non-parametric approaches.

Findings

Provides uniform confidence intervals for RDD parameters.

Extends naturally to multiple running variables.

Demonstrates improved estimation accuracy.

Abstract

The increasing popularity of regression discontinuity methods for causal inference in observational studies has led to a proliferation of different estimating strategies, most of which involve first fitting non-parametric regression models on both sides of a treatment assignment boundary and then reporting plug-in estimates for the effect of interest. In applications, however, it is often difficult to tune the non-parametric regressions in a way that is well calibrated for the specific target of inference; for example, the model with the best global in-sample fit may provide poor estimates of the discontinuity parameter. In this paper, we propose an alternative method for estimation and statistical inference in regression discontinuity designs that uses numerical convex optimization to directly obtain the finite-sample-minimax linear estimator for the regression discontinuity parameter, subject to bounds on the second derivative of the conditional response function. Given a bound on the second derivative, our proposed method is fully data-driven, and provides uniform confidence intervals for the regression discontinuity parameter with both discrete and continuous running variables. The method also naturally extends to the case of multiple running variables.

Tables

Table 1. Table 1 : Confidence interval ( α = 95 % 𝛼 percent 95 \alpha=95\% ) for the effect of raising the minimum school-leaving age on average log-earnings, as given by local linear regression with a rectangular kernel, local linear regression with a triangular kernel, and our optimized method ( 4 ). The confidence intervals account for curvature effects, provided the second derivative is bounded by B 𝐵 B .

$B$	rect. kernel	tri. kernel	optimized
0.003	$0.0213\pm 0.0761$	$0.0321\pm 0.0737$	$0.0302\pm 0.0716$
0.006	$0.0578\pm 0.0894$	$0.0497\pm 0.0867$	$0.0421\pm 0.0841$
0.012	$0.0645\pm 0.1085$	$0.0633\pm 0.1037$	$0.0557\pm 0.1003$
0.03	$0.0645\pm 0.1460$	$0.0710\pm 0.1367$	$0.0710\pm 0.1329$

Table 2. Table 2 : Summary statistics for a subset of the dataset of Matsudaira ( 2008 ) .

	fail reading		pass reading
	fail math	pass math	fail math	pass math
number of students	3,586	1,488	10,331	15,336
summer school attendance	69.6%	61%	52.9%	10.6%

Table 3. Table 3 : Estimates for the effect of summer school on math and reading scores on the following year’s test, using different estimators and choices of B 𝐵 B . Reported are bias-adjusted 95% confidence intervals, a bound on the maximum bias given our choice of B 𝐵 B , and an estimate of the sampling error conditional on { X i } subscript 𝑋 𝑖 \left\{X_{i}\right\} .

estimator:		unweighted CATE (5)			weighted CATE (6)
subject	$B$	conf. int.	m. bias	samp. err	conf. int.	m. bias	s. err
math	$0.5\times {40}^{-2}$	$0.037\pm 0.093$	0.030	0.038	$0.076\pm 0.037$	0.009	0.017
math	$1.0\times {40}^{-2}$	$0.013\pm 0.126$	0.041	0.052	$0.068\pm 0.043$	0.011	0.019
reading	$0.5\times {40}^{-2}$	$0.014\pm 0.098$	0.030	0.041	$0.044\pm 0.037$	0.009	0.017
reading	$1.0\times {40}^{-2}$	$-0.015\pm 0.130$	0.040	0.054	$0.047\pm 0.043$	0.011	0.019

Table 4. Table 4 : Estimates for the effect of television advertising on voter turnout, as well as placebo outcomes (age in years, black indicator, Democrat indicator, and female indicator). The upper half of the table displays estimates for the conditional average treatment effect at the point marked X in Figure 3 , whereas the lower half allows for a weighted treatment effect as in ( 6 ). We provide a point estimate for the treatment effect, a 95% confidence interval, as well as the worst case bias and sampling error used to produce the confidence interval.

	outcome	pt. estimate	conf. interval	max. bias	std. err.
$\tau$ at point	2008 turnout	0.130	(-0.072, 0.331)	0.071	0.079
	age	8.605	(-8.949, 26.158)	6.552	6.682
	black	-0.006	(-0.080, 0.069)	0.032	0.026
	democrat	0.019	(-0.080, 0.118)	0.035	0.039
	female	0.065	(-0.060, 0.191)	0.045	0.049
weighted $\tau$	2008 turnout	0.045	(-0.107, 0.197)	0.052	0.060
	age	3.105	(-8.702, 14.913)	4.147	4.649
	black	0.016	(-0.043, 0.075)	0.024	0.021
	democrat	-0.027	(-0.110, 0.055)	0.026	0.034
	female	0.019	(-0.083, 0.120)	0.033	0.042

Equations

$$\mu_{w}(x)=\mathbb{E}\left[Y_{i}(w)\,\big{|}\,X_{i}=x\right],\ \ \tau(x)=\mu_{1}(x)-\mu_{0}(x).$$

$$\tau(c)=\lim_{x\downarrow c}\mathbb{E}\left[Y_{i}\,\big{|}\,X_{i}=x\right]-\lim_{x\uparrow c}\mathbb{E}\left[Y_{i}\,\big{|}\,X_{i}=x\right].$$

$$\hat{\tau}=\operatorname{argmin}\left\{\frac{1}{nh_{n}}\sum_{i=1}^{n}K\left(\left\lvert\Delta_{i}\right\rvert\big{/}{h_{n}}\right)\left(Y_{i}-a-\tau W_{i}-\beta_{-}\left(\Delta_{i}\right)_{-}-\beta_{+}\left(\Delta_{i}\right)_{+}\right)^{2}\right\},$$

$$\begin{array}{rl}& \hat{\tau }=\sum _{i=1}^n{\hat{\gamma }}_i{Y}_i,\hat{\gamma }={\mathrm{argmin}}_{\gamma }\left\{\sum _{i=1}^n{\gamma }_i^2{\sigma }_i^2+I_B^2(\gamma )\right\},\\ & I_B\left(\gamma \right):=\underset{{\mu }_0(\cdot ),{\mu }_1(\cdot )}{\sup }\left\{\sum _{i=1}^n{\gamma }_i{\mu }_{W_i}(X_i)-\left({\mu }_1(c)-{\mu }_0(c)\right):\left|{\mu }_w^{\prime \prime }(x)\right|\le B\text{ for all }w,x\right\}.\end{array}$$

$${\hat{\gamma }}_c={\mathrm{argmin}}_{\gamma }\left\{\sum _{i=1}^n{\gamma }_i^2{\sigma }_i^2+{\left(\underset{\Vert {\nabla }^2{\mu }_w(x)\Vert \le B}{\sup }\left\{\sum _{i=1}^n{\gamma }_i{\mu }_{W_i}(X_i)-\left({\mu }_1(c)-{\mu }_0(c)\right)\right\}\right)}^2\right\}$$

$$\begin{array}{rl}{\hat{\gamma }}_{\ast }& ={\mathrm{argmin}}_{\gamma }\left\{\sum _{i=1}^n{\gamma }_i^2{\sigma }_i^2+{\left(\underset{\Vert {\nabla }^2{\mu }_0(x)\Vert \le B}{\sup }\left\{\sum _{i=1}^n{\gamma }_i{\mu }_0(X_i)\right\}\right)}^2:\sum _{i=1}^n{\gamma }_i{W}_i=1\right\}.\end{array}$$

$$\left\lvert\mathbb{E}\left[\hat{\tau}_{*}\,\big{|}\,\left\{X_{i}\right\}\right]-\bar{\tau}(\hat{\gamma}_{*})\right\rvert\leq\sup_{\left\lVert\nabla^{2}\mu_{0}(x)\right\rVert\leq B}\left\{\sum_{i=1}^{n}\hat{\gamma}_{*,i}\mu_{0}(X_{i})\right\},\ \ \ \bar{\tau}(\hat{\gamma}_{*}):=\sum_{i=1}^{n}W_{i}\,\hat{\gamma}_{*,i}\tau(X_{i}).$$

$$\hat{\gamma }={\mathrm{argmin}}_{\gamma }\left\{\sum _{i=1}^n{\gamma }_i^2{\sigma }_i^2+A_B^2(\gamma )\right\},$$

$$A_B(\gamma )=\underset{{\mu }_0(\cdot ),{\mu }_1(\cdot )}{\sup }\left\{\sum _{i=1}^n{\gamma }_i{\mu }_{W_i}(X_i)-\tau (c):\left|{\mu }_w(x)-{\mu }_w(c)-{\mu }_w^{\prime }(c)(x-c)\right|\le \frac{B}{2}(x-c{)}^2\right\}.$$

$$\underset{n\to \infty }{\mathrm{liminf}}\inf \left\{\mathbb{P}\left[{\mu }_1(c)-{\mu }_0(c)\in {\mathcal{I}}_{\alpha }\right]:\Vert {\nabla }^2{\mu }_w(x)\Vert \le B\text{ for all }w,x\right\}\ge 1-\alpha .$$

$$\max_{1\leq i\leq n}\left\{\hat{\gamma}_{i}^{2}\right\}\,\Big{/}\,\sum_{i=1}^{n}\hat{\gamma}_{i}^{2}\rightarrow_{p}0.$$

$$\left(\hat{\tau}-b\left(\hat{\gamma}\right)\right)\,\big{/}\,s\left(\hat{\gamma}\right)\Rightarrow\mathcal{N}\left(0,\,1\right),\ \ b(\hat{\gamma})=\sum_{i=1}^{n}\hat{\gamma}_{i}\mu_{W_{i}}(X_{i})-\tau(c),\ \ s^{2}\left(\hat{\gamma}\right):=\sum_{i=1}^{n}\hat{\gamma}_{i}^{2}\sigma_{i}^{2},$$

$$\tau (c)\in \hat{\tau }\pm l_{\alpha },l_{\alpha }=\min \left\{l:\mathbb{P}\left[\left|b+s\left(\hat{\gamma }\right)Z\right|\le l\right]\ge \alpha \text{ for all }\left|b\right|\le \hat{t}\right\},Z\sim \mathcal{N}\left(0,1\right),$$

$$\begin{split}&\tau(c)\in\mathcal{I}^{*},\ \ \mathcal{I}^{*}=\text{range}\Big{\{}\mu_{(1)}(c)-\mu_{(0)}(c):\left\lVert\nabla^{2}\mu_{(w)}(x)\right\rVert\leq B,\text{ and }\\ &\ \ \ \ \mu_{(w)}(x)=\mathbb{E}\left[Y_{i}\,\big{|}\,X_{i}=x,\,W_{i}=w\right]\text{ for all }x,\,w\in\text{supp}\left\{X_{i},\,W_{i}\right\},\Big{\}},\end{split}$$

$$\begin{array}{rl}& \underset{\gamma ,t}{\mathrm{minimize}}\sum _{i=1}^n{\gamma }_i^2{\sigma }_i^2+B^2{t}^2\text{ subject to }\\ & \sum _{i=1}^n{\gamma }_i\left(f_0(X_i)+\psi w(X_i)\left(f_1(X_i)-f_0(X_i)\right)\right)\le t,\\ & \text{for all }{f}_w(c)=0,\nabla f_w(c)=0,\Vert {\nabla }^2{f}_w(x)\Vert \le 1\text{ with }w\in \left\{0,1\right\}\\ & \sum _{i=1}^{n}w(X_i){\gamma }_i=1,\sum _{i=1}^n(1-w(X_i)){\gamma }_i=-1,\\ & \sum _{i=1}^n{\gamma }_i(X_i-c)=0,\psi \sum _{i=1}^n\left(2w(X_i)-1\right){\gamma }_i(X_i-c)=0,\end{array}$$

$$\begin{split}&\underset{f(\cdot),\,\lambda}{\text{maximize}}\ \inf_{\gamma,\,t}\Bigg{\{}\sum_{i=1}^{n}\gamma_{i}^{2}\sigma_{i}^{2}+B^{2}t^{2}\\ &\ \ \ \ \ \ +\lambda_{1}\left(\sum_{i=1}^{n}\gamma_{i}\left(f_{0}(X_{i})+\psi\,w(X_{i})\left(f_{1}(X_{i})-f_{0}(X_{i})\right)\right)-t\right)\\ &\ \ \ \ \ \ +\lambda_{2}\left(\sum_{i=1}^{n}\gamma_{i}w(X_{i})-1\right)+\lambda_{3}\left(\sum_{i=1}^{n}\gamma_{i}(1-w(X_{i}))+1\right)\\ &\ \ \ \ \ \ +\sum_{i=1}^{n}\gamma_{i}\left(\lambda_{4}+\psi\,\lambda_{5}(2w(X_{i})-1)\right)(X_{i}-c)\Bigg{\}}\\ &\text{subject to }f_{w}(c)=0,\ \nabla f_{w}(c)=0,\ \left\lVert\nabla^{2}f_{w}(x)\right\rVert\leq 1\text{ for }w\in\left\{0,\,1\right\},\\ &\ \ \ \ \ \ \lambda_{1},\geq 0,\ \lambda_{2},\,\lambda_{3}\in\mathbb{R},\ \lambda_{4},\,\lambda_{5}\in\mathbb{R}^{k},\end{split}$$

$$-2{\sigma }_i^2{\hat{\gamma }}_i={\hat{\lambda }}_1\left({\hat{f}}_0(X_i)+\psi w(X_i)\left({\hat{f}}_1(X_i)-{\hat{f}}_0(X_i)\right)\right)+{\hat{\lambda }}_{2}w(X_i)+\dots ,$$

$$\begin{array}{rl}& \underset{\tilde{f}(\cdot ),\lambda }{\mathrm{minimize}}\frac{1}{4}\sum _{i=1}^n\frac{G_i^2}{{\sigma }_i^2}+\frac{1}{4}\frac{{\lambda }_1^2}{B^2}+{\lambda }_2-{\lambda }_3\\ & \text{subject to }{G}_i=\tilde{f}(X_i)+{\lambda }_{2}w(X_i)+{\lambda }_3(1-w(X_i))\\ & +{\lambda }_4(X_i-c)+\psi {\lambda }_5(2w(X_i)-1)(X_i-c)\\ & \tilde{f}(x)={\tilde{f}}_0(x)+\psi w(x)\left({\tilde{f}}_1(x)-{\tilde{f}}_0(x)\right),{\lambda }_1\ge 0,{\lambda }_2,{\lambda }_3\in \mathbb{R},{\lambda }_4,{\lambda }_5\in {\mathbb{R}}^k,\\ & {\tilde{f}}_w(c)=0,\nabla {\tilde{f}}_w(c)=0,\Vert {\nabla }^2{\tilde{f}}_w(x)\Vert \le {\lambda }_1\text{ for }w\in \left\{0,1\right\},\end{array}$$

$$\hat{s}^{2}\left(\hat{\gamma}\right)=\sum_{i=1}^{n}\hat{\gamma}_{i}^{2}\left(Y_{i}-\hat{\mu}_{W_{i}}(X_{i})\right)^{2},\ \ \hat{s}^{2}\left(\hat{\gamma}\right)\,\big{/}\,\sum\hat{\gamma}_{i}^{2}\sigma_{i}^{2}\geq 1-o_{P}(1).$$

$$\hat{\tau}_{LATE}=\sum_{i=1}^{n}\hat{\gamma}_{i}Y_{i}\ \Big{/}\ \sum_{i=1}^{n}\hat{\gamma}_{i}W_{i},$$

$$\sum _{i=1}^n{\hat{\gamma }}_i{Z}_{ij}=0\text{ for all }j=1,...,p$$

$$\mathbb{E}\left[\hat{\tau}\,\big{|}\,X_{1},\,...,\,X_{n}\right]-\tau(c)=b\left(\hat{\gamma}\right),\ \ \operatorname{Var}\left[\hat{\tau}\,\big{|}\,X_{1},\,...,\,X_{n}\right]=s^{2}\left(\hat{\gamma}\right).$$

$$\hat{\tau}-\mathbb{E}\left[\hat{\tau}\,\big{|}\,X_{1},\,...,\,X_{n}\right]=\sum_{i=1}^{n}\hat{\gamma}_{i}\varepsilon_{i},\ \ \varepsilon_{i}=Y_{i}-\mu_{W_{i}}(X_{i}).$$

$$\displaystyle\sum_{i=1}^{n}\mathbb{E}\left[\left(\hat{\gamma}_{i}\varepsilon_{i}\right)^{q}\,\big{|}\,\hat{\gamma}_{i}\right]\,\Big{/}\,s^{q}\left(\hat{\gamma}\right)\leq\sum_{i=1}^{n}C\hat{\gamma}_{i}^{q}\,\Big{/}\,\left(\sigma_{\min}^{q}\left\lVert\hat{\gamma}\right\rVert_{2}^{q}\right)$$

$$\displaystyle\ \ \ \ \ \ \leq\frac{C}{\sigma_{\min}^{q}}\sup_{1\leq i\leq n}\left\{\left\lvert\hat{\gamma}_{i}\right\rvert\big{/}\left\lVert\hat{\gamma}\right\rVert_{2}\right\}^{q-2}\rightarrow_{p}0$$

$$\limsup_{n\rightarrow\infty}\ \mathbb{P}\left[\sum_{i=1}^{n}\mathbb{E}\left[\left(\hat{\gamma}_{i}\varepsilon_{i}\right)^{q}\,\big{|}\,\hat{\gamma}_{i}\right]\,/\,s^{q}\left(\hat{\gamma}\right)>a_{n}\right]=0.$$

$$\begin{split}&\left\lVert A\right\rVert_{\mathcal{V}}\leq\left\lVert A\right\rVert\leq\left(1-\alpha^{2}(\mathcal{V})\right)^{-1}\left\lVert A\right\rVert_{\mathcal{V}},\text{ where }\\ &\left\lVert A\right\rVert_{\mathcal{V}}:=\sup\left\{v^{\top}Av\,\big{/}\,\left\lVert v\right\rVert_{2}^{2}:v\in\mathcal{V}\right\}\text{ and }\\ &\alpha^{2}(\mathcal{V})=\sup\left\{\inf\left\{1-\left(u\cdot v\right)^{2}\,/\,\left\lVert v\right\rVert:v\in\mathcal{V}\right\}:\left\lVert u\right\rVert_{2}=1\right\};\end{split}$$

$$\mathcal{L}\left({\Vert \cdot \Vert }_{\mathcal{V}},{\lambda }_1\right)\le \mathcal{L}\left(\Vert \cdot \Vert ,{\lambda }_1\right)\le \mathcal{L}\left({\Vert \cdot \Vert }_{\mathcal{V}},\left(1-{\alpha }^2(\mathcal{V})\right){\lambda }_1\right).$$

$$\mathcal{L}\left({\Vert \cdot \Vert }_{\mathcal{V}},\left(1-{\alpha }^2(\mathcal{V})\right){\lambda }_1\right)\le \left(1-{\alpha }^2(\mathcal{V})\right)\mathcal{L}\left({\Vert \cdot \Vert }_{\mathcal{V}},{\lambda }_1\right).$$

$$\mathcal{L}\left({\Vert \cdot \Vert }_{\mathcal{V}},{\lambda }_1\right)\le \mathcal{L}\left(\Vert \cdot \Vert ,{\lambda }_1\right)\le \left(1-{\alpha }^2(\mathcal{V})\right)\mathcal{L}\left({\Vert \cdot \Vert }_{\mathcal{V}},{\lambda }_1\right),$$

$$\begin{array}{rl}& \text{minimize }\frac{1}{4}\sum _{i=1}^n{\sigma }_i^{-2}{G}_i^2+\frac{{\lambda }_1^2}{4B^2}+{\lambda }_2-{\lambda }_3\\ & \text{subject to }{G}_i={\tilde{f}}^{(h)}(N_h(X_i))+{\lambda }_{2}w(X_i)+{\lambda }_3(1-w(X_i))\\ & +{\lambda }_4(X_i-c)+\psi {\lambda }_5(2w(X_i)-1)(X_i-c)\\ & {\tilde{f}}^{(h)}(x)={\tilde{f}}_0^{(h)}(x)+\psi w(x)\left({\tilde{f}}_1^{(h)}(x)-{\tilde{f}}_0^{(h)}(x)\right)\forall x\in {\mathcal{X}}_h,{\lambda }_1\ge 0,{\lambda }_2,{\lambda }_3\in \mathbb{R},{\lambda }_4,{\lambda }_5\in {\mathbb{R}}^k,\\ & {\tilde{f}}_w^{(h)}(c)=0,{\tilde{f}}_w^{(h)}(c+h{e}_j)=0\text{ for all }j=1,...,k,w\in \left\{0,1\right\}\\ & {\Vert hv\Vert }_2^{-2}\left|{\tilde{f}}_w^{(h)}(x+hv)+{\tilde{f}}_w^{(h)}(x-hv)-2{\tilde{f}}_w^{(h)}(x)\right|\le {\lambda }_1\text{ for all }x\in {\mathcal{X}}_h,v\in \mathcal{V},w\in \left\{0,1\right\}.\end{array}$$

Full text

Keywords: Convex optimization, discrete running variable, multiple running variables, partial identification, uniform asymptotic inference.

1 Introduction

Regression discontinuity designs, first developed in the sixties (Thistlethwaite and Campbell, 1960), often allow for simple and transparent identification of treatment effects from observational data (Hahn, Todd, and Van der Klaauw, 2001; Imbens and Lemieux, 2008; Trochim, 1984), and their statistical properties have been the subject of recent interest (Armstrong and Kolesár, 2018; Calonico, Cattaneo, and Titiunik, 2014; Cheng, Fan, and Marron, 1997; Kolesár and Rothe, 2018). The sharp regression discontinuity design assumes a treatment assignment generated by a running variable $X\in\mathbb{R}^{k}$, such that individuals get treated if and only $X\in\mathcal{A}$ for some set $\mathcal{A}\subset\mathbb{R}^{k}$. For example, in epidemiology, $X\in\mathbb{R}$ could be a severity index and patients are assigned a medical intervention whenever $X\geq c$ for some threshold $c$ (i.e., $\mathcal{A}=\left\{x\in\mathbb{R}:x\geq c\right\}$); or, in political science, $X\in\mathbb{R}^{2}$ could denote the latitude and longitude of a household, and $\mathcal{A}$ could be an administrative region that has enacted a specific policy.

Given appropriate assumptions, we can identify a causal effect by comparing subjects $i$ with $X_{i}$ barely falling within the treatment region $\mathcal{A}$ to those with $X_{i}$ just outside of it. Variants of this identification strategy have proven to be useful in education (Angrist and Lavy, 1999; Black, 1999; Jacob and Lefgren, 2004), political science (Caughey and Sekhon, 2011; Keele and Titiunik, 2014; Lee, 2008), criminal justice (Berk and Rauma, 1983), program evaluation (Lalive, 2008; Ludwig and Miller, 2007), and other areas. As discussed in more detail below, standard methods for inference in the regression discontinuity design rely on plugin-in estimates from local linear regression.

In this paper, motivated by a large literature on minimax linear estimation (Armstrong and Kolesár, 2018; Cai and Low, 2003; Donoho, 1994; Donoho and Liu, 1991; Ibragimov and Khas’minskii, 1985; Johnstone, 2011; Juditsky and Nemirovski, 2009), we study an alternative approach based on directly minimizing finite sample error bounds via numerical optimization, under an assumption that the second derivative of the response surface is bounded away from the boundary of the treatment region.222Of these papers, our work is most closely related to that of Armstrong and Kolesár (2018), who explicitly consider minimax linear estimation in the regression discontinuity model for an “approximately linear” model in the sense of Sacks and Ylvisaker (1978) that places restrictions on second differences relative to the response surface at the threshold. In contrast, we assume bounded second derivatives away from the threshold. An advantage of their approach is that it allows for a closed form solution. However, a disadvantage is that they allow for discontinuities in the response surface away from the threshold, which implies that given the same value for our bound on the second derivative and their bound on second differences, our confidence intervals can be substantially shorter (moreover, allowing for discontinuities in the response surface does not seem conceptually attractive given that the assumption of continuity of the conditional expectation at the threshold is fundamental to the regression discontinuity design). We discuss this comparison further in Section 1.3. This approach has several advantages relative to local regression. Our estimator is well defined regardless of the shape of the treatment region $\mathcal{A}$, whether it be a half line as in the standard univariate regression discontinuity specification or an oddly shaped region as might appear with a geographic regression discontinuity; moreover, our inference is asymptotically valid for both discrete and continuous running variables. Finally, even with univariate designs, our approach strictly dominates local linear regression in terms of minimax mean-squared error.

For simplicity of exposition, we start by presenting our method in the context of classical univariate regression discontinuity designs with a single treatment cutoff, i.e., with $X_{i}\in\mathbb{R}$ and $\mathcal{A}=\left\{x\in\mathbb{R}:x\geq c\right\}$; a solution to the more general problem will then follow by direct extension. A software implementation, optrdd for R, is available on CRAN and from github.com/swager/optrdd.

1.1 Optimized Inference with Univariate Discontinuities

We have access to $i=1,\,...,\,n$ independent pairs $(X_{i},\,Y_{i})$ where $X_{i}\in\mathbb{R}$ is the running variable and $Y_{i}\in\mathbb{R}$ is our outcome of interest; the treatment is assigned as $W_{i}=\mathbf{1}\left(\left\{X_{i}\geq c\right\}\right)$. Following the potential outcomes model (Imbens and Rubin, 2015; Neyman, 1923; Rubin, 1974), we posit potential outcomes $Y_{i}(w)$, for $w\in\left\{0,\,1\right\}$ corresponding to the outcome subject $i$ would have experienced had they received treatment $w$, and define the conditional average treatment effect $\tau(x)$ in terms of the conditional response functions $\mu_{w}(x)$:

[TABLE]

Provided the functions $\mu_{w}(x)$ are both continuous at $c$, the regression discontinuity identifies the conditional average treatment effect at the threshold $c$,333In the fuzzy regression discontinuity design where the probability of receiving the treatment changes discontinuously at $x=c$, but not necessarily from zero to one, the estimand can be written as the ratio of two such differences. The issues we address in this paper also arise in that setting, and the present discussion extends naturally to it; see Section 5 for a discussion.

[TABLE]

Given this setup, local linear regression is a popular strategy for estimating $\tau(c)$ (Hahn et al., 2001; Porter, 2003):

[TABLE]

where $K(\cdot)$ is some weighting function, $h_{n}$ is a bandwidth, $\Delta_{i}=X_{i}-c$, and $a$ and $\beta_{\pm}$ are nuisance parameters. When we do not observe data right at the boundary $c$ (e.g., when $X_{i}$ has discrete support), then $\tau(c)$ is not point identified; however, given smoothness assumptions on $\mu_{w}(x)$, Kolesár and Rothe (2018) propose an approach to local linear regression that can still be used to construct partial identification intervals for $\tau(c)$ in the sense of Imbens and Manski (2004) (see Section 2.1 for a discussion).

The behavior of regression discontinuity estimation via local linear regression is fairly well understood. When the running variable $X$ is continuous (i.e., $X$ has a continuous positive density at $c$) and $\mu_{w}(x)$ is twice differentiable with a bounded second derivative in a neighborhood of $c$, Cheng, Fan, and Marron (1997) show that the triangular kernel $K(t)=(1-t)_{+}$ minimizes worst-case asymptotic mean-squared error among all possible choices of $K$, Imbens and Kalyanaraman (2012) provide a data-adaptive choice of $h_{n}$ to minimize the mean-squared error of the resulting estimator, and Calonico, Cattaneo, and Titiunik (2014) propose a method for removing bias effects due to the curvature of $\mu_{w}(x)$ to allow for asymptotically unbiased estimation. Meanwhile, given a second-derivative bound $\left\lvert\mu_{w}^{\prime\prime}(x)\right\rvert\leq B$, Armstrong and Kolesár (2018) and Kolesár and Rothe (2018) construct confidence intervals centered at the local linear estimator $\hat{\tau}$ that attain uniform asymptotic coverage, even when the running variable $X$ may be discrete.

Despite its ubiquity, however, local linear regression still has some shortfalls. First of all, under the bounded second derivative assumption often used to justify local linear regression (i.e., that $\mu_{w}(x)$ is twice differentiable and $\left\lvert\mu_{w}^{\prime\prime}(x)\right\rvert\leq B$ in a neighborhood of $c$), local linear regression is not the minimax optimal linear estimator for $\tau(c)$—even with a continuous running variable. Second, and perhaps even more importantly, all the motivating theory for local linear regression relies on $X$ having a continuous distribution; however, in practice, $X$ often has a discrete distribution with a modest number of points of support. When the running variable is discrete there is no compelling reason to expect local linear regression to be particularly effective in estimating the causal effect of interest.444One practical inconvenience that can arise in local linear regression with discrete running variables is that, if we use a data-driven rule to pick the bandwidth $\hat{h}$ (e.g., the one of Imbens and Kalyanaraman (2012)), we may end up with no data inside the specified range (i.e., there may be no observations with $|X_{i}-c|\leq h$); the practitioner is then forced to select a different bandwidth ad-hoc. Ideally, methods for regression discontinuity analysis should be fully data-driven, even when $X$ is discrete. In spite of these limitations, local linear regression is still the method of choice, largely because of its intuitive appeal.

As discussed above, the goal of this paper is to show that we can systematically do better. Regardless of the shape of the kernel $K(\cdot)$ in (3), local linear regression yields a linear estimator555We note the unfortunate terminological overlap between “local linear regression” estimators of $\tau$ and “linear” estimators of type $\hat{\tau}=\sum_{i=1}^{n}\hat{\gamma}_{i}Y_{i}$. The word linear in these two contexts refers to different things. All “local linear regression” estimators are “linear” in the latter sense, but not vice-versa. for $\tau$, i.e., one of the form $\hat{\tau}=\sum_{i=1}^{n}\hat{\gamma}_{i}Y_{i}$ for weights $\hat{\gamma}_{i}$ that depend only on the distances $X_{i}-c$. Here, we find that if we are willing to rely on numerical optimization tools, then minimax linear estimation of $\tau(c)$, i.e., minimax among all estimators of the form $\hat{\tau}=\sum_{i=1}^{n}\hat{\gamma}_{i}Y_{i}$ conditionally on $X_{i}$, is both simple and methodologically transparent under natural assumptions on the regularity of $\mu_{(w)}(x)$. If we know that \operatorname{Var}\left[Y_{i}\,\big{|}\,X_{i}\right]=\sigma_{i}^{2} and that $\left\lvert\mu^{\prime\prime}_{w}(x)\right\rvert\leq B$, we propose estimating $\tau(c)$ as follows:

[TABLE]

Because this estimator is minimax among the class of linear estimators, it is at least as accurate as any local linear regression estimator in a minimax sense over all problems with \operatorname{Var}\left[Y_{i}\,\big{|}\,X_{i}\right]=\sigma_{i}^{2} and $\left\lvert\mu^{\prime\prime}_{w}(x)\right\rvert\leq B$. For further discussion of related estimators, see Cai and Low (2003), Donoho (1994), Donoho and Liu (1991), and Juditsky and Nemirovski (2009).

When $X_{i}$ has a discrete distribution, the parameter $\tau(c)$ is usually not point identified because there may not be any observations $X_{i}$ in a small neighborhood of $c$. However, we can get meaningful partial identification of $\tau(c)$ thanks to our bounds on the second derivative of $\mu_{w}(x)$. Moreover, because our approach controls for bias in finite samples, the estimator (4) is still justified in the partially identified setting and, as discussed further in Section 2.1, provides valid confidence intervals for $\tau(c)$ in the sense of Imbens and Manski (2004). We view the fact that our estimator can seamlessly move between the point and partially identified settings as an important feature.

To motivate the above estimator qualitatively, note that the first term in the minimization problem corresponds to the conditional variance of $\hat{\tau}$ given $\left\{X_{i}\right\}$, while the second term is the worst-case conditional squared bias given that $|\mu^{\prime\prime}_{w}(x)|\leq B$; thus, given our regularity assumptions, the estimator $\hat{\tau}$ minimizes the worst-case conditional mean-squared error among all linear estimators. Because no constraints are placed on $\mu_{w}(c)$ or $\mu^{\prime}_{w}(c)$, the optimization in (4) also automatically enforces the constraints $\sum_{i}W_{i}\hat{\gamma}_{i}=1$, $\sum_{i}(1-W_{i})\hat{\gamma}_{i}=-1$, $\sum_{i}W_{i}(X_{i}-c)\hat{\gamma}_{i}=0$, and $\sum_{i}(1-W_{i})(X_{i}-c)\hat{\gamma}_{i}=0$.666At values of $\gamma$ for which these constraints are not all satisfied, we can choose $\mu_{1}(x)$ and $\mu_{0}(x)$ with second derivative bounded by $B$ such as to make the conditional bias arbitrarily bad, i.e., $I_{B}(\gamma)=+\infty$. Thus, the solution $\hat{\gamma}$ to (4) must satisfy the constraints. This is a convex program, and can be efficiently solved using readily available software described in, e.g., Boyd and Vandenberghe (2004).

Although this estimator depends explicitly on knowledge of $\sigma_{i}^{2}$ and $B$, we note that all practical methods for estimation in the regression discontinuity model, including Calonico et al. (2014), Imbens and Kalyanaraman (2012) and Kolesár and Rothe (2018), require estimating related quantities in order to tune the algorithm. Then, once estimators for these parameters have been specified, the procedure (4) is fully automatic—in particular, there is no need to ask whether the running variable is discrete or continuous, as the optimization is conditional on $\left\{X_{i}\right\}$—whereas the baseline procedures still have other choices to make, e.g., what weight function $K(\cdot)$ to use, or whether to debias the resulting $\smash{\hat{\tau}}$-estimator. See Sections 3 and 4 for discussions on how to choose $B$ in practice.

Figure 1 compares the weights $\hat{\gamma}_{i}$ obtained via (4) in two different settings: one with a discrete, asymmetric running variable $X$ depicted in the left panel of Figure 2, and the other with a standard Gaussian running variable. We see that, for $n=1,000$, the resulting weighting functions look fairly similar, and are also comparable to the implicit weighting function generated by local linear regression with a triangular kernel. However, as $n$ grows and the discreteness becomes more severe, our method changes both the shape and the scale of the weights, and the discrepancy between the optimized weighting schemes for discrete versus continuous running variables becomes more pronounced.

In the right panel of Figure 2, we also compare the worst-case conditional mean-squared error of our method relative to that of optimally tuned local linear regression, both with a rectangular and triangular kernel. For the smallest sample size we consider, $n=333$, the discreteness of the running variable has a fairly mild effect on estimation and—as one might have expected—the triangular kernel is noticeably better than the rectangular kernel, while our method is slightly better than the triangular kernel. However, as the sample size increases, the performance of local linear regression relative to our method ebbs and flows rather unpredictably.777As a matter of intellectual curiosity, it is intriguing to ask whether there exist discrete distributions for which the rectangular kernel may work substantially better than the triangular kernel, or whether additional algorithmic tweaks—such as using different bandwidths on different sides of the threshold—may have helped (in the above example, we used the same bandwidth for local linear regression on both sides of the boundary). However, from a practical perspective, the estimator (4) removes the need to consider such questions in applied data analysis, and automatically adapts to the structure of the data at hand.

1.2 Optimized Inference with Generic Discontinuities

The methods presented above extend naturally to the general case, where $X_{i}\in\mathbb{R}^{k}$ may be multivariate and $\mathcal{A}$ is unrestricted. The problem of regression discontinuity inference with multiple running variables is considerably richer than the corresponding problem with a single running variable, because an investigator could now plausibly hope to identify many different treatment effects along the boundary of the treated region $\mathcal{A}$. Most of the existing literature on this setup, including Papay et al. (2011), Reardon and Robinson (2012) and Wong et al. (2013), have focused on these questions of identification, while using some form of local linear regression for estimation.

In the multivariate case, however, questions about how to tune local linear regression are exacerbated, as the problems of choosing the kernel function $K(\cdot)$ and the bandwidth $h$ are now multivariate. Perhaps for this reason, it is still popular to use univariate methods to estimate treatment effects in the multivariate setting by, e.g., using shortest distance to the boundary of the treatment region $\mathcal{A}$ as a univariate running variable (Black, 1999), or only considering a subset of the data where univariate methods are appropriate (Jacob and Lefgren, 2004; Matsudaira, 2008).

Here, we show how our optimization-based method can be used to side-step the problem of choosing a multivariate kernel function by hand. In addition to providing a simple-to-apply algorithm, our method lets us explicitly account for the curvature of the mean-response function $\mu_{w}(x)$ for statistical inference, thus strengthening formal guarantees relative to prior work.

Relative to the univariate case, the multivariate case has two additional subtleties we need to address. First, in (4) it is natural to impose a constraint $|\mu_{w}^{\prime\prime}(x)|\leq B$ to ensure smoothness; in the multivariate case, however, we have more choices to make. For example, do we constrain $\mu_{w}(x)$ to be an additive function, or do we allow for interactions? Here, we opt for the more flexible specification, and simply require that $\left\lVert\nabla^{2}\mu_{w}(x)\right\rVert\leq B$, where $\left\lVert\cdot\right\rVert$ denotes the operator norm (i.e., the largest absolute eigenvalue of the second derivative).

Moreover, as emphasized by Papay et al. (2011), whereas the univariate design only enables us to identify the conditional average treatment effect at the threshold $c$, the multivariate design enables us to potentially identify a larger family of treatment effect functionals. Here, we focus on the following two causal estimands. First, writing $c$ for a focal point of interest, let $\hat{\tau}_{c}=\sum_{i=1}^{n}\hat{\gamma}_{c,i}Y_{i}$ with

[TABLE]

denote an estimator for the conditional average treatment effect at $c$. The upside of this approach is that it gives us an estimand that is easy to interpret; the downside is that, when curvature is non-negligible, (5) can effectively only make use of data near the specified focal point $c$, thus resulting in relatively long confidence intervals.

In order to potentially improve precision, we also study weighted conditional average treatment effect estimation with weights greedily chosen such as to make the inference as precise as possible: In the spirit of Crump et al. (2009), Li et al. (2017) or Robins et al. (2008), we consider $\hat{\tau}_{*}=\sum_{i=1}^{n}\hat{\gamma}_{*,i}Y_{i}$, with

[TABLE]

In other words, we seek to pick weights $\gamma_{i}$ that are nearly immune to bias due to curvature of the baseline response surface $\mu_{0}(x)$. By construction, this estimator satisfies

[TABLE]

Because $\sum W_{i}\hat{\gamma}_{*,i}=1$, we see that $\bar{\tau}(\hat{\gamma}_{*})$ is in fact a weighted average of the conditional average treatment effect function $\tau(\cdot)$ over the treated sample. If we ignored the curvature of $\tau(\cdot)$, we could interpret $\hat{\tau}_{*}$ as an estimate for the conditional average treatment effect at $x_{*}=\sum\hat{\gamma}_{*,i}W_{i}X_{i}$.

In some cases, it is helpful to consider other interpretations of the estimand underlying (6). It we are willing to assume a constant treatment effect $\tau(x)=\tau$, then $\bar{\tau}=\tau$, and $\hat{\tau}_{*}$ is the minimax linear estimator for $\tau$. Relatedly, we can always use the confidence intervals from Section 2.1 built around $\hat{\tau}_{*}$ to test the global null hypothesis $\tau(x)=0$ for all $x$.

To gain intuition for the multivariate version of our method, we outline a simple example building on the work of Keele and Titiunik (2014) on the effect of television advertising on voter turnout in presidential elections. To estimate this effect, Keele and Titiunik (2014) examine a school district in New Jersey, half of which belongs to the Philadelphia media market and the other half to the New York media market. Before the 2008 presidential elections, the Philadelphia half was subject to heavy campaign advertising whereas the New York half was not, thus creating a natural experiment for the effect of television advertising assuming the media market boundary didn’t coincide with other major boundaries within the school district. Keele and Titiunik (2014) use this identification strategy to build a regression discontinuity design, comparing sets of households straddling the media market boundary.

However, despite the multivariate identification strategy, Keele and Titiunik (2014) then reduce the problem to a univariate regression discontinuity problem for estimation: They first compute Euclidean distances $D_{i}=\left\lVert X_{i}-c\right\rVert_{2}$ to a focal point $c$, and then use $D_{i}$ as a univariate running variable. In contrast, our approach allows for transparent inference without needing to rely on such a reduction. Figure 3 depicts $\gamma$-weights generated by our optimized approach; the resulting treatment effect estimator is then $\sum_{i}\gamma_{i}Y_{i}$. Qualitatively, we replicate the “no measurable effect” finding of Keele and Titiunik (2014), while directly and uniformly controlling for spatial curvature effects. We discuss details, including placebo diagnostics and the choice of tuning parameter, in Section 4.2.

We also see that, at least here, standard heuristics used to reduce the multivariate regression discontinuity problem to a univariate one are not sharp. In the setup of the left panel of Figure 3, where we seek to estimate the treatment effect at a focal point $c$, some treated points due west of $c$ get a positive weight, whereas points the same distance south from $c$ get a mildly negative weight, thus violating the heuristic of Keele and Titiunik (2014) that weights should only depend on $D_{i}=\left\lVert X_{i}-c\right\rVert_{2}$. Meanwhile, we can compare the approach in the right panel of Figure 3, where we allow for averaging of treatment effects along the boundary, to the popular heuristic of using shortest distance to the boundary of the treatment region as a univariate running variable (Black, 1999). But this reduction again does not capture the behavior of our optimized estimator: There are some points at the eastern edge of the treated region that are very close to the boundary, but get essentially zero weight (presumably because there are no nearby units on the control side of the boundary).

1.3 Related Work

The idea of constructing estimators of the type (4) that are minimax with respect to a regularity class for the underlying data-generating process has a long history in statistics. In early work, Legostaeva and Shiryaev (1971) and Sacks and Ylvisaker (1978) independently studied inference in “almost” linear models that arise from taking a Taylor expansion around a point; see also Cheng et al. (1997). For a broader discussion of minimax linear estimation over non-parametric function classes, see Cai and Low (2003), Donoho (1994), Ibragimov and Khas’minskii (1985), Johnstone (2011), Juditsky and Nemirovski (2009), and references therein. An important result in this literature that, for many problems of interest, minimax linear estimators are within a small explicit constant of being minimax among all estimators (Donoho and Liu, 1991).

Armstrong and Kolesár (2018) apply these methods to regression discontinuity designs, resulting in an estimator of the form (4), except with weights888Armstrong and Kolesár (2018) also consider a more general setting where we assume accuracy of the $k$-th order Taylor expansion of $\mu_{w}(x)$ around $c$; and, in fact, our method also extends to this setting. Here, however, we focus on second-derivative bounds, which are by far the most common in applications.

[TABLE]

Now, although this class of functions is cosmetically quite similar to the bounded-second-derivative class used in (4), we note that the class of weights allowed for in (8) is substantially larger, even if the value of $B$ is the same. This is because the functions $\mu_{w}(\cdot)$ underlying the above weighting scheme need not be continuous, and can in fact have jumps of magnitude $B(x-c)^{2}/2$. Given that the key assumption underlying regression discontinuity designs is continuity of the conditional means of the potential outcomes at the threshold for the running variable, it would appear to be reasonable to impose continuity away from the threshold as well. Allowing for jumps through the condition (8) can make the resulting confidence intervals for $\tau(c)$ substantially larger than they are under the smoothness condition with bounded second derivatives. One key motivation for the weighting scheme (8) rather than our proposed one (4) appears to be that the optimization problem induced by (8) is substantially easier, and allows for closed-form solutions for $\hat{\gamma}_{i}$. Conversely, we are aware of no closed-form solution for (4), and instead need to rely on numeric convex optimization.

In the special case where the running variable $X$ is assumed to have a continuous density around the threshold $c$, there has been a considerable number of recent proposals for asymptotic confidence intervals while imposing smoothness assumptions on $\mu_{w}(x)$. Calonico, Cattaneo, and Titiunik (2014) propose a bias correction to the local linear regression estimator that allows for valid inference, and Calonico, Cattaneo, and Farrell (2018) provide further evidence that such bias corrections may be preferable to undersmoothing. Meanwhile, Armstrong and Kolesár (2016) show that when $\mu_{w}(x)$ is twice differentiable and $X$ has a continuous density around $c$, we can use local linear regression with a bandwidth chosen to optimize mean-squared error as the basis for bias-adjusted confidence intervals, provided we inflate confidence intervals by an appropriate, universal constant (e.g., to build 95% confidence intervals, one should use a critical threshold of 2.18 instead of 1.96). Gao (2017) characterizes the asymptotically optimal kernel for the regression discontinuity parameter under the bounded second derivative assumption with a continuous running variable. As discussed above, the value of our approach relative to this literature is that we allow for considerably more generality in the specification of the regression discontinuity design: the running variable $X$ may be discrete and/or multivariate, and the treatment boundary may be irregularly shaped.

Optimal inference with multiple running variables is less developed than in the univariate case. Papay et al. (2011) and Reardon and Robinson (2012) study local linear regression with a “small” bandwidth, but do not account for finite sample bias due to curvature. Zajonc (2012) extends the analysis of Imbens and Kalyanaraman (2012) to the multivariate case, and studies optimal bandwidth selection for continuous running variables given second derivative bounds; the inference, however, again requires undersmoothing. Keele, Titiunik, and Zubizarreta (2015) consider an approach to geographic regression discontinuity designs based on matching. To our knowledge, the approach we present below is the first to allow for uniform, bias-adjusted inference in the multivariate regression discontinuity setting.

Finally, although local methods for inference in the regression discontinuity design have desirable theoretical properties, many practitioners also seek to estimate $\tau(c)$ by fitting \mathbb{E}\left[Y_{i}\,\big{|}\,X_{i}=x\right] using a global polynomial expansion along with a jump at $c$; see Lee and Lemieux (2010) for a review and examples. However, as argued by Gelman and Imbens (2017), this approach is not recommended, as the model with the best in-sample fit may provide poor estimates of the discontinuity parameter. For example, high-order polynomials may give large influence to samples $i$ for which $X_{i}$ is far from the decision boundary $c$, and thus lead to unreliable performance.

Another approach to regression discontinuity designs (including in the discrete case) builds on randomization inference; see Cattaneo, Frandsen, and Titiunik (2015), Cattaneo, Idrobo, and Titiunik (2017) and Li, Mattei, and Mealli (2015) for a discussion. The problem of specification testing for regression discontinuity designs is considered by Cattaneo, Jansson, and Ma (2016), Frandsen (2017) and McCrary (2008).

2 Formal Considerations

2.1 Uniform Asymptotic Inference

Our main result verifies that optimized designs can be used for valid asymptotic inference about $\tau(c)$. We here consider the problem of estimating a conditional average treatment effect at a point $c$ as in (4) or (5); similar arguments extend directly to the averaging case as in (6). Following, e.g., Robins and van der Vaart (2006), we seek confidence intervals $\mathcal{I}_{\alpha}$ that attain uniform coverage over the whole regularity set under consideration:

[TABLE]

As in Armstrong and Kolesár (2018), our approach to building such confidence intervals relies on an explicit characterization of the bias of $\hat{\tau}$ rather than on undersmoothing. Our key result is as follows.

Theorem 1**.**

Suppose that we have a moment bound \mathbb{E}[(Y_{i}-\mathbb{E}\left[Y_{i}\,\big{|}\,X_{i}\right])^{q}\,\big{|}\,X_{i}=x]\leq C uniformly over all $x\in\mathbb{R}^{k}$, for some exponent $q>2$ and constant $C\geq 0$. Suppose, moreover, that $0<\sigma_{\min}\leq\sigma_{i}$ for all $i=1,\,...,\,n$ for a deterministic value $\sigma_{\min}$, and that none of the weights $\hat{\gamma}_{i}$ derived in (4) or (5) dominates all the others, i.e.,999The bound on the relative contribution of any single $\hat{\gamma}_{i}$ is needed to obtain a Gaussian limit distribution for $\hat{\tau}$. In related literature, Armstrong and Kolesár (2018) follow Donoho (1994), and side-step this issue by assuming Gaussian errors $Y_{i}(w)-\mu_{w}(X_{i})$, in which case no central limit theorem is needed. Conversely, Athey, Imbens, and Wager (2018) adopt an approach more similar to ours, and explicitly bound $\hat{\gamma}_{i}$ from above during the optimization.

[TABLE]

Then, our estimator $\hat{\tau}$ from (4) or (5) is asymptotically Gaussian,

[TABLE]

where $b\left(\hat{\gamma}\right)$ denotes the conditional bias, and $s^{2}\left(\hat{\gamma}\right)\rightarrow_{p}0$.

Now, in solving the optimization problem (4), we also obtain an explicit bound $\hat{t}$ on the conditional bias, $b(\hat{\gamma})\leq\hat{t}$, and so can use the following natural construction to obtain confidence intervals for $\tau(c)$ (Imbens and Manski, 2004),

[TABLE]

where $\alpha$ is the significance level. These confidence intervals are asymptotically uniformly valid in the sense of (9), i.e., for any $\alpha^{\prime}<\alpha$, there is a threshold $n_{\alpha^{\prime}}$ for which, if $n\geq n_{\alpha^{\prime}}$, the confidence intervals (12) achieve $\alpha^{\prime}$-level coverage for any functions $\mu_{w}(\cdot)$ in our regularity class.

Finally, whenever $X_{i}$ does not have support arbitrarily close to $c$, e.g., in the case where $X_{i}$ has a discrete distribution, the parameter $\tau(c)$ is not point identified. Rather, even with infinite data, the strongest statement we could make is that

[TABLE]

where $\text{supp}\left\{X_{i},\,W_{i}\right\}$ denotes the support of $(X_{i},\,W_{i})$. In this case, our confidence intervals (12) remain valid for $\tau(c)$; however, they may not cover the whole optimal identification interval $\mathcal{I}^{*}$. In partially identified settings, this type of confidence intervals (i.e., ones that cover the parameter of interest but not necessarily the whole identification interval) are advocated by Imbens and Manski (2004). From the perspective of the practitioner, an advantage of our approach is that intervals for $\tau(c)$ have the same interpretation whether or not $\tau(c)$ is point identified, i.e., uniformly in large samples, $\tau(c)$ will be covered with probability $1-\alpha$. Then, asymptotically, intervals (12) will converge to a point if and only if $\tau(c)$ is point identified. For a further discussion of regression discontinuity inference with discrete running variables, see Kolesár and Rothe (2018).

2.2 Implementation via Convex Optimization

In our presentation so far, we have discussed several non-parametric convex optimization problems, and argued that solving them was feasible given advances in the numerical optimization literature over the past few decades (e.g., Boyd and Vandenberghe, 2004). Here, we present a concrete solution strategy via quadratic programming over a discrete grid, and show that the resulting discrete solutions converge uniformly to the continuous solution as the grid size becomes small.101010In the case where $X$ is univariate, the resulting optimization problem is a classical one, and arguments made by Karlin (1973) imply that the weights $\hat{\gamma}_{i}$ can be written as $\hat{\gamma}_{i}=g(X_{i})$ where $g$ is a perfect spline; and our proposed optimization strategy reflects this fact. However, in the multivariate case, we are not aware of a similar simple characterization.

To do so, we start by writing the optimization problems underlying (4), (5) and (6) in a unified form. Given a specified focal point $c$, we seek to solve

[TABLE]

where $w(x)$ denotes the treatment assignment function, and $\psi$ lets us toggle between different problem types. If we want to estimate the CATE at $c$ as in (5) we set $\psi=1$, whereas if we want an optimally weighted CATE estimator as in (6) we set $\psi=0$.

To further characterize the solution to this problem, we can use Slater’s constraint qualification (e.g., Ponstein, 2004, Theorem 3.11.2) to verify that strong duality holds, and that the optimum of (14) matches the optimum of the following problem:

[TABLE]

where $k$ is the number of running variables. Here, we also implicitly used von Neumann’s minimax theorem to move the maximization over $f$ outside the $\inf_{\gamma,\,t}\{\}$ statement.

The advantage of this dual representation is that, by examining first order conditions in the $\inf_{\gamma,\,t}\{\}$ term, we can analytically solve for $\gamma$ and $t$ in the dual objective, e.g.,

[TABLE]

where $\hat{f}_{0}(\cdot)$, $\hat{f}_{1}(\cdot)$, $\hat{\lambda}_{1}$, etc., are the maximizers of (15). Carrying out the substitution results in a more tractable optimization problem over the space of twice differentiable functions $f$, along with a finite number of Lagrange parameters $\lambda_{j}$:

[TABLE]

where $\tilde{f}_{w}(x)$ in the above problem corresponds to $\lambda_{1}f_{w}(x)$ in (15), and we can recover our weights of interest via $\hat{\gamma}_{i}=-\sigma_{i}^{-2}\widehat{G}_{i}/2$ and $\hat{t}=\hat{\lambda}_{1}/(2B^{2})$. The upshot of these manipulations is that (17) can be approximated via a finite-dimensional quadratic program. In our software implementation optrdd, we use this type of a finite dimensional approximation to obtain $\hat{\gamma}_{i}$ following the construction described in the proof of Proposition 2.

Proposition 2**.**

Suppose that $X_{i}\in\mathcal{X}$ belong to some compact, convex set $\mathcal{X}\subset\mathbb{R}^{k}$. Then, for any tolerance level $\eta>0$, there exists a finite-dimensional quadratic program that can recover the solution $\hat{\gamma}$ to (17) with $L_{2}$-error at most $\eta$.

2.3 Minimizing Confidence Interval Length

As formulated in (4), our estimator seeks to minimize the worst-case mean-squared error over the specified bounded-second-derivative class. However, in some applications, we may be more interested in making the confidence intervals (12) as short as possible; and our approach can easily be adapted for this objective. To do so, consider the minimization objective in (14). Writing $\hat{v}^{2}=\sum_{i=1}^{n}\hat{\gamma}_{i}^{2}\sigma_{i}^{2}$, we see that both the worst-case mean-squared error, $\hat{v}^{2}+B^{2}\hat{t}^{2}$, and the confidence interval length in (12) are monotone increasing functions of $\hat{v}$ and $\hat{t}$; the only difference is in how they weight these two quantities at the optimum.

Now, to derive the full Pareto frontier of pairs $(\hat{v},\,\hat{t})$, we can simply re-run (14) with the term $B^{2}t^{2}$ in the objective replaced with $\lambda B^{2}t^{2}$, for some $\lambda>0$. A practitioner wanting to minimize the length of confidence intervals could consider computing this whole solution path to (14), and then using the value of $\lambda$ that yields the best intervals; this construction provides minimax linear fixed-length confidence intervals (Donoho, 1994). Since this procedure never looks at the responses $Y_{i}$, the inferential guarantees for the resulting confidence intervals remain valid.

In our applications, however, we did not find a meaningful gain from optimizing over $\lambda$ instead of just minimizing worst-case mean-squared error as in (14), and so did not pursue this line of investigation further. This observation is in line with the analytic results of Armstrong and Kolesár (2016) who showed that, when $X$ has a continuous density and $\mu_{w}(x)$ is twice differentiable, using the mean-squared error optimal bandwidth for local linear regression is over 99% efficient relative to using a bandwidth that minimizes the length of bias-adjusted confidence intervals.

Finally, although it is beyond the scope of the present paper, it is interesting to ask whether we can generalize our approach to obtain asymptotically quasi-minimax estimators for $\tau(c)$ when the per-observation noise-scale $\sigma_{i}$ needs to be estimated from the data. The resulting question is closely related to the classical issue of when feasible generalized least squares can emulate generalized least squares; see Romano and Wolf (2017) for a recent discussion.

3 Univariate Optimized Designs in Practice

To use this result in practice, we of course need to estimate the sum $\sum\hat{\gamma}_{i}^{2}\sigma_{i}^{2}$ and choose a bound $B$ on curvature. Estimating the former is relatively routine; and we recommend the following. First, we estimate $\mu_{w}(x)$ globally, or over a large plausible relevant interval around the threshold, and average the square of the residuals $R_{i}=Y_{i}-\hat{\mu}_{W_{i}}(X_{i})$ to obtain an estimate $\hat{\sigma}^{2}$ of the average value of $\sigma_{i}^{2}$. Then, we optimize weights $\hat{\gamma}_{i}$ using (4), with $\sigma_{i}^{2}\leftarrow\hat{\sigma}^{2}$. Finally, once we have chosen the weights $\gamma_{i}$, we estimate the sampling error of $\hat{\tau}$ as below, noting that the estimator will be consistent under standard conditions

[TABLE]

Conceptually, this strategy is comparable to first running local linear regression without heteroskedasticity adjustments to get a point estimate, but then ensuring that the uncertainty quantification is heteroskedasticity-robust (White, 1980). We summarize the resulting method as Procedure 1. We always encourage plotting the weights $\hat{\gamma}_{i}$ against $X_{i}$ when applying our method.

Conversely, obtaining good bounds on the curvature $B$ is more difficult, and requires problem specific insight. In particular, adapting to the true curvature $\mu_{w}(x)$ without a-priori bounds for $B$ is not always possible; see Armstrong and Kolesár (2018), Bertanha and Moreira (2016), and references therein. In applications, we recommend considering a range of plausible values of $B$ that could be obtained, e.g., from subject-matter expertise or from considering the mean-response function globally. For example, we could estimate $\mu_{w}(x)$ using a quadratic function globally, or over a large plausible relevant interval around the threshold, and then multiply maximal curvature of the fitted model by a constant, e.g., 2 or 4. The larger the value of $B$ we use the more conservative the resulting inference. In practice, it is often helpful to conduct a sensitivity analysis for the robustness of confidence intervals to changing $B$; see Figure 6 for an example.111111An interesting wrinkle is that if we are able to bound $B$ in large samples—but not uniformly—then confidence intervals built using estimated values of $B$ will have asymptotic but not uniform coverage.

1111footnotetext: Here, the algorithm assumes that all observations are of roughly the same quality (i.e., we do not know that $\sigma_{i}^{2}$ is lower for some observations than others). If we have a-priori information about the relative magnitudes of the conditional variances of different observations, e.g., some pairs outcomes $Y_{i}$ are actually aggregated over many observations, then we should run steps 2 and 3 below using appropriate inverse-variance weights. Our software allows for such weighting.

3.1 Application: The Effect of Compulsory Schooling

In our first application, we consider a dataset from Oreopoulos (2006), who studied the effect of raising the minimum school-leaving age on earnings as an adult. The effect is identified by the UK changing its minimum school-leaving age from 14 to 15 in 1947, and the response is log-earnings among those with non-zero earnings (in 1998 pounds). This dataset exhibits notable discreteness in its running variable, and was used by Kolesár and Rothe (2018) to illustrate the value of their bias-adjusted confidence intervals for discrete regression discontinuity designs. For our analysis, we pre-process our data exactly as in Kolesár and Rothe (2018); we refer the reader to their paper and to Oreopoulos (2006) for a more in-depth discussion of the data.

As in Kolesár and Rothe (2018), we seek to identify the effect of the change in minimum school-leaving age on average earnings via a local analysis around the regression discontinuity; our running variable is the year in which a person turned 14, with a treatment threshold at 1947. Kolesár and Rothe (2018) consider analysis using local linear regression with a rectangular kernel and a bandwidth chosen such as to make to make their honest confidence intervals as short as possible (recall that we can measure confidence interval length without knowing the point estimate, and so tuning the interval length does not invalidate inference). Here, we also consider local linear regression with a triangular kernel, as well as our optimized design.141414 Oreopoulos (2006) analyze the dataset using a global polynomial specification with clustered random variables, following Lee and Card (2008). However, as discussed in detail by Kolesár and Rothe (2018), this approach does not yield valid confidence intervals.

In order to obtain confidence intervals, it remains to choose a bound $B$. Following the discussion in Section 3, a 2nd-order polynomial fit with a “large” bandwidth of either 12 or 18 has a curvature of $0.006$ (the estimate is insensitive to the choice of large bandwidth); thus, we tried $B=0.006$ and $B=0.012$. We also consider the more extreme choices $B=0.003$ and $B=0.03$. For $\sigma_{i}^{2}$, we proceed as discussed in Section 3. Figure 4 shows the effective $\hat{\gamma}(X_{i})$ weighting functions for all 3 considered methods, with $B=0.012$.

We present results in Table 1. Overall, these results are in line with those presented in Figure 2. The optimized method yields materially shorter confidence intervals than local linear regression with a rectangular kernel: for example, with $B=0.03$, the rectangular kernel intervals are 11% longer. In comparison, the triangular kernel comes closer to matching the performance of our method, although the optimized method still has shorter confidence intervals. Moreover, when considering comparisons with the triangular kernel, we note that the rectangular kernel is far more prevalent in practice, and that the motivation for using the triangular kernel often builds on the optimality results of Cheng et al. (1997). And, once one has set out on a quest for optimal weighting functions, there appears to be little reason to not just use the actually optimal weighting function (4).

Finally, we note that a bound $B$ on the second derivative also implies that the quadratic approximation (8) holds with the same bound $B$. Thus, we could in principle also use the method of Armstrong and Kolesár (2018) to obtain uniform asymptotic confidence intervals here. However, the constraint (8) is weaker than the actual assumption we were willing to make (i.e., that the functions $\mu_{w}(\cdot)$ have a bounded second derivative), and so the resulting confidence intervals are substantially larger. Using their approach on this dataset gives confidence intervals of $0.0518\pm 0.0969$ with $B=0.006$ and $0.0682\pm 0.1760$ with $B=0.03$; these intervals are not only noticeably longer than our intervals, but are also longer than the best uniform confidence intervals we can get using local linear regression with a rectangular kernel as in Kolesár and Rothe (2018). Thus, the use of numerical convex optimization tools that let us solve (4) instead of (8) can be of considerable value in practice.

4 Applications with Multivariate Discontinuities

4.1 A Discontinuity Design with Two Cutoffs

We next consider the behavior of our method in a specific variant of the multiple running variable problem motivated by a common inference strategy in education. Some school districts mandate students to attend summer school if they fail a year-end standardized test in either math or reading (Jacob and Lefgren, 2004; Matsudaira, 2008), and it is of course important to understand the value of such summer schools. The fact that students are mandated to summer school based on a sharp test score cutoff suggests a natural identification strategy via regression discontinuities; however, standard univariate techniques cannot directly be applied as the regression discontinuity now no longer occurs along a point, but rather along a surface in the bivariate space encoding both a student’s math and reading scores.

We illustrate our approach using the dataset of Matsudaira (2008). As discussed above, the goal is to study the impact of summer school on future test scores, and the effect of summer school is identified by a regression discontinuity: At the end of the school year, students need to take year-end tests in math and reading; then, students failing either of these tests are mandated to attend summer school. Here, we focus on the 2001 class of graduating 5th graders, and filter the sample to only include the $n=30,741$ students whose 5th-grade math and reading scores both fall between 40 points of the passing threshold; this represents $44.7\%$ of the full sample. Matsudaira (2008) analyzed this dataset with univariate methods, by using reading score as a running variable and only consider the subset of students who passed the math exam, etc. This allows for a simple analysis, but may also result in a loss of precision.151515In order to make a formal power comparison, we need to compare two estimators that target the same estimand. In the simplest case where $\tau(x)=\tau$ is constant, our estimator (6) presents an unambiguous gain in power over those considered in Matsudaira (2008).

We present some summary statistics in Table 2. Clearly, not all students mandated to attend summer school in fact attend, and some students who pass both tests still need to attend for reasons discussed in Matsudaira (2008). That being said, the effect of passing tests on summer school attendance is quite strong and, furthermore, the treatment effect of being mandated to summer school is interesting in its own right, so here we perform an “intent to treat” analysis without considering non-compliance.

We consider both of our optimized estimators, (5) and (6), and compare weight functions $\hat{\gamma}$ learned by both methods in Figure 5. The estimator $\hat{\tau}_{c}$ is in fact quite conservative, and only gives large weights to students who scored close to $c$. Our choice of estimating the conditional average treatment effect at $(0,\,0)$ may have been particularly challenging, as it is in a corner of control-space and so does not have particularly many control neighbors.

In contrast, the weighted method $\hat{\tau}_{*}$ appears to have effectively learned matching: It constructs pairs of observations all along the treatment discontinuity, thus allowing it to use more data while canceling out curvature effects due to $\mu_{0}(x)$. As seen in Table 2, in this sample, it is much more common to fail math and pass reading than vice-versa; thus, the mean of the samples used for “matching” lies closer to the math pass/fail-boundary than the reading one.

In order to proceed with our approach, we again need to choose a value for $B$. Running a 2nd order polynomial regression on the next year’s math and reading scores for both treated and control students separately, we find the largest curvature effect among the reading score of control students; roughly a curvature of $0.46\times 40^{-2}$ along the $(1,2)$ direction. Thus, we run our algorithm with both an optimistic choice of $B=0.5\times 40^{-2}$ and a more conservative choice $B=1.0\times 40^{-2}$ (we report curvatures on the “scale” of the plots in Figure 5, such that a curvature of $1.0\times 40^{-2}$ results in a worst-case bias of $1$ in the corners of the plot).

Results are given in Table 3. As expected, the confidence intervals using the weighted method (6) are much shorter than those obtained using (5), allowing for a 0.95-level significant detection in the first case but not in the second. Since the weighting method allows for shorter confidence intervals, and in practice seems to yield a matching-like estimator, we expect it to be more often applicable than the unweighted estimator (5).

Figure 6 presents some further diagnostics for our result. The first two panels depict a sensitivity analysis for our weighted CATE result: We vary the maximum bound $B$ on the second derivative, and examine how our confidence intervals change.161616We note that, if the CATE function $\tau(\cdot)$ is not constant, then our weighted CATE estimand $\bar{\tau}_{*}=\sum_{i}W_{i}\hat{\gamma}_{i}\tau(X_{i})$ may vary with $B$. The result in Figure 6 should thus formally be interpreted as either a sensitivity analysis for the constant treatment effect parameter $\tau$ if we are willing to assume constant treatment effects, or as a robustness check for our rejection of the global null hypothesis $\tau(x)=0$ for all $x$. The result on the effect of summer school on math scores appears to be fairly robust, as we still get significant bias-aware 95% confidence intervals at $B=2\times 40^{-2}$, which is 4 times the largest apparent curvature observed in the data. The third panel plots a measure of the effective size of the treated and control samples used by the algorithm, $\text{ESS}_{w}=1/\sum_{\left\{i:W_{i}=w\right\}}\hat{\gamma}_{i}^{2}$. Although our analysis sample has almost exactly the same number of treated and control units ($\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{W}=0.501$), it appears that our algorithm is able to make use of more control than treated samples, perhaps because the treated units are “wrapped around” the controls.

Finally, it is of course natural to ask whether the bivariate specification considered here gave us anything in addition to the simpler approach used by Matsudaira (2008), i.e., of estimating the treatment effect of summer school on the next year’s math exam by running a univariate regression discontinuity analysis on only those students who passed the reading exam, and vice-versa to the effect on the reading exam.171717The estimator of Matsudaira (2008) is not exactly comparable to the two we consider. For example, when only focusing on students who passed the reading exam, his estimator effectively averages treatment effects over the math pass/fail boundary but not the reading pass/fail boundary. In contrast, we either estimate the conditional average treatment effect at a point $c$ (5), or allow for averaging over the full boundary (6). It is unclear whether the restriction of Matsudaira (2008) to averaging over only one of the two boundary segments targets a meaningfully more interpretable estimand than (6). We ran both of these analyses using our method (14), again considering bounds $B=0.5\times 40^{-2},\ 1\times 40^{-2}$ on the second derivative. For math, we obtained 95% confidence intervals of $0.083\pm 0.040$ and $0.079\pm 0.047$ for the smaller and larger $B$-bounds respectively; for reading, we obtained $0.037\pm 0.075$ and $0.030\pm 0.090$. In both cases, the corresponding bounds for the weighted estimator (6) in Table 3 are shorter, despite accounting for the possibility of bivariate curvature effects. The difference is particularly strong for the reading outcome, since our estimator $\hat{\tau}_{*}$ can also use students near the math pass/fail-boundary for improved precision.181818The corresponding headline numbers from Matsudaira (2008) are a 95% confidence interval of $0.093\pm 0.029$ for the effect on the math score, and $0.046\pm 0.045$ for the reading score; see Tables 2 and 3, reduced form estimates for 5th graders. These confidence intervals, however, do not formally account for bias. They estimate the discontinuity parameter using a global cubic fit; such methods, however, do not reliably eliminate bias (Gelman and Imbens, 2017).

4.2 A Geographic Discontinuity Design

Our final application expands on the example of Keele and Titiunik (2014), briefly discussed in Section 1.2, the goal of which is to estimate the effect of political advertising on participation in presidential elections by comparing households straddling a media market boundary. As discussed in Keele and Titiunik (2014), inference hinges on the fact that they found a media market boundary that appears not to coincide with any other major administrative boundaries, thus making it more plausible that discontinuous responses across the geographic boundary are in fact caused by varying media exposure.

The way Keele and Titiunik (2014) approach this problem is that, given a focal point $c$, they first measure Euclidean distance of each sample to it, $D_{i}=\left\lVert X_{i}-c\right\rVert_{2}$, and then use $D_{i}$ as a running variable in a univariate regression discontinuity analysis. They establish consistency of their approach via a local regression argument following Hahn et al. (2001), and also discuss advantages of their approach relative to one that considers shortest distance to the treatment boundary (rather than $c$) as input to a univariate analysis. Beyond consistency arguments, however, the resulting estimator has no guarantees in terms of statistical optimality. For example, the approach of Keele and Titiunik (2014) collapses all treated (and respectively control) observations at distance $D_{i}$ from $c$ to the same point in the univariate representation, which may increase bias.

We apply our method to this problem and, following Keele and Titiunik (2014), seek to measure the effect of media exposure on our main outcome of interest (participation in the 2008 presidential elections), as well as placebo outcomes (age, black race, political party affiliation, and gender). There were $n=24,460$ observations in this school district. Relative to the other examples considered in our paper, the responses here have a much richer dependence on the geographic running variables, thus making the problem of bounding $B$ much more delicate. To illustrate this issue, Figure 7 depicts non-parametric estimates of expected turnout and age as a function of spatial location, and reveals strong local structure. For example, we see two regions in the south-eastern and northern edges of the school district with a notably higher average age than elsewhere.

We did not find simple quadratic methods for bounding $B$ to be sufficiently sensitive to this strong local structure, and thus chose instead to set $B$ as the worst-case local curvature of a non-parametric global fit. Specifically, we ran a cross-validated ridge regression with interacted 7th-order natural splines as features on each side of the media market boundary, and then set $B$ at the 95-th percentile of the estimated operator norm curvature over all the training points. Cross-validation was carried out via the 1.se rule in the R-package glmnet (Friedman et al., 2010). We caution that this type of local choice of $B$ is only possible on fairly large datasets with strong local effects. With smaller datasets or weaker signals, cross-validation may over-smooth the estimated response function, thus resulting in an over-optimistically small value for $B$.

Table 4 shows results from our analysis, both in terms of pointwise estimates following (5) and weighted estimates (6). We set the curvature parameter $B$ separately for each different outcome, following the procedure outlined above.191919Keele and Titiunik (2014) also had a placebo outcome indicating Hispanic ethnicity. However, this outcome does not appear to have much local structure, and cross-validated ridge regression learned a constant function (resulting in $B=0$). Thus, we omit this placebo check from our present analysis; another alternative would have been to fall back on our previous approach and to select $B$ via a global quadratic fit. The $\gamma$-weighting functions underlying our analysis for the primary outcome are shown in Figure 3. Overall, we replicate the finding of Keele and Titiunik (2014), in that the media market boundary does not appear to have a substantial association with either the primary outcome of the placebo outcomes. Relative to existing methods, the value of our approach is that it allows for transparent inference of causal parameters using the original untransformed data (rather than relying on a hand-constructed univariate running variable), and uniformly accounts for spatial curvature effects.

5 Discussion

In this paper, we introduced an optimization-based approach to statistical inference in regression discontinuity designs. By using numerical convex optimization tools, we explicitly derive the minimax linear estimator for the regression discontinuity parameter under bounds on the second derivative of the conditional response surface. Because any method based on local linear regression is also a linear estimator of this type, our approach dominates local linear regression in terms of minimax mean-squared error. We also show how our approach can be used to build uniformly valid confidence intervals.

A key advantage of our procedure is that, given bounds on the second derivative, estimation of the regression discontinuity parameter is fully automatic. The proposed algorithm is the same whether the running variable is continuous or discrete, and whether or not $\tau(c)$ is point identified. Moreover, it does not depend on the shape of treatment assignment boundary when $X$ is multivariate. We end our discussion with some potential extensions of our approach.

Fuzzy regression discontinuities

In the present paper, we only considered sharp regression discontinuities, where the treatment assignment $W_{i}$ is a deterministic function of $X_{i}$. However, there is also considerable interest in fuzzy discontinuities, where $W_{i}$ is random but \mathbb{P}\left[W_{i}=1\,\big{|}\,X_{i}=x\right] has a jump at the threshold $c$; see Imbens and Lemieux (2008) for a review. In this case, it is common to interpret the indicator $\mathbf{1}\left(\left\{X_{i}\geq c\right\}\right)$ as an instrument, and then to estimate a local average treatment effect via two-stage local linear regression (Imbens and Angrist, 1994). By analogy, we can estimate treatment effects with fuzzy regression discontinuity via two-stage optimized designs as

[TABLE]

where the $\hat{\gamma}_{i}$ are obtained as in (14) with an appropriate choice penalty on the maximal squared imbalance $t^{2}$. This approach would clearly be consistent based on results established in this paper; however, deriving the best way to trade off bias and variance in specifying $\hat{\gamma}_{i}$ and extending the approach of Donoho (1994) for uniform asymptotic inference is left for future work.

Balancing auxiliary covariates

In many applications, we have access to auxiliary covariates $Z_{i}\in\mathbb{R}^{p}$ that are predictive of $Y_{i}$ but unrelated to the treatment assignment near the boundary $c$. As discussed in, e.g., Imbens and Lemieux (2008), such covariates are not necessary for identification; but controlling for them can increase robustness to hidden biases. One natural way to use such auxiliary covariates in our optimized designs is to require the weights $\hat{\gamma}_{i}$ to balance these covariates, i.e., to add a constraint

[TABLE]

to the optimization problem (4).202020The constraint (20) is a linear constraint, and so the optimization problem (14) remains a quadratic program with this constraint. In principle, if the distribution of $Z_{i}$ is in fact independent of $X_{i}$ when $X_{i}$ is near the threshold $c$, we would expect the balance conditions (20) to hold approximately even if we do not enforce them; however, explicitly enforcing such balance may improve robustness.212121A related idea would be to use the covariates $Z_{i}$ for post-hoc specification testing as in Heckman and Hotz (1989) or Imbens and Lemieux (2008). Their strategy is to obtain weights $\hat{\gamma}_{i}$ without looking at the $Z_{i}$, and then to reject the modeling strategy if (20) does not hold approximately. If we have an additive, linear dependence of $Y_{i}$ on $Z_{i}$, then enforcing balance as in (20) would also result in variance reduction, as the conditional variance of our estimator $\hat{\tau}$ would now depend on \operatorname{Var}\left[Y_{i}\,\big{|}\,X_{i},\,Z_{i}\right], which is always smaller or equal to \operatorname{Var}\left[Y_{i}\,\big{|}\,X_{i}\right].

Working with generic regularity assumptions

Following standard practice in the regression discontinuity literature, we focused on minimax linear inference under bounds on the second derivative of $\mu_{w}(\cdot)$ (e.g., Kolesár and Rothe, 2018; Imbens and Kalyanaraman, 2012). However, our conceptual framework can also be applied with higher order smoothness assumptions via bounds on the $k$-th derivative of $\mu_{w}(\cdot)$, and can easily be combined with other forms of structural information about the conditional response functions (e.g., perhaps we know from theory that the functions $\mu_{w}(\cdot)$ must be concave). Thanks to the flexibility of our optimization-based approach, acting on either of these ideas would simply involve implementing the required software using standard convex optimization libraries.

6 Appendix: Proofs

6.1 Proof of Theorem 1

By construction, we already know that

[TABLE]

Thus, to establish our desired result, it suffices to establish asymptotic Gaussianity of

[TABLE]

To do so, we use the Lyapunov central limit theorem. Thanks to our bound on the $q$-th moment of the $\varepsilon_{i}$, we know that

[TABLE]

by assumption (10). Thus, in particular, there exists a sequence $a_{n}$ such that $a_{n}\rightarrow 0$ and

[TABLE]

Now, define weights $\tilde{\gamma}_{i}$ such that $\tilde{\gamma}=\hat{\gamma}$ whenever \sum_{i=1}^{n}\mathbb{E}[(\hat{\gamma}_{i}\varepsilon_{i})^{q}\,\big{|}\,\hat{\gamma}_{i}]\,\Big{/}\,s^{q}(\hat{\gamma})\leq a_{n}, and $\tilde{\gamma}_{i}=1$ for all $i=1,\,...,\,n$ else. By Lyapunov’s central limit theorem, $\sum_{i=1}^{n}\tilde{\gamma}_{i}\varepsilon_{i}\,/\,(\sum_{i=1}^{n}\gamma_{i}^{2}\sigma_{i}^{2})^{1/2}$ is asymptotically standard normal. Moreover, we know that $\tilde{\gamma}=\hat{\gamma}$ with probability tending to 1, and so our estimator $\hat{\tau}$ must also be conditionally standard normal as claimed.

6.2 Proof of Proposition 2

Our proposed approximation to (17) has two main components. First, we need to handle the second derivative constraints $\lVert\nabla^{2}\tilde{f}_{w}(x)\rVert\leq\lambda_{1}$ for all values of $x\in\mathcal{X}$. One approach would be to write this as a positive semidefinite matrix constraint, i.e., $\nabla^{2}\tilde{f}_{w}(x)\preceq\lambda_{1}\mathcal{I}_{k\times k}$, and then rely on semidefinite programming tools; here, however, we find that we can effectively approximate $\lVert\nabla^{2}\tilde{f}_{w}(x)\rVert\leq\lambda_{1}$ via a finite number of linear constraints. Second, we need address the effect of finite differencing on our optimization.

For simplicity, throughout the following discussion, we take the curvature bound $\lambda_{1}$ to be fixed, and only consider minimization over $\lambda_{2}$, …, $\lambda_{5}$, $G$, and the $\tilde{f}_{w}$. Given this assumption, as well as the fact that $\mathcal{X}$ is bounded, we can verify that our objective is Lipschitz, which is helpful in bounding approximation errors. Then, after establishing that our approximation scheme is valid for any fixed value of $\lambda_{1}$, it follows that the approximation is also valid when we optimize over $\lambda_{1}$ because our objective is strongly convex in $\lambda_{1}$.

We start by approximating the curvature constraint. Let $A$ be any positive semidefinite matrix, and let $\mathcal{V}=\left\{v_{1},\,...,\,v_{l}\right\}$ be a set of vectors. Then, the spectral norm of $A$ is bounded by

[TABLE]

note that $\alpha(\mathcal{V})$ measures the sine of the worst angular error we may make by approximating a vector $u$ with an element from $\mathcal{V}$. The first inequality in (21) follows from the definition of the spectral norm, while the second is immediate by expressing elements in $\mathcal{V}$ in the spectral basis of $A$.

In our finite dimensional approximation to (17), we propose replacing the spectral norm constraints $\left\lVert\cdot\right\rVert$ with constraints of the type $\left\lVert\cdot\right\rVert_{\mathcal{V}}$ described above. To motivate this strategy, write $\mathcal{L}\left(\left\lVert\cdot\right\rVert,\,\lambda_{1}\right)$ for the objective value of (17) given a fixed value of $\lambda_{1}$, and let $\mathcal{L}\left(\left\lVert\cdot\right\rVert_{\mathcal{V}},\,\lambda_{1}\right)$ denote the related quantity with the approximating constraint.

Because the value of the objective cannot get worse as we increase the size the size of the constraint set, (21) implies that for any value of $\lambda_{1}$,

[TABLE]

Moreover, any feasible solution to (17) with bound $\lambda_{1}$ on the constraint can be turned into a feasible solution with bound $c\lambda_{1}$ by simply multiplying $\tilde{f}(\cdot)$ and $\lambda_{2},\,...,\,\lambda_{5}$ by $c$. When we do so, the positive term $\sum\sigma_{i}^{-2}G_{i}^{2}$ term scales quadratically in $c$, while the negative term $\lambda_{2}-\lambda_{3}$ scales linearly in $c$ (the second term must be negative, since otherwise we could improve on the objective by setting everything to 0); in particular, this implies that

[TABLE]

Chaining these two inequalities together, we find that

[TABLE]

and so, given a small enough value of $\alpha^{2}(\mathcal{V})$, the values of the two objectives match. Finally, because the objective of (17) is strongly convex in $\tilde{f}(X_{1}),\,...,\,\tilde{f}(X_{n})$, $\lambda_{2},\,\lambda_{3},\,\lambda_{4},\,\lambda_{5}$, we find that if $\mathcal{L}\left(\left\lVert\cdot\right\rVert_{\mathcal{V}},\,\lambda_{1}\right)-\mathcal{L}\left(\left\lVert\cdot\right\rVert,\,\lambda_{1}\right)$ is goes to zero, then the solutions to the two optimization problems must also converge to each other.

The remaining question is then to exhibit constructions of approximating sets $\mathcal{V}$ that make $\alpha^{2}(\mathcal{V})$ small. With a two-dimensional running variable, there exist grid-aligned sets $\mathcal{V}$ with 4 and 8 elements respectively, depicted in Figure 8, that attain values of $\alpha^{2}(\mathcal{V})$ of 0.15 and 0.05 respectively. This construction naturally extends to arbitrary $k$, resulting in sets $\mathcal{V}$ whose size scales polynomially in the inverse of the tolerance $\alpha^{2}$. In practice, though, this other computational strategies may be needed when $k$ is larger than 2 or 3.

Our second task is in showing that we can solve the $\left\lVert\cdot\right\rVert_{\mathcal{V}}$-norm optimization problem using a finite difference approximation. Let $h$ be a bandwidth, and let $\mathcal{X}_{h}$ be the set of grid points of a lattice with edge-length $h$ that fall within $h$, and assume that: the set $\mathcal{V}$ consists only of vectors whose end points lie on the unit lattice (as in, e.g., Figure 8), the basis vectors $e_{j}$ belong to $\mathcal{V}$ for $j=1,\,...,\,k$, and that the center point $c$ is in $\mathcal{X}_{h}$. Then, writing $N_{h}(x)$ for a mapping that takes all points in $\mathcal{X}$ to one of the nearest points in $\mathcal{X}_{h}$, we consider the following finite-dimensional approximation to (17):

[TABLE]

By analogy to our earlier notation, write $\mathcal{L}_{h}\left(\left\lVert\cdot\right\rVert_{\mathcal{V}},\,\lambda_{1}\right)$ for the objective value in the above problem for a fixed value of $\lambda_{1}$.

We first bound $\mathcal{L}_{h}\left(\left\lVert\cdot\right\rVert_{\mathcal{V}},\,\lambda_{1}\right)$ in terms of $\mathcal{L}\left(\left\lVert\cdot\right\rVert_{\mathcal{V}},\,\lambda_{1}\right)$. To do so, simply note that any feasible solution to the continuous problem can be turned into a solution to the discrete problem by setting

[TABLE]

and making appropriate corrections to $\lambda_{2},\,...,\,\lambda_{5}$. The goal of this affine correction is to ensure that the discrete derivate constraints $\tilde{f}^{(h)}_{w}(c+he_{j})=0$ are satisfied. Now, because $e_{j}\in\mathcal{V}$ and $\nabla\tilde{f}_{w}(c)=0$, we know that $\tilde{f}\left(c+he_{j}\right)\leq\lambda_{1}h^{2}/2$; thus, because $\mathcal{X}$ is compact, we find that $\mathcal{L}_{h}\left(\left\lVert\cdot\right\rVert_{\mathcal{V}},\,\lambda_{1}\right)\leq\mathcal{L}\left(\left\lVert\cdot\right\rVert_{\mathcal{V}},\,\lambda_{1}\right)+\mathcal{O}\left(h\lambda_{1}\right)$.

Now, to go the other direction, we start we a feasible solution $\tilde{f}^{(h)}_{w}(x)$ to the discrete solution, and need to construct a feasible solution to the continuous problem without making the objective much worse. To do so, let $\varphi_{k}(\cdot)$ for the standard Gaussian density in $\mathbb{R}^{k}$, and define

[TABLE]

with $\beta_{h}=h^{1/4}$, where $a_{w}$ is a linear correction term used to ensure that $\nabla\tilde{g}_{w}(c)=0$, and we again adjust the parameters $\lambda_{2},\,...,\,\lambda_{5}$ appropriately. Then we set $\tilde{f}_{w}(z):=\tilde{g}_{w}(z)$ for all points $z$ for which a ball of radius $h^{1/8}$ is contained inside $\mathcal{X}$, and use Whitney’s extension theorem to extend this function to the rest of $\mathcal{X}$ without increasing the supremum of the second derivative. In order to verify that $\mathcal{L}\left(\left\lVert\cdot\right\rVert_{\mathcal{V}},\,\lambda_{1}\right)\leq\mathcal{L}_{h}\left(\left\lVert\cdot\right\rVert_{\mathcal{V}},\,\lambda_{1}\right)+o(1)$, it suffices to show that the new $\tilde{f}_{w}(z)$ almost satisfy the desired curvature constraints with almost as good an objective value as before; the rest of the argument then follows as usual.

As the a full discussion is rather technical, we here focus on simply showing how to bound the second derivative of $\tilde{g}_{w}(z)$ along vectors $v\in\mathcal{V}$ at a point $z$ in the interior of $\mathcal{X}$. We emphasize that the our construction of $\tilde{g}_{w}(z)$ and $\tilde{f}_{w}(z)$ from $\tilde{f}^{(h)}_{w}(x)$ is merely a “proof of concept” used to establish asymptotic equivalence of two optimization problems, and that this construction plays no role in the actual algorithm we use.

Given a point of interest $z$, it is helpful to partition the space $\mathcal{X}_{h}$ into chords indexed by $r\in\mathcal{R}$, such that any point $x\in\mathcal{X}_{h}$ can uniquely be written as $x=a_{r}+\left(\delta_{r}+mh\right)v$ for some $r\in\mathcal{R}$ and $m\in\mathbb{Z}$. Here, $a_{r}$ is the intersection of the chord with the normal space of $v$ through $z$, and $|\delta_{r}|\leq h/2$ is the smallest correction term enabling such a representation. Now, we see that

[TABLE]

and, decomposing this quantity into a sum over chords, we get that

[TABLE]

where $\mathcal{S}_{r}\subset\mathbb{Z}$ denotes the set of valid $m$ indices in the chord $r$.

Now, the reason we used a Gaussian kernel $\varphi$ is that it has independent orthogonal components, and so we can split the $k$-dimensional Gaussian density $\varphi_{k}$ into a $(k-1)$-dimensional density over the space normal to $v$ and a $1$-dimensional density for the contribution along $v$:

[TABLE]

We are then ready to take a “Taylor expansion” where the first summand below captures the contributions of a linear effect in each chord, and the second chord has the resulting second-order terms:

[TABLE]

where the term $\Delta(r,\,m)$ simply denotes the residual from the 1-st order approximation.

If we didn’t have any discretization, the $\sum_{m\in\mathcal{S}_{r}}$ sum inside the term due to first-order effects would be exactly 0 for each $r$ individually. Here, the discretization causes error on the order of at most $h$ at each sampling point while each chord has at most on the order of $h^{-1}$ points, meaning that the total sum is bounded on the order of

[TABLE]

given our chosen growth rate on $\beta_{h}$. We also recall that $z$ is in the interior of $\mathcal{X}$, such that a ball of radius $\beta_{h}^{1/2}$ centered at $z$ is contained inside $\mathcal{X}$; thus, edge effects from $\mathcal{X}_{h}$ being finite are negligible.

Meanwhile, by feasibility of $\tilde{f}^{(h)}_{w}(\cdot)$, we know that $|\Delta(r,\,m)|\leq\lambda_{1}\left\lVert hv\right\rVert_{2}^{2}m(m-1)/2$. From here, we can verify that

[TABLE]

in the limit where $h\rightarrow 0$. To verify the last result, note that the density of points $a_{r}$ in the $k-1$ dimensional normal space to $v$ is $\left\lVert v\right\rVert_{2}h^{1-k}$; we also recall that $\left(\mathbb{E}\left[Z^{4}\right]-\mathbb{E}\left[Z^{2}\right]\right)/2=1$ for a standard Gaussian random variable $Z$.