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.
78,552 characters · 24 sections · 23 citation commands
Potential Outcome Modeling and Estimation in DiD Designs with Staggered Treatments
{
} }
{\it Keywords: Bayesian inference, Bernstein-von Mises theorem, Causal inference, Regularization, Random effects, Shrinkage prior.}
\spacingset{1.8}
Difference-in-Differences (DiD) methods are now routinely applied in settings with multiple time periods and staggered treatment adoption. In such designs, standard two-way fixed effects (TWFE) regressions generally fail to recover causal treatment effects of interest. When treatment timing varies across units and effects are heterogeneous, TWFE averages across comparisons that include already-treated units serving as controls for later-treated units, yielding coefficients that lack a causal interpretation deChaisemartin_dHaultfoeuille2020AER,Goodman2021JoE,SunAbraham2021dynamic,BorusyakJaravelSpiess2024Restud,Roth2023JoE.
In recent work, estimators of the group-time Average Treatment Effect on the Treated (ATT) have been developed based on proper comparisons between the treated and untreated units. These effects are identified under parallel trends and no anticipation, and are estimable by nonparametric or doubly robust methods CallawaySantAnna2021did,SantAnnaZhao2020JoE. These methods have clarified identification and estimation in staggered designs, but they can be unstable when cohorts are small. A related issue is that when ATTs are defined at finer categories, the effective number of observations is smaller.
To address this problem, we propose a novel probabilistic model for potential outcomes that enforces parallel trends and no anticipation conditions directly at the model level. The model provides a precise definition of group-time ATTs and yields a likelihood-based framework for estimation and inference. We show that the model provides a common foundation for Bayesian and frequentist analysis. When treated cohort sizes are small, prior information on the ATTs can be incorporated through thick-tailed t-priors that shrink treatment effects of small magnitude toward zero. In addition, hierarchical priors can be used to enhance estimation when ATTs are defined not only at treatment cohort level, but also at finer sub categories. In other settings, one can use an iterated feasible GLS (IFGLS) estimator for the ATTs that mirrors the posterior sampling steps.
Unlike some of the recent works that estimates an untreated potential outcome model on untreated observations and then imputes counterfactual untreated outcomes for treated units to construct ATTs (BorusyakJaravelSpiess2024Restud; GardnerThakralToYap2025), we specify a unified model for both untreated and treated potential outcomes. Under this framework, the ATT vector is a deterministic function of the model parameters, which facilitates incorporating regularization and prior information.
In addition, our modeling offers a natural approach to assessing pre-treatment restrictions. Rather than relying on conventional pre-trend tests, which can have low power and induce pre-testing distortions FreyaldenhovenChristianShapiro2019AER,Roth2022AERI,Roth2023JoE, we can compare models with and without pre-treatment restrictions using marginal likelihoods and report results under the selected model. Posterior model probabilities quantify uncertainty about pre-trends that is propagated into inference on the ATTs.
On the theory side, we establish the large-sample behavior of the posterior distribution under correct specification and derive a Bernstein-von Mises theorem for the posterior of the parameters that justifies posterior inference for the treatment effects.
The remainder of the paper proceeds as follows. Section (ref) presents the model and reparametrization. Sections (ref) and (ref) develop extensions for covariate-dependent heterogeneity and for evaluating pre-treatment parallel trends via marginal likelihoods. Sections (ref)--(ref) describe estimation and inference, including posterior simulation, prior specification, asymptotic theory, and IFGLS. Section (ref) report simulations and an empirical application.
We consider an event-study setting in which units are treated at different times. We observe a random sample of $n$ units over $T$ time periods. For each unit $i\in\{1,\ldots,n\}$ and each period $t\in\{1,\ldots,T\}$, we observe an outcome $Y_{it}$ and a treatment indicator $X_{it}\in\{0,1\}$, where $X_{it}=1$ if unit $i$ is treated in period $t$ and $X_{it}=0$ otherwise. We assume that $\{(Y_{it},X_{it})\}_{t=1}^T$ are independently distributed across units.
Treatment is absorbing: once a unit is treated, its treatment status remains treated thereafter. That is, for all $1\le t<t'\le T$, $X_{it}\le X_{it'}$. Each unit's treatment path is therefore uniquely characterized by its first treatment period. Define the first treatment time (or cohort index) \[ S_i:=\min\{t\in\{1,\ldots,T\}:X_{it}=1\}, \] with the convention $S_i=1$ for never-treated units. We assume that never-treated units exist. Let $\mathcal{S}\subset\{2,\ldots,T\}$ denote the set of possible treatment timings among treated units. Then $S_i\in\{1\}\cup\mathcal{S}$ and there is a one-to-one mapping between $S_i$ and the treatment path $X_i=(X_{i1},\ldots,X_{iT})$:
We denote baseline covariates by $\bm w_i \in \mathbb{R}^{d_w}$. The vector $\bm w_i$ does not include a constant term for identification purposes. We do not use post-treatment covariates in order to avoid potential endogeneity.
Each unit $i$ is randomly assigned (ex-ante researcher's perspective) to one of the sequences in $\{1\}\cup \mathcal{S}$ defined above. Denote the POs (potential outcomes) for the sequence $s\in \{1\}\cup \mathcal{S}$, at time $t$ in two states: the treated PO, $Y_{s,it}^{(1)}$, and the untreated PO, $Y_{s,it}^{(0)}$. Let $Y_{s,it}$ be the observed outcome in sequence $s$ for unit $i$ at time $t$. These respect the standard notion of consistency. For the never-treated units, we observe untreated potential outcomes, and for the treated units, we observe untreated potential outcomes before adoption and treated potential outcomes after adoption. In other words, for $t=1,\ldots,T$,
Note that consistency is a mapping rule. It specifies which potential outcome is observed at each time for each unit. In contrast, the no-anticipation (NA) assumption, which we describe below, is a restriction on the potential outcomes themselves. Equation (ref) does not impose the NA restriction. It does not assert that $Y^{(1)}_{s,it}=Y^{(0)}_{s,it}$ for $t<s$. Rather, it states only that before treatment adoption, the observed outcome coincides with the untreated potential outcome, while after adoption, the observed outcome coincides with the treated potential outcome. The no-anticipation assumption is imposed separately below as a restriction on the potential outcomes.
We make the following standard assumptions for identification of the ATTs.
This assumption rules out anticipation effects for units that are eventually treated and is likely to hold when treatment timing is not chosen by the units themselves. Under this assumption, $\tau_{\mathrm{ATT}}(s,t)=0$ for all $t<s$.
This assumption states that, after treatment adoption, the counterfactual (untreated) outcome increments for treated sequences evolve in parallel with those of the never-treated sequence. That is, absent treatment, treated and never-treated units would have experienced the same expected changes in outcomes in each post-treatment period.
Relaxations of the parallel trends assumption are discussed in ye2024negative. In this paper, we maintain the standard parallel trends assumption and impose it explicitly in the specification of the potential outcomes model below.
We now describe our model for potential outcomes. In this model, we explicitly impose the parallel trend and no-anticipation assumptions (Assumptions (ref) and (ref)). In Section (ref), we extend this model to let the no-anticipation and parallel trend assumptions hold conditional on pre-treatment covariates. We begin with the never treated sequence. We let
The intercept $\beta_{11}$ is the initial value. The parameter $b_{12}$ measures the increment in the expected potential outcome from period 1 to 2. In general, for $k>2$, $b_{1k}$ is the increment from period $k-1$ to $k$ (the slope over the interval $[k-1,k)$). Next, for each treated sequence $s\in\mathcal{S}$, we let, for $t=1,\ldots,T$,
The initial value $\beta_{s1}$ is specific to the sequence $s$. The parameter $b_{sk}$ measures the expected increment of the treated potential outcome from period $k-1$ to $k$ in the sequence $s$.
By construction, this model satisfies NA and PT. NA is satisfied because $\E[Y^{(0)}_{s,it}]$ and $\E[Y^{(1)}_{s,it}]$ coincide for all $t<s$. PT is satisfied since the increments of the untreated outcomes in the treated sequences, $\E[Y^{(0)}_{s,it}-Y^{(0)}_{s,i,t-1}]$, for $t\ge s$, match those of the never-treated sequence, $\E[Y^{(0)}_{1,it}-Y^{(0)}_{1,i,t-1}]$. In this way, the model explicitly incorporates the restriction that, absent treatment, post-treatment outcome dynamics evolve identically across sequences.
This model has an elegant vector-matrix form. Define the $T\times 1$ vectors $\bm Y^{(0)}_{1,i} =(Y^{(0)}_{1,i1},\ldots,Y^{(0)}_{1,iT})'$, and $\bm Y^{(0)}_{s,i} =(Y^{(0)}_{s,i1},\ldots,Y^{(0)}_{s,iT})'$ and $\bm Y^{(1)}_{s,i} =(Y^{(1)}_{s,i1},\ldots,Y^{(1)}_{s,iT})'$, for $s\in\mathcal{S}$. Also, let \[ \bm \beta_s = \big(\beta_{s1}, b_{s2}, \cdots, b_{sT}\big)', \quad s\in \{1\} \cup \mathcal{S}. \] Finally, let $\bm L_T$ denote the $T\times T$ lower triangular matrix of ones, $\bm L_s^{\text{pre}}$ the matrix obtained by setting columns $s$ to $T$ of $\bm L_T$ to zero, and $\bm L_s^{\text{post}}$ the matrix obtained by setting columns 1 to $s-1$ of $\bm L_T$ to zero. For example, when $T = 5$ and $s = 4$, these matrices are \[ \bm{L}_T =
, \quad \bm{L}_s^{pre} =
, \quad \bm{L}_s^{post} =
. \] It is easily seen that $$ \bm L_T = \bm L_s^{\text{pre}}+\bm L_s^{\text{post}}, $$ which we use repeatedly. Then, the models in (ref), (ref), and (ref) can be written as
In this model, $\bm L_s^{\text{pre}}$ and $\bm L_s^{\text{post}}$ pick up the increments for the pre-treatment periods and post-treatment periods, respectively. We illustrate our model of potential outcomes in Figure (ref) for $\mathcal{S} = \{2,3,4\}$ and $T = 5$, along with the corresponding ATT’s. \FloatBarrier
\FloatBarrier
The models (ref), (ref), and (ref) allow for arbitrary evolution of the POs but are over-parameterized. We can achieve a more parsimonious formulation by letting \[ \bm \beta_s = \bm \beta_1 + \bm \delta_s, \quad s \in \mathcal{S}, \] where $\bm \delta_s$ represents the difference in increments between sequence $s$ and the never-treated sequence. Substituting $\bm \beta_s$ in (ref), and (ref), we get
For the purpose of constructing the likelihood, only equations (ref) and (ref) will be relevant.
The enhanced re-parametrization also simplifies the selection of priors. For instance, if one believes that the ATT is zero for every $t \geq s$, this can be encoded by centering the prior distribution of $\delta_{st}$ around zero. Similarly, ex-ante beliefs that reflect optimism or pessimism about the ATT can be constructed by concentrating the prior distribution around positive or negative values, respectively. Moreover, knowledge about the ATT from related studies can also be utilized.
We now extend the baseline potential-outcomes framework to allow for unit-level heterogeneity that depends on observed pre-treatment covariates. We first introduce heterogeneity in outcome levels through random intercepts. We then generalize the model further by allowing covariates to shift the entire never-treated counterfactual path through random slopes. In both cases, the identifying assumptions—no anticipation (NA) and parallel trends (PT)—are preserved, either unconditionally or conditionally on covariates.
When baseline covariates $\bm w_i$ are available, we incorporate heterogeneity in outcome levels through random intercepts whose distribution depends on $\bm w_i$. These random intercepts shift the entire outcome path for each unit without altering the identifying restrictions. Let $\alpha_{s,i}$ denote a unit-specific intercept for unit $i$ in sequence $s$. We model heterogeneity in levels by allowing $\alpha_{s,i}$ to depend on pre-treatment covariates according to
Given the random intercepts, the vectors of observed outcomes can be written as
where $\bm 1$ denotes the $T$-dimensional column vector of ones. Because the intercept enters both potential outcomes symmetrically, NA and PT continue to hold, and the ATT expressions derived earlier are unchanged.
The intercept-only specification allows covariates to affect the level of the outcome path but not its shape. A natural extension is to allow covariates to shift the entire never-treated counterfactual path. In this case, NA and PT are imposed conditionally on pre-treatment covariates $\bm w_i$ (see, for example, Assumptions 3 and 4 in CallawaySantAnna2021did). Our framework accommodates this extension by introducing covariate-dependent random slopes within the potential-outcomes model.
Throughout this subsection, expectations are understood as conditional on $\bm w_i$. The key point is that the extended specification preserves NA and PT conditionally on covariates. Treatment affects outcomes only through post-treatment increments, while covariate-dependent random effects shift the entire counterfactual path symmetrically across cohorts and, therefore, do not alter the identifying restrictions. Let $\bm a_{s,i}=(a_{s,i1},a_{s,i2},\ldots,a_{s,iT})'$, for $s\in\{1\}\cup\mathcal{S}$. We express the conditional mean vectors as
where $\bm a_{s,i}$ follows the covariate-dependent distribution
and $\bm \Gamma$ is a $T\times d_w$ coefficient matrix and $\bm D$ is a $T\times T$ covariance matrix. The corresponding unit-specific deviation in counterfactual levels is $\bm L_T\bm a_{s,i}$. This formulation nests the intercept-only random effects model as a special case. If $\bm a_{s,i}=(\alpha_{s,i},0,\ldots,0)'$, then $\bm L_T\bm a_{s,i}=\bm 1\,\alpha_{s,i}$ and the preceding expressions reduce to those in Section (ref). Using $\bm L_T=\bm L_s^{\text{pre}}+\bm L_s^{\text{post}}$ and writing $\bm \beta_s=\bm \beta_1+\bm \delta_s$, the conditional means can be expressed as
Note first that $\bm \Gamma$ is common across cohorts. We impose this restriction to ensure that conditional parallel trends hold by construction. Second, $\delta_{sk}$ are not functions of covariates. Relaxing this assumption yields conditional average treatment effects of the form $\tau_{\mathrm{ATT}}(s,t;\bm w)$. While such heterogeneity may be empirically relevant, recovering unconditional ATT parameters from these conditional objects is nontrivial, as it requires integration over the joint distribution of covariates and treatment timing. We leave the details of the MCMC algorithm for fitting this model for future work.
We now develop a Bayesian framework to evaluate whether the parallel trends (PT) assumption holds in the pre-treatment periods. While support for pre-treatment PT is often taken as evidence that the assumption is credible post-treatment CallawaySantAnna2021did, such assessments can induce frequentist pre-test bias when model specification depends on test outcomes. Also, the conventional implementation of the pre-test does not provide clear guidance when the hypothesis is rejected (Roth2023JoE). Even in such a case, the researcher may still want to estimate the treatment effect of interest.
We instead frame the question as one of model comparison. We estimate models that impose and relax the pre-treatment PT restriction and compare them using marginal likelihoods computed by the method of chib1995. This strategy measures the empirical support for the identifying PT assumption directly, without conditioning inference on a preliminary test. Recall that the parallel trends condition (Assumption (ref)) used to identify the ATT parameters is given by \[ \E[ Y^{(0)}_{s,it}-Y^{(0)}_{s,i,t-1} ] = \E[ Y^{(0)}_{1,it}-Y^{(0)}_{1,i,t-1} ]=b_{1t}, \; t\geq s. \] Suppose now this condition also holds for the pre-treatment periods (i.e.\ for $t<s$). Together with the no-anticipation condition (Assumption (ref)), this leads to \[ \E[ Y^{(1)}_{s,it}-Y^{(1)}_{s,i,t-1} ]=b_{1t}, \text{ for } t<s, \] which implies that $\delta_{st}=0$ for $t<s$. We call the model with these additional restriction `reduced' because of the fewer effective number of parameters in the vector $\bm \delta_s$.
Therefore, under the additional restriction of parallel trends on pre-treatment periods, the models (ref) and (ref) become
for $s\in \mathcal{S}$ and $t=1,\ldots,T$. See Appendix for the matrix form of the model.
It may be noted that if the researcher has a prior probability that the parallel trend condition holds (or does not hold) in the pre-treatment periods, then the corresponding posterior probabilities can be obtained based on the computed marginal likelihoods. The probabilistic statements on the ATT's can then be weighted according to the posterior probabilities. This appropriately reflects the ex-post uncertainty of the pre-trends, which is in contrast to the traditional approach that conditions on passing the pre-trends test and is known to have low power.
We focus on the main model because the fitting of the reduced model follows as a special case. We first proceed under the Bayesian approach to facilitate uncertainty quantification and the incorporation of prior information. This framework also provides valid inference without relying on asymptotic approximations, which is particularly advantageous in small samples, as demonstrated in the simulation exercises in Section (ref). Importantly, our model also serves as a foundation for frequentist inference, as we discuss in Section (ref).
For simplicity, we suppose that the errors are mean zero and Gaussian $\bm \varepsilon_{s,i} \sim N_T(\bm 0, \bm \Sigma_s), \ s\in \{1\}\cup \mathcal{S}$, where each of the $T \times T$ covariance matrices are assumed to be diagonal: $\bm \Sigma_s =\text{diag}( \sigma^2_{s1},\ldots, \sigma^2_{sT})$ for $s\in \{1\}\cup \mathcal{S}$. For further generality, one can allow the errors to be autocorrelated and follow a thick-tailed distribution. The estimation procedure described below would then be modified by steps detailed in chib1993bayes.
Let $N_s=\{i:S_i=s\}$ for $s=1,\ldots,T$, the set of units in the $s$th sequence and let $n_s$ be the cardinality of $N_s$. Let $\bm y_i = \bm Y_{S_i}$ be the $T$-dimensional vector of observed outcomes for unit $i$. If $S_i=1$, the likelihood contribution of unit $i$ is
If $S_i=s$ for $s\in \mathcal{S}$, the likelihood contribution of unit $i$ is
We can integrate out $\alpha_{si}$ to get
The collection of parameters to be estimated is given by \[ \bm \theta = \left( \bm \beta_1, \bm \Sigma_1, \{ \bm \delta_s, \bm \Sigma_s , s\in \mathcal{S}\}, \{\bm \gamma_s, D_s: s\in \{1\}\cup \mathcal{S} \} \right). \] Let $\bm Y=\{\bm y_i:i=1,\ldots,n \}$ and $\bm W=\{\bm w_i: i=1,\ldots,n \}$ be the observed outcomes and covariates. Given the observations, the likelihood is defined as \[ p(\bm Y \vert \bm \theta, \bm W) = \prod_{i \in N_1} p_1\left( \bm y_i \vert \bm \beta_1, \bm \Sigma_1, \bm \gamma_1, D_1\right) \cdot \prod_{s\in \mathcal{S}} \prod_{i \in N_s} p_s\left( \bm y_i \vert \bm \beta_1,\bm \delta_s, \bm \Sigma_s, \bm \gamma_s, D_s \right). \]
The priors are independently specified as $\bm \beta_1\sim N_T(\bm \mu_{\beta_1}, \bm V_{\beta_1})$, $\bm \delta_s\sim N_{T}(\bm \mu_{\delta_s}, \bm V_{\delta_s}), \ s\in \mathcal{S}$, and for each $s\in \{1\} \cup \mathcal{S}$, $\bm \gamma_s\sim N_{d_w}(\bm \mu_{\gamma_s}, \bm V_{\gamma_s})$, $ D_s\sim \text{InvGam}(a_{D_s}/2,b_{D_s}/2)$, and $\sigma^2_{st}\sim \text{InvGam}(a_{st}/2,b_{st}/2), \ t=1,\ldots,T$. Section (ref) discusses hyperparameters selection. \
Let $\bm \alpha=(\alpha_{S_11},\ldots,\alpha_{S_nn})'$ be the vector of random effects. The joint distribution of the observed and unobserved variables is
where $\pi(\bm \theta)$ is the prior density. Integrating out the random-effects, the joint distribution can be written as
A MCMC sampler is used to efficiently sample the posterior distribution. Posterior inferences of objects of interest (e.g.\ ATT's) are based on the sample of draws produced by the algorithm. The posterior sample consists of \[ \bm \theta^{(g)} = \left( \bm \beta_1^{(g)}, \bm \Sigma_1^{(g)}, \{ \bm \delta_s^{(g)}, \bm \Sigma_s^{(g)} , s\in \mathcal{S}\}, \{\bm \gamma_s^{(g)}, D_s^{(g)}: s\in \{1\}\cup \mathcal{S} \} \right), \ g=1,\ldots, G, \] where $G$ is the number of MCMC draws (beyond a suitable burn-in). The algorithm repeats the following steps (details are supplied the Appendix (ref)):
The prior on the increment differences $\delta_{st}$ needs to be carefully chosen, as the objects of interest, such as the ATT's, are functions of these parameters. We recommend prior selection in two cases that empirical researchers are likely to encounter in practice: (1) pre-training the prior hyperparameters based on a subsample when the sample size is sufficiently large, and (2) an adaptive shrinkage prior for small sample sizes.
When the sample size is sufficiently large, we recommend the following pre-training approach to set the priors. The idea is to use a random sample of the available data to train the prior, ie., to estimate the hyperparameters, similar to an empirical Bayes approach. We use $15\%$ of the data for this purpose in the simulations and application below. Estimation and inference is then on the remaining data.
When the sample is small, the training prior approach may not be feasible. A similar issue occurs in the existing frequentist approaches. For example, the did package of CallawaySantAnna2021did returns an error message if the number of observations in a group is less than or equal to 5. The advantage of the proposed model-based approach is that the prior can be chosen to reflect an ex-ante belief of the researcher that the ATT's should not be too far off from zero. We suggest use of a student-t prior, which is given hierarchically as follows:
Note that the prior variance $V_{\delta_st}$ of $\delta_{st}$ has its own prior specified as an inverse-gamma distribution, instead of being fixed. Integrating over the $V_{\delta_st}$'s, one can show that the marginal prior on $\delta_{st}$ is a student-t distribution centered around zero. It is natural to ex-ante believe that the increment differences $\delta_{st}$ (and hence the ATT's) are centered around zero. This prior distribution encourages $\delta_{st}$ of small magnitudes to shrink towards zero while the thick tails allow for possible larger effects. See, for example, ArmaganZaretzki2010 for more discussion on the adaptive shrinkage with the t-priors and their relation to ridge regression. A review of shrinkage priors and their applications in economics can be found, for example, in KorobilisShimizu2022.
In this case, $V_{\delta_st}$ can be easily updated from the inverse-gamma distribution as an additional step in the MCMC: \[ V_{\delta_st}\vert \bullet \sim \text{InvGam}((\rho+1)/2,(\xi+\delta^2_{st})/2) \quad t=1,\ldots,T, s\in \mathcal{S}. \] We specify the priors for the other parameters as: $\beta_{1t}\sim N(0,10)$, $\sigma^2_{st}\sim \text{InvGam}(1/2,1/2),$ $\gamma_{s}\sim N(0,10),$ and $D_s\sim \text{InvGam}(1/2,1/2)$, independently.
In many applications, units can be grouped into not only high-level categories that determine treatment timings (i.e.\ cohorts), but also sub-categories. In such cases, sub-category ATTs might be of interest. For instance, in the minimum wage example we visit below, the units (i.e.\ counties) belong to different states. The treatment cohort $s$ for a specific county is determined by when the state it belongs to adapts the new policy. Multiple states can adapt the policy in a same year. Intuitively, the ATTs might be similar within a cohort, but for policy makers, state-specific ATTs can be useful information (see, e.g.,\ KarimWebb2024). However, an empirical challenge is that there would be fewer observations as we work with finer categories.
Our framework combined with hierarchical priors can be used for overcoming this challenge. To see this, let $r\in \mathcal{R}$ be a state and $s(r)$ be the cohort the state $r$ belongs to, where $\mathcal{R}$ is the collection of states. Our model for potential outcomes can be defined in the same way as before but with state level increments: for each $r\in \mathcal{R}$,
which leads to the vector of state level ATTs: $\bm \tau^{(r)}_{ATT}=\bm L_{s(r)}^{\text{post}} \bm \delta_{r}$. A hierarchical prior can be defined, for example, as follows:
where $\pi$ denotes a prior for cohort specific means $\bm \mu_{\delta_s}$ and variances $\bm V_{\delta_s}$. Such prior reflects the researcher's ex-ante belief that the state-specific ATTs are similar within the cohort and helps the estimation by pooling information at the cohort level rather than estimating the state ATTs completely separately. This approach can naturally extends to even finer subcategories (e.g.\ cities) by adding extra layers into the hierarchical structure. We leave the implementation of the hierarchical priors for sub-category ATTs for future work.
We note that the ability to incorporate regularization on ATTs, whether it is shrinkage towards zero or hierarchical structure for sub-category ATTs, is a unique feature of our unified model for both treated and untreated potential outcomes. In a framework such as BorusyakJaravelSpiess2024Restud and GardnerThakralToYap2025 that estimates a model for untreated potential outcomes and imputes the missing untreated potential outcomes for the treated, introducing such regularization is not straightforward.
In this section, we derive the frequentist properties of the proposed Bayesian approach. In particular, we establish a Bernstein-von Mises type theorem for the ATTs. This indicates that the Bayesian credible sets of the ATTs, which can be obtained easily from the draws of the simple MCMC algorithm, have frequentist interpretations.
For notational simplicity, we establish theoretical results under the original parametrization. All the results apply to the reparamatrized model. Let $\bm \beta$ and $\bm \gamma$ be the vectors with elements $\{\bm \beta_s\}_{s \in \{1\} \cup \mathcal{S}}$ and $\{\bm \gamma_s\}_{s \in \{1\} \cup \mathcal{S}}$, respectively. Let $\bm \sigma^2$ and $\bm D$ be the vectors with elements $\{ (\sigma^2_{s1},\ldots,\sigma^2_{sT}) \}_{s \in \{1\} \cup \mathcal{S}}$ and $\{ D_{s} \}_{s \in \{1\} \cup \mathcal{S}}$, respectively. Recall that $\bm \Sigma_s = \text{diag}((\sigma^2_{s1},\ldots,\sigma^2_{sT}))$. For this section, we define the vector that collects all the parameters as $\bm \theta=( \bm \beta, \bm \gamma, \bm \sigma^2, D) \in \bm \Theta$. Let $\bm \theta^*$ be the true value of $\bm \theta$.
We consider the asymptotic framework with a fixed $T$ and an increasing $n$. Conditional on covariates $\bm w_i \in \mathbb{R}^{d_w}$ and treatment assignment $S_i=s\in \{1\} \cup \mathcal{S}$, a sequence of outcomes $\bm y_i=(y_{i1},\ldots,y_{iT})'\in \mathbb{R}^T$ is generated from the model $p_{\theta^*}$ defined as: \[ p_{\theta}(\bm y_i \vert \bm w_i, S_i=s)=N(\bm y_i \vert \bm X_i \bm \phi_s, \bm \Lambda_s), \] where $ \bm X_i=(\bm 1 \bm w_i', \bm L_T)$, $\bm \phi_s=(\bm \gamma_s',\bm \beta_s')'$, and $\bm \Lambda_s = \bm \Sigma_s + D_s \bm 1 \bm 1'$. The data contains outcomes $\bm y_i$, covariates $\bm w_i$, and treatment assignment $S_i$: $\bm D^n=\{\bm D_i=(\bm y_i, \bm w_i, S_i): i=1,\ldots,n\}$. The covariates $\bm w_i$ are iid and generated from $g^*$. The treatment assignments are iid and generated from $h^*=\{h^*_s: s\in \{1\} \cup \mathcal{S} \}$ with $h^*_s=\Pr(S_i=s)>0$ for $\forall s$. The positivity of $h^*_s$ is ensured by the definition of $\mathcal{S}$. We do not model $g^*$ and $h^*$. The joint probability measure implied by $p_{\theta^*}$, $g^*$, and $h^*$ is denoted by $F_0$. We first establish Bernstein-von Mises theorem for the model parameters $\bm \theta$.
Now we study the asymptotic behavior of the posterior distribution of the ATT's. Let $P$ denote the prior for $\bm \theta$ and $\bm \theta_n^P$ denote the random variable with the law equal to the posterior distribution of $\bm \theta$ given a sample $\bm D^n$ of size $n$. By Theorem 1, it converges in distribution to $N(\bm \theta^* + \frac{1}{\sqrt{n}}\bm \Delta_{n,\theta^*}, (n \bm I_{\theta^*})^{-1} )$. Let $\bm \tau_{ATT}$ be the vector with elements $\{ \bm \tau_{ATT}^{(s)} \}_{s\in \mathcal{S}}$, where $ \bm \tau_{ATT}^{(s)} = \bm L^{post} (\bm \beta_s - \bm \beta_1)$. Denote the mapping $\bm \theta \mapsto \bm \tau_{ATT}$ by $f$; i.e. $f(\bm \theta) = \bm \tau_{ATT}$. Define the random variable $\bm \tau_{ATT,n}^P=f(\bm \theta_n^P)$ and, similarly, $\hat{\bm \tau}_{ATT, n} = f(\hat{\bm \theta}_n)$. Note that the plug-in estimator is consistent, i.e, $\hat{\bm \tau}_{ATT, n} \overset{p}{\to} \bm \tau_{ATT}^*$, the true value of $\bm \tau_{ATT}$.
By applying the Bayesian delta-method (e.g.\ BernardoSmith1994bayesian, Section 5.3), we obtain the following as a consequence of Theorem 1.
In other words, the posterior of ATT's is asymptotically normal. Hence, Bayesian credible sets have asymptotically correct nominal coverage and are valid confidence sets. The priors introduced in Section (ref), including the Student-t prior that we advocate for improved robustness in small samples, satisfy the assumptions for the theoretical results.
The model in this paper also provides the basis for frequentist estimation and inference. For its simplicity, we consider an iterated generalized least squares approach. In this approach, one iterates between (i) joint generalized least squares updates of the mean parameters and (ii) closed-form updates of the variance components. The procedure is a direct frequentist analog of the conditional Bayes updates used in Appendix A, and it preserves the same identifying structure and interpretation of the ATT parameters.
Fix a balanced panel length $T$. Units are partitioned into groups indexed by $s\in\{1,\dots,S\}$. Let $n_s$ be the number of units in group $s$. Stack outcomes in the $T\times n_s$ matrix $ \bm Y_s = (\bm y_{s,1},\dots,\bm y_{s,n_s}), \bm y_{s,i}\in\mathbb{R}^T $, and let $\bm w_{s,i}\in\mathbb{R}^{d_w}$ denote the baseline covariate vector for unit $i$ in group $s$. Collect covariates in the matrix $\bm W_s=(\bm w_{s,1}',\dots,\bm w_{s,n_s}')'$, so $\bm W_s$ is $n_s\times d_w$. Let $\bm 1_T$ denote a $T\times1$ vector of ones. Under intercept heterogeneity, the outcome equation for unit $i$ in group $s$ is
where $\bm\beta_1\in\mathbb{R}^T$ is common across groups, $\bm\gamma_s\in\mathbb{R}^{d_w}$ is group-specific, $\alpha_{s,i}$ is a scalar unit-level random effect, and $\bm\varepsilon_{s,i}$ is a $T$-vector of idiosyncratic errors. We assume \[ \alpha_{s,i}\sim (0,D_s), \qquad \bm\varepsilon_{s,i}\sim (0,\bm\Sigma_s), \qquad \bm\Sigma_s=\mathrm{diag}(\sigma_{s,1}^2,\dots,\sigma_{s,T}^2), \] independently across $i$ and independently of $\alpha_{s,i}$. Marginalizing out $\alpha_{s,i}$ yields the covariance $\mathrm{Var}(\bm y_{s,i}) = \bm\Sigma_s + D_s \bm 1_T\bm 1_T' $. Define the precision matrix
which is the inverse of $\bm\Sigma_s + D_s\bm 1_T\bm 1_T'$ via the Sherman–Morrison formula. Now let \[ \bm\beta = (\bm\beta_1', \bm\gamma_1',\dots,\bm\gamma_S', \bm\delta_2',\dots,\bm\delta_S')' \] collect all mean parameters. Conditional on $\{\bm\Lambda_s^{-1}\}_{s=1}^S$, the joint GLS normal equations are obtained by stacking the group-wise quadratic forms \[ \sum_{i=1}^{n_s} (\bm y_{s,i}-\bm\mu_{s,i}(\bm\beta))' \bm\Lambda_s^{-1} (\bm y_{s,i}-\bm\mu_{s,i}(\bm\beta)), \] where $\bm\mu_{s,i}(\bm\beta)$ denotes the conditional mean implied by (ref). The resulting first-order conditions define a linear system $ \bm A(\bm\Lambda^{-1})\,\bm\beta = \bm b(\bm\Lambda^{-1}) $, where $\bm A(\cdot)$ and $\bm b(\cdot)$ aggregate contributions across all groups. Solving yields the joint GLS estimator \[ \hat{\bm\beta}(\bm\Lambda^{-1}) = \bm A(\bm\Lambda^{-1})^{-1}\bm b(\bm\Lambda^{-1}). \] This update is the exact frequentist analog of the Gaussian conditional posterior mean updates for $(\bm\beta_1,\bm\gamma_s,\bm\delta_s)$ in the Bayesian sampler.
Next, given the current mean estimates, define the residual matrices \[ \bm E_s = \bm Y_s - (\bm L_1\hat{\bm\beta}_1)\bm 1_{n_s}' - \bm 1_T(\bm W_s\hat{\bm\gamma}_s)' - \mathbb{I}\{s\ge2\}(\bm L_1\hat{\bm\delta}_s)\bm 1_{n_s}'. \] The best linear unbiased predictor (BLUP) of the random effects is \[ \hat{\alpha}_{s,i} = D_s \frac{\bm 1_T'\bm\Sigma_s^{-1}\bm e_{s,i}} {1+D_s\,\bm 1_T'\bm\Sigma_s^{-1}\bm 1_T}, \] where $\bm e_{s,i}$ is the $i$th column of $\bm E_s$. Variance components are updated by
with lower bounds imposed to ensure numerical stability. This procedure converges rapidly in practice and produces estimates that numerically match the posterior means from the Bayesian sampler given in Appendix A.
Finally, ATT estimation and inference proceeds as follows. For treated groups $s\ge2$, the vector of average treatment effects on the treated is $ \bm\tau_s = \bm L_s\bm\delta_s $. Standard errors are obtained by the delta method using the covariance matrix of $\hat{\bm\beta}$ implied by the final GLS system, and confidence intervals are constructed in the usual way. Unlike existing frequentist procedures that estimate group-time ATTs sequentially or through sample splitting, the proposed iterated feasible GLS estimator exploits the full system of mean restrictions implied by the potential outcomes model and estimates all parameters jointly. This yields a transparent estimator that targets the same ATT parameters as the Bayesian approach under identical identifying assumptions. Thus, our model provides a unified frequentist/Bayesian framework for estimation and inference in staggered DiD designs.
We conduct simulation studies to illustrate the performance of the proposed approach in both large and small sample cases. In summary, the proposed approach is comparable to the existing frequentist method of CallawaySantAnna2021did when the sample is large, and it outperforms the existing approach in small samples. The latter case underscores the advantage of the proposed model-based approach: shrinkage priors can be used to improve estimation when the existing approaches may struggle with small sample sizes.
We generate the data sets that, roughly speaking, mimic the real data that we use in the empirical application in Section (ref). The data set contains five periods (i.e. $T=5$), and units are treated in the second, fourth, and fifth periods (i.e. $\mathcal{S}=\{2,4,5\}$). Thus, there are 7 causal parameters of interest, as listed in the first column of Table (ref). We consider two sample sizes: a medium $n=250$ and moderately large $n=500$ values. The $n$ units are randomly allocated into the 4 treatment sequences with probabilities $(0.4,0.2,0.2,0.2)$. We generated $w_{i}\sim N(3.3,1)$ to roughly match the moments of the log population in the real data in Section (ref). To mimic the log teen employment in the real data, we use estimated parameter values to generate the observed outcomes based on our model. We make $|\delta_{ts}|$ slightly larger than the estimates for the pre-treatment periods in order to ensure that the DGP is realistic in the sense that the parallel trends do not hold in pre-treatment periods.
We compare 4 approaches: (1) CallawaySantAnna2021did (CS21), (2) our baseline approach (Bayes), (3) our approach that imposes the parallel trend condition on pre-treatment periods as in Section (ref) (Bayes-PrePT), and (4) our approach chosen by the marginal likelihood (Bayes-ML). We use the R package, did, for CS21. For each approach, we report Bias, root-mean-squared-errors (RMSE), empirical coverage, and length of the 95% confidence and credible intervals (Cov and IL). In addition, we report the average of the p-values in CS21 for testing the null hypothesis that the parallel trend holds in pre-treatment periods, as well as the number of times each of our proposed models is selected by the marginal likelihood (ML). We repeat the experiments 500 times.
The upper and lower panels of Table (ref) present the simulation results for $n=250$ and $n=500$, respectively. First, regarding the estimation of the ATT's, comparing the results between CS21 and our Bayes model, we see that the two approaches are comparable in most measures, especially when the sample is large enough (i.e. $n=500$), confirming our theoretical finding in Section (ref).
Second, the p-value of the null hypothesis of parallel trends in pre-treatment periods is, on average, quite small, as expected. In contrast to the conventional implementation of the pre-test which does not provide clear guidance when the hypothesis is rejected, we compare the two models with and without the restriction. With a medium sample size (i.e. $n=250$), ML selects the baseline model (Bayes) over the restrictive model (Bayes-PrePT) about 77% of the time. The correct model is the former. Although mis-specified, the Bayes-PrePT model is sometimes preferred, probably because it has fewer parameters than the Bayes model. When the sample is large enough (i.e. $n=500$), ML selects the correct Bayesian model almost always, as expected.
In Table (ref), we show the results when the DGP actually imposes the parallel trend for pre-treatment periods. In this case, ML selects the restrictive model almost always. The bias and RMSE of this model are smaller compared to CS21 and Bayes, especially under $n=500$ due to the model having fewer parameters.
Recall that when $n_s$'s are sufficiently large, we recommend to `pre-train' the hyperparameters. When they are small, this approach is not feasible. Similarly, CS21's did package can not be implemented when $n_s<6$. In this simulation study, we examine the performance of the proposed estimator when sample size is small.
In small samples, our recommendation (Section (ref)) is to use a student-t prior centered around zero. Such prior reflects the ex-ante belief of the researcher that the ATT's should not be too far off from zero but still allowing for thick tails for possible outliers. In small samples, we do not recommend computing ML due to known issues.
We generate data sets just as in the previous simulation setting, but we let $n=24$ with $n_s=6$ for $s=1,2,4,5$. This is the smallest sample size that the CS21 package can take. We compare CS21 with our Bayesian approach under two priors on $\delta_{st}$: (1) default prior and (2) t-priors. In the default prior, we set $\delta_{st} \sim N(0,10)$ and in the t-priors, we let $\delta_{st} \sim N(0,V_{\delta_st})$ and $ V_{\delta_st}\sim \text{InvGam}(\rho,\xi)$ with $(\rho,\xi)=(1,1)$. Here we focus on the performance of the point estimators. The asymptotic coverage of confidence/credible intervals are not guaranteed for the given sample size.
Table (ref) shows the results. CS21 and Bayes-default both suffer from large biases and RMSEs of similar sizes, as expected. Clearly, Bayes-t prior gives smaller biases/RMSEs than these alternatives as the benefit of the shrinkage. On average, for Bayes-t, $\bm V_{\delta_2}^{-1}=(14.5,77.8,20.1,18.5,28.7)$, $\bm V_{\delta_4}^{-1}=(75.1,34.8,18.2,78.1,8.5)$, $\bm V_{\delta_5}^{-1}=(35.4,159.7,30.2,53.3,65.8)$, introducing substantial shrinkage while the thick tails of the t-distribution allow for outliers. Note that Bayes-default fixes these elements at 0.1.
Somewhat surprisingly, Bayes-default gives slightly smaller RMSE than CS21. A potential reason for this is that even the prior with a large fixed variance $(=10)$ can help regularize the parameter space in small samples. Note that both priors are centered around zero and yet provide sensible performance when the true ATT's are non-zero.
This simulation study shows the benefits of our model-based formulation of staggered DiD designs. It allows us to incorporate prior knowledge (shrinkage) that enables efficient estimation of the ATT's and valid uncertainty quantification, even when the sample size is too small for existing frequentist methods to deliver reliable inference.
We illustrate the proposed approach on real data regarding minimum wage policy and teen employment in the U.S. that we obtained from CallawaySantAnna2021did, (CS21)\footnote{\url{https://bcallaway11.github.io/did/}}. This is a subset of the dataset analyzed in their paper. During the period between 2001 and 2007, the federal minimum wage was flat at \$5.15 per hour. They focus on county-level teen employment in states where the minimum wage was equal to the federal minimum wage at the beginning of the period. Some of these states increased their minimum wage over this period. The counties in these states are the treated units. In particular, the treatment sequences are defined by the time period when a state first increased its minimum wage. Other states did not increase their minimum wage and the counties in these states belong to the never-treated sequence. The data includes 500 counties (i.e.\ $n=500$): 309 counties were never-treated, 20 counties were treated in 2004, 40 were treated in 2006, and 131 were treated in 2007. The outcome variable is the log of county-level teen employment. We use log population as a pre-treatment baseline covariate (i.e.\ $w_i$). The data spans from 2003 to 2007 (i.e.\ $T=5$). We construct the prior from 15% of the data. This data is randomly selected from each sequence. The remaining 85% of the data are used for estimation. See Appendix for the trained hyperparameters. We compare three methods: (1) CS21, (2) our baseline approach (Bayes), and (3) our approach with pre-treatment parallel trends (Bayes-PrePT).
Table (ref) and Figure (ref) show estimation results. First, as expected, our baseline model (Bayes) produces results similar to CS21. The estimated sign and magnitude of the ATT's are all similar between the two methods. For both models, $\text{ATT}(2,4)$ and $\text{ATT}(2,5)$ are `significant' in the sense that the confidence/credible intervals exclude zero. $\text{ATT}(4,5)$ is additionally significant in the proposed model. Last, in the Bayes-PrePT model, the posterior standard deviations tend to be smaller for the sequences that the pre-treatment periods PT would be relevant, namely sequences 4 and 5. This is because, under this condition, fewer parameters are required to be estimated in these sequences. In this model, $\text{ATT}(5,5)$ is additionally found to be significant. The posterior credible intervals of the ATT's from the two models are given in blue in Figure (ref).
Second, the table shows the difference in differences (DiDs) for the pre-treatment periods, denoted as $\text{PreDiD}(s,t)\equiv\sum_{k=2}^{t}\delta_{sk}$, for $t<s$. It simply measures the pre-treatment period DiDs, which is zero if the trends are exactly parallel before the treatment\footnote{When PT holds in the pre-treatment periods, we have
where the first equality is due to no-anticipation. It states that the difference in differences (DiDs) of the observed outcomes over subsequent periods is zero. This implies that $ \E[Y^{(1)}_{st}-Y^{(1)}_{s1}]=\E[Y^{(0)}_{1t}-Y^{(0)}_{11}], \ t<s. $ i.e. the DiDs over period 1 and period $t$ are zero. In our model, this means that $ \text{PreDiD}(s,t)\equiv \sum_{k=2}^t \delta_{sk}=0, \ t<s. $ } and is often used in empirical contexts to check the validity of the parallel trends. Figure (ref)(a) shows them in green. This quantity is zero in Bayes-PrePT. All the DiDs for pre-treatment periods are found to be non-significant except for $\text{PreDiD}(5,2)$ whose 95% credible interval slightly includes zero (the lower bound is very close to zero). This finding is actually similar to CS21 which finds this object to be only `weakly' non-significant in the sense that the confidence interval almost misses zero. In any event, the decision rule that we recommend is easy and clear: to look at the comparison between the models with and without the parallel trends in pre-treatment periods. The log marginal likelihood, computed by the method of chib1995, is substantially larger for the model that assumes the parallel trends over pre-treatment periods i.e.\ Bayes-PrePT. Indeed, with equal prior model probabilities of 0.5, the posterior probability for this smaller model is approximately 1.
This paper develops a unified model for treated and untreated potential outcomes for DiD designs with multiple periods and staggered treatment adoption. We model the evolution of potential outcomes while respecting the key identifying assumptions of parallel trends and no anticipation. The formulation clarifies the structure of causal comparisons and yields parameters that are directly interpretable as differences-in-differences.
We also incorporate unobserved heterogeneity through random effects. In one specification, we allow for sequence-specific random intercepts, and in the other, random effects whose distributions depend on baseline covariates. This structure provides a flexible way to capture persistent differences across units and treatment cohorts.
Our model provides a common foundation for both frequentist and Bayesian inference. The latter is especially attractive in applications where some group-time cells may contain few observations. Prior information on ATTs can be incorporated through training-sample priors, adaptive shrinkage priors, or hierarchical priors. We establish a Bernstein–von Mises theorem showing that, under correct specification, the posterior distribution of the ATTs is asymptotically normal and that Bayesian credible sets have correct frequentist coverage.
We also show that we can obtain an iterated generalized least squares estimator of the ATT parameters. By combining explicit modeling of the potential outcomes, independent or crossed random effects, and unified inferential tools, this framework offers a promising new approach for the analysis of DiD designs with staggered treatments.
Our approach is implemented in a user-friendly software package, bdid, to support its use by practitioners. It is available for MATLAB, R, and Stata.