*! version 3.1 2023-11-22

Title

did_imputation - Treatment effect estimation and pre-trend testing in event studies: difference-in-differences designs with staggered adoption of treatment, using the imputation approach of Borusyak, Jaravel, and Spiess (April 2023)

Syntax

did_imputation Y i t Ei [if] [in] [estimation weights] [, options]

Y outcome variable

i variable for unique unit id

t variable for calendar period

Ei variable for unit-specific date of treatment (missing = never-treated)

Note: These main parameters imply:

- the treatment indicator: D=1[t>=Ei];

- "relative time", i.e. the number of periods since treatment: K=(t-Ei) (possibly adjusted by the shift and delta options described below).

Note: did_imputation requires a recent version of reghdfe. If you get error messages (e.g. r(123) or "verbose must be between 0 and 5"), please (re)install reghdfe to make sure you have the most recent version.

Note: Before emailing the authors about errors, please read the Bug reporting section of this helpfile.

Description

did_imputation estimates the effects of a binary treatment with staggered rollout allowing for arbitrary heterogeneity and dynamics of causal effects, using the imputation estimator of Borusyak et al. (2023).

The benchmark case is with panel data, in which each unit i that gets treated as of period Ei stays treated forever;
some units may never be treated. Other types of data (e.g. repeated cross-sections) and other designs (e.g. triple-diffs) are also allowed;
see Usage examples.

Estimation proceeds in three steps:

1. Estimate a model for non-treated potential outcomes using the non-treated (i.e. never-treated or not-yet-treated)
observations only. The benchmark model for diff-in-diff designs is a two-way fixed effect (FE) model: Y_it = a_i + b_t + eps_it,
but other FEs, controls, etc., are also allowed.

2. Extrapolate the model from Step 1 to treated observations, imputing non-treated potential outcomes Y_it(0),
and obtain an estimate of the treatment effect tau_it = Y_it - Y_it(0) for each treated observation. (See What if imputation is not possible)

3. Take averages of estimated treatment effects corresponding to the estimand of interest.

A pre-trend test (for the assumptions of parallel trends and no anticipation) is a separate exercise
(see the pretrends option). Regardless of whether the pre-trend test is performed, the reference group
for estimation is always all pre-treatment (or never-treated) observations.

To make "event study" plots, please use the accompanying command event_plot.

What if imputation is not possible

The imputation step (Step 2) is not always possible for all treated observations:

- With unit FEs, imputation is not possible for units treated in all periods in the sample;

- With period FEs, it is impossible to isolate the period FE from the variation in treatment effects
in a period when all units have already been treated (and if there are never-treated units);

- If you include group#period FEs, imputation is further impossible once all units in the group have been treated;

- Similar issues arise with other covariates in the model of Y(0).

This is a fundamental issue: the model you specified does not allow to find unbiased estimates of treatment effects
for those observations (without restrictions on treatment effects; see Borusyak et al. 2023).

If this problem arises (i.e. there is at least one treated observation that enters one of the estimands of interest with a non-zero weight
and for which the treatment effect cannot be imputed),
the command will throw an error and generate a dummy variable cannot_impute which equals to one for those observations.

You have two ways to proceed:

- Modify the estimand, excluding those observations manually: via the if clause or by setting the weights on them to zero
(via the wtr option);

- Specify the autosample option that will do this automatically in most cases. But we recommend that you still review cannot_impute
to understand what estimand you will be getting.

Options

Model of Y(0)

fe(list of FE): which FE to include in the model of Y(0). Default is fe(i t) for the diff-in-diff (two-way FE) model.
But you can include fewer FEs, e.g. just period FE fe(t) or, for repeated cross-sections at the individual level, fe(state t).
Or you can have more FEs: e.g. fe(i t#state) (for state-by-year FE with county-level data)
or fe(i#dow t) (for unit by day-of-week FE). Each member of the list has to look like
v1#v2#...#vk. If you want no FE at all, specify fe(.).

controls(varlist): list of continuous time-varying controls. (For dummy-variable controls, e.g. gender, please use the
fe option for better convergence.)

unitcontrols(varlist): list of continuous controls (often unit-invariant) to be included interacted with unit dummies.
E.g. with unitcontrols(year) the regression includes unit-specific trends.
(For binary controls interacted with unit dummies, use the fe option.)

Use with caution: the command may not recognize that imputation is not possible for some treated observations.
For example, a unit-specific trend is not possible to estimate if only one pre-treatment observation is available for the unit, but it is not
guaranteed that the command will throw an error.

timecontrols(varlist) list of continuous controls (often time-invariant) to be included interacted with period dummies.
E.g. with timecontrols(population) the regression includes i.year#c.population.
(For binary controls interacted with period dummies, use the fe option.)

Use with caution: the command may not recognize that imputation is not possible for some treated observations.

Estimands and Pre-trends

wtr(varlist): A list of variables, manually defining estimands of interest by storing the weights on the treated observations.

- If estimation weights (aw/iw/fw) are used, wtr weights will be applied in addition
to those weights in defining the estimand.

- If nothing is specified, the default is the simple ATT across all treated observations (or, with horizons or allhorizons, by horizon). So wtr=1/number of the relevant observations.

- Values of wtr for untreated observations are ignored (except as initial values in the iterative procedure for computing SE).

- Values below 0 are only allowed if sum is also specified (i.e. for weighted sums and not weighted averages).

- Using multiple wtr variables is faster than running did_imputation for each of them separately, and produces a joint variance-covariance matrix.

sum: if specified, the weighted sum, rather than average, of treatment effects is computed (overall or by horizons).
With sum specified, it's OK to have some wtr<0 or even adding up to zero;
this is useful, for example, to estimate the difference between two weighted averages of treatment effects
(e.g. across horizons or between men and woment).

horizons(numlist): if specified, weighted averages/sums of treatment effects will be reported for each of these horizons separately
(i.e. tau0 for the treatment period, tau1 for one period after treatment, etc.). Horizons which are not specified will be ignored.
Each horizon must a be non-negative integer.

allhorizons: picks all non-negative horizons available in the sample.

hbalance: if specified together with a list of horizons, estimands for each of the horizons will be based
only on the subset of units for which observations for all chosen horizons are available
(note that by contruction this means that the estimands will be based on different periods).
If wtr or estimation weights are specified, the researcher needs to make sure that the weights are constant over time
for the relevant units---otherwise proper balancing is impossible and an error will be thrown.
Note that excluded units will still be used in Step 1 (e.g. to recover the period FEs) and for pre-trend tests.

hetby(varname): reports estimands separately by subgroups defined by the discrete (non-negative integer or string)
variable provided. This is the preferred option for treatment effect heterogeneity analyses (but see also project).

project(varlist): projects (i.e., regresses) treatment effect estimates on a set of numeric variables
and reports the constant and slope coefficients. The variables should not be collinear. (To analyse effect heterogeneity
by subgroup, option hetby is preferred. Note that standard errors may not agree exactly between hetby and
project.)

minn(#): the minimum effective number (i.e. inverse Herfindahl index) of treated observations,
below which a coefficient is suppressed and a warning is issued.
The inference on coefficients which are based on a small number of observations is unreliable. The default is minn(30).
Set to minn(0) to report all coefficients nevertheless.

autosample: if specified, the observations for which FE cannot be imputed will be automatically dropped from the sample,
with a warning issued. Otherwise an error will be thrown if any such observations are found.
autosample cannot be combined with sum or hbalance;
please specify the sample explicitly if using one of those options and you get an error that imputation has failed (see What if imputation is not possible).

shift(integer): specify to allow for anticipation effects.
The command will pretend that treatment happened shift periods earlier for each treated unit.

- Do NOT use this option for pre-trend testing;
use it if anticipation effects are expected in your setting.
(This option can be used for a placebo test but we recommend a pretrend test instead; see Section 4.4 of Borusyak et al. 2023.)

- The command's output will be labeled relative to the shifted treated date Ei-shift.
For example, with horizons(0/10) shift(3) you will get coeffieints _b[tau0]..._b[tau10] where tauh
is the effect h periods after the shifted treatment. That is, tau1 corresponds to the average anticipation effect 2 periods before
the actual treatment, while tau8 to the average effect 5 periods after the actual treatment.

pretrends(integer): if some value k>0 is specified, the command will performs a test for parallel trends,
by a separate regression on nontreated observations only: of the outcome on the dummies for 1,...,k periods before treatment,
in addition to all the FE and controls. The coefficients are reported as pre1,...,prek.
The F-statistic (corresponding to the cluster-robust Wald test), corresponding pvalue, and degrees-of-freedom are reported in e(pre_F),
e(pre_p), and e(pre_df) resp.

- Use a reasonable number of pre-trends, do not use all of the available ones unless you have a really large never-treated group. With too many pre-trend coefficients, the power of the joint test will be lower.

- The entire sample of nontreated observations is always used for pre-trend tests, regardless of hbalance and other options that restrict the sample for post-treatment effect estimation.

- The number of pretrend coefficients does not affect the post-treatment effect estimates, which are always computed under the assumption of parallel trends and no anticipation.

- The reference group for the pretrend test is all periods more than k periods prior to the event date (and all never-treated
observations, if available).

- Because of this reference group, it is expected that the SE are the largest for pre1 (opposite from some conventional tests).

- This is only one of many tests for the parallel trends and no anticipation assumptions. Others are easy to implement manually; please see
the paper for the discussion.

Standard errors

cluster(varname): cluster SE within groups defined by this variable. Default is i.

avgeffectsby(varlist): Use this option (and/or leaveout) if you have small cohorts of treated observations, and after reviewing
Section 4.3 of Borusyak et al. (2023). In brief, SE computation requires averaging the treatment effects by groups of treated observations.

- These groups should be large enough, so that there is no downward bias from overfitting.

- But the larger they are, the more conservative SE will be, unless treatment effects are homogeneous within these groups.

- The varlist in avgeffectsby defines these groups.

- The default is cohort-years avgeffectsby(Ei t), which is appropriate for large cohorts.

- With small cohorts, specify coarser groupings: e.g. avgeffectsby(K) (to pool across cohorts) or
avgeffectsby(D) (to pool across cohorts and periods when computing the overall ATT).

- The averages are computed using the "smart" formula from Section 4.3, adjusted for any clustering and any choice of avgeffectsby.

leaveout: Recommended option. In particular, use it (and/or avgeffectsby) if you have small cohorts of treated observations.
The averages of treatment effects will be computed excluding the own unit (or, more generally, cluster).
See Section 4.3 of Borusyak et al. (2023) for details.

alpha(real): confidence intervals will be displayed corresponding to that significance level. Default is 0.05.

nose: do not produce standard errors (much faster).

Miscellaneous

saveestimates(newvarname): if specified, a new variable will be created storing the estimates of treatment effects for each observation
are stored. The researcher can then construct weighted averages of their interest manually (but without SE).

Note: Individual estimates are of course not consistent. But weighted sums of many of them typically are, which is what did_imputation generally reports.

saveweights: if specified, new variables __w_* are generated, storing the weights corresponding to (each) coefficient.
Recall that the imputation estimator is a linear estimator that can be represented as a weighted sum of the outcomes.

- These weights are applied on top of any estimation weights.

- For treated observations these weights equal to the corresponding wtr - that's why the estimator is unbiased
under arbitrary treatment effect heterogeneity.

- If a weighted average is estimated (i.e. sum is not specified) and there are no estimation weights,
the weights add up to one across all treated observations.

- With unit and period FEs, weights add up to zero for every unit and time period (when weighted by estimation weights).

loadweights(varlist): use this to speed up the analysis of different outcome variables with an identical specification
on an identical sample. To do so, provide the set of the weight variables (__w*, but can be renamed),
saved using the saveweights option when running the analysis for the first outcome.
[Warning: the validity of the weights is assumed and not double-checked.] [Currently works only if the varnames are not __w*]

saveresid(name): if specified, a new variable name_* is generated for each estimand (e.g. name_tau0) to store
model residuals used in the computation of standard errors. For untreated observations, they are residuals from the estimation step. For
treated observations, they are the epsilon-tildes from the BJS Theorem 3 (based on the avgeffectsby option) or the leave-out
versions from Appendix A.6. The residuals may be heterogeneous across estimands in general --- see equation (8) in the paper
which depends on the estimator weights, and thus on the estimand. This option can be helpful for reproducing the standard errors manually.

delta(integer): indicates that one period should correpond to delta steps of t and Ei.
Default is 1, except when the time dimension of the data is set (via tsset or xtset);
in that case the default is the corresponding delta.

tol(real), maxit(integer): tolerance and the maximum number of iterations.
This affects the iterative procedure used to search for the weights underlying the estimator (to produce SE).
Defaults are 10^-6 and 100, resp. If convergence is not achieved otherwise, try increasing them.

verbose: specify for the debugging mode.

Estimation weights

Estimation weights (only aw or iw are allowed) play two roles:

1. Step 1 estimation is done by a weighted regression. This is most efficient when the variance
of the error terms is inversely proportionate to the weights.

2. In Step 3, the average that defines the estimand is weighted by these weights.
If wtr is specified explicitly, estimation weights are applied on top of them.

- If sum is not specified, the estimand is the average of treatment effects weighted by wtr*aw

- If sum is specified, the estimand is the sum of treatment effects multiplied by wtr*aw

Do NOT specify regression weights if you only have the second motivation to use the weights,
e.g. if you want to measure the ATT weighted by county size but you have no reason to think that the outcomes of larger counties have less noise.
Instead specify the estimand of your interest via wtr.

All weight types (aw/iw/fw) produce identical results. Weights should always be non-negative.

Stored results

did_imputation stores the following in e():

Matrices

e(b) A row-vector of (i) the estimates, (ii) pre-trend coefficients, and (iii) coefficients on controls from Step 1:

- If horizons is specified, the program returns tauh for each h in the list of horizons.

- If multiple wtr are specified, the program returns tau_v for each v in the list of wtr variables.

- Otherwise the single returned coefficient is called tau.

- If hetby is specified, the coefficient names are appended with underscore and the values of the grouping variable. For instance, with hetby(female) where female takes values 0 and 1, the command will return tauh_0 and tauh_1 for each horizon, corresponding to average treatment effects by horizon and sex.

- If project is specified, the coefficient names are appended with underscore and the constant plus the list of slope coefficients. For instance, with project(female income), the command will return tauh_cons, tauh_female, and tauh_income.

- In addition, if pretrends is specified, the command returns preh for each pre-trend coefficient h=1..pretrends.

- And if controls is specified, the command returns the coefficients on those controls. (Estimated fixed effects, unitcontrols, and timecontrols are not reported in the e(b).)

e(V) Corresponding variance-covariance matrix

e(Nt) A row-vector of the number of treated observations used to compute each estimator

Scalars

e(Nc) the number of control observations used in imputation (scalar)

e(pre_F), e(pre_p), e(pre_df) if pretrends is specified, the F-statistic, pvalue, and dof for the joint test for no pre-trends

e(Niter) the # of iterations to compute SE

Macros

e(cmd) did_imputation

e(droplist) the set of coefficients suppressed to zero because of insufficient effective sample size (see the minn option)

e(autosample_drop) the set of coefficients suppressed to zero because treatment effects could not be imputed for any observation (if autosample is specified)

e(autosample_trim) the set of coefficients where the sample was partially reduced because treatment effects could not be imputed for some observations (if autosample is specified)

Functions

e(sample) Marks the estimation sample: all treated observations for which imputation was successful and the weights are non-zero for at least one coefficient + all non-treated observations used in Step 1

Usage Examples

Conventional panels

1) Estimate the single average treatment-on-the-treated (ATT) across all treated observations, assuming that FE can be imputed for all treated observations (which is rarely the case)
. did_imputation Y i t Ei

2) Same but dropping the observations for which it cannot be imputed. (After running, investigate that the sample is what you think it is!)
. did_imputation Y i t Ei, autosample

3) Estimate the ATT by horizon
. did_imputation Y i t Ei, allhorizons autosample

4) Estimate the ATT at horizons 0..+6 only
. did_imputation Y i t Ei, horizons(0/6)

5) Estimate the ATT at horizons 0..+6 for the subset of units available for all of these horizons (such that the dynamics are not driven by compositional effects)
. did_imputation Y i t Ei, horizons(0/6) hbalance

6) Include time-varying controls:
. did_imputation Y i t Ei, controls(w_first w_other*)

7) Include state-by-year FE
. did_imputation Y county year Ei, fe(county state#year)

8) Drop unit FE
. did_imputation Y i t Ei, fe(t)

9) Additionally report pre-trend coefficients for leads 1,...,5. The estimates for post-treatment effects will NOT change.
. did_imputation Y i t Ei, horizons(0/6) pretrends(5)

10) Estimate the difference between ATE at horizons +2 vs +1 [this can equivalently be done via lincom after estimating ATT by horizon]
. count if K==1
. gen wtr1 = (K==1)/r(N)
. count if K==2
. gen wtr2 = (K==2)/r(N)
. gen wtr_diff = wtr2-wtr1
. did_imputation Y i t Ei, wtr(wtr_diff) sum

11) Reduce estimation time by using loadweights when analyzing several outcomes with identical specifications on identical samples:
. did_imputation Y1 i t Ei, horizons(0/10) saveweights
. rename __* myweights* // optional
. did_imputation Y2 i t Ei, horizons(0/10) loadweights(myweights*)

12) Treatment effect heterogeneity:

To estimate heterogeneity by individual's sex female, you can obtain individual treatment effect estimates
(via saveestimates) and simply run a second-step regression of the estimates on female:

. did_imputation Y i t Ei, saveestimates(tau)

. reg tau female (DON'T DO THIS!)

HOWEVER, standard errors will be incorrect. Instead, use hetby or project:

. did_imputation Y i t Ei, hetby(female)

will produce the ATT for males (tau_0) and females (tau_1). You can further use

. lincom tau_1-tau_0

to compute the difference between them with a SE. Alternatively,

. did_imputation Y i t Ei, project(female)

will produce the ATT for males (tau_cons) and the difference in ATTs between females and males (tau_female).

Both options can be combined with horizons to do heterogeneity analysis within each horizon.

If you also want to allow different period effects for the two groups, you can use the fe() option as usual, e.g.:

. did_imputation Y i t Ei, hetby(female) fe(i t#female)

21) Repeated cross-sections:

When in each period you have a different sample of individiuals i in the same groups (e.g. regions),
replace individual FEs with group FEs and consider clustering at the regional level:

. did_imputation Y i t Ei, fe(region t) cluster(region) ...

Note the main parameters still include i, and not region, as the unit identifier.

22) Triple-diffs:

When observations are defined by i,g,t when, say, i are counties and g are age groups,
specify a variable ig identifying the (i,g) pairs as the unit identifier, add appropriate FEs, and choose your clustering level, e.g.:

. did_imputation Y ig t Eig, fe(ig i#t g#t) cluster(i) ...

Note that the event time Eig should be specific to the i,g pairs, not to the i. For instance, Eig is missing for a never-treated age group in a county where other groups are treated at some point.

Missing Features

- Save imputed Y(0) in addition to treatment effects

- Making hbalance work with autosample

- Throw an error if imputation is not possible with complicated controls

- Allow to use treatment effect averages for SE computation which rely even on observations
outside the estimation sample for the current wtr

- Estimation when treatment switches on and off

- More general interactions between FEs and continuous controls than with timecontrols and unitcontrols

- Frequency weights

- Verify that the unit ID variable is numeric

- In Stata 13 there may be a problem with the df option of test

- loadweights doesn't work when the weights are saved with the default names __w*

- Allow for designs in which treatment is not binary

- Add a check that ranges of t and Ei match

If you are interested in discussing these or others, please contact me.

Bug reporting

If you get an error message, please:

- Reinstall reghdfe and ftools to have the most recent version

- Double check the syntax: e.g. make sure that Ei is the event date and not the treatment dummy, and that
the treatment dummy can be obtained as t>=Ei

- If it's a message with an explanation, read the message carefully.

If this doesn't help:

- Rerun your command adding the verbose option and save the log-file

- If possible, create a version of the dataset in which you can replicate the error (e.g. with a fake outcome variable).

- If you can't share a fake dataset, summarize all the relevant variables in the log-file before calling did_imputation.

- Report this on github or
email me with all of this.

References

Borusyak, Kirill, Xavier Jaravel, and Jann Spiess (2023). "Revisiting Event Study Designs: Robust and Efficient Estimation," Working paper.

Acknowledgements

We thank Kyle Butts for the help in preparing this helpfile.

Author

Kirill Borusyak (UC Berkeley), k.borusyak@berkeley.edu