EconBase
← Back to paper

Local-Polynomial Estimation for Multivariate Regression Discontinuity Designs

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

57,751 characters · 9 sections · 52 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Local-Polynomial Estimation for Multivariate Regression Discontinuity Designs.

\def\spacingset#1{ {#1}} \spacingset{1}

\if11 \fi

\if01 {

center[center omitted — 108 chars of source]

} \fi

abstractWe study a multivariate regression discontinuity design in which treatment is assigned by crossing a boundary in the space of multiple running variables. We document that the existing bandwidth selector is suboptimal for a multivariate regression discontinuity design when the distance to a boundary point is used for its running variable, and introduce a multivariate local-linear estimator for multivariate regression discontinuity designs. Our estimator is asymptotically valid and can capture heterogeneous treatment effects over the boundary. We demonstrate that our estimator exhibits smaller root mean squared errors and often shorter confidence intervals in numerical simulations. We illustrate our estimator in our empirical applications of multivariate designs of a Colombian scholarship study and a U.S. House of representative voting study and demonstrate that our estimator reveals richer heterogeneous treatment effects with often shorter confidence intervals than the existing estimator.

{\it Keywords:} Causal Inference, Multiple Running Variables, Distance Running variable

\spacingset{1.9}

Introduction

The regression discontinuity (RD) design takes advantage of a particular treatment assignment mechanism that is set by the running variables. \footnote{See Imbens.Lemieux2008, Lee.Lemieux2010, DiNardo.Lee2011, and \citeauthor*{Cattaneo.Idrobo.Titiunik2019} (Cattaneo.Idrobo.Titiunik2019,Cattaneo.Idrobo.Titiunik2023) for extensive surveys of RD literature} An example of such a mechanism is a scholarship that is awarded to applicants whose scores are above a threshold. The eligibility sometimes involves an additional requirement. For example, the applicants' poverty scores must be below another threshold to be eligible. These RD designs are multivariate in their running variables because a student must exceed a policy boundary in the space of multivariate running variables to be treated.

Existing approaches often handle multivariate designs as if they are univariate designs. \footnote{There are a few studies which tackled the multivariate problem as multivariate. For example, Papay.Willett.Murnane2011 and reardonRegressionDiscontinuityDesigns2012 are early exception which consider extensions of the classical polynomial based estimation of Imbens.Lemieux2008.} The most popular approach aggregates observations over the boundary to handle multivariate RD designs. For example, Matsudaira2008 considers participation in a program based on either a failure in language or math exams. Matsudaira2008 reduces the multivariate design by aggregating the language-passing students who are at the boundary of the math exam.\footnote{wongAnalyzingRegressionDiscontinuityDesigns2013 consider a decomposition of the boundary average effects into a weighted average of the boundary specific estimate of the similar strategy.} While there is no theoretical issue with the aggregation strategy, one may wish to estimate heterogeneous treatment effects across the policy boundary. \footnote{If we segment the boundary into a few intervals, then we may estimate heterogeneous effects separately for each segment. Nevertheless, finding an appropriate set of segments can be challenging and one cannot easily take its limit of this strategy to estimate the heterogeneous effect at each boundary point.}

To estimate heterogeneous treatment effects over the boundary, another popular approach constructs a running variable as the Euclidean distance from a boundary point. For example, Keele.Titiunik2015 propose a procedure to conduct the ordinary univariate regression discontinuity estimation with the Euclidean distance from a particular boundary point. \footnote{The distance approach dates back to Black1999, for example, which computes the closest boundary point for each unit and compares units of the same closest boundary point to achieve the mean effect across the boundary. In this paper, we focus on estimating the heterogeneous effects across the boundary points.} The distance approach produces a valid estimate with a valid inference under the procedure of Calonico.Cattaneo.Titiunik2014 because of its self-normalizing property of the t-statistic. \footnote{We thank an anonymous referee for this point.} The estimator is straightforward to implement, and available as Stata and R packages, rdrobust or its wrapper rdmulti (\citealp*{Cattaneo.Titiunik.Vazquez-Bare2020}).

However, the distance strategy selects bandwidth for the incorrect rate of convergence for the underlying multivariate design: the existing estimators select the optimal bandwidth for a univariate problem, but the underlying design is multivariate. As a result, the existing bandwidth selectors are suboptimal and hence their estimations are inefficient.

In this study, we document that the existing bandwidth selectors including Calonico.Cattaneo.Titiunik2014 are suboptimal when they are applied to a multivariate design with the distance from a boundary point as a running variable. We further propose a multivariate RD estimator with a Mean-Squared Error (MSE) optimal bandwidth selector. We demonstrate preferable properties of our estimator in simulation and empirical analyses.

Our estimator demonstrates favorable performances with smaller MSEs and shorter confidence intervals in most of designs. We demonstrate our estimator in two empirical contexts to compare with rdrobust. First, we apply our estimates to the multivariate RD design data of \citet*{Londono-Velez.Rodriguez.Sanchez2020} who study the impact of a Colombian scholarship program on the college attendance rate. Second, we consider a pseudo-multivariate RD design for the Lee2008 data with continuous covariates to study heterogeneous treatment effects across different values of the covariates. In the first application, our estimates exhibit shorter or comparable confidence intervals and better stability in the choice of scaling in two runnnig variables. In the second application, our estimates reveal shorter confidence intervals or richer heterogeneity.

We contribute to the literature on the estimation for RD designs. For a scalar running variable, the local-linear estimation of Calonico.Cattaneo.Titiunik2014 is the first choice for estimating treatment effects. Its statistical package, rdrobust (\citealp*{Calonico.Cattaneo.Titiunik2014a}, \citealp*{Calonico.Cattaneo.Farrell.Titiunik2017}, \citealp*{Calonico.Cattaneo.Farrell2022}), is the dominant and reliable package for a uni-variate RD design with a large sample. However, we demonstrate that the rdrobust bandwidth selector is suboptimal for a multivariate RD design when the distance from a boundary point is its univariate running variable. We further provide an alternative local-linear estimator with an optimal bandwidth selector.

We note that multivariate estimations are also available in a non-kernel bias-aware procedure such as Imbens.Wager2019 and Kwon.Kwon2020a which are derived from Armstrong.Kolesar2018, for example. These bias-aware methods may fully adapt the underlying distribution of the running variable. The bias-aware approach is a valid but different alternative to the kernel procedure because they employ the worst-case second derivative as the tuning parameter instead of the bandwidth. Given that the two approaches are in different principle, we contribute to fill a missing piece in the kernel procedure for a multivariate RD design with an optimal bandwidth selector.

The most related study is the recent work by cattaneoEstimationInferenceBoundary2025 which has reported an important boundary bias in the distance approach when the evaluation point is at the corner of the policy boundary. Combined with our arguments about the suboptimality of the distance approach, both contributions jointly alert that the univariate distance approach should not be used for the multivariate designs not just at the corner but also at any boundary points. Contributions in our estimator is also complementary. On the one hand, we allow for selecting dimension specific bandwidths which can differ substantially when the scaling of the running variables differ as demonstrated in our simulation. On the other hand, cattaneoEstimationInferenceBoundary2025 provide a uniform inference across evaluation points. Hence, both contributions are complementary both in terms of studying the distance approach as well as developing an appropriate estimator.

The remainder of the paper is organized as follows. We document the problem of the existing approach and introduce our estimator in Section 2. In Section 3, we evaluate the proposed estimator in Monte Carlo simulations and in empirical studies by Londono-Velez.Rodriguez.Sanchez2020 and a modification of Lee2008. Finally, we conclude the paper and discuss future challenges in Section 4.

Methods

Set up and identification

Consider a multivariate RD design for a student with a pair of test scores $(R_1,R_2)$. For example, we consider a program that accepts students whose scores exceed their corresponding thresholds $(c_1,c_2)$. In this program, the eligibility is set by a treatment region $\mathcal{T} = \{(R_1,R_2) \in \mathbb{R}^2: R_1 \geq c_1, R_2 \geq c_2 \}$ (Figure (ref) (a)). For another example, consider a program that accepts students whose total score exceeds a single threshold $c_1 + c_2$. The eligibility is set by another region $\mathcal{T} = \{(R_1,R_2) \in \mathbb{R}^2: R_1 + R_2 \geq c_1 + c_2 \}$ (Figure (ref) (b)). In general, we consider a binary treatment $D \in \{0,1\}$ and associated pair of potential outcomes $\{Y(1),Y(0)\}$ such that $ Y = D Y(1) + (1 - D)Y(0) $ for an observed outcome $Y \in \mathbb{R}$. We consider a sharp RD design with a vector of running variables $R \in \mathcal{R} \subseteq \mathbb{R}^d$ for some integer $d \geq 1$. Specifically, let $\mathcal{T}$ be the treatment region, which is an open subset of the support, $\mathcal{R}$. Let $\mathcal{T}^C$ be the complement of the closure of $\mathcal{T}$. This $\mathcal{T}^C$ is the control region, and both $\mathcal{T}$ and $\mathcal{T}^C$ have non-zero Lebesgue measures, and $D = 1\{R \in \mathcal{T}\}$.

figure[figure omitted — 541 chars of source]

We consider the i.i.d. sample of $(Y,D,R)$, $(Y_i,D_i,R_i)_{i \in \{1,\ldots, n\}}$, where $R_i = (R_{i,1},R_{i,2})$ and $R = (R_1,R_2)$. Let $c$ be a particular point on the boundary of the closure of $\mathcal{T}$. Our target parameter is $\theta(c) := \lim_{r \rightarrow c, r \in \mathcal{T}} E[Y(1)-Y(0)|R=r] - \lim_{r \rightarrow c, r \in \mathcal{T}^C} E[Y(1)-Y(0)|R=r]$. In the following section, we focus on the issues in estimating the given identified parameter $\theta(c)$. Under the following assumptions (\citealp*{Hahn.Todd.Klaauw2001}; Keele.Titiunik2015), $\theta(c)$ is the average treatment effect (ATE) at each point of the boundary $c$:

proposition(Keele.Titiunik2015, Proposition 1) If $E[Y(1)|R=r]$ and $E[Y(0)|R=r]$ are continuous in $r$ at all points $c$ of the boundary of the closure of $\mathcal{T}$; $P(D_i = 1) = 1$ for all $i$ such that $R_i \in \mathcal{T}$; $P(D_i = 1) = 0$ for all $i$ such that $R_i \in \mathcal{T}^C$, then, $\theta(c) = E[Y(1) - Y(0)|R = c]$ for all $c$ in the boundary.

Issues in Conventional Estimators

To estimate heterogeneous treatment effects over the boundary points, one often employs the distance strategy which explicitly reduces a multivariate running variable to a scalar distance measure. A frequent choice is the Euclidean distance from a point or the closest boundary Keele.Titiunik2015. The distance strategy can be easily implemented in most designs via the local-linear estimation Fan.Gijbels1992 for the uni-variate RD designs with a MSE optimal bandwidth selection. However, the existing bandwidth selector is not rate optimal when it uses for the distance strategy for a multivariate design.

Our first observation is the property of the density function of the distance running variable at a boundary point. Let $Z_i$ be the scalar running variable as a distance from a boundary point $c$. Then its density $f_Z(z)$ shrinks to zero as it approaches the boundary when the distance $\tilde{d}$ bounds the Euclidean distance with some constant:

propositionLet $\tilde{d}(\cdot, \cdot)$ be a distance on $\mathbb{R}^d$ such that $\bar{c} \|a - b\| \leq \tilde{d}(a,b)$ for any $a, b \in \mathbb{R}^d$ and some constant $\bar{c} > 0$. Here $\|a - b\|$ is the Euclidean distance between $a = (a_1,\ldots, a_d)'$ and $b = (b_1,\ldots, b_d)'$. Define $Z_i = \tilde{d}(R_i,c)$ with $c = (0,\ldots, 0)'$ and assume that $R_i$ and $Z_i$ have density functions $f$ and $f_Z$, respectively. \footnote{The boundary point value $c$ is set to zero for illustration. The same argument applies in general by normalizing the running variables with respect to the boundary point. The distance $\tilde{d}$ includes the Euclidean norm $\|a - b\|$, $\ell^{\infty}$-norm $\|a - b\|_{\infty} = \max_{1 \leq j \leq d} |a_j - b_j| \geq (1/d) \|a - b\|$, and $\ell^1$-norm $\|a - b\|_1 = \sum_{j=1}^d |a_j - b_j| \geq \|a - b\|$.} Assume that $f$ and $f_Z$ are continuous. Then we have $f_Z(z) \rightarrow 0$ as $z \rightarrow 0$.
proofBy construction of $Z_i$, for $z > 0$, \begin{align*} \int_0^z f_Z(r) dr =& P(Z_i \leq z) = P(\tilde{d}(R_i,0) \leq z) \leq P(\bar{c} \|R_i\| \leq z) = P(\|R_i\| \leq z/\bar{c})\\ =& \int_0^{z/\bar{c}} t\left(\int_0^{2\pi} f(t \cos \theta, t \sin \theta) d\theta \right)dt\\ =& \int_0^z (1/\bar{c})^{2} r \left(\int_0^{2\pi} f(\bar{c}^{-1} r \cos \theta , \bar{c}^{-1} r \sin \theta) d\theta \right) dr \end{align*} where the last equality uses the change of variable $r = \bar{c} t$. If $f$ is continuous, $f_Z(0) = 0$ by using the above inequality. Since $f_Z$ is continuous, the statement follows.

To illustrate the proposition in an example, consider $R_i = (R_{1i}, R_{2i})$ where $R_{1i}$ and $R_{2i}$ independent each other, and $R_{1i} \sim U[-1,1]$ and $R_{2i} \sim U[0,1]$. The distribution function of $Z_i = \|R_i\|$ is $P(Z_i \leq z) = P(R_{1i}^2 + R_{2i}^2 \leq z^2) = (\pi/4) z^2$. The half-circle area shrinks to zero at the order of $z^2$ as $z$ approaches the value $0$ at the boundary point $(0,0)$.

This zero-density problem itself may appear to be not an immediate problem for the Calonico.Cattaneo.Titiunik2014 (henceforth, CCT or rdrobust) estimator as its bandwidth selector does not estimate the density directly. Nevertheless, it leads to another problem that makes CCT bandwidth selector suboptimal when it is used for the distance strategy.

To demonstrate its mechanism, first we consider the simpler Imbens.Kalyanaraman2012 (IK) bandwidth selector with the Euclidean distance running variable $Z_i = D_i \|R_i\| - (1-D_i) \|R_i\|$. The IK bandwidth selector for the distance strategy takes the following form \[ \hat{h}_{IK} = C \cdot \left( \frac{\hat{V}_{IK} / \hat{f}_Z(0)}{\hat{B}_{IK}} \right)^{1/5} n^{-1/5} \] where $\hat{B}_{IK}$ depends on the regularization term and the estimator of the second derivative of $E[Y_i|Z_i=z]$, $\hat{V}_{IK}$ depends on the estimator of the conditional variance $V(Y_i|Z_i=z)$, $\hat{f}_Z(0) = \frac{1}{nh_{\mathrm{pilot}}} \sum_{i=1}^n K(Z_i/h_{\mathrm{pilot}})$ with some kernel function $K$ and the pilot bandwidth $h_{\mathrm{pilot}}$. In Online Appendix (ref), we show that $\hat{f}_{Z}(0)$ converges to zero while $h_{\mathrm{pilot}}^{-1} \cdot \hat{f}_{Z}(0)$ converges to a positive constant when $f(0)$ is strictly positive. Hence, if $\hat{V}_{IK}$ and $\hat{B}_{IK}$ converge to strictly positive constants, then the variance term $\hat{V}_{IK}/\hat{f}_{Z}(0)$ diverges while $h_{pilot} \hat{V}_{IK}/\hat{f}_{Z}(0)$ converges to a strictly positive constant. Hence, we obtain $ \hat{h}_{IK} = O_p(h_{\mathrm{pilot}}^{-1/5} n^{-1/5}), $ and the rate of $\hat{h}_{IK}$ depends on the pilot bandwidth. For instance, if $h_{\mathrm{pilot}} = O_p(n^{-1/5})$, then $\hat{h}_{IK} = O_p(n^{-4/25})$, which is suboptimal in a two-dimensional estimation problem.

This diverging variance term problem arises in CCT bandwidth selector as well even though it avoids the density estimation directly. The CCT bandwidth selector has the form

equation[equation omitted — 140 chars of source]

where $\tilde{B}_{CCT}$ depends on the regularization term and the estimator of the second derivative of $E[Y_i|Z_i=z]$ and $\tilde{V}_{CCT}$ depends on the estimators of the conditional variance $V(Y_i|Z_i=z)$ and $f_Z(0)$ but they are estimated in a sandwich form so that the density estimation does not arise explicitly. Specifically, in Online Appendix F, we study the variance term $ \tilde{V}_{CCT} = n h_{\mathrm{initial}} \left\{ \hat{V}_{+}(h_{\mathrm{initial}}) + \hat{V}_{-}(h_{\mathrm{initial}}) \right\} $ for $h_{\mathrm{initial}} = O_p(n^{-1/5})$ where the variance term is the sum of the elements of sandwich forms: $\hat{V}_{+}(h) = e_1' \Gamma_{+}(h)^{-1} \Psi_{+}(h) \Gamma_{+}(h)^{-1} e_1 /n$, and $\hat{V}_{-}(h) = e_1' \Gamma_{-}(h)^{-1} \Psi_{-}(h) \Gamma_{-}(h)^{-1} e_1 /n$ for $r_1(z) = (1,z)'$ and $e_1 = (1,0)'$. \footnote{See Online Appendix F for the formal definitions for $\Gamma_+, \gamma_-, \Psi_+,$ and $\Psi_-$. We consider a simplified version for the $\Psi_+$ and $\Psi_{-}$ matrices which take a known variance function instead of the original formula with the plug-in estimates of the residual variance function.} As Proposition F.1, we show that $h^{-1}\Gamma_+(h)$ and $\Psi_+(h)$ converges to constant matrices, instead of $\Gamma_+(h)$ and $h \Psi_+(h)$ convergence as required in the original procedure. Hence, the positive side of the original variance term $n h_{\mathrm{initial}} \hat{V}_{+}(h)$ diverges while \[ n h_{\mathrm{initial}}^2 \hat{V}_{+}(h_{\mathrm{initial}}) = e_1' \left\{ h_{\mathrm{initial}}^{-1}\Gamma_{+}(h_{\mathrm{initial}}) \right\}^{-1} \Psi_{+}(h_{\mathrm{initial}}) \left\{ h_{\mathrm{initial}}^{-1}\Gamma_{+}(h_{\mathrm{initial}}) \right\}^{-1} e_1 \] converges to a positive constant. If $\tilde{B}_{CCT}$ also converges to a positive constant, we have \[ \hat{h}_{CCT} = C \cdot \left( \frac{h_{\mathrm{initial}} \tilde{V}_{CCT}}{\tilde{B}_{CCT}} \right)^{1/5} h_{\mathrm{initial}}^{-1/5} n^{-1/5} = O_p\left( h_{\mathrm{initial}}^{-1/5} n^{-1/5} \right). \] Hence, if $h_{\mathrm{initial}} = O_{p}(n^{-1/5})$, the convergence rate of $\hat{h}_{CCT}$ is $n^{-4/25}$ which is the same suboptimal rate as the IK bandwidth for the multivariate distance bandwidth selector. See the complete discussion for the Online Appendix F.

Our Estimator

Given the suboptimality of the distance strategy, we propose a multivariate RD estimator for the heterogeneous treatment effect over the boundary with the MSE optimal bandwidths.

We demonstrate our estimator in a special case of two-dimensional running variables. Consider the following local-linear estimator $\hat{\beta}^+(c)= (\hat{\beta}^+_0(c), \hat{\beta}^+_1(c),\hat{\beta}^+_2(c))'$ \[ \hat{\beta}^+(c) = \mathop{\rm arg~min}\limits_{(\beta_0,\beta_1,\beta_2)' \in \mathbb{R}^{3}}\sum_{i=1}^{n}(Y_i - \beta_0 - \beta_1(R_{i,1} - c_1) - \beta_{2}(R_{i,2}-c_2))^2K_h\left(R_i - c\right)1\{R_i \in \mathcal{T}\} \] where $ K_{h}(R_i - c) = K\left(({R_{i,1} - c_1) / h_1},{(R_{i,2} - c_2) / h_2}\right) $ and each $h_j$ is a sequence of positive bandwidths such that $h_j \to 0$ as $n \to \infty$. Similarly, let $\hat{\beta}^-(c)$ be the estimator using $1\{R_i \in \mathcal{T}^c\}$ subsample. Hence, our multivariate RD estimator at $c$ is $\hat{\beta}^+_0(c) - \hat{\beta}^-_0(c)$. Our estimator uses the theoretical results of Ruppert.Wand1994, Masry1996, and guMultivariateLocalPolynomial2015 for the multivariate RD designs. Specifically, we employ guMultivariateLocalPolynomial2015 with a slightly extended result such as non-product kernels and explicit higher-order expressions to allow us to conduct the Calonico.Cattaneo.Titiunik2014 style bias-correction procedure.

As we consider a random sample, the treated sample is independent of the control sample. Without the loss of generality, we consider the following nonparametric regression models for each sample: $Y_i = m_+(R_i) + \varepsilon_{+,i},\ E[\varepsilon_{+,i}|R_i] = 0,\ i \in \{1,\dots, n: R_i \in \mathcal{T}\}$ and $Y_i = m_-(R_i) + \varepsilon_{-,i},\ E[\varepsilon_{-,i}|R_i] = 0,\ i \in \{1,\dots, n: R_i \in \mathcal{T}^C\}.$

For the asymptotic normality, we impose the following regularity conditions that are standard in kernel regression estimations. We provide the conditions under its general possible form. In Online Appendix (ref), we present the general results for $p$th order local-polynomial estimation with $d$-dimensional running variables. The general results in the Online Appendix are the basis of the bias correction procedure of our estimator.

In Assumption (ref), we assume the existence of a continuous density function for the running variable $R$. Assumption (ref) is the regularity conditions for a kernel function. We select a particular set of kernel functions for our subsequent analysis. Assumption (ref) imposes a set of smoothness conditions for the conditional mean functions $m_+$ and $m_-$ and for the conditional moments of residuals $\varepsilon_{+,i}$ and $\varepsilon_{-,i}$. Assumption (ref) specifies the rate of convergence of the vector of bandwidths $\{h_1,\ldots, h_d\}$ relative to the sample size $n$.

assumptionLet $U_{r}$ be a neighborhood of $r=(r_1,\dots, r_d)' \in \mathcal{R}$. \begin{itemize} • The vector of random variables $R_i$ has a probability density function $f$. • The density function $f$ is continuous on $U_r$ and $f(r)>0$. • For each $r \in \mathcal{R}$ on the boundary of the treated region $\mathcal{T}$, \textcolor{black}{there exists $\delta>0$ such that $\mathcal{T} \cap \prod_{j=1}^d[r_j -\delta, r_j + \delta] = [r_2,r_2+\delta] \times \prod_{j \neq 2} [r_j -\delta, r_j + \delta]$.} \end{itemize}

Condition (c) means that, in a sufficiently small neighborhood of the point $r$ of interest on the boundary of $\mathcal{T}$, the boundary is linearly separated. \footnote{This condition (c) excludes the evaluation of the corner point. The following implementation and the actual numerical simulation and empirical analyses avoid the evaluation exactly at the corner of the boundary. Adapting the finding of cattaneoEstimationInferenceBoundary2025, the same statement should follow with a relaxed condition (c) which allows for the corner point. One may allow for more complex boundary structures in the neighborhood of $r$ in a different setting such as those in cattaneoEstimationInferenceBoundary2025; however, extending our framework to their setting is beyond the scope of this paper.}

assumptionLet $K :\mathbb{R}^{d} \to [0,\infty)$ be a kernel function such that \textcolor{black}{ \begin{itemize} • $\int K_\pm(z)dz = 1$ where $K_+(z)=K(z)1\{z_2\geq 0\}$, and $K_-(z) = K(z)1\{z_2 <0\}$. • The kernel function $K$ is bounded and there exists a constant $C_{K}>0$ such that $K$ is supported on $[-C_{K},C_{K}]^d$. • Define $\check{z}:= (1,(z)'_1,\dots,(z)'_p)',\ (z)_L = \left(\prod_{\ell=1}^{L}z_{j_\ell}\right)'_{1\leq j_1 \leq \dots \leq j_L\leq d},\ 1\leq L \leq p. $ The matrix $S_\pm=\int K_\pm(\bm{z})\left( \begin{array}{c} 1 \\ \check{\bm{z}} \end{array} \right) (1\ \check{\bm{z}}')d\bm{z}$ is non-singular. \end{itemize} }
assumptionLet $U_{r}$ be a neighborhood of $r \in \mathcal{R}$. \begin{itemize} • The mean function $m_+$ is $(p+1)$-times continuously partial differentiable on $U_r$ and define $\partial_{j_1\dots j_L}m_+(r):= \frac{\partial m_+(r)}{\partial r_{j_1}\dots r_{j_L}}$, $1 \leq j_1,\dots, j_L \leq d$, $0 \leq L \leq p+1$. When $L=0$, we set $\partial_{j_1 \dots j_L}m_+(r) = \partial_{j_0}m_+(r)= m_+(r)$. The parallel restriction holds for the mean function $m_-$. • The variance function $\sigma_+^2(z) = E[\varepsilon_{+,i}^2|R_i = z]$ and $\sigma_-^2(z) = E[\varepsilon_{-,i}^2|R_i = z]$ are continuous at $r$. • There exists a constant $\delta>0$ such that $\sup_{z \in U_r}E[|\varepsilon_{+,1}|^{2+\delta}|R_1=z]\leq U(r) < \infty$ and $\sup_{z \in U_r}E[|\varepsilon_{-,1}|^{2+\delta}|R_1=z]\leq U(r) < \infty$ \end{itemize}
assumptionAs $n \to \infty$, \begin{itemize} • $h_j \to 0$ for $1 \leq j \leq d$, • $nh_1 \cdots h_d \times h_{j_1}^2 \dots h_{j_p}^2 \to \infty$ for $1 \leq j_1 \leq \dots \leq j_{p} \leq d$, • $nh_1 \cdots h_d \times h_{j_1}^2 \dots h_{j_p}^2h_{j_{p+1}}^2 \to c_{j_1\dots j_{p+1}} \in [0,\infty)$ for $1 \leq j_1 \leq \dots \leq j_{p+1} \leq d$. \end{itemize}

Under these assumptions, we establish the asymptotic normality of our estimator $\hat{\beta}^+_0(c) - \hat{\beta}^-_0(c)$. \textcolor{black}{The result follows from Theorem (ref) in Appendix (ref).}

theorem[Asymptotic normality of local-linear estimators] Under Assumptions (ref), (ref), (ref) and (ref) for $r=c$ with $d = 2$ and $p = 1$, as $n \to \infty$, we have \begin{align*} \sqrt{nh_1 h_2} &\left[\right.H^{ll}\left((\hat{\beta}^+_0(c) - \hat{\beta}_0^-(c)) - e'_1(M_+(c) - M_-(c))\right)\\ -& \textcolor{black}{e'_1 \{S_+^{-1}B^{(2,1)}M_{+,n}^{(2,1)}(c) - S_-^{-1}B^{(2,1)}M_{-,n}^{(2,1)}\}}\left.\right] \\ &\textcolor{black}{\stackrel{d}{\to} N\left( \0 , e'_1\left\{\frac{\sigma_{+}^2(c)}{f(c)} S_+^{-1}\mathcal{K}_+S_+^{-1} +\frac{\sigma_{-}^2(c)}{f(c)} S_-^{-1}\mathcal{K}_-S_-^{-1}\right\}e_1 \right), } \end{align*} for $e_1 = (1,0,0)'$ and $H^{ll} = {\rm diag}(1,h_1,h_2) \in \mathbb{R}^{3 \times 3}$ where \begin{align*} M_+(c) &= \left(m_+(c),\partial_1 m_+(c), \partial_2 m_+(c)\right)', M_-(c) = \left(m_-(c),\partial_1 m_-(c), \partial_2 m_-(c)\right)',\\ M_{+,n}^{(2,1)}(r) &= \left(\frac{\partial_{11}m_+(c)}{2}h_1^{2},\partial_{12}m_+(c)h_1 h_2,\frac{\partial_{22}m_+(c)}{2}h_2^{2}\right)',\\ M_{-,n}^{(2,1)}(r) &= \left(\frac{\partial_{11}m_-(c)}{2}h_1^{2},\partial_{12}m_-(c)h_1 h_2,\frac{\partial_{22}m_-(c)}{2}h_2^{2}\right)', and \\ B^{(2,1)} &= \int \left( \begin{array}{c} 1 \\ \check{\bm{z}} \end{array} \right) (\bm{z})'_{2}d\bm{z},\ \textcolor{black}{\mathcal{K}_\pm= \int K_\pm^2(\bm{z})\left( \begin{array}{c} 1 \\ \check{\bm{z}} \end{array} \right) (1\ \check{\bm{z}}')d\bm{z}.} \end{align*}

Given the bias and variance expressions in Theorem (ref), we may find the common bandwidth $h = h_1 = h_2$ that minimizes the following asymptotic expansion of the mean-squared error (MSE) of $\hat{m}_+(c) - \hat{m}_{-}(c)$: for $e_1 = (1,0,0)'$, \textcolor{black}{

align*[align* omitted — 533 chars of source]

} In general, it would be more intuitive and reasonable to consider heterogeneous bandwidth $h_1 \neq h_2$, and our main numerical illustrations are based on the heterogeneous bandwidths. In one of our empirical analysis dataset, two running variables take quite different ranges of values because one of the running variable has twice or three times larger scale than the other. If we use a common bandwidth, a possibly awkward squared area will be used for the estimation while it may be too large for one dimension and too small for the other dimension. One can avoid such an awkward situation by rescaling the running variables appropriately, but the results may change substantially by rescaling. The heterogeneous bandwidths allows users to avoid such a difficult rescaling task and use the original scaling for the estimation. See Appendix (ref) for the details for the heterogeneous case.

Numerical Results

We demonstrate the numerical properties of our estimator in Monte Carlo simulations and empirical applications. Numerical simulations use the first empirical context of a Colombian scholarship, \citet*{Londono-Velez.Rodriguez.Sanchez2020,londono_data_2020}. Specifically, we evaluate the performances of our estimator in simulations which take higher-order approximations of the Colombian data as true data generating processes, and in empirical application with the actual dataset. In their application, the scholarship of interest is primarily determined by two thresholds: merit-based and need-based. As a result, a policy boundary exists instead of a single cutoff. Figure (ref) is a scatter plot of two running variables with $30$ boundary points. The $30$ points are selected from taking $15$ points from the maximum value across the boundary to $0$. In the simulations below, we found that a largest point in SISBEN boundary is challenging for all methods evaluated, and we evaluate $28$ points denoted as red filled circles after removing the extreme boundary points denoted as blank black circles in the empirical application later. We explain the institutional details further in Section (ref).

figure[figure omitted — 852 chars of source]

Simulation Results

Given the dataset, we constructed four designs which are all two-dimensional saturated higher-order polynomial approximations of the conditional expectation functions at four boundary points. Specifically, we use the fully saturated polynomials up to fourth orders plus the fifth order terms for $X$ and $Y$ each. The four boundary points are at a higher SISBEN (need-based) boundary $(7)$, an intermediate SIBEN boundary $(13)$, an intermediate SABER11 (merit-based) boundary $(19)$ and a higher SABER11 (merit-based) boundary $(25)$. Figures (ref) show the two-dimensional plots of the mean functions. For each draw of a simulation sample, we draw a random sample of two-dimensional running variables as $R_1 \sim U[-1,1]$ and $R_2 \sim 2 \times Beta(2,4)-1$ independent of each other over a rescaled rectangular support, and generate the outcome variable as $m(R_{i1},R_{i2}) + \epsilon_i$ where $\epsilon_i \sim N(0,0.1295^2)$.

figure[figure omitted — 1,155 chars of source]

We compare the quality of our estimator relative to the rdrobust estimation. Figure (ref) shows histograms of realized estimates of $10,000$ times replications for the primary data generating process. The light-colored histograms of our 2D local poly estimates tend to have thinner shapes than the dark-colored histograms of rdrobust estimates.

figure[figure omitted — 959 chars of source]
table[table omitted — 1,766 chars of source]

We report the detailed results in Table (ref). Our first observation is that estimation with heterogeneous bandwidths $h_1 \neq h_2$ matters. The common bw estimator is a version of our 2D local poly estimator that imposes $h_1 = h_2$. For all designs, our preferred 2D local poly - diff bw has smaller or approximately equal bias than common bw. The better bias correction with heterogeneous bandwidths selection appears to induce smaller root MSEs for most designs while our 2D local poly - diff bw estimator is stable and maintaining the coverage rates above $95\%$ in all four designs.

Greater differences appear in comparison of our preferred estimator with rdrobust. The RMSE of our estimator is smaller than that of the rdrobust for all designs. In particular, the RMSE is less than the half of the RMSE in rdrobust estimates for the Designs 1 and 2. Furthermore, the confidence intervals of our estimator are also shorter than that of the rdrobust for most designs. Hence, our estimates are more efficient than the rdrobust estimates and the efficiency conveys its greater performance in the inferences. Interestingly, the bias can be smaller in rdrobust than in our estimator while its RMSE is always greater than in our estimator and their coverages are always below $95\%$. This result of the \textit{rdrobust} estimator is consistent with our earlier methodological analyses. The \textit{rdrobust} estimator chooses its bandwidth as if it is a univariate design; hence, their bandwidth selector chooses a suboptimal bandwidth which overly reduce bias relative to variance. \footnote{We report the summary statistics of the bandwidths used in Table (ref). We also conduct a parallel simulation study with a binary response via a linear probability model of the same polynomial in Table (ref) and (ref).}

Finally, we conduct a parallel exercise across all $30$ points. \footnote{Note that the underlying sampling supports are different from the earlier simulation results for the four points. Unlike the four points, which are relatively center in the support, some of the $30$ boundary points are outside of the originally constructed rectangular supports for the four designs.} Table (ref) summarizes the performance comparisons across $30$ points. Except for an extreme behavior that appears in the max among $30$ points, our estimator performs favorably relative to rdrobust. \footnote{The low performing one is at point 1 where all three estimators are poorly performed. We realize that the largest boundary points (points 1 and 30) are too extreme. Hence, we exclude the extreme points from the boundary points to evaluate in the empirical analysis. See Online Appendix (ref) Table (ref) for the results of all $30$ points.}

table[table omitted — 2,243 chars of source]

Applications

We first illustrate our estimator through an empirical application of a Colombian scholarship, \citet*{Londono-Velez.Rodriguez.Sanchez2020, londono_data_2020}. From 2014 to 2018, the Colombian government operated a large-scale scholarship program called Ser Pilo Paga (SPP). The scholarship loan covers “the full tuition cost of attending any four-year or five-year undergraduate program in any government-certified `high-quality' university in Colombia” Londono-Velez.Rodriguez.Sanchez2020. The eligibility of the SPP program is based on two thresholds. The first threshold is merit-based, determined by the nationally standardized high school graduation exam, SABER 11. In 2014 of Londono-Velez.Rodriguez.Sanchez2020's study period, the cutoff was the top 9% of the score distribution. The second threshold is need-based, and is determined by the eligibility of the social welfare program, SISBEN. SISBEN-eligible families are roughly the poorest 50 percent. \footnote{Students must be also accepted by an eligible college in Colombia to receive the scholarship. Hence, the impact of exceeding both thresholds is not the impact of the program itself owing to noncompliance. The estimand is the impact of the program eligibility, which is the intention-to-treat effect.} The sample consists of $347,673$ observations of the control units and $15,423$ observations of the treated units.

The aggregation approach is the empirical strategy of Londono-Velez.Rodriguez.Sanchez2020. They run rdrobust separately for two boundaries: the merit-based criterion (SAVER11) and the need-based criterion (SISBEN) as in Figure (ref). They report the effect of exceeding the merit-based (SABER11) threshold on enrollment in any eligible college is $0.32$ with a standard error of $0.012$ for the need-based (SISBEN) eligible subsample, and the effect of exceeding the need-based (SISBEN) threshold on enrollment in any eligible college is $0.274$ with a standard error of $0.027$ for the merit-based (SABER11) eligible subsample. Students with the need eligibility in the $x$-axis boundary of Figure (ref) have a slightly higher effect than students with the merit eligibility in the $y$-axis boundary of Figure (ref). Indeed, their strategy captures certain heterogeneity in the sub-populations, albeit with richer heterogeneity within. The SISBEN threshold students are heterogeneous in their SABER11 scores; the SABER11 threshold students are heterogeneous in their SISBEN scores.

We estimate the heterogeneous effects over the entire boundary. We summarize our results in Figure (ref) with Panel (a) of the SABER = 0 boundary and Panel (b) of the SISBEN = 0 boundary. The dark-colored intervals are the pointwise $95\%$ confidence intervals from our estimates at each boundary point value, and the light-colored intervals are the pointwise $95\%$ confidence intervals from the rdrobust estimates. For most points, the two estimates show similar patterns across the boundary points with a notable difference in the length of the confidence intervals. Our estimates exhibit shorter confidence intervals than rdrobust when there are enough neighboring observations around the boundary points (such as SISBEN values from 2 through 24 in (a) and SABER values from 7 through 21 in (b)). On the other hand, our confidence intervals widen when there are only a few neighboring observations around the boundary points (such as SISBEN values at 44 and 48 in (a) and SABER values from 70 or more in (b)). Hence, our estimates are more stable for various designs and efficient at least when the effective sample size is large enough.

Both estimates suggest substantial heterogeneity in the effects among the merit-eligible students (Panel (b)) but not among the need-eligible students (Panel (a)). Specifically, the program has similar effects among the majority of students, but has no impact on extremely capable students. The null effect for extremely capable students is reasonable because they would have received other scholarships to attend college anyway. Consequently, the program could have benefited from accepting a larger number of students with higher household incomes because their impact is expected to be similar.

figure[figure omitted — 703 chars of source]

We further assess the stability of our estimates relative to rdrobust by changing the scalings of the two running variables. Figure (ref) compares estimates with and without scaling by the absolute maximum values of each axis. Compared with the Panel (a) and (c) which exhibit substantial changes in the estimated confidence intervals of rdrobust, our estimates in Panel (b) and (d) show the stability in the underlying (relative) scale of the running variables. An appropriate relative scaling of the two axes is hardly known. Hence, our approach is superior in handling the relative scaling of the two-dimensional data as is because our estimator is more robust against the choice of scaling.

figure[figure omitted — 1,045 chars of source]

We further apply our procedure to the dataset used in Lee2008 (also in caugheyDataCitation, caugheyDataCitation, Caughey.Sekhon2011) that studies the U.S. House Elections and finds the positive significant incumbent margin. There are a few baseline covariates with continuous variations as reported in Caughey.Sekhon2011. We use four baseline covariates: percentages of black voters, foreign born voters, government worker voters, and of urban areas for each electoral district. See their scatter plots and evaluation points in Online Appendix (ref) Figure (ref). Among the four covariates, three covariate designs exhibit shorter confidence intervals of our estimates relative to rdrobust. The confidence intervals were larger among the Government Worker Percentage design, however, our estimates capture more distinct heterogeneity that the higher government worker is related to higher incumbent margin.

figure[figure omitted — 1,065 chars of source]

Conclusion

We document that the existing bandwidth selectors are suboptimal when they are used for a multivariate RD design when they take the distance from a boundary point as the running variable. We further provide an alternative estimator for a multivariate RD design to estimate the heterogeneous treatment effects. In numerical simulations, we demonstrate the favorable performance of our estimator against a frequently used rdrobust procedure with the distance from a point as the scalar running variable. We apply our estimator to the study of Londono-Velez.Rodriguez.Sanchez2020 who study the impact of a scholarship program that has two eligibility requirements and a quasi-multivariate design for Lee2008 dataset with a baseline covariate to study the heterogeneous effects across the covariate values. In these application, our estimates are consistent with the original estimates, often produce shorter confidence intervals, and reveal a richer heterogeneity in the program impacts over the policy boundary than the original estimates.

Hence, we contribute to the RD estimation literature in two ways. We provide a detailed argument that the distance approach is suboptimal for a multivariate design and we provide a remedy for the problem with a dimension-specific bandwidths selector. Combined with the recent work by cattaneoEstimationInferenceBoundary2025 which documents another problem of the distance approach for designs with a corner or kink and provides an alternative estimator with a uniform inference, we provide the reason why the distance from a boundary point should not be used for a multivariate RD design to estimate heterogeneous effects across the boundary as well as an appropriate estimator to remedy the estimation problem.

Some theoretical and practical issues remain. First, our consideration is limited to a random sample; hence, spatial RD designs are excluded from our consideration. We defer our focus to spatial design because of its theoretical and conceptual complexity. Nevertheless, we aim to propose a spatial RD estimation based on newly developed asymptotic results of Kurisu.Matsuda2022 in a separated study. Second, our theoretical results can be applied to any finite-dimensional RD design; however, the practical performance of estimators with more than two dimensions is limited. Third, our approach requires a sufficiently large sample over the boundary, and its performance with an extremely small sample size is limited. For a smaller sample, an explicit randomization approach is a compelling alternative. \citet*{Cattaneo.Frandsen.Titiunik2015}, \citet*{Cattaneo.Titiunik.Vazquez-Bare2016} and \citet*{Cattaneo.Titiunik.Vazquez-Bare2017} propose the concepts and a randomization inference. Their approach requires a substantially stronger assumption but is applicable to a geographical RD design as well \citep*{Keele.Titiunik.Zubizarreta2015}. Fourth, covariates are often incorporated in the estimation procedures in RD designs. For the efficiency gain, Frolich.Huber2019 propose a method with a multi-dimensional non-parametric estimation; Calonico.Cattaneo.Farrell.Titiunik2019 develop an easy-to-implement augmentation; and recently \citet*{Noack.Olma.Rothe2021} considers flexible and efficient estimation including machine-learning devices and several studies such as \citet*{Kreiss.Rothe2021} and \citet*{Arai.Otsu.Seo2021} explore augmentation with high-dimensional covariates. We defer these analyses to theoretical and conceptual complications for a companion study for a geographic RD design. Fifth, we provided the optimal bandwidths for multivariate RD estimation; however, the optimal kernel for this class of estimators is unknown. Exploring the optimal kernel for a multivariate estimator is a topic for future research. Finally, we do not provide any procedure to aggregate heterogeneous estimates over the set of boundary points. For example, a major feature of the rdmulti package, \citet*{Cattaneo.Titiunik.Vazquez-Bare2020}, is averaging over multiple boundary points; \citet*{Cattaneo.Keele.Titiunik.Vazquez-Bare2016} offers a target pooling parameter; and \citet*{Cattaneo.Keele.Titiunik.Vazquez-Bare2021} uses a different policy in Columbia with multiple cutoffs to extrapolate the missing part of the support. These ideas can be a benchmark to consider averaging and extrapolation when the support has holes in the boundary.