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.
128,233 characters · 37 sections · 177 citation commands
Fused Extended Two-Way Fixed Effects for Difference-in-Differences With Staggered Adoptions
\affil{University of Southern California}
Keywords: Difference-in-differences, staggered adoptions, treatment effect heterogeneity, bridge regression, fused lasso.
We investigate difference-in-differences estimation under staggered adoptions. Suppose \(N\) units are observed at \(T \geq 2\) times, where our theoretical results will assume \(N\) grows asymptotically while \(T\) is fixed. For convenience we assume a balanced panel where all units are observed at all times, so we have \(NT\) observations.
Suppose \(R \leq T - 1\) cohorts start receiving a treatment at times \(r \in \mathcal{R} =\{r_1, r_2, \ldots, r_R\} \subseteq \{2, \ldots, T\}\). For example, we might want to investigate the effect of a law or policy that is implemented state-by-state across time callaway2021difference, goodman2021difference, borusyak2024revisiting, or the effect of a product or service that is released across different regions over time de2020two. We consider settings where the treatment is an absorbing state, so that once a unit starts being treated, they continue receiving treatment thereafter, but this can be relaxed with little complication. Since no cohort enters at time \(t = 1\) and we assume the probability of never selecting into treatment is strictly positive for all units (see Assumption F2 in Section (ref)), for sufficiently large \(N\) we have untreated units available as a baseline comparison for all treated cohorts at all times.
For \(i \in \{1, \ldots, N\} =: [N]\) and \(t \in [T]\), we encode the treatment status for each unit in a random variable \(W_i \in \{0\} \cup \mathcal{R}\), where \(W_i = r\) if unit \(i\) is in cohort \(r \in \mathcal{R}\) or \(W_i = 0\) if unit \(i\) is never treated. Let \(\tilde{y}_{(it)}(r)\) be the potential outcome for unit \(i\) at time \(t\) if they were in cohort \(r\). We will be interested in estimating treatment effects of the form
as well as conditional treatment effects that depend on the value of unit covariates; see (ref), to come. It is well-known that in the \(T =2\), \(\mathcal{R} = \{2\}\) case, under no anticipation and parallel trends assumptions
and this treatment effect can be estimated consistently by the canonical two-way fixed effects regression
where \(\alpha_i^*\) are fixed effects corresponding to units, \(\gamma_t^*\) are fixed effects corresponding to time, \(\mathbbm{1}\{ \cdot \} \) is an indicator variable that equals 1 on the event \(\cdot\) and 0 otherwise, \(\epsilon_{it} \) is a noise term, and for conciseness we denote \(\tilde{y}_{it} := \tilde{y}_{it}(W_i)\), the observed response for unit \(i\) at time \(t\).
In the more general \(T \geq 2\) and arbitrary \(\mathcal{R}\) setting, one might hope that regression (ref) would aggregate the treatment effects in a sensible way to estimate some form of an average treatment effect, but it is now well-known that this isn't the case, and many alternative estimators have been proposed borusyak2018revisiting, de2020two, sun2021estimating, goodman2021difference, callaway2021difference.
Most of the alternative estimators proposed in these and other works depart from ordinary least squares (OLS) estimation, but wooldridge2021two\footnote{Along with sun2021estimating.} shows that linear regression can still be unbiased provided that enough parameters are estimated. He first shows that (ref) is also the estimand for \(\tau\) in the regression
where \( \eta \) can be interpreted as the expected outcome for the never-treated units at \(t = 1\), \(\nu_2\) is the selection bias for the treated units at \(t = 1\), and \(\gamma_2\) is the trend in the untreated potential outcomes. wooldridge2021two proposes generalizing this to the arbitrary \(T\) case and allowing for each parameter to be linear in covariates rather than constant. This leads to a linear model on cohort dummy variables, time dummies, covariates, treatment indicators, and interactions of those terms: for all \(i \in [N]\), \(t \in [T]\), and \(r \in \mathcal{R}\),
Here \(\boldsymbol{X}_{i} \in \mathbb{R}^d\) is a set of time-invariant (pre-treatment) covariates, \( \eta^* + \boldsymbol{X}_{i}^\top \boldsymbol{\kappa}^* \) is the conditionally expected outcome of the never-treated units at \(t = 1\), the \( \gamma_t^* + \boldsymbol{X}_{i}^\top \boldsymbol{\xi}_t^* \) terms are the conditional trends in the untreated potential outcomes, \(\nu_r^* + \boldsymbol{X}_{i}^\top \boldsymbol{\zeta}_r^*\) is the conditional selection bias for cohort \(r\), and \(\tau_{rt}^* + \boldsymbol{\dot{X}}_{(ir)}^\top \boldsymbol{\rho}_{rt}^* \) is the conditional treatment effect for cohort \(r\) at time \(t\). (See Assumption (LINS) in Section (ref) for more details.)
wooldridge2021two calls (ref) the extended two-way fixed effects (ETWFE) estimator wooldridge2021two and shows that in this regression \(\tau_{rt}^*\) is unbiased for \(\tau_{\text{ATT}} (r, t)\) under conditional no-anticipation and parallel trends assumptions; see Assumptions (CNAS) and (CCTS) in Section (ref). The ETWFE estimator can also be used to estimate conditional average treatment effects on the treated units
As wooldridge2021two points out, the ETWFE estimator (ref) may contain a large number of parameters to estimate if \(R\), \(T\), or \(d\) is moderately large relative to \(N\). In this case, ordinary least squares estimation of (ref) may lead to very imprecise estimates of at least some of the parameters, trading the bias problem that model (ref) has for a low-efficiency (or high-variance) problem.
wooldridge2021two suggests some possible ways to mitigate this problem by imposing restrictions in order to reduce the number of parameters to be estimated. For example, one might assume that the treatment effects only vary across cohorts and not with time, so that \(\tau_{rt} = \tau_r\) for all \(r \in \mathcal{R}\), \(t \geq r\). (Other works that consider similar restrictions include borusyak2018revisiting and borusyak2024revisiting.) If many of these restrictions hold, then imposing them results in more efficient estimation of the treatment effects than regression (ref). However, imposing these restrictions in an ad hoc way could reintroduce bias if the wrong restrictions are imposed, or leave unnecessary inefficiency if not enough restrictions are imposed.
The first goal of this paper is to propose a novel estimator that imposes these restrictions in an automatic, data-driven way. We propose the fused extended two-way fixed effects (FETWFE) estimator, a novel machine learning method for difference-in-differences estimation. Before we describe FETWFE, we discuss sparse fusion penalties in linear regression, which play a key role in our estimator.
Given a \(N \times p\) real-valued design matrix \(\boldsymbol{X}\) and a response \(\boldsymbol{y} \in \mathbb{R}^N\), the well-known lasso Tibshirani1996 estimator is defined by the optimization problem
where \(\lambda_N\) is a tuning parameter that can be chosen from a set of candidate values by cross-validation or a criterion like BIC, and recall that the \(\ell_q\) norm\footnote{It is a norm for finite \(q \geq 1\), though we also consider \(q \in (0, 1)\), in which case \(\lVert \cdot \rVert_q \) does not satisfy the triangle inequality.} of an \(n\)-vector \(\boldsymbol{v} = (v_1, v_2, \ldots, v_n)^\top \) is defined as \(\lVert \boldsymbol{v} \rVert_q := \left( \sum_{i=1}^n |v_i|^q \right)^{1/q}\) . The \(\ell_1\) penalty term leads to sparse estimated coefficients \(\boldsymbol{\hat{\beta}}(\lambda_N)\) with entries set exactly equal to 0.
The fused lasso tibshirani2005sparsity adds a second tuning parameter and \(\ell_1\) penalty term:
The fused lasso tends to yield solutions that both set some individual coefficients exactly equal to 0 (due to the first \(\ell_1\) penalty term) and also set the differences between some adjacent pairs of coefficients exactly equal to 0 (due to the second \(\ell_1\) penalty term). That is, some adjacent pairs of coefficients are set exactly equal to each other, or fused together.
FETWFE uses \(\ell_q\) fusion penalties for \(q \in (0, 2]\) (bridge regularization) to shrink certain sets of parameters in (ref) towards each other in a way that aligns with intuition and previously proposed restrictions. See Figure (ref) for a visual depiction of which of the marginal average treatment effects \(\hat{\tau}_{rt}\) we propose penalizing towards each other, and see further details about the FETWFE penalty term (ref) in Section (ref). A lasso estimator is within our framework if one chooses \(q = 1\). Choices of \(q \in (0, 1]\) will result in sparse solutions, and choices of \(q \in (0,1)\) will allow for convergence at \(1/\sqrt{N}\) rates and asymptotic normality of the estimated coefficients. These fusion penalties allow the restrictions to be selected in a hands-off, automatic way that comes with theoretical guarantees.
The second goal of this paper is to provide theoretical guarantees for the performance of FETWFE. We prove that under a suitable sparsity assumption on the true coefficients and regularity conditions, FETWFE learns the correct restrictions with probability tending to one (Theorem (ref)). FETWFE then estimates only the needed parameters, improving efficiency. We also prove that for \(q \in (0, 2]\) our proposed estimators of two classes of heterogeneous marginal average treatment effects are consistent under either conditional or marginal parallel trends, and two classes of conditional average treatment effect estimators are consistent under conditional parallel trends (Theorem (ref)). To the best of our knowledge, this is the first work establishing the consistency of a difference-in-differences estimator for marginal average treatment effects under either form of parallel trends. For estimating propensity-weighted conditional average treatment effects, the generalized propensity scores can be estimated by an arbitrary nonparametric estimator with theoretical convergence guarantees, so our theory applies to propensity score estimators as simple as multinomial logit or as flexible as deep neural networks farrell2021deep or random forests gao2022towards.
We also show that our treatment effect estimators are asymptotically normal for \(q \in (0, 1)\) and we provide feasible, provably consistent finite-sample variance estimators (Theorems (ref) and (ref)). For probability-weighted treatment effect estimates, we propose both an asymptotically normal test statistic if the generalized propensity scores are estimated on a separate data set from FETWFE (Theorem (ref)b) and an asymptotically subgaussian test statistic if both models are estimated on the same data set (Theorem (ref)c). The data-driven restrictions FETWFE learns improve parameter estimation while preserving asymptotic unbiasedness and consistency, allowing for the construction of asymptotically valid confidence intervals.
Under the assumptions of Theorem (ref) we show that FETWFE has an oracle property fan2001variable, fan2004nonconcave: the asymptotic variance of its estimates depend only on the parameters of the smaller restricted model, and it converges at the same \(1/\sqrt{N}\) rate as an OLS-estimated model on the true restricted model even if the number of covariates grows to infinity with \(N\). Our theoretical guarantees leverage extensions of the bridge regression theory of kock2013oracle (Theorems (ref) and (ref) in the appendix) that may be of independent interest.
Finally, a third goal of this paper is to show demonstrate and facilitate the practical use of FETWFE. An open-source R package implementing FEWFE is available on the Comprehensive R Archive Network Faletto:2025aa. We demonstrate that FETWFE is useful in practice through simulation studies that we conduct in Section (ref) comparing the FETWFE estimator to competitor estimators. We also illustrate FETWFE in an empirical application in Section (ref).
There has been an abundance of recent research on difference-in-differences and it is not possible to cite all of the work in this stream; see de2023tworeview, sun2022linear, roth2023s, and zimmermann2020handbookcallaway for recent reviews. lechner2011estimation reviews earlier literature. chiu2023and review the literature from a political science perspective. de2023two and arkhangelsky2023causal review causal inference in panel data more broadly. Beyond those reviews and the papers we have already cited, throughout the paper we cite papers in the difference-in-differences literature and relate them to FETWFE; see Sections (ref) and (ref) and Appendix (ref), in particular.
Besides wooldridge2021two, there are many other estimators that target similar estimands to ours, and we discuss several of these in Section (ref), though none of these are machine learning estimators. It is also natural to compare FETWFE to other machine learning and non- and semi-parametric difference-in-differences estimators. Generally these methods are more flexible in the kinds of models than they can estimate than FETWFE (which is a regularized linear model, though our theoretical results allow for nonparametric estimation of the generalized propensity scores), but much less flexible in the estimands they target. Further, these methods tend to require assumptions that an estimator is available that achieves consistency at particular rates, so applying these results would involve further checking the required assumptions for a particular estimator. In contrast, we provide “end-to-end" assumptions for the consistency, asymptotic normality, and oracle property of three classes of estimators (ref), (ref), and (ref).
The parametric and semi-parametric estimators proposed by zimmert2020efficient, sant2020doubly, and chang2020double yield theoretical results for a wide variety of machine learning estimators with theoretical guarantees under conditional parallel trends assumptions, but only in the \(T = 2\) setting and only for a single average treatment effect estimand (specifically, our estimand in Equation (ref)). blackwell2024estimating accomplish similar goals with different estimands under different assumptions, focusing mainly on a setting with \(T = 2\) and a single mediator of interest. nie2021nonparametric propose related non-parametric estimators that achieve \(1/\sqrt{N}\) convergence and asymptotic normality in the \(T = 2\) repeated cross-sections setting only for an average treatment effect, and consistency for conditional average treatment effects at much slower rates, in contrast to the \(1/\sqrt{N}\) rate for conditional average treatment effects estimated by FETWFE (see Theorem (ref)). hatamyar2023machine provide similar semi-parametric estimators for both marginal and conditional average treatment effects in the staggered adoptions setting with arbitrary \(T\), but they do not provide any theoretical guarantees.
gavrilova2023dynamic propose an approach using random forests that is somewhat more flexible in its estimands. Their approach allows for arbitrary \(T\), though they focus on the case with only one treated cohort that starts an absorbing treatment at time \(t= 2\), providing a brief sketch as to how their method might be extended to allow for staggered adoptions. Further, their method requires splitting the data and estimating completely separate models to estimate the treatment effect of the treated cohort at each time \(2, 3, \ldots, T\), while FETWFE allows for the full sample to be used to jointly estimate all treatment effects simultaneously. Indeed, a central idea of our method is that the fusion penalties we use allow FETWFE to borrow strength across cohorts and time in order to improve estimation efficiency. (None of the above works mentions fusion penalties.)
Our work also contributes to the much broader topic of estimating conditional average treatment effects. In Section (ref) we mention some prominent works in this stream of literature.
Outside of causal inference, our approach is clearly closely related to both the fused lasso tibshirani2005sparsity and the generalized lasso tibshirani2011solution; see (ref). Our work also draws much from the bridge regression theory of kock2013oracle, and, in turn, Huang2008. kaddoura2022estimation propose a linear panel data model with a fused lasso-like penalty, though not for causal inference. park2023change propose a related Bayesian approach. heiler2021shrinkage also propose a penalized panel data model that combines categories, but not in the causal inference setting. Similarly, kwon2021optimal proposes a shrinkage estimator for linear panel data models (though not for causal inference) and cites more related work in this stream.
FETWFE also relates more specifically to the generalized fused lasso hofling2010coordinate, xin2014efficient. Like FETWFE, the generalized fused lasso involves applying \(\ell_q\) penalties to the absolute differences between a variety of pairs of coefficients, though the generalized fused lasso only uses the \(q =1\) case. Some works that prove theoretical results for the generalized fused lasso include viallon2013adaptive and viallon2016robustness. viallon2013adaptive prove consistency for the generalized fused lasso under different assumptions from ours, including a restriction to the fixed \(p\) case. They prove selection consistency and asymptotic normality only for an adaptive adalasso version of their method. viallon2016robustness prove similar results (with similar limitations) for a broader class of generalized linear models.
After we introduce some notation, we establish our setup and basic assumptions in Section (ref) and present our difference-in-differences assumptions in Section (ref). In particular, we provide a testable sufficient condition (CIUN) for several kinds of estimators, including FETWFE, to be consistent under either marginal or conditional parallel trends. In Section (ref) we introduce causal estimands that FETWFE can be used to estimate, and in Section (ref) we propose the fused extended two-way fixed effects model, as well as our estimators for our proposed causal estimands. We present our theoretical guarantees for FETWFE in Section (ref). In Section (ref) we demonstrate the superior performance of FETWFE in our setting against extended two-way fixed effects (ETWFE), bridge-penalized ETWFE, and a slightly more flexible model than (ref) in simulation studies. We illustrate the usefulness of FETWFE in an empirical application using data from stevenson2006bargaining in Section (ref). Finally, Section (ref) concludes.
Throughout we will use subscripts like \(\tilde{y}_{(it)}\) to refer to the \(i^\text{th}\) unit at the \(t^{\text{th}}\) time, \(\boldsymbol{\tilde{y}}_{(i \cdot)} \in \mathbb{R}^T\) to refer to all response observations corresponding to unit \(i\), and so on. For matrices, \(\boldsymbol{Z}_{(it) \cdot}\) refers to the row corresponding to unit \(i\) at time \(t\), \(\boldsymbol{Z}_{(i\cdot) j}\) refers to the \(T\) observations of the \(j^\text{th}\) feature in unit \(i\), and so on. We will assume that unit-specific covariates are time invariant, and it will be more convenient to refer to the covariates for unit \(i\) with the notation \(\boldsymbol{X}_i\).
Recall that for a sequence of random variables \(\{X_N\}\) and a sequence of positive real numbers \(\{a_N\}\), we say that \(X_N = \mathcal{O}_{\mathbb{P}}(a_N)\) if for every \(\epsilon > 0\) there exists a finite \(M_\epsilon > 0\) and \(N_\epsilon > 0\) such that for all \(N > N_\epsilon\), \[ \mathbb{P} \left( \left| \frac{X_N}{a_N} \right| > M_\epsilon \right) < \epsilon . \] If \(a_N \to 0\), \(X_N = \mathcal{O}_{\mathbb{P}}(a_N)\) is sufficient for \(X_N \xrightarrow{p} 0\), and \(a_N\) can be interpreted as the rate of convergence dasgupta2008asymptotic.
We will expand on the setup from the introduction. For each \(i \in [N]\), assume \[ \boldsymbol{\tilde{y}}_{(i\cdot)} = \mathbb{E} [ \boldsymbol{\tilde{y}}_{(i\cdot)} \mid W_i, \boldsymbol{X}_{i} ] + c_i \boldsymbol{1}_T + \boldsymbol{u}_{(i\cdot)} , \] where \(c_i \in \mathbb{R}\) and \(\boldsymbol{u}_{(i\cdot)} \in \mathbb{R}^T\) are both random variables and \(\boldsymbol{1}_T\) is a \(T\)-vector with 1 in every entry.
Assumption (F1): \(\mathbb{E} [ c_i \mid W_i, \boldsymbol{X}_{i} ] = 0\), \( \mathbb{E} [\boldsymbol{u}_{(i\cdot)} \mid c_i, W_i, \boldsymbol{X}_{i} ] = \boldsymbol{0},\) \( \mathrm{Var} ( c_i \mid W_i, \boldsymbol{X}_{i} ) = \sigma_c^2\), and \( \mathrm{Var} ( \boldsymbol{u}_{(i\cdot)} \mid c_i, W_i, \boldsymbol{X}_{i} ) = \sigma^2 \boldsymbol{I}_T\) almost surely, where \(\sigma_c^2\) and \(\sigma^2\) are assumed to be known and finite and \(\boldsymbol{I}_T \in \mathbb{R}^{T \times T}\) is the identity matrix.
With Assumption (F1) we can use a generalized least squares transform to create a data set with conditionally uncorrelated rows. Let \(\boldsymbol{\epsilon}_{(i\cdot)} := c_i \boldsymbol{1}_T + \boldsymbol{u}_{(i\cdot)}\), and observe that for every \(i \in [N]\) almost surely \(\mathrm{Var} ( \boldsymbol{\epsilon}_{(i\cdot)} \mid W_i, \boldsymbol{X}_{i} )= \sigma^2 \boldsymbol{I}_T + \sigma_c^2 \boldsymbol{1}_T \boldsymbol{1}_T^\top =: \boldsymbol{\Omega} \in \mathbb{R}^{T \times T}\). Under (F1), for every \(i \in [N]\) almost surely \(\mathbb{E} [ \sigma \boldsymbol{\Omega}^{-1/2} \boldsymbol{\epsilon}_{(i\cdot)} \mid W_i, \boldsymbol{X}_{i} ] = \boldsymbol{0}\) and
We can carry out this transformation for any known finite and symmetric positive definite \(\boldsymbol{\Omega} = \mathrm{Var} ( \boldsymbol{\epsilon}_{(i\cdot)} \mid W_i, \boldsymbol{X}_{i} ) \). All of our results can be easily adapted to apply for any such \(\boldsymbol{\Omega}\), though we focus on the specified setting for concreteness.
Next we make an assumption on the joint distribution of the potential outcomes, treatment, and covariates. Let
be the conditional probability of a unit's treatment status, a propensity score.
Assumption (F2): The \(N\) joint random variables \((W_i, \boldsymbol{X}_{i}, c_i, \boldsymbol{u}_{(i\cdot)})_{i=1}^N\) are independent and identically distributed (iid). Further, assume that for all \(r \in \{0\} \cup \mathcal{R}\), almost surely \(0 < \pi_r(\boldsymbol{X}_i) < 1\).
Notice that for each \(\boldsymbol{X}_i\) all of the potential outcomes \(\tilde{y}_{it}(r)\) are well-defined under Assumption (F2) since treatment probabilities are always strictly between 0 and 1. So far we have made no assumptions on the relationship between the treatment assignments and the potential outcomes for a given unit, but we do so in the next section.
We will use the generalizations of the no anticipation and common trends assumptions described in the introduction---Assumptions (CNAS) and (CCTS), respectively---that were proposed by wooldridge2021two, in addition to some variations and other assumptions.
Assumption (CNAS) (“conditional no-anticipation with staggered interventions"): For treatment cohorts \(r \in \mathcal{R} \), almost surely \[ \mathbb{E} [ \tilde{y}_{(it)}(r) - \tilde{y}_{(it)}(0) \mid W_i = r, \boldsymbol{X}_{i} ] = 0 \qquad \forall i \in [N], t < r. \] This requires that before treatment, in conditional expectation each individual's expected outcome is identical to what it would have been if the unit were never treated. Some works (two recent examples are caetano2022difference and ghanem2023selection) define the potential outcomes of the units at each time in terms of their treatment status at that time rather than in terms of their cohort membership. In our notation, this amounts to assuming that \(\tilde{y}_{(i t)}(r) = \tilde{y}_{(i t)}(0)\) almost surely for \( t < r\). As wooldridge2021two points out, this stronger assumption is sufficient for (CNAS). Such a framework may not be troubling, but we only need the weaker assumption (CNAS). See further discussion in wooldridge2021two.
We will discuss and use a variety of parallel or common trends assumptions.
Assumption (CCTS) (“conditional common trends with staggered interventions"): For all \(r \in \mathcal{R}\), \( i \in [N], t \in \{2, \ldots, T\}\), almost surely
Assumption (CCTS) requires that each unit's untreated potential outcomes \(\tilde{y}_{(it)} (0)\) follow the same trend over time in conditional expectation regardless of when or if they are treated. Assumption (CCTS) is used by zimmermann2020handbookcallaway, and also by gavrilova2023dynamic in the \(\mathcal{R} = \{2\}\) case. Assumption (CCTS) has become widely adopted in the \(T = 2\) setting; early examples include heckman1997matching, heckman1998characterizing, and abadie2005semiparametric, and some more recent examples include sant2020doubly, zimmert2020efficient, chang2020double, roth2023s, and ghanem2023selection.
We also consider versions of (CCTS) that only apply before or after treatment begins for each cohort.
Assumption (CCTSB) (“conditional common trends before treatment"): (ref) holds almost surely for all \(r \in \mathcal{R} \setminus \{2\}\), \( i \in [N]\), and \(t \in \{2, \ldots, r - 1\}\).
Assumption (CCTSA) (“conditional common trends after treatment"): (ref) holds almost surely for all \(r \in \mathcal{R}\), \( i \in [N]\), and \(t \in \{r, \ldots, T\}\).
Assumption (CCTSB) requires that (CCTS) holds only until the period before each unit is treated, while (CCTSA) requires (CCTS) to hold after treatment begins. Many previous works have distinguished between the parallel trends before treatment versus after treatment, noting that (CCTSB) can be directly tested under assumption (CNAS) while (CCTSA) cannot, since untreated potential outcomes are not observed for units after treatment, and noting that parallel trends holding before treatment is no guarantee of holding after treatment. See bilinski2018seeking, kahn2020promise, dette2020difference, sun2021estimating, ban2022generalized, callaway2022pre, roth2022pretest, henderson2022complete, borusyak2024revisiting, rambachan2023more, and zimmermann2020handbookcallaway. See wooldridge2021two for an example of how (CCTSB) can be tested under (CNAS).
We can contrast (CCTS) with an unconditional or marginal version, which wooldridge2021two calls Assumption (CTS).
Assumption (CTS) (“common trends with staggered interventions"): Almost surely for all \(r \in \mathcal{R}\)
We will also refer to a version of (CTS) that holds only after treatment:
Assumption (CTSA) (“common trends with staggered interventions after treatment"): (ref) holds almost surely for all \( i \in [N]\), \(r \in \mathcal{R}\), and \( t \in \{r, \ldots, T\}\).
Because (CTS) does not account for covariates, (CCTS) is often thought of as more plausible than (CTS). We generally agree, though in Appendix (ref) we add some nuance by pointing out that (CTS) can hold when (CCTS) does not hold (Theorem (ref)). In such a setting it is not possible to estimate conditional average treatment effects using regression (ref); we show this formally in Theorem (ref)(a) in the appendix. Further, it turns out that even marginal average treatment effect estimates from models like FETWFE that rely on (CCTS) will in general be inconsistent.
However, (CTS) holding is enough to estimate marginal average treatment effects consistently using many estimators, so it is reasonable to wonder whether regression (ref) might be able to estimate treatment effects consistently under (CTS) and some additional assumptions. The following assumption will turn out to allow FETWFE to consistently estimate marginal average treatment effects if (CTS) holds, even if (CCTS) does not.
Assumption (CIUN) (“conditional independence of trends of the untreated potential outcomes for never-treated units"): Almost surely for all \(i \in [N]\) and any \(t \in \{2, \ldots, T\}\), \[ \mathbb{E} \left[ \tilde{y}_{(i t)}(0) - \tilde{y}_{(i1)}(0) \mid W_i = 0 \right] = \mathbb{E} \left[ \tilde{y}_{(i t)}(0) - \tilde{y}_{(i1)}(0) \mid W_i = 0, \boldsymbol{X}_{i} \right] . \]
Assumption (CIUN) requires that the trends in the untreated potential outcomes are mean-independent of the covariates conditional on treatment status for units that are untreated. Since (CIUN) depends only on the observed potential outcomes, it can be tested, as we discuss in the next section.
We present an assumption that is slightly different than one introduced by wooldridge2021two on the functional forms of the conditional distributions of the untreated potential outcomes, the trends \(\tilde{y}_{(i t)}(0) - \tilde{y}_{(i1)}(0) \), and the difference-in-differences estimands.
Assumption (LINS): For each \(N \in \mathbb{N}\), assume we observe \(d_N\) covariates, where \(d_N\) may be fixed or increasing in \(N\). There exist fixed parameters \(\eta^* \in \mathbb{R}, \{\nu_r^*\}_{r \in \mathcal{R}} \in \mathbb{R}\), \(\{\gamma_t^*\}_{t \in \{2, \ldots, T\}} \in \mathbb{R}\), and \(\{ \tau_{rt}^* \}_{r \in \mathcal{R}, t \geq r}\in \mathbb{R} \), along with fixed sequences of parameters \(\{\boldsymbol{\kappa}_N^*\}_{N \in \mathbb{N}} \in \mathbb{R}^{d_N}, \{ \boldsymbol{\zeta}_{rN}^*\}_{r \in \mathcal{R}, N \in \mathbb{N}} \in \mathbb{R}^{d_N}, \{\boldsymbol{\xi}_{tN}^*\}_{t \in \{2, \ldots, T\}, N \in \mathbb{N}} \in \mathbb{R}^{d_N}\), and \( \{ \boldsymbol{\rho}_{rtN}^* \}_{r \in \mathcal{R}, t \geq r, N \in \mathbb{N}} \in \mathbb{R}^{d_N}\) such that for each \(N \in \mathbb{N}\), for all \(i \in [N]\) and any fixed \(\boldsymbol{x} \in \mathbb{R}^{d_N}\) in the support of the random vector \(\boldsymbol{X}_i\) it holds that
where \[ \boldsymbol{\dot{x}}_r := \boldsymbol{x} - \mathbb{E} \left[ \boldsymbol{X}_{i} \mid W_i = r \right] \] are covariate vectors that are centered with respect to their cohort means. (For conciseness, hereafter we omit \(N\) from the subscripts of the parameters defined here that may vary with \(N\).)
The version of this assumption presented by wooldridge2021two uses a slightly different expression for (ref) and requires \(d_N = d\) to be fixed, since ETWFE is inconsistent in the high-dimensional setting we consider. Notice the similarity of the quantity on the left side of (ref) to the canonical two-period difference-in-differences estimand (ref). wooldridge2021two shows that Assumption (LINS), along with Assumptions (CNAS) and (CCTS), make (ref) correctly specified for the potential outcomes and the the left side of (ref) equal to the conditional average treatment effects. See also our derivation (ref) and Theorem (ref)(a) in the appendix.
Under Assumption (CIUN) and (LINS), the covariates \(\boldsymbol{\xi}_t^*\) all equal \(\boldsymbol{0}\). So a test of this null hypothesis (for example, a Wald test) serves as a test for the null hypothesis that both (CIUN) and (LINS) hold.
Next we will discuss several classes of causal estimands we hope to estimate using FETWFE. The building block of our estimands will be the conditional average treatment effects on the treated units \(\tau_{\text{CATT}} (r, t, \boldsymbol{x}) \) from (ref), which matches definition (6.23) in wooldridge2021two. Estimating \( \tau_{\text{CATT}} (r, t, \boldsymbol{x}) \) is useful if we are interested in identifying which units are most affected by the treatment at a specific time. Besides economics florens2008identification, djebbari2008heterogeneous, domains where conditional average treatment effect estimation have been found useful include personalized medicine feinstein1972estimating, lesko2007personalized,frankovich2011evidence, foster2011subgroup, education morgan2001counterfactuals, brand2010benefits, murphy2016handbook, targeted marketing bottou2013counterfactual, athey-yelp-2018, hitsch2018heterogeneous, ascarza2018retention, yang2020targeting, ellickson2022estimating, psychology bolger2019causal, winkelbeiner2019evaluation, sociology xie2012estimating, breen2015heterogeneous, and political science imai2011estimation, grimmer2017estimating.
Many methods for estimating conditional average (or heterogeneous\footnote{Within the differences-in-differences literature, the descriptor “heterogeneous" has often been used in reference to treatment effects that vary across cohorts or across time; e.g., de2020two, sun2021estimating, wooldridge2021two. Here we refer to heterogeneity in treatment effects across the distribution of \(\boldsymbol{X}_i\), though our model also accommodates heterogeneity in treatment effects across cohorts and time. To avoid ambiguity, we will generally use the term “conditional average treatment effects."}) treatment effects have been developed within econometrics over the last decade; see abrevaya2015estimating,athey_imbens_2016, Wager2018, nie2021quasi, semenova2021debiased, and fan2022estimation for a few examples.
We can use (ref) to construct a variety of causal estimands. For any set of finite constants \(\{\psi_{rt}\}_{r \in \mathcal{R}, t \in \{r, \ldots, T\}}\), consider the estimand
a weighted average of the cohort- and time-specific conditional average treatment effects. This is similar to the framework of estimation targets of borusyak2024revisiting and callaway2021difference, though they study average treatment effects (like our estimand (ref), to come) rather than conditional average treatment effects. Refer to these papers for many examples of choices of \(\{\psi_{rt}\}\). Similar aggregations are also discussed by RePEc:qed:wpaper:1495 and wooldridge2021two.
To give one example among many possible estimands of interest using specific \(\{\psi_{rt}\}\), we could choose for some \(r' \in \mathcal{R}\)
to define the conditional average treatment effect across time for a unit in cohort \(r'\) with covariates \(\boldsymbol{x}\):
A version of (ref) for average (rather than conditional average) treatment effects was considered by wooldridge2021two, and a slightly more general form (also accommodated by our framework) was considered by goodman2021difference.
A second class of causal estimands can be generated by weighting each \(\tau_{\text{CATT}} (r, t, \boldsymbol{x}) \) by the probability a unit is in cohort \(r\) given that they are treated conditional on \(\boldsymbol{x}\),
In particular, for any set of finite constants \(\{\psi_{rt}\}\), consider the estimands
One example of an estimand we can form in this class arises from choosing
which yields a conditional average treatment effect by marginalizing (ref) across treatment status:
This estimand is closely related to many previously considered heterogenous treatment effects; one of many works that considers a similar estimand (outside of the specific context of difference-in-differences) is Wager2018.
Outside of conditional average treatment effects, we can also marginalize these conditional treatment effects over the distribution of \(\boldsymbol{X}_i\) to get average treatment effects. Since these estimands are "coarser," it may be possible to get more precise estimates of these estimands given a fixed sample size.
For \(\tau_{\text{ATT}} (r, t) = \mathbb{E} \left[ \tau_{\text{CATT}} ( r, t, \boldsymbol{X}_i) \mid W_i = r \right] \) from (ref), consider the classes of estimands
and
where \(\tilde{\pi}_r := \mathbb{P}(W_i = r \mid W_i \neq 0) = \mathbb{E}[ \tilde{\pi}_r (\boldsymbol{X}_{i}) \mid W_i \neq 0 ] \) is the marginal probability of a unit being in cohort \(r\) given that it is treated.
Estimand (ref) was used by callaway2021difference and hatamyar2023machine, and it matches the estimand from wooldridge2021two and RePEc:qed:wpaper:1495. It is also closely related to the cohort-specific average treatment effect on the treated (3) from sun2021estimating (though we use the acronym CATT to refer to conditional, not cohort, effects) and estimands from goodman2021difference.
As we mentioned earlier, the class of estimands (ref) is the same framework of estimation targets considered by borusyak2024revisiting and is similar to aggregations considered by callaway2021difference, wooldridge2021two, RePEc:qed:wpaper:1495, and hatamyar2023machine.
The probability-weighted class of estimands (ref) bears a resemblance to the probability-weighted estimands proposed by callaway2021difference. Theorem (ref) in the appendix extends our consistency result to accommodate classes of estimands with much broader forms of estimated probability weights than (ref) and (ref), like estimand (26) of sun2021estimating or the broad class of estimands used by callaway2021difference. Theorems (ref) and (ref) extend our asymptotic normality results to these more general estimators as well.
As a particular choice of \(\{\psi_{rt}\}\), we will again highlight (ref) for any \(r' \in \mathcal{R}\), which for estimand (ref) marginalizes \(\tau_{\text{CATT}} ( r', \boldsymbol{X}_{i}) \) over the distribution of \(\boldsymbol{X}_i\):
Estimand (ref) is identical to estimand (3.7) in callaway2021difference, and is also the estimand of estimator (6.41) in wooldridge2021two.
Lastly, one of the most basic causal parameters of interest is the average effect of treatment on the treated group (ATT). By using the \(\psi_{rt}\) from (ref) with (ref), we can express this estimand as
Estimand (ref) is similar to estimand (6.40) proposed by wooldridge2021two (a similar estimand is also mentioned by borusyak2024revisiting), though we weight the cohort average treatment effects \(\tau_{\text{CATT}} (r, \boldsymbol{X}_{i}) \) by the probabilities \(\tilde{\pi}_r(\boldsymbol{X}_{i}) \) rather than taking an average where each average treatment effect in each time period is weighted equally (we could also take an equal average with our estimand from Equation (ref)). This probability-weighted average is closely related to estimand (26) from sun2021estimating.
The virtues of estimand \(\tau_{\text{ATT}}\) include summarizing the effect of the treatment in a single number and potentially being a less noisy quantity to estimate. Estimating \(\tau_{\text{ATT}}\) will involve taking a weighted average over many estimated coefficients rather than relying on a single, less precise, estimated coefficient. In contrast, more specific estimands like \( \tau_{\text{CATT}} (r, t, \boldsymbol{x}) \) will typically be difficult to estimate with precision because their estimators will depend heavily on a small subset of the observed data (though this problem will be mitigated by our fusion penalty, discussed in the next section, which allows us to borrow information from nearby estimators).
Notice that \( \tau_{\text{ATT}} \) can be obtained not just by marginalizing \( \tau_{\text{CATT}} ( \boldsymbol{X}_{i} )\) over the distribution of \(\boldsymbol{X}_i\), but also by marginalizing \(\tau_{\text{ATT}} ( W_i)\) over the distribution of cohort assignments:
which shows that estimand (ref) is in the class defined by (ref) and matches estimand (3.11) from callaway2021difference. This formulation will also be useful for estimation.
We will carry out regression on the same design matrix \(\boldsymbol{\tilde{Z}}\) proposed by wooldridge2021two to estimate regression (ref) by OLS, but we will add fusion penalties to leverage our belief that some of the parameters are equal to each other. This design matrix contains \(R\) columns of cohort dummies, \(T - 1\) columns of time dummies for times \(\{2, \ldots, T\}\), \(d_N\) covariates, \(\mathfrak{W} := \sum_{r \in \mathcal{R}} (T - r + 1)\) treatment dummies for each possible base treatment effect \(\tau_{rt}\), and \(d_N(R + T - 1 + \mathfrak{W})\) interactions between covariates and the cohort dummies, time dummies, and treatment dummies. As is typical in penalized regression, we will center \(\boldsymbol{\tilde{y}}\) and the columns of \(\boldsymbol{\tilde{Z}}\) before estimating our regression. It is not necessary to estimate the intercept term \(\eta^*\) in order to estimate the treatment effects.
Let \(p_N\) be the number of columns of \(\boldsymbol{\tilde{Z}}\); that is, \(\boldsymbol{\tilde{Z}} \in \mathbb{R}^{NT \times p_N}\). We will assume \(p_N \leq N T\) (see Assumption R2 later), but as we discussed in the introduction, it is clear that \(p_N = R + T - 1 + \mathfrak{W} + d_N(1 + R + T - 1 + \mathfrak{W})\) may be fairly large compared to \(NT\), particularly if \(d_N\) is not very small. This could lead to imprecise treatment effect estimates if (ref) is estimated by OLS. wooldridge2021two suggests alleviating this problem by assuming restrictions in the model---assuming some of these parameters equal each other. Some of the possibilities wooldridge2021two suggests include assuming that the treatment-covariate interactions \(\boldsymbol{\rho}_{rt}\) are fixed across time or assuming the within-cohort treatment effects \(\tau_{rt}\) are fixed across time. goodman2021difference also considers the possibility of treatment effects that are fixed across time. Another weaker assumption wooldridge2021two proposes is assuming that different treatment effects do not need to be estimated for every individual time period since treatment starts; instead, these times could be condensed into “early treated" and “late treated" effects, for example.
These suggestions seem sensible, but making these simplifying assumptions without strong justification risks re-introducing the bias that was removed by adding these parameters to begin with. In settings where a practitioner believes some of these restrictions hold but has no clear way to choose which particular parameters may be equal, we instead propose a data-driven approach to estimate which, if any, of these parameters may be equal.
As described in the introduction, we will assign \(\ell_q\) fusion tibshirani2005sparsity penalties to the differences between estimated parameters which may be equal, and then carry out a penalized linear regression. For example, for the treatment effects, we propose penalizing the absolute differences between treatment effects within the same cohort at adjacent times, \(|\tau_{r_k,t} - \tau_{r_k,t-1}|\), and we propose the same penalty structure for interaction effects, \(|\rho_{r_k,t,j} - \rho_{r_k,t-1,j}|\) separately for each \(j \in [d_N]\). This leverages our belief that some of these effects may be fixed across time while allowing the optimization problem to learn from the data which restrictions to choose.
We also propose penalizing the absolute differences between initial treatment effects for each cohort, \(|\tau_{r_k,r_k} - \tau_{r_{k-1},r_{k-1}}|\), and likewise we penalize the absolute differences between adjacent initial interaction effects \(|\rho_{r_k,r_k,j} - \rho_{r_{k-1},r_{k-1},j}|\) for each \(j\). This is motivated by the idea that perhaps time since beginning treatment could be the only dimension in which treatment effect changes. It is also natural to penalize the adjacent cohort fixed effects absolute differences \(|\nu_{r_k} - \nu_{r_{k-1}}|\) and time fixed effects absolute differences \(|\gamma_t - \gamma_{t-1}|\) towards each other, with analogous penalties for the interaction terms \(|\zeta_{r_kj} - \zeta_{r_{k-1},j}|\) and \(|\xi_{tj} - \xi_{t-1,j}|\).
Finally, we may also want to regularize the coefficients directly, as is common in regressions with large numbers of parameters to estimate relative to the number of observations. Because of the fusion penalties we employ, regularizing only a handful of base terms allows for this regularization. We can also penalize all of the base covariate coefficients directly, allowing us to screen out irrelevant covariates as \(d_N\) grows.
In summary, we propose estimating a model like (ref) but with a penalty term
where \(\boldsymbol{\beta} \in \mathbb{R}^{p_N}\) collects all of the coefficients. (Refer again to Figure (ref) for a graphical depiction of the penalty terms on the estimated treatment effects \(\tau_{rt}\).) We can more compactly express this as
where \(\boldsymbol{D}_N \in \mathbb{R}^{p_N \times p_N}\) is a suitably defined differences matrix which we later define explicitly in (ref).
It is much easier to see that these candidate restrictions are plausible due to the time-structured nature of the parameters than it is to know exactly which restrictions to choose. FETWFE allows practitioners to leverage this structure to improve estimation efficiency without requiring them to choose restrictions by hand, which could introduce bias if too many restrictions are chosen or compromise efficiency if too few restrictions are chosen.
Finally, we take care of some details. First we account for the random effects discussed in Section (ref). As suggested by (ref), in practice we will arrange the rows of \(\boldsymbol{\tilde{Z}}\) to contain consecutive blocks of \(T\) rows for each observation and carry out a regression on \(\boldsymbol{\tilde{y}}\) and \(\boldsymbol{\tilde{Z}}\) transformed by
where \(\otimes\) denotes the Kronecker product (recall that \(\boldsymbol{\Omega}\) was defined in Section (ref)). We will also center both the response \(\boldsymbol{G}_N \boldsymbol{\tilde{y}}\) and each column of the covariates \(\boldsymbol{G}_N \boldsymbol{\tilde{Z}}\), as is common in penalized regression. Denote these centered objects by \(\boldsymbol{y}\) and \(\boldsymbol{Z}\). Then our optimization problem is
where \(\lambda_N\) is a tuning parameter. For brevity, we will hereafter write \(\boldsymbol{\hat{\beta}}_N^{(q)}(\lambda_N) = \boldsymbol{\hat{\beta}}^{(q)}\).
As shown in Theorem (ref)(a), under (CNAS), (LINS), and (CCTS) \[ \tau_{\text{CATT}} (r, t, \boldsymbol{x}) = \tau_{rt}^* + ( \boldsymbol{x} - \mathbb{E} \left[ \boldsymbol{X}_{i} \mid W_i = r \right] )^\top \boldsymbol{\rho}_{rt}^* \] for any fixed \(\boldsymbol{x}\). Since under these assumptions (and other regularity conditions) FETWFE consistently estimates \( \tau_{rt}^* \) and \( \boldsymbol{\rho}_{rt}^* \), we can estimate this quantity by
where \(\hat{\tau}_{rt}^{(q)}\) and \(\hat{\rho}_{rt}^{(q)}\) are the estimators of \(\tau_{rt}\) and \(\rho_{rt}\) from \(\boldsymbol{\hat{\beta}}^{(q)}\) as defined in (ref) and \[ \boldsymbol{\overline{X}}_r = \frac{1}{N_r} \sum_{\{i: W_i = r\}} \boldsymbol{X}_{i} \] is the sample mean of the covariates for units observed in cohort \(r\), with \(N_r := \sum_{i=1}^N \mathbbm{1}\{W_i = r\}\) the total number of units in cohort \(r\). The forms of (ref) and (ref) then naturally suggest the estimators
and
where \(\hat{\pi}_r(\boldsymbol{x})\) is an estimator of the propensity scores \(\pi_r(\boldsymbol{x})\). We will discuss what properties this estimator \(\hat{\pi}_r(\boldsymbol{x})\) should have in more detail in Section (ref).
Finally we discuss estimating the average treatment effects marginalized over \(\boldsymbol{X}_i\). As we showed in Theorem (ref)(a) and (b), because the treatment effects are interacted with the covariates centered with respect to their cohort means, under either (CCTS) or (CTS) and (CIUN) we can estimate the average treatment effects using
as proposed by wooldridge2021two. In particular, we propose the estimators
and
for estimands (ref) and (ref), respectively. Notice that estimator (ref) avoids requiring the specification of an estimator \(\boldsymbol{\hat{\pi}}(\boldsymbol{x})\) of the generalized propensity scores by instead using the observed cohort counts to estimate the marginal cohort assignment probabilities, as we alluded to at the end of Section (ref).
Before we present theoretical guarantees for FETWFE, we will state a sparsity assumption that we will require along with some additional regularity conditions. We will require fewer assumptions to prove consistency alone than we will to prove restriction selection consistency and asymptotic normality. We emphasize that the conceptually most important assumptions are (CCTSB), (CCTSA), and (CTSA), which validate the difference-in-differences approach; (LINS), which justifies a model that is linear in the covariates \(\boldsymbol{X}_i\); and Assumption S(\(s_N\)), defined below, which assumes sparsity in a sense that matches our penalty (ref). (See Remark (ref) for some comments on how this can be made more flexible.)
Assumption S(\(s_N\)): Assume the vector \(\boldsymbol{D}_N \boldsymbol{\beta}_N^*\) (where the differences matrix \(\boldsymbol{D}_N\) associated with our penalty (ref) is defined in Equation (ref)) is sparse, with indices \(\mathcal{S} \subseteq [p_N]\) nonzero and the remaining indices in \( [p_N] \setminus \mathcal{S}\) equal to 0, where \( |\mathcal{S}| = s_N\).
Assumption S(\(s_N\)) allows us to have sparsity in the appropriate sense: the differences between parameters that we penalized in (ref) are sparse, and irrelevant covariates have true coefficient equal to 0. That is, all but \(s_N\) of the restrictions that we expect to hold will in fact hold. Crucially, Assumption S(\(s_N\)) does not require us to know which restrictions to select, or even how many restrictions to select---the value of \(s_N\) does not need to be known. Our consistency and selection consistency theorems will allow \(s_N\) to grow with \(N\) subject to certain constraints in later assumptions.
The remaining regularity conditions for Theorem (ref) are more routine. Before we describe them in detail, we briefly characterize them informally.
In particular, these assumptions allow the number of covariates \(d_N\) and the number of relevant coefficients \(s_N\) to grow to infinity with \(N\) and allow for unbounded distributions of the covariates.
We now state and describe the assumptions.
Assumption (R1): For each \(i \in [N]\), \(\mathbb{E} [ X_{ij}^4 ] < \infty\) for all \(j \in [d_N]\) and \(\mathbb{E} [ u_{(it)}^4 ] < \infty\) for all \(t \in [T] \). Further, for some \(K < \infty\), \[ \sup_{1 \leq N \leq \infty} \left\{ \frac{1}{p_N} \sum_{t=1}^T \sum_{j=1}^{p_N} \mathbb{E} \left[ ( \boldsymbol{Z} \boldsymbol{D}_N^{-1})_{(1t)j}^2 \right] \right\} \leq K . \] We show later (Lemma (ref)) that \(\boldsymbol{D}_N\) is invertible. The matrix \(\boldsymbol{Z} \boldsymbol{D}_N^{-1}\) is of interest because it turns out this is the design matrix we will actually use for estimation; see the proof of Theorem (ref) for details. Because the singular values of \(\boldsymbol{D}_N^{-1}\) are bounded away from 0 and from above for all \(N\) even if \(d_N \to \infty\) (see Lemma (ref)), this is qualitatively similar to making the same statement about \(\boldsymbol{Z}\) itself. We discuss a closely related point in more detail below.
Before we state Assumption (R2), we will define some needed notation. For any matrix \(\boldsymbol{A} \in \mathbb{R}^{N \times p_N}\), denote the empirical Gram matrix by \(\boldsymbol{\hat{\Sigma}}(\boldsymbol{A})\):
Notice that if \(\boldsymbol{A}\) has centered columns then \(\boldsymbol{\hat{\Sigma}}(\boldsymbol{A})\) is the estimated covariance matrix. Let \(e_{1N}\) and \(e_{2N}\) be the smallest and largest eigenvalues of \(\boldsymbol{\hat{\Sigma}}(\boldsymbol{Z})\).(Since \(\boldsymbol{\hat{\Sigma}}(\boldsymbol{Z})\) is random, \(e_{1N}\) and \(e_{2N}\) are as well.)
Assumption (R2): Assume \(e_{1N}\) is almost surely upper-bounded by a constant for all \(N\) and \(h_N' \xrightarrow{a.s.} 0\), where
and \(\xrightarrow{a.s.}\) denotes almost sure convergence. Also, assume \( \lim_{N \to \infty} \lambda_N \sqrt{ \frac{s_N }{N}} = 0 \).
kock2013oracle points out that (ref) requires the design matrix to be full rank, so this assumption requires that \(p_N\) grows slower than \(N \), \(\lambda_N\) does not grow too quickly, and \(e_{1N}\) does not vanish too quickly. (In practice, in finite samples \(\lambda_N\) can be selected by cross-validation or a criterion like BIC, similarly to the lasso.) Under the qualitative assumption that the eigenvalues of \(\boldsymbol{\hat{\Sigma}}(\boldsymbol{Z})\) and \(\boldsymbol{\hat{\Sigma}}(\boldsymbol{G}_N\boldsymbol{\tilde{Z}})\) are “close,"\footnote{This can be made formal by noting that centering the columns of \(\boldsymbol{G}_N\boldsymbol{\tilde{Z}}\) can be characterized as multiplication on the left by a centering matrix \(\boldsymbol{C}\), then bounding the singular values of \(\boldsymbol{C}\). See, for example, Section 8.1 of ding2024linear for further details.} (R2) roughly requires that \(\boldsymbol{G}_N\boldsymbol{\tilde{Z}}\) is full rank, and we can provide sufficient conditions for this to hold. First, the singular values of \(\boldsymbol{G}_N\) are bounded from above and below by positive constants for all \(N\) (see Lemma (ref) in the appendix), so \(\boldsymbol{\tilde{Z}}\) being full rank is enough for \(\boldsymbol{G}_N\boldsymbol{\tilde{Z}}\) to be full rank. For \(\boldsymbol{\tilde{Z}}\) to be full rank, we require at least \(d_N + 1\) observed units in each cohort (see Lemma (ref) in the appendix), all of the marginal cohort probabilities must be large enough relative to the growth rate of \(d_N\) as \(N \to \infty\), and the distribution of \(\boldsymbol{X}_i\) must be well-behaved as \(N \to \infty\).
Assumption (R3): There exist constants \(0 < b_0 < b_1 < \infty\) such that \(b_0 \leq \min_{j \in \mathcal{S}} \left\{ \left| ( \boldsymbol{D}_N \boldsymbol{\beta}_N^*)_j \right| \right\} \leq \max_{j \in \mathcal{S}} \left\{ \left| ( \boldsymbol{D}_N \boldsymbol{\beta}_N^*)_j \right| \right\} \leq b_1\).
The lower bound in Assumption (R3) is essentially a signal strength assumption.
Note that we allow \(s_N \to \infty\) and \(p_N \to \infty\) as long as the above assumptions are satisfied.
We are now prepared to state our consistency theorem.
Theorem (ref) shows that all of our causal estimators are consistent and characterizes their rates of consistency. In particular, it also shows that the classes of marginal parallel trends estimators (ref) and (ref) are consistent under either (CCTSA) or (CIUN) and (CTSA) for the marginal average treatment effects (ref) and (ref). Although we require (CCTSB) to hold regardless of whether (CIUN) holds, (CCTSB) also depends only on the observed outcomes and is testable. In the next section, we show that under mild assumptions FETWFE automatically provides a valid estimate of whether (CIUN) holds (by estimating whether the \(\boldsymbol{\xi}_t^*\) coefficients all equal \(\boldsymbol{0}\)).
To interpret the rate of convergence \( \min\{h_N, h_N'\}\), consider the upper bound \(h_N\) from (ref). Suppose the minimum eigenvalue \(e_{1N}\) is bounded away from 0 with high probability, which happens under relatively mild assumptions in our \(p_N \leq NT\) setting (see our more detailed discussion later in Section (ref) about Assumption R6). Then \( \mathcal{O}_\mathbb{P} \left( \min\{h_N, h_N'\}\right) = \mathcal{O}_\mathbb{P} \left( \sqrt{p_N/N} \right)\); see Lemma (ref) in the appendix.
In general we will focus more on \(q \in (0, 1)\) as required by Theorems (ref), (ref), and (ref), but we note that Theorem (ref) holds for any \(q > 0\). This is useful because the convex lasso (\(q = 1\)) and ridge (\(q = 2\)) optimization problems, for example, may be more tractable on very large data sets than the nonconvex bridge estimator with \(q < 1\).
To give a concrete example of an estimator \(\boldsymbol{\hat{\pi}}(\boldsymbol{x})\), in the setting of Theorem (ref)(b) where \(d_N = d\) is fixed, the multinomial logit model is \(1/\sqrt{N}\)-consistent for the generalized propensity scores by standard maximum likelihood theory under correct specification. If this model is overly simplistic or \(d \) is large relative to \(N\), regularized machine learning methods with theoretical convergence guarantees---like deep neural networks farrell2021deep, random forests gao2022towards, the \(\ell_1\)-penalized classifiers of levy2023generalization, or the group lasso multinomial logit model of farrell2015robust---would also work. On the other hand, since the cohorts have a natural ordering, we could add more structure by using an ordinal response estimator like the proportional odds model mccullagh1980regression.
It may be possible to relax Assumption (LINS) so that linearity only has to hold approximately, along the lines of e.g. Condition ASTE in belloni2014inference. (At a high level, we might expect that if approximate linearity holds so that the error due to misspecification is on the same order as the estimation error of the population least squares model, we can achieve convergence at the same rate as when the linearity assumption holds exactly.)
It may also be possible to relax the fixed \(d\) assumption in part \((b)\) by using high-dimensional central limit theorems like those in lopes2022central and chernozhukov2023nearly, to name two recent examples among many papers containing such results.
Lastly, we point out that under Theorem (ref) the generalized propensity scores and treatment effects can be estimated on the same data sets; that is, the full data set can be used for both. Intuitively, this is because these estimators are used for the distinct tasks of estimating the effects \(\tau_{\text{CATT}} (r, t, \boldsymbol{x})\) and estimating the weights used to combine them in a weighted average. For technical details, see the proof of Theorem (ref) in Appendix (ref).
In the appendix, we present an extension of Theorem (ref), Theorem (ref), which proves consistency for a set of estimators that allows for much broader functions of the estimated cohort membership probabilities than the classes (ref) and (ref). This result includes, for example, all of the estimands from callaway2021difference as special cases.
Again, we first state some additional needed regularity conditions and then present our selection consistency result, Theorem (ref). Theorem (ref) shows that with probability tending towards one FETWFE identifies the correct restrictions, resulting in fusing together parameters that equal each other, screening out irrelevant covariates, and improving estimation efficiency. At a high level, Assumptions (R4) and (R5) add the following requirements:
Assumption (R4): There exists a positive constant \(e_{\text{max}} < \infty\) such that \( e_{2N} \leq e_{\text{max}}\) almost surely (recall from Section (ref) that \(e_{2N}\) is the largest eigenvalue of \(\boldsymbol{\hat{\Sigma}}(\boldsymbol{Z})\)).
A sufficient condition for Assumption (R4) is that the distribution of \(\boldsymbol{X}_i\) has bounded support in the sense that \(\lVert \boldsymbol{X}_i \rVert_2 \leq B \) almost surely for some \(B < \infty\).
Assumption (R5): \[ \lambda_N \frac{e_{1N}^{2- q}}{ \sqrt{N^q p_N^{2 - q }}} \xrightarrow{a.s.} \infty . \] We will briefly examine how Assumption (R5) interacts with Assumption (R2). For simplicity, assume the minimum eigenvalue \ \(e_{1N} \geq c\) for some \(c > 0\) with probability tending towards 1, so that we can ignore \(e_{1N}\) asymptotically (using a similar argument to the one in the proof of Lemma (ref) in the appendix). Then Assumption (R5) is equivalent to the deterministic condition \[ \lambda_N \frac{1}{ N^{q/2} p_N^{1 - q/2 }} \to \infty . \] If we also assume \(s_N = s\) is fixed, we have from Assumption (R2) that
We will require \(q \in (0,1)\), and we see that under these assumptions for (R5) and (R2) to hold simultaneously we require \(N^{q/2}p_N^{1-q/2}\) to grow slower than \(N^{1/2}\); that is, \(p_N\) must grow slower than \[ \left( \frac{N^{1/2}}{N^{q/2}} \right)^{2/(2-q)} = N^{(1-q)/(2-q)} . \]
Even if \(e_{1N}\) is not bounded away from 0 with high probability, Assumption (R5) requires that \(e_{1N}\) does not vanish too quickly. As we discussed in reference to Assumption (R2), this requires that the probability of being in any one cohort does not vanish too quickly as \(N \to \infty\). It also requires that \(\lambda_N\) not decrease too quickly.
Notice that other than these relatively mild regularity conditions on the eigenvalues of \(\boldsymbol{\hat{\Sigma}}(\boldsymbol{Z})\), we do not require a condition like the irrepresentable condition Zhao2006 or neighborhood stability condition meinshausen2006high that are required for selection consistency of the lasso. However, our restriction selection consistency result will only work for \(q < 1\); we do not prove selection consistency of FETWFE with a lasso (\(q =1\)) penalty.
Theorem (ref) shows that FETWFE fulfills its promise to identify the correct restrictions in accordance with our \(\ell_q\) penalty structure, as discussed in Section (ref). This serves our goal of improving estimation by avoiding unnecessarily estimating parameters that are in fact equal.
FETWFE also screens out irrelevant covariates and interactions. As mentioned in Section (ref), if Assumptions (LINS) and (CIUN) both hold than the \(\boldsymbol{\xi}_t^*\) coefficients must all equal 0. If this is true, in large enough samples FETWFE will estimate \(\boldsymbol{\xi}_t^*\) to equal exactly \(\boldsymbol{0}\) with high probability.
The proof of Theorem (ref) relies on an extension of Theorem 2 of kock2013oracle (Theorem (ref) in the appendix) to show that not only does the bridge estimator exclude false selections, it also selects all of the relevant components of \( \boldsymbol{D}_N \boldsymbol{\beta}_N^*\) with probability tending to 1. Theorem (ref) requires no added assumptions from Theorem 2 of kock2013oracle.
One consequence of Theorem (ref) is that if \(\tau_{\text{ATT}} (r, t) = 0\) for any \(r\), then \[ \lim_{N \to \infty} \mathbb{P} \left( \hat{\tau}_{\text{ATT}} ( r, t) = 0\right) = 1. \] The same is true of all of the other causal estimators. This is stronger than convergence in probability to 0, and this stronger notion of convergence turns out to be a key ingredient in the proof of Theorems (ref) and (ref), to come.
Our asymptotic convergence results (Theorems (ref) and (ref) below) require one additional regularity condition, an eigenvalue condition which does not require much more than Assumptions (R2) and (R5).
Assumption (R6): The minimum eigenvalue of
is greater than or equal to a fixed \(\delta > 0\) for all \(N \) sufficiently large, and
Further, for all \(N\) and all \(r \in \{0\} \cup \mathcal{R}\), all of the eigenvalues of \( \mathrm{Cov} ( \boldsymbol{X}_i \mid W_i = r)\) are bounded between \(\lambda_{\text{min}} > 0\) and \(\lambda_{\text{max}} < \infty\).
The assumption that the minimum eigenvalue of (ref) is bounded away from 0 is fairly mild since we already require \(\boldsymbol{\hat{\Sigma}}(\boldsymbol{Z})\) to be full rank and for its minimum eigenvalue to not vanish too quickly, and the minimum singular value of \(\boldsymbol{D}_N\) is bounded from below by a positive constant; see Lemma (ref) in the appendix. Further, this assumption is stronger than necessary for the sake of interpretability; we only require a bound on the minimum eigenvalue of a different matrix (ref) for all \(N\). This condition holding for (ref) is sufficient but not necessary.
Condition (ref) is another assumption that requires \(e_{1N}\) to not vanish too quickly, which is qualitatively similar to Assumptions (R2) and (R5). Holding \(s = s_N\) fixed (as is required in Theorems (ref) and (ref)) and ignoring the distinction between almost sure convergence and convergence in probability, (ref) is more restrictive on the growth rate of \(\lambda_N\) than (R2) since (R2) requires \[ \frac{\lambda_N }{N e_{1N} } \to 0 \] while (ref) requires \[ \frac{\lambda_N}{ \sqrt{N} e_{1N}} \to 0 . \] However, if \(e_{1N}\) is bounded away from 0 with high probability, using a similar argument to the one in the proof of Lemma (ref) in the appendix, (ref) is essentially equivalent to the other requirement in (R2) that \(\lim_{N \to \infty} \lambda_N \sqrt{s_N/N} = 0\) (in the setting of Theorem (ref) where \(s_N = s\) is fixed). If the smallest eigenvalue of \(\mathbb{E} \left[ \boldsymbol{\hat{\Sigma}}(\boldsymbol{Z}) \right]\) is bounded away from 0, one can show that the asymptotic distribution for the minimum eigenvalue of \( \boldsymbol{\hat{\Sigma}}(\boldsymbol{Z}) \) has support bounded away from 0 in our setting where \(\boldsymbol{X}_i\) has finite fourth moments and \(p_N \leq NT\) marchenko1967distribution. If we also assume that \(\boldsymbol{X}_i\) is subgaussian (a weaker assumption than requiring \(\boldsymbol{X}_i\) to have bounded support), under Theorem 5.39 from vershynin_2012 \(e_{1N}\) is bounded away from 0 with high probability for finite \(N\) large enough relative to \(p_N\).
Next we show that under Assumptions (R1) - (R6) the estimators (ref) and (ref) are both oracle estimators fan2001variable, fan2004nonconcave, converge at a \(1/\sqrt{N}\) rate (faster than the guarantee from Theorem (ref)), and are asymptotically Gaussian.
Theorems (ref) and (ref) leverage novel (to the best of our knowledge) extensions of Theorem 2 from kock2013oracle that establish finite-sample variance estimators, including when the weights in a linear combination of the estimated coefficients are themselves estimated; see Theorem (ref) in the appendix.
Theorem (ref) (and, more directly, Theorem (ref)(d) and (f) in the appendix) show that even if \(p_N \to \infty\) FETWFE converges at a \(1/\sqrt{N}\) rate and its asymptotic covariance matrix depends only on the \(s\) parameters of the correct model. That is, FETWFE estimates the model with the same asymptotic efficiency as an OLS-estimated ETWFE estimator that knows the \(p_N - s\) correct restrictions. FETWFE leverages the possible restrictions in the ETWFE estimator to improve efficiency without requiring the practitioner to make assumptions on which restrictions to choose, which could induce asymptotic bias if the wrong restrictions are chosen or compromise asymptotic efficiency if not enough restrictions are chosen. This addresses the need to consider both bias and variance when estimating treatment effects gelman2021slamming.
In Remark (ref) we mentioned that the inclusion of the parameters (ref) in our regression allows for unbiased estimation of treatment effects under both (CCTS) and unconfoundedness of the untreated potential outcomes, but under unconfoundedness the parameters equal exactly 0 and are inefficient to include. However, Theorem (ref) shows that when coefficients that equal 0 are included, FETWFE sets these coefficients equal to exactly 0 with high probability and estimates the nonzero coefficients with the same asymptotic efficiency as if those coefficients were omitted from the model. FETWFE is therefore a good model to choose if one wants to maximize robustness without compromising asymptotic efficiency. The same is true if (CIUN) holds and the regression estimands \(\boldsymbol{\xi}_t^*\) all equal \(\boldsymbol{0}\), or if any irrelevant covariates are included in the regression. So practitioners do not need to hesitate to include a possibly relevant covariate that could be helpful for the plausibility of Assumption (CCTSB) or (CCTS) as long as \(\boldsymbol{Z}\) has full column rank.
Theorem (ref) is conceptual and yields insight on the kind of behavior we can expect from FETWFE in large samples. The following result, which has the same assumptions as Theorem (ref), is more practical and provides explicit formulas for variance estimators for asymptotically valid confidence intervals of treatment effect estimates.
Theorem (ref) shows that we can form asymptotically valid \(1 - \alpha\) confidence intervals with estimator \(\hat{T}_N\) and variance estimator \(\hat{v}_N\) with
where \(\Phi(\cdot)\) is the distribution function of the standard normal distribution.
Theorem (ref)(b) requires data splitting for asymptotic normality of the estimator (ref). However, the data used for estimating the marginal cohort probabilities do not need to include observed responses. In some domains like public health and medicine where difference-in-differences has many potential uses rothbard2023tutorial, doi:10.1146/annurev-publhealth-040617-013507, such unlabeled data may be more widely available than labeled data Candes2018. In a setting like this, asymptotic normality can be achieved without reducing the amount of data used to estimate the model. (It is straightforward to adjust the variance estimator appropriately if the sample sizes of the two data sets are not equal.)
Theorem (ref)(c) does not require sample splitting, at the price of sacrificing asymptotic normality for mere asymptotic subgaussianity and necessitating a conservative variance estimator. Because a subgaussian random variable can be defined as a variable \(X\) satisfying \(\mathbb{P}(|X| \geq x) \leq 2 \exp\{-x^2/K^2\}\) for some \(K > 0\) for all \(x \geq 0\) vershynin2018high, the tails of such a variable decay at least as quickly as a Gaussian random variable, so we might expect \(95\%\) confidence intervals to have similar asymptotic validity to the confidence intervals from part (b). In practice, we show in the simulation studies in Section (ref) that these confidence intervals do seem to retain validity.
We cannot test the null hypothesis of treatment effects equalling 0 under Theorem (ref), since these test statistics are not asymptotically normal if the estimated coefficients all equal 0; instead they converge in probability to 0 (even after scaling by \(\sqrt{NT}\)).
Theorem (ref) in the appendix is a generalization of both Theorems (ref) and (ref) for the classes of conditional average treatment effects estimators (ref) and (ref) in the fixed \(d_N = d\) setting under the assumption that an asymptotically Gaussian estimator of the generalized propensity scores (like the multinomial logit or proportional odds models) is available and an extra split of the data is available to estimate the cohort covariate means.
To test the efficacy of FETWFE under our assumptions, we conduct simulation studies in R using the fetwfe package Faletto:2025aa and the simulator package bien2016simulator. We choose parameters that bear resemblance to the empirical application from Section III of goodman2021difference (\(N = 51, T = 33, R = 12, d_N = 3\)), which we explore in Section (ref). We use fewer cohorts for simplicity of presentation, more covariates to increase \(p_N\) in compensation for the smaller \(R\), and more units to avoid realizations where no units are assigned to a cohort by random chance. We generate data with \(N = 120\) units, \(T = 30\) time periods, \(R = 5\) cohorts entering at times \(\{2, \ldots, 6\}\), and \(d_N = 12\) features, which results in a total of \(p_N = 2209\) covariates and \(NT = 3600\) observations. To generate a \(\boldsymbol{\beta}_N^*\) such that \(\boldsymbol{D}_N \boldsymbol{\beta}_N^*\) is sparse, we generate a single random sparse \(\boldsymbol{\theta}_N^* \in \mathbb{R}^{p_N}\), then transform this into a coefficient vector \(\boldsymbol{\beta}_N^* = \boldsymbol{D}_N^{-1} \boldsymbol{\theta}_N^*\) to use across all simulations. We generate \(\boldsymbol{\theta}_N^*\) by taking a \(p_N\)-vector of all 0 entries and setting each entry equal to 2 randomly with probability 0.1. We set the sign of each term randomly, but since each entry of \(\boldsymbol{\beta}_N^*\) is a sum of terms in \(\boldsymbol{\theta}_N^*\), we bias the signs to avoid treatment effects that are too close to 0: individual nonzero terms in \(\boldsymbol{\theta}_N^*\) are positive with probability 0.6 and negative otherwise.
On each of 700 simulations, we randomly generate \(N\) independent realizations of \(\boldsymbol{X}_{(i)} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{I}_{d_N})\) for the time-invariant covariates. We randomly assign treatments with probabilities \((\pi_0, \pi_1, \ldots, \pi_R) = (1/(R+1), \ldots, 1/(R+1))\), where
are the marginal probabilities of treatment assignments. To ensure the model is estimable, in the rare instances where there are no untreated units or one cohort has no observed units, we draw another set of random assignments.
After drawing covariates and treatment assignments, we construct \(\boldsymbol{\tilde{Z}}\). Then we draw \(N\) independent random effects \(c_i \sim \mathcal{N}(0, 5)\) as well as \(NT\) independent noise terms \(u_{(it)} \sim \mathcal{N}(0, 5)\) to generate \(\boldsymbol{\epsilon}_{(i\cdot)} = c_i \boldsymbol{1}_T + \boldsymbol{u}_{(i\cdot)}\) and finally generate \(\boldsymbol{\tilde{y}} = \boldsymbol{\tilde{Z}} \boldsymbol{\beta}_N^* + \boldsymbol{\epsilon}_{(\cdot \cdot)}\).
Besides FETWFE, we also consider three competitor methods. The first is ETWFE, estimated by OLS on \(\boldsymbol{\tilde{Z}}\). We also consider bridge regression on \(\boldsymbol{Z}\) directly penalizing \(\boldsymbol{\beta}\) rather than penalizing \(\boldsymbol{D}_N \boldsymbol{\beta}\) (BETWFE). For both FETWFE and BETWFE, we follow the experiments of kock2013oracle and use \(q = 0.5\), selecting the penalty \(\lambda_N\) over a grid of 100 values (equally spaced on a logarithmic scale) using BIC. We implement bridge regression for both BETWFE and FETWFE using the R grpreg package breheny2009penalized. Finally, we also consider a slightly more flexible version of (ref), where we estimate the model \[ \tilde{y}_{it} = \nu_r + \gamma_t + \boldsymbol{X} \boldsymbol{\kappa} + \sum_{r \in \mathcal{R}} \tau_r \cdot \mathbbm{1}\{W_i = r \} \mathbbm{1}\{t \geq r\} + \epsilon_{it} \] by OLS (TWFE_COVS). (Notice that adding unit-specific covariates to Model (ref) induces collinearity. Fitting the model instead on cohort fixed effects, like ETWFE, also matches, e.g., Equation 3.2 in callaway2021difference. See also the discussion around Equation 2.5 in sant2020doubly.) The inclusion of covariates might make it seem more plausible that Assumptions (CNAS) and (CCTS) could hold in this model, though wooldridge2021two points out that adding time-invariant covariates does not change the estimates of the treatment effects. Separate treatment effects for each cohort, while insufficient to avoid bias and capture the variation in treatment effects, should improve estimation relative to (ref) and allows for the estimation of cohort-specific treatment effects to compare to the other methods.
We start by evaluating the error of each method in estimating the average treatment effect (ref). The assumptions of Theorem (ref) are satisfied, so we expect FETWFE to estimate the treatment effects more accurately than ETWFE, which will suffer from high variance; BETWFE, which does not assume the correct form of sparsity; and TWFE_COVS, which is asymptotically biased. We estimate the average treatment effect using the FETWFE estimator (ref) with \(\psi_{rt}\) from (ref), as well as analogous estimators using the competitor estimated regression coefficients and observed cohort counts.
We calculate the squared error of each average treatment effect estimate for each method on each iteration. Boxplots of the results are displayed in Figure (ref). We also provide the means and standard errors for the squared errors of each method in Table (ref). Table (ref) contains \(p\)-values for paired one-tailed \(t\)-tests of the alternative hypothesis that the squared error for FETWFE is less than the squared error for each competitor method; all results show significantly better performance for FETWFE at the 0.05 significance level.
These results show that even if most of the restrictions we would naturally consider hold and we are only interested in overall average treatment effect (ref), a simple model like TWFE_COVS is an unsuitable estimator. This aligns with the conclusions of previous works, as we discussed in the introduction. We see that ETWFE has better estimation error, though in practice it still does not perform as well as FETWFE because it estimates about ten times as many parameters as are needed, leading to imprecise estimates.
In Appendix (ref), we also present estimation error results for the cohort-specific average treatment effects \( \tau_{\text{ATT}} (r)\) as defined in (ref) for each of the five cohorts, and estimation error for the coefficients of \(\boldsymbol{\beta}_N^*\) corresponding to the interaction effects, \(\boldsymbol{\rho}^*\). These results show that FETWFE outperforms the competitor methods for these estimands as well.
Next, we examine the extent to which we can trust FETWFE to select the correct restrictions, as we know it does with asymptotically high probability under Theorem (ref). In the same experiment as before, on each simulation we calculate the percentage of restrictions FETWFE successfully identifies among the treatment parameters. In particular, we calculate the percentage of entries of \(\boldsymbol{D} \boldsymbol{\beta}_N^*\) that FETWFE successfully estimates as either zero or nonzero.
On average across the 700 simulations, FETWFE successfully decides whether to use a restriction or not in \(95.0\%\) of the \(p_N\) possible cases, with a standard error of \(0.025\%\). Of the true restrictions (that is, the entries of \(\boldsymbol{D} \boldsymbol{\beta}_N^*\) that equal 0), on average \(97.5\%\) (\(0.021\%\)) are correctly identified. This suggests that FETWFE is effective at identifying restrictions in practice, which contributes to FETWFE's success in estimating treatment effects that we saw in the previous section. (Figure (ref) in Appendix (ref) shows a boxplot of the proportion of restrictions successfully identified across all 700 simulations.)
To get a better sense of the asymptotics, we conduct an additional simulation study in the same way as the first one, but with \(T = 5\),\footnote{Notice that increasing \(T\) increases the number of observations while also increasing the number of treatment effects and interactions to estimate. So we might expect larger \(T\) to have a roughly neutral effect on the convergence properties, while greatly increasing the computational cost of simulations with large \(N\).} \(R = 3, d_N = 2\) (resulting in \(p_N = 50\)), and \(N = 1200\) observations (resulting in 300 observations per cohort in expectation). The sparse \(\boldsymbol{\theta}_N^*\) is generated as before, but with entries each equalling 0 with probability 0.5. Again, we conduct 700 simulations. In this second simulation study, FETWFE successfully decides whether to use a restriction or not in \(99.3\%\) (\(0.042\%\)) of the \(p_N\) possible cases across the 700 simulations, and on average FETWFE identifies \(98.5\%\) (\(0.096\%\)) of the true restrictions.
Finally we investigate the asymptotic distribution of the FETWFE estimators. On each simulation, we estimate nominal asymptotic \(95\%\) confidence intervals for the cohort average treatment effects using (ref) and the variance estimator from Theorem (ref)(a). In the first simulation study, we find coverage rates of 0.791 for the \(r = 2\) cohort (standard error 0.015), 0.884 for the \(r = 3\) cohort (standard error 0.012), \(0.854\) (0.013) for the \(r = 4\) cohort, \(0.803\) (0.015) for the \(r = 5\) cohort, and \(0.850\) (0.014) for the \(r = 6\) cohort. We see that the coverage is unsatisfactory for FETWFE in this relatively small sample size, where in expectation 20 units are observed in each cohort.
We also estimate conservative confidence intervals for the average treatment effect with the conservative variance estimator from Theorem (ref)(c), yielding a coverage rate of 0.993 (0.001). That is, the conservative confidence intervals seem to achieve the desired coverage level in practice, even though the asymptotic distribution of the test statistic is not known to be Gaussian (merely subgaussian).
Lastly, we examine the non-conservative estimator from Theorem (ref)(b): on each simulation study, we generate \(N\) independent unlabeled samples where only cohort assignments (drawn from the same distribution as the \(N\) labeled units) are observed, and we use these to estimate asymptotically exact confidence intervals for the average treatment effect. The resulting coverage is 0.986 (0.001). The coverage may be better for the average treatment effect than the cohort average treatment effects because the number of observations that are relevant for estimating the overall average treatment effect is much larger than that for each of the cohort average treatment effects.
In the second simulation study, described previously in Section (ref), we find coverage rates of 0.944 (0.009) for the \(r = 2\) cohort, 0.941 (0.009) for the \(r = 3\) cohort, and 0.930 (0.010) for the third cohort. The conservative confidence intervals for the average treatment effects have coverage 0.997 (0.001), and the asymptotically exact confidence intervals have coverage 0.973 (0.002). Overall, the coverage rates are close to the nominal \(95\%\) level in this finite \(N\) setting. This is especially promising because we are using heuristics to solve the FETWFE optimization problem, as discussed in Remark (ref).
We use the R fetwfe package Faletto:2025aa to analyze the data from stevenson2006bargaining, as collected in the divorce data set from bacondecomp. stevenson2006bargaining estimate the causal effect of unilateral or “no-fault" divorce laws on suicide rates among women. The data set consists of observations from \(N=51\) states (including the District of Columbia) for \(T = 33\) years, along with demographic covariates. The response variable is the aggregate state suicide rate among women by year (as in the first column of Table 1 of stevenson2006bargaining). See stevenson2006bargaining and goodman2021difference for more details on the data.
We remove states treated at time \(T = 1\) from the data, leaving \(N = 42\) observations. 5 states are never treated, and \(R = 12\) cohorts are treated at various times. Similarly to the analysis in goodman2021difference, we consider the time-varying covariates murderrate (state female homicide rate), lnpersinc (natural log of state personal income per capita), and afdcrolls (state welfare participation rate), including the values from time \(t = 1\) as time-invariant controls. However, murderrate is missing for some observations at \(t = 1\) so is omitted, leaving \(d_N = 2\) covariates. In total, we have \(NT = 1386\) observations and \(p_N = 908\) coefficients to estimate.
We estimate the variances \(\sigma_c^2\) (the variance of the unit random effects) and \(\sigma^2\) (the variance of the added noise) using the estimators from pesaran-2015-text with ridge regression. We use FETWFE with \(q = 0.5\) to estimate the marginal average treatment effects, again leveraging the R grpreg package as in the simulation studies and selecting \(\lambda_N\) by BIC. We estimate that the overall average percent change to the suicide rate for women is \(-3.76\%\), which is similar to the estimate of \(-2.52\%\) (standard error \(1.09\%\)) that goodman2021difference calculates when using similar controls. Applying Theorem (ref)(c), we estimate a conservative standard error of \(4.58\%\), yielding a conservative \(95\%\) confidence interval of \((-12.73\%, 5.21\%)\) for the average treatment effect.
In Table (ref) we present point estimates for the estimated average treatment effects for each of the cohorts. (We omit standard errors since the individual cohorts have very small sample sizes.) We estimate that four cohorts have non-zero treatment effects, with the 1970 cohort (consisting of California and Iowa) having by far the largest treatment effect estimate. We emphasize that these individual estimates are likely less reliable than the average treatment effect estimate since most cohorts have a small number of units.
In this work, we have proposed the fused extended two-way fixed effects estimator for difference-in-differences with staggered adoptions. Our estimator leverages the fact that the extended two-way fixed effects estimator includes coefficients that are likely to be similar (for example, because of their proximity in time), but it selects appropriate restrictions automatically rather than relying on hand-selected restrictions, which could be biased. FETWFE retains the asymptotic unbiasedness (and normality) of ETWFE but leverages sparsity for large efficiency gains. In particular, we showed that our estimator is consistent, selection consistent, asymptotically Gaussian, and has an oracle property for several classes of heterogeneous marginal and conditional average treatment effects. Finally, we demonstrated the usefulness of FETWFE in practice in both simulation studies and an empirical application.
Throughout the paper we have mentioned places where our work might be relatively straightforward to extend. Here we specify some other possible extensions. It would be useful to develop theory allowing for FETWFE to be estimated when the design matrix \(\boldsymbol{Z}\) is deficient in column rank; for example, when \(p_N > NT\). One could extend our method to this setting by suitably adapting the marginal bridge estimator of kock2013oracle, which involves a screening step to remove irrelevant covariates before estimating a regression on the reduced design matrix in a second step. This allows for consistency and asymptotic normality of the second-step estimator under similar assumptions to those that we have used in this work. Alternatively, consistency (but not asymptotic normality) could be established under weaker assumptions in the \(q = 1\) case by suitably adapting theory for lasso estimators like Corollary 2 from negahban2012unified.
Another potentially useful extension would be relaxing the parallel trends assumption to allow for heterogeneous trends as described in wooldridge2021two. This requires adding extra parameters to the ETWFE estimator, but FETWFE is well-poised to adapt to this setting because it can automatically restrict to 0 any parameters that are added unnecessarily.
I thank Jacob Bien for very helpful discussions that shaped the direction of this work and edits to an early draft.
This work was supported by the University of Southern California.