Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
64,553 characters · 17 sections · 36 citation commands
Improved Inference for CSDID Using the Cluster Jackknife
The difference-in-differences (or DiD) model has been a workhorse of empirical economics for several decades. Until about five years ago, treatment effects for DiD models were typically estimated using two-way fixed effects (TWFE) regressions. Recently, a series of highly influential papers \citep*{de2020twott,callaway2021difference,goodman2021difference} has shown that TWFE regressions only recover the intended treatment effects when some very strong assumptions hold. Several of these assumptions fail to hold in a great many empirical applications. TWFE regressions often yield inconsistent estimates of the desired treatment effects because they involve “forbidden comparisons” and “negative weights.” New methods have therefore been proposed, and thousands of papers have been written using them. The focus of this new literature has been on estimating treatment effects correctly, but relatively little attention has been paid to inference. Notable exceptions are weiss2024much and mizushima2025inference. Perhaps the most widely-used of these new methods is the one proposed by callaway2021difference, which is often referred to as CSDID.
It has been known for many years that inference in TWFE regressions can be seriously unreliable when the disturbances are correlated within clusters. \citet*{BDM_2004} popularized this issue, identified serial correlation as one source of the problem, and suggested using cluster-robust standard errors to improve inferences. \citet*{CGM_2008} then showed that these standard errors can be unreliable when the number of clusters in a sample is not large and proposed the wild cluster bootstrap. The literature on cluster-robust inference for conventional DiD models and other regression models has exploded since then. \citet*{MNW-guide} is a recent guide.
Recently, it has been recognized that the cluster jackknife often provides a convenient and reliable way to obtain cluster-robust standard errors \citep*{MNW-bootknife,Hansen-jack,Hansen_2025}. The first of these papers also introduces a jackknife-inspired transformation to the wild cluster bootstrap, which can greatly improve inferences.
Given the difficulty of making reliable cluster-robust inferences with traditional TWFE regressions, it seems very likely that modern DiD methods, including the CSDID estimator, will also encounter difficulties. In (ref), we show that they do. To improve inference in finite samples, we therefore propose using the cluster jackknife. The jackknife substantially improves the reliability of inferences, but the problem of few treated clusters remains in some settings. To facilitate the estimation of jackknife standard errors, we provide both Stata and R packages.
The remainder of the paper is organized as follows. The first part of (ref) introduces the notation, while (ref) reviews the problem of staggered adoption for TWFE models. (ref) briefly reviews the CSDID estimator and the associated inference procedures. (ref) describes the two existing methods for cluster-robust inference with CSDID, and (ref) proposes the cluster jackknife. (ref) describes the Monte Carlo design to compare the various procedures, and (ref) discusses the results. (ref) briefly examines two empirical examples. (ref) describes the Stata and R packages. Finally, (ref) concludes.
The primary objective of most analyses using the difference-in-differences approach is to identify the average treatment effect on the treated ($\mbox{ATT}$). Following the potential outcomes framework, the $\mbox{ATT}$ is defined as:
where $Y_{i,r,t}(1)$ and $Y_{i,r,t}(0)$ represent potential outcomes for person $i$ in region $r$ in period $t$, with and without treatment, respectively, and $D_r$ denotes the binary treatment indicator, which is assigned at the region level. In the data, we only observe realized outcomes $Y_{i,r,t}$. While \hyperref[{eq:defATT}]{\tagform@{\ref*{eq:defATT}}} defines the target parameter, the appropriate estimator depends on the temporal distribution of treatment across units; specifically, on whether treatment occurs simultaneously or is staggered over time.
Consider a panel setting with regions $r \in \{1,\ldots,R\}$ observed over time periods $t \in \{1,\ldots,T\}$. We partition regions into those that receive treatment ($r \in R_1$) and those that are never treated ($r \in R_0$). To formalize the timing of these interventions, let ${\bm{S}} = (S_1, \ldots, S_R )$ be a $1 \times R$ vector where each entry $S_r$ denotes the first treatment period for region $r$. For never-treated regions, $S_r = 0$. We assume that treatment is an absorbing state; once a region is treated, it remains treated for all subsequent periods $t \geq S_r$.
The complexity of the estimation strategy is determined by the number of unique non-zero values in ${\bm{S}}$. If ${\bm{S}}$ contains only a single unique non-zero value, then the design involves common adoption. In that case, it is standard practice to use the two-way fixed effects estimator (TWFE) discussed in (ref). If, however, ${\bm{S}}$ contains multiple unique non-zero entries, the design involves staggered adoption.
In the case of staggered adoption, it is computationally and conceptually useful to group regions into treatment timing cohorts indexed by $g$. These cohorts are defined by the period of initial treatment adoption. Let $\mathcal{G}$ be the set of unique non-zero values in ${\bm{S}}$. For each $g \in \mathcal{G}$, we define a $1 \times R$ binary vector ${\bm G}_g$, where the $r$-th entry $G_{g,r} = 1$ if and only if $S_r = g$. Similarly, for the never-treated group, we define ${\bm G}_\infty$ such that $G_{\infty,r} = 1$ if and only if $S_r = 0$.
This cohort-based structure allows for the estimation of group-time average treatment effects (on the treated), denoted as $\mbox{ATT}(g,t)$. These parameters represent the effect of treatment for a specific cohort $g$ at time $t$. For simplicity, references to the cohort ${\bm G}_g$, or cohort $g$, imply the set of regions with values of 1 for that particular vector.
As discussed in (ref), the $\mbox{ATT}(g,t)$ serve as the building blocks for aggregate parameters, such as cohort-specific effects $\mbox{ATT}(g)$ or a single weighted average $\mbox{ATT}$ for the entire study period. This decomposition is particularly relevant when comparing standard TWFE models against heterogeneous-robust estimators like those proposed by callaway2021difference.
Major innovations in DiD methods have come from taking the potential outcomes framework seriously and determining what assumptions are necessary for a model to recover the relevant treatment effect. \citet*{Roth_2023} provides a recent survey of these developments. In general, DiD models are used to estimate the $\mbox{ATT}$ in \hyperref[{eq:defATT}]{\tagform@{\ref*{eq:defATT}}} as a summary measure of treatment effects.
As only one of the potential outcomes can be realized for any individual, care must be taken when constructing the counter-factual potential outcome for individuals. While there are many assumptions needed for DiD to recover the ATT, the two most critical are the parallel trends assumption and the no-anticipation assumption. The parallel trends assumption implies that the trends of untreated potential outcomes are parallel across treated and control regions roth_2023PTs, while the no-anticipation assumption implies that treatment has no effect on treated units prior to treatment Roth_2023. To fix ideas, consider a simple $2\times2$ setup where there are two regions, treatment and control, and two time periods 1 and 2. Treatment occurs in the treated region in period 2.
Following callaway2021difference, under the parallel trends and no-anticipation assumptions, the ATT in \hyperref[{eq:defATT}]{\tagform@{\ref*{eq:defATT}}} can be expressed in terms of observable outcomes, with the unit $i$ suppressed for concision, as
In this equation observations with $D_r=1$, or regions $r \in R_1$, are treated, while those with $D_r=0$, or regions $r \in R_0$, are controls. Here $Y_{r,1}$ are observations from the first time period, and $Y_{r,2}$ are observations from the second time period.
For a long time, researchers believed that, with these two primary assumptions satisfied, a simple DiD regression would yield the ATT. With just two regions (one treatment and one control) and at least two time periods, the simple DiD regression is often specified as
Here, ${\rm post}_t$ is an indicator variable set equal to 1 for observations in the post intervention period, $D_r$ was previously defined, and ${\rm did}_{r,t} = {\rm post}_t \times D_r$. The coefficient $\beta$ on ${\rm did}_{r,t}$ will recover the ATT in this simple setting.
In settings where there are multiple control regions and multiple treatment regions, it is natural to extend the simple model to accommodate the additional regions and time periods. Specifically, \hyperref[{eq:2by2}]{\tagform@{\ref*{eq:2by2}}} is modified by replacing ${\rm post}_t$ with a set of year fixed effects ($\delta_t$) and $D_r$ with a set of region fixed effects ($\lambda_r$). This yields the TWFE model,
The ${\rm did}_{r,t}$ variable is conceptually the same as before, but it has to be hand coded to equal 1 for treated time periods within treated regions. This requires accounting for the treatment timing in each region. Accordingly, ${\rm did}_{r,t}$ is set to 0 if $r\in R_0$, and set to 1 if $r \in R_1 $ and $t \geq S_r$. Until about five years ago, it was believed that $\beta$ from the TWFE model would consistently estimate an ATT; see the discussion in the popular textbook MHE_2008 as an example.
In the past five or so years, however, several extremely influential papers have shown that, among other things, the TWFE estimator is inconsistent in some cases de2020twott, goodman2021difference,callaway2021difference. This failure is particularly acute in staggered adoption settings (where different regions adopt policies at different times) when treatment effects are heterogeneous across groups. This research also demonstrates that, with the TWFE estimator, certain treated region $\times$ year cells, which should get a positive weight for the average treatment effect, are instead weighted negatively. The TWFE estimator can also include forbidden comparisons, where the estimated treatment effect is a function of sub-effects calculated based on contrasting newly treated units with previously treated ones.
The callaway2021difference (CSDID) procedure, which is implemented in the Stata and Python packages csdid and the R package did, applies a two-step approach to estimate an ATT. The first step estimates a series of $2 \times 2$ treatment effects, and the second step aggregates them. Specifically, the first step enumerates all the two-group, two-period differences that could estimate a treatment effect, excluding any forbidden comparisons. Each such comparison includes a treated and a control group for one period before and one period after treatment. The estimators focus attention on the year immediately before treatment for all comparisons as the pre-intervention period. The CSDID procedure then estimates a $2 \times 2$ ATT using only the observations for those four group $\times$ year cells.
In contrast to the TWFE estimator, the “groups” are defined by treatment timing, rather than by regions. These cohorts are defined by the time when treatment was first adopted, if ever. Using the values of $S_r$ and ${\bm G}_g$ defined in (ref), CSDID estimates time-by-cohort treatment effects.
The $2 \times 2$ ATT based on four group-year cells is referred to as an $\mbox{ATT}(g,t)$, where $g$ indexes the treatment timing cohort, and $t$ indexes the post-treatment year. The control group can vary from one $\mbox{ATT}(g,t)$ to the next. It also depends on the choice of never-treated or not-yet-treated comparison units. For a particular $\mbox{ATT}(g,t)$, we can denote the control group as $\mathcal{G}_{{\rm comp}(g,t)}$. It is sometimes useful to differentiate between the never-treated control group $\mathcal{G}_{{\rm comp}(g,t),{\rm NT}}$ and the not-yet-treated control group $\mathcal{G}_{{\rm comp}(g,t),{\rm NY}}$. In both cases, the time periods used for the control group are the same as for the treated group, specifically, $t$ and $g-1$. Recall that $g-1$ is the last pre-adoption period for cohort $g$.
The never-treated control group $\mathcal{G}_{{\rm comp}(g,t),{\rm NT}}$ is comprised of the cohort ${\bm G}_\infty$, that is, the regions with $S_r=0$. In contrast, the not-yet-treated control group $\mathcal{G}_{{\rm comp}(g,t), {\rm NY}}$ is comprised of the ${\bm G}_\infty$ cohort and cohorts with $g > t$, that is, regions that were treated after period $t$, for which $S_r > t$. Sometimes, all groups are eventually treated. In that case, $R_0$ is empty, and only the not-yet-treated control group is used; see (ref) below. For convenience, the comparison group is simply referred to as $\mathcal{G}_{\rm comp}$ in what follows.
Without covariates, a definition of the $\mbox{ATT}(g,t)$ is
The argument $t$ is the specific post-treatment time period for which the $\mbox{ATT}(g,t)$ is calculated. For each $g \in \mathcal{G}$, there is an $\mbox{ATT}(g,t)$ for all $g \leq t \leq T$. Thus there is an $\mbox{ATT}(g,t)$ for each timing cohort, for all post-treatment periods in which there exists an untreated comparison group. The $\mbox{ATT}(g,t)$ can be estimated in various ways. If treatment is exogenous and we have a fully-balanced panel, then a straightforward estimator, analogous to \hyperref[{eqn: callaway estimand}]{\tagform@{\ref*{eqn: callaway estimand}}}, is
where $g$ indexes the group or cohort for which the $\mbox{ATT}(g,t)$ is estimated. This is the ${\bm G}_g$ cohort, or the cohort which is first treated at time $g$. Recall that this cohort is comprised of the regions with $S_r=g$. Here, $Y_{i,r,t}$ is the post-treatment period outcome for person $i$ in region $r$, and $Y_{i,r,g-1}$ is the pre-treatment period outcome. The two double summations in \hyperref[{eqn: callaway sample analog}]{\tagform@{\ref*{eqn: callaway sample analog}}} indicate that we are summing over all observations $i$ from all regions $r$ that are part of the specified cohort, either ${\bm G}_g$ or $\mathcal{G}_{\rm comp}$. If $N_{r,t}$ is the number of observations in region $r$ in time period $t$, then $N_{g,t}$ is the total number of observations in cohort $g$ in time period $t$, so that $N_{g,t} = \sum_{r : S_r = g}N_{r,t}$. Similarly, $N_{\mathcal{G}_{{\rm comp},t}}$ is the total number of observations in time period $t$ for the comparison group, so that $N_{\mathcal{G}_{{\rm comp},t}} = \sum_{r \in \mathcal{G}_{\rm comp}}N_{r,t}$.
With a fully-balanced panel, every unit is observed in every period, and there are no missing observations. Thus $N_{g,t}$ is also the total number of observations in cohort $g$ in time period $g-1$. In repeated cross-sectional datasets, or in panel datasets with attrition, the number of observations in each cohort can differ across time periods, and \hyperref[{eqn: callaway sample analog}]{\tagform@{\ref*{eqn: callaway sample analog}}} would have to be replaced by a more complicated equation.
The second step involves taking a weighted average of the $\widehat\mbox{ATT}(g,t)$ to obtain an overall estimate of the ATT,
The weights $w_{g,t}$ sum to one and are chosen by the user, but the default is equal weighting (referred to as “simple” weighting), so that $\widehat\mbox{ATT}$ is simply the average of the $\widehat\mbox{ATT}(g,t)$.
In the remainder of this subsection, we consider three examples which illustrate how the CSDID procedure works for staggered adoption.
Clustering at the region level causes difficulties for inference. The estimated ATT is normally a weighted sum of several $\widehat\mbox{ATT}(g,t)$ terms; see, for example, \hyperref[{eq:cssampleagg}]{\tagform@{\ref*{eq:cssampleagg}}}. If we want to cluster at the region level, we have to take into account correlations of the scores across all observations that are used to estimate the ATT. These correlations arise because observations from a single region may be used to estimate multiple $\mbox{ATT}(g,t)$ terms. For instance, in (ref), A,1 and C,1 are used in the estimation of both $\mbox{ATT}(2,2)$ and $\mbox{ATT}(2,3)$, and C,3 is used in the estimation of both $\mbox{ATT}(2,3)$ and $\mbox{ATT}(3,3)$. This implies that data from region C are used to estimate all three $\mbox{ATT}(g,t)$ terms, and that data from region A are used to estimate both $\mbox{ATT}(2,2)$ and $\mbox{ATT}(2,3)$. A valid routine for cluster-robust inference has to account for the resulting correlations.
Difficulties with inference when the observations are clustered is not unique to the procedures of callaway2021difference, as two recent papers discuss. mizushima2025inference investigates the size properties of several modern DiD estimators. It finds that the default, or asymptotic, procedures are often over-sized when there are few treated clusters. However, the wild cluster bootstrap and a randomization inference procedure can both greatly improve test size for imputation-type tests. weiss2024much considers the inference properties of many DiD estimators and finds that all of them can lack power, notably CSDID, TWFE, imputation \citep*{borusyak2024revisiting}, and DCDH de2020twott. In fact, power can even decline when additional time periods are included in the analysis.
Neither of these papers considers the jackknife for improving inference, as we do in (ref), although weiss2024much does consider the so-called $\mbox{CV}_{\kern -0.08333em2}$ variance estimator; see MNW-guide.
The default method for obtaining standard errors in callaway2021difference uses the contributions to the influence functions by each of the observations. The expressions for the influence functions are quite long, so we omit them; see Theorem 2 of the paper. Here we are concerned not with how the RIF variance estimator is derived but with how well it performs in practice.
A simplified version of the second result in callaway2021difference can be written as
Here, following the notation in santanna2023did, $\textbf{ATT}_{g \leq t}$ and $\widehat{\textbf{ATT}}_{g \leq t}$ denote the vectors of $\mbox{ATT}(g, t)$ and $\widehat\mbox{ATT}(g, t)$, respectively, for all $g = 2, \ldots, t$ and $t = 2, \ldots, T$. The variance matrix ${\bm{\Sigma}}$ is defined by
where ${\bm{W}}_i$ is the data vector for observation $i$, and ${\bm{\Psi}}_{g\leq t}({\bm{W}}_i)$ are the influence functions. The result \hyperref[{ATTvar}]{\tagform@{\ref*{ATTvar}}} is used to compute asymptotic RIF standard errors based on $\hat{\bm{\Sigma}}$, the sample analog of ${\bm{\Sigma}}$. Note that constructing $\hat{\bm{\Sigma}}$ can be a bit tricky.
These RIF standard errors are then used with the standard normal distribution to calculate $P$ values and confidence intervals. This makes sense when there is no clustering, which is assumed, as the presence of the factor $\sqrt{n}$ in \hyperref[{ATTvar}]{\tagform@{\ref*{ATTvar}}} makes evident. When the observations are clustered, however, it seems like an odd thing to do. Almost all empirical work that uses cluster-robust $t$-statistics for linear regression models assumes that these are distributed as $t(H-1)$, where $H$ is the number of clusters. As \citet*{MNW-logit} illustrates, using the ${\rm N}(0,1)$ distribution instead of the $t(H-1)$ distribution increases rejection frequencies substantially when $H$ is small. Thus, it seems highly unlikely that the default RIF method for inference will yield reliable results when there is clustering and $H$ is not large. Indeed, callaway2021difference points out that their approach requires a large number of clusters.
As an alternative way to obtain standard errors, which may or may not be cluster-robust, callaway2021difference suggests using a form of multiplier bootstrap. For each of $B$ bootstrap samples, indexed by $b$, a vector of bootstrap estimates is constructed as
Here $\hat{\bm{\Psi}}_{g\leq t}({\bm{W}}_i)$ is the sample analog of the influence function ${\bm{\Psi}}_{g \leq t}({\bm{W}}_i)$. The $V_{b,i}^*$ are random draws from an auxiliary distribution with mean 0 and variance 1, typically the Mammen distribution. When the data are clustered, these draws are constant within each cluster, as for the wild cluster bootstrap; see santanna2023did for additional details. The bootstrap estimates are centered around the original estimate, with perturbations of the influence functions at the observation or cluster level added to them. The $B$ vectors of bootstrap $\widehat\mbox{ATT}^*_b(g,t)$ values given by \hyperref[{ATTboot}]{\tagform@{\ref*{ATTboot}}} are then used to form confidence intervals for the individual $\mbox{ATT}(g,t)$ or for the ATT.
This multiplier bootstrap is conceptually simple and should be computationally efficient if it is implemented appropriately. Using it avoids the need to construct $\hat{\bm{\Sigma}}$ analytically. It differs from the widely used wild cluster bootstrap \citep*{CGM_2008,DMN_2019} in one important way. With the wild cluster bootstrap, each bootstrap replication involves re-estimating the coefficients and test statistics of interest. With the multiplier bootstrap, however, these are not re-estimated. The bootstrap draws simply add noise to the original estimate, and the $\widehat\mbox{ATT}^*_b(g,t)$ that result are then used to estimate how much the $\widehat\mbox{ATT}(g,t)$ vary. Since the model is never re-estimated, this is just a convenient way to estimate ${\bm{\Sigma}}$. Thus it is not surprising that the multiplier bootstrap and asymptotic approaches to inference yield very similar answers; see (ref).
As an alternative to the methods suggested in callaway2021difference, we propose to use the cluster jackknife to calculate standard errors for $\widehat\mbox{ATT}$. The cluster jackknife has been shown to perform well in many cases, especially in settings with few and unbalanced clusters \citep*{MNW-bootknife,Hansen_2025}. Another attractive feature of the jackknife is that no analytical expressions for the influence functions are needed. Instead, the cluster jackknife simply omits each cluster $h$ in turn, estimates each of the $\mbox{ATT}(g,t)$, and then uses those estimates to re-estimate the ATT. The standard error is then simply the square root of
where $\widehat\mbox{ATT}$ is the original CSDID estimate of the ATT, and $\widehat\mbox{ATT}^{(h)}$ is the estimate when cluster $h$ is omitted. Here, following much of the literature, we use CV$_{\kern -0.08333em3}$ to denote the cluster-jackknife variance estimator, because it generalizes the jackknife HC$_3$ variance estimator of MW_1985.
Although the cluster jackknife has proven to work well in many settings, a few caveats are in order. The first is that the jackknife is not amenable to all staggered adoption settings. weiss2024much discusses this issue and does not investigate the jackknife for that reason. For instance, consider (ref) with clustering at the region level. In this case, some of the $\mbox{ATT}(g,t)$ terms cannot be calculated when certain regions are dropped. The $\mbox{ATT}(2,2)$ and $\mbox{ATT}(2,3)$ terms cannot be calculated when region A is dropped, the $\mbox{ATT}(3,3)$ term cannot be calculated when region B is dropped, and no ATT terms can be calculated when region C is dropped. In contrast, all of the $\mbox{ATT}(g,t)$ terms can be estimated for each jackknife sample for a setup like the one in (ref). For that reason, we designed the Monte Carlo simulations described in (ref) so that there are always multiple treated clusters and untreated clusters. However, in simulations with one early adopter and one late adopter, there are not multiple regions in each cohort.
The second caveat is that some staggered adoption settings will alter the weights given to different $\mbox{ATT}(g,t)$ across jackknife replications. Consider a modification of (ref), where there are now two control groups, C and D. In this case, it would be possible to estimate an ATT for each jackknife replication. Dropping either region C or D would still allow all three $\mbox{ATT}(g,t)$ terms to be calculated. However, dropping region A would mean that the ATT would be based solely on $\mbox{ATT}(3,3)$. Likewise, dropping region B would mean that the ATT would just be a weighted average of $\mbox{ATT}(2,2)$ and $\mbox{ATT}(2,3)$. If the ATT is sensitive to the regions used to compute it, this variation in weights may highlight that fact and thus provide valuable information. However, the $\widehat\mbox{ATT}^{(h)}$ may now vary across jackknife samples in ways that would not happen for other applications of the jackknife. Thus, unlike with other inference procedures, care must be taken to assess whether the jackknife is a sensible tool to use for a given treatment design.
The third caveat is that the jackknife procedure can be slow. Many of the modern DiD routines, including CSDID, can be slow to calculate when samples are large and there are many $\mbox{ATT}(g,t)$ terms to estimate. The jackknife we propose needs to call whatever command implements CSDID $H+1$ times, once for the full sample estimate $\widehat\mbox{ATT}$, and $H$ additional times for each of the $\widehat\mbox{ATT}^{(h)}$ estimates. Fortunately, the results of our simulations in (ref) suggest that this increase in computational time can result in large reductions in over-rejection by hypothesis tests or under-coverage by confidence intervals.
Our Monte Carlo simulations use Current Population Survey (CPS) data. The objective is to compare the rejection frequencies for several different tests of $\mbox{ATT}=0$. We conduct a set of “placebo laws” experiments inspired by \citet*{BDM_2004}, in which we pretend that treatment occurred in various states at various times. Independently, mizushima2025inference ran very similar experiments, comparing several DiD estimators. Their simulations examine state-by-year unemployment rates, while ours investigate earnings at the individual level. Our simulations were conducted in Stata, using both the csdid package and our new csdidjack package. In our experiments, we compare the rejection frequencies using the default cluster option and the bootstrap option wboot within csdid with those using jackknife standard errors from our new post-estimation package \texttt{csdidjack} discussed in (ref).
We employ the same dataset used for the placebo-laws simulations in MW-JAE. Our dataset is extracted from the Merged Outgoing Rotation Group (MORG) of the CPS. It includes individual-level information on weekly earnings, state, and year, covering the period from 1979 to 1999. The dependent variable is the log of women's earnings. We exclude observations with earnings less than \$20, which may be erroneous. The total sample, from all years and states, has 547,818 observations. Across replications, we randomly pick 8 consecutive years of data and $R$ states from the full sample. To mimic staggered adoption, treatment is assigned in two cohorts, an early cohort of $J$ states treated in year 4, and a late cohort of $L$ states treated in year 6. This results in $J+L$ total treated states, and we set $J=L$.
Note that treatments apply to all individuals in a state. Treatment is also “absorptive,” so that once a state is treated in a particular replication, it is always treated. For the TWFE estimator, there can be severe inferential problems when there are few clusters or few treated clusters. To see whether similar problems arise for the CSDID procedure, we consider different values for both the number of states $R$ and the numbers of early and late adopters. (ref) shows the various values of $R$, $J$, and $L$ that we consider.
In each simulation, after obtaining a subsample, we estimate the ATT using the csdid package, both with and without covariates. In all simulations, we use simple aggregation (i.e. unweighted averaging) to estimate the ATT, and we test the hypothesis that $\mbox{ATT}=0$. The covariates are age, age squared, and a set of four education dummies. These covariates, and the outcome variable, have all been demeaned at the state level using the full set of data. This demeaning should strengthen the null of no treatment effects.
We calculate three rejection frequencies at the 5% level, always clustering at the state level. The first is for $t$-statistics based on the “cluster” option of csdid, which uses the recentered influence functions discussed in (ref). The second is for the multiplier bootstrap discussed in (ref). Since the bootstrap procedure in csdid does not provide a $P$ value, we calculate the rejection rate based on whether the 95% confidence interval covers the true value of 0. The final rejection frequency is for $t$-statistics based on cluster-jackknife standard errors discussed in (ref) and calculated using our csdidjack package. Since all three inference methods use states as clusters, $H=R$ in all our experiments.
Our experiments have 2,400 replications. This number is quite a bit smaller than we would like, but the simulations are very slow, particularly for the jackknife. We hope to speed them up in the future using a methodology similar to the one in the fastdid package in R. This package is up to 100x faster than the did package in R, but it requires panel data rather than cross-sectional data. It can be found at \url{https://cran.r-project.org/web/packages/fastdid/index.html}.
(ref) shows rejection frequencies for the Monte Carlo simulations with covariates. There are three panels, one for each of the inference methods that we study. The results for the two procedures proposed in callaway2021difference in Panels A and B generally resemble each other greatly, while those for our cluster-jackknife procedure in Panel C are always much closer to the nominal level of 0.05.
For the asymptotic and multiplier bootstrap methods, there are clearly severe problems with inference when there are few clusters. The ideal case for each value of $H$ is when half the clusters are treated; see MW-EJ and many other papers. This case occurs when $H=2(J+L)=4J$. Except with $H=32$, the methods in Panels A and B work best for that case. For $H=32$, they work slightly better when $J=10$ than when $J=8$, but that may be due to the small number of replications. It is clear that both methods improve as $H$ increases, but they never perform very well. In contrast, the jackknife, in Panel C, always rejects around 5% of the time for several of the better cases.
For every method, there are serious inferential problems when there are few treated clusters. The top row of each panel, where $J=L=1$, always contains the largest rejection frequencies. The asymptotic and bootstrap procedures reject between 30% and 40% of the time. In contrast, the jackknife procedure never rejects much more than 16% of the time. Increasing the number of early and late adopters reduces rejection rates substantially, at least up to the point where $H=4J$.
The jackknife is always reliable when $J=L$ is large enough. Interestingly, the number of treated clusters that are needed for the jackknife to be reliable depends on the total number of clusters. Similar results were found for TWFE with CV$_1$ inference in MW-JAE and with jackknife inference in \citet*{MNW-bootknife}.
The similarity of the asymptotic and bootstrap rejection frequencies is immediately evident from the table. At most, they differ by about 2% (0.3333 vs 0.3458), but by much smaller amounts when the magnitude of the over-rejection is smaller. It is clear that the multiplier bootstrap does not improve inference much at all, as the discussion in (ref) suggested. For that reason, we will not focus on the multiplier bootstrap in the empirical examples.
(ref) shows rejection frequencies for the regressions without covariates. These results are broadly similar to the results with covariates. Again, the asymptotic and bootstrap tests are always very similar and always over-reject, in some cases very severely. The jackknife tests generally provide much more reliable inferences, but they can still over-reject substantially when there are few treated clusters.
Overall, there are four major takeaways from these experiments:
We now present two brief empirical examples to demonstrate that there can be meaningful differences between inferences based on asymptotic and cluster-jackknife standard errors. Both examples use American data clustered by state. The first example uses individual-level data, and the second uses county-level data. Both population and the number of counties vary considerably across states.
We first consider an empirical example from \citet*{flex_2026}, which uses the variation explored in \citet*{meinhofer2022prenatal}. The latter paper examines whether the enactment of punitive prenatal substance-use policies affects outcomes of children, while the former considers whether these policies had an impact on the mental health of mothers of child-bearing age. We thank the authors of flex_2026 for making their replication files easily available.
We use data from the Behavioral Risk Factor Surveillance System (BRFSS) for 2005--2018, a large annual survey conducted by the Centers for Disease Control and Prevention. We restrict attention to women aged 18--44 with at least one child, resulting in a sample of 440,446 observations from 34 U.S.\ states. Our outcome of interest is the number of days in the past month that a respondent reports having good mental health. This variable is highly left-skewed, with a sample mean of approximately 25.5 days.
Thirteen states adopted punitive prenatal substance-use policies during the sample period. The timing of adoption varies considerably, with Idaho adopting in 2007, South Carolina in 2008, Arizona in 2009, and a number of states enacting policies in 2012--2018. The never-treated group consists of 21 states that never adopted such policies during the sample period.
The analysis includes a limited set of covariates: age, race (white, black, hispanic), education (no degree, high school, some college, college degree), and the log of income. The mean age of respondents is 34 years.
We estimate the ATT for the CSDID estimator with simple aggregation and clustered at the state level, using the csdid package. As shown in (ref), this standard approach yields a statistically significant estimate at the 5% level. However, the cluster jackknife (using csdidjack) leads to a larger standard error and a $P$ value greater than 0.05. The sample is relatively unbalanced across the 34 states, with the largest states having nearly 5 times as many observations as the smallest ones. This is a setting in which we might expect the jackknife to be more reliable \citep*[Section 6]{MNW-bootknife}.
We also consider an empirical example adapted from callaway2021difference. The analysis examines the impact of minimum-wage increases on teen employment using county-level data from the Quarterly Workforce Indicators for 2001--2007. During this time, the federal minimum wage remained flat, allowing us to define treatment based on when a state first increased its minimum wage. We focus on 29 states that changed their minimum wage during the sample period. The outcome is the log of teen employment (lemp), and the covariates include demographic and economic controls: region indicators, percent white, percent with a high school diploma, poverty rate, and both county population and median income as quadratics.
Estimation is performed using the csdid command with simple aggregation, clustered at the state level. Results are in Panel A of (ref). The ATT is $-0.0490$, with an asymptotic standard error of $0.0171$ ($P = 0.004$), suggesting a statistically very significant reduction in teen employment. In contrast, the cluster-jackknife standard error is $0.1960$ ($P = 0.8045$).
Given the striking difference between the asymptotic and jackknife standard errors, we investigate the influence of each state on the aggregate effect. In (ref), we plot these estimates against the number of observations in the state/cluster. The red dots indicate control states, and the blue dots indicate treated states. It is evident from this figure that one control state is an extreme outlier. The outlier state is ID 6 (Indiana). Although Indiana has an intermediate number of observations, dropping it results in a subsample ATT estimate of $0.1476$, whereas all other leave-one-out ATT estimates are negative and close in magnitude to the full sample estimate.
The extreme influence of Indiana suggests that it is a high-leverage outlier. We therefore drop Indiana and repeat the analysis with the remaining 28 states. These results are presented in Panel B of (ref). With Indiana excluded, the ATT is $0.1476$. The asymptotic standard error increases to 0.1385, and the jackknife standard error increases to $0.5159$; in both cases, the coefficient is no longer statistically significant. Evidently, as is the case when using the jackknife for cluster-robust inference in least squares regression, it is highly instructive to inspect the full vector of omit-one-cluster estimates to ensure that results are not driven by a single cluster. Our summclust package \citep*{MNW-influence} provides these estimates for TWFE and other least squares regression models. Our csdidjack package (see below) makes them available for CSDID estimation of the ATT.
We provide open-source software implementations of the cluster jackknife (CV$_3$) inference procedure for both Stata and R users.
The csdidjack package is a post-estimation command for use with the csdid command \citep*{rios2021csdid} and with the built-in Stata command hdidregress. The package supports the agg(simple), agg(group), and \texttt{agg(calendar)} aggregation schemes.
To install or update csdidjack in Stata, run:
Once a model has been estimated using csdid or hdidregress and estat aggregation, users can run:
This command returns CV$_3$ standard errors, $t$-statistics, $P$ values, and confidence intervals. For further documentation, type help csdidjack in Stata. Source code is available at \url{https://github.com/liu-yunhan/csdidjack}.
The didjack package provides CV$_3$ jackknife inference for users of the did package by Callaway and Sant’Anna in R.
To install the package from GitHub, run:
An example of its usage is:
Full documentation and source code are available at \url{https://github.com/liu-yunhan/didjack}.
Obtaining reliable cluster-robust inferences for DiD estimates of treatment effects has been the subject of much research for over two decades. A great deal of progress has been made, but there are still cases (such as few clusters and/or few treated clusters) where inference can be quite unreliable; see \citet*{MNW-guide}.
In recent years, well-founded concerns over staggered adoption have led to new ways to estimate treatment effects using the so-called “modern DiD” methods. Unfortunately, for these new estimators, the concern over obtaining unbiased estimates of the ATT has shifted focus away from the question of statistical inference. This short paper highlights the fact that the long-standing problems of few clusters and few treated clusters for TWFE estimation can be severe with the CSDID estimator of callaway2021difference. We propose to use cluster-jackknife standard errors instead of the ones suggested in that paper. Simulations strongly suggest that this can lead to much more reliable inferences.
\addcontentsline{toc}{section}{\refname}