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.
97,464 characters · 25 sections · 61 citation commands
Correcting invalid regression discontinuity designs with multiple time period data
Regression Discontinuity (RD) is a widely used design for causal inference that exploits an observable discontinuity in treatment assignment at a cutoff imbens2008regression, lee2010regression. A common identification assumption in RD is that the mean potential outcomes are continuous at the cutoff hahn2001identification. However, this assumption may fail when other policies or treatments are also implemented at the same cutoff. In such cases, researchers often incorporate identification assumptions similar to difference-in-differences (DID) designs, but without sufficient theoretical justification. This lack of formal grounding may lead to incorrect empirical conclusions. In this paper, we develop a general identification and estimation framework for such settings, which we term RD-DID.
To illustrate these challenges, consider the following two applications where the continuity assumption might fail. bertrand2021improving studied an educational reform in Norway, where children born after January 1, 1978, entered a reformed schooling system. However, birth timing before and after January 1 also affects school cohort placement, creating an additional discontinuity at the cutoff. Similarly, bennedsen2022firms examined Danish pay transparency laws, where firms above a 35-employee threshold faced new reporting requirements. A potential concern is that other regulations also take effect at the 35-employee threshold. In both of these examples, there is a concern the continuity assumption does not hold and the authors utilized time periods prior to the reform, where no units were treated, to correct the RD design.
We first analyze sharp RD designs when the continuity assumption does not hold, and show that the observed outcome discontinuity at the cutoff can be decomposed into a causal effect plus a bias term. The causal effect is the average treatment effect on the treated (ATT) or untreated (ATU) local to the cutoff, while the bias term is the discontinuity in treated or untreated potential outcomes at the cutoff. In fuzzy RD, we find that the bias terms are a weighted average of the treated and untreated potential outcome discontinuities.
To identify the local ATT and ATU when the the continuity assumption is violated, we develop a general identification framework for sharp RD treatment assignment using multiple time periods. Specifically, we incorporate time periods in which all units are treated or all units are untreated. In time periods where all units are treated (untreated) the discontinuity in treated (untreated) potential outcome means is identified. Assuming how these discontinuities change over time--e.g., assuming they remain constant or follow a linear trend--we can recover the bias in time periods with RD treatment assignment, enabling identification of the local ATT and ATU.
Building on the general RD-DID framework, we study three extensions to the identification framework which are common in applications. First, we analyze RD-DID with fuzzy RD treatment assignment. Unlike single-period fuzzy RD, which requires additional assumptions such as local conditional independence or monotonicity (cattaneo2022regression), we show that fuzzy RD-DID also demands more data: observations from both time periods where all units are treated and where all units are untreated. Second, we show that carry-over effects complicate identification in RD-DID, especially if periods where all units share the same treatment status occur after RD assignment periods. We leverage an assumption similar to no-anticipation in DID callaway2021difference, to obtain identification. Lastly, we consider settings where the running variable changes over time. We show this introduces composition effects. While removing treatment switchers might mitigate composition effects, this practice alters the population near the cutoff, creating an analysis similar to ”donut-hole” RD (bajari2011regression).
A common alternative in applied research when under sharp RD assignment uses DID methods (e.g., schonberg2014expansions, lalive2014parental, danzer2018paid). We highlight two differences between RD-DID and this alternative.\footnote{The comparison between RD-DID and DID with an RD treatment assignment relates to literature on the justification and testing of the parallel trends assumption (e.g., ghanem2022selection,rambachan2023more,roth2023parallel).} First, DID identifies a global causal effect, whereas RD-DID target local effects. Moreover, the parallel trends assumption in DID is much stronger than the local assumptions required for RD-DID, making it less plausible in many RD settings. Therefore, we argue that researchers should not apply DID methods in RD contexts without careful justification.
Building on the contemporary single-period RD estimation framework (for a recent review see cattaneo2022regression), we develop an estimation approach suited for multiple periods in RD-DID. Our approach applies local linear regression estimators on both sides of the cutoff at each time period, and then aggregates the estimated outcome discontinuities according to their assumed trend over time. We implement this approach for both conventional and bias-corrected estimators (calonico2014robust). A R package \href{https://github.com/dorleventer/rddid}{rddid} that implements the proposed estimation framework is provided for general use.
We derive variance estimators for different sampling schemes. We consider repeated cross-section (CS), where different units are observed in each period, and panel data, where the same units are tracked over time. For panel data, we distinguish between two cases: time-constant running variable (PC), and time-varying running variable (PV). We show that the variance differs across sampling schemes due to the form of the covariance of outcome discontinuity estimators across periods. Through Monte Carlo simulations, we evaluate the finite-sample performance of our variance estimators. Coverage rates align with known results in single-period RD, with bias-corrected confidence intervals outperforming conventional ones. Finally, we show that correctly specifying the sampling scheme is crucial--assuming CS sampling while the data is PC leads to overestimated standard errors.
To illustrate our results, we revisit grembi2016fiscal, which studies the impact of fiscal rules on financial outcomes in Italy. The 2001 reform imposed deficit targets on municipalities above 5,000 residents, creating a sharp RD design. However, other regulations also change discontinuously at the same cutoff, calling for an RD-DID design. We test the assumption of time-constant potential outcome discontinuities using equivalence tests (hartman2018equivalence), and find the data reject this identification assumption. Estimating the ATU under a linear-in-time discontinuity instead, we find larger treatment effects in absolute terms compared to estimates under constant discontinuities, suggesting underestimation of the effects fiscal rule in the original study. Finally, we find that 11% of municipalities switch treatment status, potentially introducing composition effects. Direct estimates confirm the composition effects are non-negligible. However, removing switchers distorts the running variable density, creating a drop in the density at the cutoff.
This paper contributes to the literature on identification in confounded RD designs (mealli2012evaluating,wing2013strengthening,grembi2016fiscal,eggers2018regression,galindo2018fuzzy,picchetti2024difference). grembi2016fiscal highlight the role of additional time periods in RD settings, while our approach generalizes RD-DID to accommodate time-varying discontinuities, moving beyond the standard two-period differencing approach. We further extend the framework to fuzzy RD-DID, carry-over effects, and time-varying running variables, expanding the range of applications where RD-DID can be applied. A related issue arises in studies that implement RD-DID using differenced outcomes (e.g., bagues2021can), as formalized in picchetti2024difference. However, when the running variable changes over time, such regressions lack a clear interpretation. Our framework remains well-defined under both constant and time-varying running variables, ensuring broader applicability and clearer causal interpretation.
Beyond identification, this paper also contributes to the literature on estimation in RD designs (see, e.g., imbens2012optimal,calonico2014robust,gelman2019high). We extend the contemporary local linear framework to multiple time periods, providing a systematic approach for estimation and inference in such designs. These contributions are particularly relevant for empirical studies that estimate outcome discontinuities using pooled regressions across multiple periods, a common practice in applied work (e.g., bertrand2021improving,baltrunaite2019let,bennedsen2022firms,avdic2018modern).\footnote{By pooled regression, we refer to a linear regression that combines multiple time periods while allowing intercepts and slopes to vary across periods.} Unlike pooled regressions, which implicitly impose strong assumptions on the evolution of potential outcome discontinuities, our approach allows for greater flexibility in modeling how discontinuities evolve over time. Additionally, our framework facilitates bias correction and variance estimation, while clarifying the role of fully treated and untreated periods, which pooled approaches often implicitly combine rather than use to identify distinct parameters.
The paper is structured as follows. Section (ref) introduces the setup and notation. Section (ref) examines bias in single-period RD under continuity violations. Section (ref) develops the RD-DID framework and extends it to fuzzy RD, carry-over effects, and time-varying running variables. Section (ref) compares RD-DID to DID with RD treatment assignment. Section (ref) discusses estimation and inference, while Section (ref) evaluates finite-sample performance via simulations. Section (ref) illustrates the results in an application. Section (ref) concludes. Proofs of all theorems are presented in Appendix (ref). A R package \href{https://github.com/dorleventer/rddid}{rddid} implements the estimation framework.
Let $\mathcal{T}=\{1,..,T\}$ be a finite set of $T$ time periods. For each time period $t\in\mathcal{T}$, data on $n_t$ units is collected. Let $R_{it}$ and $W_{it}$ be the running variable and treatment value for unit $i$ at time $t$, respectively. We assume all units are $iid$, so we occasionally omit the index $i$ to improve clarity. We denote by $c$ a time-invariant cutoff.
Let $p_t(r)=\Pr\left(W_t=1\mid R_t=r\right)$ be the probability of receiving the treatment conditionally on the value of the running variable. We denote limits of the conditional treatment probabilities above and below the cutoff $c$ by $p_{t,(+)}=\lim\limits_{r\rightarrow c^{+}}p_{t}\left(r\right)$ and $p_{t,(-)}=\lim\limits_{r\rightarrow c^{-}}p_{t}\left(r\right)$, respectively. In any time period $t$ when treatment is allocated according to an RD design, $p_{t,(+)}\neq p_{t,(-)}$. In a “sharp RD”, one of the limits equals one and the other zero, and WLOG in such designs $p_{t,(+)}=1$ and $p_{t,(-)}=0$. All other RD designs are termed “fuzzy RD”.
Using the potential outcomes framework, let $Y_{i,t}(0)$ and $Y_{i,t}(1)$ be the potential outcomes under non-treatment and treatment, respectively, for unit $i$ at time $t$. We take the standard assumptions that the observed and potential outcomes are connected via the realized treatment at time $t$, i.e., $Y_{i,t}=Y_{i,t}(W_{i,t})$. This setup assumes no carry-over effects of the treatment, an assumption we relax in Section (ref). Let $\mu_{w,t}(r)= \mathbb{E}[Y_{t}(w)|R_{t}=r]$ and $\mu_t(r)=\mathbb{E}[Y_{t}|R_{t}=r]$ be the mean potential and observed outcomes, respectively, at $R_{i,t}=r$ and time $t$. Similar to the conditional probability of the treatment, we use plus and minus subscripts to denote limits of the conditional expectations of the outcome above and below $c$, respectively. That is, $\mu_{w,t,(+)}=\lim\limits_{r\rightarrow c^{+}}\mu_{w,t}\left(r\right)$ and $\mu_{w,t,(-)}=\lim\limits_{r\rightarrow c^{-}}\mu_{w,t}\left(r\right)$ for the mean potential outcomes, and, $\mu_{t,(+)}=\lim\limits_{r\rightarrow c^{+}}\mu_t\left(r\right)$ and $\mu_{t,(-)}=\lim\limits_{r\rightarrow c^{-}}\mu_t\left(r\right)$ for the observed outcome mean.
A leading approach for identification in a single time period RD designs relies on the assumption of continuity of the conditional expectations of the potential outcomes at the cutoff (hahn2001identification).\footnote{An alternative identification framework is local randomization (see cattaneo2022regression).} In this section, we examine the bias arising when the continuity assumption, defined formally below, does not hold in a classical single-period RD.\footnote{The bias discussed in this section is due to identification. It is inherently different from the widely discussed estimation bias analyzed in the RD literature.} We return to the multiple-period setup in Section (ref), where we present remedies for this bias. We therefore omit in this section the subscript $t$ so, for example, $\mu_{w}(r)=\mu_{w,t}(r)$, and similarly $\mu_{w,(+)}, \mu_{w,(-)}, \mu_{(+)}$ and $\mu_{(-)}$ all represent the respective quantities at a specific time period with an RD assignment.
In RD designs, researchers center their efforts around the average treatment effect (ATE) local to the cutoff, $ATE(c)= \mu_{1}(c) - \mu_{0}(c)$. As discussed above, a leading approach for identification relies on the following continuity assumption.
Under Assumption (ref), the $ATE(c)$ is identified in a sharp RD design by $ATE(c)=\mu_{(+)}-\mu_{(-)}$. To provide intuition about violations\footnote{It is important to distinguish between a discontinuity in potential outcomes at the cutoff due to a confounded treatment or policy, and a discontinuity due to manipulation of the running variable (mccrary2008manipulation,lee2010regression). In the analysis that follows, we assume away the case of manipulation. We return to this point in Section (ref) which considers the case of time-varying running variables.} of Assumption (ref), we visualize in Figure (ref) a stylized example. As a benchmark, the dashed dark-blue line represents a case where $\mu_1(r)$ is continuous at $r=c$, and hence Assumption (ref) is not violated. If, however, the mean treated potential outcomes is as described by one of the dotted light-blue lines, then there is a discontinuity at the cutoff, and Assumption (ref)(i) does not hold. A similar discussion can be made with respect to the mean potential outcomes under no treatment, by comparing the dashed and dotted orange lines.
To formalize the impact of wrongfully assuming continuity of the mean potential outcomes at the cutoff, we consider the following assumption, which is considerably weaker than Assumption (ref).
This assumption is mainly technical, and is expected to hold whenever the potential outcomes and the functions $\mu_{w}(r)$ are well-defined around the cutoff. Under Assumption (ref), the difference between the limits of the mean potential outcomes under no treatment and treatment are well-defined for each side of the cutoff. We denote these differences by $\alpha_{0} =\mu_{0,(+)} - \mu_{0,(-)}$ and $\alpha_{1} =\mu_{1,(+)} - \mu_{1,(-)}$. Assumption (ref) can be viewed as a special case of Assumption (ref), with $\alpha_0=\alpha_1=0$.
We can now define two new causal estimands of interest which are well-defined under Assumption (ref). Let $ATT(c)=E[Y(1)-Y(0)|W=1, R=c]$ be the average treatment effect on the treated (ATT) local to the cutoff and $ATU(c)=E[Y(1)-Y(0)|W=0, R=c]$ be the average treatment effect on the untreated (ATU) local to the cutoff. Under Assumption (ref), these causal estimands can be written as $ATT(c) = \mu_{1,(+)} - \mu_{0,(+)}$ and $ATU(c) = \mu_{1,(-)} - \mu_{0,(-)}$.
The following theorem provides a characterization of the bias in the single time-period sharp RD design under the more general Assumption (ref). Let $D^Y = \mu_{(+)}-\mu_{(-)}$ be the discontinuity of the mean observed outcome at the cutoff. Under Assumption (ref), $D^Y$ can be decomposed into a causal effect and a bias term.
Note that under the assumption that $\alpha_0=0$, the outcome discontinuity $D^Y$ can be interpreted as the $ATT(c)$. An analogous statement about the $ATU(c)$ can be made when researchers believe $\alpha_1=0$. When $\alpha_0=\alpha_1=0$, and hence Assumption (ref) holds, $ATE(c)=ATT(c)=ATU(c)$.\footnote{The case where only one of the $\alpha$'s is zero and the other is not, is technically possible, but presumably quite rare. For example, we will have that $\alpha_0=0$ and $\alpha_1 \ne0$ if at the cutoff, a second policy is implemented (in addition to $W$), however this policy only affects the outcome if $W=1$.}
Returning to Figure (ref), the letters depict a hypothetical example for the true limit points of the potential outcome conditional means at the cutoff $c=50$. In this example, the $ATT(c)$ is equal to the difference in the y-coordinates between point B (representing $\mu_{1,(+)}$) and point C ($\mu_{0,(+)}$), while the $ATU(c)$ is the difference between point A ($\mu_{1,(-)}$) and point D ($\mu_{0,(-)}$). Lastly, the observed discontinuity $D^Y$ is the difference between points B and D. In this example, the observed discontinuity is larger than the $ATT(c)$ and smaller than the $ATU(c)$. In terms of the result in Theorem (ref), these differences between $D^Y$ and $ATT(c)$ or $ATU(c)$ are due to $\alpha_1<0$, since A is above B, and $\alpha_0>0$, since C is above D.
We now characterize the bias in fuzzy RD designs. Let $D^W=p_{(+)} - p_{(-)}$ be the discontinuity in the observed treatment probability at the cutoff. In fuzzy RD, the $ATE(c)$ is identified by $\frac{D^Y}{D^W}$ under continuity and the additional conditional independence assumption (CIA) between treatment and potential outcomes at the cutoff (hahn2001identification).\footnote{Because a local conditional independence assumption might be too restrictive in practice, alternative assumptions have been considered, such as monotonicity, which lead to identification of alternative causal estimands, such as the local average treatment effect of compliers (see cattaneo2022regression and references therein).} For clarity, we present here the characterization under CIA, defined formally as $Y(0),Y(1)\perp \!\!\! \perp W\mid R=r$, for $r$ near $c$. The following theorem shows that whenever Assumption (ref) does not hold, but Assumption (ref) does, then $\frac{D^Y}{D^W}$ equals to a causal effect plus a bias term.
A main difference between sharp and fuzzy RD designs when the continuity assumption is violated, is that in the fuzzy setting the bias terms for the $ATT(c)$ and $ATU(c)$ each depends on both $\alpha_0$ and $\alpha_1$. Therefore, even if researchers believe, for example, that $\alpha_0=0$, neither the $ATT(c)$ nor the $ATU(c)$ are identifiable from the data without the additional assumption that $\alpha_1=0$. This property has consequences on identification in the multiple time periods setup, as we show in Section (ref).
In the previous section, we derived the bias of the RD estimator when the continuity assumption is violated. The results of Theorems (ref) and (ref) direct us to target the $ATT(c)$ or $ATU(c)$, by the means of identifying of $\alpha_0$ and $\alpha_1$. To this end, we present in this section a general identification framework that utilizes data from multiple periods. Similar to the difference-in-differences (DID) design, we formulate assumptions on how mean potential outcomes are related across these periods. Therefore, we term this identification framework RD-DID. We further refer to the canonical RD-DID design as the design with two time periods, where all units are untreated in the first time period and a sharp RD treatment assignment takes place in the second time period. This creates two cohorts: a never-treated cohort (not treated in the RD), and a treated cohort (treated in the second period RD). In terms of the notations introduced in Section (ref), in a canonical RD-DID design $T=2$, $\mathcal{T}_0=\left\{1\right\}$, $\mathcal{T}_1=\emptyset$, and $\mathcal{T}_{\mathrm{RD}}=\{2\}$.
Bringing back the $t$ notation, and following Section (ref), we define the following period-specific causal estimands. Let $ATE(c,t) = \mu_{1,t,(+)} - \mu_{0,t,(-)}$, $ATT(c,t) = \mu_{1,t,(+)} - \mu_{0,t,(+)}$ and $ATU(c,t) = \mu_{1,t,(-)} - \mu_{0,t,(-)}$ denote the ATE, ATT and ATU local to the cutoff $c$ at time period $t$, respectively.
We begin by generalizing Assumption (ref) to the setting of multiple periods.
Differently from Assumption (ref), Assumption (ref) also assumes the existence of the limits for the periods in either $\mathcal{T}_0$ or $\mathcal{T}_1$. Nevertheless, it is, like Assumption (ref), a technical and rather weak assumption. Whenever Assumption (ref)(i) holds, let also $\alpha_{1,t}=\mu_{1,t,(+)}-\mu_{1,t,(-)}$ be the period-specific discontinuity of $\mu_{1,t}$ at the cutoff $c$. Similarly, define $\alpha_{0,t}=\mu_{0,t,(+)}-\mu_{0,t,(-)}$ at each $t\in\mathcal{T}$ Assumption (ref)(ii) holds for.
Let $t^\star\in \mathcal{T}_{\mathrm{RD}}$ be a time period of interest. Assume the RD assignment follows a sharp RD design, and that Assumption (ref) holds. According to Theorem (ref), the observed outcome discontinuity in period $t^\star$ is a biased estimator of $ATT(c,t^\star)$ (or of $ATU(c,t^\star)$) with a bias term $\alpha_{0,t^\star}$ ($\alpha_{1,t^\star}$). If data from other periods can be used to identify and estimate $\alpha_{0,t^\star}$ and $\alpha_{1,t^\star}$, those parameters can be used to identify $ATU(c,t^*)$ and $ATT(c,t^*)$. The following assumption formalizes the general assertion that $\alpha_{0,t^\star}$ and $\alpha_{1,t^\star}$ can be learned from data obtained at other periods.
Assumption (ref) asserts that $\alpha_{0,t^\star}$ and $\alpha_{1,t^\star}$ can be written as functions of similar potential outcome discontinuities in periods where no unit is treated or all units are treated, respectively. We give two specific examples of such functions below.
We are now ready to present our main identification result. Let $D^Y_t=\mu_{t,(+)}-\mu_{t,(-)}$ denote the discontinuity in the mean observed outcome at the cutoff at time $t$.
The intuition behind Theorem (ref)(i) is as follows. In time periods $t\in\mathcal{T}_0$, observed outcomes from both sides of the cutoff are equal to the untreated potential outcomes, and hence $\alpha_{0,t}$ is identifiable. Under Assumption (ref), the bias parameter $\alpha_{0,t^\star}$ is obtained by applying $g_0$ to the $\alpha_{0,t}$'s in time periods $t\in\mathcal{T}_0$. By Theorem (ref), the $ATT(c,t^\star)$ is then identified as the sum of $D^Y_{t^\star}$ and the bias parameter $\alpha_{0,t^\star}$. The intuition behind Theorem (ref)(ii) is similar.
Although without the formal framework presented here, it is common in practice to assume the discontinuities $\alpha_{0,t}$ or $\alpha_{1,t}$ are constant across time. The following corollary underpins this approach, and extends it to linear-in-time trends for $\alpha_{0,t}$ or $\alpha_{1,t}$.
The results for constant $\alpha_{w,t}$ hold for any set of weights $\pi_t$. However, as we discuss in Section (ref), for estimation certain weighting schemes may be preferable. Nevertheless in certain contexts, the constant $\alpha_{w,t}$ assumption might be plausible only for a subset of periods in $\mathcal{T}_0$ or $\mathcal{T}_1$. In these situations, one can simply use these subsets of time periods instead of the whole set of time periods. Note also that Corollary (ref) states that identification in the canonical RD-DID is equivalent to a difference between two outcome discontinuities (grembi2016fiscal).
In some applications, the constant $\alpha_{w,t}$ assumption cannot be justified, and researchers might opt for a linear trend across time in $\alpha_{0,t}$ or $\alpha_{1,t}$, as formally presented Corollary (ref)(ii). The parameters $m_0$ or $m_1$ can be identified using any two or more time periods in $t\in\mathcal{T}_0$ or $t\in\mathcal{T}_1$, respectively.
For a specific dataset, one can test whether the chosen functional forms of $g_0$ and/or $g_1$ are supported by the data. In Section (ref), we consider an equivalence testing procedure for the constant discontinuities assumption in the application, and show the data does not support this assumption for the considered outcomes. We hence opt for a linear-in-time discontinuities assumption for that outcome. An alternative approach to pre-specifying $g_0$ and/or $g_1$, which may be useful in applications with a large number of time periods, is to avoid strong assumptions on these functions and learn them from the data in a non- or semi-parametric manner.
Similar to the discussion in callaway2021difference, we can consider aggregation of $ATT(c,t)$s or $ATU(c,t)$s. Due to the time-invariance of the cutoff there is only one treatment cohort. Hence the aggregations by calendar time and length of exposure are equivalent. Also, the aggregation to an overall treatment effect is equivalent whether by group, time or length of exposure.
The setup in Theorem (ref) is quite common in practice. Nevertheless, it is overly simplistic with respect to certain applications. The following subsections consider more complex setups, including fuzzy RD designs, carry-over effects, and time-varying running variable.
As with sharp RD, the identification framework for the causal estimands $ATT(c,t^\star)$ and $ATU(c,t^\star)$ for $t^\star\in\mathcal{T}_{\mathrm{RD}}$ under a fuzzy RD treatment design is motivated by Theorem (ref). Let $D^W_t = p_{t,(+)} - p_{t,(-)}$ denote the discontinuity in treatment probability at the cutoff at time $t$.
The proof is similar to the proof of Theorem (ref). If Assumption (ref) holds, then $\alpha_{0,t_0}$ and $\alpha_{1,t_1}$ are identified in time periods $t_0\in\mathcal{T}_0$ and $t_1\in\mathcal{T}_1$, respectively. If Assumption (ref) also holds then $\alpha_{0,t^\star}$ and $\alpha_{1,t^\star}$ are identified through the functions $g_0$ and $g_1$. From Theorem (ref), the $ATU(c,t^\star)$ and the $ATT(c,t^\star)$ are hence identified.
Theorem (ref) reveals that unlike identification in the sharp RD-DID design (Theorem (ref)), identification under the fuzzy RD-DID design requires both sets $\mathcal{T}_0$ and $\mathcal{T}_1$ to be non-empty, since both $\alpha_{0,t^\star}$ and $\alpha_{1,t^\star}$ are needed to identify the causal estimands of interest. Put differently, when a fuzzy RD suffers from a violation of continuity (with respect to potential outcomes), both time periods when all units are treated and time periods when all units are untreated are needed for identification.
An important implication of the above result is that in the two-period canonical RD-DID, the target causal estimands are not identifiable. Because in a canonical RD-DID $\mathcal{T}_1$ is empty, it is impossible to identify the bias parameter $\alpha_{1,2}$. Thus, stronger assumptions are needed for identification. For example, if one is willing to assume that $\alpha_{0,t}=\alpha_{1,t}=\alpha$ for $t=1,2$, i.e., that the discontinuities for untreated and treated potential outcomes are equal and constant across time, then the single time period $t=1$ suffices for identification, and $ATT(c,2)=ATU(c,2)$.
Carry-over effects arise when potential outcomes depend on treatment status in prior time periods. To accommodate carry-over effects, potential outcomes are defined as functions of the treatment path, i.e., of the vector of treatment statuses at each time period. In a setup with a time-invariant running variable and cutoff, only two treatment paths are possible, since treatment is the same for all units at time periods $t\in\mathcal{T}_0$ or $t\in \mathcal{T}_1$, and is different between units at time periods $t\in\mathcal{T}_{\mathrm{RD}}$. Consequentially, the causal estimands and bias parameters are defined using these two treatment paths. The following example illustrates that carry-over effects introduce a problem for identification in RD-DID designs.
The example illustrates that in the presence of carry-over effects, identification of the mean potential outcome discontinuities $\alpha_{w,t}$ requires additional assumptions. We now present such an assumption. Let $tp_1$ ($tp_0$) denote the treatment path where a unit is treated (untreated) in $\mathcal{T}_{\mathrm{RD}}$. In the above example, $tp_1= (0,1,1)$ and $tp_0=(0,0,1)$. Let $Y_{i,t}(tp)$ denote the potential outcome of unit $i$ at time $t$ under treatment path $tp$, and let $\mu_{tp,t}(r)=\mathbb{E}[Y_t(tp)\mid R=r]$. Similar to before, let $\mu_{tp,t,(+)}$ and $\mu_{tp,t,(-)}$ represent the limits above and below the cutoff, respectively. The following assumption states that the limits of the potential outcome means are equivalent under both treatment paths, below the cutoff for time periods $t\in\mathcal{T}_1$ and above the cutoff for time periods $t\in\mathcal{T}_0$.
Assumption (ref)(ii) states that local to the cutoff $c$, the potential outcome means under treatment and non-treatment are equivalent for units above the cutoff, i.e., the treatment group in RD, in time periods where no unit is treated. Assumption (ref)(i) makes a similar claim for units below the cutoff, i.e., the untreated group in RD, in time periods where all units are treated. An alternative interpretation of Assumption (ref)(i) is that $ATU(c,t)=0$ for $t\in\mathcal{T}_1$.\footnote{Allowing carry-over effects, the ATU and ATT are defined as $ATU(c,t)=\mu_{tp_{1},t,(-)}-\mu_{tp_{0},t,(-)}$ and $ATT(c,t)=\mu_{tp_{1},t,(+)}-\mu_{tp_{0},t,(+)}$. } Similarly, Assumption (ref)(ii) can be interpreted as assuming $ATT(c,t)=0$ for $t\in\mathcal{T}_0$. Finally, Assumption (ref)(ii) is a local version of the no-anticipation assumption in the DID identification framework (roth2023s). In DID designs, researchers are usually willing to assume Assumption (ref)(ii) for time periods before treatment takes place. Whether the researcher can make such claims on time periods after treatment takes place depends on the context. Returning to the above example, $ATU(c,2)$ is not identifiable if Assumption (ref)(i) does not hold in $t=3$. However, since in $t=1$ all units are untreated, and assuming differences in future treatment status do not affect current potential outcomes, then it is plausible that Assumption (ref)(ii) holds. This implies that $D_{1}^{Y}=\alpha_{0,1}$, and hence $ATT(c,2)$ is identifiable under Assumption (ref).
For completeness, Appendix (ref) presents a version of the bias characterization theorem and identification theorem when carry-over effects are present. The results are essentially unchanged, except that the studied causal estimands are those defined in this section, and that the additional Assumption (ref) is required for identification.
We now extend the setup to scenarios where the running variable $R_{i,t}$ may change over time within units. Consider the setup with no carry-over effects. If $R_{i,t}$ may change over time, then units can cross above and below the cutoff in different time periods. These changes in the treatment allocation resulting from changes in the running variable can introduce composition effects. Furthermore, such changes and effects may be the result of manipulation of the running variable. We first consider composition effects in general, and then discuss testing procedures for manipulation using multiple time periods.
The following example illustrates how time-varying running variables can introduce composition effects. By composition effects, we mean changes in expectations of potential outcomes due to changes in the composition of units that the expectation is taken over.
The above example illustrates that, in a setup with a time-varying running variable, assumptions placed on how $\alpha_{w,t}$ changes (or does not change) over time, i.e., specifying the functions $g_0$ or $g_1$ in Assumption (ref), entails an implicit strong assumption on the effect of a composition change on potential outcomes, namely on $\Delta R_{s,t} $. This suggests, that in applications with time-varying running variables, researchers should test for composition effects. Using the notation of the above example, this amounts to testing whether $\Delta R_{s,t} $ is equal to zero. This is a viable approach, since in time periods where no units is treated $\Delta_{s,t} R$ is identifiable from the data. Such a testing procedure is shown in the application in Section (ref).
As indicated above, assuming that $\Delta Y_{s,t} = 0$ immediately implicates that $\alpha_{w,t}$ is constant, a viable identification assumption when the running variable does not change over time (Corollary (ref)). If, however, $\Delta R_{s,t} \ne 0$, the composition around the cutoff varies across periods. In that case, assuming $\Delta Y_{s,t} = 0$ effectively means that any change in discontinuities between periods is driven by the groups composition above and below the cutoff, and not by changes in the outcome due to the treatment. Therefore, assuming $\alpha_{w,s}-\alpha_{w,t}=\Delta R_{s,t}$ is a more subtle assumption, because we are allowing composition effects but ruling out shifts in how potential outcomes evolve for any given composition.
Due to these potential composition effects, researchers are faced with a trade-off when the running variable may change over time. Assuming away these composition effects might be an implausible assumption. A possible alternative is to omit units which cross the cutoff $c$ between time periods due to a change in their running variable. Such units are often termed treatment switchers. If the probability of crossing the cutoff is higher for units that are initially close to the cutoff, omitting treatment switchers creates a drop in the density of the running variable at the cutoff, hence mechanically creating a situation similar to “donut-hole” RD (bajari2011regression,noack2023donut). In RD-DID we expect this to be the common case, since minor changes in the running variable may cause units that are close to the cutoff to cross the cutoff. Such a drop can create problems in both the interpretation and estimation of the estimand. Removing units close to the cutoff affects the composition of units near the cutoff and thus the estimand's population. For estimation, since RD leverages local regression around the cutoff and gives more weight to units closest to the cutoff, omitting treatment switchers introduces bias and increases variance.
When deciding whether to include or omit treatment switchers, one must consider how units crossing the cutoff differ from units who did not. If in some sense these units are different in their potential outcomes, both omitting (selecting on units who did not cross) and including (composition effects) makes interpretation of the results harder. We illustrate this trade-off in the application in Section (ref).
Changes in the running variable over time may stem from manipulation, where units strategically affect their running variable to influence treatment assignment. While not the focus of this paper, multiple period data can help to detect manipulation. Existing approaches to detect manipulation in a multiple time-period setting extend methods developed for single-period RD designs. One method examines the continuity of covariates and outcomes around the cutoff in pre-treatment periods (cellini2010value), complementing traditional tests of covariate continuity (lee2008randomized). A second approach evaluates the continuity of the running variable’s density at the cutoff (mccrary2008manipulation), both within and across time periods (grembi2016fiscal).
We propose two new methods to detect manipulation with multiple time-period data. First, we propose to examine treatment switchers. Treatment switching may reflect manipulation if such changes correlate with covariates. Therefore, we propose to model treatment switching using pre-treatment covariates. If pre-treatment covariates are predictive of treatment switching, this may suggest the presence of manipulation.
Our second suggestion is motivated by the definition of manipulation given by mccrary2008manipulation. Let $\widetilde{R}_{i,t}$ be the value of the running variable had the RD treatment not taken place. We say that unit $i$ manipulated her running variable at time $t$ if $R_{i,t}\neq \widetilde{R}_{i,t}$. With this definition in mind, we propose to model the running variable non-manipulated values. First, fit a model for the running variable using observations from the pre-treatment time periods. Then, predict the running variable in the RD period. Under the assumption of no manipulation in the pre-treatment period, the model predictions can be interpreted as estimates of $\widetilde{R}_{i,t}$. If the model is correctly specified, the difference between observed and predicted running variable values provides a measure of manipulation, which can be further used. For example, for testing whether the mean difference at the cutoff is zero, or, characterizing units with a large difference measure.
So far, we have considered identification from an RD perspective. In some empirical studies with multiple time period data and a sharp RD treatment assignment, researchers employ a standard DID analysis, ignoring the RD treatment assignment (for recent reviews on DID see roth2023s,de2023two). Such analyses are carried out, for example, in studies of parental leave reforms (e.g., schonberg2014expansions, lalive2014parental, danzer2018paid). In this section, we compare between DID, RD and RD-DID when the underlying treatment mechanism is a sharp RD, for simplicity focusing on the case of two time periods. In Appendix (ref) we also analyze a model-based imputation approach (e.g., borusyak2024revisiting) that targets the same estimands as RD-DID.
Assume the canonical RD-DID setup with $\mathcal{T}_0=\{1\}$ and $\mathcal{T}_{\mathrm{RD}}=\{2\}$. We adopt the classical DID notation that $Y_{i,t}(0,w)=Y_{i,t}(w)$ (roth2023s). The target causal estimand in a DID design with two time periods is the global average treatment effect on the treated in the second time period, defined as $$ ATT(2)=\mathbb{E}\left[Y_2(1) - Y_2(0)\mid W_2 = 1,R\geq c\right]. $$ In RD designs, the target causal estimand is the local average treatment effect in the second time period, $$ATE(c,2)=\mathbb{E}[Y_2(1) - Y_2(0)\mid R = c].$$ In RD-DID designs, the target causal estimand is the local average treatment effect on the treated in the second time period, $$ATT(c,2)=\mathbb{E}[Y_2(1) - Y_2(0)\mid W_2=1, R = c].$$ Contrasting RD and RD-DID with DID, the main difference is in the global vs. local population of the target estimand. Whether a global or local estimand should be preferred, depends on the subject matter of the respective study and the research question at hand. Intuitively, the global estimand provides information on the entire (treated) population and so might be more relevant in certain cases. Put differently, the local estimand lacks external validity, if we are interested in the effect on the general population (lee2008randomized). Nevertheless, some research questions focus on local populations. For example, goldstein2023learning studies how crossing a round test score boosts the self-esteem of young and initially low-scoring prospective students and is beneficial in the long-term for their careers. The target estimand for a relevant policy in this context would not include older or initially high scoring prospective students, and hence a local estimand is preferred.
Identification in DID without covariates relies on a parallel trends assumption, namely that the mean trend in untreated potential outcomes is equal between the treated and untreated groups. Under the considered setup, this assumption can be written as
Conversely, the constant potential outcome discontinuity version of Assumption (ref), as presented in Corollary (ref)(i), implies that $\alpha_{0,2}=\alpha_{0,1}$, which can be rewritten as
Contrasting the DID and RD-DID identification assumptions, they are similar yet different. Similar, in that they make an assumption on the mean trend of untreated potential outcomes. Different, in that one is global in its nature and the other is local. The rationale behind RD designs is to replace indefensible global assumptions (i.e. assumptions on the entire population) with more plausible assumptions on the units that are local to the cutoff. Therefore, the main drawback we see in the parallel trends assumption is that it goes against this rationale and might be too ambitious in many applications due to its global nature. That is, under a sharp RD treatment assignment, we see the DID identification assumptions as stronger than the RD-DID identification assumptions.
To illustrate this point, Appendix Figure (ref) presents an illustrative example of a canonical RD-DID design. We compare two scenarios: first where the association between $Y_t(0)$ and $R$ is constant over time, and second where it is not. The identification assumptions underlying RD-DID hold in both scenarios, as RD-DID is agnostic to changes not at the cutoff. By contrast, the parallel trends is violated in the second scenario.
We propose a new estimation framework for RD with multiple time periods. Our proposal estimates the outcome discontinuity at each time period separately, and then constructs the estimator of the target estimand using the assumed $g_0$ and/or $g_1$ (Assumption (ref)).
A common approach in the literature is to pool multiple time periods into a single regression (e.g., grembi2016fiscal,avdic2018modern). A second, more recent approach implements a single RD regression to differenced outcomes between two time periods (picchetti2024difference).\footnote{However, it is not clear how to generalize the differenced outcomes approach when there are more than two time periods and/or under fuzzy RD design.} Our approach has three distinct advantages over these estimation frameworks. First, our estimation framework is directly motivated by the identification formula in Theorem (ref). Second, our approach allows us to use contemporary RD estimation methods, e.g., bias-correction procedures (calonico2014robust). Third, as we show below, variances under different sampling schemes can be studied. While in previous frameworks it is unclear how to estimate the variance under different sampling schemes, in our framework the distinction is made explicit and easy for the researcher to specify.
We focus here on the widely used constant potential outcome discontinuity case under sharp RD, formally presented in Corollary (ref)(i); we study a simplified version of the linear case (Corollary (ref)(ii)) in Appendix (ref). We present estimators for $ATT(c,t)$; the derivations for $ATU(c,t)$ are analogous. As is typical when analyzing RD estimators, the entire analysis is conditioned on the values of the running variable $R$. Detailed calculations are given in Appendix (ref).
For each period $t\in\mathcal{T}$, we propose to estimate the outcome discontinuity $D^Y_t$ by the local linear RD estimator for a single time period (hahn2001identification,porter2003estimation). Let $K(\cdot)$ be a kernel function, e.g., uniform or triangular, $h_n$ a bandwidth sequence, and $\boldsymbol{X}_p(r)=\boldsymbol{X}_p(r,c)=\left[1,(r-c)^1,...,(r-c)^p\right]^\prime$ a vector of polynomial terms with respect to the running variable centered at $c$.\footnote{One can use different kernel functions $K(\cdot)$ and bandwidths $h_n$ for the regression above $(+)$ and below $(-)$ the cutoff and at each time period $t$. Current practices for single time periods commonly choose the same bandwidth and kernel for both regressions (imbens2012optimal). Therefore, we also analyze an estimator with the same bandwidth and kernel.} For a given choice of $K(\cdot)$, $h_n$ and polynomial order $p$, the weighted least squares polynomial regression estimators on each side of the cutoff are
where the minus and plus subscripts denote the regressions below and above the cutoff, respectively. Let $\widehat{\beta}^{(v)}_{t,p,(+)}(h_n)$ and $\widehat{\beta}^{(v)}_{t,p,(-)}(h_n)$ be the $v$-entries in the coefficient estimators above and below the cutoff, respectively. The estimator for the outcome discontinuity at time $t$ is $\widehat{D}^Y_{t,p}(h_n)=\widehat{\beta}^{(0)}_{t,p,(+)}(h_n) - \widehat{\beta}^{(0)}_{t,p,(-)}(h_n)$. Following Corollary (ref), for a set of non-negative weights $\pi_t$, such that $\sum_{t \in t \in \mathcal{T}_0}\pi_t=1$, the ATT is estimated by
A natural question is how to choose the weights? A simple choice is to set the weight of the time period closest to $t^\star$ to one, and all others to zero. This choice will likely be more robust to deviations from the constant discontinuities assumption, if in practice there is some trend as we move to time periods further away from $t^\star$. On the other hand, this choice will result in a larger standard error compared to other choices, because it does not utilize all available data. A second option is to set all the weights according to a kernel, e.g., the uniform kernel which assigns equal weights. One consideration when choosing the weights is the sample size in each period. If, for example, more data is available from one time period compared to other time periods, one might opt for a weighting scheme that puts higher weight on that time period. If the sample size is equal in all time periods, a uniform kernel can be used. On the other hand, if there is a slight deviation over time from the constant discontinuities assumption, then the triangular kernel might be more robust. Hence the kernel choice may balance efficiency and robustness.
A third potential approach is to weight time periods according to their similarity to the time period of interest $t^\star$ with respect to covariates. A similar weighting scheme based on synthetic controls is proposed for DID estimation by arkhangelsky2021synthetic and is shown to have good properties. Two more possible weighting schemes are based on bias or variance considerations; we discuss them below.
We now derive the first-order asymptotic bias of $\widehat{ATT}(c,t^{\star};h_n,p)$ and propose a bias-corrected (BC) estimator for $ATT(c,t^{\star})$ (calonico2014robust). We assume the existence of one-sided $v$ derivatives of the mean outcomes above and below the cutoff, denoted by $\mu_{t,(+)}^{(v)}$ and $\mu_{t,(-)}^{(v)}$, respectively. Combining known results on the asymptotic bias of $\widehat{\beta}_{t,p,(+)}(h_n)$ and $\widehat{\beta}_{t,p,(-)}(h_n)$, presented in Appendix (ref), we can write the bias in the multiple time period estimator as
where $\mathtt{B}_{t,p}(h_n)=\frac{h_n^{2}}{2}\Big[\mu_{t,(+)}^{\left(2\right)}B_{t,(+),0,p,2}(h_n)-\mu_{t,(-)}^{\left(2\right)}B_{t,(-),0,p,2}(h_n)\Big]$ and the terms $B_{t,(+),0,p,2}$ and $B_{t,(-),0,p,2}$ are known expressions defined explicitly in Appendix (ref), and $h_n\rightarrow 0$ as $n$ goes to infinity. Returning to choice of $\pi_{t}$, a fourth option is to choose weights that minimizes the first-order bias.
Using the above bias formula, we can construct a BC estimator, computed in two steps. Note that in $\mathtt{B}_{t,p}(h_n)$, the only unknown quantities are $\mu_{t,(+)}^{\left(2\right)}$ and $\mu_{t,(-)}^{\left(2\right)}$. In the first step, $\mu_{t,(+)}^{\left(2\right)}$ and $\mu_{t,(-)}^{\left(2\right)}$ are estimated by $\widehat{\mu}_{t,(+)}^{\left(2\right)}=2\widehat{\beta}^{(2)}_{t,(+),q}\left(b_{n}\right)$ and $\widehat{\mu}_{t,(-)}^{\left(2\right)}=2\widehat{\beta}^{(2)}_{t,(-),q}\left(b_{n}\right)$, where $\widehat{\beta}_{t,(+),q}(b_n)$ and $\widehat{\beta}_{t,(+),q}(b_n)$ are estimated coefficient vectors, as in (ref), with bandwidth $b_n$ and polynomial order $q\geq2$. Denote the obtained estimated bias by $\widehat{\mathtt{B}}_{t,p,q}(h_n,b_n)$. In the second step, for each period we calculate the BC estimator of the outcome discontinuity by $\widehat{D}^Y_{\mathtt{BC},t,p,q}(h_n,b_n)=\widehat{D}^Y_{t,p}(h_n)-\widehat{\mathtt{B}}_{t,p,q}(h_n,b_n)$. Finally, the BC estimator of $ATT(c,t^{\star})$ is
We characterize the variance of the estimators, with and without the bias correction. An important implication of the multiple time period setup is that the data sampling scheme impacts the estimators' variance. We consider three different sampling types: a repeated cross-section (CS), a panel where the running variable is constant across time (PC), and a panel where the running variable varies over time (PV).
Let $V_{p}^{\mathrm{CS}}\left(t^\star;h_{n}\right)$, $V_{p}^{\mathrm{PC}}\left(t^\star;h_{n}\right)$, and $V_{p}^{\mathrm{PV}}\left(t^\star;h_{n}\right)$ be the variance of $\widehat{ATT}(c,t^\star;h_n,p)$ under CS, PC and PV sampling, respectively. The CS, PC and PV cases differ by the covariance of estimators above and below the cutoff across time periods. In the CS case, there is zero covariance between estimators at different time periods. In the PC case, the covariance is zero between estimators from different sides of the cutoff. In the PV case, all covariances can be non-zero. As we show in Appendix (ref), the variance of the estimator in each of these cases can be written as
where
The variances of the outcome discontinuity estimators and the covariance of outcome discontinuity estimators between two time periods are given explicitly in Appendix (ref). Similar to the discussion in the previous subsection, a fifth option to set the weights $\pi_t$ is obtained by minimizing the variance of the estimator under the constraint that the weights sum to one. If the variances $V\big(\widehat{D}_{t,p}^{Y}(h_n)\big)$ are similar across $t$ and the covariances are negligible, the optimal choice is to take equal weights for all time periods.
Following previous results on RD estimation (hahn2001identification,calonico2014robust), to test the null hypothesis that $H_0:ATT(c,t^\star)=0$, one can use the statistic $\widehat{ATT}(c,t^\star;h_n,p) / \sqrt{\widehat{V}_{p}(t^\star;h_{n}})$ and similarly calculate confidence intervals $\big[\widehat{ATT}(c,t^\star;h_n,p)\pm z_{1-\alpha/2}\sqrt{\widehat{V}_{p}(t^\star;h_{n})}\big]$.
As thoroughly studied by (calonico2014robust), BC procedures with desired asymptotic confidence interval coverages are possible and often preferred over the standard estimator. We therefore derive the variance of the BC estimator in (ref). Since the final form is similar in spirit to the variance of the standard estimator but the derivation requires further technical arguments, we refer to Appendix (ref) for the explicit formulas.
We conduct a Monte-Carlo simulation study to assess the finite-sample performance of the proposed estimation framework and compare the different variance estimators. We consider three data generating processes (DGPs), corresponding to CS, PC and PV as defined in Section (ref). For simplicity, we consider the canonical sharp RD-DID.
We empirically motivate the DGP of the simulation study by taking parameters based on estimated values from the application in Section (ref), as done in imbens2012optimal,calonico2014robust. We summarize the main parts of the DGP here, and describe it in more detail in Appendix (ref).
For each unit $i$, we generate the running variable by $R_i=\left(B_i-0.375\right)\times5000$, with $B_i$ sampled from Beta(2,4) distribution. Under the CS DGP, a sample of different $n$ units is generated at each time period, and for each unit we draw $R_{i,t}$ from the above distribution. For PC DGP, we draw a sample of $n$ units, $R_{i,1}$ is drawn from the above distribution, and $R_{i,2}=R_{i,1}$. For PV DGP, we draw a sample of $n$ units, $R_{i,1}$ is drawn from the above distribution, and $R_{i,2}=0.97\times R_{i,1}+\tau_i$, with $\tau_i \sim N(153,410)$. In all three DGPs, the outcome model is $Y_{i,t}=f_{t}\left(R_{i,t}\right)+1_{\{R_{i,t}\geq0\}}\times D_{t}+\delta_{i}+\delta_{t}+\varepsilon_{i,t}$, where $f_t(.)$ is a fifth order polynomial, $1_{\{.\}}$ is the indicator function, $D_{t}$ is the outcome discontinuity at the cutoff, $\delta_i$ is a unit fixed effect, $\delta_t$ is a time fixed effect, and $\varepsilon_{i,t}$ is an idiosyncratic error for each unit and time period. The values and distributions of these parameters are discussed in Appendix (ref). We take constant potential outcome discontinuities (as in Corollary (ref)(i)). Hence, the true effect is set as the difference of the outcome discontinuities at the cutoff between the two time periods, i.e., $ATT(c,2)=D_2-D_1$.
For each combination of DGP type (CS, PC, and PV) and sample size $n\in\{500,1000\}$, we simulate 1,000 samples. In each sample, we estimate the ATT using both the ”conventional” estimator given by (ref) and the BC estimator given by (ref), with $p=1$ and $q=2$. We estimate the variances assuming each of the three sampling types, for the conventional estimator using (ref)-(ref) and for the BC estimator using (ref) presented in Appendix (ref). We consider $h$ values between 200 to 1,000, and $b=2h$. We calculate confidence intervals using $z=1.96$ (targeting 95% coverage) as explained in Section (ref). Finally, we calculate the empirical coverage rate as the proportion of confidence intervals containing the true parameter across simulations.
The results are reported in Table (ref). Starting from the CS DGP, the variance estimator corresponding to the correct DGP is $\widehat{V}^{\mathrm{CS}}$. However, the variance estimator assuming the DGP is PC or PV, $\widehat{V}^{\mathrm{PC}}$ or $\widehat{V}^{\mathrm{PV}}$, have only negligible differences, since the estimated covariances over time periods between estimators are very small. The empirical coverage rates of the BC estimator are closer to the desired 95% than the coverage of the conventional estimator, similar to results previously obtained for the single time period RD setup (calonico2014robust). A similar result on the difference between the empirical coverage rates of the conventional and BC estimators can be observed for the PC and PV DGPs. A stark difference between the variance estimators is found in the PC DGP, where coverage rates are too high when estimated using variance that assumes CS data. This is due to the high covariance of estimators across time, which is correctly subtracted from the standard error estimator under the PC and PV specifications.
We illustrate our theoretical results by revisiting grembi2016fiscal, who examined how fiscal laws affect fiscal outcomes at the municipality level. grembi2016fiscal answer this question by analyzing a reform, enacted in Italy in 2001, which required municipalities with a population ($R$) larger than 5,000 ($c$) to abide by an annual deficit growth target ($W$). Therefore, treatment is assigned by a sharp RD. Here, we study the two primary outcomes ($Y$), fiscal gap (total expenditures minus total revenues net of transfers and debt services, in euros) and deficit (total expenditures minus total revenues, in euros). A challenge in employing RD in this application is that numerous regulations at the municipality level are set in motion depending on population size. One of these is the salary of the mayor, which also changes discontinuously at the 5,000 population threshold. Therefore, Assumption (ref) is not defensible. Fortunately, data from multiple time periods is available, which calls for the RD-DID design. Below we present the estimand and identification assumptions in this study, and briefly summarize the main results. We then turn to falsification tests in the RD-DID context. We finish with the issue of time-varying running variable (Section (ref)), and its impact on identification and estimation in this application.
The sample consists of 1,375 municipalities, observed across $1997,...,2004$.\footnote{The sample size is around 1,200 each year, with a minimum of 1213 and a maximum of 1246, since not all municipalities have available data in all years. The sample construction process mimics grembi2016fiscal, in that we do not use years 1997 and 1998 for the main analysis, and we drop observations with $R_{i,t}\leq3,500$ and $R_{i,t}\geq7,000$.} These time periods are classified into the following sets: $\mathcal{T}_1=\{1999,2000\}$, $\mathcal{T}_{\mathrm{RD}} =\{2001 , ... , 2004 \}$, and $\mathcal{T}_0=\emptyset$. Therefore, $t_1<t^\star$ for all $t_1\in\mathcal{T}_1$ and $t^\star\in\mathcal{T}_{\mathrm{RD}}$. Our target estimand is the $ATU(c,t^\star)$ for $c=5,000$ and $t^\star\in\mathcal{T}_{\mathrm{RD}}$.\footnote{As discussed in Section (ref), when allowing for carry-over effects, the definition of the ATU depends on the possible treatment paths. Therefore, in this application, for any year $t^\star\in \{2001,2002,2003,2004\}$ the ATU is defined as $ATU\left(c,t^{\star}\right)=\lim_{r\rightarrow c^{-}}\mathbb{E}\left[Y_{t^{\star}}\left(1,1,1,1,1,1\right)\mid R=r\right]-\lim_{r\rightarrow c^{-}}\mathbb{E}\left[Y_{t^{\star}}\left(1,1,0,0,0,0\right)\mid R=r\right].$} The ordering of time periods makes Assumption (ref) (Section (ref)) likely to hold. Consequently, we can use $t_1\in\mathcal{T}_1$ to identify $\alpha_{1,t^\star}$. To this end, by Theorem (ref), to identify $ATU(c,t^\star)$, we need to posit the functional form of $g_1$, the function connecting the unidentifiable bias $\alpha_{1,t^\star}$ with the identifiable $\alpha_{1,t}$ for $t\in\mathcal{T}_1$.
Appendix Figure (ref) presents the estimated yearly outcome discontinuities for the fiscal gap and the deficit. It presents visual evidence that fiscal rules lower fiscal gaps of municipalities with population size around 5,000. Table (ref) panel A reports the estimated $ATU(5000,t^\star)$, for $t^\star \in \mathcal{T}_{\mathrm{RD}}$, assuming constant potential outcome discontinuities (i.e., constant $\alpha_{1,t}$). Estimates are presented both for the conventional estimator (ref) and the BC estimator (ref). Standard errors under CS, PV and PV are calculated according to (ref)-(ref) for the conventional estimator and according to (ref) for the BC estimator, noting that in the application the data is of type PV. The estimated bias and variance were approximately the same in 1999 and 2000 so we took the weights $w_{1999}=w_{2000}=0.5$. Similar to the results from grembi2016fiscal, four years post implementation, the conventional (BC) estimated average treatment effect on the untreated is -23 (-25) with SE of 13 (14) under PV sampling for deficit, and -103 (-113) with SE of 36 (40) under PV sampling for fiscal gap, suggesting that fiscal rules improve fiscal outcomes at the municipality level.
We now study the plausibility of the constant discontinuity assumption and consider the alternative linear-in-time discontinuities.
For the deficit outcome, the estimated discontinuity $D^{Y}_{t}$ increases from 3 in 1999 to 23.2 in 2000 (Table (ref)A), an increase of 673%. In comparison, for the fiscal gap we observe an increase from 63.4 to 88.9, an increase of 40%. If the discontinuity in either outcome that is unrelated to the fiscal law policy increases over time, assuming a constant $g_1$ will result in underestimating (in absolute terms) the ATU. A formal testing procedure can help assessing the constant $g_1$ assumption. Denote $\Delta = \alpha_{1,2000} - \alpha_{1,1999}$. The null hypothesis of constant $g_1$ is $H_{0}:\Delta = 0$. The calculated t-statistic is $1.46$, which does not provide sufficient evidence to reject $H_{0}$ at the 10% significance level. The similar time-constant discontinuity hypothesis for the fiscal gap outcome is also not rejected in the 10% level, with a t-statistic of $0.62$.
As recently discussed in other designs (hartman2018equivalence,hartman2021equivalence,bilinski2018nothing) not rejecting such falsification tests (here, tests for $H_{0}$) does not provide evidence for the assumption in question (here, the time-constant discontinuity). These tests control the type I error -- the error of falsely stating the constant discontinuity assumption does not hold -- while we would like to control the type II error -- the error of falsely not rejecting the constant discontinuity assumption. To alleviate these concerns, we introduce an equivalence testing procedure, which tests the hypothesis $H_{\delta,0}:|\Delta|>\delta$ for a specified non-negative value $\delta$. We test $H_{\delta,0}$ using two one sided t-tests (TOST) (hartman2018equivalence). Rejecting $H_{\delta,0}$ means there is evidence in the data that $-\delta\leq \Delta \leq \delta$, namely that the identification assumption of constant discontinuity is approximately (up to $\delta$) correct. We consider the following equivalence testing procedure (hartman2018equivalence). We test $H_{\delta,0}$ for $\delta=0.36\times\sigma$, where $\sigma$ is the standard deviation of the outcome, using two one sided t-tests (TOST) each at $\alpha$ significance level.\footnote{The estimated standard deviation of deficit and fiscal gap in 1999 and 2000 ranges from 38.6 to 46.7. For comparison, we also study what equivalence does the data support, by calculating the minimal $\delta$ that will reject $H_{\delta,0}$. For deficit the minimal $\delta$ is 43.1, and for fiscal gap the minimal delta is 93.8. A similar procedure is proposed to suggest at equivalence of pre-trends in a DID setting by liu2024practical.} If both TOST are rejected, we conclude that the data supports equivalence of $\delta$ at $1-2\alpha$ significance levels. We find that for 10% significance level $H_{\delta,0}$ is not rejected for $\delta=0.36\sigma$ for both the deficit and fiscal gap outcomes. Therefore, using the aforementioned equivalence testing procedure, the $\alpha_{1,t}$ values are not equivalent between 1999 and 2000 for both deficit and fiscal gap using 10% significance level.
As an alternative to time-constant discontinuities, we use a simple linear approximation, as discussed in Corollary (ref)(ii). The formulation of the conventional and BC estimators, and their variances, for the linear-in-time assumption is presented in Appendix (ref).\footnote{Due to the data limitations of the application, with only two periods observed in $\mathcal{T}_1$, the linear-in-time discontinuity assumption is not testable, and the slope is calculated using only two data points. Hence we see these results as suggestive, and present them for illustration purposes. If more time periods were available, a linear model could be fit to the estimated outcome discontinuities of each period in $\mathcal{T}_1$.} Panel B of Table (ref) reports the results of this analysis. The estimated treatment effects under linear-in-time discontinuity are stronger (more negative). For example, the estimated $ATU(5000, 2004)$ for deficit became about 400% larger in absolute terms.
As discussed in Section (ref), composition effects may impact the analysis when the running variable is time varying. Consider the target estimand $ATU(5000,2004)$. The parameter $\alpha_{1,1999}$ is identifiable from the 1999 data, and assuming constant $\alpha_{1,t}$, such that $\alpha_{1,2004}=\alpha_{1,1999}$, identification of the ATU follows from Corollary (ref). Following Section (ref) and Equations (ref)--(ref), the constant $\alpha_{1,t}$ assumption can be re-written as $\Delta Y_{2004,1999}+\Delta R_{2004,1999}=0$.
Such an assumption is not easily defensible if either $\Delta Y_{2004,1999}$ or $\Delta R_{2004,1999}$ are non-zero, since this implies that they cancel out. A more plausible assumption is that $\Delta Y_{2004,1999}=\Delta R_{2004,1999}=0$. Note that $\Delta Y_{2004,1999}$ is unobserved in the data, since $\lim_{r\rightarrow c^{-}}\mathbb{E}[Y_{2004}(1)\mid R_{2004}=r]$ is unobserved. However, $\Delta R_{2004,1999}$ is identifiable from the data, since in $t=1999$ for all units $Y_{i,t}=Y_{i,t}(1)$. We estimate $\Delta R_{2004,1999}$ as follows. Let $D_{s,t}=\lim_{r\rightarrow c^{+}}\mathbb{E}\left[Y_{t}\mid R_{s}=r\right]-\lim_{r\rightarrow c^{-}}\mathbb{E}\left[Y_{t}\mid R_{s}=r\right]$. We estimate $D_{2004,1999}$ and $D_{1999,1991}$ using local linear regressions, presented in Section (ref), denoted $\widehat{D}_{2004,1999,p}(h_n)$ and $\widehat{D}_{1999,1991,p}(h_n)$. The estimator for the composition effect is $\widehat{\Delta} R_{2004,1999,p}(h_n)=\widehat{D}_{2004,1999,p}(h_n)-\widehat{D}_{1999,1991,p}(h_n)$. If $\widehat{\Delta} R_{1999,2004,p}(h_n)$ is non-zero and statistically significant this is suggestive for composition effects. Replacing 1999 as the baseline year with 2000, the composition effect is $\Delta R_{2000,2004}$ and its estimator is $\widehat{\Delta} R_{2004,2000,p}(h_n)$. Since the estimator is a difference between two single time period outcome discontinuities, it is equivalent to the proposed estimation approach under constant potential outcome discontinuities, and hence estimation and variance estimation can be carried out by our conventional and BC procedures.
Table (ref) reports point estimates and standard errors of estimated composition effects, for both 1999 and 2000 as baseline years.\footnote{In the considered application, since the observed $R_{i,t}$ does not change between $t^\star\in\mathcal{T}_{\mathrm{RD}}$, the estimated composition effects are the same for all target periods. That is, repeating the analysis for $t^\star\in\{2001,...,2004\}$ will produce equivalent $\widehat{\Delta}R$.} For deficit, focusing on the conventional estimator, taking 1999 as the baseline year, the composition effects are of small magnitude, $-4$ (SE of 12). However, using 2000 the estimated composition effects are larger, $-21$ (SE of 20). Compared to the estimated $\widehat{ATU}(c,2004)=-23$ (Table (ref)), the estimated composition effects are non-negligible in magnitude, although they are not significant at the 5% level. For fiscal gap we estimate $-65$ (SE of 36) using 1999 and $-93$ (SE of 39) using 2000. Again, compared to $\widehat{ATU}(c,2004)=-103$, the magnitude of the composition effects is non-negligible. Also, in 2000, the estimate is statistically different from zero. These results suggest that composition effects pose a problem for identification of $\alpha_{1,t^\star}$ for $t^\star\in\mathcal{T}_{\mathrm{RD}}$, and hence of $ATU(c,t^\star)$, in this context. To summarize, we found evidence of composition effects that are not small in magnitude but mostly not significant. And, if we assume composition effects enter the identification assumption, e.g., $\alpha_{1,t^\star}-\alpha_{1,2000}=\Delta R_{t^\star,2000}$, the estimated ATUs become much smaller (in absolute terms).
Given the above results, researchers might be tempted to drop treatment switchers from the sample. In our application, 137 (11%) municipalities changed their treatment status. The running variable, population size, is measured using a census prior 2001 and a census post 2001, and hence for each municipality there is one observed value before the reform, and one observed value after the reform. Figure (ref) visualizes the distribution of the population size across municipalities in each census. Comparing the distributions before (in blue) and after (orange) removal of treatment switchers, the density at each census after omitting treatment switchers displays a drop around the cutoff. Hence, omitting treatment switchers may cause the researcher to interpret such evidence as manipulation, as it will plausibly not pass a Mc'Crary test (mccrary2008manipulation). Furthermore, omitting treatment switchers produces a type of analysis similar to “donut-hole” RD, due to the drop in the density of the running variable around the cutoff.
The drop in the distribution before the reform (left panel of Figure (ref)) appears slightly to the left of the cutoff, whereas the dip after the reform (right panel of Figure (ref)) appears slightly to the right of the cutoff. This pattern suggests that most treatment-switching municipalities experienced population growth. Figure (ref) shows that most of the 137 municipalities that switched treatment status are initially between 4,500 and 5,500 population size, with some outliers below the 4,500 mark. Out of the 137 treatment-switching municipalities, 102 municipalities moved from below the cutoff to above the cutoff, while only 35 municipalities moved in the opposite direction, confirming the above suspicion.
In this paper, we studied identification and estimation of designs combining RD without the continuity assumption, with multiple time periods, termed RD-DID. We formulated the bias in RD studies when continuity does not hold, developed a general identification framework for RD-DID and discussed several key extensions. We then compared the proposed RD-DID framework with DID, and derived estimators for the causal estimands and variances of the estimators under several sampling schemes. Finally, we studied the finite-sample performance of the estimation framework in a simulation study, and illustrated the utility of the identification and estimation approach in an application.
Our framework can be applied when the data can be classified into three sets analogous to $\mathcal{T}_0$, $\mathcal{T}_{\mathrm{RD}}$ and $\mathcal{T}_{1}$, even when these sets are not defined based on time periods. For example, one can consider comparing between a hospital where units are treated according to some sharp RD assignment and a hospital where the RD is not implemented, and no unit, or all units, are treated.
Furthermore, while our formulation focused on the continuity approach to RD, a second common approach for identification in RD designs is the local randomization framework (lee2010regression). Extending our framework to the local randomization approach is left for future research.
Our work assumed the cutoff is fixed across time periods. Future research may consider time-varying cutoffs, and how they might be used to identify causal parameters. For example, multiple cutoffs might allow identification of treatment effects at different points across the support of the running variable, as discussed in cattaneo2021extrapolating. In addition, we briefly discussed the possible use of model-based imputation combined with the multiple time-period data to identify treatment effects which are far from the cutoff. Although promising, a more rigorous study is warranted to understand the implications of such assumptions on identification and estimation.
This study has shown the benefits of incorporating multiple time period data into an RD design. We believe that this paper will equip researchers with a principled framework on how to approach RD designs when the continuity assumption is violated.
\printbibliography