Global Testing in Multivariate Regression Discontinuity Designs

Global Testing in Multivariate Regression Discontinuity Designs

Introduction

Regression discontinuity (RD) designs are widely regarded as one of the "...leading quasi-experimental empirical strategies in economics, political science, education, and many other social and behavioral sciences" (calonico2014robust). This paper develops a global testing framework for multivariate regression discontinuity designs, allowing researchers to determine whether discontinuities in conditional expectations or densities exist anywhere along a multivariate treatment boundary, even in settings with limited sample sizes. From a theoretical perspective, the paper establishes a global testing procedure whose validity does not rely heavily on precise estimation of local multivariate discontinuities. Practically, the proposed procedure exhibits good finite-sample size control even when the sample size is modest. \\

While traditional RD designs use a single running variable, recent research has increasingly explored settings with multiple running variables (papay2011extending, reardon2012regression, wong2013analyzing). For example, keele2015geographic examine geographic RD designs, where treatment and control groups are separated by geographic boundaries in two-dimensional space. matsudaira2008mandatory studies the effects of mandatory summer school, where treatment assignment depends on both math and reading scores. londono2020upstream explore a government tuition subsidy in Colombia where eligibility was determined by both merit (via a test score) and economic need (via an index of wealth). For more applications, see papers such as frey2019cash, narita2021algorithm, elacqua2016short, egger2015impact, evans2017smart, cohodes2014merit, becht2016does. \\

In response to these empirical applications, several methodological papers have been written to estimate treatment effects in multivariate RD settings. Notable contributions include imbens2019optimized, who use convex optimization, cheng2023estimation, who applies thin plate splines, and gunsilius2023free, who develop methods for cases where the treatment boundary is unknown. More recently, cattaneo2025estimationdist and cattaneo2025estimation explore using local polynomial regressions to estimate treatment effects along two-dimensional boundaries, as well as aggregated effects. \\

Despite these methodological developments, relatively little attention has been devoted to global testing: the task of determining whether a discontinuity exists at any point along a multivariate treatment boundary. A straightforward strategy is to estimate treatment effects at many boundary locations and construct uniform confidence intervals at each point, assessing whether zero lies outside any of the resulting bands. When the sample size is large, this approach is appealing: it not only detects whether a discontinuity is present but also reveals where along the boundary it occurs. As a result, uniform confidence bands provide a rich and granular view of underlying treatment effect patterns. However, their performance deteriorates when the number of running variables increases or when the sample size is modest. In such settings, the number of observations effectively informing each boundary point becomes small, leading to unreliable asymptotic approximations. Through a simple simulation in the next subsection, this standard approach is shown to exhibit poor finite-sample performance. \\

To combat this issue, this paper proposes a new approach to global testing where a multivariate estimator is used to pool individual treatment effects together, which allows for more robust global inference even with smaller sample sizes. In simulations, a test size closer to a target 5% is achieved even with a modest sample size. Even with two or more running variables, this approach can be used to test for discontinuities in multivariate conditional expectation functions and multivariate densities. This method is intended to complement existing multivariate RD estimators by providing a robust tool for assessing the reliability of estimated effects and confidence intervals. \\

The proposed global testing procedure works in two stages. First, a multivariate estimator is used to estimate whether the treatment effects along the boundary are positive or negative. The local linear forest estimator of friedberg2020local is used to estimate treatment effect heterogeneity, while a variation of the random forest density estimator of wen2022random is used for the density discontinuity test. Afterwards, these estimates are incorporated into a univariate test statistic that is used in the final hypothesis test. This approach has a number of advantages. First, by pooling together estimates along the boundary, the asymptotics for the final test statistic are made more reliable than they would be for each estimate individually. Second, the random forest estimators help maintain reasonable performance even when the number of running variables becomes relatively large. Third, the final test statistic only uses the signs of the multivariate estimates rather than the complete point estimates. This makes the final test statistic robust even when the initial treatment effects are estimated poorly, which translates to powerful theoretical properties. \\

This paper makes two contributions. First, it contributes to the literature on multivariate regression discontinuity designs by proposing a global test that remains robust even in relatively small samples. Second, it contributes to the literature on bunching and manipulation testing in regression discontinuity designs. Several papers have been written on univariate manipulation tests (cattaneo2020simple, mccrary2008manipulation, cattaneo2024local). crippa2025manipulation extends these types of tests to accommodate many running variables. However, this test only applies to boundaries like the one in Figure 1 in the next subsection. The test outlined in this paper can be applied to a wide range of boundary shapes, making it applicable even in settings with complex treatment rules, as in geographic regression discontinuities. \\

The remainder of this paper is organized as follows. Subsection 1.1 discusses a simulation example that demonstrates where the uniform confidence bands approach fails, while subsection 1.2 discusses notation and assumptions. Section 2 discusses identification for both the conditional expectation and density cases. Section 3 discusses estimation and inference for the final test statistics used for global testing. Section 4 shows simulation results that demonstrate the properties of the outlined global testing procedure. Section 5 discusses an application, and section 6 concludes.

Simple Simulation Example

To illustrate the limitations of the uniform confidence band approach in finite samples, consider the following process: $$ Y_i = \frac{X_{1,i} + X_{2,i}}{3} + \epsilon_i $$ where $\epsilon_i \sim \mathcal{N}(0, 0.05)$, $X_{1,i}, X_{2,i}$ come from independent uniform distributions on the interval $[-1, 1]$, and all observations are i.i.d across $i$. Further, define the treatment indicator as follows: $$ T_i =

$$ For this treatment rule, the treatment region looks as follows:

The goal is to test whether there is a discontinuity in $\mathbb{E}[Y_i | X_i = x]$ for some $x$ along the boundary between the treated and non-treated regions. To do this, the local polynomial estimator of cattaneo2025estimation is used to estimate treatment effect at 79 evenly-spaced points along the boundary and apply 95% uniform confidence bands for inference. Afterwards, rejection rates are computed over 5,000 simulations. Since there are no discontinuities in this example, the rejection rate should be approximately 5%. To assess robustness, the rejection rates from the sum of treatment effects and sum of squared treatment effects are also reported. Additionally, these three approaches are performed only on "safe" points $x = (0.5, 0)$ and $x = (0, 0.5)$ to avoid estimation at edges and kink points. Using number of observations $n = 20,000$ and $n = 1,000$, the rejection rates are as follows:

In the case where $n = 20,000$, the rejection rates are fairly close to 0.05. However, with a smaller number of observations, the rejection rate becomes much larger than the expected rejection rate of 0.05. Summation mitigates this problem only partially, because the aggregated statistic remains sensitive to inaccuracies in the individual local estimates used to construct it.

Setup and Notation

Let $Y_i$ denote the outcome variable for observation $i$ and let $X_i$ be a vector of running variables with support $\mathbb{X}$. Let $\Omega$ denote the set of treated units and let $\partial \Omega$ be the boundary of $\Omega$. Let $T_i = 1$ for $i \in \Omega$ and $T_i = 0$ for $i \notin \Omega$. Also, let $Y_i(t)$ denote the potential outcome for unit $i$ with treatment status $t$. Let $g_\Omega(X_i)$ be the signed (Euclidean) distance function relative to $\partial \Omega$, with $f_{g_\Omega}(g)$ being the density of $g_\Omega(X_i)$. For simplicity, define $g_\Omega(X_i) \equiv g_i$. With that, here are some necessary assumptions:

Identification

This section develops the identification framework for two proposed global test statistics, covering both conditional expectation discontinuities (treatment effect heterogeneity) and density discontinuities (manipulation).

Treatment Effect Heterogeneity

Define $\tau(x) = \mathbb{E}[Y_i(1) - Y_i(0) | X_i = x]$. Here, the hypothesis of interest is as follows:

To test the hypothesis above without uniform confidence bands, the running vector $X_i$ is first aggregated into a scalar $g_i$. Using $g_i$ as a running variable leads to the following:

Proofs of Theorem 1 and all subsequent theorems are deferred to the appendix. Intuitively, Theorem 1 implies that when the signed distance function is used as a running variable, the standard RD difference identifies a weighted average across $\tau(x)$ for all $x \in \partial \Omega$. Although this quantity may be of interest to researchers, it is not sufficient to test the main hypothesis because $\mathbb{E}[Y_i(1) - Y_i(0) | g_i = 0]$ may equal zero even when there exists points $x \in \partial \Omega$ where $\tau(x) \neq 0$. This can happen when the positive discontinuities mix with negative discontinuities to produce a null result. To address this issue, let $\gamma(X_i) \equiv \gamma_i$ denote the closest boundary point to $X_i$ and let $\Gamma(X_i) \equiv \Gamma_i = 1$ if $\tau(\gamma(X_i)) \geq 0$ and $\Gamma_i = 0$ otherwise. This leads to the following result:

Here, $\tau_\Gamma$ represents the weighed average of $|\tau(x)|$ over $x \in \partial \Omega$. Intuitively, $\Gamma_i$ serves as a sign‐normalization factor: it reverses the sign of $Y_i$ whenever the corresponding treatment effect is negative, which aligns all effects in the same direction. In other words, the original hypothesis can be rewritten as:

To make estimation easier and more efficient, notice that $\tau_\Gamma$ can be rewritten as: $$ \tau_\Gamma = \lim_{\epsilon \rightarrow 0^+} \mathbb{E}[(2\Gamma_i - 1)Y_i | g_i = \epsilon] - \lim_{\epsilon \rightarrow 0^-} \mathbb{E}[(2\Gamma_i - 1)Y_i | g_i = \epsilon] $$ A key advantage of this approach is that it only uses the sign of $\tau(\gamma(X_i))$. As a result, $\Gamma_i$ can be estimated reliably even if the underlying $\tau(\gamma(X_i))$ is noisy, which is an important property in small samples or in settings with many running variables. In contrast, procedures that aggregate raw treatment effects along the boundary are fragile: a single poorly estimated effect can materially distort the resulting test statistic. This phenomenon is evident in Table 1.

Density Manipulation

Consider the problem of testing whether the joint density of the running variables is continuous at the boundary. Such a test serves as a diagnostic for endogenous manipulation of the assignment variables. In this sense, the density-discontinuity test of mccrary2008manipulation is being extended to a multivariate setting. Define the following:

In this case, the hypothesis of interest is as follows:

Using a change-of-variables for densities: $$ f_{g_\Omega}(g) = \int_{x: g_\Omega(X) = g} f_X(x) dS_\Omega \implies \lim_{\epsilon \rightarrow 0^+} f_{g_\Omega}(\epsilon) - \lim_{\epsilon \rightarrow 0^-} f_{g_\Omega}(\epsilon) = \int_{x: g_\Omega(x) = 0} \tau_f(x)dS_\Omega $$ Let $\Lambda(X_i) \equiv \Lambda_i = 1$ if $\tau_f(\gamma(X_i)) \geq 0$ and $\Lambda_i = 0$ otherwise. With that, the result below follows:

Here, $\tau_\Lambda^f$ is equivalent to the sum of $|\tau_f(x)|$ over all $x \in \partial \Omega$. Like in the boundary heterogeneity case, a multivariate density test can be written through the following hypothesis:

To make this expression easier to estimate, $\tau_\Lambda^f$ can be rewritten as follows: $$ \tau_\Lambda^f = \lim_{\epsilon \rightarrow 0^+} f_{(2\Lambda-1)g}(\epsilon) - \lim_{\epsilon \rightarrow 0^-} f_{(2\Lambda-1)g}(\epsilon) $$

In other words, $\tau_\Lambda^f$ can be computed as the discontinuity in the density of $g_i^* \equiv (2\Lambda_i - 1)g_i$ at $g_i^* = 0$.

Estimation and Inference

The following subsections describe estimation and inference procedures for both the treatment effect heterogeneity test and the density manipulation test. Although the two procedures are closely related, there are several important distinctions between them.

Treatment Effect Heterogeneity

Based on the identification results of the previous section, the quantity of interest is the following: $$ \tau_\Gamma = \lim_{\epsilon \rightarrow 0^+} \mathbb{E}[Y_i^* | g_i = \epsilon] - \lim_{\epsilon \rightarrow 0^-} \mathbb{E}[Y_i^* | g_i = \epsilon] $$ where $Y_i^* = (2\Gamma(X_i) - 1)Y_i$. However, $\Gamma(X_i)$ is not observed, so it will need to be estimated as $\Gamma_n(X_i)$. Now define $\hat{Y}_i = (2\Gamma_n(X_i) - 1)Y_i$. With that, notice the following:

In other words, $\hat{\eta}_i$ can be thought of as a measurement error term, where $\Gamma_i$ is a higher dimensional nuisance parameter. As in chernozhukov2018double, estimating $\Gamma_i$ and $\tau_\Gamma$ using the same data will cause problems for estimation and inference. To address these issues, $\hat{Y}_i$ will be generated as follows:

where the index $k$ denotes the fold of observation $i$. To estimate $\Gamma(\cdot)$, the Local Linear Forest of friedberg2020local will be used. The advantage of this method is that it easily allows for efficient boundary point estimation with any number of running variables.

Once the outcome variable is defined, the fold-level version of $\tau_\Gamma$ (defined as $\tau_{\Gamma_{k,p}}$) will be estimated as follows:\footnote{The superscript $\pm$ is used as shorthand, with $+$ indexing quantities computed for observations with $g_{i,k} \ge 0$ and $-$ indexing quantities computed for observations with $g_{i,k} < 0$.}

where $K(\cdot)$ is some bounded, symmetric kernel function, $r_p(\cdot)$ is a vector of polynomial terms up to order $p$, and $h_k$ is the bandwidth. In this case, $K(\cdot)$ will be an Epanechnikov kernel and $p=1$. Additionally, the multiplier bootstrap approach of chiang2019robust will be used for inference and bandwidth selection, with polynomial order $q=2$ used for bias correction. More specifically, bootstrap estimates will be created as follows: $$ \phi_{b,q}^\pm = \sum_{i=1}^n u_{i,b} \frac{e_1'(\Gamma_q^\pm)^{-1} (\hat{Y}_{i,k} - r_q(g_{i,k})\hat{\beta}_{k,q}^\pm)r_q\left( \frac{g_{i,k}}{h_k} \right) K\left( \frac{g_{i,k}}{h_k} \right) \delta_{i,k}^\pm}{n h_k \hat{f}_g(0)} $$ where $u_{i,b}$ are independent and take values in $\{-1, 1\}$ (with equal probability), and $\Gamma_p^\pm$ is defined as follows: $$ \Gamma_p^+ = \int_0^1 r_p(u) r_p(u)' K(u) du, \quad \Gamma_p^- = \int_{-1}^0 r_p(u) r_p(u)' K(u) du $$

The variance of $\hat{\tau}_{\Gamma_{k, q}}$ can be computed using the variance of $\phi_{b,q}^+ - \phi_{b,q}^-$. With these components, the result below follows:\footnote{The theorem below can be trivially modified to accommodate two different bandwidths on either side of the cutoff.}

An important observation is that asymptotic normality holds even though the local linear forest estimator converges more slowly than the second-stage local polynomial regression. The key reason is that, when a true discontinuity exists, inference depends only on correctly determining the signs of the treatment effects. Because the data lie on a bounded support, the probability of sign errors can be made arbitrarily small, regardless of the slower convergence rate of the first-stage estimator. As a result, the asymptotic properties of the standard local polynomial regression remain valid.

Density Test

Estimation and inference for the aggregated univariate density are carried out using a bias-corrected kernel density estimator, with sampling variability assessed via a multiplier bootstrap. A standard kernel density estimator is given as follows: $$ \hat{f}_{g^*}(0)^+ - \hat{f}_{g^*}(0)^- = \sum_{i=1}^n \left( \frac{1}{nh_k} \right) K\left( \frac{\hat{g}_{i,k}}{h_k} \right) (2\delta_i-1) $$ where $K(\cdot)$ is a symmetric and bounded kernel function, and $\delta_i = \mathbbm{1}(\hat{g}_{i,k} \geq 0)$. Since the estimator is known to exhibit boundary bias, an approach similar to hazelton2009linear is adopted, whereby the base kernel is adjusted to achieve unbiasedness at the boundary. The resulting estimator is given by: $$ \hat{f}_{g^*}(0)^+ - \hat{f}_{g^*}(0)^- = \sum_{i=1}^n \left( \frac{1}{nh_k} \right) \left(\alpha_1 + \alpha_2 \left| \frac{\hat{g}_{i,k}}{h_k} \right| \right) K\left( \frac{\hat{g}_{i,k}}{h_k} \right) (2\delta_i-1) $$ where $\alpha_1, \alpha_2$ solves the following system of equations: $$

=

$$ With this construction of $\alpha_1$ and $\alpha_2$, the boundary bias inherent in the standard kernel density estimator is eliminated, as the modified kernel integrates to one on the relevant side of the cutoff. These coefficients also remove the first-order bias, yielding behavior analogous to that of a kernel density estimator evaluated at an interior point. The optimal bandwidth is selected using the following MSE-optimal expression: $$ h_k = \left( \frac{V_k}{B_k^2} \right)^\frac{1}{5} $$ where $V_k$ is the variance term and $B_k$ is the bias term. Their exact expressions are given as follows:

Here, the density and second derivative of the density will be estimated via preliminary local polynomial regressions (cattaneo2020simple, cattaneo2024local). For consistency of this aggregated estimator, a variation of the random forest density estimator used to construct $\hat{g}_{i,k}$ must be consistent. When joint density estimation is conducted at boundary points, the volume of each cell generated by the partitioning algorithm of wen2022random must be adjusted, particularly when the support of the running variables is non-rectangular. The effective volume of a cell is obtained by multiplying its original volume by the proportion of the cell that lies within the support of the running variables. This proportion is approximated by uniformly sampling points within the cell and computing the fraction that fall inside the support. With that, the following lemma holds:

As with the treatment effect heterogeneity case, a K-fold cross-fitting procedure is used to separate the estimation of $\hat{g}_{i,k}$ and the final test statistic. Also, as with other kernel-based nonparametric estimators, the MSE-optimal bandwidth cannot be used directly for inference because it does not adequately control for bias. To address this issue, the kernel bias-correction procedure of calonico2018effect is employed. Specifically, inference is conducted using the following re-centered estimator: $$ \frac{1}{n} \sum_{i=1}^n \frac{1}{h_k} \left[ \left(\alpha_1 + \alpha_2 \left| \frac{\hat{g}_{i,k}}{h_k} \right| \right) K\left( \frac{\hat{g}_{i,k}}{h_k} \right) - L^{(2)}\left( \frac{\hat{g}_{i,k}}{h_k} \right) \int_0^1 \frac{u^2 (\alpha_1 + \alpha_2 u)K(u)}{2!} du \right] (2\delta_i-1) $$ where $$ L(u) = 30(1-|u|)^2 u^2 \mathbbm{1}(|u| \leq 1) $$ Here, $L(u)$ is designed to estimate the curvature of the given density at the boundary point $g = 0$. With that, the following will be used to generate bootstrap replications to compute standard errors:

where $$ \hat{f}_g(0)^\pm =

, \quad \hat{f}_g”(0)^\pm =

$$

Now, the result below follows:\footnote{The theorem below can be trivially modified to accommodate two different bandwidths on either side of the cutoff.}

As in the treatment effect heterogeneity setting, the asymptotic results for the bias-corrected kernel density estimator continue to hold. This is because the probability of assigning an incorrect sign to the estimated discontinuity can be driven to zero as the sample size grows, ensuring that such errors do not affect the limiting distribution.

Accounting for Random Split Selection

Asymptotically, the specific split used in the $K$-fold cross-fitting procedure is inconsequential. In finite samples, however, different splits can lead to non-negligible variation in both point estimates and standard errors. To account for this variability, the $K$-fold cross-fitting procedure is repeated $S$ times, and the resulting main estimators are averaged across the $S$ splits. For inference, influence functions are computed for each split; multiplier weights are then applied to the average influence function, with the weights for each observation $i$ held fixed across splits. This construction parallels the clustered multiplier bootstrap extension of chiang2019robust and incorporates the randomness of the sample split into the standard error calculation.

Simulations

This section illustrates the implementation of the treatment effect heterogeneity test and the density manipulation tests using simulated data. For each case, two data-generating processes (DGPs) are considered: one that is discontinuous almost everywhere, and one that is continuous everywhere.

Treatment Effect Heterogeneity Simulations

The first DGP will be constructed as follows: $$ Y = \frac{1}{3}

+ \epsilon $$ where $\epsilon \sim N(0, \sigma_\epsilon^2)$, $X_1, X_2$ take values in $[-1, 1]$, and $f_X(x) = 0.25$ with $X_1 \perp \!\!\! \perp X_2$. Here, $\tau_\Gamma = \frac{1}{6}$. Below is a plot of this function:

For the second DGP, the data will be generated as follows: $$ Y = \frac{X_1 + X_2}{3} + \epsilon $$ where $X_1$, $X_2$, and $\epsilon$ are constructed in the same way as in the first DGP. For each of these DGPs, the following parameters are used: sample size $n = 1000$, noise variance $\sigma_{\epsilon}^2 = 0.05$, number of folds $K = 2$, number of splits $S \in \{1,10,20\}$, and 5{,}000 simulation replications. Also, two different coverage-error optimal bandwidths are used on either side of the one-dimensional cutoff. The resulting simulation outcomes are presented below.

These results reveal several patterns. First, the bias under the first DGP exceeds that under the second. This occurs because, in the second DGP, the sign of the outcome is immaterial in the absence of discontinuities. In contrast, under the first DGP, misclassification of the treatment effect sign has a non-trivial influence on the test statistic, leading to greater bias. Consequently, DGPs featuring discontinuities are expected to exhibit slightly larger bias. This bias does not seem to change on average as the number of splits increases. \\

Second, the standard errors exhibit substantial sensitivity to the number of splits. For both DGPs, the standard error is approximately 0.09 with a single split, but decreases markedly when the number of splits is increased to 10. Although the initial reduction is pronounced, further increases in the number of splits yield progressively smaller declines. For the discontinuous DGP, this reduction improves the rejection rate, whereas the rejection rate remains relatively stable for the continuous DGP.

Density Manipulation Simulations

For the first density test example, the following bivariate density will be used: $$ f_X(X_1, X_2) = \frac{1}{3}

$$ Note that this is the same function that was used in DGP 1 for the boundary heterogeneity simulations. Here, $\tau_{\Lambda}^f = 1/3$. For the second density, the following joint density will be used: $$ f_X(X_1, X_2) = 1/4 $$ where $X_1, X_2$ are independent and take values on $[-1, 1]$. Here, the first density is discontinuous almost everywhere along the treatment boundary, whereas the second density is continuous everywhere. For bandwidth selection, a single MSE optimal bandwidth is used on either side of the one-dimensional cutoff. For $n = 1000$, $K = 2$ and $S = 1,10,20$, notice the following simulation results: \\

In this scenario, the simulation results exhibit a pattern similar to that observed in the treatment effect heterogeneity case. Specifically, the bias is slightly greater under the discontinuous DGP, and the standard errors decrease up to a certain point as the number of splits increases.

Application

The empirical application in this paper follows the setting of frey2019cash. In 2003, the conditional cash transfer (CCT) program Bolsa Família (BF) was launched across Brazil. Beyond serving as a major source of income stabilization for extremely poor households, there is substantive speculation that BF also reduced the scope for political influence by diminishing incentives to engage in clientelism.\footnote{Here, clientelism refers to the act of "...replacing public good distribution with private transfers targeting groups or individual voters" (frey2019cash).} To evaluate this claim, frey2019cash exploit cross-municipal variation in CCT coverage induced by the Family Health Program (FHP). Beginning in August 2004, municipalities with fewer than 30,000 residents and a Human Development Index (HDI) below 0.7 became eligible for a 50% increase in FHP funding. This expansion influenced BF coverage by improving households’ access to information about the program. According to a survey of 10,000 poor households, more than 10% of BF beneficiaries learned about the program through their family doctor. Consequently, the political effects of BF can be identified using variation in BF coverage generated by the FHP expansion through an instrumental variables or regression discontinuity analysis. \\

In the original frey2019cash study, the author implements a multivariate regression discontinuity design along the FHP funding threshold to assess the strength of the first stage and to estimate the effects of BF on various political outcomes. However, with only a few thousand municipalities, a full boundary-based treatment effect analysis may lack reliability due to limited effective sample size near each cutoff point. To evaluate the robustness of the original findings, the methods of this paper will be applied to test for the presence of any treatment effects along the FHP boundary for the same outcomes examined in frey2019cash. For simplicity, only municipalities from 2008 are included in this analysis. \\

For the first set of results, the outcomes will be health funds received and CCT coverage. These outcomes measure the extent to which the intervention affects the type of treatment that was received. The first column reports the weighted average absolute treatment effect, with p-values in brackets. The second column presents the weighted average treatment effect based on simple signed distances (as in Theorem 1). The third column indicates whether a discontinuity is detected using the uniform bands of cattaneo2025estimation. The fourth column reports whether both positive and negative discontinuities are detected using the local polynomial regression approach of cattaneo2025estimation. All running variables are standardized using the standard deviation.

Unlike frey2019cash, the analysis presented here does not uncover evidence of a strong first stage. This conclusion holds both when applying the uniform confidence bands of cattaneo2025estimation and when using the aggregation methods developed in this paper. Together, these results indicate that the available data do not provide sufficient support for a strong first-stage effect. At the same time, there may be scope to improve the power of these tests without compromising size. As noted by crippa2025manipulation, wong2013analyzing, the scaling of the running variables can meaningfully affect the performance of distance-based aggregators. Intuitively, rescaling alters the relative emphasis placed on different segments of the boundary. A common remedy is to standardize each running variable so that all coordinates receive equal weight. However, there is no reason to expect equal weighting to be optimal for power. For instance, if discontinuities occur primarily along one portion of the boundary, assigning greater weight to that region may deliver more powerful tests than treating all boundary segments symmetrically. One possible way to address this issue is to treat the scaling parameters as tuning choices and select them to maximize power. Although this extension lies beyond the scope of this paper, it may improve the power of the proposed test and potentially reveal first-stage discontinuities that remain undetected under equal weighting. \\

For the second set of results, the analysis will use political outcomes that will help determine whether BF had any impact on clientelism in Brazil.

In this case, all but one of the political outcomes retain their significance from frey2019cash. In particular, pro-poor spending exhibits a clear discontinuity at the FHP boundary. Notably, both the local polynomial estimator and the ML-based distance aggregation approach detect this discontinuity, whereas the standard distance-based aggregation estimator does not. A closer examination of the local polynomial estimates suggests that pro-poor spending increases at some boundary locations and decreases at others. Below is a figure showing heterogeneity in the treatment effects estimated by the local polynomial regression.

The patterns observed in pro-poor spending are intuitive: larger and less-developed municipalities tend to exhibit higher pro-poor expenditures than smaller and more-developed municipalities. Under the standard distance-based aggregation method, the positive and negative discontinuities offset one another, attenuating the aggregated statistic toward zero. In contrast, the aggregation method proposed in this paper avoids this cancellation issue and is therefore able to detect a discontinuity despite the heterogeneous pattern of effects along the boundary. \\

To check for robustness, a density test will be conducted to ensure that municipalities did not manipulate their running variables to receive more funding. More specifically, a density test is conducted for four different subsets of the data, each of which were used in the main analysis. For the first subset, only municipalities in which the candidate was the incumbent mayor were retained. The second subset includes municipalities where the incumbent mayor was eligible to run for re-election. The third subset restricts the sample to municipalities in which the candidate was not the incumbent but the incumbent mayor was running for re-election. Finally, the fourth subset applies the same criteria as the third, further restricting the sample to candidates whose rank was less than three. These tests show no signs of manipulation.

Conclusion

This paper develops a testing procedure for detecting discontinuities along multidimensional regression discontinuity (RD) boundaries. A central advantage of the proposed approach is its strong performance in relatively small samples, which is a setting where many existing multivariate RD methods struggle due to the curse of dimensionality. The simulation evidence demonstrates that the procedure yields reliable inference even with limited data, and the empirical application illustrates how it can complement and strengthen insights obtained from multivariate RD estimators. \\

Several promising directions for future research remain. One important extension concerns the choice of distance metric used to aggregate local estimates along the boundary. Because the power of the global test depends on how distances are scaled, developing principled data-driven methods for tuning or learning the metric may lead to substantial power gains. Another avenue involves identifying where along the boundary discontinuities occur. Combining the proposed global test with localization techniques could provide a more refined understanding of heterogeneity in treatment effects or sorting, while limiting the disadvantages of standard multivariate RD estimators.

\printbibliography