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.
89,619 characters · 28 sections · 49 citation commands
Difference-in-Discontinuities: Estimation, Inference and Validity Tests
{\bf Keywords:} Difference-in-discontinuities, regression discontinuity design, difference-in-differences
The difference-in-discontinuities design (DiDC) aims to address the limitations of both regression discontinuity designs (RDD) and difference-in-differences designs (DiD) by combining temporal and discontinuity-based sources of variation from the data-generating process (DGP). grembi2016fiscal and eggers2018regression proposed this quasi-experimental approach, using differences between the pre- and post-treatment periods around a threshold for both the treated and the untreated groups. Despite its increasing use in applied microeconomics \citep*{azuaga2017violencia, chicoine2017homicides, garcia2022plant, albright2024hidden}, the econometric theory of DiDC (identification, estimation, and asymptotic behavior) remains limited. leventer2025correcting is one of the few papers that study the identification of regression discontinuity in the case of violation of the continuity of potential outcomes assumption and multiple periods. Their estimator can be an alternative to the DiDC.
The method offers more flexibility over the standard RDD in some settings and can often be used in contexts where neither the RDD nor the DiD is applicable. In particular, the DiDC can handle cases where the control and treatment groups differ significantly and do not satisfy the parallel trends assumption of the DiD, or when multiple time-invariant confounders are at or near the threshold in an RDD setting. Additionally, by incorporating more information into the estimation, the DiDC eliminates bias in RDD estimates, providing more accurate and reliable estimates of treatment effects under certain assumptions.
The advantages of the DiDC have been explored in several applied microeconomic studies. For example, azuaga2017violencia use the DiDC to analyze the impact of introducing legislation on domestic violence in Brazil. butts2023geographic explore the use of Diff-in-Disc in geographic settings, chicoine2017homicides investigates the expiration of the Assault Weapon Ban (AWB) by comparing municipalities in which the incumbent mayoral party wins a close election with those where the incumbent is defeated, and albright2024hidden isolates the causal effects of algorithmic recommendations on decision-makers using DiDC. These real-life examples demonstrate how effective this approach can be in various research settings, highlighting its value as a practical analysis method.
In this paper, we develop the econometric theory for the Difference-in-Discontinuity design. While galindo2021fuzzy develop an identification theory for the fuzzy difference-in-discontinuity design based on the difference of RDD estimations, this work will focus on the sharp design, developing identification, inference and validity tests for the DiDC design. leventer2025correcting also provide an alternative method to the DiDC. In fact, they use the time to address violations of the RDD's continuity assumption (e.g., the confounding effect at the threshold). However, they do not specifically study the DiDC estimator and compare it to the traditional RDD.
Drawing from standard assumptions in cross-sectional RDD and assumptions specific to the DiDC framework, we establish conditions for reliable identification. These new assumptions, specific to the DiDC design, constrain how confounding effects behave around the threshold and over time.
We show that the parameter of interest can be recovered using local polynomial estimation of the differences in the outcomes employing the methodology proposed by calonico2014robust. We derive the asymptotic properties of the DiDC estimator and highlight scenarios in which its asymptotic bias can be smaller than that of the RDD. In addition, we introduce simple tests to assess the testable implications of the identifying assumptions. We explore the conditions under which DiDC mitigates bias more effectively than RDD, offering practical guidance for unbiased treatment-effect estimation. We also provide partial identification results for settings in which the identifying assumptions fail to hold.
Through Monte Carlo simulations, we evaluate the finite-sample properties of the estimator, comparing its performance to that of the local linear RDD estimator proposed by calonico2014robust and the nonparametric DiD regression estimator proposed by santanna2020doubly. By comparing DiDC's performance with DiD and RDD methods, we identify the scenarios where DiDC outperforms scenarios where DiDC outperforms RDDs and contexts for each method. Finally, we apply our estimator to the study by grembi2016fiscal on the impact of fiscal rules on municipal deficits in Italy.
This work is closely connected to calonico2014robust, whose contributions have been instrumental in robust nonparametric RDD estimation and whose estimation methods are used for our DiDC approach. frolich2019impact discusses the potential intersection between RDD and DiD, but does not provide a formal derivation of the econometric properties, offering only high-level remarks on identification and estimation. This work also aligns with the broader literature on the intersection of RDD and panel data. pettersson2012does explore settings where RDD is combined with fixed effects to address small sample issues and violations of the continuous support assumption. lemieux2008incentive, use a first-difference RD approach to eliminate individual-specific fixed effects by capitalizing on the longitudinal nature of the Finnish Census data. Lastly, cellini2010value introduce "dynamic RD" models, accommodating scenarios with multiple treatment opportunities and examining the dynamics of treatment effects.
The remainder of this article is organized as follows. Section (ref) presents the DiDC as the discontinuity of differences in a potential outcomes model and shows the main identification results, along with the necessary assumptions. Estimation procedures for the treatment effect in the sharp setting are presented in Section (ref), along with the derivation of large-sample properties, optimal bandwidths and robust confidence intervals. In Section (ref) we present tests for 2 of the identifying assumptions and in Section (ref) we consider the partial identification of treatment effects when the assumptions described in Section (ref) are violated. Monte Carlo simulations are conducted in Section (ref) to examine the properties of the estimators. Section (ref) provides an empirical illustration and Section (ref) concludes. In the \hyperref[app:estimation_setup]{Appendix} we provide detailed notation, proofs and other methodological results.
Our model is best suited for panel data. We consider a setting with two time periods ($t\in\left\{0,1\right\}$) and $n$ units, indexed by $i$. For each unit, we observe an outcome $Y_{i,t}$ at each time period and a running variable $Z_{i}$, which is assumed to be time-invariant.
We are interested in measuring the effect of a binary treatment introduced between time periods 0 and 1, denoted as an indicator function $D_{i,1}$. Treatment assignment is a function of the running variable $Z_{i}$. Units with the running variable above a cutoff $z_{0}$ are assigned to treatment, whereas units with the running variable below $z_{0}$ are assigned to control. That is, the assignment rule for treatment $D_{i,1}$ can be formalized as $D_{i,1}=\mathbf{1}\left\{Z_{i}\geq z_{0},t=1\right\}$.
Our setting differs from the canonical RD setting due to the presence of a confounding policy, denoted by $D_{i,0}$, which is introduced before the treatment of interest, but following the same assignment rule around the cutoff $z_{0}$. The assignment rule for the confounding policy is $D_{i,0}=\mathbf{1}\left\{Z_{i}\geq z_{0},t\geq 0\right\}$.
For concreteness, and in anticipation of the empirical analysis, let $Y_{it}$ denote a fiscal outcome of municipality $i$ in period $t$. The running variable $Z_{i}$ is the population of municipality $i$, the confounding policy $D_{i,0}$ is a high wage for the mayor of municipality $i$, and the treatment of interest $D_{i,1}$ is a relaxation in the fiscal rule for municipality $i$.
Potential outcomes are defined as a function of both the treatment of interest and the confounding policy. We define $Y_{i,t}(d_{0},d_{1})$ as the potential outcome for unit $i$ in period $t$ if the confounder and the treatment are set to $d_{0}$ and $d_{1}$, where $(d_{0},d_{1})\in\left\{0,1\right\}^{2}$.
The observed outcomes at each time period are related to the potential outcomes through
and
Since the policy of interest is only introduced between periods 0 and 1, at period 0 we only observe two of the four possible potential outcomes, whereas in period 1 we can observe all four possible potential outcomes.
The presence of the confounding policy at the threshold is a violation of the continuity of mean potential outcomes assumption invoked in the canonical RD setting (Assumption 1 from hahn2001identification). In the DiDC setting, the continuity assumption is modified to accommodate the confounding policy:
Assumption (ref) is a straightforward modification of the standard continuity assumption. It can be interpreted as stating that, other than the confounder $D_{0}$, there is no systematic difference in mean potential outcomes around the threshold $z_{0}$. In the empirical application, Assumption 1 states that, aside from the wages of mayors and the relaxation of fiscal rules, there are no other policies that change around the population threshold $z_{0}$.
The second identification assumption is the random treatment assignment around the cut-off that is known as non-manipulation at the thresold:
A fundamental assumption in the canonical RD setting is that the probability of receiving treatment is discontinuous at the threshold $z_{0}$. In the DiDC, this assumption is extended for the treatment of interest and the confounding policy:
Assumption (ref) states that there is perfect compliance towards the assignment rule around the threshold $z_{0}$. Borrowing the jargon from the standard RD setting, we call the case with perfect compliance the sharp discontinuity setting. In the sharp setting, we observe $Y_{i,0}=Y_{i,0}(1,0)$ and $Y_{i,1}=Y_{i,1}(1,1)$ for units above the threshold $z_{0}$, whereas we observe $Y_{i,0}=Y_{i,0}(0,0)$ and $Y_{i,1}=Y_{i,1}(0,0)$ for units below the threshold.
Assumptions (ref) and (ref) modify the traditional RDD assumptions stated at hahn2001identification to accommodate the confounding policy. However, to identify the effect of the treatment of interest, an additional assumption regarding the evolution of the confounder effect is required:
Assumption (ref) states that the effect of confounding policy $D_{i,0}$ remains constant over time. That is, any differences between the treatment and control groups at the threshold in $t=1$, not caused by the treatment, should have existed in $ t=0$ before the treatment was introduced.
Assumptions (ref)-(ref) are the assumptions invoked by grembi2016fiscal in order to derive a causal interpretation for the DiDC estimand. In the next section, we discuss the Difference-in-Discontinuities estimand, as well as the Discontinuity-in-Differences estimand.
Before discussing identification, we define the relevant target parameters for the setting. There are two relevant causal parameters for understanding the effect of the policy of interest. We are interested in the identification of the average difference between potential outcomes of treated units at the cutoff and the potential outcomes of untreated units, holding the confounder fixed.
Thus, the first parameter we define is $\tau_{c}\equiv\mathbb{E} \left[Y_{i,1}(1,1)-Y_{i,1}(1,0)|Z_{i}=z_{0}\right]$, which is the average treatment effect for individuals at the cutoff exposed to the confounder in period 1. The second parameter of interest is the average treatment effect for individuals at the cutoff not exposed to the confounder in period 1, defined as $\tau_{uc}\equiv\mathbb{E} \left[Y_{i,1}(0,1)-Y_{i,1}(0,0)|Z_{i}=z_{0}\right]$. We now discuss the point identification of these parameters.
We define $\tau_{t}^{RD}$ as the regression discontinuity estimand at period $t$. Let $Y_{t}^{+}=\lim_{\varepsilon \rightarrow 0} E[Y_{i,t} |Z_{i}=z_0 + \epsilon]$ and $Y_{t}^{-}=\lim_{\varepsilon \rightarrow 0} E[Y_{i,t} |Z_{i}=z_0 - \epsilon]$. Thus, the RD estimand at period $t$ is simply $\tau_{t}^{RD}=Y_{t}^{+}-Y_{t}^{-}$. grembi2016fiscal define the DiDC estimand as the difference between the RD estimand in period 1 and the RD estimand in period 0: $\tau^{DiDC}=\tau_{1}^{RD}-\tau_{0}^{RD}$.
The next lemma shows that the DiDC estimand can be alternatively defined as an RD estimand evaluated at the difference between outcomes over time:
The result in Lemma (ref) is intuitive and presents an attractive feature for the DiDC setting: the estimand can be implemented via a single RD estimand rather than taking the difference across estimands. We now turn to the causal interpretation of the estimand, which is already well established in the literature:
The difference-in-discontinuities estimand identifies the effect of the treatment of interest in period 1 for units at the threshold that are exposed to the confounding policy. To identify a more general causal effect, an additional assumption is required:
Assumption (ref) states that the effect of the treatment of interest does not depend on the confounding policy. Such an assumption can be justified by a potential outcomes model in which the treatment effect of interest and the effect of the confounder are linearly separable. Assumption (ref) might be overly restrictive, as one might expect the treatment of interest and the confounding policy to interact; nevertheless, it allows the DiDC estimand to identify a more general causal effect:
Assumption (ref) can be overly restrictive, and sometimes there could be scenarios where the confounding policy and the treatment of interest might interact. Rather than assuming additive effects at time $t=1$, we can relax this assumption with one that allows the inclusion of multiplicative effects:
This assumption allows for the possibility that the combined effect of the confounding policy and the treatment of interest is not simply the sum of their individual effects, but also includes an additional component reflecting their interaction. This additional component is determined by multiplying the effects together. By incorporating this concept, we can derive a new estimand that uses both the RDD at time 0 and the DiDC to identify the causal effect of the treatment of interest.
This estimand can capture more complex relationships between the treatment and the confounding factors. Our future work will focus on the properties and robustness of this estimator.
Following standard practice in the regression discontinuity literature, we propose a local polynomial estimation to recover the parameter of interest. This nonparametric method involves fitting a polynomial to the data near the threshold and using the estimated function to calculate differences in outcomes between the treatment and control groups at that threshold.
To implement the estimation at $z_0$, we rely on the methodologies proposed by calonico2014robust for local polynomial estimation of the RDD. In our case, we estimate a local polynomial regression of the differences in outcomes over time ($\Delta Y_{i}$).
Under a mild continuity condition, hahn2001identification showed that the average treatment effect at the threshold is nonparametrically identifiable as the difference of two conditional expectations evaluated at the (induced) boundary point $z_0=0$. Similarly, the sharp DiDC parameter can be identified as the difference of the difference (in time) of two conditional expectations evaluated at $z_0=0$ at each side of $z_0$:
Appendix (ref) states the assumptions underlying the nonparametric local polynomial regression estimation. These assumptions impose restrictions on the kernel function, require the existence of certain moments, ensure continuity of the running variable in the relevant region, impose smoothness conditions on the regression functions, and bound the conditional variance of the observed outcome.
For a given $\nu \leq p \in \mathbf{N}$, define $\Delta \mu^{(\nu)}$ as the $\nu$th-order derivatives of the $p$th-order local polynomial of the difference. We are interested in the limits of this function around the threshold, $\Delta \mu^{(\nu)}_{+}$ and $\Delta \mu^{(\nu)}_{-}$. The general estimand of interest is $\tau^{DiDC} = \Delta \mu_{+} - \Delta \mu_{-}$, with $\Delta \mu=\Delta \mu^{0}$. The $p$th-order local polynomial estimators of the $\nu$th-order derivatives $\Delta \mu^{(\nu)}_{+,p}$ and $\Delta \mu^{(\nu)}_{-,p}$ are:
where $e_\nu$ is a conformable $(\nu +1)$ unit vector, $K_h(u) = K(u/h)/h$, $h_n$ is a positive bandwidth sequence, $r_p(x) =
'$ and $ \Delta Y =
'$. Therefore, for a positive bandwidth $h_n$, the nonparametric estimator of $\tau_{\nu,p}$ is
To perform local polynomial estimation, it is necessary to choose an appropriate bandwidth $h_n$. This parameter determines the range of observations used for the estimation and impacts the trade-off between bias and variance in the estimated treatment effect $(\tau^{DiDC})$. In point estimation, the standard approach is to select the bandwidth that minimizes the asymptotic Mean Squared Error (MSE) of the estimator. Let $\chi_n =
'$, the MSE is:
We now discuss the asymptotic properties of the estimator. Using Lemma (ref) found in the Appendix (ref), it is possible to recover the leading asymptotic bias, expressed as:
where $\mathcal{B}_{+,\nu,p,r} = \nu!e_\nu'\Gamma_{+,p}^{-1}(h_n)\vartheta_{+,p,r}(h_n)$ and $\mathcal{B}_{-,\nu,p,r} = \nu!e_\nu'\Gamma_{-,p}^{-1}(h_n)\vartheta_{-,p,r}(h_n)$ are asymptotically bounded. Further notation is available in Appendix (ref) and a detailed proof for the statement above can be found in Appendix (ref) under Lemma (ref).
This is where we believe one of the main contributions of this research lies: the asymptotic bias of the DiDC can be zero if we include an assumption similar to that of parallel trends, or if the shapes of the data-generating processes for both groups are time-invariant. In the first case, imposition of “parallel trends“ restricts the functional form in such a way that $\Delta \mu_{+}^{(r)}\mathcal{B}_{+,\nu,p,r}(h_n) = \Delta \mu_{-}^{(r)}\mathcal{B}_{-,\nu,p,r}(h_n)$. It's important to note that the symmetry of the kernel function, as imposed in Assumption (ref), plays a significant role in this result. In the case of time-invariant data-generating processes, both $\Delta \mu_+^{(\nu)}$ and $\Delta\mu_-^{(\nu)}$ equate to zero.
We present Claim (ref), demonstrating the bias of difference-in-discontinuities in relation to RDDs:
Equation (ref) demonstrates how incorporating additional data can help reduce the bias in RD analysis. Traditional RD analysis of interventions typically uses only a cross-section of post-intervention data. However, by incorporating pre-intervention data into the analysis, bias can be substantially reduced and, under specific conditions, eliminated. Even when the bias is not completely eliminated, it can still be reduced if the RD estimation of pre-treatment data exhibits bias in the same direction as the post-treatment RD and if its magnitude is not too large to overcome the original bias.
Using the MSE-optimal bandwidth for point inference can result in bandwidths that are “too large“, which may seem attractive for reducing variance but can introduce first-order asymptotic bias\footnote{calonico2014robust}. To address this issue, we employ robust bias-corrected confidence intervals (CIs) proposed by calonico2014robust. These intervals take into account the asymptotic bias of the point estimate by 1) estimating the bias and recentering the CI, and 2) incorporating the additional variance from estimating the bias for bias correction into the CI. This process requires estimating a separate local polynomial of order $q$, with $q > p \geq \nu$. The bias-corrected estimator is defined as:
with $\Delta \hat{\mu}^{(p+1)}_{+,q}(b_n)$ and $\Delta \hat{\mu}^{(p+1)}_{-,q}(b_n)$ being local polynomial estimations as described in Appendix (ref). $\hat{\mathbf{B}}_{\nu,p,q} \left(h_n, b_n \right)$ is the estimationg we get of the bias from the $q$th-order local polynomial.
To determine the MSE-optimal bandwidth for estimating the bias, we need to conduct a separate local polynomial estimation of order $q$ with $q > p \geq \nu$. Once again, we aim to minimize the MSE. Following Lemma (ref) and applying it to the bias estimate, we can find the MSE-optimal bandwidth for the bias estimation.
Following calonico2014robust, we derive a large-sample distributional approximation that accounts for the added variability introduced by the bias estimate. The large-sample approximation for the standardized t-statistic is formalized in the theorem below:
This motivates the following confidence interval:
Just like in the standard RD setting, this confidence interval have better properties compared to conventional bias-corrected intervals. They are more robust to bandwidth selection, feature coverage error decays at a faster rate, and offer shorter interval lengths, as explained by calonico2014robust.
As highlighted in Section (ref), the importance of Assumption (ref) - that the confounding effect is constant over time - cannot be overstated. Without this assumption, all estimations are inherently flawed, leading to biased treatment effect estimates. Therefore, an initial step when considering the applicability of a differences-in-discontinuity approach is to assess the validity of this assumption.
To accomplish this, we propose a simple test to check violations of this assumption. In essence, it involves using available pre-treatment periods to estimate stacked RDs. The goal is to examine whether the RD coefficients remain consistent across multiple periods. This involves a set of $k$ periods where it is known that no changes occurred at the threshold.
It is crucial to note that this estimation method is not suitable if any alterations occurred at the threshold between the initial period in this sample and the period just preceding the introduction of the treatment of interest. To ensure the validity of this approach, it is necessary to select a period during which no events occurred at this threshold that could influence the observed outcome.
The procedure involves considering the following stacked RDs regression model:
where $T_{-k}$ is a dummy indicating the period to which the data belong, that is, the RD of which period. The hypothesis to be tested is:
A Wald test is sufficient to test these hypotheses. Rejection of the null hypothesis indicates that the differences-in-discontinuity design may not be suitable for estimating the effect of treatment in this setting.
A question arises about concerning bandwidth to use when conducting this test. There are several options available, and Appendix (ref) provides details on the simulations used to assess the better bandwidth. In summary, in practice, any bandwidth derived from the data that is optimal for the specific RD should work well.
In Section (ref), we briefly discussed an important aspect of the DiDC method: its potential to achieve zero asymptotic bias when choosing DiDC over the standard RD, provided that the functional forms of the data-generating processes for both groups are time-invariant. Essentially, this would mean that the derivatives of the functions on each side of the threshold at each $t=\{0,1\}$ ($\mu_{+,1}(z)$ and $\mu_{+,0}(z)$, $\mu_{-,1}(z)$ and $\mu_{-,0}(z)$) are equal at each point of the running variable $z$, resulting in zero bias in the local polynomial estimation.
When $\Delta \mu_+^{(\nu)} = \Delta \mu_-^{(\nu)} = 0$, the bias term $\mathbf{B}_{\nu,p,r} (h_n) = \frac{\Delta \mu_{+}^{(r)}\mathcal{B}_{+,\nu,p,r}(h_n) - \Delta \mu_{-}^{(r)}\mathcal{B}_{-,\nu,p,r}(h_n) }{ r!}$ is also zero. This condition implies that $\Delta \mu(z)$ is a constant, allowing unbiased point estimates through linear estimations on both sides of the threshold. As a result, two simple OLS regressions would be sufficient for accurately estimating the treatment effect.
Therefore, it can be very beneficial for researchers to know if the shapes of the data-generating functions remain stable over time. This knowledge would allow them to determine whether they can use simpler, less biased estimation methods. We propose a simple nonparametric Two-Sample Kolmogorov-Smirnov (KS) test for time-invariance of the functional forms of conditional means. The procedure is an adaptation of the KS test for conditional moment restrictions from Whang (2001).
In order to implement the test, it is required to have data from more than one time period before the implementation of the treatment of interest ($t\in\left\{-1,0,1\right\}$). Using the notation from Appendix A, we write the nonparametric regression model for the outcome in time $t$ as
where $\mu_{t}(Z_{i})=\mathbb{E}\left[Y_{i,t}|Z_{i}\right]$ and $\varepsilon_{i,t}=Y_{i,t}-\mathbb{E}\left[Y_{i,t}|Z_{i}\right]$. Let $\mathcal{S}^{+}(Z)$ and $\mathcal{S}^{-}(X)$ denote, respectively, the support of the running variable above and below the threshold $z_{0}$. We want to test whether $\mathbb{E}\left[\mu_{0}(Z_{i})-\mu_{-1}(Z_{i})|Z_{i}\right]=0$ almost surely for $z\in\mathcal{S}^{+}(Z)$ and whether $\mathbb{E}\left[\mu_{0}(Z_{i})-\mu_{-1}(Z_{i})|Z_{i}\right]=0$ almost surely for $z\in\mathcal{S}^{-}(Z)$. To do so, we assume the conditional mean of $Y_{i,t}$ can be approximated by a $K\times1$ vector of approximating functions $p^{K}(.)=(p_{1}(.),...,p_{K}(.))^{'}$. Using this vector, we write the regression model for the outcome in period $t$ above and below the threshold, respectively, as
Under the null hypothesis of time-invariant conditional means, we have
which motivates the following KS type test statistics:
where $n^{+}=\sum_{i=1}^{n}\mathbf{1}\left\{Z_{i}\geq z_{0}\right\}$, $n^{-}=\sum_{i=1}^{n}\mathbf{1}\left\{Z_{i}< z_{0}\right\}$, $\widehat{\gamma}^{+}_{t}$ is the vector of least square estimates of the regression of $Y_{i,t}$ on $p^{K}(Z_{i})^{'}$ for units above the threshold and $\widehat{\gamma}^{-}_{t}$ is the vector of least square estimates of the regression of $Y_{i,t}$ on $p^{K}(Z_{i})^{'}$ for units below the threshold.
Confidence intervals and p-values for the test statistics can be obtained using a recentered bootstrap hall1996bootstrap in which the null hypothesis is imposed for the bootstrapped distribution. In Appendix F, we provide the conditions under which the bootstrap is consistent, and the test statistics are powerful against local alternatives.
In this section, we consider the partial identification of treatment effects when Assumption (ref) is violated. In many settings, assuming that confounding effects remain constant over time is implausible. If that is the case, then the effect of the treatment of interest is not point identified.
We replace the time-invariance assumption with bounded variations assumptions, in the spirit of manski2018how, and derive identified sets for the causal effect of interest as a function of the sensitivity parameters that bound variations in potential outcomes across time. Formally, we invoke the following assumption:
The first inequality in Assumption (ref) states that the difference between the mean potential outcome $Y_{i,0}(0,1)$ (observed) and the mean potential outcome $Y_{i,1}(1,0)$ (not observed) at the threshold is no greater than the scalar $c_{1}$. Thus, the parameter $c_{1}$ can be interpreted as a time-trend in the evolution of the potential outcome associated to receiving the confounder, but not the treatment of interest.
The second inequality states that the difference between the mean confounding effect in period 0 (point identified under Assumptions (ref) and (ref)), and the mean confounding effect in period 1 (not observed nor identified) at the threshold is no greater than the scalar $c_{2}$. In that sense, $c_{2}$ can be interpreted as a time-trend in confounding effects. In the next lemma, we derive the identified set for $\tau_{c}$:
Lemma (ref) shows that average treatment effect for confounded individuals at the cutoff is partially identified as a function of the sensitivity parameters $c_{1}$ and $c_{2}$. It is straightforward to identify breakdown values for this identified set. That is, the largest values of $c_{1}$ and $c_{2}$ for which a particular conclusion holds. For instance, researchers might be interested in the largest values of bounded variation under which one can conclude that the causal effect of interest is positive. If the lower bound is above 0 for large values of $c_{1}$ and $c_{2}$, then one can conclude that the qualitative conclusions of the empirical exercise are robust to violations of the identifying assumptions.
Estimation of the bounds for $\tau$ is straightforward. We a nonparametric estimator in which the local polynomial estimates described in Section 3 are plugged in the max and min operators. The estimators for the lower and the upper bound are, respectively,
We use the Delta method for Hadamard directionally differentiable mappings fang2018inference to derive the asymptotic properties of the plug-in estimators for the bounds. The next theorem shows that the estimators of the bounds converge to a non-Gaussian limiting process at the $\sqrt{n\min\left\{h_{n},b_{n}\right\}}$ rate:
Inference for the bounds estimates is particularly challenging. Using the limiting process to which the estimates converge to obtain analytical asymptotic confidence bands is difficult. An alternative would be a bootstrap. Although the Delta Method is valid for Hadamard directionally differentiable functions, the standard nonparametric bootstrap is not. Instead, we recommend researchers to use the bootstrap procedure proposed in Section(ref) in fang2018inference\footnote{Note, however, that the mapping is fully differentiable if the quantities within the min and max operators are different, in which case the standard bootstrap is valid}.
An alternative approach to partial identification in the Diff-in-Disc setting is to partially identify treatment effects by exploiting assumptions regarding the interaction between the treatment of interest and the confounding policy. The first step is to assume that potential outcomes are bounded up to known values:
In some settings the values $y^{min}$ and $y^{max}$ readily justified by the nature of the outcome, and in others these values are sensitivity parameters that reflect the beliefs about the smallest and largest possible values of potential outcomes kline2025finite.
Assumption (ref) motivates the worst-case bounds manski1989anatomy. Without further assumptions, it implies that the treatment effects lie in the interval $\left[y^{min}-y^{max},y^{max}-y^{min}\right]$. Such bounds are usually uninformative as the bound covers zero, and hence the sign of the treatment effect is not identified.
Theoretical knowledge of the applied setting can be used to tighten the bounds. For instance, researchers can draw insights from economic theory to assume that the treatment of interest and the confounding policy are complementary:
Assumption (ref) is the supermodularity assumption from twinam2017complementarity for the case of two binary treatments. It formalizes the notion that the magnitude of the causal effect of interest increases with the confounding policy. Mapping the assumption to our empirical setting, one can assume, for instance, that the effect of relaxing fiscal rules on public finance outcomes is greater for municipalities with less-skilled incumbents. The proposition below derives the sharp identified sets for $\tau_{c}$ and $\tau_{uc}$ under complementarity:
Alternatively, researchers can assume that the treatment of interest and the confounder are substitutes:
Assumption (ref) is the submodularity assumption from twinam2017complementarity for the case of two binary treatments. It states that the treatment effect of interest decreases in the confounding policy. The proposition below derives the identified sets for $\tau_{c}$ and $\tau_{uc}$ under substitutability:
So far in this section, we have described two approaches to partially identify the target parameters in the DiDC setting under different sets of assumptions.
The bounded variation approach partially identifies $\tau_{c}$ as a function of a measure of deviation from Assumption (ref). The procedure has several desirable features. It allows for researchers to partially identify the parameters under different values for the deviation, which also allows researchers to assess the robustness of empirical findings. However, it does not allow for the partial identification of $\tau_{uc}$, not at least without further assumptions regarding the homogeneity of treatment effects.
The modularity approach, on the other hand, allows for the partial identification of $\tau_{c}$ ad $\tau_{uc}$, by invoking assumptions on the interaction of the confounding policy and the treatment of interest which can be rooted in economic theory or background knowledge of the empirical setting. However, the procedure does not allow researchers to exploit the temporal structure of the data. Nevertheless, it is (particularly) useful for partial identification in RD settings in which there is a confounding policy at the threshold, but no data from periods before the implementation of the treatment of interest.
In the lemmas below, we show that combining bounded variation and modularity assumptions can strengthen the bounds on $\tau_{uc}$ and $\tau_{c}$:
Lemmas (ref) and (ref) derive the identified sets for $\tau_{c}$ and $\tau_{uc}$ when bounded variation, bounded outcomes and modularity assumptions are combined. When it comes to $\tau_{c}$, the combination of assumptions strengthens the lower and the upper bound under both modularity assumptions. When it comes to $\tau_{uc}$, however, only one of the bounds is strengthened by the combination of assumptions, whereas the other bounds remains the worst-case bound. When it comes to estimation and inference, the bootstrap procedure outlined in Theorem 2 remains valid for bounds under alternative sets of assumptions as well.
We now analyze the finite-sample properties of the Diff-in-Disc estimator through Monte Carlo simulations and compare the performance of the proposed estimator to that of the local linear RDD estimator proposed by calonico2014robust, the nonparametric DiD regression estimator proposed by santanna2020doubly, and the Two-Way Fixed Effects (TWFE) estimator. We consider Data Generating Processes (DGP) based on model 3 from calonico2014robust, with small modifications that will be described.
We conduct our simulation studies in four distinct settings: with time-invariant confounding factors at the threshold and without any confounding factors, which is the typical scenario for RDDs, for both time-invariant and time-varying functional forms of the conditional means of potential outcomes. For each simulation, we conduct 1000 replications, and for each replication, we consider a sample size $n=1000$, with $Z_{i}\sim(2\mathcal{B}(2,4)-1)$ where $\mathcal{B}(p_1,p_2)$ is a beta distribution with parameters $p_1$ and $p_2$. We also consider $\varepsilon_{it}\sim N(0,\sigma^{2}_{\varepsilon})$, $\sigma_{\varepsilon}=0.1295$ and the outcome generated is $Y_{i,t}=\mu_{it}(Z_{i})+\varepsilon_{i,t}$. The detailed specifications and additional functional forms of both models are provided in Appendix (ref) for clarity. In both scenarios, we observe that our DiDC estimator performs better than the RDD estimator, with a smaller bias and improved coverage. Additionally, we find that the DiDC estimator has smaller bias and better coverage than the DiD estimator when the DGPs change over time.
In the first simulation, we mimic a scenario where one or more confounding discontinuities are present at the threshold $z_0=0$, but the functional forms for conditional means of potential outcomes are the same in both periods. This setting is comparable to that of grembi2016fiscal, where the treatment of interest was introduced at some point between $t=0$ and $t=1$ and was given to units whose running variable values $Z_i$ are above the threshold $z_0=0$ however there were other pre-existing treatments determined by the same threshold $z_0=0$ on the same running variable.
The second model we consider is a scenario with identical functional forms over time, the only distinction being that in period $t=1$, there is an effect of treatment $\tau$ for units with $Z_i \geq 0$. Importantly, there are no other sources of discontinuity in the outcome at the threshold $z_0=0$, making this an ideal scenario for estimating the effect of the treatment using a regular RDD.
We present the results of these simulation studies in Tables 1 and 2. We estimate the Diff-in-Disc along with the RDD, DiD and TWFE, and compare the average bias, median bias, root-mean-squared errors, 95% coverage probability, and the length of the 95% confidence interval for each estimator when the treatment effect $\tau$ is equal to 0.
Table (ref) shows a significant improvement over the RD estimator when a confounding factor is present at the threshold. The result is not surprising, as the presence of the confounder violates the continuity assumption that underlies the validity of cross-sectional RD. Note also that the mean size of the bandwidths, both for the estimator and for the bias-correction, is greater for the Diff-in-Disc estimator than the RD estimator. In terms of bias, both the nonparametric DiD and the TWFE estimators exhibit desirable finite-sample properties. However, when it comes to coverage of the confidence interval, we find that those from the TWFE estimator are severely biased.
Table (ref) shows the results for the case where there is no confounding policy at the threshold. In that case, the cross-sectional RD is valid, as evidenced by the small finite-sample bias in the simulation. Once again, the Diff-in-Disc estimator has smaller finite-sample bias and larger bandwidths, which illustrates the point in Section (ref) that, even in settings where the standard RD is valid, there might be gains in using the Diff-in-Disc approach.
The results show that Difference-in-Discontinuites estimator has an apparent improvement over the standard RD, yielding a smaller bias and better coverage. This is noteworthy as it highlights that even in cases where the RDD would traditionally be regarded as suitable, the differences-in-discontinuity approach can give more desirable results by incorporating more data into the estimation. The results hold when the functional format exhibits little temporal variation, as shown in the next section, with the additional period contributing to more reliable estimates than the RD whenever the bias from the extra period $t=0$ does not exceed that from period $t=1$. In Appendix (ref), we show that the results for the simulations are robust even in the case where the variance of potential outcomes vary over time.
Next, we introduce scenarios where the functional forms for conditional means change between periods, contrasting to the prior section where functional forms were identical across time. This setup allows us to analyze how changes in the conditional mean of potential outcomes affect the performance of the considered estimators.
Again, we consider simulations with and without confounders at the threshold, but now we alter the model for $t=1$ to be a linear model derived from the original. The functional form for time period $t=0$ is identical to that of models 1 and 2. Results are shown in Tables 3 and 4.
Table 3 shows that in the presence of changes in functional forms over time and confounding policies at the threshold, only the Diff-in-Disc approach is valid. The finite-sample bias of the standard RD estimator is close to the one presented in Table 1. The main difference in Table 3 in comparison to Table 1 is the poor finite-sample properties of DiD methods, as both the nonparametric DiD and the TWFE estimator exhibit larger finite-sample bias than the standard RD.
Table 4 displays the results for the case where functional forms change over time, but there is no confounding policy at the threshold. Once again, the standard RD and the Diff-in-Disc exhibit desirable finite-sample properties, with the Diff-in-Disc showing better coverage and smaller bias. However, unlike when functional forms are constant over time, when functional forms change, the standard RD method yields larger optimal bandwidths. For the DiD estimators, the results show severe bias, and the performance of the nonparametric DiD and the TWFE are similar to what we observed in Table 3. In Appendix D.2, we show that the simulation results are robust to changes in the time-varying variance.
Overall, the Monte Carlo results show several desirable properties of the Diff-in-Disc estimator. Not only it remains unbiased when the alternative approaches are not valid, but it also exhibits smaller finite-sample bias than the standard RD even in the absence of confounders. In the next section, we revisit a well-known political economy setting to analyze the estimator's performance on a real dataset.
We illustrate the use of our estimator by revisiting the empirical application in grembi2016fiscal, which analyzes the impact of fiscal rules on Italian municipal finances by exploiting a 2001 fiscal rule relaxation as a natural experiment. The relaxation applied to municipalities with fewer than 5,000 inhabitants, which therefore form the treatment group, while municipalities above this threshold serve as controls.
Our objective is to replicate their setting using the Difference-in-Discontinuities (DiDC) estimator, compare the results to the original findings, and assess the validity of the identifying assumption using the constant-confounder test from Section (ref).
The dataset comprises data from Italian municipalities, focusing on the period surrounding the government's relaxation of fiscal rules in 2001. Municipalities with fewer than 5,000 inhabitants experienced a relaxation of fiscal rules and were treated, while those with more than 5,000 inhabitants served as controls.
We implement a 2×2 DiDC design, computing before–and–after differences and estimating local linear RDDs as detailed in Section (ref). grembi2016fiscal's original study utilized a large panel dataset, a rectangular kernel and a polynomial of degree one as in the model below: {
} where $S_i$ is a dummy variable for cities below 5,000 (treatment indicator), $T_i$ is a dummy variable for the post-treatment period and $\beta_0$ is the parameter of interest.
Table 5 compares our DiDC estimates (difference of RDDs), along with the effective bandwidths and sample sizes to grembi2016fiscal's estimations. Due to differences in specifications, the results must be interpreted with caution.
The results are qualitatively similar, yet, our specification is underpowered due to the smaller sample size, and thus estimates are not statistically significant.
We conduct the test from Section (ref) to evaluate the assumption of time-invariant confounding effects. To test this, we estimate stacked RDs for the years 1998--2000, interacting the running variable with treatment and period dummies. Table (ref) reports the interaction coefficients ($z\times D$) for the smallest bandwidth, and Table (ref) summarizes the corresponding joint Wald tests for equality of these coefficients across years.
The joint test of $H_0:\theta_{2000}=\theta_{1999}=\theta_{1998}$ is implemented using a Wald $F$-test based on the pooled regressions. The results are shown in Table (ref).
For taxes, the null hypothesis of time-invariant confounding is rejected at conventional significance levels (p $\in$ [0.01, 0.05]), indicating that the RD discontinuity varied across pre-treatment years. In contrast, the test fails to reject $H_0$ for both deficit and fiscal balance (p-values between 0.75 and 0.87).
In summary, our estimated treatment effects are similar in magnitude to the difference-in-discontinuities approach utilized in grembi2016fiscal study. Notably, our difference-in-RDDs approach yields smaller confidence intervals, suggesting it is the most powerful. However, our validity test indicates that confounding effects vary over time, suggesting that the use of DiDC in this setting leads to biased estimates. We also replicate our estimates using the bandwidths in ludwig2007does, with similar results. These estimates are shown in Table (ref) in Appendix (ref).
We also implement the KS-type test from Section (ref) to assess whether the shapes of the conditional mean functions remain stable across pre-treatment years. Applying this test to each outcome, we find strong evidence of time-varying functional forms for taxes and fiscal balance: on both sides of the cutoff, the bootstrap p-values for imposte and saldo are essentially zero (0.001 on both sides), indicating clear violations of the time-invariance condition. For the deficit, however, the evidence is mixed. The left-hand side yields a p-value of 0.266, consistent with time-invariant functional forms, whereas the right-hand side again produces a very small p-value (0.001), suggesting that the relationship between the running variable and the deficit outcome changed over time for municipalities above the 5,000-inhabitant threshold. Overall, these results show that the shape of the outcome–running-variable relationship is not stable across pre-treatment years for most outcomes, reinforcing the need for careful consideration when applying the method.
Table (ref) reports the identified sets for the treatment effect on Taxes under the bounded-variation sensitivity framework. Each cell shows the lower and upper bound for the identified set $\left[\tau_{c}^{LB},\,\tau_{c}^{UB}\right]$ as a function of the sensitivity parameters $c_1$ and $c_2$. The bounds are constructed by allowing the conditional potential-outcome functions and the confounding effect to deviate from the time-invariance baseline up to $c_1$ and $c_2$, respectively (Section (ref)). The same grid procedure was used to produce the analogous tables for Fiscal balance (Tables (ref)--(ref) in Appendix (ref)) and Deficit (Tables 12 and 13 from Section (ref)).
Under Assumptions (ref) and (ref), we construct worst-case modularity bounds using $y_{min}$ and $y_{max}$ equal to the minimum and maximum observed values of each outcome in the data (Table (ref)). As expected in this worst-case setting, the resulting identified sets are very wide: for example, the bounds for taxes span roughly $[-157.59,\; 1203.19]$ for $\tau_{c}$ and $[-1203.19,\; 108.08]$ for $\tau_{uc}$, with similarly large intervals for the fiscal gap and deficit. Consequently, these bounds don't provide much information. To obtain informative conclusions about the direction or magnitude of the treatment effect, researchers must impose stronger and more economically grounded restrictions.
When the two approaches are combined, identification becomes substantially sharper. Together, these restrictions eliminate many of the extreme scenarios allowed by modularity alone and narrow the identified sets relative to either assumption in isolation.
The Deficit outcome clearly illustrates the value of this combination. Under worst-case modularity, the identified sets for Deficit include both large negative and positive effects, but once bounded-variation constraints are imposed, the intervals become notably narrower and remain strictly negative across all $c_1$ and $c_2$ values considered. Table (ref) shows that the bounds for $\tau_c$ fall roughly between $-27$ and $-19$, and Table (ref) similarly reports negative upper endpoints for $\tau_{uc}$. The combined assumptions, therefore, identify not only the sign of the effect but also restrict its magnitude for $\tau_c$ to a reasonably small range.
The results in Table (ref) show that the bounds on the effects of the fiscal rules are overall uninformative, suggesting that the results are not robust to violations of the time-invariance assumptions. Tables 12 and 13, on the other hand, show that the bounds on the effects on deficit are robust both to violations on time-invariance assumptions and complementarity between treatment and confounding effects. The results in the table show that the evolution of the mean confounding effect and the mean confounded deficit at the threshold could exceed 6 euros per capita (roughly 30.000 euros given the 5.000 population rule), and the true effect of relaxing fiscal policy would still be positive. However, the setting is severely underpowered, which means that confidence intervals are uninformative regarding the true signal of the treatment effects\footnote{In Tables (ref) and (ref), which show informative identified sets, the lower bound of the confidence interval for the lower bound is always negative.}.
Finally, Table (ref) shows that modularity assumptions alone are not sufficient to yield informative identified sets for the treatment effects.
The difference-in-discontinuities (DiDC) design is emerging as a promising method for estimating causal inference, addressing the limitations of both regression discontinuity (RDD) and difference-in-difference (DiD) approaches. This paper lays the theoretical groundwork for DiDC, examines its identification assumptions, estimation procedures, and asymptotic properties. We showcase its advantages through Monte Carlo simulations and an empirical application.
DiDC can handle scenarios in which the control and treatment groups differ significantly, violating the parallel trends assumption of DiD, or when RDD encounters confounding factors at the threshold. By incorporating more information, DiDC eliminates bias in RDD estimates under specific assumptions about the data-generating processes.
However, it is important to give due attention to the identification assumptions, particularly the time-invariance of confounding effects. We propose a test based on stacked RDDs to assess its validity in practice. Additionally, DiDC requires the treatment effect to be independent of confounding policy, though we introduce a possible relaxation for potential interaction effects.
We find that the DiDC method can eliminate bias completely if the function shapes remain stable on both sides of the threshold over time. This suggests it could offer significant advantages over standard RDD estimators, even in settings where no other confounding variables are present at the threshold. We also propose a test to compare the derivatives of estimated functions on either side of the threshold, allowing researchers to evaluate the stability of data-generating processes over time.
Monte Carlo simulations demonstrate DiDC's potential to improve upon RDD, yielding lower bias and better coverage. It can provide more desirable results by incorporating more data, especially when the functional form exhibits minimal temporal variation. Notably, DiDC is the only viable approach when confounding factors render both RDD and DiD unsuitable. The empirical application highlights the importance of the time-invariance assumption.
Future research directions include developing robust alternative estimators that are robust to violations of identification assumptions, as well as exploring other confidence interval methods tailored to the DiDC design. Overall, the DiDC method offers a valuable addition to the causal inference toolkit. It is applicable in settings where no other methods were previously available and shows potential to reduce bias in estimation in other settings.
\nocite{imbens2008regression} \nocite{rdrobust}