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.
67,965 characters · 16 sections · 45 citation commands
Sequential Synthetic Difference in Differences
Keywords: event studies, synthetic control, difference in differences, interactive fixed effects, panel data, sequential analysis.
JEL Codes: C21, C23, C38.
Event study designs, where a researcher observes units before and after an event alongside a comparison group, are a cornerstone of modern applied economics currie2020technology,goldsmith2024tracking. These studies routinely employ the difference-in-differences (DiD) strategy angrist2008mostly,Bertrand2004did,card1990impact,card1994minimum. The credibility of DiD, however, hinges on a parallel trends assumption, requiring that counterfactual outcomes for the treated group would have evolved in parallel with the outcomes of the comparison group. In the common case of staggered treatment adoption, most modern estimators, while flexible in accommodating heterogeneous treatment effects, still rely on this fundamental assumption borusyak2024revisiting,CallawaySantAnna,de2020two,sun2020estimating.
We address this limitation by proposing the Sequential Synthetic Difference-in-Differences (Sequential SDiD) estimator, a method for staggered adoption designs that is robust to violations of the parallel trends assumption. Our proposal adapts the Synthetic DiD estimator of arkhangelsky2021synthetic, which combines ideas from the Synthetic Control (SC) and DiD literature abadie2003,abadie2010synthetic. The procedure is iterative: it first averages outcomes for units sharing the same adoption date, then sequentially estimates treatment effects. After each estimation step, it uses the result to impute the missing counterfactual for the treated cohort. This imputed data is then used in the analysis of subsequent cohorts, preventing bias from cascading through the estimation.
We analyze Sequential SDiD within a linear model with interactive fixed effects (IFE), which explicitly models the unobserved confounders that drive violations of parallel trends arellano2011identifying,bai2009panel,chamberlain1992efficiency,Freyberger,holtz1988estimating,pesaran2006estimation, and has been routinely used to establish statistical properties for the SC-type estimators abadie2010synthetic,arkhangelsky2021synthetic,ben2021augmented,ben2022synthetic,ferman2019synthetic,ferman2021synthetic. Our analysis departs from much of the SC literature by employing an asymptotic regime with a fixed number of periods and a large number of units. This allows us to connect the SC literature to the classic moment-based econometrics panel data literature (going back at least to chamberlain1984panel). This framework is suitable for our setting because the initial averaging of outcomes within large cohorts effectively reduces the noise level in the problem, connecting our work to the "low-noise" regime in the SC literature Hirshberg.
Our analysis is built around an infeasible "oracle" OLS estimator. This oracle knows the unobserved interactive fixed effects and uses them as regressors---a benchmark that captures what applied researchers attempt to do with proxies like unit-specific trends. We prove that our feasible Sequential SDiD estimator is asymptotically equivalent to this infeasible oracle. This equivalence is the paper's key theoretical result, as it immediately implies that our estimator is asymptotically normal and unbiased, providing a basis for standard inference. It also delivers, to our knowledge, the first formal efficiency guarantees for an SC-type method, showing that Sequential SDiD is asymptotically as efficient as the oracle OLS estimator.
A key intermediate step in our proof, which is of independent interest, is that we derive an alternative representation for this oracle regression. We show that it can be implemented via a sequential algorithm, where at each step the oracle constructs a weighted DiD estimator and uses the result for imputation. This finding provides new insight into the mechanics of modern imputation estimators, such as that of borusyak2024revisiting for the two-way fixed effects model. More broadly, this sequential representation of a correctly specified oracle opens up possibilities for developing other new methods that relax standard DiD assumptions. Our work focuses on one such relaxation (allowing for IFE), but other routes, such as relaxing selection assumptions, could be built on a similar sequential logic.
Our work offers contributions to several strands of the applied and theoretical econometrics literature.
First, we introduce a new, transparent, and robust estimator for event studies with staggered adoption. A key feature is its sequential algorithm, where estimates for early-adopting cohorts are used to impute counterfactuals that, in turn, inform the analysis of later ones. This iterative structure provides a cohesive analysis of the entire event study panel, distinguishing it from related SC approaches for staggered designs that effectively treat different adoption cohorts as separate problems (e.g., ben2022synthetic, cattaneo2021prediction).
Second, we provide an answer to a long-standing critique of SC methods. A common question is why one should rely on synthetic control weights instead of directly estimating the underlying factor model that motivates the procedure. Our results resolve this tension by formally demonstrating that our SC-based method is asymptotically equivalent to an oracle OLS estimator that directly uses the unobserved factors. This clarifies that, in this setting, the choice is not between balancing and direct estimation, but rather how to feasibly approximate the same ideal oracle benchmark. Crucially, this equivalence allows us to establish efficiency guarantees for an SC-type method. These results complement the optimal regret properties of the SC method under adversarial sampling derived by chen2023synthetic.
Finally, we advance the literature on difference-in-differences with staggered adoption. While recent methods have addressed the issues with heterogeneous treatment effects, they still effectively operate within a standard two-way fixed effects model that implies the parallel trends assumption. Our estimator also accommodates unrestricted treatment effect heterogeneity, but it remains robust when the parallel trends assumption is violated due to the presence of unobserved interactive fixed effects.
We demonstrate our method's performance with an empirical application and two data-based simulation studies. First, we revisit the analysis of Community Health Centers by bailey2015war and show that Sequential SDiD produces results that are reassuringly similar to standard DiD in a setting where the parallel trends assumption is believed to hold. We then use two sets of simulations to assess performance when this assumption is violated by confounding interactive fixed effects. The first is calibrated to the county-level data from the application. The second, which extends the simulation design of arkhangelsky2021synthetic to the staggered case and is inspired by the seminal work of Bertrand2004did, uses a state-level wage panel where treatment assignment is correlated with the underlying factors. In both simulation designs, we find that while standard DiD becomes severely biased, our Sequential SDiD estimator provides reliable estimates and valid inference.
Our analysis has two main limitations. First, our theoretical framework requires that adoption cohorts are relatively large. This size is crucial as it allows us to leverage the law of large numbers; by averaging outcomes within a cohort, we can substantially reduce the influence of idiosyncratic noise. This approach makes our statistical results possible under otherwise mild assumptions and provides a reasonable approximation for many empirical settings. It is also commonly used in the modern DiD literature (e.g., abadie,callaway2020difference).
Second, our main results rely on the assumption that idiosyncratic errors are independent across units, a condition that ensures they concentrate around zero upon aggregation. The appropriateness of this assumption depends on the context. While some influential work has relied on it for inference card1994minimum, other research highlights that common shocks can induce error correlation that survives aggregation Bertrand2004did. Our model is designed to address this concern. The interactive fixed effects capture precisely these common shocks, and recent evidence suggests that this structure can account for the vast majority of variation in aggregated data, leaving a much smaller, plausibly independent residual error arkhangelsky2021synthetic. Nonetheless, in applications where this separation is not clean---for instance, where idiosyncratic shocks themselves have heavy tails and drive aggregate outcomes---our statistical guarantees may not apply.\footnote{This type of granularity, where idiosyncratic shocks affect aggregates, is documented by gabaix2011granular in the context of the US economy.}
The rest of the paper proceeds as follows. Section (ref) presents the econometric setup and our estimator. Section (ref) establishes the theoretical results. Section (ref) discusses covariates. Section (ref) contains our empirical illustration and simulation experiments, and Section (ref) concludes.
\paragraph{Notation:} We use standard notation for expectations and variance operators, $\mathbb{E}[\cdot]$ and $\mathbb{V}[\cdot]$, respectively. For a sequence of random variables $X_n, Y_n$ we write $X_n = o_p(Y_n)$ if $\frac{X_n}{Y_n}$ converges to zero in probability. For two sequences $a_n$ and $b_n$ we write $a_n \lesssim b_n$ if $\frac{a_n}{b_n}$ is bounded, and $a_n \ll b_n$ if $\frac{a_n}{b_n}$ converges to zero. We write $a_n\sim b_n$ if $a_n \lesssim b_n$ and $b_n \lesssim a_n$. We use the same notation with subscript $p$ for the corresponding concepts for random sequences.
This section lays out our proposed methodology. We begin by describing the econometric framework, including our assumptions about potential outcomes and the data-generating process with interactive fixed effects. We then formally present the Sequential Synthetic Difference-in-Differences (Sequential SDiD) estimator. We conclude by outlining procedures for statistical inference and model validation.
We observe $n$ units over $T$ periods, with $i$ being a generic unit and $t$ being a generic period. In our theoretical analysis, we treat $T$ as fixed and $n$ as going to infinity -- the asymptotic regime that provides a reasonable approximation for a large class of empirical applications. For each unit $i$ and period $t$, we observe a real-valued outcome $Y_{i,t}\in \mathbb{R}$ and a binary treatment indicator $W_{i,t}\in \{0,1\}$. Following most of the applied work, we focus on settings with staggered adoption, thus assuming $W_{i,t} \ge W_{i,t-1}$.
We formalize causality by interpreting the observed outcomes using potential outcomes neyman1923,rubin1974estimating,imbens2015causal.
This assumption incorporates two separate restrictions. The first is no anticipation: only treatments realized by period $t$ can affect the outcomes in that period. The second is the absence of cross-unit spillovers, a key component of the Stable Unit Treatment Value Assumption (SUTVA), meaning potential outcomes for unit $i$ only vary with that unit's own treatment assignment. See arkhangelsky2024causal for a discussion of these assumptions.
Given our focus on staggered adoption designs, we re-index potential outcomes by adoption time. We begin by defining the set of all possible staggered treatment paths, $\mathbb{W}:= \{\mathbf{w}\in \{0,1\}^T: w_t\ge w_{t-1}\}$. For any such path $\mathbf{w} \in \mathbb{W}$, we can define its corresponding adoption time as $a(\mathbf{w}):= \inf \{t: w_t =1\}$. This mapping creates a one-to-one correspondence between treatment paths in $\mathbb{W}$ and adoption times in the set $\mathbb{A}:= \{1,\dots, T, +\infty\}$. This allows us to define potential outcomes indexed by adoption time, $Y_{i,t}(a) := Y_{i,t}(\mathbf{w}^t(a))$, where $\mathbf{w}^t(a)$ is the history of the unique treatment path corresponding to adoption time $a$. For each unit $i$, we denote its observed adoption time by $A_i := \inf\{t \le T: W_{i,t} = 1\}$. The observed outcome can then be written as a function of the potential outcome corresponding to the observed adoption time: $Y_{i,t} = Y_{i,t}(A_i)$. The internal consistency of this representation relies on Assumption (ref). In what follows, we use both representations interchangeably.
Our next assumption specifies the data-generating process for the potential outcomes. We model potential outcomes as a combination of four components: (i) a standard unit fixed effect, $\alpha_i$; (ii) a standard time fixed effect, $\beta_t$; (iii) a low-rank interactive fixed effect, $\theta_i^\top \psi_t$, which captures unobserved confounding factors that violate the parallel trends assumption; and (iv) the treatment effect itself, $\tau_{i,a,k}$, which is allowed to be heterogeneous across units, adoption cohorts, and time since treatment.
The interpretation of Assumption (ref) depends on the underlying sampling scheme and treatment assignment protocol. We now discuss two scenarios that justify this assumption and clarify the meaning of the treatment effect parameters $\tau_{i,a,k}$.
\paragraph{Example 1: Latent Unconfoundedness.} Suppose that $\{Y_i, A_i\}_{i=1}^n$ are $n$ i.i.d. draws from a population. For each unit, suppose there is a latent characteristic $U_i$ such that a conditional independence (unconfoundedness) assumption holds:
Further, suppose the never-treated potential outcomes satisfy the IFE structure:
We can then define the treatment effect parameter as the average causal effect conditional on $U_i$, and define the error term $\epsilon_{i,t}$ as the residual component:
Under these definitions, the model for the observed outcome $Y_{i,t}$ conforms to the structure in Assumption (ref).
Independent sampling is standard in the panel data literature and is natural for environments where units of observation correspond to economic agents, such as individuals or firms. It is also common in the theoretical DiD literature (e.g., abadie2005semiparametric, callaway2020difference, sant2020doubly). The latent unconfoundedness restriction (ref) is a form of strict exogeneity: conditional on $U_i$, the adoption date is as good as random. Strict exogeneity has a long tradition in panel data econometrics chamberlain1984panel and underlies analyses based on parallel trends assumptions (see ghanem2022selection for a discussion).
\paragraph{Example 2: Aggregate and Idiosyncratic Shocks.} Consider a fixed set of $n$ units over $T$ periods. In each period $t$, unit $i$ is exposed to an unobserved aggregate shock $H_t$ and an idiosyncratic shock $\nu_{i,t}$. Suppose the potential outcomes have the following structure:
where $\boldsymbol{h}^{t}$ and $\boldsymbol{u}^t$ are the histories of aggregate and idiosyncratic shocks. This decomposition is general; the term $\tau_{i,a,k}(\cdot)$ is simply the difference $Y_{i,a+k}(a, \cdot) - Y_{i,a+k}(+\infty, \cdot)$, implying that shocks do not restrict the causal effects.
Let $\boldsymbol{\nu}_i$ and $\boldsymbol{H}$ be the realized shock vectors. We make two independence assumptions: (i) idiosyncratic shocks are independent of adoption times conditional on aggregate shocks, and (ii) unit-level data is independent across units conditional on aggregate shocks:
Furthermore, suppose the conditional expectation of the untreated component has an IFE structure:
Finally, we define the parameters and the error term by conditioning on the aggregate shocks $\boldsymbol{H}$:
The resulting model for $Y_{i,t}$ again satisfies Assumptions (ref) and (ref).
This second example may seem less familiar, but it captures features common in applications where aggregate shocks drive outcomes and treatment assignment. First, it allows treatment trajectories to be correlated across units. Second, it permits the errors $\epsilon_{i,t}$ to have a factor structure (i.e., be correlated), since they are only independent “conditional” on $\boldsymbol{H}$. This structure is prevalent in large panel data models (e.g., bai2009panel,pesaran2006estimation) and is related to research on shift-share designs adao2019shift.
In both examples, the parameter of interest $\tau_{i,a,k}$ represents a specific conditional average treatment effect. In Example 1, we condition on the latent characteristic $U_i$ that drives selection. In Example 2, we average over idiosyncratic shocks $\boldsymbol{\nu}_i$ while conditioning on aggregate shocks $\boldsymbol{H}$. Our analysis will focus on estimating averages of $\tau_{i,a,k}$ for subpopulations of adopters. These examples illustrate that the interpretation of these averages depends on the underlying probability model.
In this section, we introduce the new estimator, which we call Sequential SDiD. As the name suggests, it is based on sequential application of a version of the SDiD estimator introduced in arkhangelsky2021synthetic. The key difference, though, is that we apply SDiD principles to aggregated data. Let $\mathcal{A}$ be the support of $A_i$; for each adoption cohort $a \in \mathcal{A}$, we define aggregate outcomes:
where $n_a := \sum_{i=1}^n \mathbf{1}\{A_i = a\}$ is the total number of units in cohort $a$. We also define the corresponding shares, $\pi_a := \frac{n_a}{n}$. Assumption (ref) guarantees that these aggregate outcomes follow a similar interactive fixed effects model:
where $\epsilon_{a,t} :=\frac{\sum_{i : A_i = a} \epsilon_{i,t}}{n_a}$, and other variables are the cohort-level averages of their unit-level counterparts. Representation (ref) is the key to our algorithm, which we present formally in Algorithm (ref).
\RestyleAlgo{boxruled} \LinesNumbered
Algorithm (ref) estimates the treatment effects $\tau_{a,k}$ sequentially, starting with the contemporaneous effect ($k=0$) and proceeding to longer horizons. For each cohort $a$ and horizon $k$, the procedure applies the principles of the SDiD estimator. It first constructs unit-specific weights $\hat{\omega}^{(a,k)}$ and time-specific weights $\hat{\lambda}^{(a,k)}$ (Line 3). The unit weights create a synthetic control for cohort $a$ using a weighted average of not-yet-treated cohorts ($j>a$). The time weights create a synthetic match for the event period $a+k$ using a weighted average of pre-event periods ($l < a+k$). A key difference from the original SDiD estimator is that we do not impose non-negativity constraints on these weights. The estimator $\hat{\tau}_{a,k}^{SSDiD}$ is then formed using these weights in a generalized DiD formula (Line 4).
The sequential nature of our method is driven by the imputation step in Line 5. After estimating $\hat{\tau}_{a,k}^{SSDiD}$, we use it to update the outcome $Y_{a,a+k}$. This step replaces the observed outcome with an estimate of the missing average counterfactual outcome. The imputation in Line 5 is the engine of the sequential procedure. When we move to estimate a longer-run effect for cohort $a$ (i.e., for $k' > k$), or any effect for a later cohort $a' > a$, the period $a+k$ becomes part of the pre-treatment history. However, the observed outcome $Y_{a,a+k}$ is contaminated by treatment. The imputation step $Y_{a,a+k} := Y_{a,a+k} - \hat \tau_{a,k}^{SSDiD}$ replaces this treated outcome with its estimated untreated counterfactual. This ensures that the weights computed in subsequent steps are based on a panel where the parallel trends assumption is enforced on the best available estimate of the counterfactual data, preventing the cascading bias that would arise from using raw treated outcomes as controls. This imputed value is then carried forward and used to construct weights and estimates for longer horizons (i.e., for $k' > k$). This process ensures that information from early estimates is efficiently incorporated into later ones.
The algorithm depends on three user-specified parameters: the regularization parameter $\eta$, the maximum horizon $K$, and the range of cohorts $(a_{\min}, a_{\max})$. The parameter $\eta \ge 0$ controls the penalty on the weights; its role is discussed further in our theoretical analysis. In practice, it serves to prevent the synthetic control from overfitting to idiosyncratic noise in the pre-treatment periods, especially when the number of control units or pre-treatment periods is small. The choice of $K$ and $(a_{\min}, a_{\max})$ determines which treatment effects are estimated and must satisfy the data constraint $a_{\max} +K \le T$, since estimating $\hat \tau_{a_{\max}, K}$ requires observing outcomes in period $a_{\max} + K$.
We use the estimated cohort-specific effects to construct an average effect for each horizon $k$:
where the weights $\mu$ are user-specified. In our analysis, we focus on weights that are proportional to the cohort shares, $\mu_{a} = \frac{\pi_a}{\sum_{a' \in \{a_{\min},\dots, a_{\max}\}}\pi_{a'}}$, and use $\hat \tau_k^{SSDiD}$ to denote the resulting estimator.
To conduct inference on the estimated treatment effects, we rely on the Bayesian bootstrap rubin1981bayesian,chamberlain2003nonparametric. The procedure involves drawing a set of $n$ independent weights $\boldsymbol{\xi} := \{\xi_i\}_{i=1}^n$ from an Exponential distribution, $\xi_i \sim \text{Exp}(1)$. These weights are used to construct a bootstrapped version of the aggregate outcomes:
We then apply Algorithm (ref) to these re-weighted outcomes $Y_{a,t}(\boldsymbol{\xi})$ to produce a bootstrap replicate of our estimator, $\hat \tau_k^{SSDiD}(\mu, \boldsymbol{\xi})$. By repeating this process many times, we can simulate the conditional distribution of the estimator given the data, and use the quantiles of this distribution to form confidence intervals.
A simpler alternative, which we use in our simulations, is to construct normal-approximation confidence intervals. This involves computing the standard deviation of the bootstrap replicates,
and using it to form a standard Wald-type interval:
where $q_{1-\alpha/2}$ is the $1-\frac{\alpha}{2}$ quantile of the standard normal distribution. While we expect such intervals may be conservative in some settings hahn2021bootstrap, we find that they perform reasonably well in our simulations.
We also adapt our method for placebo validation, analogous to testing for pre-trends in a standard DiD analysis. For a chosen placebo horizon $P > 0$, we define a set of placebo adoption times by shifting the true adoption times backward:
We then construct aggregate outcomes based on these placebo cohorts $A_i^P$ and apply Algorithm (ref) to this new dataset, setting the maximum horizon to $K = P-1$. Under the model's identifying assumptions, the resulting placebo estimates, $\hat{\tau}_{a,k}^{SSDiD}$ for $k \in \{0, \dots, P-1\}$, should be centered at zero. Significant deviations from zero would cast doubt on the validity of the model.
In this section, we establish the theoretical properties of our estimator. Our analysis proceeds in two main steps: we first introduce an infeasible "oracle" OLS estimator that serves as a natural benchmark, and then we show that our Sequential SDiD estimator is asymptotically equivalent to this oracle. This approach has two distinct advantages. First, the oracle estimator we consider is grounded in established empirical practice, making the link between our procedure and what applied researchers might ideally want to do explicit and immediately valuable. Second, the analysis of the oracle uncovers new algorithmic properties of OLS in staggered adoption settings, a result of independent interest. All proofs are collected in the Appendix.
We begin our theoretical analysis with an infeasible "oracle" OLS estimator constructed from the aggregate data, $Y_{a,t}$. This oracle knows the factor loadings $\theta_a$ and factors $\psi_t$ and uses them as regressors. Specifically, we consider the solution to the following weighted least squares problem:
The optimization in (ref) treats the factor components $\theta_{a}$ and $\psi_t$ as known regressors to be interacted with unknown cohort-specific and period-specific coefficients, $\nu_a$ and $\phi_t$. This can be viewed as a tractable linearization of the fully nonlinear interactive fixed effects problem. As we will show, the solution to this problem is intimately connected to our proposed Sequential SDiD estimator.
A priori, it is not obvious whether the parameters of interest, $\hat \tau^{OLS}_{a,k}$, are uniquely defined by (ref). In a standard two-way fixed effects model, the conditions for identification are straightforward, requiring only that a valid comparison group exists. The situation is more challenging here. In principle, it is possible for an effect at a longer horizon, $\hat \tau^{OLS}_{a,l}$ (for $l>k$), to be identified even if the shorter-horizon effect $\hat \tau^{OLS}_{a,k}$ is not. Our next assumption provides a sufficient condition to ensure orderly identification.
Assumption (ref) has a natural interpretation. For the time factors $\{\psi_t\}$, it requires that factors far enough in the future ($t \ge t^{\star}$) can be perfectly predicted by an affine combination of past factors. The requirement for the unit-specific loadings $\{\theta_a\}$ is similar: the loadings of early-adopting cohorts can be represented as an affine combination of the loadings of sufficiently late-adopting cohorts.
This assumption allows us to state our first main result, which is not only crucial for our subsequent analysis but also holds independent interest.
Proposition (ref) is important because it provides an explicit, sequential algorithm for computing the oracle OLS estimates. The structure of Algorithm (ref) deliberately mirrors that of our main proposal, Sequential SDiD (Algorithm (ref)). For each $(a,k)$, it constructs unit- and time-specific weights, $\tilde \omega^{(a)}$ and $\tilde \lambda^{(a,k)}$, and uses them to form a weighted DiD estimator. It then concludes by imputing the counterfactual outcome and proceeds sequentially.
\RestyleAlgo{boxruled} \LinesNumbered
The sequential representation of $\hat \tau_{a,k}^{OLS}$ in Algorithm (ref) is key to establishing the statistical connection between the oracle OLS and the feasible Sequential SDiD estimators in the next section. However, this result is also interesting in its own right as it reveals the underlying mechanics of OLS in this setting. For instance, if we consider a standard two-way model (i.e., $\theta_a \equiv \psi_l \equiv 0$), Assumption (ref) holds trivially. In this case, Algorithm (ref) provides a sequential implementation of the imputation estimator proposed by borusyak2024revisiting. The weights $\tilde{\lambda}$ become uniform, and the weights $\tilde{\omega}$ are inversely proportional to cohort size, reducing the procedure to a sequence of standard DiD estimators applied to imputed data.
The sequential nature of Algorithm (ref) also has practical implications, showing that $\hat \tau_{a,k}^{OLS}$ can, in principle, be computed "online" using only information available up to period $a+k$. While less critical for the small-scale problems we consider, this property is valuable in large-scale industrial applications. Finally, this representation opens several paths for generalizations of OLS, with our Sequential SDiD estimator being one such path. Other routes could include adding further regularization to the weights (e.g., a simplex constraint) or restricting the information set to allow for weaker exogeneity assumptions.
In this section, we formally connect our feasible Sequential SDiD estimator (Algorithm (ref)) to the infeasible oracle OLS estimator (Algorithm (ref)). As the procedural similarities suggest, the two are closely related. To establish this statistical relationship, however, we require some mild restrictions on the data-generating process.
We begin by imposing a regularity condition on the aggregate errors, $\boldsymbol{\epsilon}_a := (\epsilon_{a,1},\dots, \epsilon_{a,T})$.
Assumption (ref) is a mild condition that ensures the cohort-level noise diminishes at a standard $\frac{1}{\sqrt{n}}$ rate. It implicitly requires that the share of each adoption cohort, $n_a/n$, is non-vanishing as $n \to \infty$, which is central to our large-$n$, fixed-$T$ asymptotic framework.
To state our main result, we introduce additional notation. For each relevant cohort-horizon pair $(a,k)$, we define a matrix $L^{(a,k)}$ that captures the demeaned interactive fixed effects:
for cohorts $j >a$ and periods $l <a+k$. Here, $\overline{\theta}^{(a)}$ and $\overline{\psi}^{(k)}$ are the averages of the respective factor components over the relevant control units and pre-treatment periods. Assumption (ref) guarantees that for the $(a,k)$ pairs we consider, this matrix has full rank $r$. We use $\tilde \sigma_{(a,k)}$ to denote the smallest non-zero singular value of $L^{(a,k)}$.
We are now ready to state the main theorem, which establishes the asymptotic equivalence between the Sequential SDiD and the oracle OLS estimator.
Theorem (ref) provides the key theoretical justification for our method. The conditions are relatively mild. We require the oracle OLS estimator to be well-behaved (Condition 2), which is natural for a benchmark. The key restriction is on the strength of the factors (Condition 3). This condition requires the identifying variation from the interactive fixed effects to be stronger than the statistical noise, but it allows the factors to be "weak" in the sense that their explanatory power can vanish as $n \to \infty$. Finally, we require the regularization parameter to be chosen appropriately (Condition 4), balancing the bias from regularization with the variance from overfitting.
As long as these conditions hold, our feasible Sequential SDiD estimator mimics the infeasible oracle. This result is particularly relevant because the "weak factor" case is common in applications. While standard two-way fixed effects often explain a large share of the variance in aggregated data, the incremental contribution of interactive fixed effects can be small, yet still larger than the noise and thus important for valid inference arkhangelsky2021synthetic. Our result, which allows the weakest factor to be only marginally stronger than the noise level ($n^{-1/2}$), is therefore especially appealing for empirical work. This contrasts with estimators designed for "strong" factors, which are typically assumed in the large panel data literature (e.g., bai2009panel), and avoids the complexities of bias-aware inference required when factors are "at the noise level" armstrong2022robust.
A key practical implication of Theorem (ref) is that it provides a direct justification for using standard inference procedures, as formalized in the following corollary.
The asymptotic equivalence established in Theorem (ref) allows us to analyze the efficiency of our Sequential SDiD estimator by studying the well-understood properties of the oracle OLS estimator. Because $\hat{\tau}_{a,k}^{SSDiD}$ is first-order equivalent to $\hat{\tau}_{a,k}^{OLS}$, it inherits its asymptotic efficiency properties.
The most direct efficiency result stems from the Gauss-Markov theorem. If the aggregate errors are homoskedastic and serially uncorrelated for each cohort (i.e., $\mathbb{V}[\boldsymbol{\epsilon}_a|n_a] = \frac{\sigma^2}{n\pi_a}\mathcal{I}_{T}$ for some $\sigma^2 > 0$), then the oracle estimator in (ref) is the Best Linear Unbiased Estimator (BLUE). This guarantee is analogous to the one established by borusyak2024revisiting for their imputation procedure under similar conditions. Our results therefore extend this type of efficiency guarantee to the Sequential SDiD estimator. To our knowledge, this provides the first formal efficiency result for an SDiD-type estimator and, more broadly, for any synthetic control-type procedure.
The connection to OLS also clarifies the limits of our estimator's efficiency. While OLS is BLUE under spherical errors, it is not, in general, the most efficient estimator in the presence of heteroskedasticity or serial correlation; Generalized Least Squares (GLS) would be preferred. While we do not formally establish this, we anticipate that under the i.i.d model (Example 1 in Section (ref)) the relevant limit experiment for estimating $\tau_{a,k}$ is a normal model with a mean structure given by (ref) but with a non-spherical error covariance. In such cases, the semiparametric efficiency bound would be attained by a GLS-type estimator. This implies that while our proposed estimator is efficient within a class that mirrors common empirical practice, we do not expect it to be semiparametrically efficient without further assumptions or modifications.
Our analysis thus far has abstracted from covariates. In practice, however, researchers usually observe both time-invariant characteristics, which we denote by $X_i$, and time-varying characteristics, $Z_{i,t}$. In this section, we discuss how to incorporate such information into our framework.
Our focus is on time-invariant covariates $X_i$ (assumed to be discrete), for two reasons. First, from a practical standpoint, interacting time fixed effects with time-invariant characteristics is a powerful and common way to control for unobserved heterogeneity. Time-varying controls $Z_{i,t}$, in our experience, often have limited additional predictive power in many applications.\footnote{This is the case in the empirical example we analyze in Section (ref).} Second, there is a theoretical challenge with time-varying covariates: if $Z_{i,t}$ has a substantial, dynamic causal effect on the outcome, it essentially behaves like an additional treatment variable. Properly analyzing such multi-treatment settings requires strong assumptions about effect homogeneity and is beyond the scope of this paper (see, e.g., de2023two).
We therefore focus on three distinct strategies for incorporating discrete, time-invariant covariates $X_i$ into our framework, each suited to a different empirical context.
\paragraph{Strategy 1: Full Stratification.} The most flexible approach is to allow all model parameters to vary with the value of the covariate, $x \in \mathcal{X}$. This corresponds to a fully-interacted model:
This approach is equivalent to stratifying the data by each value of $x$ and applying our main algorithm separately to each stratum. While conceptually straightforward, this is often impractical. If the set of covariates $\mathcal{X}$ is large, the resulting subsamples for each $(a,x)$ pair may be too small for our asymptotic arguments to hold, leading to noisy and unstable estimates. We therefore do not recommend this strategy unless every stratum is large.
\paragraph{Strategy 2: A Practical Hybrid Model (Recommended Approach).} Our recommended strategy for most applications is to adopt a more parsimonious, hybrid model. We allow the additive unit- and time-specific fixed effects to depend on $X_i$, but assume the multiplicative factors $\psi_t$ are common across all units:
where we assume $\mathbb{E}[\epsilon_{i,t}|A_i, X_i] = 0$. This specification is more restrictive than the fully stratified model in (ref), but it remains more general than the standard conditional parallel trends assumption made in much of the literature (e.g., abadie2005semiparametric,sant2020doubly). It is also particularly natural if we view $\psi_t$ as representing pervasive aggregate shocks, as in Example 2 of Section (ref).
The key advantage of model (ref) is that it behaves well under aggregation. A researcher can first compute cohort- and covariate-specific averages, $Y_{a,t}(x)$, and then average these over $x$ to produce a single aggregate series, $Y_{a,t}$, for each cohort $a$. The resulting series will have exactly the same structure as in our main model (ref), allowing for the direct application of Algorithm (ref). This approach is robust to arbitrary heterogeneity in the treatment effects across covariate groups and is our recommended strategy in the common case where adoption times vary within covariate groups (i.e., $A_i$ is not a deterministic function of $X_i$).
A special consideration arises for the never-treated cohort ($a=\infty$). Since we do not estimate treatment effects for this group, the outcomes $Y_{\infty, t}(x)$ for different $x$ can be aggregated using flexible, data-driven weights. In practice, if the number of never-treated units is large, we suggest including each $Y_{\infty, t}(x)$ as a separate control series in the synthetic control algorithm.
\paragraph{Strategy 3: The Case of Group-Level Treatment Assignment.} Finally, we consider a different approach for the specific but important case where adoption time is a deterministic function of the covariates, i.e., $A_i = A(X_i)$. This is common in applications where $X_i$ represents a geographic unit like a state or county. In this setting, the model with conditional parallel trends is too flexible to allow for identification, and researchers often rely on the stronger assumption of unconditional parallel trends. As an attractive alternative, our framework allows one to strengthen the exogeneity assumption using $X_i$ without including it directly in the functional form for the outcome:
In this case, rather than averaging outcomes within adoption cohorts, it is natural to average them within the groups defined by $X_i$. We illustrate this strategy in one of the simulations in the next section.
We first apply our method to reevaluate the findings of bailey2015war, who study the effect of Community Health Centers (CHCs) on mortality rates. The original dataset consists of county-year level observations from 1959-1988. We follow the authors' main sample construction choices and exclude the three most populous counties, which were all treated in the same year, to ensure a more balanced design. The main outcome, $Y_{i,t}$, is the adjusted mortality rate, and the treatment, $W_{i,t}$, is the presence of a CHC in a given county and year.
To implement our procedure, we first aggregate the data. For treated counties, we define cohorts based on their CHC adoption date ($A_i$). For never-treated counties, we create five distinct control cohorts based on their 1960 urban population percentage, a key baseline covariate. This follows the logic of our Strategy 2 from Section (ref). All cohort-level outcomes, $Y_{a,t}$, are calculated as weighted averages, using 1960 county populations as weights.
We then compute our Sequential SDiD estimates using Algorithm (ref), setting the regularization parameter $\eta^2 = \hat{\sigma}^2/n^{0.9}$, where $\hat{\sigma}^2$ is a preliminary estimate of the error variance from a standard two-way fixed effects model. For comparison, we also compute standard DiD estimates by setting $\eta = \infty$. Standard errors are produced via the Bayesian bootstrap with $1000$ replications.
Figure (ref) plots the resulting estimates. As shown, the Sequential SDiD estimates closely track the standard DiD estimates. This is expected, as the original study's pre-trends analysis suggests that the two-way fixed effects model is a good approximation for this application. We view this as a successful "proof of concept," demonstrating that our method produces sensible results in settings where simpler methods are believed to perform well.
We now use a simulation study to assess our estimator's performance when the parallel trends assumption fails due to unobserved interactive fixed effects. Our simulation design is based on the bailey2015war data. First, we use the matrix completion method of athey2021matrix on the original data to create a complete panel of potential outcomes under no treatment, $Y_{i,t}(\infty)$. From this completed panel, we estimate and extract the interactive fixed effects component, which we find to have a rank of 5.
We then generate simulated data where the strength of these interactive fixed effects varies. We define a "signal" strength parameter, where a 0% signal corresponds to a standard two-way fixed effects model (no interactive effects), and an 80% signal corresponds to a model where the variance of the interactive fixed effects component is four times larger than the variance of the idiosyncratic noise. To ensure our cohorts are large enough for the asymptotics to be relevant, we expand the sample size by replicating each unit's data (including its treatment status and factor loadings) four times.
For each signal level, we run 1000 simulations. In each simulation, we apply both our Sequential SDiD estimator and the standard DiD estimator, and we compute a $t$-statistic for the estimated treatment effect using 100 bootstrap replications to calculate the standard error.
Figure (ref) displays the distribution of these t-statistics for the contemporaneous treatment effect, $\tau_0$. In the 0% signal case (top row), where the DiD model is correctly specified, both estimators perform well, though the t-statistics for SSDiD are slightly more concentrated, suggesting our bootstrap standard errors are modestly conservative. In the 80% signal case (bottom row), the standard DiD estimator is severely biased; its t-statistics are centered far below zero, implying its confidence intervals would have near-zero coverage. In sharp contrast, the t-statistics for our Sequential SDiD estimator remain centered at zero, demonstrating that it provides valid inference even in the presence of strong confounding factors.
Figure (ref) repeats this exercise for the treatment effect four years after adoption, $\tau_4$. Again, our estimator performs well in both scenarios. The standard DiD estimator remains biased in the 80% signal case, although its bias is less extreme. This is likely because the estimation of longer-run effects is inherently noisier, making the bias a smaller component of the overall estimation error. Nonetheless, only the Sequential SDiD estimator provides reliable inference across all scenarios.
In our second experiment, we assess the performance of our estimator in a setting calibrated to real-world aggregated panel data where interactive fixed effects are empirically important. Inspired by the work of Bertrand2004did, the simulation design uses a state-by-year panel of average log wages for women, constructed from the March Current Population Survey (CPS). The panel consists of 50 states over 40 years. Crucially, this design simulates data at the aggregate (state) level, mirroring many empirical applications where researchers do not have access to the underlying microdata. Following the simulation design of arkhangelsky2021synthetic, we decompose the observed log-wage panel to estimate its structural components to use for the data-generating process. We find that a standard two-way fixed effects model explains about 94% of the variation, and an interactive fixed effects component with a rank of $r=4$ explains an additional 5%. For the simulation, we treat these estimated components as the true, fixed parameters and generate new data by drawing random shocks from the fitted AR(2) process that models the remaining idiosyncratic error.
A key feature of this experiment is that we induce a correlation between treatment timing and the interactive fixed effects, creating a direct violation of the parallel trends assumption required by standard DiD estimators. We implement a two-stage procedure to assign adoption dates. First, we designate states as "ever-treated" with a probability that depends on the state's factor loading on the first principal component of the interactive-effects model. Second, for these ever-treated states, we assign an adoption date drawn from a normal distribution. The mean of this distribution is also a function of the state's factor loading, ensuring that states with different unobserved trends are systematically treated at different times. The adoption period begins at $t=20$, and a subset of states is always left untreated. Figure (ref) plots the resulting empirical cumulative distribution function of the simulated adoption dates.
The results demonstrate that Sequential SDiD provides valid inference while the standard DiD estimator fails. Figure (ref) shows that the $t$-statistics for the standard DiD estimator are biased for both contemporaneous ($k=0$) and lagged ($k=4$) effects, whereas the $t$-statistics for our estimator remain correctly centered at zero. Table (ref) quantifies this finding: Sequential SDiD has a consistently lower Root Mean Squared Error (RMSE) and its 95% confidence intervals achieve the nominal coverage rate. In contrast, the coverage for the standard DiD estimator falls to around 70% due to its bias, rendering it unreliable for inference in this setting.
We propose a new method, Sequential Synthetic Difference-in-Differences (Sequential SDiD), for estimating treatment effects in event studies with a staggered rollout. Our estimator applies the principles of the original SDiD estimator sequentially to aggregated data, using an iterative imputation procedure where estimates for early cohorts inform those for later ones. We establish the estimator's theoretical properties by proving its asymptotic equivalence to an oracle OLS estimator. This result is significant as it delivers, to our knowledge, the first formal efficiency guarantees for a synthetic control-type method. An empirical application and data-based simulations demonstrate that our estimator is competitive with traditional DiD when its assumptions hold and provides robust inference when they fail.