The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
88,886 characters
rd2d: Causal Inference in Boundary Discontinuity Designs
\maketitle
\begin{abstract}
Boundary discontinuity designs---also known as Multi-Score Regression Discontinuity (RD) designs, with Geographic RD designs as a prominent example---are often used in empirical research to learn about causal treatment effects along a continuous assignment boundary defined by a bivariate score. This article introduces the \texttt{R} package \texttt{rd2d}, which implements and extends the methodological results developed in \cite{Cattaneo-Titiunik-Yu_2025_BoundaryRD} for boundary discontinuity designs. The package employs local polynomial estimation and inference using either the bivariate score or a univariate distance-to-boundary metric. It features novel data-driven bandwidth selection procedures, and offers both pointwise and uniform estimation and inference along the assignment boundary. The numerical performance of the package is demonstrated through a simulation study.
\end{abstract}
\textit{Keywords}: treatment effects, regression discontinuity designs, nonparametric regression.
\section{Introduction}\label{sec: Introduction}
Regression Discontinuity (RD) designs are commonly used for treatment effect estimation and causal inference in quantitative sciences \citep[see][and references therein]{Cattaneo-Titiunik_2022_ARE}. In their canonical form, each unit $i\in \{1,2,\cdots,n\}$ is assigned to control ($T_i=0$) or treatment ($T_i=1$) according to the discontinuous rule $T_i = \mathds{1}(X_i \geq c)$, where $X_i$ denotes a scalar score variable, $c$ denotes a scalar cutoff, and $\mathds{1}(\cdot)$ is the indicator function. The key idea underlying all RD designs is that units with a score near the cutoff determining treatment assignment are comparable in terms of all pretreatment observables and unobservable characteristics, the only difference being that some units are assigned to control ($X_i < c$) while other are assigned to treatment ($X_i \geq c$). Therefore, in the absence of score manipulation, units having a score near (but on different sides of) the cutoff can be used as counterfactual of each other to learn about causal treatment effects.
Boundary discontinuity designs generalize the canonical RD design to allow for a multi-dimensional score variable \citep{Papay-Willett-Murnane_2011_JoE,Reardon-Robinson_2012_JREE,Keele-Titiunik_2015_PA}. The most common case is a bivariate score $\mathbf{X}_i = (X_{1i},X_{2i})^\top$ with an assignment boundary curve $\mathcal{B}$ on its support. For example, \cite{LondonoVelezRodriguezSanchez_2020_AEJ} study the effect of a Colombian social program where $X_{1i}$ is a poverty index and $X_{2i}$ is an exam grade, for each student $i$, and where eligibility for receiving the treatment required either a minimum poverty index ($X_{1i} \geq c_1$) or a minimum exam grade ($X_{2i} \geq c_2$). Therefore, in their application, the boundary determining treatment is $\mathcal{B}=\{(X_{1i},X_{2i}): \{X_{1i} \geq c_1 \text{ and } X_{2i} = c_2\} \cup \{X_{1i}=c_1 \text{ and }X_{2i}\geq c_2\}\}$. This type of bivariate RD designs, where each score component has its own cutoff, is illustrated in Figure \ref{fig:fig1a}. Another prototypical class of boundary discontinuity designs are the Geographic RD designs: for instance, \cite{Keele-Titiunik_2015_PA} study the effect of political advertisements on voter turnout during a presidential campaign by leveraging sharp discontinuities in exposure to presidential ads induced across geographic media market boundaries. Figure \ref{fig:fig1b} illustrates a generic example of a Geographic RD design employing the boundary separating two US states. \cite{jardim2024local} gives another recent empirical application of a boundary discontinuity design in the context of labor markets, and provides further references. See \cite{Cattaneo-Idrobo-Titiunik_2020_CUP,Cattaneo-Idrobo-Titiunik_2024_CUP} for a two-part practical introductory monograph.
\begin{figure}
\centering
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{inputs/fig1a.pdf}
\caption{Two-Score RD Design.}
\label{fig:fig1a}
\end{subfigure}
\quad
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{inputs/fig1b.pdf}
\caption{Geographic RD Design.}
\label{fig:fig1b}
\end{subfigure}
\caption{Two Examples of Boundary Discontinuity Designs}
\label{fig:fig1}
\end{figure}
While classical RD designs based on a scalar score are well-understood in the literature, boundary discontinuity designs are surprisingly less studied. Methodological developments have lagged empirical practice for a while, leading to different approaches in practice, but without a foundational understanding of their relative merits. \cite{Cattaneo-Titiunik-Yu_2025_BoundaryRD} address this gap in the literature by studying the properties of two leading approaches often used in empirical research leveraging boundary discontinuity designs:
\begin{itemize}
\item the \textit{location-based approach} employs bivariate local polynomial regression analysis based on the bivariate location score relative to each point on the boundary $\mathcal{B}$; and
\item the \textit{distance-based approach} employs univariate local polynomial regression analysis based on a scalar score constructed as distance to each point on the boundary $\mathcal{B}$.
\end{itemize}
\cite{Cattaneo-Titiunik-Yu_2025_BoundaryRD} study the two methodologies, and establish novel identification, estimation and inference results, both pointwise for each point on $\mathcal{B}$ and uniformly over $\mathcal{B}$. Importantly, they demonstrate that the distance-based approach can exhibit a large bias near kinks or other irregularities in the assignment boundary $\mathcal{B}$, while the location-based approach remains valid even in those cases. Based on their findings, it is recommended to employ the bivariate location-based approach whenever possible, but their results provide foundational theoretical guidance for both empirical approaches.
This article introduces the \texttt{R} software package \texttt{rd2d}, which expands and implements the methodological
results in \cite{Cattaneo-Titiunik-Yu_2025_BoundaryRD}, thereby offering data-driven general-purpose methods for the analysis and interpretation of boundary discontinuity designs. The package includes the following four functions.
\begin{itemize}
\item \texttt{rd2d()}. This function implements location-based local polynomial regression analysis for estimation and inference of causal treatment effects in boundary discontinuity designs. For the $n$ units of analysis, the function takes as inputs their outcomes $\mathbf{Y} = (Y_i: i=1,\cdots,n)^\top$, their bivariate location scores $\mathbf{X} = [\mathbf{X}_i = (X_{1i},X_{2i})^\top:i=1,\cdots,n]^\top$, their treatment assignment indicators $\mathbf{T} = (T_i: i=1,\cdots,n)^\top$, and a collection of cutoffs $\mathbf{b} = [\mathbf{b}_j = (b_{1j},b_{2j})^\top : j=1,\cdots,J]^\top$ on the boundary $\mathcal{B}$ determining treatment assignment. The function then implements estimation and (robust bias-corrected) inference via bivariate local polynomial regression, both pointwise and uniformly over the cutoffs $\mathbf{b}_j\in\mathcal{B}$. As it is customary in nonparametric regression settings, the function requires specifying a bandwidth (localization) parameter: if not provided by the user, it is selected via the companion function \texttt{rdbw2d()} for data-driven bandwidth selection.
\item \texttt{rdbw2d()}. This function employs mean square error (MSE) approximations to implement (approximate) MSE-optimal bandwidth selection for treatment effect estimation and inference in boundary discontinuity designs. It provides second-generation direct plug-in (DPI) rules \citep{hardle2004nonparametric,wand1994kernel}, incorporating several regularization schemes.
\item \texttt{rd2d.dist()}. This function implements distance-based local polynomial regression analysis for estimation and inference of causal treatment effects in boundary discontinuity designs. For the $n$ units of analysis, the function takes as inputs their outcomes $\mathbf{Y} = (Y_i: i=1,\cdots,n)^\top$, their scalar distance scores $\mathbf{D} = [[D_i(\mathbf{b}_j): j=1,\cdots,J]^\top:i=1,\cdots,n]^\top$, and a collection of cutoffs $\mathbf{b} = [\mathbf{b}_j = (b_{1j},b_{2j})^\top : j=1,\cdots,J]^\top$ on the boundary $\mathcal{B}$. That is, $T_i(\mathbf{b}_j) = \mathds{1}(D_i(\mathbf{b}_j) \geq 0)$ denotes the treatment assignment for unit $i$ relative to cutoff $\mathbf{b}_j\in\mathcal{B}$. This function also implements estimation and (robust bias-corrected) inference via univariate local polynomial regression, both pointwise and uniformly over the cutoffs $\mathbf{b}_j \in \mathcal{B}$. If the bandwidth (localization) parameter is not provided by the user, then it is chosen via the companion function \texttt{rdbw2d.dist()} for data-driven bandwidth selection.
\item \texttt{rdbw2d.dist()}. This function implements rule-of-thumb (ROT) bandwidth selection rules, depending on the specific assumptions on $\mathcal{B}$ imposed. More precisely, depending on whether the assignment boundary $\mathcal{B}$ is assumed to be smooth or not, different ROT bandwidth selectors are implemented as supported by the underlying theoretical results in \cite{Cattaneo-Titiunik-Yu_2025_BoundaryRD}. Unfortunately, when $\mathcal{B}$ exhibits kinks or other irregularities, it is difficult to develop MSE-optimal bandwidth selection, in which case the function reverts back to a simple ROT implementation based on a rate-optimality criteria.
\end{itemize}
The methods \texttt{print()} and \texttt{summary()} are supported for objects returned by \texttt{rd2d()}, \texttt{rdbw2d()}, \texttt{rd2d.dist()}, and \texttt{rdbw2d.dist()}. We also demonstrate how to use the outputs to generate useful plots for empirical work, depicting treatment effect estimation, confidence intervals, and confidence bands, along the treatment assignment boundary $\mathcal{B}$.
In addition, the four functions in the package \texttt{rd2d} offer several practically relevant features, including (i) heteroskedasticity-robust and cluster-robust variance estimation, (ii) explicit regularization for the presence of mass points in the bivariate location score $\mathbf{X}_i$ or the univariate distance score $D_i(\mathbf{b}_j)$, and (iii) explicit regularization accounting for specific extreme shape features of the underlying unknown conditional expectation functions for bandwidth selection. In the case of \texttt{rd2d()}, the bivariate location-based approach, the data-driven point estimator is approximately MSE-optimal, while in the case of \texttt{rd2d.dist()}, the univariate distance-based approach, the rate-optimality of the data-driven point estimator depends on the underlying geometry of the assignment boundary $\mathcal{B}$. Putting aside the induced bias by the possibly non-smooth $\mathcal{B}$, both functions offer pointwise (for each $\mathbf{b}_j\in\mathcal{B}$) and uniform (over $\mathcal{B}$) robust bias-corrected inference \citep{Calonico-Cattaneo-Farrell_2018_JASA,Calonico-Cattaneo-Farrell_2022_Bernoulli}.
The main contribution of this article is to introduce and discuss the first general-purpose software implementation of (MSE-optimal) treatment effect estimation and (pointwise and uniform) uncertainty quantification methods for boundary discontinuity designs, given by the \texttt{R} package \texttt{rd2d}. To this end, the article develops second-generating DPI rules for bandwidth selection, along with principled regularization schemes for specific empirically relevant settings (e.g., mass points in $\mathbf{X}_i$) and practically relevant variance estimators (e.g., cluster-robust).
The rest of the paper proceeds as follows. Section \ref{sec: Methods and Implementation} reviews the main methodological contributions in \cite{Cattaneo-Titiunik-Yu_2025_BoundaryRD}, and presents additional results related to bandwidth selection and regularized implementation. Section \ref{sec: Numerical Illustration} demonstrates the performance of the \texttt{R} package \texttt{rd2d} using a simulation study with a data generated process calibrated using the data in \cite{LondonoVelezRodriguezSanchez_2020_AEJ}. Section \ref{sec: Conclusion} concludes. Replication codes, background references, and other information related the software package \texttt{rd2d} can be found at: \url{https://rdpackages.github.io/rd2d/}.
\section{Methods and Implementation}\label{sec: Methods and Implementation}
We employ standard potential outcomes notation. Suppose that $(Y_i(0), Y_i(1), \mathbf{X}_i^\top)^\top$, $i = 1,2,\dots, n$, is a random sample, where $Y_i(0)$ and $Y_i(1)$ denote the scalar potential outcomes for unit $i$ under control and treatment assignment, respectively. Units are assigned to control group or treatment group according to their bivariate location score $\mathbf{X}_i= (X_{1i},X_{2i})^\top$ relative to a known one-dimensional boundary $\mathcal{B}$ splitting the support of $\mathbf{X}_i$ in two disjoint regions: $\mathcal{A}_0$ denotes the control region, and $\mathcal{A}_1$ denotes the treatment region. Thus, $\mathcal{B} = \mathtt{bd}(\mathcal{A}_0) \cap \mathtt{bd}(\mathcal{A}_1)$, where $\mathtt{bd}(\mathcal{A}_t)$ denotes the boundary of the set $\mathcal{A}_t$. The observed response variable is $Y_i = (1-T_i) \cdot Y_i(0) + T_i \cdot Y_i(1)$, where $T_i = \mathds{1}(\mathbf{X}_i \in \mathcal{A}_1)$. Without loss of generality, we assume that the boundary belongs to the treatment group, that is, $\mathtt{bd}(\mathcal{A}_1)\subset\mathcal{A}_1$ and $\mathcal{B}\cap\mathcal{A}_0=\emptyset$. Figure \ref{fig:fig1} gives two graphical examples.
The causal parameter of interest is the \textit{average treatment effect curve along the boundary}:
\begin{align*}
\tau(\mathbf{x}) = \mathbb{E}[Y_i(1) - Y_i(0) | \mathbf{X}_i = \mathbf{x}],
\qquad \mathbf{x}\in\mathcal{B}.
\end{align*}
Identification follows directly by the usual continuity assumptions invoked in canonical RD designs \citep{Hahn-Todd-vanderKlaauw_2001_ECMA}:
\begin{align*}
\tau(\mathbf{x}) = \lim_{\mathbf{u}\to\mathbf{x},\mathbf{u}\in\mathcal{A}_1} \mu_1(\mathbf{u}) - \lim_{\mathbf{u}\to\mathbf{x},\mathbf{u}\in\mathcal{A}_0} \mu_0(\mathbf{u}),
\qquad \mathbf{x}\in\mathcal{B},
\end{align*}
where $\mu_t(\mathbf{x})=\mathbb{E}[Y_i(t)|\mathbf{X}_i = \mathbf{x}]=\mathbb{E}[Y_i|\mathbf{X}_i = \mathbf{x}, T_i=t]$, $t\in\{0,1\}$, are assumed to be smooth functions. See \cite{Papay-Willett-Murnane_2011_JoE}, \cite{Reardon-Robinson_2012_JREE}, \cite{Keele-Titiunik_2015_PA}, and references therein.
For implementation, the continuous assignment boundary $\mathcal{B}$ is first discretized into $J$ cutoff points $\mathbf{b}=(\mathbf{b}_1,\cdots,\mathbf{b}_J)^\top$ with $\mathbf{b}_j\in\mathcal{B}$ for all $j=1,\dots,J$. Then, the empirical analysis is conducted pointwise for each cutoff or uniformly over all cutoffs, employing either the bivariate location score $\mathbf{X}_i$ directly, or an induced univariate distance to each cutoff point. See \citet[Section 5]{Cattaneo-Idrobo-Titiunik_2024_CUP} for an introductory discussion.
Following the methodological recommendations in \cite{Cattaneo-Titiunik-Yu_2025_BoundaryRD}, most of the discussion focuses on location-based methods via bivariate local polynomial regression based on the data $(Y_i,\mathbf{X}_i^\top)^\top$, $i=1,\cdots,n$, which are implemented in the functions \texttt{rd2d()} and \texttt{rdbw2d()}. However, given their predominance in empirical work, Section \ref{sec: Distance-Based Methods} also discusses distance-based methods via univariate local polynomial regression, which require a user-chosen scalar distance score to each cutoff point on the assignment boundary, and are implemented in the functions \texttt{rd2d.dist()} and \texttt{rdbw2d.dist()}. We omit assumptions and other technical details in the remaining of the paper, which can be found in the references given.
\subsection{Location-Based Methods}
The location-based treatment effect curve estimator of $\tau(\mathbf{x})$ is
\begin{align*}
\widehat{\tau}(\mathbf{x}) = \mathbf{e}_1^{\top}\widehat{\boldsymbol{\beta}}_1(\mathbf{x}) - \mathbf{e}_1^{\top}\widehat{\boldsymbol{\beta}}_0(\mathbf{x}),
\qquad \mathbf{x} \in \mathcal{B},
\end{align*}
where, for $t \in \{0,1\}$,
\begin{align*}
\widehat{\boldsymbol{\beta}}_t(\mathbf{x})
= \operatorname*{arg\,min}_{\boldsymbol{\beta} \in \mathbb{R}^{\mathfrak{p}_p+1}} \sum_{i=1}^n \big(Y_i - \mathbf{R}_p(\mathbf{X}_i - \mathbf{x})^{\top}\boldsymbol{\beta} \big)^2 K_h(\mathbf{X}_i - \mathbf{x})\mathds{1}(\mathbf{X}_i \in \mathcal{A}_t),
\end{align*}
with $\mathfrak{p}_p = (2+p)(1+p)/2-1$, $\mathbf{R}_p(\mathbf{u}) = (1,u_1,u_2,u_1^2,u_2^2,u_1u_2,\cdots,u_1^p,u_2^p)^\top$ denotes the $p$th order polynomial expansion of the bivariate vector $\mathbf{u}=(u_1,u_2)^\top$, $K_h(\mathbf{u})=K(u_1/h,u_2/h)/h^2$ for a bivariate kernel function $K(\cdot)$, and a bandwidth parameter $h$.
In practice, it is often important to first standardize each dimension of the bivariate location score $\mathbf{X}_i = (X_{1i},X_{2i})^\top$, and in some applications it may also be useful to allow for different bandwidths for control and treatment groups. Thus, the function \texttt{rd2d()} allows for four different bandwidths: $h_{01}$ is used for $X_{1i}$ and control units, $h_{02}$ is used for $X_{2i}$ and control units, $h_{11}$ is used for $X_{1i}$ and treatment units, and $h_{02}$ is used for $X_{2i}$ and treatment units. The data-driven bandwidth selection function \texttt{rdbw2d()} also allows for both standardization of each dimension of $\mathbf{X}_i$ (via the option \texttt{stdvar=TRUE}) and different bandwidth selection for control and treatment regions (via the options \texttt{bwselect="msetwo"} or \texttt{bwselect="imsetwo"}). Thus, the function \texttt{rdbw2d()} can report up to four distinct estimated bandwidths: $(\widehat{h}_{01},\widehat{h}_{02},\widehat{h}_{11},\widehat{h}_{11})$ corresponding to $(h_{01},h_{02},h_{11},h_{11})$. The discussion in this article focuses on a single common bandwidth $h$ for simplicity, but we explain how different MSE-optimal bandwidth are estimated as appropriate. In addition, while the notation does not explicitly reflect clustering, which is a common feature in geographic and other multidimensional RD designs, the package \texttt{rd2d} allows for cluster-robust inference as explained below.
\subsubsection{Point Estimation, MSE Expansions, and Bandwidth Choices}\label{sec: Point Estimation and MSE Expansions}
Under minimal regularity conditions, the bivariate location-based treatment effect estimator $\widehat{\tau}(\mathbf{x})$ is pointwise and uniform consistent for the treatment effect curve $\tau(\mathbf{x})$ along the assignment boundary: $\widehat{\tau}(\mathbf{x}) \to_\mathbb{P} \tau(\mathbf{x})$ for each $\mathbf{x}\in\mathcal{B}$, and $\sup_{\mathbf{x} \in \mathcal{B}} |\widehat{\tau}(\mathbf{x}) - \tau(\mathbf{x})| \to_\mathbb{P} 0$, where $\to_\mathbb{P}$ denotes convergence in probability as $h\to0$ and $nh^2\to\infty$. Furthermore, under similar regularity conditions, precise (conditional) MSE expansions can be established along $\mathcal{B}$, which can then be used for principled bandwidth selection.
The pointwise (conditional) MSE expansion is
\begin{align*}
\mathbb{E} [(\widehat{\tau}(\mathbf{x}) - \tau(\mathbf{x}))^2 | \mathbf{X}]
\approx_\mathbb{P} \frac{1}{n h^2} V_{\mathbf{x}} + h^{2(p+1)} B_{\mathbf{x}}^2,
\end{align*}
where $V_{\mathbf{x}} = V_{1,\mathbf{x}} + V_{0,\mathbf{x}}$ and $B_{\mathbf{x}} = B_{1,\mathbf{x}} - B_{0,\mathbf{x}}$ denote the fixed-$n$ conditional variance and the leading conditional bias of the treatment effect estimator, respectively, and $\approx_\mathbb{P}$ denotes equality in probability up to vanishing higher-order terms. More precisely, using standard multi-index notation and least squares algebra, the fixed-$n$ conditional variance and leading conditional bias for each group $t\in \{0,1\}$ are
\begin{align*}
V_{t,\mathbf{x}} = \mathbf{e}_{1}^{\top} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} \boldsymbol{\Sigma}_{t,\mathbf{x},\mathbf{x}} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1}\mathbf{e}_{1}
\qquad\text{and}\qquad
B_{t,\mathbf{x}}
= \sum_{|\mathbf{k}| = p + 1} \frac{\mu_t^{(\mathbf{k})}(\mathbf{x})}{\mathbf{k}!} \mathbf{e}_1^{\top} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} \widehat{\boldsymbol{\vartheta}}_{t,\mathbf{x}}(\mathbf{k}),
\end{align*}
respectively, where
\begin{align*}
\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}
&= \frac{1}{n} \sum_{i=1}^n \mathbf{r}_p\big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\big) \mathbf{r}_p \big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\big)^{\top} K_h(\mathbf{X}_i - \mathbf{x}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t),\\
\boldsymbol{\Sigma}_{t,\mathbf{x},\mathbf{x}}
&= \frac{h^2}{n} \sum_{i=1}^n \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big) \mathbf{r}_p\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)^{\top} K_h(\mathbf{X}_i - \mathbf{x})^2 \varepsilon_{i,t,\mathbf{x}}^2 \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t),\\
\widehat{\boldsymbol{\vartheta}}_{t,\mathbf{x}}(\mathbf{k})
&= \frac{1}{n} \sum_{i=1}^n \mathbf{r}_p \Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)\Big(\frac{\mathbf{X}_i - \mathbf{x}}{h}\Big)^{\mathbf{k}}K_h(\mathbf{X}_i - \mathbf{x}) \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t),
\end{align*}
and $\varepsilon_{i,t,\mathbf{x}} = Y_i - \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) \mu_t(\mathbf{x})$.
Therefore, for each $\mathbf{x}\in\mathcal{B}$, and noting that $\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}$, $\boldsymbol{\Sigma}_{t,\mathbf{x},\mathbf{x}}$, and $\widehat{\boldsymbol{\vartheta}}_{t,\mathbf{x}}(\mathbf{k})$ converge in probability (as $nh^2\to\infty$) to well-defined limits independent on the bandwidth $h$, an MSE-optimal bandwidth choice is
\begin{align*}
h_{\operatorname{MSE},\mathbf{x}} = \Big(\frac{2 V_{\mathbf{x}} }{(2p+2) B_{\mathbf{x}}^2} \frac{1}{n} \Big)^{1/(2p+4)},
\end{align*}
provided that $B_{\mathbf{x}} \neq 0$.
Similarly, given a weighting function $w(\mathbf{x})$, an integrated MSE expansion is
\begin{align*}
\int_\mathcal{B} \mathbb{E} [(\widehat{\tau}(\mathbf{x}) - \tau(\mathbf{x}))^2 | \mathbf{X}] w(\mathbf{x}) d\mathbf{x}
\approx_\mathbb{P} h^{2(p+1)} \int_\mathcal{B} B_{\mathbf{x}}^2 w(\mathbf{x}) d \mathbf{x} + \frac{1}{n h^2} \int_\mathcal{B} V_{\mathbf{x}} w(\mathbf{x}) d\mathbf{x},
\end{align*}
using the notation already introduced. Therefore, along $\mathcal{B}$, an integrated MSE-optimal bandwidth choice is
\begin{align*}
h_{\operatorname{IMSE}} = \Big(\frac{2 \int_{\mathcal{B}} V_{\mathbf{x}} w(\mathbf{x})d \mathbf{x}}{(2p+2) \int_{\mathcal{B}} B_{\mathbf{x}}^2 w(\mathbf{x}) d \mathbf{x}} \frac{1}{n} \Big)^{1/(2p+4)},
\end{align*}
provided that $\int_{\mathcal{B}} B_{\mathbf{x}}^2 w(\mathbf{x}) d \mathbf{x} \neq 0$.
The basic MSE-optimal and IMSE-optimal bandwidth choices, $h_{\operatorname{MSE},\mathbf{x}}$ and $h_{\operatorname{IMSE}}$, can be extended to accommodate different selections for each coordinate of the bivariate location score and/or for control and treatment groups separately. Different bandwidths for each coordinate in $\mathbf{X}_i = (X_{1i},X_{2i})^\top$ are obtained by first standardizing each component, then applying the basic bandwidth rules, and finally removing the standardization: letting $\tilde{\mathbf{X}}_i = (X_{1i}/\sigma_{X_1},X_{2i}/\sigma_{X_2})^\top$, where $\sigma_{X_l}^2 = \mathbb{V}[X_{li}]$ for $l\in\{1,2\}$, then for each coordinate $l\in\{1,2\}$ the bandwidth selectors are
\begin{align*}
h_{l,\operatorname{MSE},\mathbf{x}} = \sigma_{X_l} \cdot h_{\operatorname{MSE},\mathbf{x}}
\qquad\text{and}\qquad
h_{l,\operatorname{IMSE}} = \sigma_{X_l} \cdot h_{\operatorname{IMSE}},
\end{align*}
where $h_{\operatorname{MSE},\mathbf{x}}$ and $h_{\operatorname{IMSE}}$ are computed using the standardized bivariate location score $\tilde{\mathbf{X}}_i$ (instead of using the original score $\mathbf{X}_i$). Different bandwidth selection for control and treatment groups are obtained by implementing $h_{\operatorname{MSE},\mathbf{x}}$ and $h_{\operatorname{IMSE}}$ for each group $t\in\{0,1\}$ separately:
\begin{align*}
h_{\operatorname{MSE},t,\mathbf{x}} = \Big(\frac{2 V_{t,\mathbf{x}} }{(2p+2) B_{t,\mathbf{x}}^2} \frac{1}{n} \Big)^{1/(2p+4)}
\quad\text{and}\quad
h_{\operatorname{IMSE},t} = \Big(\frac{2 \int_{\mathcal{B}} V_{t,\mathbf{x}} w(\mathbf{x})d \mathbf{x}}{(2p+2) \int_{\mathcal{B}} B_{t,\mathbf{x}}^2 w(\mathbf{x}) d \mathbf{x}} \frac{1}{n} \Big)^{1/(2p+4)},
\end{align*}
assuming the denominators are not zero.
A combination of the two ideas gives four distinct MSE-optimal and IMSE-optimal bandwidth selection rules: $(\sigma_{X_1} \cdot h_{\operatorname{MSE},0,\mathbf{x}},\sigma_{X_2} \cdot h_{\operatorname{MSE},0,\mathbf{x}},\sigma_{X_1} \cdot h_{\operatorname{MSE},1,\mathbf{x}},\sigma_{X_2} \cdot h_{\operatorname{MSE},1,\mathbf{x}})$ and $(\sigma_{X_1} \cdot h_{\operatorname{IMSE},0},\sigma_{X_2} \cdot h_{\operatorname{IMSE},0},\sigma_{X_1} \cdot h_{\operatorname{IMSE},1},\sigma_{X_2} \cdot h_{\operatorname{IMSE},1})$, respectively. For implementation, the package \texttt{rd2d} employs the following options:
\begin{align*}
\begin{tabular}{clclclclc}
\hline\hline
\multicolumn{3}{c}{Option} & & \multicolumn{5}{c}{Bandwidth Selection} \\ \cline{1-3} \cline{5-9}
\texttt{stdvar} & & \texttt{bwselect} & & Score Coordinates & & Treatment Groups & & Default \\ \cline{1-1} \cline{3-3} \cline{5-5} \cline{7-7} \cline{9-9}
\texttt{TRUE} & & \texttt{mserd} & & Distinct & & Same & & Yes \\
\texttt{FALSE} & & \texttt{mserd} & & Same & & Same & & No \\
\texttt{TRUE} & & \texttt{msetwo} & & Distinct & & Distinct & & No \\
\texttt{FALSE} & & \texttt{msetwo} & & Same & & Distinct & & No \\
\hline\hline
\end{tabular}
\end{align*}
To implement the bandwidth selection procedures it is necessary to (i) estimate the unconditional variances $\sigma^2_{X_l}$ for standardization of each coordinate of $\mathbf{X}_i=(\mathbf{X}_{li}:l=1,2)^\top$, as needed; (ii) estimate residuals $\varepsilon_{i,t,\mathbf{x}}$ entering the variances $V_{t,\mathbf{x}}$ for each group $t\in\{0,1\}$; (iii) estimate higher-order derivatives $(\mu_t^{(\mathbf{k})}(\mathbf{x}) : |\mathbf{k}| = p + 1)$ entering the bias $B_{t,\mathbf{x}}$ quantities for each group $t\in\{0,1\}$; and (iv) select (preliminary) bandwidth entering the matrices $\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}$, $\boldsymbol{\Sigma}_{t,\mathbf{x},\mathbf{x}}$, and $\widehat{\boldsymbol{\vartheta}}_{t,\mathbf{x}}(\mathbf{k})$ for each group $t\in\{0,1\}$. (A data-driven regularization term is also added to the denominators to avoid near-zero bias, as explained below.) These unknown quantities are estimated as follows.
\begin{itemize}
\item The unconditional variances used for standardization are estimated using their sample variance counterparts: $\widehat{\sigma}^2_{X_l} = \frac{1}{n-1} \sum_{i=1}^n (X_{li} - \frac{1}{n}\sum_{i=1}^n X_{li})^2$, for each coordinate $l\in\{1,2\}$.
\item Given a (preliminary) bandwidth, the package \texttt{rd2d} implements several variance estimators for each group $t\in\{0,1\}$: $V_{t,\mathbf{x}}$ is replaced by $\widehat{V}_{t,\mathbf{x}} = \mathbf{e}_{1}^{\top} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1} \widehat{\boldsymbol{\Sigma}}_{t,\mathbf{x},\mathbf{x}} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1}\mathbf{e}_{1}$, where $\widehat{\boldsymbol{\Sigma}}_{t,\mathbf{x},\mathbf{x}}$ is either a heteroskedasticity-consistent (HC) or cluster-consistent (CR) variance estimator based on replacing the unknown residuals $\varepsilon_{i,t,\mathbf{x}}$ with the plug-in residuals estimates
\begin{align*}
\widehat{\varepsilon}_{i,t,\mathbf{x}} = Y_i - \mathbf{R}_p(\mathbf{X}_i - \mathbf{x})^\top\widehat{\boldsymbol{\beta}}_t(\mathbf{x}).
\end{align*}
More precisely, the package allows for the following options:
\begin{align*}
\resizebox{0.9\textwidth}{!}{
\begin{tabular}{cccc}
\hline\hline
\multicolumn{2}{c}{Option} & & Variance Implementation \\ \cline{1-2}
\texttt{vce} & \texttt{cluster} & & HC-estimate (\texttt{cluster=NULL}) or CR-estimate (\texttt{cluster=}\textit{varname}) \\ \cline{1-2} \cline{4-4}
\texttt{hc0} & \texttt{NULL} or \textit{varname} & & no weighting of estimated residuals\\
\texttt{hc1} & \texttt{NULL} or \textit{varname} & & degrees-of-freedom weighting of estimated residuals\\
\texttt{hc2} & \texttt{NULL} or \textit{varname} & & inverse-diagonal-projection weighting of estimated residuals\\
\texttt{hc3} & \texttt{NULL} or \textit{varname} & & inverse-diagonal-projection-squared weighting of estimated residuals\\
\hline\hline
\end{tabular}
}
\end{align*}
See \cite{Zeileis_2004_JSS} and \cite{Zeileis-etal_2020_JSS} for more discussion.
\item Given a (preliminary) bandwidth, the package \texttt{rd2d} estimate the higher-order curvature of the unknown conditional expectations, $\mu_t^{(\mathbf{k})}(\mathbf{x})$ and $\mu_1^{(\mathbf{k})}(\mathbf{x})$ for each $\mathbf{k}\in\mathbb{N}^2$ such that $|\mathbf{k}| = p + 1$, using a higher-order polynomial approximation. By default, a local polynomial regression of order $q=p+1$ is used.
\item The preliminary bandwidth(s) needed to construct the matrices $\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}$, $\boldsymbol{\Sigma}_{t,\mathbf{x},\mathbf{x}}$, and $\widehat{\boldsymbol{\vartheta}}_{t,\mathbf{x}}(\mathbf{k})$, as well as the plug-in residuals and higher-order derivative estimates, are selected using a combination of ROT and DPI-2 methods. Specifically, two bandwidth are sequentially constructed as follows:
\begin{enumerate}
\item[Step 1.] Using a Gaussian distribution reference model, construct a plug-in IMSE-optimal ROT bandwidth selector for the canonical kernel density estimator of $f_{\mathbf{X}}(\mathbf{x})$, the Lebesgue density of score. The resulting data-driven bandwidth choice is:
\begin{align*}
\widehat{c} = \widehat{C} n^{-1/6},
\end{align*}
where $\widehat{C}$ is a function of the variance of $\mathbf{X}_i$ and known constants determined by the kernel function used. This preliminary bandwidth choice is motivated by the fact that $\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}} \to_\mathbb{P} f_{\mathbf{X}}(\mathbf{x})\boldsymbol{\Gamma}_{t,\mathbf{x}}$ and $\widehat{\boldsymbol{\vartheta}}_{t,\mathbf{x}}(\mathbf{k}) \to_\mathbb{P} f_{\mathbf{X}}(\mathbf{x})\boldsymbol{\vartheta}_{t,\mathbf{x}}(\mathbf{k})$, where $\boldsymbol{\Gamma}_{t,\mathbf{x}}$ and $\boldsymbol{\vartheta}_{t,\mathbf{x}}(\mathbf{k})$ are non-random matrices, only function of $\mathcal{B}$, $K$, $p$, $|\mathbf{k}| = p + 1$, and $t\in\{0,1\}$.
\item[Step 2.] Construct an MSE-optimal bandwidth choice for estimating the linear combination given by $\sum_{|\mathbf{k}| = p + 1} \frac{\mu_t^{(\mathbf{k})}(\mathbf{x})}{\mathbf{k}!} \mathbf{e}_1^{\top} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1}\widehat{\boldsymbol{\vartheta}}_{t,\mathbf{x}}(\mathbf{k})$, using a $q$th order local polynomial estimator and the bandwidth $\widehat{c}$ for the coefficients $\mathbf{e}_1^{\top} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}^{-1}\widehat{\boldsymbol{\vartheta}}_{t,\mathbf{x}}(\mathbf{k})$. To implement this bandwidth choice, $\widehat{c}$ is used for variance estimation, and a preliminary nearest-neighbor-based polynomial regression approximation is used for bias estimation. The resulting data-driven bandwidth choice is:
\begin{align*}
\widehat{b}_t = \widehat{C}_t n^{-1/(2q+4)},
\qquad t\in\{0,1\},
\end{align*}
where $\widehat{C}_t$ depends on the variance and bias estimates for the target linear combination, and $K$ and $q$.
\item[Step 3.] Construct the final variance and bias constants using the preliminary bandwidth estimates $(\widehat{c},\widehat{b}_0,\widehat{b}_1)$. Specifically, the basic MSE-optimal (or IMSE-optimal) bandwidth choice is implemented as follows: for each $t\in\{0,1\}$,
\begin{itemize}
\item $V_{t,\mathbf{x}}$ is replaced by $\widehat{V}_{t,\mathbf{x}}$, where $\widehat{c}$ is used to construct the matrices, and a $p$th order local polynomial regression is used for residual estimation; and
\item $B_{t,\mathbf{x}}$ is replaced by $\widehat{B}_{t,\mathbf{x}}$, where $\widehat{c}$ is used to construct the matrices, and $\widehat{b}_t$ is used to estimate the derivatives of the regression function for each group.
\end{itemize}
\end{enumerate}
See \cite{wand1994kernel} and \cite{hardle2004nonparametric} for technical details.
\end{itemize}
\subsubsection{Statistical Inference}\label{sec: Statistical Inference}
To assess uncertainty in the estimation of the causal effects $\tau(\mathbf{x})$ along the boundary $\mathcal{B}$, we consider the usual Wald-type test statistic
\begin{align*}
\operatorname{T}(\mathbf{x})
= \frac{\widehat{\tau}(\mathbf{x}) - \tau(\mathbf{x})}
{\sqrt{\widehat{V}_{\mathbf{x}}/(nh^2)}},
\qquad \mathbf{x}\in\mathcal{B},
\end{align*}
where $\widehat{\tau}(\mathbf{x})$ and $\widehat{V}_{\mathbf{x}}$ are constructed using the polynomial order $p$ and bandwidth $h$, after choosing the appropriate HC and CR variance estimator, as explained above in the context of bandwidth selection.
In practice, the bandwidth is chosen to be (I)MSE-optimal, and thus the sampling distribution of the statistic satisfies the following pointwise distributional approximation:
\begin{align*}
\operatorname{T}(\mathbf{x}) - \text{Bias}(\mathbf{x}) \overset{a}{\thicksim} \mathsf{N}(0, 1),
\qquad \text{Bias}(\mathbf{x}) = \frac{h^{p+1} B_\mathbf{x}}{\sqrt{\widehat{V}_{\mathbf{x}}/(nh^2)}},
\qquad \mathbf{x}\in\mathcal{B},
\end{align*}
where $\overset{a}{\thicksim}$ denote an approximation in distribution as $nh^2\to\infty$, $\mathsf{N}(0,1)$ denotes the standard Gaussian distribution, and $\text{Bias}(\mathbf{x})$ denotes the standardized leading bias emerging whenever a ``large'' bandwidth is used (i.e., when the (I)MSE-optimal bandwidth is used, or any other bandwidth choice such that $\text{Bias}(\mathbf{x}) \not\to_\mathbb{P} 0$). The standard confidence interval estimator with $100(1-\alpha)\%$ nominal coverage is
\begin{align*}
\operatorname{CI}(\mathbf{x})
= \left[\widehat{\tau}(\mathbf{x}) \pm \Phi_{1-\alpha/2} \sqrt{\widehat{V}_{\mathbf{x}}/(nh^2)} \right],
\qquad \mathbf{x}\in\mathcal{B},
\end{align*}
where $\Phi_\alpha$ is the $\alpha$th quantile of the standard Gaussian distribution. However, for ``large'' bandwidths such as the (I)MSE-optimal choice, $\operatorname{CI}(\mathbf{x})$ will be invalid due to the bias $\text{Bias}(\mathbf{x})$, thereby delivering empirical coverage well below its nominal target. A solution to this problem is to employ ad-hoc undersmoothing, that is, to implement $\operatorname{CI}(\mathbf{x})$ with a ``smaller'' bandwidth $h$ relative to the (I)MSE-optimal choice. \cite{Calonico-Cattaneo-Farrell_2018_JASA,Calonico-Cattaneo-Farrell_2022_Bernoulli} showed that undersmoothing is sub-optimal (possibly invalid) under standard assumptions, while the robust bias-correction (RBC) methodology introduced by \cite{Calonico-Cattaneo-Titiunik_2014_ECMA} enjoys validity and better (in some cases optimal) higher-order distributional properties. The core idea behind the RBC approach can be summarized as follows: (i) employ the (I)MSE-optimal bandwidth for constructing the point estimator $\widehat{\tau}(\mathbf{x})$, (ii) de-bias (bias correct) the numerator of the statistic $\operatorname{T}(\mathbf{x})$, and (iii) adjust the variance estimate to account for the variability introduced by the debiasing of $\widehat{\tau}(\mathbf{x})$.
The implementation of the RBC inference methodology is straightforward: given a chosen (I)MSE-optimal bandwidth for the $p$th order local polynomial point estimator $\widehat{\tau}(\mathbf{x})$, an adjsted test statistic is constructed using a $q$th order local polynomial point estimator $\widehat{\tau}_q(\mathbf{x})$ and its associated variance estimate $\widehat{V}_{\mathbf{x},q}$, with $q>p$. Thus, the approach employs the RBC statistic $T_{q}(\mathbf{x}) = (\widehat{\tau}_q(\mathbf{x}) - \tau(\mathbf{x}))/\sqrt{\widehat{V}_{\mathbf{x},q}/(nh^2)}$, instead of the original statistic $\operatorname{T}(\mathbf{x})$ above. The resulting RBC confidence interval estimator is
\begin{align*}
\operatorname{CI}_q(\mathbf{x})
= \left[\widehat{\tau}_q(\mathbf{x}) \pm \Phi_{1-\alpha/2} \sqrt{\widehat{V}_{\mathbf{x},q}/(nh^2)} \right],
\qquad \mathbf{x}\in\mathcal{B},
\end{align*}
which is constructed with an (I)MSE-optimal bandwidth choice for the $p$th order local polynomial point estimator $\widehat{\tau}(\mathbf{x})$. (This amounts to a form of robust bias correction because $\widehat{\tau}_q(\mathbf{x}) = \widehat{\tau}(\mathbf{x}) - h^{p+1} \widehat{B}_\mathbf{x}$, where $\widehat{B}(\mathbf{x})$ is an ``estimate'' of $B_\mathbf{x}$.) The package \texttt{rd2d} employs $p=1$ and $q=2$ as defaults, which are standard choices for implementation.
Uniform inference and confidence bands along the boundary $\mathcal{B}$ also employ RBC methodology. Specifically, given an (I)MSE-optimal bandwidth for the $p$th order local polynomial point estimator $\widehat{\tau}(\mathbf{x})$, the associated RBC confidence band estimate is
\begin{align*}
\operatorname{CB}_q(\mathcal{B}) =
\left\{\left[\widehat{\tau}_q(\mathbf{x}) \pm \mathsf{q}_{\alpha} \sqrt{\widehat{V}_{\mathbf{x},q}/(nh^2)} \right]: \mathbf{x} \in \mathcal{B} \right\},
\end{align*}
where $\mathsf{q}_{\alpha}$ denotes a suitably chosen quantile to control false rejections uniformly over $\mathcal{B}$. In practice, the continuous assignment boundary is discretized to consider the $J$ cutoff in $\mathbf{b} = (\mathbf{b}_1,\cdots,\mathbf{b}_J)^\top$ jointly. Then, a feasible quantile choice is
\begin{align*}
\mathsf{q}_{\alpha}
= \inf\left\{ u\geq 0:\ \mathbb{P}\Big[ \max_{1\leq j \leq J} \big|\mathbf{e}_j^\top\widehat{\mathbf{C}}_{\mathbf{b},q}^{1/2} \mathbf{Z} \big| \geq u \Big| \text{Data} \Big] \leq \alpha \right\},
\end{align*}
where the $J$-dimensional standard Gaussian vector $\mathbf{Z}\thicksim\mathsf{N}(\mathbf{0}_J, \mathbf{I}_J)$ is independent of the data, and the $(J\times J)$ covariance matrix $\widehat{\mathbf{C}}_{\mathbf{b},q} = \big[\widehat{V}_{\mathbf{b}_j,\mathbf{b}_k} / (\widehat{V}_{\mathbf{b}_j,\mathbf{b}_j} \widehat{V}_{\mathbf{b}_k,\mathbf{b}_k})^{1/2}: 1\leq j,k \leq J \big]$ where $\widehat{V}_{\mathbf{b}_j,\mathbf{b}_k} = \widehat{V}_{1,\mathbf{b}_j,\mathbf{b}_k} + \widehat{V}_{0,\mathbf{b}_j,\mathbf{b}_k}$ with
\begin{align*}
V_{t,\mathbf{b}_j,\mathbf{b}_k}
= \mathbf{e}_{1}^{\top} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{b}_j}^{-1} \widehat{\boldsymbol{\Sigma}}_{t,\mathbf{b}_j,\mathbf{b}_k} \widehat{\boldsymbol{\Gamma}}_{t,\mathbf{b}_k}^{-1}\mathbf{e}_{1}
\end{align*}
and
\begin{align*}
\widehat{\boldsymbol{\Sigma}}_{t,\mathbf{b}_j,\mathbf{b}_k}
= \frac{h^2}{n} \sum_{i=1}^n \mathbf{r}_q \Big(\frac{\mathbf{X}_i - \mathbf{b}_j}{h}\Big) \mathbf{r}_q \Big(\frac{\mathbf{X}_i - \mathbf{b}_k}{h}\Big)^{\top} K_h(\mathbf{X}_i - \mathbf{b}_j) K_h(\mathbf{X}_i - \mathbf{b}_k) \widehat{\varepsilon}_{i,\mathbf{b}_j}\widehat{\varepsilon}_{i,\mathbf{b}_k} \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t)
\end{align*}
with $\widehat{\varepsilon}_{i,t,\mathbf{x}} = Y_i - \mathbf{R}_q(\mathbf{X}_i - \mathbf{x})^\top\widehat{\boldsymbol{\beta}}_t(\mathbf{x})$,
for $t\in\{0,1\}$. This formulation is robust to unknown conditional heteroskedasticity, while the clustered-robust analogue formula is omitted to save space; see \cite{Zeileis_2004_JSS} and \cite{Zeileis-etal_2020_JSS}.
The RBC method produces confidence intervals/bands that are not centered at the treatment effect point estimator because different polynomial orders are used for estimation and inference. As a result, the point estimates may lie outside the RBC confidence intervals/bands, particularly if the underlying treatment effect curve $\tau(\mathbf{x})$ exhibits high curvature at certain evaluation points. One possible solution is to increase the polynomial orders $p$ and $q$, or to use a bandwidth smaller than the (I)MSE-optimal one.
\subsubsection{Regularization Strategies}\label{sec: Implementation and Regularization Details}
The package \texttt{rd2d} implements several regualization schemes to ensure robustness in applications.
\begin{itemize}
\item \textit{Small bias regularization}. Ignoring the asymptotically constant and higher-order terms, the approximate MSE-optimal and IMSE-optimal bandwidth choices require $B_{\mathbf{x}}\neq0$ and $\int_{\mathcal{B}} B_{\mathbf{x}}^2 \omega(\mathbf{x})d\mathbf{x} \neq 0$, respectively. Thus, a small estimated bias can result in a bandwidth that is too large. To avoid this problem, a regularization term is added to the term of estimated bias, leading to the regularized MSE-optimal bandwidth choice,
\begin{align*}
h_{\operatorname{MSE},\mathbf{x}}
= \Big(\frac{2 \widehat{V}_{\mathbf{x}} }{(2p+2) (\widehat{B}_{\mathbf{x}}^2 + s\cdot \mathbb{V}[\widehat{B}_{\mathbf{x}}]) } \frac{1}{n} \Big)^{1/(2p+4)}, \qquad \mathbf{x} \in \mathcal{B},
\end{align*}
and the regularized IMSE-optimal bandwidth choice,
\begin{align*}
h_{\operatorname{IMSE}} = \Big(\frac{2 \int_{\mathcal{B}} \widehat{V}_{\mathbf{x}} d \mathbf{x}}{(2p+2) (\int_{\mathcal{B}} \widehat{B}_{\mathbf{x}}^2 + s \cdot \mathbb{V}[\widehat{B}_{\mathbf{x}}] d \mathbf{x}) } \frac{1}{n} \Big)^{1/(2p+4)}.
\end{align*}
where the regularization terms account for variance of the bias estimator, and are estimated as discussed previously. The factor $s$, defaulted to $3$, controls the degree of regularization
\item \textit{Minimum sample size}. A sample size of at least \texttt{bwcheck} is required by (possibly) enlarging the selected/provided bandwidth until \texttt{bwcheck} number of observations are included in the estimation region. The default is $\text{\texttt{bwcheck}} = 50 + (2+p)(1+p)/2-1$. When kernel type is \texttt{"prod"}, a smallest rectangle centered at the evaluation point with two edges proportional to $(\sigma(X_{i1}), \sigma(X_{i2}))$, $\sigma$ stands for the standard deviation, is found, and the bandwidth for local polynomial fitting should allow the smallest rectangle to be contained in its resulting kernel. When kernel type is \texttt{"rad"}, a smallest ball centered at evaluation point with \texttt{bwcheck} number of data points is found, and the bandwidth is increased until its resulting kernel contains the smallest ball.
\item \textit{Mass points in $\mathbf{X}_i$}. The \texttt{masspoint} option checks for unique number of points in the data. The default is \texttt{masspoint = "check"}, where unique number of data points is reported, and a warning is issued if duplication exceeds 20\% of the data. When \texttt{masspoint = "adjust"}, bandwidths are regularized so that the resulting kernels contain a minimal number of unique observations (see \textit{minimum sample size}). When \texttt{masspoint = "off"}, the potential presence of mass points is ignored.
\end{itemize}
\subsection{Aggregated Average Treatment Effects Along the Boundary}
For a weight function $w: \mathcal{B} \rightarrow [0,\infty)$, the \textit{aggregated average treatment effect} (AATE) along the boundary $\mathcal{B}$ is
\begin{align*}
\tau_{\text{AATE},\mathcal{B}} = \frac{\int_\mathcal{B} \tau(\mathbf{b}) w(\mathbf{b}) \, d\mathbf{b}}{\int_\mathcal{B} w(\mathbf{b}) \, d\mathbf{b}}
\end{align*}
A choice of weight function recovers a scalar causal effect, aggregating potential heterogeneous treatment effects along the boundary. Estimation and inference methods for this class of causal parameters can be deduced from our results. For example, consider the generic plug-in ``estimator'' of $\tau_{\text{AATE},\mathcal{B}}$ given by
\begin{align*}
\widehat{\tau}_{\text{AATE},\mathbf{b}} = \frac{\sum_{j=1}^J \widehat{\tau}(\mathbf{b}_j) w(\mathbf{b}_j)}{\sum_{j=1}^J w(\mathbf{b}_j)},
\end{align*}
where the cutoff points $\mathbf{b}=(\mathbf{b}_1,\dots,\mathbf{b}_J)^\top$ along the assignment boundary are assumed to be evenly-spaced over $\mathcal{B}$, and $\widehat{\tau}(\mathbf{b}_j)$ is constructed using $p$th order local polynomial fits. It follows that the IMSE-optimal bandwidth choice is more natural. Feasible RBC confidence intervals are
\begin{align*}
\operatorname{CI}_{q}(\text{AATE})
= \bigg[\widehat{\tau}_{\text{AATE},\mathbf{b},q} \pm \Phi_{1-\alpha/2}
\sqrt{\frac{\mathbf{w}^\top \widehat{\mathbf{V}}_{\mathbf{b},q} \mathbf{w}}{n h^2}}\, \bigg],
\end{align*}
where $\widehat{\tau}_{\text{AATE},\mathbf{b},q}$ is the point estimator $\widehat{\tau}_{\text{AATE},\mathbf{b}}$ constructed using a $q$th order polynomial basis, $\widehat{\mathbf{V}}_{\mathbf{b},q} = \big[\widehat{V}_{\mathbf{b}_j,\mathbf{b}_k} : 1\leq j,k \leq J \big]$ denotes its associated variance estimator, as explained previously, and $\mathbf{w} = (w(\mathbf{b}_1), \dots, w(\mathbf{b}_J))^\top/(\sum_{j = 1}^J w(\mathbf{b}_j))$ is the user-chosen vector of weights. In \texttt{rd2d}, AATE estimation and inference are implemented via the optional argument \texttt{AATE = w} in the \texttt{summary()} method for \texttt{rd2d} objects.
\subsection{Distance-Based Methods}\label{sec: Distance-Based Methods}
For each unit $i=1,\ldots,n$, their scalar distance-based score to the boundary point $\mathbf{x}=(x_1,x_2)^\top\in\mathcal{B}$ is
$D_i(\mathbf{x}) = (2T_i-1)\mathcal{d}(\mathbf{X}_i, \mathbf{x})$, where $\mathcal{d}(\cdot,\cdot)$ denotes a distance function such as the Euclidean distance $\mathcal{d}(\mathbf{X}_i, \mathbf{x}) = \|\mathbf{X}_i - \mathbf{x}\| = \sqrt{(X_{1i} - x_1)^2 + (X_{2i} - x_2)^2}$. Therefore, for each $\mathbf{x}\in\mathcal{B}$, the setup reduces to a standard univariate RD design with distance score $D_i(\mathbf{x})\in\mathbb{R}$ and cutoff $c=0$, the observed data now being $(Y_1,D_1(\mathbf{x})),\cdots,(Y_n,D_n(\mathbf{x}))$ for each point on the assignment boundary $\mathbf{b}=(\mathbf{b}_1,\cdots,\mathbf{b}_J)^\top$.
The distance-based local polynomial treatment effect curve estimator of $\tau(\mathbf{x})$ is
\begin{align*}
\widehat{\tau}_{\operatorname{dis}}(\mathbf{x}) = \mathbf{e}_1^{\top} \widehat{\boldsymbol{\gamma}}_1(\mathbf{x}) - \mathbf{e}_1^{\top} \widehat{\boldsymbol{\gamma}}_0(\mathbf{x}),
\qquad \mathbf{x} \in \mathcal{B},
\end{align*}
where, for $t \in \{0,1\}$,
\begin{align*}
\widehat{\boldsymbol{\gamma}}_{t}(\mathbf{x})
= \operatorname*{arg\,min}_{\boldsymbol{\gamma} \in \mathbb{R}^{p+1}} \frac{1}{n} \sum_{i=1}^n \big(Y_i - \mathbf{r}_p(D_i(\mathbf{x}))^{\top} \boldsymbol{\gamma} \big)^2 k_h(D_i(\mathbf{x})) \mathds{1}(D_i(\mathbf{x})\in\mathcal{I}_t),
\end{align*}
with $\mathbf{r}_p(u)=(1,u,u^2,\cdots,u^p)^\top$ the usual univariate polynomial basis, $k_h(u)=k(u/h)/h^2$ for univariate kernel function $k(\cdot)$ and bandwidth parameter $h$, and $\mathcal{I}_0 = (-\infty,0)$ and $\mathcal{I}_1 = [0,\infty)$. \citet{Cattaneo-Titiunik-Yu_2025_BoundaryRD} studied the statistical properties of the distance-based approach in boundary discontinuity designs, and obtained the following main results (under regularity conditions).
\begin{enumerate}
\item \textit{Consistency}. As $h\to0$ and $nh^2\to\infty$, $\widehat{\tau}_{\operatorname{dis}}(\mathbf{x}) \to_\mathbb{P} \tau_{\operatorname{dis}}(\mathbf{x})$ for all $\mathbf{x}\in\mathcal{B}$, and $\sup_{\mathbf{x} \in \mathcal{B}} |\widehat{\tau}_{\operatorname{dis}}(\mathbf{x}) - \tau_{\operatorname{dis}}(\mathbf{x})| \to_\mathbb{P} 0$, where $\tau_{\operatorname{dis}}(\mathbf{x}) = \lim_{r\downarrow0} \theta_{1,\mathbf{x}}(r) - \lim_{r\uparrow0} \theta_{0,\mathbf{x}}(r)$ with
\begin{align*}
\theta_{t,\mathbf{x}}(r)
= \mathbb{E}\big[Y_i \big| D_i(\mathbf{x}) = r , D_i(\mathbf{x}) \in \mathcal{I}_t \big]
= \mathbb{E}\big[Y_i(t) \big| \mathcal{d}(\mathbf{X}_i,\mathbf{x}) = |r|\big],
\end{align*}
for $\mathbf{x}\in\mathcal{B}$ and $t\in\{0,1\}$. The functions $\theta_{t,\mathbf{x}}(r)$ are the univariate induced conditional expectations based on distance to the boundary point $\mathbf{x}\in\mathcal{B}$ for each group $t\in\{0,1\}$.
\item \textit{Identification}. $\tau(\mathbf{x}) = \tau_{\operatorname{dis}}(\mathbf{x})$ for all $\mathbf{x}\in\mathcal{B}$, thereby showing that the distance-based estimator is a valid treatment effect estimator.
\item \textit{Bias}. If the assignment boundary $\mathcal{B}$ is non-smooth, then the uniform bias of the estimator $\widehat{\tau}_{\operatorname{dis}}(\mathbf{x})$ along the boundary is no better than of order $h$, regardless of the polynomial order $p$ used. In other words, the distance-based estimator exhibits a ``large'' bias near kinks or other irregularities of the assignment boundary $\mathcal{B}$. On the other hand, if the assignment boundary $\mathcal{B}$ is smooth enough, then the bias of $\widehat{\tau}_{\operatorname{dis}}(\mathbf{x})$ is of order $h^{p+1}$ as expected in local polynomial regression settings.
\item \textit{Mean Square Convergence and Bandwidth Choice}. Due to unknown form of distance function $\mathcal{d}(\cdot,\cdot)$ and the assignment boundary $\mathcal{B}$, it is not possible to obtain valid (I)MSE expansions and precise bandwidth selection rules. At this level of generality, only bandwidth selection in terms of rates can be established:
\begin{itemize}
\item If $\mathcal{B}$ is smooth, then $h \asymp n^{-1/(2p+4)}$ is (I)MSE rate-optimal, where $\asymp$ denotes up to a proportional constant.
\item If $\mathcal{B}$ is non-smooth, then $h \asymp n^{-1/4}$ is (I)MSE rate-optimal, regardless of the polynomial order $p$ used in constructing $\widehat{\tau}_{\operatorname{dis}}(\mathbf{x})$.
\end{itemize}
\item \textit{Statistical Inference}. Putting aside the issue of ``large'' bias whenever $\mathcal{B}$ is non-smooth, valid confidence intervals/bans can be developed based on the distance-based estimator. The same inference results outlined for the location-based approach are available for the distance based approach, with some important caveats:
\begin{itemize}
\item If $\mathcal{B}$ is smooth, then RBC inference is possible. Thus, first the (I)MSE-rate-optimal bandwidth $h \asymp n^{-1/(2p+4)}$ is used for point estimation (i.e., $\widehat{\tau}_{\operatorname{dis}}(\mathbf{x})$) using $p$th order local polynomial regression, and then inference proceeds using $q$th order local polynomial regression. This is implemented using the option \texttt{kink = "off"}, and is the default for \texttt{rd2d.dist()}.
\item If $\mathcal{B}$ is non-smooth, then the RBC inference is not possible because the leading bias is unknown and increasing the polynomial order does not reduce bias. In this case, point estimation employs the (I)MSE-rate-optimal bandwidth $h \asymp n^{-1/4}$, and then inference employs the undersmoothed bandwidth choice $h \asymp n^{-1/3}$ following the results in \cite{Calonico-Cattaneo-Farrell_2018_JASA,Calonico-Cattaneo-Farrell_2022_Bernoulli}. As a result, point estimation and inference employ the same polynomial order ($p=q$).
\end{itemize}
\end{enumerate}
Other implementation and regularization methods follow the same logic as for bivariate location-based estimation, taking into account the distance variable explicitly. In particular, $\text{\texttt{bwcheck}} = 50 + p + 1$ is used as default. See \cite{Cattaneo-Titiunik-Yu_2025_BoundaryRD} for omitted technical and methodological details.
\section{Numerical Illustrations}\label{sec: Numerical Illustration}
We illustrate the capabilities of the general-purpose \texttt{R} software package \texttt{rd2d} with a synthetic dataset of size $n = 20,000$ calibrated using the \textit{Ser Pilo Paga} (SPP) dataset \citep{LondonoVelezRodriguezSanchez_2020_AEJ}. We set $T_i = \mathds{1}(X_{1i} \geq 0, X_{2i} \geq 0)$, where $T_i = 1$ indicates unit $i$ is in the treatment group, and $T_i = 0$ indicates unit $i$ is in the control group. Covariates $\mathbf{X}_i = (X_{1i}, X_{2i})^\top$ are drawn from the product distribution $(100\mathsf{Beta}(3,4)-25,\mathsf{B}_2 = 100\mathsf{Beta}(3,4)-25)$ with independent components. Potential outcomes are generated by
\begin{align*}
Y_{i}(t) = \beta_{t,0} + X_{1i} \beta_{t,11} + X_{2i} \beta_{t,12}
+ X_{1i}^2 \beta_{t,21} + X_{2i}^2 \beta_{t,22} + X_{1i}X_{2i} \beta_{t,23}
+ \varepsilon_{t,i},
\end{align*}
where $\mathbf{X}_i$, $\varepsilon_{0,i}$ and $\varepsilon_{1,i}$ are mutually independent, and $\varepsilon_{t,i}\thicksim \mathsf{N}(0,\sigma_t^2)$, for $i = 1,2, \cdots, n$ and $t = 0,1$. We consider two DGPs as in Table~\ref{tab:true-params}, where coefficients are estimated from the \textit{Ser Pilo Paga} (SPP) dataset \citep{LondonoVelezRodriguezSanchez_2020_AEJ}, and scaled by a factor of 2 to enhance signal strength.
\begin{table}[ht]
\centering
\begin{tabular}{@{}rrrrr@{}} \hline\hline
& \multicolumn{2}{c}{DGP 1 (Linear)} & \multicolumn{2}{c}{DGP 2 (Quadratic)} \\ \cmidrule(lr){2-3} \cmidrule(lr){4-5}
& \multicolumn{1}{c}{$t=0$} & \multicolumn{1}{c}{$t=1$} & \multicolumn{1}{c}{$t=0$} & \multicolumn{1}{c}{$t=1$} \\ \midrule
$\beta_{t,0}$ & $\mathbf{2} \times 3.35 \times 10^{-1}$ & $\mathbf{2} \times 6.98 \times 10^{-1}$ & $\mathbf{2} \times 3.72 \times 10^{-1}$ & $\mathbf{2} \times 7.435 \times 10^{-1}$ \\
$\beta_{t,11}$ & $\mathbf{2} \times 2.52 \times 10^{-3}$ & $\mathbf{2} \times 2.74 \times 10^{-3}$ & $\mathbf{2} \times 4.23 \times 10^{-3}$ & $\mathbf{2} \times 2.29 \times 10^{-3}$ \\
$\beta_{t,12}$ & $- \mathbf{2} \times 1.72 \times 10^{-3}$ & $- \mathbf{2} \times 6.05 \times 10^{-4}$ & $- \mathbf{2} \times 2.45 \times 10^{-3}$ & $- \mathbf{2} \times 5.85 \times 10^{-3}$ \\
$\beta_{t,21}$ & $0$ & $0$ & $\mathbf{2} \times 1.25 \times 10^{-5}$ & $- \mathbf{2} \times 1.33 \times 10^{-7}$ \\
$\beta_{t,22}$ & $0$ & $0$ & $- \mathbf{2} \times 4.92 \times 10^{-6}$ & $\mathbf{2} \times 2.14 \times 10^{-5}$ \\
$\beta_{t,23}$ & $0$ & $0$ & $\mathbf{2} \times 3.12 \times 10^{-5}$ & $\mathbf{2} \times 1.04 \times 10^{-4}$ \\
$\sigma_{t}$ & $3.32 \times 10^{-1}$ & $4.35 \times 10^{-1}$ & $3.31 \times 10^{-1}$ & $4.35 \times 10^{-1}$ \\
\hline\hline
\end{tabular}
\caption{True parameter values for DGP 1 (Linear) and DGP 2 (Quadratic) under treatment statuses $t=0$ and $t=1$, fitted from the \textit{Ser Pilo Paga} (SPP) dataset \citep{LondonoVelezRodriguezSanchez_2020_AEJ} and coefficients boosted by $2$.}
\label{tab:true-params}
\end{table}
Figure~\ref{fig:fig2a} presents a scatterplot of a synthetic dataset constructed using the data generating process described above, which we use for numerical illustration of the main capabilities of the package \texttt{rd2d}. The gray assignment boundary separates the control (blue) and treatment (red) groups. The plot also includes forty grid points along the boundary, where the cutoff $\mathbf{b}_{21}$ is a kink point. Figure~\ref{fig:fig2b} presents the population treatment effects corresponding to the two data generating processes. To demonstrate the variability in the outcome variable, the figure includes synthetic data points derived from 300 uniform draws along the boundary.
\begin{figure}
\centering
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{inputs/fig2a.png}
\caption{Scatter Plot of $(X_{1i}, X_{2i})$}
\label{fig:fig2a}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{inputs/fig2b-cond-mean.png}
\caption{Treatment Effects Along the Boundary}
\label{fig:fig2b}
\end{subfigure}
\caption{Estimation and Inference (Simulations).}
\label{fig:fig2}
\end{figure}
\subsection[Function rd2d()]{Function \texttt{rd2d()}}
The function \texttt{rd2d()} provides point estimation, robust confidence intervals, and robust uniform confidence bands for (derivatives of) treatment effect function $\tau(\mathbf{x})$ based on bivariate location-based local polynomial regression. It takes as input an outcome vector \texttt{y}, a bivariate location score matrix \texttt{X}, a treatment indicator vector \texttt{t}, and a grid of evaluation points \texttt{b} along the boundary $\mathcal{B}$.
Optional arguments include bandwidth choices \texttt{h}, degrees of polynomial for point estimation (\texttt{p}) and inference (\texttt{q}), the partial derivative of treatment effect to be estimated \texttt{deriv}, and confidence level \texttt{level}. When optional arguments are not provided, the function defaults to estimate the value of treatment effect \texttt{deriv = c(0,0)} using the MSE-optimal bandwidth, with \texttt{p = 1} degree polynomial for point estimation and \texttt{q = 2} degree polynomial for robust bias-corrected confidence interval (and bands if requested), and \texttt{level = 95} percentage points confidence level. Additionally, \texttt{kernel\_type} indicates whether a product kernel (\texttt{"prod"}) or a radial kernel (\texttt{"rad"}) is used for weighting. The default is \texttt{kernel\_type = "prod"}.
Below is a demonstration of \texttt{rd2d()} applied on the synthetic dataset, with results for selected indices printed using the \texttt{summary()} method.
\begin{lstlisting}[style=Rstyle, caption={Functionality of \texttt{rd2d()}}]
> result.rd2d <- rd2d(y, X, t, eval)
> summary(result.rd2d, subset = c(1,5,10,15,21,25,30,35,40))
rd2d
Number of Obs. 20000
BW type. mserd-dpi-std
Kernel triangular-prod
VCE method hc1
Masspoints check
Number of Obs. 6191 13809
Estimand (deriv) 0 0
Order est. (p) 1 1
Order rbc. (q) 2 2
Unique Obs. 6191 13809
====================================================================
ID b1 b2 Est. z P>|z| 95
====================================================================
1 0.000 50.000 0.8143 9.8761 0.0000 [0.6191, 0.9257]
5 0.000 40.000 0.7298 11.8383 0.0000 [0.6125, 0.8556]
10 0.000 27.500 0.6752 11.3927 0.0000 [0.5487, 0.7767]
15 0.000 15.000 0.6456 13.4458 0.0000 [0.5448, 0.7307]
21 0.000 0.000 0.6890 8.8464 0.0000 [0.5294, 0.8307]
25 10.000 0.000 0.6647 14.1815 0.0000 [0.5492, 0.7254]
30 22.500 0.000 0.6116 11.5752 0.0000 [0.5049, 0.7107]
35 35.000 0.000 0.5722 8.6064 0.0000 [0.4110, 0.6534]
40 47.500 0.000 0.5120 6.2842 0.0000 [0.3354, 0.6393]
====================================================================
\end{lstlisting}
The first part of the output provides basic information on the options specified in the function. For example, the default estimand is the value of treatment effect, indicated by \texttt{deriv = (0,0)}. The rest of the output gives estimation results, including (i) \texttt{b1} and \texttt{b2}: First and second coordinate of the evaluation points; (ii) \texttt{Coef.}: Point estimation of (derivative) of treatment effect using $p=$ \texttt{p} polynomial order; and (iii) t-statistics, (iv) p-value, and (v) \texttt{level\%} confidence intervals using $q=$ \texttt{q} polynomial order. When $q>p$, the resulting inference procedures correspond to robust bias correction \citep{Calonico-Cattaneo-Titiunik_2014_ECMA,Calonico-Cattaneo-Farrell_2018_JASA,Calonico-Cattaneo-Farrell_2022_Bernoulli}, which is the default and recommended method; $p=q$ corresponds to standard least squares methods. Point estimates, standard errors, and other information can be easily extracted for further statistical analysis. The output is stored in a standard matrix, and can be accessed with the following command,
\begin{lstlisting}[style=Rstyle, caption={Output matrices of \texttt{rd2d()}}]
> result.rd2d$results
> result.rd2d$results.A0
> result.rd2d$results.A1
\end{lstlisting}
where \texttt{result.rd2d\$main.A0} contains results for the control group, \texttt{result.rd2d\$main.A1} contains results for the treatment group, and \texttt{result.rd2d\$main} contains results for both.
The \texttt{summary()} method allows for three optional arguments. First, the option \texttt{subset} takes the indices of evaluation points to be presented, which should be a subset of \texttt{c(1:nrow(eval))}. The default is \texttt{NULL}, and thus all evaluation points are presented. Second, \texttt{CBuniform} is boolean variable for confidence bands construction, where \texttt{FALSE} indicates that pointwise confidence intervals are provided, and \texttt{TRUE} indicates uniform confidence bands are provided. The default is \texttt{CBuniform = FALSE}.
\begin{lstlisting}[style=Rstyle, caption={Functionality of \texttt{rd2d()} with \texttt{CBuniform = TRUE}}]
> summary(result.rd2d, subset = c(1,5,10,15,21,25,30,35,40), CBuniform = TRUE)
rd2d
Number of Obs. 20000
BW type. mserd-dpi-std
Kernel triangular-prod
VCE method hc1
Masspoints check
Number of Obs. 6191 13809
Estimand (deriv) 0 0
Order est. (p) 1 1
Order rbc. (q) 2 2
Unique Obs. 6191 13809
====================================================================
ID b1 b2 Est. z P>|z| 95
====================================================================
1 0.000 50.000 0.8143 9.8761 0.0000 [0.5322, 1.0126]
5 0.000 40.000 0.7298 11.8383 0.0000 [0.5436, 0.9245]
10 0.000 27.500 0.6752 11.3927 0.0000 [0.4841, 0.8414]
15 0.000 15.000 0.6456 13.4458 0.0000 [0.4920, 0.7834]
21 0.000 0.000 0.6890 8.8464 0.0000 [0.4439, 0.9161]
25 10.000 0.000 0.6647 14.1815 0.0000 [0.4993, 0.7753]
30 22.500 0.000 0.6116 11.5752 0.0000 [0.4465, 0.7690]
35 35.000 0.000 0.5722 8.6064 0.0000 [0.3423, 0.7221]
40 47.500 0.000 0.5120 6.2842 0.0000 [0.2492, 0.7255]
====================================================================
\end{lstlisting}
When the optional argument \texttt{AATE} $=\bigl(w_j\bigr)_{j=1}^J$, representing the weights for the evaluation points $\{\boldsymbol{b}_j\}_{j=1}^J$, is supplied to \texttt{summary()}, the output gains an extra row labeled “AATE.”
This row reports the aggregated average treatment effect estimator
\[
\widehat{\tau}_{\text{AATE}}
=\frac{\sum_{j=1}^{J} w_j\,\widehat{\tau}\!\bigl(\boldsymbol{b}_j\bigr)}
{\sum_{j=1}^{J} w_j},
\]
along with its inference statistics.
\begin{lstlisting}[style=Rstyle, caption={Functionality of \texttt{rd2d()} for AATE}]
> summary(result.rd2d, subset = c(1,5,10,15,21,25,30,35,40), AATE = rep(1,neval))
rd2d
Number of Obs. 20000
BW type. mserd-dpi-std
Kernel triangular-prod
VCE method hc1
Masspoints check
Number of Obs. 6191 13809
Estimand (deriv) 0 0
Order est. (p) 1 1
Order rbc. (q) 2 2
Unique Obs. 6191 13809
====================================================================
ID b1 b2 Est. z P>|z| 95
====================================================================
1 0.000 50.000 0.8143 9.8761 0.0000 [0.6191, 0.9257]
5 0.000 40.000 0.7298 11.8383 0.0000 [0.6125, 0.8556]
10 0.000 27.500 0.6752 11.3927 0.0000 [0.5487, 0.7767]
15 0.000 15.000 0.6456 13.4458 0.0000 [0.5448, 0.7307]
21 0.000 0.000 0.6890 8.8464 0.0000 [0.5294, 0.8307]
25 10.000 0.000 0.6647 14.1815 0.0000 [0.5492, 0.7254]
30 22.500 0.000 0.6116 11.5752 0.0000 [0.5049, 0.7107]
35 35.000 0.000 0.5722 8.6064 0.0000 [0.4110, 0.6534]
40 47.500 0.000 0.5120 6.2842 0.0000 [0.3354, 0.6393]
--------------------------------------------------------------------
AATE 0.6534 28.9598 0.0000 [0.5987, 0.6856]
====================================================================
\end{lstlisting}
Finally, the \texttt{summary()} method allows for presenting the underlying bandwidths used and associated effective sample sizes via the option \texttt{output = "bw"}.
\begin{lstlisting}[style=Rstyle, caption={Functionality of \texttt{rd2d()} with \texttt{output = "bw"}}]
> summary(result.rd2d, subset = c(1,5,10,15,21,25,30,35,40), output = "bw")
rd2d
Number of Obs. 20000
BW type. mserd-dpi-std
Kernel triangular-prod
VCE method hc1
Masspoints check
Number of Obs. 6191 13809
Estimand (deriv) 0 0
Order est. (p) 1 1
Order rbc. (q) 2 2
Unique Obs. 6191 13809
====================================================================
Bdy Points BW Control BW Treatment Eff. N
ID b1 b2 h01 h02 h11 h12 Nh0 Nh1
====================================================================
1 0.000 50.000 16.687 16.717 16.687 16.717 689 1295
5 0.000 40.000 14.782 14.808 14.782 14.808 1019 1810
10 0.000 27.500 12.736 12.759 12.736 12.759 1252 2108
15 0.000 15.000 12.045 12.067 12.045 12.067 1295 2167
21 0.000 0.000 14.740 14.766 14.740 14.766 2169 1552
25 10.000 0.000 12.455 12.477 12.455 12.477 1467 1980
30 22.500 0.000 11.476 11.496 11.476 11.496 1107 1843
35 35.000 0.000 11.595 11.615 11.595 11.615 824 1418
40 47.500 0.000 14.961 14.987 14.961 14.987 638 1242
====================================================================
\end{lstlisting}
\subsection[Function rdbw2d()]{Function \texttt{rdbw2d()}}
The function \texttt{rdbw2d()} is used for MSE (or IMSE) optimal bandwidth implementation, and is used internally in \texttt{rd2d()} when user does not specify bandwidth choices manually. The function takes the same input data as \texttt{rd2d()}, that is, an outcome vector \texttt{y}, a bivariate location score matrix \texttt{X}, a treatment indicator vector \texttt{t}, and a grid of evaluation points \texttt{b} along the boundary $\mathcal{B}$. In addition, the option \texttt{bwselect} encodes four options of bandwidth type: (i) \texttt{"mserd"} finds the MSE-optimal bandwidth for estimating (derivatives of) treatment effect, (ii) \texttt{"imserd"} finds the integreated MSE optimal bandwidth for estimating (derivatives of) treatment effect, (iii) \texttt{"msetwo"} finds the MSE optimal bandwidth for estimation (the derivatives of) conditional means of two potential outcome variables, (iv) \texttt{"imsetwo"} finds the integrated MSE optimal bandwidth for estimation (the derivatives of) conditional means of two potential outcome variables. The default is \texttt{bwselect = "mserd"}. An additional Boolean argument \texttt{stdvar} indicates whether the covariates are first standardized to unit standard deviation in each coordinate, in which case the optimal bandwidth is estimated and then converted back to the original scale. The default is \texttt{stdvar = TRUE}.
\begin{lstlisting}[style=Rstyle, caption={Functionality of \texttt{rdbw2d()}}]
> bws.rd2d <- rdbw2d(y, X, t, eval)
> summary(bws.rd2d,subset = c(1,5,10,15,21,25,30,35,40))
Call: rdbw2d
Number of Obs. 20000
BW type. mserd-dpi
Kernel triangular-prod
VCE method hc1
Masspoints check
Standardization on
Number of Obs. 6191 13809
Estimand (deriv) 0 0
Order est. (p) 1 1
Unique Obs. 6191 13809
Bandwidth Selection
====================================================
Bdy Points BW Control BW Treatment
ID b1 b2 h01 h02 h11 h12
====================================================
1 0.000 2.855 20.552 20.588 20.552 20.588
5 0.000 2.284 18.668 18.701 18.668 18.701
10 0.000 1.570 16.476 16.505 16.476 16.505
15 0.000 0.856 14.098 14.123 14.098 14.123
21 0.000 0.000 19.115 19.149 19.115 19.149
25 0.572 0.000 15.737 15.765 15.737 15.765
30 1.287 0.000 14.000 14.025 14.000 14.025
35 2.002 0.000 13.562 13.586 13.562 13.586
40 2.717 0.000 18.798 18.832 18.798 18.832
====================================================
\end{lstlisting}
The first part of the \texttt{summary} output lists the options used for bandwidth selection. The second part of the \texttt{summary} output gives bandwidth selection results, including: (i) Boundary points, \texttt{b1} for the first coordinate and \texttt{b2} for the second coordinate; (ii) Bandwidths for control group, \texttt{h01} for the first coordinate and \texttt{h02} for the second coordinate; (iii) Bandwidths for treatment group, \texttt{h11} for the first coordinate and \texttt{h12} for the second coordinate.
\subsection[Function rd2d.dist()]{Function \texttt{rd2d.dist()}}
The function \texttt{rd2d.dist()} provides point estimation and inference for boundary treatment effects using distance-based univariate local polynomial regression. It takes as input an outcome vector \texttt{y}, and a signed distance matrix \texttt{D} of distance to each boundary point, where each column of \texttt{D} corresponds to the signed distance from all observations to one evaluation point, with a positive sign indicating the unit is in the treatment group and a negative sign indicating the unit is in the control group.
Optional arguments include evaluation points \texttt{b}, bandwidth choices \texttt{h}, degrees of polynomial for point estimation (\texttt{p}) and inference (\texttt{q}), option for kink adjustment \texttt{kink}, and confidence level \texttt{level}, among other options. When not provided, the function defaults to estimate the value of treatment effect using the MSE-optimal bandwidth without kink adjustment (\texttt{kink = "off"}), using \texttt{p = 1} degree polynomial for point estimation and \texttt{q = 2} degree polynomial for robust bias-corrected confidence intervals and bands, providing \texttt{level = 95}\% confidence interval and uniform confidence bands, and without displaying of evaluation points (\texttt{b = NULL}).
\begin{lstlisting}[style=Rstyle, caption={Functionality of \texttt{rd2d.dist()}}]
> result.dist <- rd2d.dist(y,D, b = eval)
> summary(result.dist, subset = c(1,5,10,15,21,25,30,35,40))
rd2d.dist
Number of Obs. 20000
BW type mserd-rot
Kernel triangular-rad
VCE method hc1
Masspoints check
Number of Obs. 6191 13809
Estimand (deriv) 0 0
Order est. (p) 1 1
Order rbc. (q) 2 2
====================================================================
ID b1 b2 Est. z P > |z| 95
====================================================================
1 0.000 50.000 0.8371 11.5847 0.0000 [0.6514, 0.9167]
5 0.000 40.000 0.7307 9.4650 0.0000 [0.5807, 0.8840]
10 0.000 27.500 0.6572 9.7763 0.0000 [0.5123, 0.7692]
15 0.000 15.000 0.6272 7.7001 0.0000 [0.4344, 0.7311]
21 0.000 0.000 0.7426 6.8180 0.0000 [0.5056, 0.9136]
25 10.000 0.000 0.6301 11.2414 0.0000 [0.5474, 0.7785]
30 22.500 0.000 0.6109 7.3862 0.0000 [0.4337, 0.7470]
35 35.000 0.000 0.5768 6.4730 0.0000 [0.3668, 0.6855]
40 47.500 0.000 0.4654 6.7582 0.0000 [0.3602, 0.6544]
====================================================================
\end{lstlisting}
The first part of the \texttt{summary} output provides basic information on the options specified in the function. The rest of the \texttt{summary} output gives estimation and inference results, including: (i) \texttt{b1} and \texttt{b2} (when \texttt{b} is provided) report first and second coordinates of the evaluation points; (ii) \texttt{Coef.} reports the treatment effect estimate using a $p$th order polynomial regression; and (iii) the last three columns correspond to t-statistic, p-value and \texttt{level}\% confidence intervals using $q$th order polynomial regression. The \texttt{summary()} method also has the option of displaying the uniform confidence bands (instead of the confidence intervals) in the last two columns as follows (numerical results omitted to conserve space).
\begin{lstlisting}[style=Rstyle, caption={Functionality of \texttt{rd2d.dist()} with \texttt{CBuniform = TRUE}}]
> summary(result.dist, subset = c(1,5,10,15,21,25,30,35,40), CBuniform = TRUE)
\end{lstlisting}
When the optional argument \texttt{AATE} $=\bigl(w_j\bigr)_{j=1}^J$, representing the weights for the evaluation points $\{\boldsymbol{b}_j\}_{j=1}^J$, is supplied to \texttt{summary()}, the output gains an extra row labeled “AATE” for estimation and inference of aggregated average treatment effect.
\begin{lstlisting}[style=Rstyle, caption={Functionality of \texttt{rd2d()} for AATE}]
> summary(result.dist, subset = c(1,5,10,15,21,25,30,35,40), AATE = rep(1,neval))
\end{lstlisting}
In addition, The \texttt{summary()} method can also display the underlying bandwidths and effective samples sizes as follows (numerical results omitted to conserve space).
\begin{lstlisting}[style=Rstyle, caption={Functionality of \texttt{rd2d.dist()} with \texttt{output = "bw"}}]
> summary(result.dist, subset = c(1,5,10,15,21,25,30,35,40), output = "bw")
\end{lstlisting}
Point estimates, standard errors, and other information can be easily extracted for further statistical analysis. The output is stored in a standard matrix, and can be accessed with the following command.
\begin{lstlisting}[style=Rstyle, caption={Output matrices of \texttt{rd2d.dist()}}]
> result.dist$results
> result.dist$results.A0
> result.dist$results.A1
\end{lstlisting}
\texttt{result.dist\$main.A0} contains results for the control group, \texttt{result.dist\$main.A1} contains results for the treatment group, and \texttt{result.dist\$main} contains results for both.
\subsection[Function rdbw2d.dist()]{Function \texttt{rdbw2d.dist()}}
The function \texttt{rdbw2d.dist()} is used to implement MSE (or IMSE) rate-optimal ROT bandwidth selectors, and is used internally in \texttt{rd2d.dist()} when user does not provide bandwidths manually. As for \texttt{rdbw2d()}, four bandwidth types are allowed, that is, \texttt{bwselect} can be \texttt{"mserd"}, \texttt{"imserd"}, \texttt{"msetwo"} or \texttt{"imsetwo"}.
\begin{lstlisting}[style=Rstyle, caption={Functionality of \texttt{rdbw2d.dist()}}]
> bws.dist <- rdbw2d.dist(y,D, b = eval)
> summary(bws.dist,subset = c(1,5,10,15,21,25,30,35,40))
Call: rdbw2d.dist
Number of Obs. 20000
BW type mserd-rot
Kernel triangular-rad
Kink off
VCE method hc1
Masspoints check
Number of Obs. 6191 13809
Estimand (deriv) 0 0
Order est. (p) 1 1
Bandwidth Selection
================================================
Bdy Points BW Control BW Treatment
ID b1 b2 h0 h1
================================================
1 0.000 50.000 32.918 32.918
5 0.000 40.000 22.071 22.071
10 0.000 27.500 20.899 20.899
15 0.000 15.000 14.455 14.455
21 0.000 0.000 14.381 14.381
25 10.000 0.000 19.712 19.712
30 22.500 0.000 14.041 14.041
35 35.000 0.000 16.722 16.722
40 47.500 0.000 27.499 27.499
================================================
\end{lstlisting}
An additional argument \texttt{kink}, taking values \texttt{"off"} or \texttt{"on"}, indicates whether a kink adjustment is made for estimation and inference. The default is \texttt{kink = "off"}, but when \texttt{kink = "on"} is specified then the bandwidth is shrank to account for lack of smoothness of the assignment boundary $\mathcal{B}$.
\begin{lstlisting}[style=Rstyle, caption={Functionality of \texttt{rdbw2d.dist()} with \texttt{kink = "on"}}]
> bws.dist <- rdbw2d.dist(y,D, kink = "on")
> summary(bws.dist,subset = c(1,5,10,15,21,25,30,35,40))
Call: rdbw2d.dist
Number of Obs. 20000
BW type mserd-rot
Kernel triangular-rad
Kink on
VCE method hc1
Masspoints check
Number of Obs. 6191 13809
Estimand (deriv) 0 0
Order est. (p) 1 1
Bandwidth Selection
================================
BW Control BW Treatment
ID h0 h1
================================
1 14.421 14.421
5 9.670 9.670
10 9.156 9.156
15 6.333 6.333
21 6.300 6.300
25 8.636 8.636
30 6.152 6.152
35 7.326 7.326
40 12.048 12.048
================================
\end{lstlisting}
\subsection{Graphical Presentation}
The package \texttt{rd2d} provides an array of estimation and inference results that can be used for graphical presentation. Figure~\ref{fig:fig3a} compares three point estimation methods: (i) bivariate method via \texttt{rd2d}; (ii) distance-based method via \texttt{rd2d.dist}, ignoring the kink (default argument \texttt{kink = "off"}); and (iii) distance-based method via \texttt{rd2d.dist}, adjusting for kink (\texttt{kink = "on"}). Figure~\ref{fig:fig3b} plots point estimation using the bivariate method, along with its associated robust bias-corrected confidence intervals and confidence band. Figure~\ref{fig:fig3c} presents a heatmap of treatment effects along the boundary, with high to low point estimation indicated by red to blue colors. Finally, Figure~\ref{fig:fig3d} presents a heatmap of p-values along the boundary, with five colors assigned to five ranges of values. The codes for generating the graphical presentations are given in the replication \texttt{R} file.
It can be seen that distance-based estimation with \texttt{kink = "off"} overshoots compared to the bivariate estimation before the kink, and undershoots after the kink. This corresponds to the phenomena of getting a first-order bias using distance based method in the presence of a kink, despite using local polynomial regression of degree greater than or equal to $1$. See \cite{Cattaneo-Titiunik-Yu_2025_BoundaryRD} for more methodological and theoretical discussions.
\begin{figure}
\centering
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{inputs/fig3a-point-estimation-simulation.png}
\caption{Point Estimation}
\label{fig:fig3a}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{inputs/fig3b-ci-and-cb-simulation.png}
\caption{Confidence Interval and Bands}
\label{fig:fig3b}
\end{subfigure}
\vspace{0.5cm}
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{inputs/fig3c-heat-spp-simulation.png}
\caption{Treatment Effects Heatmap}
\label{fig:fig3c}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\linewidth]{inputs/fig3d-heat-pvalue-simulation.png}
\caption{p-values Heatmap}
\label{fig:fig3d}
\end{subfigure}
\caption{Estimation and Inference (Simulations).}
\label{fig:fig3}
\end{figure}
\subsection{Simulation Evidence}
The discussion so far employed one realization of the data generating process to illustrate the main features of the package \texttt{rd2d}. In this final section, we conduct a Monte Carlo experiment to assess the performance of the package in repeated sampling. We consider $m = 1,000$ simulations of the two data generating processes defined at the beginning of this section, and report the simulation results in Table~\ref{tab:sim-linear} (DGP 1: linear model) and Table~\ref{tab:sim-quadratic} (DGP 2: quadratic model). Three methods are used and compared: the bivariate method \texttt{rd2d}, the distance-based method \texttt{rd2d.dist} ignoring the presence of the kink in the boundary (\texttt{kink = "off"}), and the distance-based method \texttt{rd2d.dist} adjusting for the kink (\texttt{kink = "on"}). Bandwidths are chosen automatically by the package, and their average across simulations is reported. We also report diagnostic measures including: bias, standard deviation of point estimator, root mean-squared error of point estimator, pointwise empirical coverage, pointwise interval length, uniform empirical coverage, and uniform interval length.
For both DPG 1 and DGP 2, all of the three methods give pointwise coverage around 95\%, while the pointwise interval length for \texttt{rd2d} and \texttt{rd2d.dist (kink = "off")} are smaller compared to \texttt{rd2d.dist (kink = "on")}. This is likely due to bandwidth shrinkage for kink adjustment, which does results in a significantly smaller bias compared to the one ignoring the kink. Both \texttt{rd2d} and \texttt{rd2d.dist (kink = "off")} give around 95\% uniform coverage and a relatively shorter interval length compared to \texttt{rd2d.dist (kink = "on")}, likely due to the same reason.
\begin{table}[ht]
\centering
\begin{tabular}{@{}lrrrrrrr@{}} \toprule \toprule
Method & Index & $h$ & Bias & SD & RMSE & EC & IL\\ \midrule
\multirow{10}{*}{\texttt{rd2d}}
& 1 & 15.764 & 0.005 & 0.053 & 0.054 & 0.953 & 0.313 \\
& 5 & 13.885 & 0.003 & 0.042 & 0.042 & 0.958 & 0.243 \\
& 10 & 11.997 & 0.002 & 0.038 & 0.038 & 0.952 & 0.218 \\
& 15 & 12.706 & 0.001 & 0.033 & 0.033 & 0.949 & 0.191 \\
& 21 & 13.788 & -0.003 & 0.049 & 0.049 & 0.957 & 0.329 \\
& 25 & 12.908 & 0.001 & 0.031 & 0.031 & 0.950 & 0.177 \\
& 30 & 11.688 & -0.001 & 0.036 & 0.036 & 0.948 & 0.212 \\
& 35 & 13.012 & 0.001 & 0.039 & 0.039 & 0.957 & 0.228 \\
& 40 & 15.255 & 0.001 & 0.049 & 0.049 & 0.960 & 0.289 \\ \cline{2-8}
& Uniform & & & & & 0.948 & 0.348 \\
\midrule
\multirow{10}{*}{\shortstack{\texttt{rd2d.dist}\\\texttt{kink = "off"}}}
& 1 & 34.136 & 0.031 & 0.038 & 0.049 & 0.940 & 0.261 \\
& 5 & 26.747 & 0.016 & 0.035 & 0.038 & 0.950 & 0.240 \\
& 10 & 19.478 & 0.003 & 0.035 & 0.035 & 0.953 & 0.258 \\
& 15 & 16.905 & 0.001 & 0.038 & 0.038 & 0.957 & 0.275 \\
& 21 & 21.151 & 0.002 & 0.038 & 0.038 & 0.960 & 0.283 \\
& 25 & 18.524 & -0.021 & 0.039 & 0.044 & 0.955 & 0.255 \\
& 30 & 16.986 & -0.000 & 0.037 & 0.037 & 0.948 & 0.280 \\
& 35 & 24.532 & -0.005 & 0.033 & 0.033 & 0.956 & 0.233 \\
& 40 & 33.070 & -0.015 & 0.035 & 0.038 & 0.953 & 0.244 \\ \cline{2-8}
& Uniform & & & & & 0.946 & 0.408 \\
\midrule
\multirow{10}{*}{\shortstack{\texttt{rd2d.dist}\\\texttt{kink = "on"}}}
& 1 & 14.955 & 0.012 & 0.079 & 0.080 & 0.954 & 0.729 \\
& 5 & 11.718 & 0.003 & 0.072 & 0.072 & 0.948 & 0.646 \\
& 10 & 8.533 & 0.001 & 0.078 & 0.078 & 0.954 & 0.698 \\
& 15 & 7.406 & -0.002 & 0.083 & 0.083 & 0.948 & 0.753 \\
& 21 & 9.266 & -0.003 & 0.083 & 0.083 & 0.944 & 0.810 \\
& 25 & 8.116 & 0.005 & 0.076 & 0.077 & 0.955 & 0.697 \\
& 30 & 7.442 & -0.001 & 0.079 & 0.079 & 0.951 & 0.761 \\
& 35 & 10.748 & 0.002 & 0.067 & 0.067 & 0.948 & 0.626 \\
& 40 & 14.488 & -0.001 & 0.075 & 0.075 & 0.948 & 0.667 \\ \cline{2-8}
& Uniform & & & & & 0.89 & 1.073 \\
\bottomrule
\end{tabular}
\caption{Simulation results for DGP 1 (Linear)}
\label{tab:sim-linear}
\end{table}
\begin{table}[ht]
\centering
\begin{tabular}{@{}lrrrrrrr@{}} \toprule \toprule
Method & Index & $h$ & Bias & SD & RMSE & EC & IL \\ \midrule
\multirow{10}{*}{\texttt{rd2d}}
& 1 & 15.730 & 0.004 & 0.052 & 0.053 & 0.958 & 0.315 \\
& 5 & 13.803 & 0.005 & 0.041 & 0.041 & 0.965 & 0.244 \\
& 10 & 11.954 & 0.003 & 0.038 & 0.038 & 0.943 & 0.218 \\
& 15 & 12.600 & 0.001 & 0.033 & 0.033 & 0.944 & 0.193 \\
& 21 & 13.716 & -0.008 & 0.049 & 0.049 & 0.962 & 0.331 \\
& 25 & 12.833 & -0.003 & 0.031 & 0.031 & 0.952 & 0.177 \\
& 30 & 11.679 & 0.000 & 0.035 & 0.034 & 0.962 & 0.212 \\
& 35 & 12.966 & -0.001 & 0.041 & 0.041 & 0.954 & 0.230 \\
& 40 & 15.312 & -0.007 & 0.050 & 0.050 & 0.942 & 0.288 \\ \cline{2-8}
& Uniform & & & & & 0.942 & 0.350 \\
\midrule
\multirow{10}{*}{\shortstack{\texttt{rd2d.dist}\\\texttt{kink = "off"}}}
& 1 & 34.303 & 0.040 & 0.038 & 0.055 & 0.927 & 0.260 \\
& 5 & 24.471 & 0.016 & 0.037 & 0.041 & 0.960 & 0.262 \\
& 10 & 18.392 & -0.001 & 0.036 & 0.036 & 0.951 & 0.273 \\
& 15 & 14.179 & -0.003 & 0.042 & 0.042 & 0.944 & 0.324 \\
& 21 & 21.170 & -0.004 & 0.037 & 0.038 & 0.959 & 0.283 \\
& 25 & 16.036 & -0.020 & 0.045 & 0.049 & 0.943 & 0.293 \\
& 30 & 13.771 & -0.002 & 0.043 & 0.043 & 0.947 & 0.341 \\
& 35 & 19.017 & -0.016 & 0.043 & 0.046 & 0.954 & 0.297 \\
& 40 & 31.541 & -0.071 & 0.042 & 0.083 & 0.935 & 0.260 \\ \cline{2-8}
& Uniform & & & & & 0.949 & 0.461 \\
\midrule
\multirow{10}{*}{\shortstack{\texttt{rd2d.dist}\\\texttt{kink = "on"}}}
& 1 & 15.028 & 0.015 & 0.079 & 0.080 & 0.950 & 0.718 \\
& 5 & 10.721 & 0.005 & 0.074 & 0.074 & 0.958 & 0.710 \\
& 10 & 8.058 & 0.000 & 0.082 & 0.082 & 0.937 & 0.746 \\
& 15 & 6.212 & -0.008 & 0.097 & 0.097 & 0.942 & 0.903 \\
& 21 & 9.275 & -0.001 & 0.084 & 0.084 & 0.945 & 0.811 \\
& 25 & 7.026 & 0.001 & 0.088 & 0.088 & 0.950 & 0.810 \\
& 30 & 6.033 & 0.006 & 0.102 & 0.102 & 0.944 & 0.947 \\
& 35 & 8.331 & -0.001 & 0.090 & 0.090 & 0.957 & 0.817 \\
& 40 & 13.818 & -0.020 & 0.079 & 0.081 & 0.952 & 0.708 \\\cline{2-8}
& Uniform & & & & & 0.822 & 1.199 \\
\bottomrule
\end{tabular}
\caption{Simulation results for DGP 2 (Quadratic)}
\label{tab:sim-quadratic}
\end{table}
\section{Conclusion}\label{sec: Conclusion}
This paper introduced the \texttt{R} software package \texttt{rd2d} for causal inference in Boundary Discontinuity designs. The package provides pointwise and uniform (over the treatment assignment boundary) estimation and inference methods employing either a bivariate location score or a univariate distance score. In addition, the methods can be used for graphical presentation. From a methodological perspective, this paper introduced second generation bandwidth selection methods complementing the main results in \citet{Cattaneo-Titiunik-Yu_2025_BoundaryRD}. Simulation evidence demonstrated a good performance of the package \texttt{rd2d}. Replication codes and related information are available at: \url{https://rdpackages.github.io/rd2d/}.
\section{Acknowledgments}
We thank Xinwei Ma and Filippo Palomba for insightful discussions. Cattaneo and Titiunik gratefully acknowledge financial support from the National Science Foundation (SES-2019432 and SES-2241575). Cattaneo gratefully acknowledge financial support from the Data-Driven Social Science initiative at Princeton University.
\clearpage
\bibliographystyle{plainnat}
\bibliography{CTY_2025_rd2d--bib}