Better Understanding Triple Differences Estimators

bibunit\begin{abstract} Triple Differences (DDD) designs are widely used in empirical work to relax parallel trends assumptions in Difference-in-Differences (DiD) settings. This paper highlights that common DDD implementations---such as taking the difference between two DiDs or applying three-way fixed effects regressions---are generally invalid when identification requires conditioning on covariates. In staggered adoption settings, the common DiD practice of pooling all not-yet-treated units as a comparison group can introduce additional bias, even when covariates are not required for identification. These insights challenge conventional empirical strategies and underscore the need for estimators tailored specifically to DDD structures. We develop regression adjustment, inverse probability weighting, and doubly robust estimators that remain valid under covariate-adjusted DDD parallel trends. For staggered designs, we demonstrate how to effectively utilize multiple comparison groups to obtain more informative inferences. Simulations and three empirical applications highlight bias reductions and precision gains relative to standard approaches. A companion R package is available. JEL: C10; C14; C21; C23. \\ Keywords: Triple Differences; Difference-in-Differences; Difference-in-Difference-in-Differences; Parallel Trends; Doubly Robustness; Staggered Adoption. \end{abstract} \sloppy \section{Introduction} Over the last few years, we have seen a big Difference-in-Differences (DiD) “methodological revolution” with multiple DiD estimators being proposed to address the interpretability shortcomings associated with using simple two-way fixed-effects specifications in the presence of treatment effect heterogeneity.\footnote{See, e.g., Roth2023a and Baker_etal_2025_JEL for overviews.} Although these modern DiD estimators can capture richer notions of heterogeneity, in practice, they rely on parallel trends (PT) assumptions, and a concern relates to how plausible these PT assumptions are. When such PT assumptions are not reasonable enough approximations of reality, one may doubt the conclusions of DiD studies rambachan_roth, Chiu_etal_2025. In some setups, however, it is possible to naturally relax such DiD-type PT assumptions and retain the simplicity and empirical appeal of DiD-type analysis. This is particularly the case when a unit needs to fulfill two criteria to be treated, e.g., it belongs to (i) a group (e.g., a state) in which the treatment is enabled, and (ii) a partition of the population that qualifies (or is eligible) for treatment (e.g., women). Such setups are often referred to as Triple Differences (DDD), and allow for group-specific and partition-specific violations of parallel trends. Since its introduction by Gruber1994, DDD has gained considerable popularity among empirical researchers. Some prominent recent DDD applications include antwi_effects_2013, Walker_2013_QJE, garthwaite_public_2013, Alsan_Wanamaker_2018_QJE, Patnaik_2019, hansen_national_2023, and Bailey_Helgerman_Stuart_2024_QJE---see Olden2022 for additional documentation of DDD applications, and Section (ref) for three DDD applications. Despite its empirical popularity, little attention has been devoted to better understanding the econometric foundations of DDD setups with covariates, multiple periods, and/or staggered treatment adoption. This article aims to improve our understanding of Triple Difference (DDD) designs. Our main goal is to provide a set of clear, easy-to-use, and theoretically grounded tools that empirical researchers can utilize whenever they wish to explore DDD designs. To that end, we study identification, estimation, and inference procedures for DDD when covariates may be important for the reliability of the identification assumptions, multiple periods are available, and treatment adoption is potentially staggered over time. We tackle the DDD problem using causal inference first principles and uncover interesting results that challenge some conventional wisdom and common practices. For instance, although DDD estimators can be understood as the difference between two DiD estimators in setups with two periods and no covariates Olden2022, our results highlight that this is no longer the case when covariates are required to justify the plausibility of a DDD-type parallel trends assumption. As we illustrate through simulations, proceeding as if DDD were just the difference between two DiD estimators can lead to severely biased results. Such bias arises because this naive DDD strategy fails to integrate the covariate distribution over the correct reference group---the treated group. We show that it is straightforward to avoid these problems and propose regression adjustment, inverse probability weighting, and doubly robust DDD estimators. We also show that if one wants to cast DDD in terms of DiD, one would need three---and not two---DiD terms. Each of these DiD terms compares effectively treated units with a different type of untreated units, e.g., units in a treated state but ineligible for treatment, units that are eligible but are in an untreated state, or ineligible units in an untreated state. In setups with staggered treatment adoption, our results once again challenge the interpretation of DDD as the difference between two DiDs. In DiD with staggered treatment adoptions, it is now common to pool all not-yet-treated units at a time period and use that aggregate set of units as a valid comparison group.\footnote{See, e.g., Callaway2021, deChaisemartin2020, Wooldridge2021, Borusyaetal2024.} Thus, one may think that a similar strategy should work with DDD. However, our results suggest that this is generally not the case and that pooling all not-yet-treated units and proceeding as in staggered DiD procedures can lead to biased estimators for average treatment effect parameters, even when covariates do not play a significant role. This arises because the proportion of units eligible for treatment may change across groups that enable treatment over time. As DDD allows for group-specific and partition-specific violations of DiD-type PT, these differential trends do not average out when pooling all not-yet-treated units, leading to potentially misleading estimates. We propose DDD estimators that bypass this drawback by using any specific not-yet-treated unit as a comparison group (e.g., the set of units in groups that never enabled treatment). As one can potentially use different not-yet-treated cohorts as comparison groups, we also discuss combining these to form more precise estimators. Our proposed DDD estimator that aggregates across different comparison groups can be understood as a two-step Generalized Method of Moments (GMM) procedure based on recentered influence functions. Importantly, our staggered DDD procedures can flexibly accommodate covariates using regression adjustment, inverse probability weighting, or doubly robust methods and can also be used to form event-study estimators that highlight how average treatment effects evolve with elapsed treatment time. We illustrate our results through Monte Carlo simulations and three different empirical applications that address various DDD setups. Our Monte Carlo results highlight that ignoring the key takeaways in this paper and relying on overly rigid estimators can lead to biases and imprecise conclusions; our proposed DDD estimators bypass these limitations. In terms of empirical applications, we (a) revisit cai_insurance_2016 and analyze the effects of agricultural insurance programs on financial decisions in China; (b) build on carbon_pricing and examine the impact of the emission trading scheme on carbon emissions in China; and (c) hansen_national_2023 and assess the impact of genetically modified crop adoption on countrywide yields. Compared to the three-way fixed effects estimators used by cai_insurance_2016, our doubly robust DDD estimator yields substantial gains in precision, with their confidence intervals being up to 115% wider than ours. In our application examining the effect of the emission trading scheme on carbon emissions with staggered adoption, our doubly robust DDD estimators yield a more modest and statistically insignificant effect. In contrast, the three-way fixed estimators indicate a statistically significant effect on the share of low-carbon patents. When we apply our staggered DDD tools to hansen_national_2023's data, we find that their conclusions are robust to dropping observations that are part of the never-enabling crop-country group, which is not the case when using more standard three-way fixed effects estimators. Taken together, these results highlight that our proposed tools can indeed lead to interesting new insights that are of practical interest. Related literature: This article contributes to the rapidly expanding literature on DiD-related methods. In particular, we contribute to the scarce literature on DDD procedures. Our paper is related to Olden2022, though we cover substantially more general DDD setups with (a) multiple periods, (b) staggered treatment adoption, and (c) covariates potentially playing an important role for the plausibility of the identification assumptions. In this sense, our paper can be understood as the DDD analog of Callaway2021. However, and in sharp contrast with Callaway2021, our DDD procedures cannot pool all not-yet-treated units as an aggregated comparison group, highlighting some interesting differences between staggered DiD and DDD designs. Our paper is also related to strezhnev2023, who introduced a decomposition of the DDD estimators based on three-way fixed effects specifications, illustrating when and why it fails to recover an easily interpretable causal parameter of interest when treatment effects are heterogeneous. To some extent, strezhnev2023 can be understood as the analog of Goodman2021 to DDD setups. As such, our results complement strezhnev2023, and since our estimators do not rely on a rigid three-way fixed effect specification, they circumvent the issues highlighted in his paper. Our paper is also related to Sloczyski2022,Sloczyski2024_IVLATE in the sense that our proposed tools avoid issues related to potentially misleading weights related to model misspecifications. In this paper, we use the term “triple differences” to describe designs under which units must satisfy two criteria to be (effectively) treated. However, we note that sometimes different researchers use the term triple differences more broadly, for instance, when they are interested in analyzing treatment effect heterogeneity across subgroups. In such cases, it is important to note that the underlying identification assumptions and the parameters of interest would differ from those studied in this paper; see Caron2025 for a recent example. Those procedures should be understood as a complement to the ones we discuss in this paper, as they can be used to answer different questions of interest. Organization of the paper: The rest of the paper is organized as follows. In the next section, we present our framework. In Section (ref), we challenge some standard empirical practices for DDD analyses in terms of interpreting them as a simple extension of DiD analysis, and we also highlight some important practical takeaway messages from our main results. Section (ref) introduces our formal identification, estimation, and inference results. Section (ref) presents a Monte Carlo study to demonstrate the finite sample properties of our estimator, while Section (ref) presents three empirical illustrations. Section (ref) concludes. Detailed mathematical proofs and additional results can be found in the Supplemental Appendix. To ease the adoption of our proposed DDD tools, we provide an open-source \texttt{R} package, \texttt{triplediff}, which automates all the procedures described in this article. \section{Framework} We start our analysis by discussing the specifics of our DDD research design, including potential outcomes, parameters of interest, and identification assumptions. We consider a setup with $T$ time periods, $t=1,2,\dots, T$. Units are indexed by $i$, with $i=1,2,\dots, n$. We focus on setups where $n$ is much larger than $T$, as our inference procedures are asymptotically justified using the “fixed-$T$, large-$n$” panel data framework. Each unit may be exposed to a binary treatment in any time period $t>1$. Treatment is an absorbing state such that once a unit is treated, it remains treated for the remainder of the panel. Each unit $i$ belongs to a group (e.g., a state or a country) that enables treatment for the first time in period $g>1$. Let $S_i \in \mathcal{S} \subseteq \{2,...,T\} \cup \{\infty\}$ be a variable that indicates the first time the policy/treatment was enabled, with the notion that $S = \infty$ if the policy is not enabled by $t=T$. In addition, each unit belongs to a population partition that qualifies (or is eligible) for the treatment or not (e.g., being a woman, or an indicator for specific crops). We denote this variable by $Q_i$ with $Q_i=1$ if unit $i$ (eventually) qualifies for the treatment and $Q_i=0$ otherwise. For simplicity, we assume that this population partition that eventually qualifies for treatment is time-invariant. In our DDD setup, a unit $i$ is treated in period $t$ if $t \ge S_i$ and $Q_i=1$, i.e., if it belongs to a group that has already enabled treatment by period $t$ (i.e., $t \ge S_i$) and it qualifies for treatment (i.e., $Q_i=1$). With this notation that makes it clear that a unit $i$ is treated if it satisfies two criteria, let $D_{i,t} = 1\{t\ge S_i, Q_i = 1\}$ be an indicator for whether unit $i$ receives treatment in period $t$, and let $G_i = \min\{t : D_{i,t} = 1\}$ be the earliest period at which unit $i$ has received treatment. If $i$ is never treated during the sample, then $G_i = \infty$. Here, we have that $G_i = S_i$ if $Q_i = 1$ and $G_i = \infty$ if $Q_i = 0$.\footnote{Note that when all units are eligible for treatment, we have $G_i = S_i$, getting us back to a (staggered) DiD setup; see, e.g., Callaway2021 and Sun2021.} Let $\mathcal{G}$ denote the support of $G_i$ and $\mathcal{G}_{\text{trt}} = \mathcal{G}\setminus \{\infty\}$. We assume that a group of “never-enabled” units always exists, i.e., $S_i = \infty$ for some units. In an application where all units belong to a group that eventually enables treatment, we remove all that data from all units from the time the last cohort enabled treatment onwards, i.e., we drop all observations from periods $t>\max S_i$, and retain the remaining data as the “effective” data to be used in our analysis, where the last-to-be-eligible group becomes the “never-eligible” group.\footnote{If needed, we also update the notion of support of all variables to reflect this change. In particular, $T$ here denotes the number of available periods in the subset of the data that we will use in our analysis.} Finally, we also assume that a vector of pre-treatment covariates $X_{i}$, whose support is denoted by $\mathcal{X}\subseteq \mathbb{R}^d$ is available. Regarding potential outcomes, we adopt the potential outcome framework of Robins1986 with potential outcomes indexed by treatment sequences. Let $\mathbf{0}_s$ and $\mathbf{1}_s$ be $s$-dimensional vectors of zeros and ones, respectively, and denote the potential outcome for unit $i$ at time $t$ if unit $i$ is first treated at time $g$ by $Y_{i,t}(\mathbf{0}_{g-1},\mathbf{1}_{T-g+1})$, and denote by $Y_{i,t}(\mathbf{0}_{T})$ the outcome if untreated by time $t=T$. As we focus our attention on staggered treatment adoptions, we can simplify notation and index potential outcomes by the time treatment begins, $g$: $Y_{i,t}(g) = Y_{i,t}(\mathbf{0}_{g-1},\mathbf{1}_{T-g+1}) $ and use $Y_{i,t}(\infty) = Y_{i,t}(\mathbf{0}_{T})$ to denote never-treated potential outcomes. In practice, though, we observe, \begin{equation} Y_{i,t} = \sum_{g \in \mathcal{G}} 1\{G_i = g\} Y_{i,t}(g), \end{equation} where $1\{A\}$ represents the indicator function, which equals one if $A$ is true and zero otherwise. Additionally, we assume the observation of a random sample of $(Y_{t=1}, \dots, Y_{t=T}, X', G, S, Q)'$. \begin{namedassumption}{S}[Random Sampling] $\{(Y_{i,t=1},\dots,Y_{i,t=T}, X_i', G_i, S_i, Q_i)'\}_{i=1}^n$ is a random sample from $(Y_{t=1},\dots,Y_{t=T}, X', G,S,Q)'$. \end{namedassumption} \subsection{Parameters of interest} In this paper, we aim to gain a deeper understanding of how average treatment effects vary across periods and different groups defined by the treatment starting period. More specifically, we want to make inferences on functionals of the group-time average treatment effects, $ATT(g,t)$'s, defined as \begin{align} ATT(g,t) \equiv \mathop\!\mathbb{E}[Y_{i,t}(g)-Y_{i,t}(\infty)|G_i=g] = \mathop\!\mathbb{E}[Y_{i,t}(g)-Y_{i,t}(\infty)|S_i=g, Q_i = 1], \end{align} By exploring that in our context $G_i=g$ if and only if $S_i=g$ and $Q_i = 1$, we have that $ATT(g,t) = \mathop{}\!\mathbb{E}[Y_{i,t}(g)-Y_{i,t}(\infty)|S_i=g, Q_i = 1]$. Note that $ATT(g,t)$ captures how the average treatment effects evolve over time for each treatment group $g$ Callaway2021. As such, one can use $ATT(g,t)$ to construct group-$g$-specific event studies by analyzing how average treatment effects vary with elapsed treatment time $e=t-g$. In some setups with multiple groups $g$, researchers may want to summarize over the many $ATT(g,t)$'s into a more aggregate parameter. A natural summary parameter that still allows one to understand treatment effect dynamics with respect to elapsed treatment time is the aggregated event study parameter $ES(e)$, defined as \begin{align} ES(e) \equiv \mathop\!\mathbb{E}\big[ATT\left( G, G+e\right) \big| G+e \in [2,T]\big] = \sum_{g\in\mathcal{G}_{\text{trt}}} \mathop\!\mathbb{P}(G=g|G+e \in [2,T]) ATT(g,g+e). \end{align} One may also want to aggregate the event study coefficients further to recover a scalar summary measure. Let $\mathcal{E}$ denote the support of post-treatment event time $E=t-G$, $t\geq G$, and let $N_E$ denote its cardinality. Then, \begin{align} ES_{\text{avg}} &\equiv \dfrac{1}{N_E}\sum_{e\in \mathcal{E}} ES(e), \end{align} provides a simple average of all post-treatment event study coefficients. Many other summary parameters are possible; see Callaway2021 for discussions of several alternatives. \subsection{Identification assumptions} To identify the $ATT(g,t)$'s and their functionals $ES(e)$ and $ES_{\text{avg}}$, we impose the following assumptions. \begin{namedassumption}{SO}[Strong Overlap] For every $(g,q) \in \mathcal{S} \times \{0,1\}$ and for some $\epsilon>0$, $\mathbb{P}[S=g, Q=q |X] > \epsilon$ with probability one. \end{namedassumption} Assumption (ref) is an overlap condition that ensures that for any value of $X$, there are units with any combination $(g,q) \in \mathcal{S}$ that have comparable $X$ values. Heuristically, this condition guarantees that we cannot perfectly predict which $(g,q)$-partition a unit belongs to using information from $X$. This assumption also rules out irregular identification Khan2010.\footnote{As our focus is on $ATT(g,t)$-type parameters, it is possible to relax Assumption (ref) to hold only over $X$ in the support of the covariates among the (eventually) treated units. To simplify the discussion, we abstract from these subtle points.} We also impose the following no-anticipation assumptions. \begin{namedassumption}{NA}[No-Anticipation] For every $g \in \mathcal{G}_{\text{trt}}$, and every pre-treatment period $t<g$, $\mathop{}\!\mathbb{E}[Y_{i,t}(g)|S=g, Q = 1,X] = \mathop{}\!\mathbb{E}[Y_{i,t}(\infty)|S=g, Q = 1,X]$ with probability one. \end{namedassumption} Assumption (ref) rules out anticipatory effects among treated units as in, e.g., Abbring2003, Callaway2021, and Sun2021. This assumption is important as it allows us to consider observations in pre-treatment periods $t<g$ as effectively untreated. If units are expected to anticipate some treatments---for example, if treatment is announced in advance---it is important to adjust the definition of the treatment date to account for this; see Malani2015 for a discussion. Next, we impose our final identification assumption that restricts the evolution of average untreated potential outcomes across groups. \begin{namedassumption}{DDD-CPT}[DDD-Conditional Parallel Trends] For each $g \in \mathcal{G}_{\text{trt}}$, $g' \in \mathcal{S}$ and time periods $t$ such that $t\ge g$ and $g'>\max\{g,t\}$, with probability one, \begin{eqnarray*} \mathbb{E}\left[Y_{t}(\infty) - Y_{t-1}(\infty)| S = g, Q=1, X\right] &-& \mathbb{E}\left[Y_{t}(\infty)- Y_{t-1}(\infty) | S = g, Q=0, X \right] \\[-0.2cm] &=&\\[-0.2cm] \mathbb{E}\left[Y_{t}(\infty) - Y_{t-1}(\infty)|S = g', Q=1, X\right] &-& \mathbb{E}\left[Y_{t}(\infty)- Y_{t-1}(\infty) | S = g', Q=0, X \right]. \end{eqnarray*} \end{namedassumption} Assumption (ref) is a conditional parallel trends assumption for DDD setups that generalizes the unconditional DDD parallel trend assumption for the two-period setup of Olden2022 to setups with multiple periods, staggered treatment adoption, and when assumptions are only plausible after conditioning on $X$. Assumption (ref) can also be understood as an extension of the conditional PT assumption based on not-yet-treated units from the DiD setup of Callaway2021 to our DDD setup--- i.e., we can use any unit from groups that either never enabled treatment or those that will eventually enable treatment. Moreover, if covariates do not play any important identification role in the analysis, one can take $X=1$ for all units, so Assumption (ref) would hold unconditionally. Several remarks about Assumption (ref) are worth making. First, if all units in a group $S$ are eligible for treatment, Assumption (ref) reduces to Assumption 5 of Callaway2021 under the no-anticipation condition in Assumption (ref). However, this case is not appealing to us, as that would not qualify as a DDD design. Second, as Assumption (ref) only holds after conditioning on covariates, it does not restrict the evolution of untreated potential outcomes across covariate-strata, i.e., it allows for covariate-specific trends, which can be very important in applications. Third, and perhaps the most empirically relevant, Assumption (ref) does not impose DiD-type parallel trends among units with $S=g$---i.e., it does not impose that $\mathbb{E}\left[Y_{t}(\infty) - Y_{t-1}(\infty)| S = g, Q=1, X\right] = \mathbb{E}\left[Y_{t}(\infty)- Y_{t-1}(\infty) | S = g, Q=0, X \right]$---nor impose DiD-type parallel trends across treated groups---i.e., it does not impose that $\mathbb{E}\left[Y_{t}(\infty) - Y_{t-1}(\infty)| S = g, Q=1, X\right] = \mathbb{E}\left[Y_{t}(\infty) - Y_{t-1}(\infty)|S = g', Q=1, X\right]$. As such, Assumption (ref) allows for violations of traditional DiD-based PT, as provided that these violations are stable across groups. This observation is arguably what makes DDD appealing in settings where it can be applied. \section{Implications for Empirical Practices} Before presenting our formal results on identification, estimation, and inference for average treatment effects in DDD designs, we challenge some standard empirical practices for DDD analyses and highlight some important practical takeaways from our paper. We first start with a simple setup involving only two periods, $t=1$ and $t=2$, and two eligibility groups, $S_i = 2$ (who enabled treatment in period 2) and $S_i = \infty$ (who have not enabled treatment by period two). As before, units are either eligible ($Q_i =1$) or ineligible ($Q_i = 0$) for the treatment, and we let $D_{i,t}$ be a treatment indicator for unit $i$ in time period $t$, i.e., $D_{i,t} = 1\{t\ge S_i, Q_i = 1\}$. Since there are only two eligibility groups and two time periods, the relevant group-time ATT in such a scenario is $ATT(2,2)$. As discussed in Olden2022, when covariates are not important for the analysis, one can use ordinary least squares (OLS) based on the following three-way fixed effects linear regression specification to recover the $ATT(2,2)$: \begin{align} Y_{i,t}=& \gamma_i + \gamma_{s,t} + \gamma_{q,t} + \beta_{\text{3wfe}} D_{i,t} + \varepsilon _{i,t}, \end{align} where $\gamma_i$ are unit fixed effects, $\gamma_{s,t}$ and $\gamma_{q,t}$ are enabled-group-by-time and qualified-group-by-time fixed effects, and $\beta_{\text{3wfe}}$ is the parameter of interest. Indeed, in this particular setup, under Assumptions (ref), (ref), (ref), and (ref) with $X=1$ a.s., it is straightforward to show that \begin{align} \beta_{\text{3wfe}} =& \Bigg[\underbrace{\bigg(\mathbb{E}\left[Y_{2} - Y_{1}|S = 2, Q=1\right]\bigg)-\bigg( \mathbb{E}\left[ Y_{2} - Y_{1} | S = 2, Q=0 \right] \bigg)}_{\text{DiD estimand among } S=2}\Bigg] \nonumber\\ & - \Bigg[\underbrace{\bigg(\mathbb{E}\left[ Y_{2} - Y_{1}| S = \infty, Q=1 \right] \bigg)\Bigg.-\Bigg.\bigg(\mathbb{E}\left[ Y_{2} - Y_{1}|S=\infty, Q=0 \right]\bigg)}_{\text{DiD estimand among } S=\infty}\Bigg] \\ =& ATT(2,2) \nonumber. \end{align} The observation that $\beta_{\text{3wfe}} = ATT(2,2)$ in this simple setup has two implications: (i) one can use a simple three-way fixed effects (3WFE) regression specification and use OLS to estimate $ATT(2,2)$ in the DDD design, and (ii) DDD estimates can be understood as the difference between two DiD estimates Olden2022. Based on these, one may be tempted to extrapolate these claims to more general setups. In what follows, we highlight that, unfortunately, this is not warranted and that proceeding in this manner can lead to non-negligible biases. The solution to these issues is relatively simple and involves adopting a “forward-engineering” approach to DDD setups Baker_etal_2025_JEL, recognizing its specific characteristics. See also Mogstad_Torgovitsky_2024 for a related discussion in an instrumental variable context. \subsection{DDD setup with two periods, with covariates being important} In this section, we illustrate the challenges of leveraging simple regression-based and DiD tools to DDD setups using simple simulations in a setup where covariates are important for identification, i.e., when Assumption (ref) is satisfied only after accounting for covariates. We consider a design with four different time-invariant, unit-specific covariates, $X_i = (X_{i,1}, X_{i,2}, X_{i,3}, X_{i,4})'$, and two periods and two treatment-enabling groups. The true $ATT(2,2)$ in our simulations is zero. To ease the exposition, we abstract from further details about the DGP and refer the reader to Section (ref) and Supplemental Appendix (ref) for a more detailed discussion. Based on the discussion on DDD without covariates above, it is natural to consider three alternative ways to incorporate covariates in the analysis. The first approach would be to “extrapolate” from (ref), add the interactions of the time-invariant covariates with post-treatment dummies, \begin{align} Y_{i,t}=& \gamma_i + \gamma_{s,t} + \gamma_{q,t} + \tilde{\beta}_{\text{3wfe}} D_{i,t} + (X_i 1_{ \{t=2\}})'\theta + u_{i,t}, \end{align} and interpret the OLS estimates of $\tilde{\beta}_{\text{3wfe}}$ as estimates of $ATT(2,2)$. The second natural way to proceed is similar, but it would leverage the Mundlak device and replace unit fixed effects in (ref) with $S$-by-$Q$ fixed effects, add covariates linearly, \begin{align} Y_{i,t}=& \gamma_{s,q} + \gamma_{s,t} + \gamma_{q,t} + \check{\beta}_{\text{3wfe}} D_{i,t} + X_i'\theta + e_{i,t}, \end{align} and interpret the OLS estimates of $\check{\beta}_{\text{3wfe}}$ as estimates of $ATT(2,2)$. Both strategies leverage a presumption that it is sufficient to add covariates linearly into the 3WFE regression specification (ref) to allow for covariate-specific trends. A third strategy, which is also a priori intuitive, presumes that we can write DDD estimates as the difference between two DiD estimates: one DiD using the subset with $S=2$ and considering units treated if $Q=1$, and another DiD using the subset with $S=\infty$ and considering units treated if $Q=1$. Here, one could consider different estimation strategies. We focus on the doubly robust (DR) DiD estimators proposed by SantAnna2020, as a DR estimator is generally more resilient to model misspecifications. \begin{figure}[!htp!] \begin{center} \begin{subfigure}[t]{0.48\textwidth} \caption{3WFE with covariates interacted with post} \end{subfigure} \begin{subfigure}[t]{0.48\textwidth} \begin{center} \caption{Mundlak-based 3WFE with covariates} \end{center} \end{subfigure} \begin{subfigure}[t]{0.48\textwidth} \caption{Difference between two Doubly Robust DiDs} \end{subfigure} \begin{subfigure}[t]{0.48\textwidth} \begin{center} \caption{Doubly Robust DDD } \end{center} \end{subfigure} \caption{Density of different DDD estimates for ATT(2,2): two-period setup with covariates} \end{center} \justifying \scriptsize{Notes: Simulation designs based on DGP 1 described in Section (ref) and Supplemental Appendix (ref), with $n=5,000$ and $1,000$ Monte Carlo repetitions. True $ATT(2,2)$ is zero and is indicated in the solid vertical line in all panels. Panel (a) displays the density of OLS estimates of $\tilde{\beta}_{\text{3wfe}}$ based on (ref). Panel (b) displays the density of OLS estimates of $\check{\beta}_{\text{3wfe}}$ based on (ref). Panel (c) displays the density of the DDD estimates based on the difference between two doubly robust DiD estimators SantAnna2020. Panel (d) displays the density of the estimates based on our proposed doubly robust DDD estimator described in (ref). All densities are computed across all simulation draws. Panels have the same x-axis range but different y-axis. } \end{figure} To check if such alternative strategies recover the $ATT(2,2)$, we draw $5,000$ units in each simulation draw, compute estimates using these three alternative estimators, and repeat this $1,000$ times---we defer all details of the data generating process to Section (ref) and Supplemental Appendix (ref). Panels (a) and (b) from Figure (ref) display the density of OLS estimates for the $D_{i,t}$ coefficient in the regression specifications (ref) and (ref), while Panel (c) displays the density of the DDD estimates based on the difference between two SantAnna2020 DR DiD estimators. These three panels make it clear that when covariates are necessary to justify the plausibility of the DDD research design, using any of these three procedures can lead to substantial biases and harm policy recommendations and evaluations. In other words, these results highlight that traditional 3WFE linear regression specifications are “too rigid” to be reliable for DDD analysis. They also highlight that, in general, one should not claim that DDD is the difference between two DiD procedures. A natural question that then arises is: What should we do instead? As we discuss in Section (ref), one can form regression adjustment, inverse probability weighting, and doubly robust DDD estimators that do not suffer from the shortcomings highlighted in Panels (a) - (c) in Figure (ref). Among these, we generally favor the DR DDD estimator as it is more resilient against model misspecifications than the other alternatives. To form the DR DDD estimator for $ATT(2,2)$, we need to estimates for the outcome regression models $m_{Y_2-Y_1}^{S=g,Q=q}(X) \equiv \mathbb{E}\left[Y_2 - Y_1|S=g,Q=q ,X\right]$, and for the generalized propensity score model $p^{S=g,Q=q}(X) \equiv \mathbb{P}[S=g, Q=q | X]$. Let $\widehat{m}_{Y_1-Y_0}^{S=g,Q=q}(X)$ and $\widehat{p}^{S=g,Q=q}(X)$ be working models for these---e.g, a linear regression model and a multinomial logistic linear model, though much richer, potentially machine-learning based estimators can also be used Ahrens2025_DML_JEL. Based on these estimates, we propose the following DR DDD estimator for the $ATT(2,2)$: \begin{small} \begin{align} \widehat{ATT}_{\text{dr}}(2,2) =& \mathbb{E}_n\left[ \left(\widehat{w}^{S=2,Q=1}_{\text{trt}}(S,Q) - \widehat{w}^{S=2,Q=0}_{\text{comp}}(S,Q,X) \right)\left(Y_2- Y_1 - \widehat{m}_{Y_2-Y_1}^{S=2,Q=0}(X) \right)\right]\nonumber \\ & +\mathbb{E}_n\left[ \left(\widehat{w}^{S=2,Q=1}_{\text{trt}}(S,Q) - \widehat{w}^{S=\infty,Q=1}_{\text{comp}}(S,Q,X) \right)\left(Y_2- Y_1 - \widehat{m}_{Y_2-Y_1}^{S=\infty,Q=1}(X) \right)\right] \\ &- \mathbb{E}_n\left[ \left(\widehat{w}^{S=2,Q=1}_{\text{trt}}(S,Q) -\widehat{w}^{S=\infty,Q=0}_{\text{comp}}(S,Q,X)\right) \left(Y_2- Y_1 - \widehat{m}_{Y_2-Y_1}^{S=\infty,Q=0}(X)\right)\right] \nonumber, \end{align} \end{small} where $\mathbb{E}_n [A] = n^{-1}\sum_{i=1}^n A_i$ denotes the sample mean, and the estimated weights $\widehat{w}$ are given by \begin{small} \begin{align*} \widehat{w}^{S=2,Q=1}_{\text{trt}}(S, Q)\equiv\dfrac{1{\{S=2,Q=1\}} }{\mathbb{E}_n[1{\{S=2,Q=1\}}]}, \quad \widehat{w}^{S=g, Q=q}_{\text{comp}}(S, Q, X) \equiv \dfrac{\dfrac{1{\{S=g,Q=q\}} \cdot \widehat{p}^{S=2,Q=1}(X) }{\widehat{p}^{S=g,Q=q}(X)}} {\mathbb{E}_n\left[\dfrac{1{\{S=g,Q=q\}} \cdot \widehat{p}^{S=2,Q=1}(X) }{\widehat{p}^{S=g,Q=q}(X)}\right]}. \end{align*} \end{small} Interestingly, it is worth mentioning that although the DR DDD estimator in (ref) cannot be expressed as the difference between two DR DiD estimators, it is a function of \emph{three} DR DiD estimators, each one using a particular subset of the untreated units as a comparison group. For comparisons, we report in Panel (d) of Figure (ref) the density of the estimates using the DR DDD estimates based on (ref). Our proposed DR DDD estimator not only mitigates the biases associated with other estimation strategies but also yields substantially more precise estimates. All in all, the results in Figure (ref) highlight that common DDD practices can lead to misleading conclusions. However, it is straightforward to bypass these limitations by adopting our DR DDD estimators. \subsection{DDD setups with variation in treatment timing} The practical challenges of estimating average treatment effects in DDD setups are not confined to the presence of covariates. Even in designs without covariates, the use of too-rigid 3WFE regression specifications like (ref) can lead to misleading estimates when there is variation in treatment timing across groups strezhnev2023. In such cases, new identification and estimation concerns emerge that the recent DiD literature does not address. In particular, in this section, we highlight that, unlike in staggered DiD procedures like Callaway2021, Borusyaetal2024, and Wooldridge2021\footnote{See also deChaisemartin2020, deChaisemartin2023b-intertemporal-treatments for related procedures that also pool not-yet-treated units. Their estimators allow for treatment turning on and off. However, they impose additional assumptions that restrict how past treatments affect future outcomes. We do not consider these setups in this paper.}, pooling all not-yet-treated units and using them as a comparison group does not respect the triple-differences identification assumptions and, as such, can lead to biased estimates for the parameters of interest. We also discuss straightforward and computationally simple estimators that bypass these problems. Throughout this section, we assume that all identification assumptions discussed in Section (ref) hold without covariates, i.e., by taking $X=1$ almost surely. Our methods naturally extend to setups where covariates are necessary for identification, as discussed in Section (ref). There, we also explore extensions to event-study aggregations and treatment effect heterogeneity across groups and time. To build intuition, we begin by noting that the way the DiD literature has addressed the shortcomings of using regression specifications akin to (ref) to infer overall average treatment effects is to decompose the problem into a series of $2$-period $2$-group ($2\times 2$) DiDs; for an overview, see Roth2023a and Baker_etal_2025_JEL. A popular strategy involves using the units not yet treated by period $t$ as a comparison group when estimating $ATT(g,t)$ (Callaway2021, Borusyaetal2024, Wooldridge2021). It is thus intuitive and natural to build on these DiD papers, Olden2022's DDD procedure, and the linear regression specification (ref), and attempt to estimate $ATT(g,t)$ in a DDD setup using \begin{small} \begin{align} \widehat{ATT}_{\text{cs-nyt}}(g,t) =& \Bigg[\bigg(\mathbb{E}_n\left[Y_{t} - Y_{g-1}|S = g, Q=1\right]\bigg)-\bigg( \mathbb{E}_n\left[ Y_{t} - Y_{g-1} | S = g, Q=0 \right] \bigg)\Bigg] \nonumber\\ &- \Bigg[\bigg(\mathbb{E}_n\left[ Y_{t} - Y_{g-1}| S >t, Q=1 \right] \bigg)\Bigg.-\Bigg.\bigg(\mathbb{E}_n\left[ Y_{t} - Y_{g-1}|S>t, Q=0 \right]\bigg)\Bigg] \end{align} \end{small} in any post-treatment periods $t\ge g$.\footnote{Since there are no covariates, we do not need to use three DiDs as we discussed in Section (ref). We use the notation $\text{cs}$ in (ref) to denote the estimator discussed above, which pinpoints the baseline period at period $g-1$.} The question now is whether (ref) indeed recovers $ATT(g,t)$'s under our identification assumptions. To answer this practically relevant question, we conduct some Monte Carlo simulations for a setup with three time periods, $t=1,2,3$, three treatment-enabling groups, $S\in\{2,3,\infty\}$, and two eligibility groups $Q=1$ and $Q=0$. We focus on $ATT(2,2)$, i.e., the average treatment effect in period two of being treated in period two, among units treated in period two. The true $ATT(2,2)$ in our simulations is 10. We considered a setup with $n=5,000$ and conducted $1,000$ simulation draws. To ease the exposition, we abstract from further details about the DGP and refer the reader to Section (ref) and Supplemental Appendix (ref) for a more detailed discussion. Panel (a) from Figure (ref) displays the density of the DDD estimates for $ATT(2,2)$ based on the estimator in (ref). This result makes it clear that, in general, (ref) is not a valid estimator for the $ATT(2,2)$ in DDD setups, as it is systematically biased. In fact, in our simulations, (ref) always leads to a negative estimate while the true effect is positive. This bias arises because the DDD parallel trends assumption is more flexible than its DiD counterpart: it allows for treatment-enabling-groups- and partition-specific violations of DiD-type parallel trends. In particular, when the fraction of eligible units differs across treatment-enabling groups $S$, pooling not-yet-treated units may conflate trends across heterogeneous populations, violating the assumptions necessary to interpret differences as causal. \begin{figure}[!htp!] \begin{center} \begin{subfigure}[t]{0.48\textwidth} \caption{DDD using pooled not-yet-treated units } \end{subfigure} \begin{subfigure}[t]{0.48\textwidth} \begin{center} \caption{DDD GMM using all not-yet-treated units} \end{center} \end{subfigure} \begin{subfigure}[t]{0.48\textwidth} \begin{center} \caption{DDD using never-treated units} \end{center} \end{subfigure} \begin{subfigure}[t]{0.48\textwidth} \begin{center} \caption{DDD using never-treated and GMM with not-yet-treated units} \end{center} \end{subfigure} \caption{Density of different staggered DDD estimates for ATT(2,2), without covariates} \end{center} \justifying \scriptsize{Notes: Simulation designs based on the design described in Section (ref) and Supplemental Appendix (ref), with $n=5,000$ and $1,000$ Monte Carlo repetitions. The true $ATT(2,2)$ is ten and is indicated in the solid vertical line in all panels. Panel (a) displays the density of DDD estimates that use the pooled not-yet-treated units as a comparison group as described in (ref). Panel (b) displays the density of the estimates based on our proposed DDD GMM estimator that uses all not-yet-treated units as a comparison group described in (ref). Panel (c) displays the density of the estimates based on our proposed DDD estimator that uses the never-treated units as a comparison group described in (ref) with $g_{\text{c}}=\infty$. Panel (d) compares DDD estimates using never-treated units (yellow curve) with GMM-based DDD using not-yet-treated (green curve), on the same scale. All densities are computed across all simulation draws. Panels(a)-(c) have the same x-axis range but different y-axis. } \end{figure} The key insight to address these problems is that we should be cautious when selecting the comparison group to estimate each $ATT(g,t)$. Such comparison groups must satisfy the DDD identification assumptions, which need to be verified on a group-by-group basis. Upon close inspection of Assumption (ref), a natural solution is to avoid pooling across treatment-enabling groups and use one of them at a time. Doing so yields multiple valid comparisons for the same $(g,t)$ group, generating an over-identified model. More precisely, for each available not-yet-enabled group $g_{\text{c}}>t$, we can use the following estimator for $ATT(g,t)$, $t\ge g$: \begin{small} \begin{align} \widehat{ATT}_{g_{\text{c}}}(g,t) =& \Bigg[\bigg(\mathbb{E}_n\left[Y_{t} - Y_{g-1}|S = g, Q=1\right]\bigg)-\bigg( \mathbb{E}_n\left[ Y_{t} - Y_{g-1} | S = g, Q=0 \right] \bigg)\Bigg] \nonumber\\ &- \Bigg[\bigg(\mathbb{E}_n\left[ Y_{t} - Y_{g-1}| S = g_{\text{c}}, Q=1 \right] \bigg)\Bigg.-\Bigg.\bigg(\mathbb{E}_n\left[ Y_{t} - Y_{g-1}|S=g_{\text{c}}, Q=0 \right]\bigg)\Bigg]. \end{align} \end{small} Note that when $g_{\text{c}}=\infty$, (ref) uses the set of units that never enabled treatment $S=\infty$ as the comparison group. However, one is not restricted to this unique comparison group. In the context of our simulation, one can also use the units that enabled treatment in period three to learn about $ATT(2,2)$. In this sense, instead of choosing which comparison group to use, we propose combining all available options and forming a more precise DDD estimator for the $ATT(g,t)$s. More concretely, we propose using \begin{align} \widehat{ATT}_{\text{gmm}}(g,t) = \dfrac{\mathbf{1}' \widehat{\Omega}^{-1}}{\mathbf{1}' \widehat{\Omega}_{g,t}^{-1}\mathbf{1}}\widehat{ATT}_{\text{dr}}(g,t), \end{align} where $\widehat{ATT}_{\text{dr}}(g,t)$ is the $k_{g,t}\times 1$-dimensional vector of all possible (non-collinear) estimators for $ATT(g,t)$ that uses a valid comparison group $g_{\text{c}}>t$, $\widehat{\Omega}_{g,t}$ is a consistent estimator of their variance-covariance matrix, and $\mathbf{1}$ is a ($k_{g,t}\times 1$-dimensional) vector of ones. We show that $ \widehat{ATT}_{\text{gmm}}(g,t)$ has a GMM interpretation based on re-centered influence functions in Remark (ref). Panels (b) and (c) of Figure (ref) display the density of the DDD estimates for $ATT(2,2)$ based on our GMM-based DDD estimator (ref) and our DDD estimator that only uses the never-treated as comparison group $g_{\text{c}}=\infty$ in (ref), respectively. As it is easy to see, both estimators are correctly centered at the true $ATT(2,2)$. As Panels (a) - (c) in Figure (ref) have the same scale, it is challenging to compare our DDD estimates that combine all not-yet-treated units with our DDD estimates that use never-treated units as comparison groups. In Figure (ref)(d), we address this issue and display their densities based on the $1,000$ simulation draws. Overall, it is evident that utilizing all not-yet-treated units can yield substantial gains in precision. The results of our simulations indicate that confidence intervals based on $ \widehat{ATT}_{g_{\text{c}}=\infty}(g,t)$ are around 50% wider than those based on $\widehat{ATT}_{\text{gmm}}(g,t)$, underscoring the appeal of using our GMM-based DDD estimator in terms of power. Overall, the results in this section underscore the broader lesson that DDD designs with staggered adoption should not be treated as simple extensions of DiD methods. The interaction between timing, eligibility, and heterogeneity in group composition introduces complexities that necessitate more careful attention to the identification argument and the construction of the comparison group. \section{The econometrics of DDD designs} In this section, we discuss the econometrics of DDD designs following the framework discussed in Section (ref). We start by establishing nonparametric identification of the $ATT(g,t)$'s under the identification assumptions in Section (ref). We then discuss estimation and inference procedures for $ATT(g,t)$'s and their event-study functional $ES(e)$ as defined in (ref). Throughout this section, we focus on setups where covariates are important for identification, i.e., all the assumptions discussed in Section (ref) are only plausible after you condition on covariates. Results for unconditional DDD setups follow as special cases by taking all covariates $X=1$ for all units. We also focus on staggered treatment adoption DDD setups, as they nest DDD setups with a single treatment date. \subsection{Identification} In this section, we establish the nonparametric identification for the $ATT(g,t)$'s in all post-treatment periods $t\ge g$ under Assumptions (ref), (ref), (ref), and (ref). Furthermore, we demonstrate that one can utilize regression adjustment/outcome regression (RA), inverse probability weighting (IPW), or doubly robust estimands to recover the $ATT(g,t)$'s. We also demonstrate that one can potentially utilize different comparison groups, thereby opening the door to combining them for potential efficiency gains. Before formalizing our results, we need to introduce some additional notation. Let $m_{Y_t-Y_{t'}}^{S=g,Q=q}(X) \equiv \mathbb{E}\left[Y_t - Y_{t'}|S=g,Q=q ,X\right]$ denote the population regression function of changes in outcomes from period $t'$ to period $t$ given covariates $X$ among units that enabled treatment in period $g$ ($S=g$) that belongs to eligibility group $q$ ($Q=q$). Analogously, let $p^{S=g,Q=1}_{g', q'}(X) \equiv \mathbb{P}[S=g, Q=1 | X, (S=g, Q=1) \cup (S=g',Q=q')]$ denote the generalized propensity score. Note that $p^{S=g,Q=1}_{g', q'}(X)$ indicates the probability of a unit being observed in enabling group $S=g$ and being eligible for treatment ($Q=1)$, conditional on pre-treatment covariates $X$ and on either being in the $S=g$ group and being eligible for treatment, or being in the $S=g'$ group with eligibility to treatment $Q=q'$.\footnote{We use this notion of generalize propensity score as it allow us to focus on sequences of two-groups comparisons as in Lechner2002 and Callaway2021. One can understand these generalized propensity scores as $p^{S=g,Q=1}_{g’, q’}(X) = \frac{p^{S=g,Q=1}(X)}{p^{S=g,Q=1}(X) + p^{S=g’,Q=q’}(X)}$ with $p^{S=g,Q=q}(X) \equiv \mathbb{P}[S=g, Q=q | X]$. We favor $p^{S=g,Q=1}_{g’, q’}(X)$ as this is how we implement these when constructing our DDD estimators.} For any $g_{\text{c}}\in \mathcal{S}$ such that $g_{\text{c}}>\max\{g,t\}$, and any post-treatment period $t\ge g$, let the doubly robust DDD estimand for the $ATT(g,t)$ be given by \begin{small} \begin{align} {ATT}_{\text{dr}, g_{\text{c}}}(g,t) =& \mathbb{E}\left[ \left({w}^{S=g,Q=1}_{\text{trt}}(S,Q) - {w}^{S=g,Q=1}_{g,0}(S,Q,X) \right)\left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=g,Q=0}(X) \right)\right]\nonumber \\ & +\mathbb{E}\left[ \left({w}^{S=g,Q=1}_{\text{trt}}(S,Q) - {w}^{S=g,Q=1}_{g_{\text{c}},1}(S,Q,X) \right)\left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=g_{\text{c}},Q=1}(X) \right)\right] \\ &- \mathbb{E}\left[ \left({w}^{S=g,Q=1}_{\text{trt}}(S,Q) -{w}^{S=g,Q=1}_{g_{\text{c}},0}(S,Q,X)\right) \left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=g_{\text{c}},Q=0}(X)\right)\right] \nonumber, \end{align} \end{small} where the weights ${w}$ are given by {\fontsize{10.4pt}{10pt}\selectfont \begin{align} {w}^{S=g,Q=1}_{\text{trt}}(S, Q)\equiv\dfrac{1{\{S=g,Q=1\}} }{\mathbb{E}[1{\{S=g,Q=1\}}]}, \quad {w}^{S=g, Q=1}_{g',q'}(S, Q, X) \equiv \dfrac{\dfrac{1{\{S=g',Q=q'\}} \cdot {p}^{S=g,Q=1}_{g',q'}(X) }{1 - {p}^{S=g,Q=1}_{g',q'}(X)}} {\mathbb{E}\left[\dfrac{1{\{S=g',Q=q'\}} \cdot {p}^{S=g,Q=1}_{g',q'}(X) }{1 - {p}^{S=g,Q=q}_{g',q'}(X)}\right]}. \end{align} } Analogously, let the RA DDD estimand for the $ATT(g,t)$ be given by \begin{small} \begin{equation} \resizebox{.93\textwidth}{!}{$ \begin{aligned} {ATT}_{\text{ra}, g_{\text{c}}}(g,t) =\ \mathbb{E}\Big[{w}^{S=g,Q=1}_{\text{trt}}(S, Q)\Big(Y_t - Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=g,Q=0}(X) - {m}_{Y_t-Y_{g-1}}^{S=g_{\text{c}},Q=1}(X) + {m}_{Y_t-Y_{g-1}}^{S=g_{\text{c}},Q=0}(X) \Big)\Big] \end{aligned} $} \end{equation} \end{small} and the IPW estimand be \begin{small} \begin{align} {ATT}_{\text{ipw}, g_{\text{c}}}(g,t) =& \mathbb{E}\left[\left({w}^{S=g,Q=1}_{\text{trt}}(S,Q) - {w}^{S=g,Q=1}_{g,0}(S,Q,X) \right) \left( Y_t-Y_{g-1}\right)\right]\nonumber\\ & - \mathbb{E}\left[\left({w}^{S=g,Q=1}_{g_{\text{c}},1}(S,Q,X) - {w}^{S=g,Q=1}_{g_{\text{c}},0}(S,Q,X)\right) \left( Y_t-Y_{g-1}\right)\right]. \end{align} \end{small} \begin{theorem} Let Assumptions (ref), (ref), (ref), and (ref) hold. Then, for all $g\in \mathcal{G}_{\text{trt}}$, $t\in \{2,\dots, T\}$, and $g_{\text{c}} \in \mathcal{S}$ such that $t\ge g$ and $g_{\text{c}}>t$, \begin{equation} ATT(g,t) = ATT_{\text{dr}, g_{\text{c}}}(g,t) = {ATT}_{\text{ra}, g_{\text{c}}}(g,t)={ATT}_{\text{ipw}, g_{\text{c}}}(g,t). \end{equation} \end{theorem} Theorem (ref) is the first main result of this paper. It establishes the nonparametric identification of all post-treatment $ATT(g,t)$'s in DDD setups. It extends the DiD identification results of Callaway2021 to DDD setups. As such, it also extends the difference-in-differences identification results based on the RA approach of Heckman1997, the IPW approach of Abadie2005, and the DR approach of SantAnna2020 to DDD setups with multiple periods and variation in treatment time. Theorem (ref) also highlights that one can use different parts of the data-generating process to identify the $ATT(g,t)$'s: the RA estimand only models the conditional expectation of evolution of outcomes among untreated units, the IPW approach only models the conditional probability of being observed in a given partition of the $S$-by-$Q$ groups, whereas the DR approach exploits both components. A big advantage of the DR approach is that it is based on a Neyman-orthogonal moment condition Belloni2017, and, therefore, it is more robust against model misspecifications than the IPW and RA formulations. It is very easy to show that estimators based on ${ATT}_{\text{dr}, g_{\text{c}}}(g,t)$ enjoy a very attractive doubly-robust property SantAnna2020 that allows for some forms of (global) model misspecifications.\footnote{For an overview of doubly robust estimators in cross-sectional designs, see section 2 of Sloczynski2018, and Seaman2018.} Another important result from Theorem (ref) is that our DDD model is over-identified, as we can use multiple not-yet-treated enabling groups $g_{\text{c}}$ as valid comparison groups. For instance, in a setup with $S\in \{2,3,\infty\}$, we can set $g_{\text{c}}=3$ or $g_{\text{c}}=\infty$ to identify $ATT(2,2)$, and both will lead to the same target parameter. As a direct consequence of this result, any weighted sum of these estimands that use different $g_{\text{c}}$'s will also lead to the $ATT(g,t)$, as long as the weights sum up to one. We formalize this result in the following corollary, using the DR estimand; however, this also applies to the RA and IPW. Let $\mathcal{G}_{\text{c}}^{\text{g,t}}=\{ g_{\text{c}} \in \mathcal{S}: g_{\text{c}} > \max\{g,t\}\}$. \begin{corollary} Let Assumptions (ref), (ref), (ref), and (ref) hold. Then, for all $g\in \mathcal{G}_{\text{trt}}$ and $t\in \{2,\dots, T\}$ such that $t\ge g$, and any set of weights $w^{\text{g,t}}_{\mathbf{g}_{\text{c}}}$ that sum up to one over $\mathcal{G}_{\text{c}}^{\text{g,t}}$, \begin{equation*} ATT(g,t) = \sum_{g_{\text{c}}\in \mathcal{G}_{\text{c}}^{\text{g,t}}} w^{\text{g,t}}_{g_{\text{c}}} ATT_{\text{dr}, g_{\text{c}}}(g,t). \end{equation*} \end{corollary} As Corollary (ref) indicates that all weighted sums lead to the same $ATT(g,t)$, a natural way to choose these weights is to pick them such that we maximize precision in terms of minimizing the resulting asymptotic variance. In the next session, we will discuss this in greater detail, connecting these arguments to a formulation based on generalized methods of moments using re-centered influence functions. \begin{remark} As we discussed in Section (ref), in two-period DDD setups without covariates, one can identify $ATT(2,2)$ using the difference of two DiD estimands as in (ref) Olden2022. This equivalence breaks down when the DDD identification assumptions are only satisfied after you condition on covariates $X$---see Figure (ref). The econometric reason for this failure of equivalence is that one needs to integrate the covariates using the covariate distribution among treated units, i.e., units with $S=2$ and $Q=1$. Proceeding as if $ATT(2,2)$ were the difference of two DiD estimands would integrate $X$ using the covariate distribution of untreated units ($S=\infty$ and $Q=1$), leading to biases. The results in Theorem (ref) address this problem by guaranteeing that one integrates out covariates using the correct reference distribution, which leads to a combination of three DiD estimands, not just two. \end{remark} \begin{remark} Although Theorem (ref) and Corollary (ref) allow one to use several different not-yet-treated cohorts $g_{\text{c}}$ as the comparison group, it does not allow one to pool all not-yet-treated units and use that pooled set of units as the aggregate comparison group to identify $ATT(g,t)$ in DDD. This sharply contrasts DiD procedures such as those discussed in Callaway2021---see Figure (ref) for an illustration of the bias that can arise by following this type of procedure. The econometric reasoning for such results is that Assumption (ref) allows for both enabling-group- and eligibility-group-specific trends, and it does not impose that the proportion of units in each eligibility group $Q$ is the same across all enabling groups $S$. As such, Assumption (ref) does not guarantee that, with probability one, \begin{small} \begin{eqnarray*} \mathbb{E}\left[Y_{t}(\infty) - Y_{t-1}(\infty)| S = g, Q=1, X\right] &-& \mathbb{E}\left[Y_{t}(\infty)- Y_{t-1}(\infty) | S = g, Q=0, X \right] \\[-0.2cm] &=&\\[-0.2cm] \mathbb{E}\left[Y_{t}(\infty) - Y_{t-1}(\infty)|S >t, Q=1, X\right] &-& \mathbb{E}\left[Y_{t}(\infty)- Y_{t-1}(\infty) | S >t, Q=0, X \right], \end{eqnarray*} \end{small} as it would be required to use the pooled, not-yet-treated units as a comparison group. \end{remark} \begin{remark} As Theorem (ref) establishes nonparametric identification of the $ATT(g,t)$'s over all post-treatment periods and that $\mathop{}\!\mathbb{P}(G=g|G+e \in [1,T])$ is also nonparametrically identified, it follows that event-study parameters that aggregate across eligibility-groups, $ES(e)$ as defined in (ref), is also nonparametrically identified. For instance, it follows that for any event-time $e\ge 0$, \begin{align} ES(e) = \sum_{g\in\mathcal{G}_{\text{trt}}} \mathop\!\mathbb{P}(G=g|G+e \in [1,T]) ATT_{\text{dr}, g_{\text{c}}}(g,g+e). \end{align} One can also replace $ATT_{\text{dr}, g_{\text{c}}}(g,g+e)$ with their analogs in Corollary (ref) or with the RA or IPW estimands in Theorem (ref). One can also use Theorem (ref) to establish the identification of many other aggregate summary causal parameters discussed in Section 3 of Callaway2021. \end{remark} \begin{remark} One of the biggest appeals of DiD and DDD setups is the availability of pre-treatment periods that allow the assessment of the plausibility of PT assumptions, such as Assumption (ref). Under Assumption (ref), a very popular way to assess the plausibility of PT is to construct event-study plots based on $ES(e)$ as in (ref), consider both pre-treatment ($e<0$) and post-treatment ($e\ge 0$) event times, and check whether pre-treatment event-study coefficients are all close to zero. It is straightforward to adapt this strategy in our DDD context by fixing the statistical estimand---for example, the ${ATT}_{\text{dr}, g_{\text{c}}}(g,t)$ in (ref)---consider pre-treatment periods $t<g$, and then aggregate them using cohort-size. More specifically, for any event-time $e < 0$, \begin{align} ES(e) = \sum_{g\in\mathcal{G}_{\text{trt}}} \mathop\!\mathbb{P}(G=g|G+e \in [1,T]) ATT_{\text{dr}, g_{\text{c}}}(g,g+e). \end{align} Note that when $e=-1$, $ES(e)=0$ by construction, as we fix the baseline period at the last untreated period for group $g$, $g-1$. Based on these event-study aggregations, it is also possible to conduct sensitivity analysis for the plausibility of Assumption (ref) using the results in rambachan_roth. \end{remark} \begin{remark} The DR DDD estimand can also be understood as a DDD estimand that builds on an efficient influence function in DDD setups with two periods. See Lemma (ref) in the appendix for such results. \end{remark} \subsection{Estimation and inference} In this section, we now propose simple-to-use plug-in estimators for the $ATT(g,t)$s and $ES(e)$s parameters, and discuss how one can conduct valid inference for these parameters. We focus on the doubly robust DDD estimator; the results for the RA and IPW DDD estimators are analogous. First, notice that for any $g_{\text{c}}\in \mathcal{S}$ such that $g_{\text{c}}>\max\{g,t\}$, Theorem (ref) suggests that we can estimate $ATT(g,t)$ by using the sample analogue of (ref), \begin{small} \begin{align} \widehat{{ATT}}_{\text{dr}, g_{\text{c}}}(g,t) =& \mathbb{E}_n\left[ \left(\widehat{w}^{S=g,Q=1}_{\text{trt}}(S,Q) - \widehat{w}^{S=g,Q=1}_{g,0}(S,Q,X) \right)\left(Y_t- Y_{g-1} - \widehat{m}_{Y_t-Y_{g-1}}^{S=g,Q=0}(X) \right)\right]\nonumber \\ & +\mathbb{E}_n\left[ \left(\widehat{w}^{S=g,Q=1}_{\text{trt}}(S,Q) - \widehat{w}^{S=g,Q=1}_{g_{\text{c}},1}(S,Q,X) \right)\left(Y_t- Y_{g-1} - \widehat{m}_{Y_t-Y_{g-1}}^{S=g_{\text{c}},Q=1}(X) \right)\right] \\ &- \mathbb{E}_n\left[ \left(\widehat{w}^{S=g,Q=1}_{\text{trt}}(S,Q) -\widehat{w}^{S=g,Q=1}_{g_{\text{c}},0}(S,Q,X)\right) \left(Y_t- Y_{g-1} - \widehat{m}_{Y_t-Y_{g-1}}^{S=g_{\text{c}},Q=0}(X)\right)\right] \nonumber, \end{align} \end{small} where the estimated weights $\widehat{w}$ are given by \begin{small} \begin{align*} \widehat{w}^{S=g,Q=1}_{\text{trt}}(S, Q)\equiv\dfrac{1{\{S=g,Q=1\}} }{\mathbb{E}_n[1{\{S=g,Q=1\}}]}, \quad \widehat{w}^{S=g, Q=1}_{g',q'}(S, Q, X) \equiv \dfrac{\dfrac{1{\{S=g',Q=q'\}} \cdot \widehat{p}^{S=g,Q=1}_{g',q'}(X) }{1 - \widehat{p}^{S=g,Q=1}_{g',q'}(X)}} {\mathbb{E}_n\left[\dfrac{1{\{S=g',Q=q'\}} \cdot \widehat{p}^{S=g,Q=1}_{g',q'}(X) }{1 - \widehat{p}^{S=g,Q=q}_{g',q'}(X)}\right]}, \end{align*} \end{small} and $\widehat{m}_{Y_t-Y_{g-1}}^{S=g_{\text{c}}, Q=1}(X)$ and $\widehat{p}^{S=g, Q=1}_{g',q'}(X)$ are (potentially misspecified) working models for the outcome regression ${m}_{Y_t-Y_{g-1}}^{S=g_{\text{c}}, Q=1}(X)$ and the generalized propensity score $\widehat{p}^{S=g, Q=1}_{g',q'}(X)$. These estimators extend the DR DiD estimator of Callaway2021 to the DDD setup, and remain consistent if \emph{either} outcome regression or generalized propensity score models are correctly specified. It is also worth stressing that we do not need that \emph{all} generalized propensity score working models or \emph{all} outcome regression working models in (ref) to be correctly specified to get a consistent DDD estimator for $ATT(g,t)$; it suffices that any of the working models within each of the 3 DR DiD components of (ref) to be correctly specified, allowing a greater deal of estimation flexibility.\footnote{Some people may call this a multiply-robust estimator, as one has more than two opportunities to estimate the target parameter consistently. For simplicity, we retain the doubly robust terminology to avoid new acronyms. } As Corollary (ref) highlights, one can also combine several $\widehat{{ATT}}_{\text{dr}, g_{\text{c}}}(g,t)$ that leverage different comparison groups $g_{\text{c}}$, i.e., for any (consistently estimated) weights $\widehat{w}^{\text{g,t}}_{g_{\text{c}}}$ that sum up to one over $ \mathcal{G}_{\text{c}}$, \begin{equation} \widehat{ATT}_{\text{dr},\widehat{w}}(g,t) = \sum_{g_{\text{c}}\in \mathcal{G}_{\text{c}}^{\text{g,t}}} \widehat{w}^{\text{g,t}}_{g_{\text{c}}} \widehat{ATT}_{\text{dr}, g_{\text{c}}}(g,t) = \widehat{w}^{\text{g,t} '}\widehat{ATT}_{\text{dr}}(g,t), \end{equation} where $\widehat{ATT}_{\text{dr}}(g,t)$ is the $k_{g,t}\times 1$ vector of $\widehat{{ATT}}_{\text{dr}, g_{\text{c}}}(g,t)$ for all $g_{\text{c}} \in \mathcal{G}_{\text{c}}^{\text{g,t}}$, and $\widehat{w}^{\text{g,t}~'}$ is a $k_{g,t}\times 1$ vector of (estimated) weights that sum up to one, i.e., for a generic vector of ones $\textbf{1}$, $\textbf{1}'w^{\text{g,t}}=1$. A natural question that arises is: how should one choose these weights $\widehat{w}^{\text{g,t}}_{g_{\text{c}}}$? We propose to choose the weights that lead to the asymptotically most precise (minimum variance) estimator for $ATT(g,t)$, that is, to pick weights that solve \begin{align} \min_{w^{\text{g,t}}} w^{\text{g,t} '} \widehat{\Omega}_{g,t} w^{\text{g,t}} \text{ subject to } \textbf{1}'w^{\text{g,t}} = 1, \end{align} where $\widehat{\Omega}_{g,t}$ is a $k_{g,t}\times k_{g,t}$ consistent estimator for the variance-covariance matrix of $\widehat{ATT}_{\text{dr}}(g,t)$. Notice that the solution of (ref) admits a closed-form solution, and the optimal weights are given by \begin{align} \widehat{w}^{\text{g,t}}_{\text{gmm}} =\dfrac{\widehat{\Omega}^{-1}_{g,t}\textbf{1}}{\textbf{1}'\widehat{\Omega}^{-1}_{g,t}\textbf{1}} . \end{align} In turn, this implies that the linear combination of $\widehat{{ATT}}_{\text{dr}, g_{\text{c}}}(g,t)$ that leads to the most precise estimator for $ATT(g,t)$ is given by \begin{equation} \widehat{ATT}_{\text{dr},\text{gmm}}(g,t) = \dfrac{\mathbf{1}' \widehat{\Omega}_{g,t}^{-1}}{\mathbf{1}' \widehat{\Omega}_{g,t}^{-1}\mathbf{1}} \widehat{ATT}_{\text{dr}}(g,t). \end{equation} In many situations with multiple periods and variation in treatment time, researchers are interested in summarizing the $ATT(g,t)$'s into fewer parameters that highlight treatment effect heterogeneity with respect to the time elapsed since treatment take-up. That is, very often, researchers are interested in estimating event-study type parameters $ES(e)$ as defined in (ref). A very natural estimator for $ES(e)$ is the plug-in estimator, where we replace $ATT(g,t)$ with $\widehat{ATT}_{\text{dr},\text{gmm}}(g,t)$ (or $\widehat{ATT}_{\text{dr},g_{\text{c}}}(g,t)$), and $\mathop{}\!\mathbb{P}(G=g|G+e \in [1,T])$ by its sample analogue, that is, \begin{equation} \widehat{ES}_{\text{dr},\text{gmm}}(e) = \sum_{g\in\mathcal{G}_{\text{trt}}} \mathop\!\mathbb{P}_n(G=g|G+e \in [1,T]) \widehat{ATT}_{\text{dr},\text{gmm}}(g,g+e), \end{equation} where $\mathop{}\!\mathbb{P}_n(G=g|G+e \in [1,T]) = \sum_{i=1}^n 1\{G_i = g\} 1\{G_i + e \in [1,T]\}\big/ \sum_{j=1}^n 1\{G_j + e \in [1,T]\}$. We can define $ \widehat{ES}_{\text{dr},g_{\text{c}}}(e)$ analogously by replacing $\widehat{ATT}_{\text{dr},\text{gmm}}(g,g+e)$ with $\widehat{ATT}_{\text{dr},g_{\text{c}}}(g,g+e)$ on (ref). Based on it, we can also estimate an overall summary parameter by averaging all post-treatment event times, i.e., \begin{align} \widehat{ES}_{\text{avg},\text{gmm}} = \dfrac{1}{N_E}\sum_{e\in \mathcal{E}} \widehat{ES}_{\text{dr},\text{gmm}}(e). \end{align} \begin{remark} It is also worth noticing that $\widehat{ATT}_{\text{dr},\text{gmm}}(g,t)$ in (ref) can be interpreted as an optimal Generalized Method of Moments (GMM) estimator based on re-centered influence functions. To see this, let ${\mathbb{IF}}_{\text{dr}, g_{\text{c}}}(g,t)$ denote the influence function of $\sqrt{n}\left(\widehat{ATT}_{\text{dr}, g_{\text{c}}}(g,t) - ATT_{\text{dr}, g_{\text{c}}}(g,t)\right)$. Let ${\mathbb{RIF}}_{\text{dr}, g_{\text{c}}}(g,t) = {\mathbb{IF}}_{\text{dr}, g_{\text{c}}}(g,t) + {ATT}_{\text{dr},g_{\text{c}}}(g,t)$ denote its re-centered influence function, and denote the $k_{g,t}\times 1$ vector of all ${\mathbb{RIF}}_{\text{dr}, g_{\text{c}}}(g,t)$ for $g_{\text{c}} \in \mathcal{G}_{\text{c}}^{\text{g,t}}$ by ${\mathbb{RIF}}_{\text{dr}}(g,t)$. Since influence functions are mean zero, and that $ATT(g,t) = ATT_{\text{dr}, g_{\text{c}}}(g,t)$ for any $g_{\text{c}} \in \mathcal{G}_{\text{c}}^{\text{g,t}}$, we have the vector of moment conditions $\mathop{}\!\mathbb{E}[{\mathbb{RIF}}_{\text{dr}}(g,t) - \theta^{g,t}] = 0$, with $\theta^{g,t} = ATT(g,t)$. From standard GMM results Newey_McFadden_1994_Handbook, it follows that, under mild regularity conditions, the optimal (population) GMM estimator for $\theta^{g,t}$ is given by $$\theta^{g,t}_{\text{gmm}} = \dfrac{\mathbf{1}' {\Omega}_{g,t}^{-1}}{\mathbf{1}' {\Omega}_{g,t}^{-1}\mathbf{1}} \mathop{}\!\mathbb{E}[ {\mathbb{RIF}}_{\text{dr}}(g,t)] = \dfrac{\mathbf{1}' {\Omega}_{g,t}^{-1}}{\mathbf{1}' {\Omega}_{g,t}^{-1}\mathbf{1}} ATT_{\text{dr}}(g,t),$$ where the last equality follows from $\mathop{}\!\mathbb{E}[ {\mathbb{IF}}_{\text{dr}}(g,t)] = 0$. Thus, $\widehat{ATT}_{\text{dr},\text{gmm}}(g,t)$ in (ref) is the sample-analogy of the efficient population GMM $\theta^{g,t}_{\text{gmm}}$. \end{remark} \subsubsection{Asymptotic theory for ATT(g,t)'s} In what follows, we derive the large sample properties of our DR DDD estimators $\widehat{ATT}_{\text{dr},g_{\text{c}}}(g,t)$ and $\widehat{ATT}_{\text{dr},\text{gmm}}(g,t)$. All our results are derived for the large $n$, fixed $T$ paradigm. For a generic $Z$, let $|| Z || = \sqrt{trace(Z'Z)}$ denote the Euclidean norm of $Z$ and set $W_i=(Y_{i,t=1},\dots,Y_{i,t=T}, X_i', G_i, S_i, Q_i)'$; we will omit the index $i$ to unclutter the notation. Let $g(\cdot)$ be a generic notation for the outcome regressions $m_{Y_t-Y_{t'}}^{S=g',Q=q}(X)$ and generalized propensity scores $p^{S=g,Q=1}_{g', q'}(X)$, and, with some abuse of notation, let $g(\cdot; \gamma)$ denote a parametric model for $g(\cdot)$ that is known up to the finite-dimensional parameters $\gamma$. For a generic $\kappa^{g,t}_{g_{\text{c}}} = (\gamma^{ps~\prime}_{g,t,g_{\text{c}}}, \gamma^{reg~\prime}_{ g,t,g_{\text{c}}})'$, with $\gamma^{ps}_{g,t,g_{\text{c}}}$ and $\gamma^{reg}_{ g,t,g_{\text{c}}}$ being nuisance parameters for the generalized propensity score and outcome regressions, respectively, let \begin{small} \begin{align*} h^{g,t}_{g_{\text{c}}}(W; \kappa^{g,t}_{g_{\text{c}}}) =& \left({w}^{S=g,Q=1}_{\text{trt}}(W) - {w}^{S=g,Q=1}_{g,0}(W;\gamma^{ps}_{g,t,g_{\text{c}}}) \right)\left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=g,Q=0}(X;\gamma^{reg}_{ g,t,g_{\text{c}}}) \right)\nonumber \\ & +\left({w}^{S=g,Q=1}_{\text{trt}}(W) - {w}^{S=g,Q=1}_{g_{\text{c}},1}(W;\gamma^{ps}_{g,t,g_{\text{c}}}) \right)\left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=g_{\text{c}},Q=1}(X;\gamma^{reg}_{ g,t,g_{\text{c}}}) \right) \\ &- \left({w}^{S=g,Q=1}_{\text{trt}}(W) -{w}^{S=g,Q=1}_{g_{\text{c}},0}(W;\gamma^{ps}_{g,t,g_{\text{c}}})\right) \left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=g_{\text{c}},Q=0}(X;\gamma^{reg}_{ g,t,g_{\text{c}}})\right) \nonumber, \end{align*} \end{small} where the weights $w(W;\gamma^{ps}_{g,t,g_{\text{c}}})$ are defined similarly to those in (ref), with the difference being that the true unknown generalized propensity score models are replaced by working parametric counterparts, $p^{S=g,Q=1}_{g', q'}(X;\gamma^{ps}_{g,t,g_{\text{c}}})$, and the true unknown outcome regression models ${m}_{Y_t-Y_{g-1}}^{S=g',Q=q}(X)$ are also replaced with parametric working models, ${m}_{Y_t-Y_{g-1}}^{S=g',Q=q}(X;\gamma^{reg}_{ g,t,g_{\text{c}}})$. We denote the vector of pseudo-true parameters by $\kappa^{g,t}_{0, g_{\text{c}}}$ and let $\dot{h}^{g,t}_{g_{\text{c}}}(\kappa) = \partial h^{g,t}_{g_{\text{c}}}(W;\kappa) / \partial \kappa$. To derive our results, we make the following relatively mild assumptions. \begin{namedassumption}{WM}[Working Model Conditions] (i) $g(x;\gamma)$ is a parametric model for $g(x)$, where $ \gamma \in \Theta \subset \mathbb{R}^{d_k}$ is a compact set; (ii) the mapping $\theta \mapsto g(X ; \theta)$ is a.s. continuous; (iii) the pseudo-true parameter $\theta_0 \in \operatorname{int}(\Theta)$ satisfies that for an appropriate criterion function $Q: \Theta \rightarrow \mathbb{R}$ and for any $\epsilon > 0$, there exists some $\delta >0 $ such that $\inf _{\theta \in \Theta:\left\|\theta-\theta_0\right\| \geq \epsilon} Q(\theta)-Q\left(\theta_0\right)>\delta$; (iv) there exists an open neighborhood $\Theta_0 \subset \Theta$ containing $\theta_0$ such that $g(X;\gamma)$ is a.s. continuously differentiable in a neighborhood of $\gamma_0 \in \Theta_0$. In addition, (v) there exists some $\epsilon > 0$ such that, for all $(g,g',q') \in \mathcal{G}_{\text{trt}}\times \mathcal{G}^{g,t}_c \times \{0,1\}$, we have that $0 \leq p^{S=g,Q=1}_{g',q'}(X;\theta) \leq 1 - \epsilon$ a.s. for all $\theta \in \operatorname{int}(\Theta_{ps})$, where $\Theta_{ps}$ denotes the parameter space of $\gamma$ for the generalized propensity score working model. \end{namedassumption} \begin{namedassumption}{ALR}[$\sqrt{n}$-Asymptotically Linear Representation] Let $\hat{\theta}$ be a strongly consistent estimator of $\theta_0 \mapsto g(x ; \theta_0)$ and satisfy the following linear expansion \begin{align} \sqrt{n}\left(\widehat{\theta}-\theta_0\right)=\frac{1}{\sqrt{n}} \sum_{i=1}^n l \left(W_i ; \theta_0\right)+o_p(1) \end{align} where $l \left(\cdot ; \cdot\right)$ is a function such that $\mathop{}\!\mathbb{E}\left[l \left(W_i ; \theta_0\right)\right] = 0$; $\mathop{}\!\mathbb{E}\left[l \left(W_i ; \theta_0\right) \cdot l \left(W_i ; \theta_0\right)^{'}\right] < \infty$ and is positive definite; and $\lim_{\delta \to 0}\mathop{}\!\mathbb{E}\left[\sup_{\theta \in \Theta_0: \left\|\theta-\theta_0\right\| \leq \delta} \left\| l(W; \theta) - l(W; \theta_0) \right\|^2 \right] = 0$. \end{namedassumption} \begin{namedassumption}{IC}[Integrability Conditions] For each $g \in \mathcal{G}_{\text{trt}}$, $t \in \{2,\dots, T\}$, and $g'\in \mathcal{G}^{g,t}_c$, assume that $\mathop{}\!\mathbb{E}[\| h^{g,t}_{g_{\text{c}}}(W; \kappa^{g,t}_{0, g_{\text{c}}}) \|^{2}] < \infty$ and $\mathop{}\!\mathbb{E}\left[\sup_{\kappa \in \Gamma_0} \left| \dot{h}^{g,t}_{g_{\text{c}}}(\kappa) \right|\right] < \infty$, where $\Gamma_0$ is a small neighborhood of the pseudo-true parameter $\kappa^{g,t}_{0,g_{\text{c}}}$. \end{namedassumption} Assumptions (ref), (ref), and (ref) are standard in the literature; see e.g., Abadie2005, WOOLDRIDGE20071281, SantAnna2020, Callaway2021. Assumptions (ref) and (ref) impose a well-behaved parametric structure for the first-step estimators for the nuisance parameters. This assumption is made for statistical convenience and acknowledges that, in many DDD applications, the number of units in each group is small, making it difficult to adopt a nonparametric approach reliably. It is relatively straightforward to relax these conditions and allow for nonparametric or data-adaptive/machine-learning-based estimators; see, e.g., Ahrens2025_DML_JEL for an empirically-oriented discussion of causal double machine learning methods. Assumption (ref) imposes mild regularity constraints on the moments of the estimating equations, preventing ill-behaved variance properties and ensuring the stability of higher-order approximations. In what follows, we omit $W$ and $X$ from the weights and outcome regressions to minimize notation, and for a generic $\kappa^{g,t}_{g_{\text{c}}}$, let \begin{align} \psi^{g,t}_{g_{\text{c}}}(W; \kappa^{g,t}_{g_{\text{c}}}) = \psi^{g,t}_{S=g,Q=0}(W; \kappa^{g,t}_{g_{\text{c}}}) + \psi^{g,t}_{S=g_{\text{c}},Q=1}(W; \kappa^{g,t}_{g_{\text{c}}}) - \psi^{g,t}_{S=g_{\text{c}},Q=0}(W; \kappa^{g,t}_{g_{\text{c}}}), \end{align} where, for $(g',q') \in \{(g,0), (g_{\text{c}},1), (g_{\text{c}},0)\}$, $\psi^{g,t}_{S=g',Q=q'}(W; \kappa^{g,t}_{g_{\text{c}}})$ is an influence function for one of the three DR DiD components of the DR DDD, and is given by \begin{align} \psi^{g,t}_{S=g',Q=q'}(W; \kappa^{g,t}_{g_{\text{c}}}) = \psi^{g,t,1}_{S=g',Q=q'}(W; \kappa^{g,t}_{g_{\text{c}}}) - \psi^{g,t,0}_{S=g',Q=q'}(W; \kappa^{g,t}_{g_{\text{c}}}) - \psi^{g,t,est}_{S=g',Q=q'}(W; \kappa^{g,t}_{g_{\text{c}}}), \end{align} with \begin{align*} \psi^{g,t,1}_{S=g',Q=q'}(W; \kappa^{g,t}_{g_{\text{c}}})= & {w}^{S=g,Q=1}_{\text{trt}}\left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=a,Q=b}(\gamma^{reg}_{ g,t,g_{\text{c}}}) \right) \\ & - {w}^{S=g,Q=1}_{\text{trt}}\mathop\!\mathbb{E}\left[{w}^{S=g,Q=1}_{\text{trt}}\left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=a,Q=b}(\gamma^{reg}_{ g,t,g_{\text{c}}}) \right) \right]\\ \psi^{g,t,0}_{S=g',Q=q'}(W; \kappa^{g,t}_{g_{\text{c}}})= & {w}^{S=g,Q=1}_{g',q'}(\gamma^{ps}_{ g,t,g_{\text{c}}})\left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=a,Q=b}(\gamma^{reg}_{ g,t,g_{\text{c}}}) \right) \\ & - {w}^{S=g,Q=1}_{g',q'}(\gamma^{ps}_{ g,t,g_{\text{c}}})\mathop\!\mathbb{E}\left[{w}^{S=g,Q=1}_{g',q'}(\gamma^{ps}_{ g,t,g_{\text{c}}})\left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=a,Q=b}(\gamma^{reg}_{ g,t,g_{\text{c}}}) \right) \right] \end{align*} and \begin{align*} \psi^{g,t,est}_{S=g',Q=q'}(W; \kappa^{g,t}_{g_{\text{c}}}) = l^{g,t,reg}_{S=g',Q=q'}(\gamma^{ref}_{ g,t,g_{\text{c}}})' M^{g,t, 1}_{S=g',Q=q'}(\kappa^{g,t}_{g_{\text{c}}}) + l^{g,t,ps}_{S=g',Q=q'}(\gamma^{ps}_{ g,t,g_{\text{c}}})' M^{g,t,2}_{S=g',Q=q'}(\kappa^{g,t}_{g_{\text{c}}}) \end{align*} where $l^{g,t,reg}_{S=g',Q=q'}(\cdot)$ is the asymptotic linear representation of the outcome evolution for the group with $S=g'$ and $Q=q'$ as described in Assumption (ref), $l^{g,t,ps}_{S=g',Q=q'}(\cdot)$ is defined analogously for the generalized propensity score that uses group $S=a,Q=b$ as a comparison group, and {\fontsize{10.4pt}{10pt}\selectfont \begin{align*} M^{g,t,1}_{S=g',Q=q'}(\kappa^{g,t}_{g_{\text{c}}}) &= \mathop\!\mathbb{E}\left[ \left({w}^{S=g,Q=1}_{\text{trt}} - {w}^{S=g,Q=1}_{g',q'}(\gamma^{ps}_{g,t,g_{\text{c}}})\right) \dot{m}_{Y_t-Y_{g-1}}^{S=a,Q=b}(\gamma^{reg}_{g,t,g_{\text{c}}}) \right],\\ M^{g,t,2}_{S=g',Q=q'}(\kappa^{g,t}_{g_{\text{c}}}) & = \mathop\!\mathbb{E} \left[ {\alpha}^{S=g, Q=1}_{g',q'}(\gamma^{ps}_{ g,t,g_{\text{c}}}) \left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=a,Q=b}(\gamma^{reg}_{ g,t,g_{\text{c}}}) \right) \cdot \dot{p}^{S=g,Q=q}_{g',q'}(\gamma^{ps}_{g,t,g_{\text{c}}})\right]\\ & - \mathop\!\mathbb{E} \left[ {\alpha}^{S=g, Q=1}_{g',q'}(\gamma^{ps}_{ g,t,g_{\text{c}}}) \left(\mathop\!\mathbb{E}\left[ {w}^{S=g,Q=1}_{g',q'}(\gamma^{ps}_{g,t,g_{\text{c}}})\left(Y_t- Y_{g-1} - {m}_{Y_t-Y_{g-1}}^{S=a,Q=b}(\gamma^{reg}_{ g,t,g_{\text{c}}}) \right)\right] \right) \cdot \dot{p}^{S=g,Q=q}_{g',q'}(\gamma^{ps}_{g,t,g_{\text{c}}})\right] \end{align*} } with $\dot{m}_{Y_t-Y_{g-1}}^{S=a,Q=b}(\gamma^{reg}_{g,t,g_{\text{c}}}) = \partial{m}_{Y_t-Y_{g-1}}^{S=a,Q=b}(\gamma^{reg}_{g,t,g_{\text{c}}})\big/ \partial \gamma^{reg}_{g,t,g_{\text{c}}}$, $\dot{p}^{S=g,Q=q}_{g',q'}(\gamma^{ps}_{g,t,g_{\text{c}}}) = \partial {p}^{S=g,Q=q}_{g',q'}(\gamma^{ps}_{g,t,g_{\text{c}}}) \big/\gamma^{ps}_{g,t,g_{\text{c}}} $, and {\fontsize{10.4pt}{10pt}\selectfont \begin{align*} {\alpha}^{S=g, Q=1}_{g',q'}(\gamma^{ps}_{ g,t,g_{\text{c}}}) = \left. \dfrac{1{\{S=a,Q=b\}} }{\left(1 - {p}^{S=g,Q=1}_{g',q'}(X;\gamma^{ps}_{g,t,g_{\text{c}}})\right)^2} \right/ \mathbb{E}\left[\dfrac{1{\{S=a,Q=b\}} \cdot {p}^{S=g,Q=1}_{g',q'}(X;\gamma^{ps}_{g,t,g_{\text{c}}}) }{1 - {p}^{S=g,Q=q}_{g',q'}(X;\gamma^{ps}_{g,t,g_{\text{c}}})}\right]. \end{align*} } For each $g\in \mathcal{G}_{\text{trt}}$ and each $t\in \{2,3,\dots, \}$, let ${ATT}_{\text{dr}}(g,t)$ denote the $k_{g,t}\times 1$ vector of ${{ATT}}_{\text{dr}, g_{\text{c}}}(g,t)$ for all (non-collinear) $g_{\text{c}} \in \mathcal{G}_{\text{c}}^{\text{g,t}}$, and ${\Omega}_{g,t}$ be the asymptotic variance-covariance matrix of $\sqrt{n} \left( \widehat{ATT}_{\text{dr}}(g,t) - {ATT}_{\text{dr}}(g,t)\right)$, i.e., $ {\Omega}_{g,t} = \mathop{}\!\mathbb{E}\left[ \psi^{g,t}(W; \kappa^{g,t})\psi^{g,t}(W; \kappa^{g,t})'\right]$, with $\psi^{g,t}(W; \kappa^{g,t})$ being the $k_{g,t} \times 1$ vector that stacks all non-collinear $\psi^{g,t}_{g_{\text{c}}}(W; \kappa^{g,t}_{g_{\text{c}}})$ for $g_{\text{c}} \in \mathcal{G}_{\text{c}}^{\text{g,t}}$. Let ${ATT}_{\text{dr},\text{gmm}}(g,t) = ({\textbf{1}'{\Omega}^{-1}_{g,t}\textbf{1}})^{-1} {\textbf{1}'{\Omega}^{-1}_{g,t}}~ ATT_{\text{dr}}(g,t)$, and for a generic set of weights that sum up to one, let $ {ATT}_{\text{dr},{w}}(g,t) = \sum_{g_{\text{c}}\in \mathcal{G}_{\text{c}}^{\text{g,t}}} w^{\text{g,t}}_{g_{\text{c}}}~ ATT_{\text{dr}, g_{\text{c}}}(g,t)$, and recall that $\widehat{ATT}_{\text{dr},\widehat{w}}(g,t)$ is its empirical analogue as defined in (ref). Finally, let $\widehat{\Omega}_{g,t}$ be the empirical analogue of ${\Omega}_{g,t}$, where one replaces expectations by sample analogues and $\kappa^{g,t}$ with $\widehat{\kappa}^{g,t}$, and consider the following claim: \begin{align} &\text{For each } g\in \mathcal{G}_{\text{trt}}, t \in \{2,\dots, T\} \text{ such that }t\geq g, \text{ and each } g_{\text{c}}\in \mathcal{G}_{\text{c}}^{\text{g,t}} ,\nonumber\\ &\text{ we have that, for each } (g',q') \in \{(g,0), (g_{\text{c}},1), (g_{\text{c}},0)\},\nonumber\\ & \exists \gamma^{ps}_{0, g,t,g_{\text{c}}} \in \Theta^{ps}: \mathop\!\mathbb{P}({p}^{S=g,Q=q}_{g',q'}(X;\gamma^{ps}_{0,g,t,g_{\text{c}}}) = {p}^{S=g,Q=q}_{g',q'}(X)) = 1 \text{ or} \\ & \exists \gamma^{reg}_{0, g,t,g_{\text{c}}} \in \Theta^{reg}: \mathop\!\mathbb{P}({m}^{S=g,Q=q}_{g',q'}(X;\gamma^{reg}_{0,g,t,g_{\text{c}}}) = {m}^{S=g,Q=q}_{g',q'}(X)) = 1. \nonumber \end{align} Claim (ref) states that for each $(g,t)$-pair and each suitable comparison group $g_{\text{c}}$, either the working parametric model for the generalized propensity score is correctly specified, or the working outcome regression model for the comparison group is correctly specified for each of the three DiD components of our DDD estimator. Thus, eight possible working model combinations would lead to consistent DDD estimation of the $ATT(g,t)$ parameter. The next theorem establishes the limiting distribution of $\widehat{ATT}_{\text{dr},g_{\text{c}}}(g,t)$ and $\widehat{ATT}_{\text{dr},\text{gmm}}(g,t)$. \begin{theorem}[Consistency and Asymptotic Normality] Let Assumptions (ref), (ref), (ref), (ref), (ref), (ref), and (ref) hold. Then, for all $g\in \mathcal{G}_{\text{trt}}$, $t\in \{2,\dots, T\}$, and $g_{\text{c}} \in \mathcal{G}_{\text{c}}^{\text{g,t}}$ such that $t\ge g$, provided that (ref) is true, \begin{align*} \sqrt{n}\left(\widehat{ATT}_{\text{dr},g_{\text{c}}}(g,t) - {ATT}(g,t) \right) &= \dfrac{1}{\sqrt{n}}\sum_{i=1}^n \psi^{g,t}_{g_{\text{c}}}(W_i; \kappa^{g,t}_{0, g_{\text{c}}})+ o_p(1) \overset{d}{\rightarrow} N(0,\Omega_{g,t,g_{\text{c}}}), \end{align*} where $\Omega_{g,t,g_{\text{c}}} = \mathop{}\!\mathbb{E}\left[ \psi^{g,t}_{g_{\text{c}}}(W_i; \kappa^{g,t}_{0, g_{\text{c}}})\psi^{g,t}_{g_{\text{c}}}(W_i; \kappa^{g,t}_{0, g_{\text{c}}})'\right]$. Furthermore, \begin{align*} \sqrt{n}\left(\widehat{ATT}_{\text{dr},\text{gmm}}(g,t) - {ATT}(g,t) \right) &= \dfrac{\textbf{1}'{\Omega}^{-1}_{g,t}}{\textbf{1}'{\Omega}^{-1}_{g,t}\textbf{1}} \dfrac{1}{\sqrt{n}}\sum_{i=1}^n \psi^{g,t}(W_i; \kappa^{g,t}_{0})+ o_p(1) \overset{d}{\rightarrow} N(0,\Omega_{g,t, \text{gmm}}), \end{align*} where $\Omega_{g,t, \text{gmm}} = \left({\textbf{1}'{\Omega}^{-1}_{g,t}\textbf{1}}\right)^{-1} \leq \Omega_{g,t,g_{\text{c}}}$ for any $g_{\text{c}} \in \mathcal{G}_{\text{c}}^{\text{g,t}}$. In fact, for any set of weights $w$ that sum up to one over the $\mathcal{G}_{\text{c}}^{\text{g,t}}$, $\Omega_{g,t, \text{gmm}} \leq \Omega_{g,t, w}$, with $\Omega_{g,t, w}$ defined as the asymptotic variance of $\sqrt{n}\left(\widehat{ATT}_{\text{dr},\widehat{w}}(g,t) - {ATT}_{\text{dr},{w}}(g,t) \right)$. \end{theorem} Theorem (ref) provides the influence function for estimating each $ATT(g,t)$, using different comparison groups $g_{\text{c}}$, as well as establishes the consistency and asymptotic normality of our DR DDD estimator $\widehat{ATT}_{\text{dr},g_{\text{c}}}(g,t)$. Theorem (ref) also highlights that combining different comparison groups as our DR DDD estimator $\widehat{ATT}_{\text{dr},\text{gmm}}(g,t)$ does is effective in terms of asymptotically improving precision. That is, Theorem (ref) highlights that $\widehat{ATT}_{\text{dr},\text{gmm}}(g,t)$ is optimal in the sense that it asymptotically achieves the minimum variance across all weighted average estimators that combine multiple $\widehat{ATT}_{\text{dr},g_{\text{c}}}(g,t)$s. Importantly, Theorem (ref) also highlights the doubly (or multiply) robust property of our DDD estimators: they recover the $ATT(g,t)$ provided that each of the three DR DiD estimators has a correctly specified outcome regression or generalized propensity score working model. \begin{remark} Although Theorem (ref) provides pointwise inference results for each $ATT(g,t)$, it is straightforward to extend it to hold simultaneously across multiple $ATT(g,t)$'s. For instance, by letting $\widehat{ATT}_{\text{gmm},t\ge g}$ and ${ATT}_{\text{gmm},t\ge g}$ denote the vector of $\widehat{ATT}_{\text{dr},\text{gmm}}(g,t)$ and ${ATT}_{\text{dr},\text{gmm}}(g,t)$, respectively, for all $g\in\mathcal{G}_{\text{trt}}$, $t\in \{2,\dots, T,\}$ such that $t\ge g$, it is straightforward to show that $\sqrt{n}\left(\widehat{ATT}_{\text{gmm},t\ge g} - {ATT}_{\text{gmm},t\ge g}\right) \overset{d}{\rightarrow} N(0,\Omega)$, with $\Omega = \mathop{}\!\mathbb{E}[\psi^{t\ge g}_{\text{gmm}}(W_i; \kappa^{t\ge g}_{0}) \psi^{t\ge g}_{\text{gmm}}(W_i; \kappa^{t\ge g}_{0})']$, with $\psi^{t\ge g}_{\text{gmm}}(W_i; \kappa^{t\ge g}_{0})$ the asymptotic linear representation of $\sqrt{n}\left(\widehat{ATT}_{\text{gmm},t\ge g} - {ATT}_{\text{gmm},t\ge g}\right)$. One can then construct simultaneous confidence bands using a simple-to-use multiplier bootstrap as discussed in Theorem 3 and Algorithm 1 of Callaway2021. It is also straightforward to conduct cluster-robust inference; see Remark 10 of Callaway2021. As these results are commonly accessible, we will not include them here to conserve space. \end{remark} \subsubsection{Asymptotic theory for event-study parameters} In this section, we derive large sample properties for our event-study estimator $\widehat{ES}_{\text{dr},\text{gmm}}(e)$ as defined in (ref). Given that $\mathop{}\!\mathbb{P}_n(G=g|G+e \in [1,T])$ is an $\sqrt{n}$-consistent and asymptotically normal estimator of $\mathop{}\!\mathbb{P}(G=g|G+e \in [1,T])$, then for all $g \in \mathcal{G}_{\text{trt}}$, we have that \begin{equation} \sqrt{n}(\mathop\!\mathbb{P}_n(G=g|G+e \in [1,T]) - \mathop\!\mathbb{P}(G=g|G+e \in [1,T])) = \frac{1}{\sqrt{n}} \sum_{i=1}^{n} \xi^{g,e}(W_i) + o_p(1), \end{equation} with $\mathop{}\!\mathbb{E}[\xi^{g,e}(W)] = 0$ and $\mathop{}\!\mathbb{E}[\xi^{g,e}(W) \xi^{g,e}(W)'] < \infty$ being positive definite, and \begin{equation*} \xi^{g,e}(W) = \frac{1}{\mathop\!\mathbb{P}(G+e \in [1,T])} \cdot \bigg[ 1\{G=g, G+e \in [1,T]\} - \mathop\!\mathbb{P}(G=g| G+e \in [1,T]) \cdot 1\{G+e \in [1,T] \} \bigg] \end{equation*} The following corollary can be used to conduct asymptotically valid (pointwise) inference for the event-study type parameter $ES(e)$. \begin{corollary} Under the assumptions of Theorem (ref), for each $e$ such that $\mathop{}\!\mathbb{P}(1 \leq G+e \leq T)$, as $n\rightarrow \infty$, \begin{align*} \sqrt{n}(\widehat{ES}_{\text{dr},\text{gmm}}(e) - ES(e)) &= \frac{1}{\sqrt{n}} \sum_{i=1}^{n} l^{es,e}_{\text{gmm}}(W_i) + o_p(1)\\ & \overset{d}{\rightarrow} N(0, \mathop\!\mathbb{E}[l^{es,e}_{\text{gmm}}(W)^2]), \end{align*} with $l^{es,e}_{\text{gmm}}(W) = \sum_{g \in \mathcal{G}_{\text{trt}}} \Big( \mathop{}\!\mathbb{P}(G=g| G+e \in [1,T]) \cdot \dfrac{\textbf{1}'{\Omega}^{-1}_{g,t}}{\textbf{1}'{\Omega}^{-1}_{g,t}\textbf{1}} \psi^{g,t}(W_i; \kappa^{g,t}_{0}) + \xi^{g,e}(W_i) \cdot ATT(g,t) \Big)$. \end{corollary} The results in Corollary (ref) also apply to estimators of $ES(e)$ using $\widehat{ATT}_{\text{dr},g_c}(g,g+e)$ on (ref). Corollary (ref) focuses on pointwise inference procedures. Still, as discussed in Remark (ref), it is straightforward to extend it to hold for all event-times $e$ and conduct simultaneous-based inference. The asymptotic results for our overall summary parameter $\widehat{ES}_{\text{avg},\text{gmm}}$ as defined in (ref) follow from the delta method and are omitted. \section{Monte Carlo Simulations} In this section, we evaluate the finite sample properties of our proposed DR DDD estimators via Monte Carlo simulations. We examine two scenarios: (i) when covariates play a crucial role in identification across two time periods, and (ii) when there are multiple time periods with variation in treatment timings. For the first scenario, we have panel data for two time periods, $t=1,2$, four covariates, two enabling-groups $S \in \mathcal{S}_{des-1} \equiv \{2, \infty\}$, and there are two eligibility groups: $Q=1$ and $Q=0$. In the setup with staggered adoption, we consider the simplest possible case with three time periods, $t=1,2,3$, with $S \in \mathcal{S}_{des-2} \equiv \{2,3,\infty\}$, and we abstract from covariates in the main text. We relegate simulation results with DDD staggered adoption with covariates to the Supplemental Appendix. In the main text, we compare the performance of different DDD estimators via graphs: one that presents the density of the point estimates across the 1,000 Monte Carlo repetitions, and one that presents the length of confidence intervals in each Monte Carlo draw. In the Supplemental Appendix, we also report the traditional summary statistics for the Monte Carlo involving average bias, root mean square error (RMSE), empirical 95% coverage probability, and the average length of a 95% confidence interval under 1,000 Monte Carlo repetitions. Light-gray confidence intervals mean that they do not contain the true parameter of interest, $ATT(2,2)$ in our simulations, and are appropriately colored when they contain it. We focus on results with $n=5,000$ but report results for different sample sizes in the Supplemental Appendix (ref). \subsection{Simulations for DDD with two periods and covariates} We describe the data-generating process (DGP) for the 2-period DDD setup. For a generic four-dimensional vector $O$, the conditional probability of each unit belonging to a subgroup $(g,q) \in \{2,\infty\} \times \{0,1\}$ is \begin{equation} \mathop\!\mathbb{P}[S = g, Q=q | O] \equiv p^{S=g, Q=q}(O) = \dfrac{\exp(f^{ps}_{S=g,Q=q}(O))}{\sum_{(g,q) \in \mathcal{S}_{\text{des-1}} \times \{0,1\}}\exp(f^{ps}_{S=g,Q=q}(O))}, \end{equation} where $f^{ps}_{S=g,Q=q}(O))$ is a linear index with heterogeneous coefficients across sub-groups; we defined these in the Supplemental Appendix (ref) to save space. Of course, each unit belongs to a single subgroup, and we assigned these subgroups as follows: \begin{equation} (S,Q) := \begin{cases} (\infty,0), & \text { if } U \leqslant p^{S=\infty, Q=0}(O), \\ (\infty,1), & \text { if } p^{S=\infty, Q=0}(O) < U \leq \sum_{j=0}^1 p^{S=\infty, Q=j}(O), \\ (2,0), & \text { if } \sum_{j=0}^1 p^{S=\infty, Q=j}(O) < U \leq \sum_{j=0}^1 p^{S=\infty, Q=j}(O) + p^{S=2, Q=0}(O), \\ (2,1), & \text { if } \sum_{j=0}^1 p^{S=\infty, Q=j}(O) + p^{S=2, Q=0}(O) < U,\end{cases} \end{equation} with $U$ being a uniform random variable in $[0,1]$, independent of all other variables. The potential outcomes are defined as \begin{align} Y_{i,1}(\infty) &= f^{reg}(O_{i}, S_i) + \nu_i(O_i, S_i, Q_i) + \varepsilon_{i,1}(\infty) \nonumber\\ Y_{i,2}(\infty) &= 2 f^{reg}(O_{i}, S_i)+ \nu_i(O_i, S_i, Q_i) + \varepsilon_{i,2}(\infty) \\ Y_{i,2}(2) &= 2 f^{reg}(O_{i}, S_i)+ \nu_i(O_i, S_i, Q_i) + \varepsilon_{i,2}(2),\nonumber \end{align} where $f^{reg}(O_{i}, S_i)$ is a linear regression specification with heterogeneous coefficients across the enabling groups $S$, $\nu_i(O_i, S_i, Q_i)$ is a time-invariant unobserved heterogeneity correlated with covariates and sub-groups, and $\varepsilon_{i,1}(\infty), \varepsilon_{i,2}(\infty)$ and $\varepsilon_{i,2}(2)$ are independent standard normal random variables; we provide a precise definition of $f^{reg}(O_{i}, S_i)$ and $\nu_i(O_i, S_i, Q_i)$ in the Supplemental Appendix (ref). Note that our designs' $ATT(2,2)$ equals zero, though there is treatment effect heterogeneity across units. We observe untreated outcomes for all units in period $t=1$; in period $t=2$, we observed $Y_{i,2}(2)$ if unit $i$ belongs to group $S=2$, $Q=1$, and observe $Y_{i,2}(\infty)$ otherwise. Building on kang_schafer_2007 and SantAnna2020, we allow propensity score and/or outcome regression models to be misspecified. We consider four different types of DGP: DGP 1, where all models are correctly specified; DGP 2, where outcome models are correctly specified but the propensity score model is misspecified; DGP 3, where the propensity score is correctly specified but outcome regressions are not; and DGP 4, where all models are misspecified. The source of misspecification in these nuisance models is related to whether they depend on $X$ or $Z$, where $X$ is a nonlinear transformation of all the $Z$'s. In our simulations, the observed data is $W_i = \{Y_{i,1}, Y_{i,2}, S_{i}, Q_{i}, X_{i}\}_{i=1}^{n}$, so using $Z$ as linear covariates in these nuisance models lead to working model misspecification; we relegate to the Supplemental Appendix (ref) the definition of $X$'s and $Z$'s. We compare the performance of four different estimators for $ATT(2,2)$, just like in Section (ref): our DR DDD estimator as defined in (ref) (we label it as DRDDD), 3WFE OLS estimator for $\tilde{\beta}_{\text{3wfe}}$ based on (ref) (we label it as 3WFE), 3WFE OLS estimator of $\check{\beta}_{\text{3wfe}}$ based on (ref) (we label it as M-3WFE), and the difference of two SantAnna2020' DR DiD estimators (we label it as DRDID-DIF). We summarize the results of our simulations in Figure (ref), where we consider a sample size $n=5,000$ and conducted 1,000 Monte Carlo repetitions. The left panels display the density of the point estimates across all Monte Carlo draws, and the right panels display the 95% confidence intervals for each Monte Carlo draw. See Table (ref) in the Supplemental Appendix for additional results. \begin{figure}[htp] \begin{subfigure}[t]{0.9\textwidth} \caption{DGP 1: All working models are correctly specified} \end{subfigure} \begin{subfigure}[t]{0.9\textwidth} \caption{DGP 2: Outcome working models are correctly specified} \end{subfigure} \begin{subfigure}[t]{0.9\textwidth} \caption{DGP 3: Propensity-score working models are correctly correct} \end{subfigure} \begin{subfigure}[t]{0.9\textwidth} \caption{DGP 4: All working models are misspecified} \end{subfigure} \caption{Monte Carlo Simulation Results for DDD: two-period setup with covariates} \justifying \scriptsize{Notes: Simulation designs as discussed in text, with $n=5,000$ and $1,000$ Monte Carlo repetitions. True $ATT(2,2)$ is zero and is indicated in the solid vertical line in all panels. 3WFE corresponds to the OLS estimates of $\tilde{\beta}_{\text{3wfe}}$ based on (ref). M-3WFE corresponds to the OLS estimates of $\check{\beta}_{\text{3wfe}}$ based on (ref). DRDID-DIF corresponds to the difference between two doubly robust DiD estimators SantAnna2020. DRDDD corresponds to our proposed doubly robust DDD estimator described in (ref). All densities (left) and confidence intervals (right) are computed across all simulation draws. Light grey areas in the right plots indicate confidence intervals that exclude the true $ATT(2,2)$, where increased prominence suggests lower empirical coverage. } \end{figure} The results in Figure (ref) are self-explanatory and firmly support our theoretical results. When outcome regression or propensity score models (but not necessarily both) are correctly specified, our DR DDD estimators are appropriately centered (so they are unbiased), their confidence intervals are the narrowest across all other estimators, and they still have appropriate coverage across the first three DGPs (94.4%, 94.5%, and 94.6%, respectively). For instance, when all working models are correctly specified, the average length of the 95% confidence interval of the M-3WFE estimator is 7 times longer than our DR DDD estimator; this difference is much larger for the other considered estimators. In fact, the performance of all other estimators in all our considered DGPs is poor, as they have non-negligible bias, high RMSE, and poor coverage properties. In DGP 4, when all working models are misspecified, we note that our estimator is biased, directly affecting the confidence intervals' coverage probabilities. None of the considered DDD estimators perform well when all working models are misspecified (DGP 4), highlighting that all estimators indeed depend on modeling assumptions. \subsection{Simulations for DDD with variation in treatment timing} We now discuss the staggered DDD setup with three time periods, $t=1,2,3$, three enabling groups, $S \in \{2,3,\infty\} $, and two eligibility groups $Q \in \{0,1\}$. We abstract from covariates and defer a discussion about them to the Supplemental Appendix. Each unit $i$, we have that $p^{S=2,Q=0}=0.20$, $p^{S=2,Q=1}=0.15$, $p^{S=3,Q=0}=0.30$, $p^{S=3,Q=1}=0.20$, $p^{S=\infty,Q=0}=0.05$, and $p^{S=\infty,Q=1}=0.10$. We then randomly assign the realized value of $(S,Q)$ based on the above distribution. The potential outcomes are generated as \begin{align} Y_{i,1}(\infty) &= (1 + Q_i)\alpha + \nu_i(S_i, Q_i) + \varepsilon_{i,1}(\infty) \nonumber\\ Y_{i,2}(\infty) &= (2 + Q_i) \alpha + 1.1 \nu_i(S_i, Q_i) + \varepsilon_{i,2}(\infty) \nonumber\\ Y_{i,3}(\infty) &= (3 + Q_i) \alpha + 1.2 \nu_i(S_i, Q_i) + \varepsilon_{i,3}(\infty) \nonumber\\ Y_{i,2}(2) &= (2 + Q_i) \alpha + 1.1 \nu_i(S_i, Q_i) + ATT(2,2) Q_i + \varepsilon_{i,2}(2) \\ Y_{i,3}(2) &= (3 + Q_i) \alpha + 1.2 \nu_i(S_i, Q_i) + ATT(2,3) Q_i + \varepsilon_{i,3}(2)\nonumber\\ Y_{i,3}(3) &= (3 + Q_i) \alpha + 1.2 \nu_i(S_i, Q_i) + ATT(3,3) Q_i + \varepsilon_{i,3}(3)\nonumber, \end{align} where we set $\alpha = 278.5$, $ATT(2,2)=10$, $ATT(2,3)=20$, $ATT(3,3)=25$, all $\varepsilon_{i,t}(\cdot)$ are independent standard normal, and $\nu_i(S_i, Q_i)$ is an unobserved heterogeneity term that we formally define in the Supplemental Appendix (ref). The observed data is represented as $W_i = \{Y_{i, 1}, Y_{i,2}, Y_{i,3}, S_{i}, Q_{i}\}_{i=1}^{n}$, where the $Y_{i,1} = Y_{i,1}(\infty)$, $Y_{i,2} = Y_{i,2}(2)$ for units with $S_i=2,Q_i=1$, and $Y_{i,2} = Y_{i,2}(\infty)$ otherwise, and $Y_{i,3} = Y_{i,3}(3)$ for units with $S_i=3,Q_i=1$, $Y_{i,3} = Y_{i,3}(2)$ for units with $S_i=2,Q_i=1$, and $Y_{i,3} = Y_{i,3}(\infty)$ for all other units. In what follows, we focus our attention on estimators for $ATT(2,2)$, though we compare estimators for $ATT(2,3)$ and $ATT(3,3)$ in the Supplemental Appendix (ref). We consider $n=5,000$ here and refer to Supplemental Appendix (ref) for information on other sample sizes. We compare the performance of three staggered DDD estimators for the $ATT(2,2)$ as in Section (ref). More precisely, we consider our optimally GMM-weighted estimator $\widehat{ATT}_{\text{gmm}}(g,t)$ as defined in (ref), which is a special case of the (ref) without covariates, i.e., with $X_i=1$ for all units; we refer to this estimator as $DDD_{\text{gmm}}$. We also consider our DDD estimator that uses never-treated units as the comparison group as defined in (ref) with $ g_{\text{c}} = \infty$; we refer to this estimator as $DDD_{\text{nev}}$.\footnote{This estimator is a special case of (ref) when covariates are trivial}. Finally, we consider a DDD estimator that pools all not-yet-treated units a la Callaway2021, $\widehat{ATT}_{\text{cs-nyt}}(g,t)$, as defined in (ref); we refer to this estimator as $DDD_{\text{cs-nyt}}$. Our theoretical results indicated that $DDD_{\text{gmm}}$ and $DDD_{\text{nev}}$ should both be valid DDD estimators under our identification assumptions. At the same time, we have no statistical guarantee about the performance of $DDD_{\text{cs-nyt}}$. Our theoretical results also suggest that $DDD_{\text{gmm}}$ should be more precisely estimated than $DDD_{\text{nev}}$, as it uses more information. \begin{figure}[htp] \begin{subfigure}[t]{\textwidth} \caption{Comparison across the three DDD estimators} \end{subfigure} \begin{subfigure}[t]{\textwidth} \caption{Comparing the performance of $DDD_{\text{{nev}}}$ and $DDD_{\text{{gmm}}}$} \end{subfigure} \caption[position=bottom]{Monte Carlo Simulation Results for DDD: staggered setup without covariates} \justifying \scriptsize{Notes: Simulation designs are detailed in the text, using $n=5,000$ and 1,000 Monte Carlo repetitions. The true $ATT(2,2)$ is 10, marked by a solid vertical line in all panels. $DDD_{\text{{nev}}}$ denotes our DDD estimator with $g_c = \infty$ from Equation (ref). $DDD_{\text{{gmm}}}$ is our proposed DDD estimator from Equation (ref). $DDD_{cs-nyt}$ is the estimator pooling all not-yet-treated units as defined in (ref). In the top-left panel, we plot the densities of the $DDD_{\text{{nev}}}$, $DDD_{\text{{gmm}}}$, and $DDD_{cs-nyt}$ estimates across all simulation draws. The bottom-left panel zooms into the densities for $DDD_{\text{{nev}}}$ and $DDD_{\text{{gmm}}}$. The top-right panel shows the confidence intervals for these estimators, while the bottom-right panel focuses on $DDD_{\text{{nev}}}$ and $DDD_{\text{{gmm}}}$ only. Light grey areas in the right plots indicate confidence intervals that exclude the true $ATT(2,2)$, where increased prominence suggests lower empirical coverage. } \end{figure} Figure (ref) summarizes our simulation results; see also Table (ref) in the Supplemental Appendix. Panel (a) of Figure (ref) compares the performance of the three estimators. As it is evident, $DDD_{\text{cs-nyt}}$ is severely biased in our design, and its confidence interval never covered the true $ATT(2,2)$ in the Monte Carlo draws. This result further highlights that proceeding as if DDD is just the difference of two DiD procedures can lead to misleading conclusions. At the same time, the simulations highlight that $DDD_{\text{{nev}}}$ and $DDD_{\text{{gmm}}}$ are unbiased and have good coverage properties, $93.2\%$ and $95\%$, respectively. In Panel (b), we drop $DDD_{\text{cs-nyt}}$ and focus on comparing the performance of our two proposed DDD estimators. As the plot makes it clear, the gains in efficiency of using all not-yet-treated comparison groups as in $DDD_{\text{{gmm}}}$ are notable. The average length of the confidence intervals of $DDD_{\text{{nev}}}$ is more than $50\%$ higher than that of $DDD_{\text{{gmm}}}$. Thus, in practice, we recommend favoring $DDD_{\text{{gmm}}}$, especially when the sample size of the never-enabling comparison group is low. \section{Empirical Illustrations} In this section, we revisit (i) cai_insurance_2016 and analyze the effect of agriculture insurance provision on household financial decisions, (ii) carbon_pricing and analyze the effect of the emission trading scheme implemented in China on carbon emissions, and (iii) hansen_national_2023 and assess the impact of genetically modified crop adoption on countrywide yields. These three studies examine DDD setups with a single treatment date and covariates, as well as staggered DDD setups with and without covariates. The aim is to illustrate how our proposed estimator performs compared to the methods used in these original studies. Overall, we find that our DR DDD procedure can (a) produce substantial gains in precision, with 3WFE 95% confidence intervals being up to 115% wider than ours, (b) lead to different conclusions about the effect of a policy, or (c) provide more substantial evidence supporting the effectiveness of a policy. \subsection{Effect of insurance provision on financial decisions} cai_insurance_2016 analyzes the effects of agricultural insurance programs implemented in rural Chinese households on their production, borrowing, and saving behaviors. These rural areas are particularly vulnerable to unique weather shocks, resulting in income volatility for households engaged in agriculture. To address this, the People's Insurance Company of China (PICC) launched the first weather-indexed crop-insurance program in 2003 for tobacco farmers in selected counties of Jiangxi province. The contracts were compulsory for tobacco growers in the treatment regions, allowing cai_insurance_2016 to take advantage of the variation in insurance provision across regions and household types. By design, the policy aimed to simultaneously protect farm income and sustain county fiscal revenues that heavily depend on tobacco taxes. The insurance policy rollout introduces three sources of variation: the timing of insurance provision (pre-policy, February 2000, and post-policy, August 2003), geographic regions (counties that enabled policy vs. those that did not), and household eligibility conditions (tobacco-specialized vs. other farm or off-farm households). This setup enables us to allow for industry-specific and region-specific trends. cai_insurance_2016 utilizes an administrative panel from the Rural Credit Cooperative (RCC), encompassing data from over 5,000 households annually surveyed between 2000 and 2008. This dataset merges financial and savings records with RCC survey details on demographics, land allocations, and crop revenues. It includes more than 3,000 tobacco-growing households (approximately 1,200 of which are located in treated counties) and 2,200 non-tobacco households. The study evaluates outcomes such as the area of tobacco production, loan size, loan limits, interest rates, savings rates, and flexible-term savings\footnote{Although the original study analyzed outcomes related to production, borrowing, and saving behaviors, the DDD design specifically estimates the impact of insurance provision on borrowing and saving behavior, accounting for region-specific effects.}. We adopt a similar approach as in cai_insurance_2016 by employing a DDD design to evaluate how insurance provision affects household saving behaviors. cai_insurance_2016 considered the following 3WFE event study specification\footnote{See, e.g., Equations (4) and (6) in cai_insurance_2016.} to estimate dynamic average treatment effects of insurance provision on household saving behavior \begin{equation} Y_{i,t} = \gamma_i + \gamma_{r,t} + \gamma_{j,t} + \sum_{e \neq-1} \beta_{e} \mathbf{1}\{E_{i,t}=e\} + X'_{i,j,r}\theta + u_{i,t}, \end{equation} where where $i,j,r,t$ are household, sector, region, and years indices, respectively, $\gamma_i, \gamma_{r,t}, \gamma_{j,t}$ denote household fixed effects, region-by-time fixed effect and industry-by-time fixed effect, respectively, $E_{i,t} = t - G_i$ is the time relative to the period at which household $i$ has received the insurance\footnote{In this setup, household $i$ is considered to have received the treatment if it is located in a treated region and qualifies as a tobacco household.} and $u_{i,t}$ is an idiosyncratic error term. Covariates in $X_{i,j,r}$ include household size, head of household age, and education level. Regarding outcome variables $Y_{i,t}$ for measuring household saving decisions, cai_insurance_2016 examines net saving, defined as the annual increase in total savings; saving rate, defined as the ratio of net saving to current household income; and flexible-term saving, which is the ratio of net savings in checking accounts to the total net savings. We compare the estimates of $\beta_{e}$ with those of $ES(e)$, as specified in (ref), using our DR procedure, using the same set of covariates. In our setup, $S$ represents the geographic region that enables or not the insurance policy, and $Q$ represents the household eligibility condition---i.e., tobacco-specialized or not. \begin{figure}[pht] \begin{subfigure}[t]{0.8\textwidth} \caption{Panel A: Insurance Provision on Net saving} \end{subfigure} \begin{subfigure}[t]{0.8\textwidth} \caption{Panel B: Insurance Provision on Saving rate} \end{subfigure} \begin{subfigure}[t]{0.8\textwidth} \caption{Panel C: Insurance Provision on Flexible-term saving} \end{subfigure} \caption{Impact of Insurance Provision on Savings Decisions.} \justifying \scriptsize{Notes: This figure presents event study estimates with covariates including household size and head of household age. Panel A uses Net Saving as outcome, Panel B focuses on the Saving Rate, and Panel C examines Flexible-term Saving. In each panel, blue dots represent the point estimates from the specification in (ref), while red triangles show the $\widehat{ES}_{dr, g_{\text{c}}}(e)$ estimates aggregated from our proposed DR estimator as in (ref). The average $\widehat{ES}_{\text{avg}}$, calculated from all $\widehat{ES}_{dr, g_{\text{c}}}(e)$ values to the right of the vertical dashed line at \(e = -1\), is provided alongside its bootstrapped standard errors, based on 999 repetitions, and confidence interval. The corresponding blue and red vertical lines represent the 95% confidence intervals. } \end{figure} Figure (ref) plots the event study estimates for the 2003 weather insurance policy. Both the 3WFE and our DR DDD procedure yield similar results in terms of magnitude. However, our estimator offers a precision advantage in this particular context, with 3WFE 95% confidence intervals being up to 1.15 times wider than the one associated with our proposed DR DDD estimator. In Panel A of Figure (ref), we observe an insignificant effect on net saving, in line with the results reported in cai_insurance_2016. In Panel B, we observe a negative effect on the saving rate, consistent with previous findings by cai_insurance_2016. However, unlike earlier results, which indicated an insignificant negative effect, our proposed estimator reveals a small yet statistically significant negative effect at the 95% confidence level. This suggests a change in savings behavior after obtaining insurance. The findings in Panel C, consistent in magnitude with those from cai_insurance_2016, demonstrate that weather insurance significantly increases the proportion of flexible-term savings. This suggests that households prefer to shift their savings to more flexible instruments after acquiring insurance. Visually, results in Panel C affirm the previous findings in cai_insurance_2016, though we emphasize the enhanced precision our method provides. For example, the confidence intervals associated with $\widehat{\beta}_e$ are 16.29%, 85.53%, 64.21%, and 115.84% wider than the confidence intervals associated with $\widehat{ES}(e)$ at event times $e=2,3,4,5$, respectively. This illustrates that our DDD procedure can provide much more powerful analysis than alternative procedures. \subsection{Effect of emission trading scheme on carbon emissions} We next revisit the study by carbon_pricing, which examines the impact of the emission trading scheme (ETS) implemented by the Chinese government to regulate and reduce carbon emissions at minimal cost, focusing on its effect in promoting low-carbon innovation among firms. China tackled its dual objectives of economic growth and reducing carbon emissions by approving seven regional ETS pilot programs in 2011. These pilot programs spanned four major cities, two provinces, and one special economic zone, each developing their own allowance and enforcement regulations within national guidelines, to peak national carbon emissions by around 2030. carbon_pricing leverages three sources of variation. First, it considers firms' patent applications before and after the introduction of ETS pilots. Second, the pilots were launched in different regions whereas some regions are not impacted by the introduction of an ETS. Third, the pilots covered various manufacturing sectors within these regions such as power and heating, chemical, cement etc. Therefore, the authors utilize these variations in ETS pilots over time, across sectors, and across regions, employing a DDD approach to identify the effects of ETS on low-carbon innovation. This setup exemplifies a DDD model with variation in treatment timing and covariates. carbon_pricing uses the year of announcement (i.e., 2011) rather than the actual launch year as the post-policy indicator. The regional carbon market pilots were phased in starting from 2013, with initial launches in Shenzhen, Shanghai, Beijing, Guangdong, and Tianjin, followed by the introduction in Hubei and Chongqing the subsequent year. Thus, we adjust the construction of the treatment variable to use the official launch date instead of the announcement date, allowing for staggered treatment adoption.\footnote{We also considered the original specification using treatment announcement as the definition of the treatment. In such a case, the results obtained using our DR DDD procedure and those in carbon_pricing are very close; these results are available upon request. Thus, we stress that the empirical results in this section should be interpreted as an extension of carbon_pricing, and not a “replication”.} The dataset utilized in carbon_pricing includes information on publicly listed firms in China from 2003 to 2015, sourced from the Shanghai and Shenzhen stock exchanges. This dataset comprises both financial data, obtained through the China Stock Market and Accounting Research (CSMAR), and patent application details, collected from China's State Intellectual Property Office (SIPO). The dataset includes $18,937$ firm-year observations from $1,956$ firms spanning the years 2003 to 2015. The three outcomes evaluated by carbon_pricing are the number of low-carbon patent applications, the number of patent applications in other non-low-carbon technologies, and the ratio of low-carbon patents relative to the total number of patents. In terms of estimators, we use a 3WFE similar to carbon_pricing\footnote{See, e.g., Equation (1) in carbon_pricing}, for firm $i$ at year $t$ \begin{equation} Y_{i,t} = \gamma_i + \gamma_{r,t} + \gamma_{s,t} + \beta_{3wfe} D_{i,t} + X'_{i,s,r}\theta + \varepsilon_{i,t}, \end{equation} where $\beta_{3wfe}$ is meant to capture the impact of ETS on low-carbon technology innovation, $D_{i,t}$ is the treatment indicator post-announcement of ETS for covered sectors in pilot regions, $\gamma_i$ is a firm fixed effect, $\gamma_{r,t}$ and $\gamma_{s,t}$ are region-by-time and sector-group-by-time fixed effects, respectively. Covariates $X_{i,s,r}$ include firm-level attributes such as assets, revenue, and current liabilities. As we are also interested in treatment effect dynamics, we also consider a 3WFE event-study specification \begin{equation} Y_{i,t} = \gamma_i + \gamma_{r,t} + \gamma_{s,t} +\sum_{e \neq-1} \beta_{e} \mathbf{1}\{E_{i,t}=e\} + X'_{i,s,r}\theta + u_{i,t}, \end{equation} where $E_{i,t} = t - G_i$ denotes the time relative to when the regional ETS pilot was initiated for firm $i$, and $u_{i,t}$ is an idiosyncratic error term. In this setup, a firm $i$ is considered to have received the treatment if it is located in a region that is a carbon market pilot and its sector qualifies as a regulated sector by the ETS. We compare these 3WFE estimators with our DR DDD estimator (ref). In our setup, $S$ represents the year a region launched the ETS pilot, and $Q$ represents the manufacturing sector that is eligible for ETS. \begin{figure}[!htp] \begin{subfigure}[t]{0.75\textwidth} \caption{Panel A: Emission Trading Scheme regulation on low-carbon patents} \end{subfigure} \begin{subfigure}[t]{0.75\textwidth} \caption{Panel B: Emission Trading Scheme regulation on other non-low-carbon patents} \end{subfigure} \begin{subfigure}[t]{0.75\textwidth} \caption{Panel C: Emission Trading Scheme regulation on share of low-carbon patents} \end{subfigure} \caption{Impact of ETS on share of low-carbon patents.} \justifying \scriptsize{Notes: This figure displays event study estimates incorporating firm characteristics such as total assets, total revenue, and current liabilities as covariates. Panel A examines the number of low-carbon patents, Panel B considers other non-carbon patents, and Panel C uses the proportion of low-carbon patents relative to the total number of patents. In each panel, blue dots represent the point estimates from the specification in (ref), while red triangles show the $\widehat{ES}_{\text{dr},\text{gmm}}(e)$ estimates aggregated from our proposed DR estimator as in (ref). The average $\widehat{ES}_{\text{avg}, \text{gmm}}$, calculated from all $\widehat{ES}_{\text{dr},\text{gmm}}(e)$ values to the right of the vertical dashed line at \(e = -1\), is provided alongside its bootstrapped standard errors, based on 999 repetitions, and confidence interval. The corresponding blue and red vertical lines represent the 90% confidence intervals. } \end{figure} Figure (ref) compares the event study estimates from the specification (ref) with the estimates $\widehat{ES}_{\text{dr},\text{gmm}}(e)$, as specified in (ref), using our DR procedure. Consistent with carbon_pricing, we present results at a 90% confidence level. In Panel A, as reported by carbon_pricing, the ETS pilot appears to encourage innovation, but we cannot reject the null hypothesis of no effect at the same confidence level. Panel B's results align with carbon_pricing, indicating minimal evidence for the ETS's crowding-out effect on other patent applications. The findings in Panel C show no significant impact of ETS regulation on the proportion of low-carbon patents. This contrasts with the results associated with a 3WFE, which indicated a positive and significant effect at a 90% confidence level\footnote{See Table 1 in carbon_pricing}. For instance, our method estimates an aggregate effect of 1.7% without sufficient evidence to reject the null hypothesis of no effect, while specification (ref) estimates a significant aggregate effect of 5.3% of ETS adoption on share of low-carbon patents. Overall, our findings suggest there is no evidence to support that the ETS promotes innovation in low-carbon patents. Additionally, our tighter confidence intervals reflects that DR DDD can provide practical gains in power in empirical applications. For example, in Panel C of Figure (ref), the confidence intervals associated with $\widehat{\beta}_e$ are 20.28%, 24.26% and 1.51% wider than the confidence intervals associated with $\widehat{ES}_{\text{dr},\text{gmm}}(e)$ at event times $e=0,1,2$, respectively. This implies that $\widehat{\beta}_e$ would need a sample with 44%, 54%, and 3% more observations than our DR DDD event study estimator for event times zero, one, and two to achieve the same precision as DR DDD.\footnote{These calculations are based on asymptotic relative efficiency. For any parameter $\eta$ of a distribution $F$, and for estimators $\widehat{\eta}_{1}$ and $\widehat{\eta }_{2}$ approximately $N\left( \eta,V_{1}/n\right) $ and $N\left( \eta, V_{2}/n\right) $, respectively, the asymptotic relative efficiency of $\widehat{\eta}_{2}$ with respect to $\widehat{\eta}_{1}$ is given by $V_{1}/V_{2}$; see, e.g., Section 8.2 in VanderVaart1998.} \subsection{Effect of genetically modified crops on yields} Finally, we revisit the study by hansen_national_2023 on the impact of genetically modified (GM) crop adoption on countrywide yields. While farmers in North and South America and parts of Asia adopted these seeds almost immediately, regulators in the European Union, much of Africa, and several middle‑income economies imposed outright bans or stringent de-facto restrictions, often in response to food-safety scares and NGO pressure. By 2019 only 29 countries permitted commercial cultivation, another 42 permitted imports but not planting, and the rest maintained comprehensive bans. This sharp global policy divergence, paired with the agronomic premise that pest‑ and herbicide-tolerant GM varieties raise yields mainly where weeds and insects are prevalent, sets the stage for a compelling cross‑country analysis of aggregate production impacts. hansen_national_2023 exploits the fact that GM technologies were adopted at different times in various countries starting in the mid-1990s to examine the effects of GM adoption on outcomes such as crop yields, harvested area, and trade flows, using balanced crop-country panel data covering over 120 countries from 1986 to 2019. Different sources of variation underpin the DDD design in this context. First, there was a variation in policy timing as each country passed GM cultivation legislation in different years. This allows one to compare outcomes before and after GM commercialization is allowed, between countries that have already passed these legislations and those that have not yet. Second, crop eligibility for GM traits exists for only four field crops---cotton, maize, soybean, and rapeseed--- leaving rice, wheat, and dozens of other staples as the primary comparisons. Differencing outcomes over time, across crops, and across policy regimes therefore allows for common shocks (e.g., weather, prices) and crop‑specific global trends, isolating the causal effect of GM adoption. hansen_national_2023 utilizes publicly available country-level data on production quantities from FAOSTAT, covering yields (production per hectare), producer prices, and trade flows (exports and imports by commodity). Following their setting, our focus is also on countries with a minimum of $100,000$ hectares of cropland, thereby drawing attention to nations with significant agricultural activities. Similarly, we also include 60 field crops, including cereals, pulses, roots and tubers, oil crops, and fiber crops, as outlined by the FAO FAO2012. hansen_national_2023 collected official legislative documents, USDA Foreign Agricultural Service reports, and industry publications to determine the earliest legal and feasible year of commercial cultivation for each country-crop combination.\footnote{For further details, refer to the documentation for GM approval dates in hansen_national_2023 Online Appendix B.} In some cases, GM approval did not directly result in large-scale planting, while in others, formal bans restricted or effectively prohibited cultivation. Thus, following hansen_national_2023, we establish $G_{i,c}$ as the first year in which a GM variety of a given crop in a given country could be harvested and sold commercially for human consumption or animal feed without violating a ban.\footnote{hansen_national_2023 uses a different notation than us, and refers to our $G_{i,c}$ as $E_{ic}$.} hansen_national_2023 considers the following event-study specification\footnote{See, e.g., Equation (1) in hansen_national_2023.} to asses the dynamic average effect of GM adoption on countrywide yields \begin{equation} \ln Y_{i,c,t} = \gamma_{c,i} + \gamma_{i,t} + \gamma_{c,t} + \sum_{e\not=-1} {\beta}_{e} \cdot \mathbf{1}\{E_{i,c,t} = e\} + \varepsilon_{i,c,t} \end{equation} where $Y_{i,c,t}$ denotes the yield of crop $c$ in country $i$ for the year $t$. The variable $E_{i,c,t} = t - G_{i,c}$ represents the time since a GM variety of a crop became legally permissible for harvest and sale for human consumption, and $\varepsilon_{i,c,t}$ is an idiosyncratic error error. The specification in (ref) includes three fixed effects: country-by-year fixed effects ($\gamma_{i,t}$) to account for shocks affecting all crops in a particular country and year, crop-by-country fixed effects ($\gamma_{c,i}$) to capture local time-invariant conditions for each crop, and crop-by-year fixed effects ($\gamma_{c,t}$) to reflect global trends specific to each crop. We compare event-study estimates based on (ref) with our proposed DR DDD-based event study in (ref). In our setup, $S$ represents the time a country allowed for GM commercialization, and $Q$ represents the eligible crop. Figure (ref) illustrates the event study estimates of the impact of GM crop adoption on yields of cotton, maize, rapeseed, and soybean together. Panel A shows results aligned with those in hansen_national_2023, with a slight enhancement in the precision of the confidence intervals. For instance, at event times $e = 0$, the confidence intervals associated with $\widehat{\beta}_{e}$ is 46.47% wider than the one associated with $\widehat{ES}_{\text{dr}, \text{gmm}}(e)$. Both methods show that approving GM technology leads to a significant increase in crop yields by approximately 13%. It is not surprising that both methods yield similar results in this context. Despite the challenges associated with 3WFE specifications with variation in treatment timing strezhnev2023, there exists a very large number of never-enabling crop-country combinations, making these “forbidden comparisons” less likely to shift results. Indeed, as noted by hansen_national_2023, only about 2 percent of the units are treated during the period of analysis. \begin{figure}[hptb] \begin{subfigure}[t]{0.8\textwidth} \caption{Panel A: GM adoption on log-yields (full sample)} \end{subfigure} \begin{subfigure}[t]{0.8\textwidth} \caption{Panel B: GM adoption on log-yields (dropping all never-treated)} \end{subfigure} \caption{Impact of GM adoption on yields} \captionsetup{justification=justified} \justifying \scriptsize{Notes: This figure shows event study estimates using the logarithm of crop yields as the outcome variable, without any covariates. The estimation window spans from 1986 to 2019, excluding country-crop combinations treated in 2010 or later. Panel A considers a fully balanced sample, as in hansen_national_2023, whereas Panel B excludes units never treated and forces the last cohort to be considered never-treated, with data post-treatment being filtered out. In each panel, blue dots represent point estimates from the specification in equation (ref), while red triangles show the $\widehat{ES}_{\text{dr},\text{gmm}}(e)$ estimates aggregated from our proposed DR estimator as in (ref). The average $\widehat{ES}_{\text{avg},\text{gmm}}$, calculated from all $\widehat{ES}_{\text{dr},\text{gmm}}(e)$ values to the right of the vertical dashed line at \(e = -1\), is provided alongside its bootstrapped standard errors, based on 999 repetitions, clustered at the country-crop level and its confidence interval. The corresponding blue and red vertical lines represent the 95% confidence intervals. } \end{figure} Having said that, a potential concern that can arise is related to whether it is indeed desirable to use these never-enabling crop-country combinations as a “de facto” comparison group. It may be the case that countries that never allow harvesting and commercialization of GM crops during a very long period of analysis may already face lower yields or stricter farm rules, so they might not be a necessarily good comparison group for eventually treated crop-countries observations. To assess the robustness of the findings in Panel A of Figure (ref), we drop data for every truly never-enabling unit, and re-run the 3WFE regression specification (ref) in this subsample as well as our staggered DDD estimator in (ref). As Panel B in Figure (ref) highlights, this has important consequences for the analysis: if one were to rely on (ref), one would find much smaller effects compared to those in Panel A, with non-negligible pre-treatment trends that may cast doubt about the plausibility of the DDD design. On the other hand, when using our proposed staggered DDD estimator, we can see that the results remain positive and statistically significant for most event-times, with no very serious violation of pre-treatment trends. In fact, our estimates for $ES_{\text{avg}}$ indicate an approximate 15% increase in crop yields following GM adoption, a slightly larger (but statistically distinguishable) than those in Panel A. Of course, statistical precision is less accurate in Panel B than in Panel A, as we use a much smaller sample size. This result essentially highlights that the main findings in hansen_national_2023 are robust to dropping the never-enabling set of units from the comparison group, as long as you use DDD procedures that are meant to work well in such setups. \section{Concluding remarks} This paper studied DDD estimators, paying close attention to situations where covariates are important for identification and to setups with staggered treatment adoption. Our findings challenge the conventional wisdom that DDD can be understood as the difference between two DiDs. We showed that when DDD-type parallel trends hold after conditioning on covariates, DDD estimators cannot generally be expressed as such, even in cases with only two time periods. In addition, when treatment adoption is staggered, pooling all not-yet-treated units is not generally valid, and proceeding as such can lead to misleading conclusions even when covariates are not crucial for the DDD identification arguments. These results highlight the need for more careful consideration when applying DDD strategies. To address these challenges, we proposed DR DDD estimators that can appropriately handle covariates and can also be used in DDD setups with staggered treatment adoption. Importantly, we proposed a DR DDD estimator that leverages information across different comparison groups, and our simulation results highlighted that the gains in precision can be substantial compared to alternatives. As such, we recommend practitioners to favor our proposed DDD estimator $\widehat{ATT}_{\text{dr},\text{opt}}(g,t)$ in applications. Finally, a companion \texttt{R} package, \href{https://marcelortiz.com/triplediff/}{\texttt{triplediff}}, is freely available on GitHub to automate all these DDD estimators proposed in this article, simplifying their adoption for practitioners. We envision extending our proposed DDD tools to incorporate data-adaptive and machine-learning estimators for the nuisance parameters, with a particular emphasis on training these models to ensure robust performance in finite samples. Additionally, we aim to generalize our DDD framework to accommodate more flexible sampling schemes, such as unbalanced panel data or repeated cross-sectional data; see, for example, Abadie2005, Callaway2021, and SantAnna2023 for related developments in DiD settings. Other promising avenues for future research include exploring semiparametric efficiency bounds and efficient estimation in overidentified DDD models Chen_SantAnna_Xie_effientDiD_2025, as well as using a DDD strategy to uncover persuasion effects Jun_Lee_2024_DiD_persuasion. We leave a full treatment of these topics to future work. { {1pt plus 0.3ex} \putbib }