Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
105,342 characters · 12 sections · 105 citation commands
Regression Discontinuity Design with Distribution-Valued Outcomes
\thispagestyle{empty}
\onehalfspacing
\@startsection{section}{1}{\z@} {-3.5ex \@plus -1ex \@minus -.2ex} {2.3ex \@plus.2ex} {\center}{Introduction}
The regression discontinuity design (RDD) is a popular non-experimental method for causal inference and program evaluation. It exploits a cutoff rule in the assignment variable—often a running variable such as a test score or an income threshold—to identify sharp changes in treatment status among units just above and just below the cutoff. In recent years, it has been widely adopted in economics and political science.
The conventional RDD setup typically assumes that the running variable and outcome are measured at the same level of aggregation-—each unit has its own running variable and a single outcome measurement. In many policy and program contexts, however, the outcome of interest takes the form of an entire distribution within an aggregate unit that receives treatment, rather than a single scalar. For example, when a school district implements an educational policy based on a district-wide threshold (like a cutoff in the district’s poverty rate), one may be interested in the effect on the distribution of student test scores in each district. Similarly, when a minimum wage is implemented along a state border, the outcome of interest could be the distribution of goods prices in each establishment, rather than a single average price. These settings are marked by two layers of randomness: one across units (districts, establishments), and one within (students within districts, goods sold in an establishment). This motivates the development of a more general framework, where the local average treatment effect is defined over distributions rather than scalars.
In this article, I extend the standard RDD framework to a functional data setting that can accommodate distribution-valued outcomes, allowing one to capture how an intervention shifts entire distributions rather than just their means or fixed quantiles. I call this the Regression Discontinuity Design with Distributions (R3D). Its key distinction from classical settings is that it models the data-generating process as sampling entire distributions together with the running variable. Hence, in this setting, distributions themselves are treated as random objects. This naturally leads to a novel concept of distribution-valued treatment effects, the “local average quantile treatment effect” (LAQTE), which captures the shift in the underlying average quantile function around the cutoff, where the average is with respect to the distribution of distributions. Identification is obtained by assuming that this conditional average distribution evolves smoothly. This constitutes an intuitive generalization of the canonical RD smoothness assumption to distribution-valued outcomes. The setting is illustrated in Figure (ref). In classical RDD (bottom panel), the data points are a random point cloud (rainbow colors), and their conditional expectation (gray color) is a smooth scalar-valued function. In R3D, the data points are random distributions (rainbow colors), and their conditional expectations (gray color) are a smooth path of distributions.
To estimate these average quantile treatment effects in practice, I propose two closely related estimators. The first estimator extends the canonical local polynomial regression estimator to random quantiles. The idea is to compute the observed outcome quantile function within each aggregate unit, pick a given quantile on these quantile functions, and estimate a local polynomial regression on the resulting “random quantiles”. This process is then repeated for every point on the quantile function. Such an approach accounts for both the vertical (distribution within unit) and horizontal (distribution across units) sampling that distinguishes the R3D setting from the canonical one. However, while doing this quantile-by-quantile may be intuitive, it is suboptimal in the sense that it does not properly treat the quantile function as a functional object.
Hence, I propose a second estimator, based on local Fr\'echet regression in 2-Wasserstein space petersen2019frechet, that estimates a local polynomial regression for the entire quantile function at once. Such a functional approach is preferred because it uses all information from the entire quantile function, leading to better finite-sample performance, and only requires picking a single bandwidth, which makes it more computationally efficient. Moreover, the resulting estimate is a conditional Wasserstein barycenter (Fr\'echet mean), which has the important property of being the central tendency of the observed quantile functions in probability space agueh2011barycenters, fan2024conditional. Importantly, this Fr\'echet (second) estimator is closely linked to the local polynomial (first) estimator. In particular, the Fr\'echet estimator is equivalent to an $L^2$ projection of the local polynomial (first) estimator onto the space of quantile functions. This close link between both estimators allows me to derive uniform, debiased confidence bands for both, by leveraging general theoretical results for local polynomial estimators developed in chiang2019robust. Deriving confidence bands for local Fr\'echet regression is in general not feasible due to the absence of the required algebraic structure on general metric spaces dubey2019frechet. A notable exception is petersen2021wasserstein, who derived confidence bands for global Fr\'echet regression in 2-Wasserstein space by leveraging that space's optimal transport geometry and the linearity of the global regression model. My results complement theirs by deriving the first confidence bands for the local Fr\'echet regression estimator in 2-Wasserstein space. I do so by similarly leveraging that space's optimal transport geometry without requiring a linear response model, but instead exploiting the connection to the pointwise local polynomial estimator.
Empirically, I first validate the estimators through extensive simulations. These results show that, unlike state-of-the-art quantile RD estimators qu2019uniform, the proposed R3D estimators do not suffer from asymptotic bias. Moreover, I show that the uniform confidence bands are asymptotically valid and consistent, quickly converging to the nominal 95% level and to a power of 1.
Further, I demonstrate the estimators' use in an empirical application. The question studied is, “what is the effect of partisan control of the state governor's office on the within-state income distribution?”. To answer this question, I leverage a close-election R3D design, which compares states where the Democratic candidate narrowly won their election to states where they narrowly lost. Because each state has only a single election outcome but an entire distribution of family incomes, this is a prototypical R3D setting. Applying the proposed estimators to this setting, I estimate reductions in income for above-median earners that get stronger with income and become statistically significant for the top 10 percentiles, but no such effects for lower-income families. These results point to a classical equality--efficiency tradeoff okun1975equality, where a decrease in income inequality can only be achieved at the cost of an overall loss of income.
To conclude the introduction, I note that the setting considered here is distinct from that of the quantile RD (Q-RD) setting first developed in frandsen2012quantile. That approach estimates quantile treatment effects for scalar-valued outcomes, and thus does not apply to the distribution-valued setting considered here. Indeed, in what follows, I show that the quantile RD estimator is biased and inconsistent in the R3D setting, both theoretically and in simulations. This bias results from the fact that its scalar-valued sampling framework is inappropriate for the R3D setting, and its identifying smoothness assumption highly restrictive. This can be seen in Figure (ref). Because of random sampling, the observed distributions (rainbow) exhibit discontinuous changes, violating the Q-RD assumption. The average distributions (gray) do evolve smoothly, however. Of course, when treatment and outcome are measured at the same level -- i.e., we find ourselves in the classical RD setting -- the quantile RD estimator is preferred over the R3D estimator.
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}*{Literature}
This article contributes to several strands of literature, the primary one being the literature on regression discontinuity design thistlethwaite1960regression, hahn2001identification, see lee2010regression and cattaneo2022regression for an older and more recent overview. I contribute to this large area of research in three ways.
First, I extend the literature on quantile treatment effects in RDD to allow for distribution-valued outcomes. frandsen2012quantile first developed the framework for quantile RD and derived uniform convergence results, though they did not derive uniform confidence bands. These were developed later for different types of quantile RD estimators in qu2019uniform, qu2024inference, chiang2019robust. Further variations of the classical quantile RD were studied in jin2025identification, chiang2019causal, qu2024inference, chen2020quantile. I build on this literature, in particular the general framework of chiang2019robust, to derive uniform confidence bands for distribution-valued RD designs. This also connects to the larger literature on distributional inference chernozhukov2013inference and quantile regression and treatment effects koenker1978regression, firpo2009unconditional, firpo2007efficient, chernozhukov2005iv.
Second, I contribute to the strands of literature that have developed robust, debiased confidence bands for local polynomial estimators with mean-squared error (MSE) based bandwidth selection procedures calonico2014robust, calonico2018effect, calonico2020optimal, calonico2022coverage, armstrong2018optimal, imbens2012optimal, by extending these tools to distribution-valued settings. That places this paper in a rich literature built on the foundational contributions in local polynomial regression, particularly related to bias reduction and bandwidth selection, made by fan1992variable, fan1993local, fan1995adaptive, linton1994multiplicative.
Third, this article relates to several other papers that have considered RD designs with varying levels of aggregation. borusyak2024regression considered the opposite design, where the treatment assignment is at a lower instead of a higher level of aggregation than the outcome. cattaneo2016interpreting, cattaneo2021extrapolating, bertanha2020regression considered aggregation schemes for RD with multiple cutoffs. Relatedly, gunsilius2023free, papay2011extending, cheng2023estimation considered RD designs with multi-dimensional or multiple assignment variables.
The other main strand this paper contributes to is the literature on frechet1948elements regression, which was originally developed by petersen2019frechet for general metric spaces, with several further contributions for distribution regression in Wasserstein space chen2023wasserstein, fan2022conditional, chen2023sliced, ghodrati2022distribution, zhou2024wasserstein and for local Fr\'echet regression chen2022uniform, iao2024deep, qiu2024random. As noted above, I contribute to this literature by deriving uniform confidence bands for local Fr\'echet regression in Wasserstein space, complementing related results for global Fr\'echet regression in petersen2021wasserstein and for Wasserstein barycenters in carlier2021entropic, agueh2017vers, kroshnin2021statistical. My results also hold for general polynomial orders while the literature has mostly focused on local linear regression, with the exception of schotz2022nonparametric. More broadly, this article contributes to the large literature on functional data analysis ramsay2005functional.
Relatedly, my results leverage the fact that Fr\'echet regression in 2-Wasserstein space is an $L^2$ projection of the local polynomial estimator onto the space of quantile functions. This relates closely to isotonic regression barlow1972statistical and monotone rearrangement methods, chernozhukov2010quantile, as well as shape-constrained inference with convex projection operators chetverikov2018econometrics, groeneboom2014nonparametric, fang2021projection, dumbgen2024shape.
This article also contributes to the literature applying optimal transport tools to causal inference -- see gunsilius2025primer for a recent overview. In particular, gunsilius2023distributional considered a similar setting to mine, where treatment is at a higher level than the outcome, in the context of synthetic controls (see also van2024return for an application to firm tenure distributions). kurisu2024geodesic introduced causal inference for objects in general metric spaces using geodesics, while zhou2025geodesic applied this to the well-known difference-in-differences estimator, thus complementing the distributional estimators of athey2006identification, torous2024optimal, callaway2018quantile.
Finally, this paper’s empirical application-—to gubernatorial party control and income distributions—fits into a rich literature linking partisan control of US state governments to inequality and other economic outcomes. Building on Hibbs' partisan theory hibbs1977political and Kelly’s market conditioning kelly2009politics, research has generally argued that Democrats, allied with lower income groups, adopt policies that narrow income gaps, whereas Republicans, favoring upper and business income constituencies, may widen them. Panel studies show that Democratic legislatures raise taxes and spending reed2006democrat, implying stronger redistribution. Though recent evidence from difference-in-difference designs and close-election RDDs found no evidence that party control significantly affects most state-level economic outcomes within a governor's tenure dynes2020noisy, other close-election RDs have shown that Democratic state control often increases minimum wages and welfare caseloads, compressing the post‐tax income distribution leigh2008estimating, and leads to more liberal policies caughey2017incremental. I contribute to this literature by providing credible causal estimates of the effect of gubernatorial party control on the income distribution, using rich individual-level data within each state instead of just state-level aggregates. That way, I estimate significant declines in pre-tax income for upper-income families, which compress the income distribution.
\@startsection{section}{1}{\z@} {-3.5ex \@plus -1ex \@minus -.2ex} {2.3ex \@plus.2ex} {\center}{Regression Discontinuity with Distribution-Valued Outcomes}
In this section, I present the distribution-valued version of the canonical regression discontinuity design. First, I formally introduce the setting, before providing several concrete examples from the literature. Then, I introduce a new definition of “local average quantile treatment effects” (LAQTE) appropriate for this setting, where the average is over random quantile functions. Before presenting two consistent estimators for these LAQTEs, I briefly discuss the distinction between my R3D setting and the classical quantile RD setting of frandsen2012quantile. I conclude providing an overview of the statistical inference tools developed in Section (ref), including extensions to fuzzy RDD and empirical quantile functions.
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Setting}
First, I define and discuss the R3D setting. Let $\mathcal{Y}$ be the space of cumulative distribution functions (cdfs) $G$ on $\mathbb{R}$ with finite variance, $\int_{\mathbb{R}} x^2 \mathop{}\!\mathrm{d} G(x) < \infty$. Let $(X,Y) \sim F$ be a random element with joint distribution $F$ on $\mathbb{R} \times \mathcal{Y}$. I call $X$ the running variable and $Y$ the outcome variable. Unlike the canonical RD design, here $Y$ is a random distribution rather than a random variable. Hence, each draw $(X_i, Y_i)$ from $(X,Y)$ provides a full distribution $Y_i$ at the running variable value $X_i$, rather than a single real number. Then, denote $T \in \{0,1\}$ the treatment status. I assume that $T$ is some monotonic function of $X$ such that, \[T=
\] for some threshold $c$. That is, treatment is assigned deterministically when the running variable $X$ crosses the threshold $c$, where I assume without loss of generality that $c=0$. This is the so-called “sharp” RD design, on which I focus in the main text for expositional clarity, though I derive statistical results for the fuzzy RDD case as well (see Section (ref)).
In addition, denote the marginal distributions of $X$ and $Y$ as $F_X, F_Y$. I assume that $\mu = E[X], \Sigma = \text{var}(X)$ and the conditional distributions $F_{X|Y}, F_{Y|X}$ exist with $\Sigma$ positive definite. Here, $F_{Y|X}$ is a probability measure supported on the set of cdfs $\mathcal{Y}$, $F_{Y \mid X=x}(A) \coloneqq P(Y \in A \mid X=x)$, $A \subseteq \mathcal{Y}$ with $A$ measurable.
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Motivating Examples}
To make the setting more concrete, I now provide several prominent examples from the literature that can be viewed as R3D designs. They are instances of broader classes of settings where treatment is assigned to units at a higher level of aggregation than the outcomes.
Common to all these examples is that for any value of the running variable (distance to the border, vote share, poverty level), I observe an entire distribution of the outcome (store prices, test scores, child mortality), and these outcome distributions vary across any two units (across restaurants, schools, or counties). This implies that one needs to model the outcome as a random distribution instead of a random variable, as discussed above. Consequently, new concepts of average treatment effects and discontinuities that are appropriate for random distributions are required, which I introduce in the next section.
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Local Average Quantile Treatment Effects}
To begin, I define a new treatment effects concept for distribution-valued outcomes appropriate for the setting introduced above. Following Neyman-Fisher-Rubin notation, denote $Y^0 \in \mathcal{Y}$ the counterfactual outcome distribution in the absence of treatment and $Y^1 \in \mathcal{Y}$ the outcome distribution under treatment. Define the observed outcome \[ Y =
,\] noting that $Y$ is a cdf so I can write $Y(y)$, $y\in \mathbb{R}$ to evaluate the function at a given point $y$.
Consider, for a moment, the canonical RD setting, such that $Z^T \in \mathbb{R}$ the classical scalar-valued counterfactual outcome. Then, assuming that the treatment effects vary between units, the classical local treatment effect is hahn2001identification \[ E[Z^1 - Z^0 \mid X = 0], \] the conditional expectation of the jump in the outcome variable at the threshold.
In the R3D setting, $Y^T$ is a full distribution function. An intuitive generalization of the classical average treatment effect to settings with distribution-valued outcomes is given in the below definition. Write $Q_{Y}(q)$ for the function mapping the cdf $Y$ to quantiles, \[ Q_Y(q) \coloneqq \inf \{y \in \mathbb{R}: q \leq Y(y)\}. \] Then I get,
Observe that the expectation is taken with respect to the conditional distribution of distributions, $F_{Y^T\mid X=0}$, \[ m_{T}(q) = E[Q_{Y^T}(q) \mid X=0] = \int_{\mathcal{Y}} Q_{y}(q) \mathop{}\!\mathrm{d} F_{Y^T \mid X=0}(0,y), \quad T=0,1. \] These average quantile treatment effects (AQTE) are a compelling way to summarize random distributional treatment effects. First, they offer an intuitive generalization of average treatment effects in the Euclidean setting. In particular, they allow one to study what happens to the outcome distribution of the “average” unit when it crosses the cutoff and receives treatment. Moreover, as I discuss in more detail below, they are equivalent to a difference of conditional Wasserstein barycenters, which respect the intrinsic geometry of the underlying probability measures being averaged over. In particular, the distribution defined by the LAQTEs has the intuitive interpretation of being the unique distribution with the lowest possible cumulative “least-squares” cost of transporting its probability mass into each of the underlying distributions of the individual units. This is exactly analogous to the interpretation of the mean as the “central tendency” in the standard Euclidean setting, i.e. the unique quantity that has the lowest expected least-squares distance to all points.
Next, I show that these unobserved LAQTEs can be identified from observed $(X,Y)$.
To identify $\tau^{R3D}$ from the data, I impose two assumptions that generalize the canonical RDD requirements. First, I assume that the average quantile function is continuous in the running variable around the threshold.
Importantly, this assumption allows for the observed random distributions $Y$ to evolve discontinuously with $x$, like in the top left panel of Figure (ref).
The following example may help to clarify this point. Suppose $F_{Y^T|X=x} \sim N(N(g(x) + \tau T,1),1)$ for $T=0,1$, $\tau > 0$ and $supp(X)=[-1,1]$. In words, the counterfactual distribution functions $Y^T$ are drawn from a class of normal distributions with normally distributed means that depend on $X$ and shift with treatment $T$. The distributions in Figure (ref) are an instance of this class. The figure clearly demonstrates what it means for distributions to be drawn randomly: the densities at a given value of the running variable fluctuate, leading to a lack of pointwise continuity with respect to $X$. This directly generalizes the Euclidean setting, where samples form a random point cloud that generally also lacks continuity. By contrast, (ref) shows the conditional average distributions estimated on either side of the cutoff using the local polynomial approach set out in Section (ref). These average distributions are clearly continuous in the running variable. This demonstrates that even this simple collection of random Gaussian distributions satisfies the weaker continuity assumption in (ref) but still fails continuity in quantiles. The following example establishes this formally. I include the proof here for intuition.
The example solidifies the intuition behind Figure (ref). While the probability of drawing a certain distribution varies smoothly in $X$ the actual distributions at any two points $x, x'$ close to each other will always be different with probability 1. This follows from the distributions being random objects themselves. In Section (ref) how this setting precludes the smoothness assumption used in the classical quantile RD frandsen2012quantile.
The second assumption I need for identification is a standard RDD assumption which posits no manipulation and a non-zero mass of observations around the threshold.
Then, I obtain the following identification result.
where the lemma defines $m_{\pm}(q), m(q)$.
The weak distributional continuity assumption (ref) introduced above implies that the treatment has an effect when there is a discontinuity in the observed average distribution $E[Q_Y(q) \mid X =c]$ at the threshold $X=0$. Thus, I can define a discontinuity in our setting to occur when, for some $q \in [0,1]$ \[ \lim_{x\to 0^+} E[Q_Y(q) \mid X =x] \neq \lim_{x\to 0^-} E[Q_Y(q) \mid X =x]. \] The uniform confidence bands I derive below allow one to test for the presence of such discontinuities for a given quantile $q$. Alternatively, one can conduct inference on entire segments of the distribution at once. An overview of inference is given in Section (ref).
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Comparison to Quantile RDD}
Before developing the estimators for the LAQTEs, I briefly discuss the difference between the R3D setting and the quantile RD estimator of frandsen2012quantile. The key insight is that quantile RDs are appropriate for estimating quantile treatment effects (QTE) for scalar-valued outcomes, while the R3D estimator can estimate (average) QTEs for distribution-valued outcomes, and there is no overlap in use cases.
A comparison of the population quantities targeted by each estimator makes this point clearer. Let $Y \in \mathcal{Y}$ and $Z \in \mathbb{R}$ as before. The two population objects targeted are, \[ \text{R3D}: \lim_{x \to 0^{\pm}}E[Q_Y(q) \mid X=x] \qquad \qquad \text{Q-RDD}: \lim_{x \to 0^{\pm}} E[\mathrm{1}(Z \leq z) \mid X=x]. \] Thus, the R3D aims to estimate a conditional average quantile. The Q-RDD on the other hand, aims to estimate a fixed distribution function. Practically, they do so with the following local linear estimators, \[ \text{R3D}: \frac1n \sum_{i=1}^n s_{\pm,i}(h) Q_{Y_i}(q) \qquad \qquad \text{Q-RDD}: \frac1n \sum_{i=1}^n s_{\pm,i}(h) \mathrm{1}(Z_i \leq z). \] As can be seen, the R3D approach first estimates quantiles and only then runs a local linear regression. This properly accounts for the two-level randomness intrinsic to the R3D setting. Distribution estimation at a given $X=x$ precedes smoothing. By contrast, the Q-RDD estimator intrinsically estimates the distribution by smoothing, ignoring the randomness within units. In the presence of such randomness, the observed distributions will almost surely not very smoothly, and the Q-RD approach will be biased and inconsistent.
Underlying these arguments are three distinct differences between the R3D and the Q-RDD setting. First, as mentioned, the sampling model imposed by the Q-RDD setting does not correctly represent the underlying data-generating process. In particular, it assumes i.i.d. sampling of scalar-valued outcomes instead of distribution-valued ones, which ignores the within-unit sampling that characterizes the R3D setting. As such, the sampling framework of the Q-RD design could never result in multiple data points having the same value of the (continuous) running variable. Second, as mentioned, the quantile continuity assumption required for the identification of the estimator in frandsen2012quantile is highly restrictive in the R3D setting, requiring that two units that are both close to the threshold have essentially identical distributions. In the examples in Section (ref), this would imply that, conditional on having the same value of the running variable, two different restaurants would have the exact same distributions of product prices, two different schools the same distribution of tests, and two different counties the same distribution of child mortality. Of course, there is no reason why the cheapest product in one restaurant should have the same price as in another, or the best student in one school the same score as in another, even if their running variables did happen to take on the same value. The estimator I propose requires a much weaker continuity assumption in (ref). In particular, it only demands, for example, that the test score distributions of schools near the cutoff on average look the same, while allowing the distributions of specific schools to differ. In this way, (ref) is the direct distribution-valued analogue of the conditional mean continuity assumption originally imposed in hahn2001identification, which only requires the expectation of the random outcome variable to be continuous but leaves its distribution otherwise unrestricted. Indeed, while (ref) is consistent with the common approach of averaging the outcome variable at the level of the aggregate unit and then estimating a standard RD, the continuity assumption in frandsen2012quantile is not, because there would be no random variation left in the averages, which are assumed to evolve smoothly. Third, and similarly, the standard assumption that treatment effects vary across units automatically implies that the counterfactual distributions must be random objects themselves: the outcome is a distribution, and receiving treatment affects this distribution differently for different units. More concretely: if a policy affects the entire workforce of a company, but does so differently at Company A compared to Company B, then even if all untreated companies have identical distributions in the absence of treatment (an unrealistically strong assumption), the outcome distributions of those companies under treatment will still differ.
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Estimators}
To estimate the local average quantile treatment effects introduced above, I now propose two intuitive estimators that generalize local polynomial regression to the R3D setting with distribution-valued outcomes. The first is based on the simple idea of running local polynomial regressions on the quantile functions, separately at each quantile. The second estimator builds on this by projecting the local polynomial estimator back onto the space of quantile functions. As shown in Proposition (ref), the resulting estimator coincides with the local Fr\'echet regression estimator of petersen2019frechet, restricted to the space of cumulative distribution functions equipped with the 2-Wasserstein distance (see Appendix (ref) for an overview). In Section (ref), I derive valid uniform confidence intervals for both approaches, though the Fr\'echet estimator is preferable due to its computational advantages, superior finite-sample performance, and its more meaningful interpretation as the “average” distribution.
A simple and intuitive first approach to estimating the average distributions near the threshold is to treat quantiles as the fundamental unit of observation, and estimate their conditional expectations using the local polynomial regression approach that has become canon in RDD hahn2001identification. The intuition behind the approach is illustrated in Figure (ref): regression lines are fitted through data points that represent randomly scattered quantiles.
The local polynomial R3D estimators $\hat{m}_{\pm, p}(q)$ of order $p$ for each quantile $q$ can then be written in their standard form
where $K_h(x) \coloneqq \frac1h K(x/h)$, $\delta_i^\pm \coloneqq 1\!\bigl\{X_{i}\,\mathrel{\mathpalette\@gtr@less@eq{\geqslant<}}\,c\bigr\}$, and $r_p(x) \coloneqq (1, x, x^2, \ldots, x^p)$. The only difference with the standard local polynomial RDD estimator is that I now have i.i.d. samples $(Q_Y(q), X_i)$ instead of $(Y_i, X_i)$. Standard derivations give the following solution for the conditional mean estimator,
where $s_{+,\,i n}^{(p)}(h)$ are the usual empirical weights for a local polynomial regression of order $p$ fan1996local, which I derive explicitly in Appendix (ref).
Note that the estimator $\hat{m}_{\pm,p}(q)$ is technically a function of $x$, but I suppress this for all estimators to ease notation, since I only consider the cutoff point $X=0$. Further, observe that since the weights $s_{\pm,in}(h)$ can be negative, $\hat{m}_{\pm,p}$ need not be a quantile function. To resolve this, I use the standard monotone rearrangement from chernozhukov2010quantile.
The corresponding R3D estimator then is, for each $q \in [0,1]$,
In Section (ref) below, I show that, under some assumptions, this estimator converges uniformly to an asymptotic normal distribution centered at the true treatment effect, for $p\geq 1$. Following chiang2019robust, I build bias correction into the estimator by leveraging Remark 7 in calonico2014robust, which establishes an equivalence between explicitly bias-corrected estimators and estimators where the MSE-optimal bandwidth is chosen based on a pilot estimator of lower order -- which is the approach I will take.
Three intuitive improvements can be made to the local polynomial regression on quantiles introduced above. First, as noted, the resulting function is not guaranteed to be a quantile function because the weights $s^{(p)}_{\pm,in}(h)$ can be negative and thus introduce non-monotonicity (quantile crossing). Second, the pointwise estimation approach ignores global function information, which degrades the estimator's finite-sample performance, as confirmed in the simulations below. Third, the pointwise estimation approach also requires repeated bandwidth selection and estimation for each quantile, leading to computational overhead. To resolve these three issues at once, I consider the following extension of the estimator,
where $Q(\mathcal{Y})$ is the space of quantile functions of the cdfs in $\mathcal{Y}$, restricted to $[a,b] \subseteq [0,1]$. I define $\Pi_{\mathcal{Q}}$ as the $L^2$ projection onto that space of restricted quantile functions.\footnote{Working on $[a,b]$ instead of $[0,1]$ requires much weaker assumptions on the support of the distributions and is nearly equivalent in practice, see Section (ref). } In Proposition (ref), I show that $\hat{m}_{\pm, \oplus, p}$ is unique and exists under the stated assumptions. The estimated treatment effects are then defined as,
The augmented estimator in (ref) is an $L^2$ projection of the local polynomial estimator introduced above, with the entire function projected onto the space of quantile functions. As such, it is a form of isotonic regression Robertson1988. Indeed, the approach can be viewed as a “double regression”: a local linear regression on pointwise quantile functions, followed by a functional regression on quantile functions. More importantly, due to the deep connection between $L^2$ space and the 2-Wasserstein space, this extended estimator is equivalent to the local Fr\'echet regression estimator of petersen2019frechet, restricted to the space of finite-variance probability distributions $\mathcal{Y}$ equipped with the 2-Wasserstein distance, $d_{W_2}$ (i.e. 2-Wasserstein space). In Appendix (ref), I define these objects and provide an overview of local Fr\'echet regression. Here, the main thing to note is that $\hat{m}_{\pm, \oplus,p}$ converges to the same population quantile function $m_\pm(q)$ as the local polynomial estimator. This is established in Theorem (ref) through the insight that the projection of $m_{\pm}(q)$ onto the space of quantile functions is just an identity operator, as $m_{\pm}(q)$ is a valid quantile function. Another way to view this connection is that the local Fr\'echet estimator converges to the conditional Fr\'echet mean on $(\mathcal{Y}, d_{W_2})$, which is the “conditional Wasserstein barycenter” agueh2011barycenters, fan2024conditional -- the unique distribution that has a quantile function equal to the average of the quantile functions at the cutoff, i.e., $m_\pm$ (see the proof in Proposition (ref)). In short, the Fr\'echet estimator offers a principled functional approach to estimating the LAQTE in Definition (ref), while converging to the same object in population.
This connection to local Fr\'echet regression in Wasserstein space explains why the “double regression” approach in (ref) is preferred over monotonizing the simple local linear estimator in (ref). Similar to the population object $m_{\pm}(q)$, the estimator in (ref) can be interpreted as the unique quantile function that minimizes the “quantile least squares distance” (the 2-Wasserstein distance) to each of the quantile functions $Q_{Y_i}$ in the space of probability distributions, weighted by the local regression weights $s^{(p)}_{\pm,in}(h)$. In other words, it is the weighted central tendency of the sample quantile functions $\{ Q_{Y_i} \}_{i=1}^n$ in probability space. This interpretation makes the projection approach in (ref) preferable over the mononotization approach of chernozhukov2010quantile, as the quantile function resulting from the latter generally does not have this desirable interpretation. Another advantage is that local Fr\'echet regression more naturally leverages global function information by smoothing large deviations across quantiles to minimize the objective function. In comparison, the monotone rearrangement approach just sorts the quantiles but does not otherwise use the global function information to do so in any optimal manner. These superior theoretical qualities of the Fr\'echet regression approach express themselves in better finite-sample performance in the simulations in Section (ref).
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Overview of Inference}
For these two R3D estimators (one local polynomial, one Fr\'echet), I derive the asymptotic distribution, uniformly over $q \in [a,b]$, a compact subset of $[0,1]$, in Section (ref) below. Further, I propose estimated multiplier bootstrap processes $\hat{\mathbb{G}}^{\mathrm{R3D}}$, $\hat{\mathbb{G}}^{\mathrm{F3D}}$ for the sharp and fuzzy design, respectively, that are shown to converge to the uniform limiting law and hence can be used to construct uniform confidence bands. This allows one to determine what quantiles have a statistically significant treatment effect while accounting for multiple testing due to the functional nature of the estimands.
Moreover, the bootstrapped distributions can also be used to construct critical values for various distributional hypothesis tests. In particular, treatment nullity and homogeneity can be tested in a particular part of the distribution $[\underline{q}, \overline{q}] \subset (0,1)$ through the following tests chiang2019causal:
where the critical values can be constructed by taking the $(1-\lambda)$-th quantiles of \\ $ \left\{ \max_{q \in [\underline{q}, \overline{q}]} \left| \hat{\mathbb{G}}^{\mathrm{R/F3D}'}(q) \right| \right\}_{b=1}^B$ and $\left\{ \max_{q \in [\underline{q}, \overline{q}]} \left| \hat{\mathbb{G}}^{\mathrm{R/F3D}'}(q) - \frac{1}{\overline{q}-\underline{q}}\right.\right.$[1] $\left.\left.\int_{[\underline{q}, \overline{q}]}\hat{\mathbb{G}}^{\mathrm{R/F3D}'}(q')\,dq' \right| \right\}_{b=1}^B$ with $\lambda$ the desired level of statistical significance and $B$ the number of bootstrap repetitions.
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Extensions}
So far, I have focused on the sharp regression discontinuity design, where treatment assignment is a deterministic function of the cutoff. Now, I show how to define a fuzzy R3D design in which treatment assignment is a random function of the cutoff, so only a fraction of units are treated on either side of it hahn2001identification.
Define $T^0,$ $T^1$, the local potential treatment states as $\lim_{x \to 0^\pm} T(x)$, where $T(x)$ is the potential treatment status as a function of the running variable. Further, define the events,
The treatment effects of interest are,
To identify these, I need the following standard additional assumptions,
This gives,
This Wald estimator takes the same form as the standard fuzzy RDD one hahn2001identification, except that the outcomes are random quantiles. Note that I can work with this simpler form compared to frandsen2012quantile because I work directly with the random quantiles and hence do not need to invert the CDFs on each side.
The corresponding treatment effect estimator, using local polynomial regressions of order $p$, then becomes,
where
The corresponding Fr\'echet estimator is,
So far, I have assumed that the researcher observes entire quantile functions $Q_{Y_i}$. This is realistic in settings where an entire population of sub-units within a given aggregate unit is observed -- for example, when a researcher has access to the census so that all firms within a US county are in the data. In practice, however, there is often another sampling layer, where one only observes further i.i.d. samples $Z_{ij}, j=1,\ldots,n_i$ from these distribution functions, with $Z_{ij} \in \mathbb{R}$ distributed according to $Y_i$.\footnote{See chen2023wasserstein for an analogous setting in the context of distribution-on-distribution regression, and zhou2024wasserstein in a similar setting.} In Section (ref), I show that under standard assumptions, the empirical quantile functions converge to the true quantile functions faster than the R3D estimators and hence do not affect the asymptotic results. The corresponding sharp RD estimator is defined as,
where
with \[ \widehat{Y}_i(x) \coloneqq \frac{1}{n_i} \sum_{j=1}^{n_i} \mathrm{1}\left( Z_{ij} \leq x \right), \] the empirical distribution function. The other estimators are similarly modified by plugging in $\widehat{Q}_{Y_i}$, and denoted with a bar instead of a hat, e.g.\ $\bar{\tau}_p^{\mathrm{R3D}}$. Sampling weights can be incorporated by constructing $\widehat{Q}_{Y_i}$ as weighted quantile functions. The asymptotic results for this setting are established in Section (ref).
\@startsection{section}{1}{\z@} {-3.5ex \@plus -1ex \@minus -.2ex} {2.3ex \@plus.2ex} {\center}{Statistical Results}
In this section, I derive the asymptotic distributions of the local polynomial and Fr\'echet regression estimators. I do so in full generality for $p$-th order local polynomials, accommodating both the sharp and fuzzy RDD setting. The results for the local polynomial estimator follow from an application of the general results in chiang2019robust, extended to random distribution-valued outcomes. The corresponding results for the local Fréchet regression follow from the functional delta method and a projection argument that uses a version of Rademacher's theorem for Banach spaces preiss2014gateaux. I conclude by extending these asymptotic results to the empirical quantile setting.
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Assumptions}
Throughout, I work on the restricted set of quantiles $[a,b]$, a compact subset of $(0,1)$, and let $\underline{c} < 0 < \overline{c}$. Also, define $\mathcal{Y}_c \coloneqq \{ Y(\omega): \omega \in \Omega^x, X(\omega) \in[\underline{c}, \bar{c}]\}$ as the set of random cdfs that are realized in a small neighborhood around the cutoff.
(ref) is a standard kernel assumption and is satisfied by the commonly used triangular and uniform kernels. (ref) is a standard bandwidth assumption, with the important benefit that for local polynomial order $p > 1$, it accommodates the bandwidth rates implied by common bandwidth selection procedures, which are typically slower than $h=n^{-1/5}$ calonico2014robust. Moreover, the assumption accommodates quantile-specific bandwidths. (ref) (i) is a stronger version of the standard continuity assumption (ref) that ensures the Taylor expansions required for local polynomial regression of order $p$ are well-defined. (ref) (ii) further provides some minimal control over the functional objects $E[Q_Y(q) | X=x]$ through the covariance of the quantiles. Note that both (i) and (ii) are implied by the much stronger assumption that the random distribution $F_{Y|X}$ evolves smoothly, which would be the random-distribution equivalent of Assumption E1 in frandsen2012quantile and is imposed in petersen2019frechet. For Assumption (ref), first note that clearly, for every $Y$, there exists an $M_Y$ such that $\sup_{q \in [a,b]} Q_Y(q) < M_Y$. However, the assumption strengthens this point-wise fact into a statement that these caps cannot `blow up' too often in all possible realizations $Y$. In practice, this means that while each $Y$ can have unbounded support, the family $\mathcal{Y}$ must not produce extremely large quantiles too often around the cutoff. As such, the assumption controls the across-distribution variance, enabling uniform statistical arguments. Finally, Assumption (ref) is a standard assumption for multiplier bootstraps that can easily be satisfied in practice.
For the extension with empirical quantile functions, I further impose the following,
(ref) implies that the empirical quantile functions converge uniformly van2000asymptotic. The assumption can be relaxed for the case with discrete distributions and finite support, since then standard pointwise convergence results for local linear regression imply uniform convergence fan1996local. Also note that if the entire distribution is observed, then these assumptions are not required and the quantile functions are allowed to have discontinuities. (ref) is a weak assumption on the number of measurements per distribution that guarantees the empirical quantile functions will converge faster than the estimators (ref) and (ref).
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Asymptotic Distribution}
Under these assumptions, I can derive the asymptotic distribution of the local polynomial estimator. For that, I need a few additional pieces of notation, borrowing from chiang2019robust. For the formal results, I assume without loss of generality that the kernel $K$ is supported on $[-1,1]$. Define $g_1: (Y, T, q) \subset (\mathcal{Y}, \{0,1\}, [a,b]) \to Q_Y(q)$, $g_2 : (Y, T, q) \subset (\mathcal{Y}, \{0,1\}, [a,b]) \to T$. Further, define the population residual $\mathcal{E}_k(y, t, x, q) \coloneqq g_k(y, t, q) - E[g_k(y, t, q) | X_i=x], \: k={1,2}$ and let \[ \sigma_{k l}(q, q' \mid x)=E\left[\mathcal{E}_k\left(Y_i, T_i, X_i, q\right) \cdot \mathcal{E}_l\left(Y_i, T_i, X_i, q'\right) \mid X_i=x\right] \] with $k,l \in \{1,2\}$, $q,q' \in [a,b]$, and $\sigma_{kl}(q, q' \mid 0^{\pm}) = \lim_{x \to 0^{\pm}} \sigma_{k l}(q, q' \mid x)$. Moreover, let $e_0$ denote the $0$th standard basis vector of $\mathbb{R}^p$, $(1, 0, \ldots ,0)$, and write $\Gamma_{\pm, p} \coloneqq \int_{\mathbb{R}_{\pm}} K(u) r_p(u) r_p'(u) \mathop{}\!\mathrm{d} u$. Also, let $X_n \leadsto X$ denote weak convergence for some sequence of random variables $X_n$ and a random variable $X$, while $X_n \stackrel[\xi]{p}{\leadsto} X$ denotes conditional weak convergence. The latter is defined as $\sup_{h \in BL_1} \left| E_{\xi \mid x}\left[ h(X_n) - E[h(X)] \right] \right| \stackrel[x]{p}{\to} 0$ where $BL_1$ the set of bounded Lipschitz functions with supremum norm bounded by 1 and $\stackrel[x]{p}{\to}$ denotes convergence in probability with respect to probability measure $P^x$ vaart1996weak. Then, I first get the following preliminary result for the conditional means.
Then, a simple application of the functional delta method yields the following result.
In practice, it is easier to approximate the limiting processes in Theorem (ref) with a multiplier bootstrap, which preserves the local structure without full resampling. To that end, I use the pseudo-random samples $\{ \xi_i \}_{i=1}^n$ defined in (ref) to define the estimated multiplier process,
where $\hat{f}_X(0)$ is any uniformly consistent estimator of $f_X(0)$, and $\hat{\mathcal{E}}_k\left(Y_i, T_i, X_i, q \right)$ is any uniformly consistent first-stage estimator of the residual $\mathcal{E}_k$. In practice, I will use the first-stage estimator proposed in chiang2019robust, described in detail in Appendix (ref). The process $\hat{\nu}^\pm_{\xi, n}(q, k)$ is an estimator for the uniform Bahadur representation of the bias-corrected processes $\hat{m}_{\pm,p}(q) - m_{\pm}(q)$, $\hat{m}_{\pm,T,p}(q) - m_{\pm,T}(q)$ chiang2019robust, see the proof of Theorem (ref) for more details. Then, I obtain the uniform validity of the multiplier bootstrap,
A practical algorithm for computing the empirical bootstrap is provided in Appendix (ref). The asymptotic validity and consistency of the tests proposed in Section (ref) follow immediately from Theorem (ref).
Turning to the Fr\'echet estimator, I now show that it has the same asymptotic distribution as the local polynomial estimator. I include a proof sketch to explain the intuition behind this striking result.
Note that in this theorem the $c_1(\cdot)$ term does not appear because the Fr\'echet estimator uses a single bandwidth for all quantiles. Then, the following result again follows by a simple application of the functional delta method.
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Extensions}
\@startsection{section}{1}{\z@} {-3.5ex \@plus -1ex \@minus -.2ex} {2.3ex \@plus.2ex} {\center}{Empirical Applications}
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Simulations}
To evaluate the proposed estimators’ performance, I conduct Monte Carlo simulations under several data-generating processes. Throughout this and the next section, I use R3D estimators of quadratic order but with bandwidths that are MSE-optimal for the linear estimators. As argued in Remark 7 calonico2014robust, this is equivalent to using explicitly bias-corrected linear estimators.
In the simulations, I estimate the quantile treatment effects $\tau^{\mathrm{R3D}}$ at 10 quantiles using three estimators: 1) a local polynomial estimator for classical quantile RDDs qu2019uniform;\footnote{Computed using the rd.qte command in R qu2024qte.} 2) the local polynomial R3D estimator in Section (ref); 3) the (Fr\'echet) R3D estimator in Section (ref). The Q-RD estimator is corrected for bias using the approach in qu2024inference. The reason for using the Q-RD estimator of qu2019uniform is to give Q-RD the best possible chance, since this estimator allows for bias-corrected, uniform inference, improving on the estimator in frandsen2012quantile.
I consider two data-generating processes, where $X_i \sim \mathrm{Uniform}(-1,1)$. \\ DGP 1: Normal with Normal Means. For each \(i\), draw
and define \(Y_i = N(\mu_i,\,\sigma_i^2)\).
DGP 2: Normal--Exponential Mixture with Normal--Exponential means. Set \(\mu_i = \mathrm{Uniform}(-5,5) + 2\,X_i\) and \(\lambda_i = \mathrm{Uniform}(0.5,1.5)\). Then, generate
In both setups, I let \(\Delta\) vary across different simulations to test different treatment effect magnitudes. For the first DGP, the true treatment effects have the closed-form solution $N(\Delta, 2)$, implying constant treatment effects. The heterogeneous treatment effects in the second DGP are estimated by averaging across a large number of simulated quantile functions.
Figure (ref) shows the estimators' performance in terms of relative bias, which is the magnitude of the estimated bias at a given quantile as a proportion of the treatment effect at that quantile. I set $\Delta = 2$ but the results are similar for other values. The green line (diamonds) shows the quantile RD estimator, the orange line (triangles) the Fr\'echet estimator, and the blue line (circles) the local polynomial one. In line with theoretical expectations, the quantile RD estimator appears to be inconsistent and suffers from large finite sample bias, with a relative bias that is at least an order of magnitude higher than the R3D estimators', for some quantiles. As expected, the quantile RD estimator performs well at the median in DGP 1, because a mixture of normals approximates the average normal at the median. Similarly, due to the heavy tails of the exponential distribution, it performs much worse at the upper quantiles in DGP 2. Between the two R3D estimators, the Fr\'echet estimator has much lower bias than the local polynomial one for small sample sizes, but both converge quickly to near-zero, supporting the asymptotic theory.
Figure (ref) further supports the theoretical arguments that the Fr\'echet estimator is preferred over the local polynomial one, as the latter has much larger variance in small samples, though again both estimators quickly converge to a similarly low variance. I do not report results for the quantile RD as its inferential properties are irrelevant due to its inconsistency and bias in the R3D settings.
To study the coverage properties of the confidence bands and tests proposed in (ref), I report their acceptance probabilities for both DGPs with varying values of $\Delta$ in Table (ref). The values of $\Delta$ are chosen to reflect an average Cohen's d (treatment effect size relative to standard error) of $0, 0.5,$ and $1$ which correspond roughly to no, medium, and large treatment effects. The coverage and acceptance probabilities of the uniform confidence intervals and the homogeneity test in the first two rows are not affected by the magnitude of the treatment effect. Moreover, both the Fr\'echet and the local polynomial estimator rapidly converge to the correct nominal coverage level, with the Fr\'echet estimator exhibiting slightly better coverage. The slight undercoverage in small samples is expected insofar as the estimators are only asymptotically unbiased, as also illustrated in Figure (ref). For DGP 2, which has heterogeneous treatment effects, the homogeneity test's coverage rapidly coverges to 0, illustrating the test's consistency and sharp power in finite sample. Finally, for both DGPs, the treatment nullity test also exhibits consistency and significant finite-sample power for rejecting the null hypothesis of no effect.
\@startsection{subsection}{2}{\z@} {-3.25ex\@plus -1ex \@minus -.2ex} {1.5ex \@plus .2ex} {\center}{Empirical Illustration: State Governors and the Income Distribution}
To further illustrate the method, I estimate the effect of partisan governorship on the income distribution within US states. To that end, I deploy a classical and widely used RD design in economics and political science: the close-election design lee2008randomized. This design compares constituencies where a political party barely won an election to those where it barely lost in order to estimate the effect of that party's win on some outcome of interest. The identification assumption is that the outcome of interest evolves smoothly with the party's vote share in a small window around the 50% electoral threshold that puts the party in power. Under that assumption, any jump observed in the outcome at the threshold is induced by the party's electoral win, and thus identifies its causal effect locally for states with close election outcomes. Such a close-election design naturally leads to an R3D setting (see also Motivating Example (ref)), since many outcomes of interest are measured at the constituent level, leading to an entire distribution of outcomes within each constituency.
I use data on gubernatorial election outcomes from Congressional Quarterly's Voting and Elections Collection, collating election data from 1984 to 2010. This produces a dataset of 356 state-year combinations where a gubernatorial election took place. Restricting the sample to data before 2010 ensures a stable and clearly defined environment for estimating gubernatorial impacts on state-level income distributions. The year 2010 marked a structural breakpoint in state politics (see e.g.\ the sharp increase in state-level polarization documented in shor2022two) due to the significant Republican gains from the Tea Party wave and the subsequent implementation of the Affordable Care Act (ACA). The ACA introduced confounding by influencing state policy choices through federal incentives, while increased partisan polarization changed the nature and meaning of gubernatorial party control itself. Restricting the analysis to pre-2010 thus guarantees a stable treatment definition, ensuring clearer identification of causal effects attributable specifically to Democratic versus Republican gubernatorial control. Indeed, while the magnitude of the effects remains similar when including post-2010 data, their precision and magnitude decrease (see Figure (ref)).
I combine these data with family-level income data from the UNICON extract of the March Current Population Survey (CPS) for the final year of the state governor's tenure, in order to capture the cumulative effect of that tenure on the income distribution. Practically, this means the election data is lagged 3 years, except in New Hampshire and Vermont, which hold gubernatorial elections every 2 years.
The variables in the sample are defined as follows. The running variable $X_{it}$ is the Democratic candidate's votes in state $i$ in year $t$ as a share of the combined Democratic and Republican votes. When this threshold exceeds 50%, the Democratic candidate is elected. As such, the treatment $T_{it}$ indicates whether state $i$ elected a Democratic governor in year $t$ compared to a Republican one.
The outcome variable $Z_{ijt'}$ is real income of family $j$ in state $i$ in year $t'=t+t_j$, where $t_j$ is a state-specific offset to match the income distribution in the final year of a governor's tenure to their electoral results. Real family income is constructed as the ratio of family income in year $t'$ to the federal poverty threshold in that year. Family income is defined in the standard fashion as the combined pre-tax cash income of the family, including earnings and cash transfers, but excluding non-cash benefits or tax credits. The federal poverty threshold is adjusted yearly and depends on family size and the number of children. Normalizing income by the year-specific poverty threshold makes the units of the outcome variable comparable across years, thus accounting for growth in real income levels over time and making the i.i.d.\ assumption required for the R3D estimator more likely to hold.
The CPS data are a sample of the full census data, thus placing this application in the empirical quantile setting discussed in Section (ref). In particular, instead of observing the full population income distribution, in each state $i$ in year $t$, I observe a sample of $n_i$ families $j=1,\ldots,n_i$. Based on that, I construct the empirical income quantile functions $\widehat{Q}_{Y_{it}}$, where $Y_{it}$ is the distribution function of family income in state $i$ at time $t$ such that $Z_{ijt} \sim Y_{it}$. I use the family probability weights provided in the CPS to construct these as weighted quantile functions. Further, I winsorize the distribution at the 95th percentile to account for top-coding in the CPS. In practice, I estimate the quantile function on an equally spaced grid of 95 points between $[1 \times 10^{-6}, 0.95+1 \times 10^{-6}]$, where the $1 \times 10^{-6}$ offset ensures I work on a compact subset of $[0,1]$ as required by the theoretical results.
The data are depicted in Figure (ref), which shows a version of the classical RD plot calonico2015optimal appropriate for the R3D setting, similar to Figure (ref). In particular, it shows a scatterplot of the “data”, which are the quantile functions at various quantiles $q$, averaged within equal-width bins $B_j $ of the running variable, $\frac{1}{|B_j|} \sum_{j \in B_i} \widehat{Q}_{Y_j}(q),$ with $B_j=\left\{i: X_i \in\left[x_{j, \min }, x_{j, \max }\right)\right\}$ the $j$ bins. For 5 illustrative quantiles $q$, I then fit a second-order polynomial regression line to these data. This simple descriptive plot already suggests that there is a drop in income at the higher (average) quantiles that becomes stronger as it moves up the income distribution.
Based on these data, I use the Fréchet estimator (Section (ref)) to estimate the local average quantile treatment effects in Definition (ref), plugging in the estimated empirical quantile functions $\widehat{Q}_{Y_{it}}$. For these, 90% uniform confidence bands are constructed using the bootstrap algorithm described in (ref), where I use the 90% nominal level to follow the standard in the literature frandsen2012quantile, qu2019uniform, chiang2019causal. To address some of the small-sample undercoverage reported in the simulations above, I apply the rule-of-thumb coverage correction of calonico2018effect to the IMSE-optimal bandwidth (see Appendix (ref)). In addition, I formally test for uniform treatment nullity and homogeneity using the tests described in Section (ref).
The main results are shown in Figure (ref). The graph depicts the LAQTE estimates, with the Y-axis indicating the effect as a multiple of the federal income threshold for the quantile of the distribution indicated by the X-axis. The light blue band depicts the 90% uniform confidence band.
As shown, treatment effects are slightly positive at the lowest quantiles and become increasingly negative farther up the income distribution, with the top 10 percentiles seeing a decline in income of 1.5 times the federal poverty threshold. By contrast, the very bottom quantiles see their income increase by nearly half the poverty threshold. Only the effects for the top 10 percentiles (85th--95th) are uniformly significant at the 90% level. The p-values for the uniform treatment nullity test and the treatment heterogeneity test are 0.0376 and 0.061, respectively, suggesting the observed negative relation between income quantile and effect size is significant.
The estimated results are very similar for alternative specifications with the local polynomial estimator of Section (ref), when using a uniform instead of a triangular kernel, or when using half the IMSE-optimal bandwidth in Figures (ref), (ref), and (ref). In contrast, when estimating the baseline specification with the income distribution of the same year as the election as outcome variable, none of the quantile treatment effects are significant, nor are the nullity and homogeneity tests. This suggests the results are not driven by reverse causality, where the pre-existing income distribution drives the election outcomes. This aligns with the small effects of local economic conditions on voting behavior estimated in the literature on retrospective voting healy2013retrospective. Additionally, I check whether the results are not driven by families “voting with their feet” by moving across states. To that end, Figure (ref) demonstrates that the results are near-identical when excluding families that moved across state borders in the previous year, barring some expected loss of precision.
Finally, Table (ref) reports estimates of the local average treatment effect using the standard RD estimator with robust confidence bands calonico2014robust, using both the state-level weighted average family income and the raw family-level outcome data. The state-level treatment effect estimate is $-0.631$ and significant at the 90% level, while the family-level estimate is $-0.525$ and is very precisely estimated. As predicted theoretically, the state-level estimates are in line with the average of the LAQTEs produced by the R3D estimators, which is -0.647, while the family-level estimates are 15% less strong, because they do not account for the two-level sampling when using the disaggregated family-level data. Both standard RD estimates cloak the underlying heterogeneity, in particular the redistribution that is achieved at the cost of the estimated drop in overall income. I also report the quantile RD estimator of qu2019uniform in Figure (ref). Unfortunately, the confidence bands in the companion R package are not currently implemented. However, in line with the simulations above, the estimated effects exhibit substantial bias, effectively precluding the need for inference. Specifically, the estimated quantile treatment effects are more than twice as small as the R3D estimates, and the corresponding average effect is only -0.158, 4 times smaller than the standard RD estimates using the aggregated data.
Taken together, these results suggest a classical equality--efficiency trade--off under Democratic governorship, with some redistribution of income achieved at the cost of a loss of income for upper-income earners okun1975equality. They also highlight the practical utility of the R3D estimator in uncovering distributional heterogeneity in treatment effects, compared to standard RD methods, while producing estimates that are consistent with those standard RD estimates in the aggregate.
\@startsection{section}{1}{\z@} {-3.5ex \@plus -1ex \@minus -.2ex} {2.3ex \@plus.2ex} {\center}{Conclusion}
This paper introduces the Regression Discontinuity Design with Distributions (R3D), a novel extension of the standard Regression Discontinuity Design (RDD) framework tailored to settings where the outcome of interest is a distribution rather than a single scalar value. This generalization is motivated by the common real-world setting where treatment is assigned at a higher level of aggregation than the outcome of interest, such as firm-level policies that affect employees, county-level policies that affect inhabitants, or school-level policies that affect students. Standard RD methods do not apply to such settings since they do not account for the two-level randomness involved in these settings, which introduces sampling at the level of distributions. To address this, I define the local average quantile treatment effect (LAQTE) as the primary estimand, which quantifies the difference in average quantile functions, instead of observed ones, just above and below a treatment cutoff. This measure offers a natural and intuitive extension of the traditional RDD treatment effect to distribution-valued outcomes.
To estimate the LAQTE, I propose two complementary estimators: one based on local polynomial regression applied to random quantiles and another leveraging local Fr\'{e}chet regression in 2-Wasserstein space. The local polynomial approach adapts familiar RDD techniques to handle distribution-valued data pointwise, while the Fr\'{e}chet regression method treats the quantile function as a cohesive functional object, improving efficiency and finite-sample performance. Both estimators are developed for the sharp as well as the fuzzy R3D setting. I establish the asymptotic normality of both estimators and develop uniform, debiased confidence bands that can be estimated with a multiplier bootstrap. Additionally, I introduce a data-driven bandwidth selection procedure for functional outcomes. Simulations confirm the robustness of these theoretical properties, demonstrating good finite-sample performance and reliable coverage of the confidence bands.
The practical utility of the R3D framework is illustrated through an empirical application examining the effect of gubernatorial party control on within-state income distributions in the United States, using a close-election RDD. The findings reveal a classical equality--efficiency trade-off under Democratic governorship, with some redistribution achieved at the cost of an overall loss of income. In particular, incomes at the top of the distribution decline, while slight but not statistically significant improvements are observed at the lower end of the distribution. This evidence underscores the method’s ability to uncover nuanced distributional impacts that scalar-based approaches might overlook. Moreover, the implied average effect is in line with standard RD estimates at the aggregate state level, unlike quantile RD methods, which estimate effects that are up to 4 times smaller in magnitude.
There are several avenues for future research. The R3D framework could be extended to allow for covariates jin2025identification, frolich2019including, multiple running variables or cutoffs bertanha2020regression, gunsilius2023distributional, cheng2023estimation, cattaneo2016interpreting, or multivariate outcome distributions chen2023sliced. Further, applying these methods to empirical domains where the R3D setting occurs frequently, such as education, labor policy, or politics, promises to yield new insights into the distributional consequences of policy interventions.
In summary, the R3D framework offers a powerful and versatile new tool for causal inference with functional outcomes, making the estimation of distributional treatment effects practical in a novel but commonly occurring setting. By providing both theoretical foundations and practical estimation strategies, this article equips researchers with a new way to address pressing questions about how policies shape distributions.