did-imputation 0.1.0

Borusyak, Jaravel & Spiess (2024) difference-in-differences imputation estimator with event-study plots.

did_imputation

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, Jaravel, and Spiess (2024).

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:

• 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.
• 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. Treated observations for which imputation of the counterfactual outcome is not possible and treatment effect cannot be estimated will be dropped. A warning is displayed in that case, as the estimand should not be interpreted as the average treatment effect for all treated observations. (See What if imputation is not possible)
• Take averages of estimated treatment effects corresponding to the estimand of interest.

A pretrend test (for the assumptions of parallel trends and no anticipation) is a separate exercise. Regardless of whether the pretrend test is performed, the reference group for estimation is always all untreated (i.e., pre-treatment and never-treated) observations.

To make "event study" plots, please use the accompanying command event_plot (documented below).

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 groupXperiod 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. 2024).

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 apply the autosample option that will exclude those observations automatically in most cases.

In case this does not work, please modify the estimand, excluding those observations manually: altering the data frame or by setting the weights on them to zero (via the wtr option);

Installation

To install or update the package, use pip:

pip install did_imputation

Dependencies

pandas, numpy, pyhdfe, statsmodels, matplotlib, scipy

Usage

Import

from did_imputation import did_imputation, event_plot

Command Structure

The structure of the command mimics the one in Stata. The basic setup is:

output = did_imputation(df, "y", "i", "t", "Ei")

where:

• df: Pandas dataFrame containing the data
• y: Outcome variable name in dataframe
• i: Entity/unit identifier
• t: Time variable
• Ei: Treatment timing variable

Note: These main parameters imply:

• the treatment indicator: $D_{it}=\mathbf {1}[t\ge E_i]$;
• "relative time", i.e. the number of periods since treatment: $K_{it}=(t-E_i)$ (possibly adjusted by the shift and delta options described below).

Options

• fe (list of variable names in DataFrame): which FEs 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 include more FEs. If you want no FEs at all, specify fe = None.

• controls (list of variable names in DataFrame): list of continuous controls.

• aw (string name of variable in DataFrame): analytic weights, or estimation weights, which are used both in Step 1 (when estimating the fixed effects) and in Step 3 (applied when averaging causal effects). For instance, one may use regional population weights when larger regions have a lower variance of the error term (for efficiency) and a population-weighted average of effects is of interest. Default is None, meaning all observations have the same weight.

• wtr (list of variable names in DataFrame): A list of variables, manually defining estimands of interest by storing the weights on the treated observations.

  • If estimation weights (aw) 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).
  • Using multiple wtr variables is faster than running did_imputation for each of them separately, and produces a joint variance-covariance matrix.

• sum (bool, default False): 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 women).

• horizons (list of integers): 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.), instead of an overall ATT. Horizons which are not specified will be ignored. Each horizon must a be non-negative integer (since the pre-trend test is a separate exercise).

• allhorizons (bool, default False): picks all non-negative horizons available in the sample.

• hbalance (bool, default False): 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).

• minn (float, default 30): 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.

• shift (integer, default 0): specify to allow for anticipation effects. The command will pretend that treatment happened shift periods earlier for each treated unit. This option is not for pretrend testing; use it if anticipation effects are expected in your setting. (It can be used for a placebo test but we recommend a pretrend test instead; see Section 4.4 of Borusyak et al. 2024.) The command's output will be labeled relative to the shifted treated date Ei-shift. For example, with horizons = list(range(0, 11)) shift(3) you will get coefficients tau0...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.

• delta (integer): indicates that one period should correpond to delta steps of t and Ei. Default is the mode of difference between two consequent time periods.

• pretrends (integer, default 0): if some value k>0 is specified, the command will perform 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 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).
  • For the same reason, it is not recommended to include too many periods in the pretrend test
  • 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.

• cluster (string variable name in DataFrame, default: same as i): cluster SE within groups of units defined by this variable. Two-way clustering is currently not supported. Note that serial correlation should always be allowed.

• avgeffectsby (list of variable names in DataFrame): Use this option if you have small cohorts of treated observations, and after reviewing Section 4.3 of Borusyak et al. (2024). 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, i.e. 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.

• nose (bool, default False): do not produce standard errors (much faster).

• saveweights (bool, default False): if specified, imputation weights are included in the output as pandas DataFrame with columns 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 the analytic 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. the option sum is not specified) and there are no analytic 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 the analytic weights).

Output

A class object with attributes:

• estimates, std_errors: estimates of average effect or by horizon/wtr coefficients and corresponding SEs
• pretrends_estimates, pretrends_std_errors: pretrends coefficients and corresponding SEs (if pretrends option specified)
• controls_estimates, controls_std_errors: contains estimates of control coefficients and corresponding SEs (if the corresponding option is specified)
• n_obs: the number of observations (treated plus untreated) used in the analysis. If autosample drops some observations, they are not counted
• weights: pandas DataFrame with imputation weights if option saveweights is True in the main command
• V: the variance-covariance matrix of effects, pretrend estimates, and control coefficients

Usage Examples

• The command produces output in the form of a class object. Estimates and standard errors are stored as class attributes. To get a description of output write:

    output = did_imputation(df, "y", "i", "t", "Ei")
    output.summarise()

• To estimate dynamic ATT effects, use horizons or allhorizons options. For example:

    output = did_imputation(df, "y", "i", "t", "Ei", allhorizons = True)
    output = did_imputation(df, "y", "i", "t", "Ei", horizons = list(range(0,5)))
    output = did_imputation(df, "y", "i", "t", "Ei", horizons = [0, 1, 2, 5])

• To estimate ATT in a dynamic model on the subset of units available for all of these horizons (such that the dynamics are not driven by compositional effects):

    output = did_imputation(df, "y", "i", "t", "Ei", allhorizons = True, hbalance = True)
    output = did_imputation(df, "y", "i", "t", "Ei", horizons = list(range(0,5)), hbalance = True)

• To include a pretrend test with 3 pretrend coefficients:

    output = did_imputation(df, "y", "i", "t", "Ei", allhorizons = True, pretrends = 3)

• By default, time and entity fixed effects are included and standard errors are clustered by entity. To change it, use the corresponding options:

    output = did_imputation(df, "y", "i", "t", "Ei", fe=["other_fe1", "other_fe2"], cluster = "other_fe1")

• Add state by year fe:

    state_by_year = list(zip(df['state'], df['year']))
    df['state_by_year'] = pd.factorize(state_by_year)[0]

    output = did_imputation(df, "y", "i", "year", "Ei", fe = ["state_by_year"])

• In a repeated cross-section setting where different individuals are observed within each region in each period, adjust clustering and FEs accordingly (and note the main parameters still include i, and not region, as the unit identifier):

    output = did_imputation(df, "y", "i", "t", "Ei", fe=["region", "t"], cluster = "region")

• Triple-differences estimation: consider a dataset by state i, age group g, and time t, with treatment happening at time Eig. 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 state where other groups are treated at some point. Use unique id for state-age as entity.

    pairs = list(zip(df['state'], df['age_group']))
    df['state_age_group'] = pd.factorize(pairs)[0]
    pairs = list(zip(df['state'], df['year']))
    df['state_year'] = pd.factorize(pairs)[0]
    pairs = list(zip(df['age_group'], df['year']))
    df['age_group_year'] = pd.factorize(pairs)[0]

    output = did_imputation(df, "y", "state_age_group", "year", "Eig", cluster="state", fe = ["state_age_group", "state_year", "age_group_year"])

• If the cohorts are small, consider using avgeffectsby option, e.g. if D is the dummy of being treated:

    output = did_imputation(df, "y", "i", "t", "Ei", avgeffectsby = ["D"])

• Estimate the difference between ATE for men versus women:

    df["D"] = (df["t"] >= df["Ei"])
    df["men"] = ((df["gender"] == "man") & (df["D"] == 1))
    df["women"] = ((df["gender"] == "woman") & (df["D"] == 1))
    df["men"] = df["men"]/np.sum(df["men"])
    df["women"] = df["women"]/np.sum(df["women"])
    df["wtr_dif"] = df["men"] - df["women"]
    output = did_imputation(df, "y", "i", "t", "Ei", wtr = ["wtr_dif"], sum=True)

Current Limitations

The following Stata options are not yet implemented:

• unitcontrols (controls interacted with unit dummies)
• timecontrols (controls interacted with period dummies)
• leaveoneout (leave-out conservative standard errors)
• project, hetby (for estimation of treatment effect heterogeneity)
• saveestimates (save estimated causal effects for each treated observation)
• saveresiduals (save residuals used for constructing SEs)

---

event_plot

event_plot - Plot the staggered-adoption diff-in-diff ("event study") estimates as errorbar: coefficients post treatment ("lags") and, if available, pretrend coefficients ("leads") along with confidence intervals (CIs).

This command is used once estimates have been produced by the imputation estimator of Borusyak et al. 2024 (did_imputation), other methods robust to treatment effect heterogeneity, or conventional event-study OLS.

Syntax

The command has two syntaxes. When using output from the did_imputation command:

output_figure = event_plot(results_obj = did_imputation_output)

When using user-provided estimates, e.g. from other packages: Effects and pretrends are expected to be in the same format as in did_imputation output, i.e. a dictionary with keys "tau1" ... "tauk" for effects, "pre1"..."prek" for pretrends, and corresponding lists for standard errors:

output_figure = event_plot(pretrends={"pre1": -0.05, "pre2": 0.02}, pretrends_std={"pre1": 0.04, "pre2": 0.03}, effects={"tau0": 0.15, "tau1": 0.25, "tau2": 0.35}, effects_std={"tau0": 0.05, "tau1": 0.06, "tau2": 0.07})

Options

• results_obj (did_imputation output object, default None): Results object with results of estimation using did_imputation command with horizons/allhorizons and, optionally, pretrends options specified. If results_obj is provided, it overrides the separate pretrends, pretrends_std, effects, and effects_std parameters. The user must specify either this option or pretrends and effects estimates separately.

• pretrends (dictionary with keys pre1, ..., prek (or just 1 ... k) (in ascending order) and estimates as values, default None): The pretrend coefficients to plot. Can be None when the effects are specified if no pretrends are estimated. The option results_obj should be None if this is specified.

• pretrends_std (dictionary with keys pre1, ..., prek (or just 1 ... k) (in ascending order) and std. errors as values, default None): The standard errors for the pretrend coefficients, to go along with the pretrends option.

• effects (dictionary with keys tau0, ..., tauk (or just 0 ... k) (in ascending order) and estimates as values, default None): The post-treatment causal effects (horizon-specific coefficients) to plot. The option results_obj should be None if this is specified.

• effects_std (dictionary with keys tau0, ..., tauk (or just 0 ... k) (in ascending order) and std. errors as values, default None): The standard errors for the horizon coefficients, to go along with the effects option.

• plot_type (string, default "rcap"): The type of plot "rcap" or "rarea" to be displayed. "rcap" means that only point estimates with confidence intervals will be displayed. "rarea" fills the space between upper and lower bounds of the confidence intervals and connects them.

• together (bool, default False): Whether to plot pretrends and effects on the same series in one style without separation. The default is False to highlight that pretrends are estimated separately from the effects.

• significance_level (float, default 0.05): The significance level for calculating confidence intervals (e.g., 0.05 for 95% CI).

• pretrends_marker (string, default "o"): The marker style for pretrends points (e.g., "o" for circle). Use kwargs instead of this option if together is specified.

• effects_marker (string, default "s"): The marker style for effects points (e.g., "s" for square). Use kwargs instead of this option if together is specified.

• pretrends_color (string, default "blue"): The color for pretrends points and lines. Use kwargs instead of this option if together is specified.

• effects_color (string, default "red"): The color for effects points and lines. Use kwargs instead of this option if together is specified.

• pretrends_linestyle (string, default "-"): The line style for pretrends lines (e.g., "-" for solid). Use kwargs instead of this option if together is specified.

• effects_linestyle (string, default "-"): The line style for effects lines (e.g., "-" for solid). Use kwargs instead of this option if together is specified.

• pretrends_marker_size (int, default 8): The size of the pretrends markers. Use kwargs instead of this option if together is specified.

• effects_marker_size (int, default 8): The size of the effects markers. Use kwargs instead of this option if together is specified.

• pretrends_line_width (float, default 2): The width of the pretrends lines. Use kwargs instead of this option if together is specified.

• effects_line_width (float, default 2): The width of the effects lines. Use kwargs instead of this option if together is specified.

• pretrends_alpha (float, default 0.7): The transparency level for the pretrends (0 to 1). Use kwargs instead of this option if together is specified.

• effects_alpha (float, default 0.7): The transparency level for the effects (0 to 1). Use kwargs instead of this option if together is specified.

• figsize (tuple, default (12, 8)): The figure size as a tuple (width, height) in inches.

• title (string, default ""): The title of the plot.

• xlabel (string, default "Time Relative to Treatment"): The label for the x-axis.

• ylabel (string, default "Coefficient"): The label for the y-axis.

• show_grid (bool, default True): Whether to display grid lines on the plot.

• show_legend (bool, default True): Whether to display the legend on the plot.

• show_zero_line (bool, default True): Whether to display a horizontal zero line on the plot.

• zero_line_position (float, default 0): The y-position of the zero line.

• zero_line_style (string, default "--"): The line style for the zero line (e.g., "--" for dashed).

• zero_line_color (string, default "grey"): The color for the zero line.

• zero_line_width (float, default 1): The width of the zero line.

• zero_line_alpha (float, default 0.5): The transparency level for the zero line (0 to 1).

• show_event_line (bool, default True): Whether to display a vertical event line on the plot.

• event_line_position (float, default 0): The x-position of the event line.

• event_line_style (string, default "--"): The line style for the event line (e.g., "--" for dashed).

• event_line_color (string, default "black"): The color for the event line.

• event_line_width (float, default 1): The width of the event line.

• event_line_alpha (float, default 0.5): The transparency level for the event line (0 to 1).

• pretrends_rarea_color (string, default "blue"): The color for pretrends rarea.

• effects_rarea_color (string, default "red"): The color for effects rarea. Is used for all graph if together is specified.

• pretrends_rarea_alpha (float, default 0.3): The transparency level for the rarea plot for pretrends. Works if plot_type = "rarea" is specified (0 to 1).

• effects_rarea_alpha (float, default 0.3): The transparency level for the rarea plot for effects. Works if plot_type = "rarea" is specified (0 to 1). Is used for all graph if together is specified.

• pretrends_rarea_edgecolor (string, default None): The color of edges for pretrends rarea.

• effects_rarea_edgecolor (string, default None): The color of edges for effects rarea. Is used for all graph if together is specified.

• pretrends_rarea_hatch (string, default None): The hatch of pretrends area, e.g. "/" for diagonal hatch.

• effects_rarea_hatch (string, default None): The hatch of effects area, e.g. "/" for diagonal hatch. Is used for all graph if together is specified.

• font_size (integer, default 12): The font size for text elements like labels and title.

• legend_loc (string, default "best"): The location of the legend (e.g., "best", "upper right").

• save_path (string, default None): The file path to save the figure. If None, the figure is not saved.

• dpi (integer, default 300): The resolution in dots per inch for saving the figure.

• kwargs (dictionary): Additional keyword arguments to pass to matplotlib's errorbar. If specified, added both to pretrends and horizons errorbar or in the common errorbar if together is specified.

Usage Examples

• The command event_plot supplied in the package can be used to plot dynamic estimates directly from did_imputation command:

    output = did_imputation(df, "y", "i", "t", "Ei", pretrends = 4, allhorizons=True)
    fig = event_plot(results_obj = output)

• To plot effects and pretrends in the same style use together option:

    output = did_imputation(df, "y", "i", "t", "Ei", pretrends = 4, allhorizons=True)
    fig = event_plot(results_obj = output, together = True)

• Graphical options from matplotlib.pyplot errorbar can be added as kwargs. They apply to both effects and pretrends. For example, the following command adjusts the capsize and adds pink edges for the markers:

    fig = event_plot(results_obj = output, together=True, capsize = 2, markeredgecolor = "pink" )

• When together is False, the options described above can adjust pretrends and effects graphs separately. They can be combined with errorbar options same for effects and pretrends as kwargs:

    fig = event_plot(results_obj = output, pretrends_color = "green", effects_color = "grey", fillstyle = "left", capsize = 2 )

• Rarea plot can be displayed:

    fig = event_plot(results_obj = output, plot_type = "rarea")

• To adjust options of rarea such as color, hatch, alpha etc., use corresponding options provided separately for pretrends and effects:

    fig = event_plot(results_obj = output, plot_type = "rarea", effects_rarea_alpha = 1, pretrends_rarea_alpha = 0.8)

• To adjust options of rarea when together is True, use options of the graph corresponding to effects:

    fig = event_plot(results_obj = output, plot_type = "rarea", together = True, effects_rarea_color = "navy")

• Options for "rcap" graph still work for rarea, e.g. one can adjust marker type using the same options:

    fig = event_plot(results_obj = output, plot_type = "rarea", pretrends_marker = "D", effects_marker = "o")

event_plot — Current Limitations

• Shift option is not realized yet for event_plot

References

If using this command, please cite:

"Revisiting Event Study Designs: Robust and Efficient Estimation"
(Kirill Borusyak, Xavier Jaravel and Jann Spiess)
Review of Economic Studies, 2024, 91(6), p.3253-3285

Author: Georgii Marinichev ([e-mail]), based on the Stata command by Kirill Borusyak ([e-mail]), with help from Ilaria Dal Barco.