Statistical Power for Estimating Treatment Effects Using Difference-in-Differences and Comparative Interrupted Time Series Estimators With Variation in Treatment Timing

Introduction

When randomized controlled trials (RCTs) are not feasible and time series data are available, panel data methods can be used to estimate treatment effects on outcomes, by exploiting variation in policies and conditions over time and across locations. A feature of panel data is that treatment timing often varies across the sample, for example, due to differences across locations in treatment implementation (such as school turnaround efforts; Redding & Nguyen, 2020), adoption of laws, or natural phenomena (such as the onset of Covid-19). Only recently has this timing feature been addressed in the literature for estimating treatment effects using difference-in-differences (DID) estimators but not for other panel data estimators (Athey & Imbens, 2018; Borusyak & Jaravel, 2017; Callaway & Sant’Anna, 2021; de Chaisemartin & D’Haultfœuille, 2020; Goodman-Bacon, 2018; Sun & Abraham, 2020). Another feature of panel data is that outcomes are often autocorrelated over time, which if not addressed can lead to estimated standard errors that are seriously biased downward (Bertrand et al., 2004).

This article discusses new power formulas to assess required sample sizes that account for both of these real-world features for two classes of commonly used panel data estimators: (1) DID estimators that contrast changes in outcomes during the pre- and posttreatment periods across the treated and comparison groups and (2) comparative interrupted time series (CITS) and associated ITS estimators that rely on fitted pre- and postperiod trend lines. We adopt common specifications found in the literature for these estimators and build on the much smaller associated literature on statistical power that has focused only on settings without variation in treatment timing. We consider binary treatments—where individuals or groups who become treated remain so during the following periods—and focus on both continuous outcomes (such as student test scores) and binary outcomes aggregated as rates to the cluster level.

Another contribution of this work is that we incorporate recent approaches for adjusting for the clustering of individuals within groups—such as geographic or educational units—an issue where there has been considerable confusion in the causal inference literature for quasi-experimental designs (QEDs; Abadie et al., 2017). We consider panel designs with separate cross sections of individuals as well as longitudinal designs where the same individuals are followed over time. Further, unlike previous power analysis studies, we allow for time periods to be unevenly spaced and for the inclusion of model covariates to adjust for potential confounding bias and to improve precision. Using an event history approach, we develop power formulas for both point-in-time effects (to examine treatment dynamics) and pooled effects averaged over the full postperiod. We use a potential outcomes framework (Holland, 1986; Rubin, 1974, 2005) to facilitate comparisons across the estimators. Our approach accommodates studies using either individual- or cluster-level (aggregate) data. A free Shiny R dashboard, Power_Panel available for download in the Supplementary Materials with a .pdf documentation file, performs the sample size calculations for the considered estimators.

The rest of this article is in eight sections. The second section discusses the related power analysis literature, and the third section discusses our statistical power framework. In the fourth section, we discuss clustering, and the fifth section presents notation and target estimands. The sixth and seventh sections present sample size formulas for the DID and CITS/ITS estimators, where our main results are presented in a series of theorems. In the eighth section, we present an illustrative power analysis for the considered estimators and then conclude. For reference, Table 1 lists the equation numbers for the variance formulas presented in the text for each estimator. Online Appendix A provides more details on the power formulas.

OPEN IN VIEWER

Table 1.

Equation Numbers of Variance Formulas in the Text, by Estimator

Panel Data Estimator (Relative to the Average Preperiod and Averaged Across All Timing Groups)	Equation Number in Text
	Variance Formula	Regression Model for Estimation
Difference-in-differences (DID) estimator
Pooled (average) effect in the postperiod
Cross-sectional	Equation 14	Equation 7
Longitudinal	Online Appendix A.1	Equation 7, including autocorrelations for individual-level errors
Point-in-time effect after a certain amount of exposure to the treatment or at a specific postperiod calendar time
Cross-sectional	Equation 15	Equation 7
Longitudinal	Online Appendix A.1	Equation 7, including autocorrelations for individual-level errors
Models with covariates	Equation 18	Equation 7 with covariates
Comparative interrupted time series (CITS) estimators (ITS estimators omit the comparison group terms)
Fully interacted estimator (allows intercepts and linear slopes to differ for the pre- and postperiods)a
Pooled (average) effect in the postperioda
Cross-sectional	Equation 24	Equation 19
Longitudinal	Online Appendix A.1	Equation 19, including autocorrelations for individual-level errors
Point-in-time effect after a certain amount of exposure to the treatment or at a specific postperiod calendar timea
Cross-sectional	Equation A.2	Equation 19
Longitudinal	Online Appendix A.1	Equation 19, including autocorrelations for individual-level errors
Common-slopes estimator (assumes the same linear slopes for the pre- and postperiods)
Pooled and point-in-time effects are the same
Cross-sectional	Equation A.3	Equation 19, omitting preperiod slope terms
Longitudinal	Online Appendix A.1	Equation 19, omitting preperiod slope terms and including autocorrelations for individual-level errors

^a The discrete estimator is similar to the fully interacted estimator, except that the postperiod in Equation 19 is modeled using discrete postperiod indicators rather than linear trendlines. The pooled variance estimators are the same for the discrete and fully interacted estimators. Variances for the point-in-time discrete estimators can be obtained using Equation 24 and Note 6 in the main text.

Related Literature

We extend the previous power analysis literature for DID estimators to simultaneously allow for variation in treatment timing (not previously addressed), autocorrelated errors, clustering, both average and point-in-time effects, and uneven time intervals and covariates (both not previously addressed). Frison and Pocock (1992) provide sample size formulas for nonclustered, longitudinal DID designs with constant correlations across time but do not consider clustered designs. McKenzie (2012) considers similar nonclustered designs with more general error structures but does not provide explicit variance formulas. While Somers et al. (2013) consider clustered designs, they examine two time periods only and do not address staggered timing or autocorrelation. Burlig et al. (2019) provide sample size formulas for cross-sectional models with both clustering and autocorrelated errors for a more restrictive model specification than ours, but do not allow for staggered timing, uneven time intervals, covariates, or point-in-time effects.

To our knowledge, Bloom (1999) is the only article that has developed statistical power formulas for CITS and ITS estimators. However, this study only considers models with discrete postperiod indicators, but not more common specifications that model both pre- and postperiod trendlines that are considered here. Further, Bloom (1999) considers a single follow-up period only and does not address staggered timing, autocorrelation, or covariates. Zhang et al. (2011), Kontopantelis (2018), Hawley et al. (2019), and Liu et al. (2019) use simulations to calculate power for the ITS design but do not develop power formulas or address staggered timing. Several articles conduct within-study comparisons using RCT data to empirically compare RCT impact estimates with those obtained using CITS or ITS methods (Baicker & Svoronos, 2019; Clair et al., 2016) but do not develop variance formulas. In the theory sections, we show how our general formulas reduce to those from the literature.

Statistical Power Framework

Our power analysis relies on formulas for minimum detectable impacts, which represent the smallest true impacts that can be detected with a high probability. We scale the minimal detectable impacts into effect size (standard deviation) units—hereafter referred to as MDEs—which is common in policy research, especially for outcomes that are difficult to interpret in nominal units (although this scaling is not necessary). This approach parallels the statistical power literature for RCTs (see, e.g., Bloom, 1995; Donner & Klar, 2000; Murray, 1998; Schochet, 2008) and regression discontinuity designs (Schochet, 2009).

All the considered impact estimators are multiple regression estimators. Accordingly, relying on asymptotic normality and a classical hypothesis testing approach, we assume the use of t-statistics to test null hypotheses of the form, H 0 : β = 0 , where β is the impact parameter of interest. This approach yields the following MDE formula for a two-tailed hypothesis test:

MDE = 1 σ Factor α , λ , d f Var β ^ ,

1

where Factor α , λ , d f = T − 1 α / 2 , d f + T − 1 λ , d f , T − 1 is the inverse of the student’s t distribution function, α is the significance level, λ is the statistical power level, d f is the degrees of freedom, β ^ is the impact estimator, Var β ^ is the variance of β ^ based on the specific design features, and σ is the standard deviation of the outcome variable at a particular time point.

Alternatively, one can solve Equation 1 to calculate required sample sizes or the number of time periods to achieve a given MDE value, assuming specific values of α , λ , and the design parameters in Var β ^ . Study sample sizes enter Equation 1 primarily through Var β ^ but also through Factor α , λ , d f because of d f . Hill et al. (2008) and Lipsey et al. (2012) discuss a framework for selecting MDE targets for education evaluations, for example, by examining the natural growth in outcomes over time, policy-relevant gaps across subgroups, and observed effect sizes from previous similar evaluations.

The power formula in Equation 1 hinges critically on Var β ^ . Thus, this article focuses on calculating closed-formed expressions for Var β ^ for the considered designs. However, we begin by discussing the general issue of clustering that can have a large effect on precision.

Clustering

Clustering occurs when the outcomes of individuals (or broader units) in the data set are correlated. Consider a panel analysis using data on the same individuals (or units) over time. In this case, outcomes are likely to be correlated for the same individual (or unit) across time periods, and our variance formulas incorporate this form of clustering. However, there has been much confusion in the QED causal inference literature about how to account for clustering across different individuals (or units) in the sample (Abadie et al., 2017).

In this article, we adopt an approach developed for RCTs and inspired by design-based theory to account for clustering across sample members or units (Abadie et al., 2017; Freedman, 2008; Imbens & Rubin, 2015; Schochet, 2010, 2013, 2016, 2020). Under this approach, one must be explicit about what is random under repeated sampling. For instance, if the study units are assumed to be fixed (which occurs in finite population settings), the main source of randomness is the treatment assignment itself. Accordingly, in this setting, the presence of clustering hinges on the considered unit of treatment assignment.

To determine this unit, it is useful to consider what would be the parallel unit under a hypothetical RCT. To help fix concepts, consider a panel data study examining whether Covid-19 had a larger effect on students in states where the pandemic hit earlier than in states where it hit later. In this case, the hypothetical unit of treatment assignment under an RCT would be the state, and thus, similarly in the panel data context. However, suppose the study instead examined the effects of pandemic-related educational policies (such as distance learning) that differed by school district. In this case, we would treat districts as the unit of treatment assignment and the source of clustering.

Note that in some designs, study units are formally sampled for the study from broader populations (e.g., studies using national data sets with multistage sample designs or multisite random block designs) or are deemed to have been so. In these cases, study results are assumed to generalize beyond the study sample, and clustering from the sampling of units becomes pertinent (Schochet, 2008). We do not formally consider such forms of clustering because their variance estimators differ from those developed here due to the emergence of correlations between the outcomes of the treatment and comparison groups within the same higher level clusters. However, our variance formulas are likely to be conservative in these cases (because they ignore these correlations). Further, our approach applies to settings where the higher level clusters are treated as fixed effects rather than as random blocks.

Clustering increases the MDEs in Equation 1 for two interrelated reasons (Donner & Klar, 2000; Hedges, 2007; Murray, 1998; Schochet, 2008). First, design effects due to correlated outcomes reduce effective sample sizes, and hence precision, as can be quantified by intraclass correlation coefficients (ICCs). Second, d f s are based on the number of clusters, not the number of individuals, which increases Factor α , λ , d f values in Equation 1.

While many designs in our setting are likely to be clustered, there may be longitudinal designs where the unit of treatment assignment is the individual. For instance, consider a study that aims to examine the effects of tutoring on the achievement of low-performing students who have access to tutors at different time points. In this case, the unit of treatment assignment could be the student. While we consider more general clustered designs, our variance formulas also apply (reduce) to nonclustered, longitudinal designs by setting the pertinent ICC parameters to 0 and treating clusters as individuals (as discussed further in Section “DID Estimator”).

Definitions, Assumptions, and Target Estimands

We consider a clustered panel data setting with M total clusters defined by the unit of treatment assignment. We assume M_T clusters in the treatment group and M C untreated clusters in the comparison group. We assume the same clusters remain in the sample for each of P time periods (e.g., where time is measured in months, semesters, or years). We denote time periods by t = 1 , … , P with both pre- and posttreatment periods. As discussed later, the time periods do not necessarily have to be evenly spaced, but we avoid this notational inconvenience for now.

We assume K treatment timing groups in total, with treatment start times of S k for timing group k ∈ 1 , … , K ordered from earliest to most recent, with M T k treatment clusters in group k. We assume M C k comparison clusters are matched to timing group k ( M T k and M C k can differ). We do not consider an analysis where each timing group is compared to the full comparison group (as is sometimes done in the literature), because under our event history framework, this approach considerably complicates the variance formulas due to covariance terms that arise across the timing group estimators due to the common comparison group. It can be shown, however, that the power formulas developed here are good approximations for panel designs with an overlapping comparison group sample (e.g., if this sample is proportionally split across the timing groups). Our stratified approach may also be desirable in practice, especially if there are differences in cluster-level characteristics across the timing groups (see Daw & Hatfield, 2018, for a discussion of matching issues in panel studies).

We assume multiple time periods and the availability of both pre- and postperiod data ( S k ≥ 2 for the DID design and 4 ≤ S k ≤ P − 2 for the CITS and ITS designs for all k). We assume treatment effects are observed and measured starting in period S_k , which pertains to the first posttreatment period. Let T_j be the treatment indicator that equals 1 for ever-treated clusters and 0 for comparison (never-treated) clusters. Further, let G T j = k denote the group of treatment clusters in timing group k (with onset time S_k ), where G C j = k are the group of matched comparison clusters. Finally, let A k = P − S k + 1 denote the number of postintervention (“after”) periods, where B k = ( S k − 1 ) are the number of preintervention (“before”) periods. For reference, Table 2 provides key parameter definitions used throughout this article.

OPEN IN VIEWER

Table 2.

Key Notation and Definitions

Input Variable	Definition	Variable Range
M , M T , M C , r	Number of total, treatment, and comparison clusters, where r = M T / M is the share of treatment clusters	M T ≥ 2 ; M C ≥ 2 for DID/CITS designs and M C = 0 for ITS design
P, t	Total number of time periods, where time is indexed by t	P ≥ 2 for DID design; P ≥ 6 for CITS and ITS designs
N	Number of individuals per cluster per time period	N ≥ 1
Tj	Treatment indicator that equals 1 for ever-treated clusters and 0 for comparison (never-treated) clusters indexed by j
Time t ; t ∈ 1 , … , P	Time of outcome measurement, in elapsed calendar time since a common reference point (e.g., months or quarters). These time points could be evenly or unevenly spaced	Time t > 0
K	Number of treatment timing groups in the treatment sample	K ≥ 1
Sk ; k ∈ 1 , … , K	Treatment start period for each timing group, measured using the time label, t, not the elapsed calendar time variable, Time t . Treatment effects can be measured starting in Sk	2 ≤ S k ≤ P for DID design; 4 ≤ S k ≤ P − 2 for CITS/ITS designs
A k = P − S k + 1 ; B k = S k − 1	Number of post- and preintervention periods	Ak , B k ≥ 1 for DID design; A k , B k ≥ 3 for CITS/ITS designs
Pre k , Post k	1/0 indicators of preperiod ( t < S k ) or postperiod t ≥ S k	Pre k , Post k = 0 or 1
M T k , M C k , p T k , p C k	Number of treatment and associated comparison clusters in timing group k, where p T k = M T k / M T and p C k = M C k / M C	M T k , M C k > 0 ; ∑ k M T k = M T ; ∑ k M C k = M C ; ∑ k p . k = 1
G T j , G C j , Gj	G T j = k for treatment clusters in timing group k, G C j = k for matched comparison clusters, and G j = k for both groups	G T j , G C j , G j ∈ 1 , … , K
Y i j t 1 , Y i j t 0	Potential outcomes in the treated and untreated conditions for individual i in cluster j at time t	Continuous or binary
Y ¯ . t = q ; Y ¯ ¯ . t …	Mean potential outcomes averaged over all clusters in a timing group at time q ( Y ¯ . ) or over a specific period Y ¯ ¯ .	Means or proportions
y i j t ; y ¯ t = q G . j = k ; y ¯¯ ¯ t … G . j = k	Observed outcomes and cluster-level means at a time point or period, averaged over a timing group and treatment status	Continuous or binary
w k q , w k	Weight assigned to timing group k in postperiod q ( w k q ) or to timing group k when calculating pooled effects (wk )	w k q > 0 ; w k = ∑ q = S k P w k q
ρ	Correlation of cluster-level outcomes over time, modeled using an AR(1) structure (or assumed constant over time)	ρ < 1
ρ ¯ k Pre , ρ ¯ k Post , …	Mean correlations over specific time periods	ρ ¯ < 1
ICC θ	Intraclass correlation coefficient measuring the percentage of total variance in the outcome that is due to variation between cluster-time cells	0 ≤ ICC θ ≤ 1
ψ ; ψ ¯ k Pre , ψ ¯ k Post , …	Correlation of outcomes for the same individual over time for longitudinal designs, modeled using an AR(1) structure	| ψ | < 1
R Y X 2 ; x i j t	Regression R 2 value from the inclusion of a vector of covariates, x i j t , conditional on the other model parameters	0 ≤ R Y X 2 < 1 ; x i j t can be continuous or binary
R T X 2	Average R 2 value from regressing the treatment status interaction terms on the model covariates (treatment-covariate collinearity)	0 ≤ R T X 2 < 1

Note. DID = difference-in-differences; CITS = comparative interrupted time series.

We assume data are available for individuals in each study cluster at each time period, which allows us to fully examine sources of variation in the model error terms. However, the power formulas apply also to analyses using data averaged to the cluster level or to lower level units within the study clusters. Data aggregation (averaging) aligns with our framework because we consider linear panel data estimators, so the regression model structure is not affected by aggregation (Greene, 2012; Schochet, 2020). The analysis, however, cannot be conducted if data are available only for units at higher levels than the unit of treatment assignment.

We index individuals by i j t for individual i in cluster j in time t. For each cluster, either all individuals are treated or not. For simplicity and as is customary in the power literature, we assume a balanced design with N j t = N individuals per cluster in each time period; for unbalanced designs, the average sample size across all time periods and clusters, N ¯ ¯ = 1 P M ∑ t = 1 P ∑ j = 1 M N j t , can be used in the power formulas as an approximation (Donner & Klar, 2000; Kish, 1995; Murray, 1998; Schochet, 2008). This assumption does not materially affect the analysis as power for clustered designs is determined primarily by the number of clusters, not the number of individuals per cluster (unless this number is small). We allow for both separate cross sections of individuals (e.g., separate cohorts of fourth graders in each time period) and longitudinal data where the same persons are followed over time.

To make the large amount of notation concrete, consider the panel study conducted by Zimmer et al. (2017) that estimated the effects of school turnaround strategies in low-performing Tennessee schools on student test scores. In this study, there were P = 8 time periods with M T = 51 schools (clusters) in the treatment group and M C = 28 comparison clusters. There was variation in treatment timing, with 19 treatment schools treated at time point t = 6 , 20 schools treated at t = 7 , and 12 schools treated at t = 8 . Thus, the study had three timing groups ( k = 1 , 2, and 3) with associated treatment start times of S 1 = 6 , S 2 = 7 , and S 3 = 8 and treatment cluster sample sizes of M T 1 = 19 in timing group 1 ( G T j = 1 ), M T 2 = 20 in timing group 2 ( G T j = 2 ), and M T 3 = 12 in timing group 3 ( G T j = 3 ). The lengths of the pre- and postperiods were B 1 = 5 and A 1 = 3 for k = 1 , B 2 = 6 and A 2 = 2 for k = 2 , and B 3 = 7 and A 3 = 1 for k = 3 . Finally, there was an average of N = 230 students per school per time period. We use this running example throughout this article.

We invoke several assumptions that apply to both the DID and CITS/ITS estimators to obtain unbiased estimators of well-defined treatment effect parameters (estimands), where additional assumptions specific to each design are presented in Sections “DID Estimator” and “CITS and ITS Estimators”:

Assumption 1. The stable unit treatment value assumption (SUTVA; Rubin, 1986 ). Let Y i j T clus denote the potential outcome for an individual given the random vector of all cluster treatment assignments, T clus . Then, if T j = T j ' for cluster j, we have that Y i j T clus = Y i j T clus ' .

SUTVA allows us to express Y i j T clus as Y i j T j , so that the potential outcomes of an individual in cluster j depend only on the cluster’s treatment assignment and not on the treatment assignments of other clusters in the sample. More specifically, SUTVA allows us to define Y i j t 1 as the potential outcome for the individual in the treatment condition and Y i j t 0 as the potential outcome in the nontreated condition. Potential outcomes are assumed to be continuous variables, although our methods also approximately apply to binary outcomes that are aggregated as rates or proportions to the cluster level. Using this framework, the data generating process for the observed outcome at each time point, y i j t , can be expressed as y i j t = T j Y i j t 1 + 1 − T j Y i j t 0 . SUTVA also requires that units cannot receive different forms of the treatment.

Assumption 2. No anticipatory behavior. We assume that a future treatment does not affect past outcomes, so that y i j t = Y i j t 1 = Y i j t 0 for t < S j . This assumption is required to rule out situations where individuals or units change their behavior in anticipation of a treatment that is likely to occur in the future, which could influence preperiod outcomes.

To define the focal estimands for our analysis, we adopt an “event history” approach that aggregates treatment effect parameters across particular treatment timing groups and time periods (Callaway & Sant’Anna, 2021; Sun & Abraham, 2020). With variation in treatment timing, an important consideration for defining our focal estimands is whether to measure the postprogram impacts in calendar time or relative to exposure to the treatment. Consider first a calendar-time analysis, where we define ATT k q as the average effect of treatment on the treated (ATT) for clusters in timing group G T j = k in postperiod q (assuming q ≥ S k ):

ATT k q = E Y ¯ 1 t = q − Y ¯ 0 t = q | G T j = k .

2

Here, E ( Y ¯ 1 t = q | G T j = k ) = 1 M T k N ∑ j : G T j = k M T k ∑ i = 1 N E Y i j q 1 is the mean cluster-level outcome in the treated condition at time q for timing group k, and E ( Y ¯ 0 t = q | G T j = k ) is the mean outcome if these clusters had not been treated. Following Callaway and Sant’Anna (2021), our first impact parameter of interest averages the ATT k q estimands across timing groups:

A T T q = 1 ∑ k = 1 K w k q I q ≥ S k ∑ k = 1 K w k q I q ≥ S k ATT k q ,

3

where I q ≥ S k is an indicator that equals 1 if q is a postperiod for timing group k and 0 otherwise, and w k q are weights (which we set to 1, as discussed in Section “Cross-Sectional Analysis: Framework”). Some timing groups may not contribute to Equation 3. For instance, in our running example, q = 7 is a postperiod for timing groups k = 1 and k = 2 with S 1 = 6 and S 2 = 7 , but not for timing group k = 3 with S 3 = 8 , so timing group 3 would be excluded from Equation 3 to calculate ATT 7 .

Note that at calendar time q, exposure to the intervention, q − S k + 1 , will differ across timing groups, which could complicate the interpretation of the ATT q estimands if impacts vary over time. Using our running example, at q = 8 , the exposure time is three periods for timing group 1, compared to only two periods for timing group 2 and one period for timing group 3. Thus, an alternative approach used by Sun and Abraham (2020) is to realign the data to measure postperiod impacts relative to treatment exposure, where exposure time, t_e , can be mapped to calendar time using t e = t − S k + 1 = t − B k . This approach yields the following ATT l e estimand pertaining to impacts at exposure point l that are averaged over the timing groups:

ATT l e = 1 ∑ k = 1 K w k l e I l ≤ A k ∑ k = 1 K w k l e I l ≤ A k ATT k l e ,

4

where ATT k l e = E Y ¯ 1 t e = l − Y ¯ 0 t e = l | G T j = k is the ATT parameter for timing group k at exposure point l; I l ≤ A k is an indicator that equals 1 if the length of the postperiod (A_k ) is at least l periods and 0 otherwise; and w k l e are weights (which we set to 1). These point-in-time ATT l e estimands are useful for examining treatment dynamics (Sun & Abraham, 2020), although potential differences in impacts across timing groups could affect interpretation.

Power concerns are similar for the point-in-time ATT q and ATT l e estimands as they both rely on the same ATT k q components, differing in only how these components are organized and averaged. Thus, our analysis applies to both impact parameters. While we prefer the ATT l e estimand (since it has a clearer interpretation), for ease of exposition, we use calendar time notation because it conforms to the format of time series data that are typically available for analysis; we then transform the results into exposure time as needed. The choice of notation does not change the results.

Another key estimand for our analysis, ATT POOL , takes a weighted average of either the ATT q or ATT l e parameters over their respective postperiods to obtain pooled treatment effects averaged over the full observed postperiod (both approaches yield the same result). This pooled estimand has policy relevance as a summary measure of treatment effects, for instance, for use in benefit–cost analyses. An intuitive way to express the ATT POOL estimand is as follows:

ATT POOL = 1 ∑ k ∑ q w k q ∑ k = 1 K ∑ q = S k P w k q A T T k q ,

5

where we set w k q = 1 for our analysis. ¹ This aggregate parameter is flexible in that it allows for heterogenous treatment effects both across timing groups and across time within timing groups.

The estimation challenge for the considered treatment effect parameters is that we do not observe counterfactual outcomes for the treated groups during the posttreatment period, that is, the E Y ¯ 0 t = q | G T j = k terms in Equation 2. Thus, our considered DID and CITS/ITS panel methods estimate these counterfactuals under various identification assumptions, which then allows for consistent estimation of the ATT q , ATT l e , and ATT POOL estimands in Equations 3 –5.

DID Estimator

DID methods identify causal effects from panel data by contrasting changes in outcomes during the pre- and posttreatment periods across the treatment and comparison groups. The large literature on DID methods focuses on designs without variation in treatment timing (e.g., Angrist & Pischke, 2009; Ashenfelter, 1978; Ashenfelter & Card, 1985; Bertrand et al., 2004). However, a smaller recent literature (cited in the Introduction section) considers impact estimation with staggered treatment timing to produce unbiased estimators. This literature underlies our analysis.

The key identifying assumption for DID methods is that in the absence of treatment, the mean outcomes for the treatment and comparison groups would exhibit parallel trends over time (Abadie, 2005; Heckman et al., 1997):

E Y ¯ 0 t = q − Y ¯ ¯ 0 t < S k | G T j = k = E Y ¯ 0 t = q − Y ¯ ¯ 0 t < S k | G C j = k .

6

Here, E ( Y ¯ ¯ 0 t < S k | G T j = k ) = 1 B k ∑ b = 1 B k E ( Y ¯ 0 t = b | G T j = k ) is the mean counterfactual outcome averaged over all preperiods and all treatment clusters in timing group k (using equal preperiod weights), and similarly for E ( Y ¯ ¯ 0 t < S k | G C j = k ) for the comparison clusters.

Intuitively, this assumption allows us to obtain an unbiased estimate of the unobserved counterfactual, E Y ¯ 0 t = q | G T j = k , using y ¯ ¯ t < S k G T j = k + y ¯ t = q G C j = k − y ¯ ¯ t < S k G C j = k , where y ¯ ¯ t < S k G T j = k = 1 B k M T k N ∑ t = 1 B k ∑ j : G T j = k M T k ∑ i = 1 N y i j t is the observed preperiod mean outcome for treatments in timing group k, and similarly for y ¯ t = q G C j = k and y ¯ ¯ t < S k G C j = k ; this yields the DID estimator discussed below.

In what follows, we develop variance formulas for the DID estimator, first for the cross-sectional design and then for the longitudinal design.

Cross-Sectional Analysis: Framework

To consider DID impact and associated variance estimators for the ATT q , ATT l e , and ATT POOL estimands, we rely on the following regression model, adapted from the literature to our context, using stacked data on separate cross sections of individuals nested within study units (clusters) over time:

y i j t = α j + ∑ k = 1 K δ k t I G j = k + ∑ k = 1 K ∑ q = S k P β k q I G j = k F i j t , q T j + θ j t + ε i j t ,

θ j t = ρθ j t − 1 + η j t .

7

In this model, α j and δ k t are cluster and time fixed effects, I G j = k is an indicator that equals 1 if either G T j = k or G C j = k and 0 otherwise, and F i j t , q is a period q indicator that equals 1 if t = q and 0 otherwise. The random error terms are assumed to have mean zero and to be distributed independently of each other: θ j t captures the correlations of individuals within the same cluster and time period, and ε i j t are iid 0,σ ε 2 individual-level errors. Following the influential work of Bertrand et al. (2004), we allow θ j t to be correlated over time using an autoregressive process of order 1 (AR(1)), where ρ < 1 is the autocorrelation parameter and η j t are iid 0 , σ η 2 errors with E η j t θ j t − 1 = 0 . For our analysis, we assume the same autocorrelation structure during the pre- and posttreatment periods and for each timing group and their associated comparison group. Under the AR(1) model with long panels, σ θ 2 = Var θ j t = σ η 2 / 1 − ρ 2 , which we assume hereafter. ² This is not a restrictive assumption as Var θ j t converges quickly over time (e.g., for ρ = 0.5 , the variance stabilizes after three periods).

Under the AR(1) model, correlations are larger for cluster observations closer in time than further apart, so Corr θ j t , θ j t − p = ρ p for p = 1 , 2 , … , t − 1 . Accounting for autocorrelated errors in many settings is important because cluster-level outcomes are typically highly correlated over time, which if ignored can lead to serious underestimates of standard errors (Bertrand et al., 2004). For instance, pertinent to studies in education, when estimating Equation 7 using National Assessment of Education Progress (NAEP) public use data on 28 school districts over 9 nine years (excluding the interaction terms), we find a ρ value of 0.49 for fourth grade reading scores and 0.46 for fourth grade math scores (both are statistically significant at the 1% level). Thus, the autocorrelations are substantial, even after controlling for time fixed effects.

In DID settings, time intervals are often equally spaced, but not always. Thus, we consider general settings where time measurements can be unevenly spaced within timing groups but not across them. We follow Baltagi and Wu (1999) who discuss ordinary least squares (OLS) estimation with unequal time intervals that fully maintains the AR(1) structure. ³ Under this approach, correlations between successive observations are based on their time differences (e.g., the correlation is ρ 3 for successive observations three time units apart). To use this approach, for each t, we use the notation, Time t , to denote the elapsed time between outcome measurement and a common reference point; thus, t refers to a time label (counter), whereas Time t refers to elapsed calendar time.

The regression model in Equation 7 includes three-way interactions between indicators of timing group, time period, and treatment status; the model, however, excludes interactions for each comparison group and preperiod. Thus, the resulting OLS estimators, β ^ DID , k q , provide DID estimates for each posttreatment time point relative to the mean preperiod. Formally, the DID estimator for timing group k in time period q ≥ S k is β ^ DID , k q = y ¯ t = q G T j = k − y ¯ ¯ { t < S k } G T j = k − y ¯ t = q G C j = k − y ¯ ¯ t < S k G C j = k , where y ¯ and y ¯ ¯ are mean observed outcomes defined above. For instance, using our running example, the DID estimator at time q = 7 for timing group k = 1 with S 1 = 6 can be calculated using β ^ DID , 17 = y ¯ t = 7 G T j = 1 − y ¯ ¯ { t < 6 } G T j = 1 − y ¯ t = 7 G C j = 1 − y ¯ ¯ t < 6 G C j = 1 . However, the DID estimator, β ^ DID , 37 , would not be germane for timing group k = 3 with S 3 = 8 , because q = 7 is not a postperiod for these clusters. Note that the OLS estimates using Equation 7 are the same as those from models run separately by timing group. ⁴

We can now aggregate the β ^ DID , k q estimators across timing groups to obtain an unbiased estimator for the calendar-time ATT q estimand in Equation 3 using

β ^ DID , q = 1 ∑ k = 1 K w k q I q ≥ S k ∑ k = 1 K w k q I q ≥ S k β ^ DID , k q .

8

Similarly, we can estimate the exposure-specific ATT l e estimand in Equation 4 using

β ^ DID , l e = 1 ∑ k = 1 K w k l e I l ≤ A k ∑ k = 1 K w k l e I l ≤ A k β ^ DID , k q e ,

9

where q e = l + B k converts exposure time to calendar time to select the pertinent β ^ DID , k q e estimators. Finally, we can aggregate β ^ DID , q or β ^ DID , l e across their postperiods to obtain an unbiased estimator for the ATT POOL estimand in Equation 5 using

β ^ DID,POOL = 1 ∑ k ∑ q w k q ∑ k = 1 K ∑ q = S k P w k q β ^ DID , k q ,

10

or using the equivalent expressions in Note 1.

For weighting, we follow Sun and Abraham (2020) who use w k q = w k l e = 1 , so that each included timing group receives the same weight in calculating the point-in-time estimators, β ^ DID , k q . This weighting scheme also implies that when calculating β ^ DID,POOL , the weight for each timing group is A_k , the length of its postperiod (i.e., w k = ∑ q = S k P w k q = A k ). For instance, in our running example, timing group 1 with A 1 = 3 would receive a larger weight than both timing group 2 with A 2 = 2 and timing group 3 with A 3 = 1 , because it has more exposure periods. In this setting, β ^ DID,POOL can be interpreted as the average treatment effect observed in the sample. Other options exist for w k q and w k l e , such as weighting according to the number of clusters ( M T k and M C k ; Callaway & Sant’Anna, 2020) or inversely proportional to timing group variances.

Finally, we note that our event history approach differs from the approach of Goodman-Bacon (2018) who examines impact estimation for a variant of model (7) where the three-way interaction terms are replaced by β D j t , where D j t is an indicator that equals 1 if cluster j in the treatment group was treated at or before time t and 0 otherwise and where a common comparison group is used across timing groups. Under this specification, β is the average ATT parameter over the postperiod. With staggered treatment timing, Goodman-Bacon (2018) shows that the OLS estimator, β ^ , becomes a weighted average of separate two-by-two DID estimators, where each treatment timing group is compared not only to the common comparison group but also to each other based on treatment timing. More specifically, a later timing group serves as a comparison group for the early timing group before its treatment begins and the early group then serves as a comparison group for the later timing group after its treatment begins.

We do not adopt this specification, however, because it only recovers the ATT estimand when treatment effects are homogeneous across timing groups and over time (Callaway & Sant’Anna, 2021; de Chaisemartin & D’Haultfœuille, 2020; Sun & Abraham, 2020). These assumptions, however, are unrealistic in practice. Instead, our event history approach allows for both sources of heterogeneity and yields unbiased and more interpretable impact estimates.

Cross-Sectional Analysis: Variance Estimation

In this section, we first obtain the variance formula for the pooled DID estimator, β ^ DID,POOL in Equation 10 and then discuss how this formula reduces to the variance formulas for the point-in-time estimators, β ^ DID , q in Equation 8 and β ^ DID , l e in Equation 9.

To obtain the variance of β ^ DID,POOL , it is useful to express β ^ DID,POOL in Equation 10 as an average of treatment effects across each timing group, β ^ DID , k , as follows:

β ^ DID,POOL = 1 ∑ k = 1 K A k ∑ k = 1 K A k β ^ DID , k ,

11

where β ^ DID , k = 1 A k ∑ q = S k P β ^ DID , k q , recalling that our weighting scheme uses w k q = 1 and w k = ∑ q = S k P w k q = A k . Because of the assumed independence of the outcomes across the timing groups, the variance of β ^ DID,POOL can then be obtained using the following simple relation:

Var β ^ DID,POOL = 1 ( ∑ k = 1 K A k ) 2 ∑ k = 1 K A k 2 Var β ^ DID , k .

12

We can now calculate Var β ^ DID,POOL by first calculating Var β ^ DID , k using the variances and covariances across the pre- and postperiod means based on the model error structure in Equation 7.

To fix concepts, we first consider the variance formula for timing group k without the AR(1) structure, which can be expressed as follows:

Var β ^ DID , k = 1 M T k + 1 M C k 1 A k + 1 B k σ θ 2 + σ ε 2 N .

13

This variance is intuitive in that it gets smaller as the number of clusters ( M T k and M C k ), the number of individuals per cluster (N), and the lengths of the pre- and postintervention periods (B_k and A_k ) increase. Clustering effects arise due to σ θ 2 . If σ θ 2 = 0 , N = 1 , and M T k and M C k are redefined to be the number of individuals rather than clusters, then Equation 13 reduces to the nonclustered design.

If we now add the AR(1) structure, the calculations become considerably more complex, because the mean outcomes for a particular cluster become correlated across time periods (both within the pre- and postperiods as well as across them), where the correlations diminish over time. This leads to our first main variance result provided in the following new theorem.

Theorem 1. The variance of the pooled DID estimator, β ^ DID,POOL , in Equation 12 obtained from the model in Equation 7 that incorporates clustering, variation in intervention timing, the AR(1) error structure, and general measurement intervals is as follows:

Var β ^ DID,POOL = 1 ( ∑ k = 1 K A k ) 2 ∑ k = 1 K A k 2 1 M T k + 1 M C k σ θ 2 1 A k + 1 B k + A k − 1 A k ρ ¯ k Post + B k − 1 B k ρ ¯ k Pre − 2 ρ ¯ k Pre-Post + σ ε 2 N 1 A k + 1 B k .

14

Here, ρ ¯ k Post = 2 A k A k − 1 ∑ a = S k P ∑ a ′ > a P ρ Time a ' − Time a denotes the average autocorrelation in the postperiod, ρ ¯ k Pre = 2 B k B k − 1 ∑ b = 1 B k ∑ b ′ > b B k ρ Time b ' − Time b denotes the average autocorrelation in the preperiod, ρ ¯ k Pre − Post = 1 B k A k ∑ b = 1 B k ∑ a = S k P ρ Time a − Time b denotes the average autocorrelation between the pre- and postperiods, and other terms are defined above.

In this expression, if ρ > 0 (the typical case), we find that the ρ ¯ k Post and ρ ¯ k Pre terms increase variance whereas the ρ ¯ k Pre - Post term reduces variance. In most cases, the ρ ¯ k Post and ρ ¯ k Pre terms will be larger than the ρ ¯ k Pre - Post term as cluster observations are, on average, closer in time within a pre- or postperiod than across them. Thus, accounting for the AR(1) structure tends to reduce precision.

However, there are several cases where accounting for the autocorrelated errors improves precision. First, this will occur if we assume constant autocorrelations across time as is sometimes specified in the power literature for shorter panels (Frison & Pocock, 1992; McKenzie, 2012). In this case, Corr θ j t , θ j t − p = ρ , and the Var β ^ DID , k formula reduces to 1 M T k + 1 M C k 1 A k + 1 B k σ θ 2 1 − ρ + σ ε 2 N . Second, in short panels with A k = 1 and B k = 1 , the Var β ^ DID , k formula reduces to 1 M T k + 1 M C k 2 σ θ 2 1 − ρ ¯ k Pre − Post + 2 σ ε 2 N , so the ρ ¯ k Pre-Post term improves precision in this case (a similar result holds if A k = 1 or B k = 1 ). Finally, if ρ is close to 1, the ρ ¯ k Pre - Post term can offset the sum of the ρ ¯ k Post and ρ ¯ k Pre terms in Equation 14, leading to a reduction in variance. Note that Equation 14 reduces to the AR(1) estimator in Burlig et al. (2019) assuming no staggered treatment timing, σ ε 2 = 0 , and equal time intervals. The illustrative power analysis presented in Section “An Illustrative Power Analysis” discusses other key features of Equation 14.

Next, we present a Corollary to Theorem 1 that provides new variance formulas for the point-in-time treatment effect estimators, β ^ DID , l e and β ^ DID , q .

Corollary to Theorem 1. The variance formula for β ^ DID , l e in Equation 9 that pertains to the DID estimator after l periods of treatment exposure is

Var β ^ DID , l e = 1 ( ∑ k = 1 K I l ≤ A k ) 2 ∑ k = 1 K I l ≤ A k 1 M T k + 1 M C k σ θ 2 1 + 1 B k + B k − 1 B k ρ ¯ k Pre − 2 ρ ¯ k , l Pre − Post + σ ε 2 N ,

15

where ρ ¯ k , l Pre − Post = 1 B k ∑ b = 1 B k ρ ( Time l + S k − 1 − Time b ) is calculated at postperiod l + S k − 1 . Further, for the calendar-time DID estimator, we can calculate Var β ^ DID , q using Equation 15, replacing I l ≤ A k with I q ≥ S k and calculating ρ ¯ k , q Pre − Post at postperiod q rather than at postperiod l + S k − 1 .

The result in Equation 15 follows from Theorem 1 by averaging over only those timing groups with l periods of observed exposure, setting A k = 1 , and translating exposure time into calendar time. The expression for Var β ^ DID , q can be obtained similarly.

In practice, standard errors for clustered DID designs using Equation 7 are often estimated using OLS with cluster-robust standard errors (CRSE) using data aggregated to the cluster-time level (Cameron & Miller 2015; Liang & Zeger 1986). Based on simulation evidence in Bertrand et al. (2004) using AR(1) errors, CRSE standard errors are typically adjusted for clustering at the unit of treatment assignment only.

Cross-Sectional Analysis: Power Calculations

To facilitate power calculations, it is customary to express variance formulas in terms of ICCs of the error variances. We adopt this strategy by setting σ θ 2 = ICC θ σ 2 and σ ε 2 = 1 − ICC θ σ 2 , where σ 2 = σ θ 2 + σ ε 2 is the total error variance. Using this formulation, σ 2 no longer enters the MDE formula in Equation 1, and the power calculations can be conducted by specifying a value for ICC θ rather than for the error variances directly.

Power calculations can now be conducted by inserting Equation 14 or 15 into Equation 1 to calculate MDE values for prespecified values of P, M T k , M C k , N, ICC θ , A_k , B_k , ρ , α , and λ . Calculating the degrees of freedom ( d f ) for clustered designs is complex (Donald & Lang, 2007; Hedges, 2007). However, in panel settings, it is customary to use the number of cluster-time observations adjusted by the number of model parameters (Cameron & Miller, 2015). Using Equation 7, this approach yields d f = M P − M − K P − ∑ k A k for the pooled estimator, β ^ DID,POOL . In our running example, d f = 523 (where M = 79 , P = 8 , K = 3 , and ∑ k A k = 6 ). Similarly, we can use d f = M l P − M l − K l P − K l for the point-in-time estimator, β ^ DID , l e , where M l = ∑ k I l ≤ A k M T k + M C k and K l = ∑ k I l ≤ A k are sample sizes for the included timing groups, and similarly for β ^ DID , q .

Alternatively, Equation 1 can be solved to calculate the total number of clusters, M = M T + M C , required to attain a given MDE value. For example, for the pooled DID estimator, we can calculate the required cluster sample size using

M = Factor α , λ , d f 2 V MDE 2 ,

16

where V is the variance in Equation 14, setting M T k = r p T k and M C k = 1 − r p C k , where the inputs r = M T / M , p T k = M T k / M T , and p C k = M C k / M C are user-specified. Because Factor α , λ , d f is also a function of M, iterative methods to solve nonlinear equations can be used to solve for M ( Power_Panel uses the secant method). For nonclustered designs, one can set ICC θ equal to 0, N to 1, and M to the number of individuals.

Incorporating Covariates

In DID analyses, a vector of time-varying covariates, x i j t , with associated parameter vector, γ , is often included in the models to adjust for potential confounding bias. To examine the effects of covariates on precision, to fix concepts, we begin by assuming no confounding bias by invoking a parallel trends assumption for the covariates (which we then relax):

Assumption DID.2. Parallel trends for the covariates. For each timing group k ≥ 1 , posttreatment time period q ≥ S k , and covariate c ∈ 1 , … , υ :

E x ¯ c , t = q − x ¯ ¯ c , t < S k | G T j = k = E x ¯ c , t = q − x ¯ ¯ c , t < S k | G C j = k ,

17

where υ is the number of covariates and the x ¯ c , t = q and x ¯ ¯ c , t < S k covariate means are defined analogously to the Y ¯ 0 and Y ¯ ¯ 0 outcome means defined in Equations 2 and 6.

Because γ is assumed constant over time and across timing groups, the implication of this assumption is that mean covariate values during the pre- and posttreatment periods are independent of treatment status for each timing group. The result is that the addition of covariates to the models does not change the DID estimators in expectation.

To see this result more formally, it can be shown that OLS estimators for β DID , k q with and without covariates, β ^ DID , k q , X and β ^ DID , k q , can be related using β ^ DID , k q , X = β ^ DID , k q − x ¯ t = q G T j = k − x ¯ ¯ t < S k G T j = k − x ¯ t = q G C j = k − x ¯ ¯ t < S k G C j = k ' γ ^ , where x ¯ and x ¯ ¯ are vectors of covariate means, and γ ^ is the OLS estimator for γ (see also Schochet et al., 2021). Note that the term in parentheses is a vector containing DID estimators for each covariate. But these DID estimators each have zero expectation under Assumption DID.2. Thus, E β ^ DID , k q , X = E β ^ DID , k q , even when γ ^ ≠ 0 .

Thus, under Assumption DID.2, precision gains from covariates can be quantified by multiplying the variance terms from the model without covariates by 1 − R Y X 2 , where R Y X 2 is the R ² value from covariate inclusion. We assume the same R ² value to explain σ θ 2 and σ ε 2 . With covariates, d f is reduced by the number of covariates, υ . This approach parallels the power analysis approach for RCTs where baseline covariates are independent of treatment status due to randomization (Raudenbush, 1997; Schochet, 2008).

Assumption DID.2, however, is not likely to hold in practice. Thus, precision gains from covariates will typically be reduced by treatment-covariate collinearity. Thus, the net variance reduction due to the covariates can be approximated using the ratio 1 − R Y X 2 / 1 − R T X 2 , where R T X 2 is the average R ² value from regressing the treatment interaction terms in Equation 7 on the covariates. This ratio can be less than 1. We summarize this result in the following theorem.

Theorem 2. If covariates are added to the model in Equation 7, the variance formula for β ^ DID,POOL, X can be approximated using

Var β ^ DID,POOL , X = 1 − R Y X 2 1 − R T X 2 Var β ^ DID,POOL ,

18

and similarly for the point-in-time estimators, β ^ DID , l e and β ^ DID , q , and the longitudinal estimators.

CITS and ITS Estimators

CITS estimators are based on specifications that model trends before and after the introduction of the treatment. They quantify whether once the treatment begins, the treatment group deviates from its preintervention trend by a greater amount than does the comparison group. The more common ITS estimators pertain to designs without a comparison group, where the treatment group is compared to its own preperiod trend only. CITS and ITS estimators were popularized in the seminal works of Cook and Campbell (1979) and Shadish et al. (2002); recent applications and reviews in various disciplines are provided in Bloom (2003), Somers et al. (2013), Kontopantelis et al. (2015), Linden (2015), Clair et al. (2016), Bernal et al. (2017), and Baicker and Svoronos (2019). We assume at least three (preferably four) pre- and postperiod time points each, so that the trends can be adequately modeled (Cook & Campbell, 1979).

For our CITS and ITS power analysis, we adopt an event history approach for the reasons discussed above for the DID design. We model pre- and postperiod trends using a linear specification, where we allow the intercepts and slopes to differ across the two time periods, which we refer to as the “fully interacted” model. This approach can be extended to allow for segmented regression lines across different postperiods. We also present results for a restricted model that assumes common slopes across the pre- and postperiods. Bernal et al. (2017) and Ferron and Rendina-Gobioff (2014) discuss other CITS estimands not considered here.

We rely on the following fully interacted regression model using data for all timing groups:

y i j t = ∑ k = 1 K I G j = k [ α 0 k T j Pre t + γ 0 k Time t T j Pre t + α 1 k T j Post t + γ 1 k Time t T j Post t + δ 0 k ( 1 − T j ) Pre t + τ 0 k Time t 1 − T j Pre t + δ 1 k 1 − T j Post t + τ 1 k Time t ( 1 − T j ) Post t + θ j t + ε i j t ] ,

19

where Pre t is an indicator that equals 1 if the time period corresponds to the preperiod (i.e., if t < S k ) and 0 otherwise, and Post t = 1 − Pre t is an indicator that corresponds to the postperiod, recalling that Time t denotes when the outcome was measured in calendar time relative to a common reference point.

In this model, α 0 k and γ 0 k are preperiod intercept and slope parameters for treatments in timing group k; α 1 k and γ 1 k are postperiod intercepts and slopes for treatments in timing group k; and δ 0 k , τ 0 k , δ 1 k , and τ 1 k are corresponding parameters for the matched comparison group. In this formulation, we include the treatment-by-time indicator terms (such as T j Pre t ) rather than the cluster-level fixed effects as we did for the DID design to simplify variance estimation. Using our running example, the model in Equation 19 contains 24 parameters (eight pre- and postperiod intercept and slope parameters for each of the three timing groups, split evenly between the treatment and comparison groups).

Applying OLS to the model in Equation 19, we can estimate the ATT k q estimand in Equation 2 using

β ^ CITS , k q = α ^ 1 k + γ ^ 1 k Time q − α ^ 0 k + γ ^ 0 k Time q − δ ^ 1 k + τ ^ 1 k Time q − δ ^ 0 k + τ ^ 0 k Time q = Pred 1 k q − Pred 2 k q − Pred 3 k q − Pred 4 k q .

20

The first part of this estimator, Pred 1 k q − Pred 2 k q , measures the difference between the predicted values at time q based on the fitted pre- and postperiod trend lines for the treatment group, and the second part, Pred 3 k q − Pred 4 k q , measures the pre–post difference—the “forecast error”—for the matched comparison group. Note that Pred 2 k q + Pred 3 k q − Pred 4 k q is the estimated counterfactual outcome at time q for the treatment group in the absence of the intervention; it incorporates both the forecast at time q from the fitted treatment group preperiod trend line as well as the comparison group forecast error.

Using Equation 20, we see that the CITS estimator differs from the DID estimator in that it measures treatment-comparison differences in estimated trendlines, not just estimated intercepts; thus, the CITS estimator will have a larger variance (which we formalize below). The ITS estimator, β ^ ITS , k q = Pred 1 k q − Pred 2 k q , is based on the treatment group trendlines only.

As with the DID approach, we can now average the β ^ CITS , k q estimators over k to yield the point-in-time CITS estimators, β ^ CITS , q and β ^ CITS , l e , which can then be averaged over their postperiods to yield the pooled estimator, β ^ CITS,POOL (and similarly for the ITS estimators). These estimators are unbiased for the ATT q , ATT l e , and ATT POOL estimands under Assumptions 1 and 2, assuming the linear specification in Equation 19 is correct, and under the following added assumption:

E Y ¯ t = q 0 − Pred 2 k q | G T j = k = E Y ¯ t = q 0 − Pred 4 k q | G C j = k .

21

This condition requires that in the absence of treatment, forecast errors from the preperiod trend lines would be the same, on average, for the treatment and comparison groups. The parallel condition for the ITS estimator is perfect mean forecasts for the treatment group: E Y ¯ t = q 0 − Pred 2 k q | G T j = k = 0 .

To calculate the variance of β ^ CITS , k q , we can use Equation 20 to compute the variances and covariances of the Pred k q terms. To facilitate the calculations, note that the OLS estimates from Equation 19 are identical to those from separate models estimated using four groups of observations for each timing group k: Group 1: treatments in the postperiod (to obtain Pred 1 k q ), Group 2: treatments in the preperiod (to obtain Pred 2 k q ), Group 3: comparisons in the postperiod (to obtain Pred 3 k q ), and Group 4: comparisons in the preperiod (to obtain Pred 4 k q ).

For illustration, consider OLS estimation using Group 2 data for timing group k (with M T k clusters) aggregated to the cluster-time level for a model with an intercept and time variable. For estimation, we center Time t around its preperiod mean, Time ¯ { t < S k } , which facilitates the calculations because α ^ 0 k now becomes the preperiod mean outcome, y ¯ ¯ { t < S k } G T j = k .

Consider first the model without the AR(1) error structure that we use to highlight key features of the more general, complex variance estimator presented later in Theorem 3. For this specification, we find after applying standard OLS methods that the variance of the Pred 2 k q forecast at time q based on the fitted preperiod trend line for Group 2 is