Note

Download Jupyter notebook: https://docs.doubleml.org/stable/examples/did/py_panel.ipynb.

Python: Panel Data with Multiple Time Periods

In this example, a detailed guide on Difference-in-Differences with multiple time periods using the DoubleML-package. The implementation is based on Callaway and Sant’Anna(2021).

The notebook requires the following packages:

import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

from lightgbm import LGBMRegressor, LGBMClassifier
from sklearn.linear_model import LinearRegression, LogisticRegression

from doubleml.did import DoubleMLDIDMulti
from doubleml.data import DoubleMLPanelData

from doubleml.did.datasets import make_did_CS2021

Data

We will rely on the make_did_CS2021 DGP, which is inspired by Callaway and Sant’Anna(2021) (Appendix SC) and Sant’Anna and Zhao (2020).

We will observe n_obs units over n_periods. Remark that the dataframe includes observations of the potential outcomes y0 and y1, such that we can use oracle estimates as comparisons.

n_obs = 5000
n_periods = 6

df = make_did_CS2021(n_obs, dgp_type=4, n_periods=n_periods, n_pre_treat_periods=3, time_type="datetime")
df["ite"] = df["y1"] - df["y0"]

print(df.shape)
df.head()

(30000, 11)

	id	y	y0	y1	d	t	Z1	Z2	Z3	Z4	ite
0	0	201.501653	201.501653	203.582951	2025-04-01	2025-01-01	-0.846596	-1.537895	0.973802	-1.09159	2.081298
1	0	191.640440	191.640440	192.080043	2025-04-01	2025-02-01	-0.846596	-1.537895	0.973802	-1.09159	0.439603
2	0	181.229363	181.229363	181.717056	2025-04-01	2025-03-01	-0.846596	-1.537895	0.973802	-1.09159	0.487693
3	0	172.665800	173.525095	172.665800	2025-04-01	2025-04-01	-0.846596	-1.537895	0.973802	-1.09159	-0.859295
4	0	165.795306	163.230785	165.795306	2025-04-01	2025-05-01	-0.846596	-1.537895	0.973802	-1.09159	2.564521

Data Details

Here, we slightly abuse the definition of the potential outcomes. \(Y_{i,t}(1)\) corresponds to the (potential) outcome if unit \(i\) would have received treatment at time period \(\mathrm{g}\) (where the group \(\mathrm{g}\) is drawn with probabilities based on \(Z\)).

More specifically

\[\begin{split}\begin{align*} Y_{i,t}(0)&:= f_t(Z) + \delta_t + \eta_i + \varepsilon_{i,t,0}\\ Y_{i,t}(1)&:= Y_{i,t}(0) + \theta_{i,t,\mathrm{g}} + \epsilon_{i,t,1} - \epsilon_{i,t,0} \end{align*}\end{split}\]

where

• \(f_t(Z)\) depends on pre-treatment observable covariates \(Z_1,\dots, Z_4\) and time \(t\)

• \(\delta_t\) is a time fixed effect

• \(\eta_i\) is a unit fixed effect

• \(\epsilon_{i,t,\cdot}\) are time varying unobservables (iid. \(N(0,1)\))

• \(\theta_{i,t,\mathrm{g}}\) correponds to the exposure effect of unit \(i\) based on group \(\mathrm{g}\) at time \(t\)

For the pre-treatment periods the exposure effect is set to

\[\theta_{i,t,\mathrm{g}}:= 0 \text{ for } t<\mathrm{g}\]

such that

\[\mathbb{E}[Y_{i,t}(1) - Y_{i,t}(0)] = \mathbb{E}[\epsilon_{i,t,1} - \epsilon_{i,t,0}]=0 \text{ for } t<\mathrm{g}\]

The DoubleML Coverage Repository includes coverage simulations based on this DGP.

Data Description

The data is a balanced panel where each unit is observed over n_periods starting Janary 2025.

df.groupby("t").size()

t
2025-01-01    5000
2025-02-01    5000
2025-03-01    5000
2025-04-01    5000
2025-05-01    5000
2025-06-01    5000
dtype: int64

The treatment column d indicates first treatment period of the corresponding unit, whereas NaT units are never treated.

Generally, never treated units should take either on the value ``np.inf`` or ``pd.NaT`` depending on the data type (``float`` or ``datetime``).

The individual units are roughly uniformly divided between the groups, where treatment assignment depends on the pre-treatment covariates Z1 to Z4.

df.groupby("d", dropna=False).size()

d
2025-04-01    7674
2025-05-01    7332
2025-06-01    7134
NaT           7860
dtype: int64

Here, the group indicates the first treated period and NaT units are never treated. To simplify plotting and pands

df.groupby("d", dropna=False).size()

d
2025-04-01    7674
2025-05-01    7332
2025-06-01    7134
NaT           7860
dtype: int64

To get a better understanding of the underlying data and true effects, we will compare the unconditional averages and the true effects based on the oracle values of individual effects ite.

# rename for plotting
df["First Treated"] = df["d"].dt.strftime("%Y-%m").fillna("Never Treated")

# Create aggregation dictionary for means
def agg_dict(col_name):
    return {
        f'{col_name}_mean': (col_name, 'mean'),
        f'{col_name}_lower_quantile': (col_name, lambda x: x.quantile(0.05)),
        f'{col_name}_upper_quantile': (col_name, lambda x: x.quantile(0.95))
    }

# Calculate means and confidence intervals
agg_dictionary = agg_dict("y") | agg_dict("ite")

agg_df = df.groupby(["t", "First Treated"]).agg(**agg_dictionary).reset_index()
agg_df.head()

	t	First Treated	y_mean	y_lower_quantile	y_upper_quantile	ite_mean	ite_lower_quantile	ite_upper_quantile
0	2025-01-01	2025-04	208.701432	198.421263	218.766826	-0.052617	-2.361279	2.293900
1	2025-01-01	2025-05	210.748585	201.461408	220.826454	0.052420	-2.271705	2.290892
2	2025-01-01	2025-06	212.502824	202.279470	222.850428	0.001605	-2.289492	2.322057
3	2025-01-01	Never Treated	214.756604	204.696059	224.927327	0.006626	-2.387401	2.349855
4	2025-02-01	2025-04	208.298789	188.085353	227.899479	-0.039705	-2.393998	2.400497

def plot_data(df, col_name='y'):
    """
    Create an improved plot with colorblind-friendly features

    Parameters:
    -----------
    df : DataFrame
        The dataframe containing the data
    col_name : str, default='y'
        Column name to plot (will use '{col_name}_mean')
    """
    plt.figure(figsize=(12, 7))
    n_colors = df["First Treated"].nunique()
    color_palette = sns.color_palette("colorblind", n_colors=n_colors)

    sns.lineplot(
        data=df,
        x='t',
        y=f'{col_name}_mean',
        hue='First Treated',
        style='First Treated',
        palette=color_palette,
        markers=True,
        dashes=True,
        linewidth=2.5,
        alpha=0.8
    )

    plt.title(f'Average Values {col_name} by Group Over Time', fontsize=16)
    plt.xlabel('Time', fontsize=14)
    plt.ylabel(f'Average Value {col_name}', fontsize=14)


    plt.legend(title='First Treated', title_fontsize=13, fontsize=12,
               frameon=True, framealpha=0.9, loc='best')

    plt.grid(alpha=0.3, linestyle='-')
    plt.tight_layout()

    plt.show()

So let us take a look at the average values over time

plot_data(agg_df, col_name='y')

Instead the true average treatment treatment effects can be obtained by averaging (usually unobserved) the ite values.

The true effect just equals the exposure time (in months):

\[ATT(\mathrm{g}, t) = \min(\mathrm{t} - \mathrm{g} + 1, 0) =: e\]

plot_data(agg_df, col_name='ite')

DoubleMLPanelData

Finally, we can construct our DoubleMLPanelData, specifying

• y_col : the outcome

• d_cols: the group variable indicating the first treated period for each unit

• id_col: the unique identification column for each unit

• t_col : the time column

• x_cols: the additional pre-treatment controls

• datetime_unit: unit required for datetime columns and plotting

dml_data = DoubleMLPanelData(
    data=df,
    y_col="y",
    d_cols="d",
    id_col="id",
    t_col="t",
    x_cols=["Z1", "Z2", "Z3", "Z4"],
    datetime_unit="M"
)
print(dml_data)

================== DoubleMLPanelData Object ==================

------------------ Data summary      ------------------
Outcome variable: y
Treatment variable(s): ['d']
Covariates: ['Z1', 'Z2', 'Z3', 'Z4']
Instrument variable(s): None
Time variable: t
Id variable: id
No. Unique Ids: 5000
No. Observations: 30000

------------------ DataFrame info    ------------------
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 30000 entries, 0 to 29999
Columns: 12 entries, id to First Treated
dtypes: datetime64[s](2), float64(8), int64(1), object(1)
memory usage: 2.7+ MB

ATT Estimation

The DoubleML-package implements estimation of group-time average treatment effect via the DoubleMLDIDMulti class (see model documentation).

Basics

The class basically behaves like other DoubleML classes and requires the specification of two learners (for more details on the regression elements, see score documentation).

The basic arguments of a DoubleMLDIDMulti object include

• ml_g “outcome” regression learner

• ml_m propensity Score learner

• control_group the control group for the parallel trend assumption

• gt_combinations combinations of \((\mathrm{g},t_\text{pre}, t_\text{eval})\)

• anticipation_periods number of anticipation periods

We will construct a dict with “default” arguments.

default_args = {
    "ml_g": LGBMRegressor(n_estimators=500, learning_rate=0.01, verbose=-1, random_state=123),
    "ml_m": LGBMClassifier(n_estimators=500, learning_rate=0.01, verbose=-1, random_state=123),
    "control_group": "never_treated",
    "gt_combinations": "standard",
    "anticipation_periods": 0,
    "n_folds": 5,
    "n_rep": 1,
}

The model will be estimated using the fit() method.

dml_obj = DoubleMLDIDMulti(dml_data, **default_args)
dml_obj.fit()
print(dml_obj)

================== DoubleMLDIDMulti Object ==================

------------------ Data summary      ------------------
Outcome variable: y
Treatment variable(s): ['d']
Covariates: ['Z1', 'Z2', 'Z3', 'Z4']
Instrument variable(s): None
Time variable: t
Id variable: id
No. Unique Ids: 5000
No. Observations: 30000

------------------ Score & algorithm ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0

------------------ Machine learner   ------------------
Learner ml_g: LGBMRegressor(learning_rate=0.01, n_estimators=500, random_state=123,
              verbose=-1)
Learner ml_m: LGBMClassifier(learning_rate=0.01, n_estimators=500, random_state=123,
               verbose=-1)
Out-of-sample Performance:
Regression:
Learner ml_g0 RMSE: [[1.94846287 1.96763921 1.92277248 2.87619482 4.10529174 1.94735088
  1.93724715 1.91847971 1.9429786  2.88544865 1.92854563 1.93995659
  1.94271395 1.91905477 1.96746858]]
Learner ml_g1 RMSE: [[1.93455165 1.92955165 1.92884366 2.84860154 3.91175491 1.98806223
  1.94970825 1.93800778 1.98730267 3.01180336 1.89867196 1.9202855
  1.96939867 1.93654801 1.90468294]]
Classification:
Learner ml_m Log Loss: [[0.65159601 0.65411577 0.65133955 0.65896451 0.6554241  0.71138805
  0.70335352 0.70181861 0.70294331 0.70253816 0.71934176 0.72085116
  0.71722285 0.71510649 0.71477689]]

------------------ Resampling        ------------------
No. folds: 5
No. repeated sample splits: 1

------------------ Fit summary       ------------------
                                  coef   std err          t         P>|t|  \
ATT(2025-04,2025-01,2025-02)  0.204730  0.194930   1.050273  2.935924e-01
ATT(2025-04,2025-02,2025-03) -0.177923  0.162847  -1.092576  2.745798e-01
ATT(2025-04,2025-03,2025-04)  1.145894  0.124601   9.196478  0.000000e+00
ATT(2025-04,2025-03,2025-05)  1.905966  0.218096   8.739105  0.000000e+00
ATT(2025-04,2025-03,2025-06)  3.250239  0.388193   8.372735  0.000000e+00
ATT(2025-05,2025-01,2025-02)  0.193202  0.097500   1.981553  4.752932e-02
ATT(2025-05,2025-02,2025-03) -0.056627  0.091892  -0.616237  5.377384e-01
ATT(2025-05,2025-03,2025-04)  0.068216  0.092263   0.739362  4.596874e-01
ATT(2025-05,2025-04,2025-05)  0.800469  0.101132   7.915083  2.442491e-15
ATT(2025-05,2025-04,2025-06)  2.111392  0.156963  13.451543  0.000000e+00
ATT(2025-06,2025-01,2025-02)  0.043119  0.104344   0.413234  6.794348e-01
ATT(2025-06,2025-02,2025-03) -0.126025  0.091910  -1.371177  1.703197e-01
ATT(2025-06,2025-03,2025-04)  0.011008  0.094901   0.115994  9.076572e-01
ATT(2025-06,2025-04,2025-05) -0.112267  0.095811  -1.171758  2.412943e-01
ATT(2025-06,2025-05,2025-06)  1.089835  0.114594   9.510436  0.000000e+00

                                 2.5 %    97.5 %
ATT(2025-04,2025-01,2025-02) -0.177326  0.586786
ATT(2025-04,2025-02,2025-03) -0.497097  0.141252
ATT(2025-04,2025-03,2025-04)  0.901679  1.390108
ATT(2025-04,2025-03,2025-05)  1.478505  2.333426
ATT(2025-04,2025-03,2025-06)  2.489394  4.011084
ATT(2025-05,2025-01,2025-02)  0.002105  0.384298
ATT(2025-05,2025-02,2025-03) -0.236732  0.123478
ATT(2025-05,2025-03,2025-04) -0.112616  0.249047
ATT(2025-05,2025-04,2025-05)  0.602254  0.998684
ATT(2025-05,2025-04,2025-06)  1.803751  2.419034
ATT(2025-06,2025-01,2025-02) -0.161393  0.247630
ATT(2025-06,2025-02,2025-03) -0.306164  0.054115
ATT(2025-06,2025-03,2025-04) -0.174994  0.197010
ATT(2025-06,2025-04,2025-05) -0.300052  0.075518
ATT(2025-06,2025-05,2025-06)  0.865236  1.314435

The summary displays estimates of the \(ATT(g,t_\text{eval})\) effects for different combinations of \((g,t_\text{eval})\) via \(\widehat{ATT}(\mathrm{g},t_\text{pre},t_\text{eval})\), where

• \(\mathrm{g}\) specifies the group

• \(t_\text{pre}\) specifies the corresponding pre-treatment period

• \(t_\text{eval}\) specifies the evaluation period

The choice gt_combinations="standard", used estimates all possible combinations of \(ATT(g,t_\text{eval})\) via \(\widehat{ATT}(\mathrm{g},t_\text{pre},t_\text{eval})\), where the standard choice is \(t_\text{pre} = \min(\mathrm{g}, t_\text{eval}) - 1\) (without anticipation).

Remark that this includes pre-tests effects if \(\mathrm{g} > t_{eval}\), e.g. \(\widehat{ATT}(g=\text{2025-04}, t_{\text{pre}}=\text{2025-01}, t_{\text{eval}}=\text{2025-02})\) which estimates the pre-trend from January to February even if the actual treatment occured in April.

As usual for the DoubleML-package, you can obtain joint confidence intervals via bootstrap.

level = 0.95

ci = dml_obj.confint(level=level)
dml_obj.bootstrap(n_rep_boot=5000)
ci_joint = dml_obj.confint(level=level, joint=True)
ci_joint

	2.5 %	97.5 %
ATT(2025-04,2025-01,2025-02)	-0.351202	0.760662
ATT(2025-04,2025-02,2025-03)	-0.642356	0.286510
ATT(2025-04,2025-03,2025-04)	0.790536	1.501251
ATT(2025-04,2025-03,2025-05)	1.283964	2.527967
ATT(2025-04,2025-03,2025-06)	2.143128	4.357350
ATT(2025-05,2025-01,2025-02)	-0.084865	0.471268
ATT(2025-05,2025-02,2025-03)	-0.318699	0.205445
ATT(2025-05,2025-03,2025-04)	-0.194914	0.331345
ATT(2025-05,2025-04,2025-05)	0.512044	1.088894
ATT(2025-05,2025-04,2025-06)	1.663741	2.559044
ATT(2025-06,2025-01,2025-02)	-0.254467	0.340705
ATT(2025-06,2025-02,2025-03)	-0.388148	0.136098
ATT(2025-06,2025-03,2025-04)	-0.259645	0.281661
ATT(2025-06,2025-04,2025-05)	-0.385514	0.160981
ATT(2025-06,2025-05,2025-06)	0.763019	1.416652

A visualization of the effects can be obtained via the plot_effects() method.

Remark that the plot used joint confidence intervals per default.

dml_obj.plot_effects()

/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 2025-04'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-05'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-06'}, xlabel='Evaluation Period', ylabel='Effect'>])

Sensitivity Analysis

As descripted in the Sensitivity Guide, robustness checks on omitted confounding/parallel trend violations are available, via the standard sensitivity_analysis() method.

dml_obj.sensitivity_analysis()
print(dml_obj.sensitivity_summary)

================== Sensitivity Analysis ==================

------------------ Scenario          ------------------
Significance Level: level=0.95
Sensitivity parameters: cf_y=0.03; cf_d=0.03, rho=1.0

------------------ Bounds with CI    ------------------
                              CI lower  theta lower     theta  theta upper  \
ATT(2025-04,2025-01,2025-02) -0.211950     0.089612  0.204730     0.319848
ATT(2025-04,2025-02,2025-03) -0.553374    -0.216750 -0.177923    -0.139096
ATT(2025-04,2025-03,2025-04)  0.821165     1.018996  1.145894     1.272791
ATT(2025-04,2025-03,2025-05)  1.389366     1.732528  1.905966     2.079403
ATT(2025-04,2025-03,2025-06)  2.411606     3.011661  3.250239     3.488816
ATT(2025-05,2025-01,2025-02) -0.073142     0.086692  0.193202     0.299711
ATT(2025-05,2025-02,2025-03) -0.312092    -0.161266 -0.056627     0.048011
ATT(2025-05,2025-03,2025-04) -0.191431    -0.039642  0.068216     0.176073
ATT(2025-05,2025-04,2025-05)  0.529955     0.695003  0.800469     0.905935
ATT(2025-05,2025-04,2025-06)  1.700513     1.955007  2.111392     2.267777
ATT(2025-06,2025-01,2025-02) -0.229818    -0.059461  0.043119     0.145698
ATT(2025-06,2025-02,2025-03) -0.384576    -0.232476 -0.126025    -0.019573
ATT(2025-06,2025-03,2025-04) -0.251377    -0.094450  0.011008     0.116466
ATT(2025-06,2025-04,2025-05) -0.380411    -0.223130 -0.112267    -0.001403
ATT(2025-06,2025-05,2025-06)  0.804514     0.988712  1.089835     1.190959

                              CI upper
ATT(2025-04,2025-01,2025-02)  0.660710
ATT(2025-04,2025-02,2025-03)  0.199521
ATT(2025-04,2025-03,2025-04)  1.485967
ATT(2025-04,2025-03,2025-05)  2.455573
ATT(2025-04,2025-03,2025-06)  4.168626
ATT(2025-05,2025-01,2025-02)  0.461050
ATT(2025-05,2025-02,2025-03)  0.200254
ATT(2025-05,2025-03,2025-04)  0.328307
ATT(2025-05,2025-04,2025-05)  1.074231
ATT(2025-05,2025-04,2025-06)  2.530724
ATT(2025-06,2025-01,2025-02)  0.319177
ATT(2025-06,2025-02,2025-03)  0.131232
ATT(2025-06,2025-03,2025-04)  0.272298
ATT(2025-06,2025-04,2025-05)  0.156872
ATT(2025-06,2025-05,2025-06)  1.384389

------------------ Robustness Values ------------------
                              H_0     RV (%)    RVa (%)
ATT(2025-04,2025-01,2025-02)  0.0   5.272462   0.000445
ATT(2025-04,2025-02,2025-03)  0.0  13.017999   0.000618
ATT(2025-04,2025-03,2025-04)  0.0  23.981930  20.544369
ATT(2025-04,2025-03,2025-05)  0.0  28.337174  24.315134
ATT(2025-04,2025-03,2025-06)  0.0  33.771239  29.325556
ATT(2025-05,2025-01,2025-02)  0.0   5.374803   0.940494
ATT(2025-05,2025-02,2025-03)  0.0   1.635009   0.000454
ATT(2025-05,2025-03,2025-04)  0.0   1.908135   0.000576
ATT(2025-05,2025-04,2025-05)  0.0  20.600481  16.661749
ATT(2025-05,2025-04,2025-06)  0.0  33.529462  29.796238
ATT(2025-06,2025-01,2025-02)  0.0   1.272344   0.000652
ATT(2025-06,2025-02,2025-03)  0.0   3.541730   0.000599
ATT(2025-06,2025-03,2025-04)  0.0   0.317290   0.000580
ATT(2025-06,2025-04,2025-05)  0.0   3.037450   0.000336
ATT(2025-06,2025-05,2025-06)  0.0  27.878898  23.964444

In this example one can clearly, distinguish the robustness of the non-zero effects vs. the pre-treatment periods.

Control Groups

The current implementation support the following control groups

• "never_treated"

• "not_yet_treated"

Remark that the ``”not_yet_treated” depends on anticipation.

For differences and recommendations, we refer to Callaway and Sant’Anna(2021).

dml_obj_nyt = DoubleMLDIDMulti(dml_data, **(default_args | {"control_group": "not_yet_treated"}))
dml_obj_nyt.fit()
dml_obj_nyt.bootstrap(n_rep_boot=5000)
dml_obj_nyt.plot_effects()

/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 2025-04'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-05'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-06'}, xlabel='Evaluation Period', ylabel='Effect'>])

Linear Covariate Adjustment

Remark that we relied on boosted trees to adjust for conditional parallel trends which allow for a nonlinear adjustment. In comparison to linear adjustment, we could rely on linear learners.

Remark that the DGP (``dgp_type=4``) is based on nonlinear conditional expectations such that the estimates will be biased

linear_learners = {
    "ml_g": LinearRegression(),
    "ml_m": LogisticRegression(),
}

dml_obj_linear = DoubleMLDIDMulti(dml_data, **(default_args | linear_learners))
dml_obj_linear.fit()
dml_obj_linear.bootstrap(n_rep_boot=5000)
dml_obj_linear.plot_effects()

/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 2025-04'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-05'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-06'}, xlabel='Evaluation Period', ylabel='Effect'>])

Aggregated Effects

As the did-R-package, the \(ATT\)’s can be aggregated to summarize multiple effects. For details on different aggregations and details on their interpretations see Callaway and Sant’Anna(2021).

The aggregations are implemented via the aggregate() method.

Group Aggregation

To obtain group-specific effects one can would like to average \(ATT(\mathrm{g}, t_\text{eval})\) over \(t_\text{eval}\). As a sample oracle we will combine all ite’s based on group \(\mathrm{g}\).

df_post_treatment = df[df["t"] >= df["d"]]
df_post_treatment.groupby("d")["ite"].mean()

d
2025-04-01    1.914848
2025-05-01    1.480135
2025-06-01    0.993531
Name: ite, dtype: float64

To obtain group-specific effects it is possible to aggregate several \(\widehat{ATT}(\mathrm{g},t_\text{pre},t_\text{eval})\) values based on the group \(\mathrm{g}\) by setting the aggregation="group" argument.

aggregated_group = dml_obj.aggregate(aggregation="group")
print(aggregated_group)
_ = aggregated_group.plot_effects()

================== DoubleMLDIDAggregation Object ==================
 Group Aggregation

------------------ Overall Aggregated Effects ------------------
    coef  std err        t  P>|t|   2.5 %   97.5 %
1.561451 0.126013 12.39115    0.0 1.31447 1.808433
------------------ Aggregated Effects         ------------------
             coef   std err          t  P>|t|     2.5 %    97.5 %
2025-04  2.100699  0.225016   9.335782    0.0  1.659676  2.541722
2025-05  1.455931  0.118203  12.317200    0.0  1.224257  1.687604
2025-06  1.089835  0.114594   9.510436    0.0  0.865236  1.314435
------------------ Additional Information     ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0

/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/doubleml/did/did_aggregation.py:368: UserWarning: Joint confidence intervals require bootstrapping which hasn't been performed yet. Automatically applying '.aggregated_frameworks.bootstrap(method="normal", n_rep_boot=500)' with default values. For different bootstrap settings, call bootstrap() explicitly before plotting.
  warnings.warn(

The output is a DoubleMLDIDAggregation object which includes an overall aggregation summary based on group size.

Time Aggregation

To obtain time-specific effects one can would like to average \(ATT(\mathrm{g}, t_\text{eval})\) over \(\mathrm{g}\) (respecting group size). As a sample oracle we will combine all ite’s based on group \(\mathrm{g}\). As oracle values, we obtain

df_post_treatment.groupby("t")["ite"].mean()

t
2025-04-01    0.946131
2025-05-01    1.432781
2025-06-01    1.992559
Name: ite, dtype: float64

To aggregate \(\widehat{ATT}(\mathrm{g},t_\text{pre},t_\text{eval})\), based on \(t_\text{eval}\), but weighted with respect to group size. Corresponds to Calendar Time Effects from the did-R-package.

For calendar time effects set aggregation="time".

aggregated_time = dml_obj.aggregate("time")
print(aggregated_time)
fig, ax = aggregated_time.plot_effects()

================== DoubleMLDIDAggregation Object ==================
 Time Aggregation

------------------ Overall Aggregated Effects ------------------
   coef  std err         t  P>|t|    2.5 %   97.5 %
1.56289 0.136353 11.462053    0.0 1.295642 1.830138
------------------ Aggregated Effects         ------------------
             coef   std err          t  P>|t|     2.5 %    97.5 %
2025-04  1.145894  0.124601   9.196478    0.0  0.901679  1.390108
2025-05  1.365815  0.138669   9.849468    0.0  1.094029  1.637601
2025-06  2.176962  0.189422  11.492663    0.0  1.805702  2.548222
------------------ Additional Information     ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0

/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/doubleml/did/did_aggregation.py:368: UserWarning: Joint confidence intervals require bootstrapping which hasn't been performed yet. Automatically applying '.aggregated_frameworks.bootstrap(method="normal", n_rep_boot=500)' with default values. For different bootstrap settings, call bootstrap() explicitly before plotting.
  warnings.warn(

Event Study Aggregation

To obtain event-study-type effects one can would like to aggregate \(ATT(\mathrm{g}, t_\text{eval})\) over \(e = t_\text{eval} - \mathrm{g}\) (respecting group size). As a sample oracle we will combine all ite’s based on group \(\mathrm{g}\). As oracle values, we obtain

df["e"] = pd.to_datetime(df["t"]).values.astype("datetime64[M]") - \
    pd.to_datetime(df["d"]).values.astype("datetime64[M]")
df.groupby("e")["ite"].mean()[1:]

e
-122 days    0.017831
-92 days    -0.003497
-61 days    -0.002369
-31 days    -0.031005
0 days       0.978670
31 days      1.903669
59 days      2.950480
Name: ite, dtype: float64

Analogously, aggregation="eventstudy" aggregates \(\widehat{ATT}(\mathrm{g},t_\text{pre},t_\text{eval})\) based on exposure time \(e = t_\text{eval} - \mathrm{g}\) (respecting group size).

aggregated_eventstudy = dml_obj.aggregate("eventstudy")
print(aggregated_eventstudy)
aggregated_eventstudy.plot_effects()

================== DoubleMLDIDAggregation Object ==================
 Event Study Aggregation

------------------ Overall Aggregated Effects ------------------
    coef  std err         t  P>|t|    2.5 %   97.5 %
2.090005 0.197902 10.560791    0.0 1.702123 2.477886
------------------ Aggregated Effects         ------------------
               coef   std err          t     P>|t|     2.5 %    97.5 %
-4 months  0.043119  0.104344   0.413234  0.679435 -0.161393  0.247630
-3 months  0.035773  0.069842   0.512204  0.608508 -0.101114  0.172660
-2 months  0.055756  0.088993   0.626522  0.530973 -0.118667  0.230179
-1 months -0.075255  0.080544  -0.934326  0.350136 -0.233118  0.082609
0 months   1.013438  0.081955  12.365831  0.000000  0.852810  1.174066
1 months   2.006338  0.163146  12.297821  0.000000  1.686578  2.326098
2 months   3.250239  0.388193   8.372735  0.000000  2.489394  4.011084
------------------ Additional Information     ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0

/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/doubleml/did/did_aggregation.py:368: UserWarning: Joint confidence intervals require bootstrapping which hasn't been performed yet. Automatically applying '.aggregated_frameworks.bootstrap(method="normal", n_rep_boot=500)' with default values. For different bootstrap settings, call bootstrap() explicitly before plotting.
  warnings.warn(

(<Figure size 1200x600 with 1 Axes>,
 <Axes: title={'center': 'Aggregated Treatment Effects'}, ylabel='Effect'>)

Aggregation Details

The DoubleMLDIDAggregation objects include several DoubleMLFrameworks which support methods like bootstrap() or confint(). Further, the weights can be accessed via the properties

• overall_aggregation_weights: weights for the overall aggregation

• aggregation_weights: weights for the aggregation

To clarify, e.g. for the eventstudy aggregation

print(aggregated_eventstudy)

================== DoubleMLDIDAggregation Object ==================
 Event Study Aggregation

------------------ Overall Aggregated Effects ------------------
    coef  std err         t  P>|t|    2.5 %   97.5 %
2.090005 0.197902 10.560791    0.0 1.702123 2.477886
------------------ Aggregated Effects         ------------------
               coef   std err          t     P>|t|     2.5 %    97.5 %
-4 months  0.043119  0.104344   0.413234  0.679435 -0.161393  0.247630
-3 months  0.035773  0.069842   0.512204  0.608508 -0.101114  0.172660
-2 months  0.055756  0.088993   0.626522  0.530973 -0.118667  0.230179
-1 months -0.075255  0.080544  -0.934326  0.350136 -0.233118  0.082609
0 months   1.013438  0.081955  12.365831  0.000000  0.852810  1.174066
1 months   2.006338  0.163146  12.297821  0.000000  1.686578  2.326098
2 months   3.250239  0.388193   8.372735  0.000000  2.489394  4.011084
------------------ Additional Information     ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0

Here, the overall effect aggregation aggregates each effect with positive exposure

print(aggregated_eventstudy.overall_aggregation_weights)

[0.         0.         0.         0.         0.33333333 0.33333333
 0.33333333]

If one would like to consider how the aggregated effect with \(e=0\) is computed, one would have to look at the corresponding set of weights within the aggregation_weights property

# the weights for e=0 correspond to the fifth element of the aggregation weights
aggregated_eventstudy.aggregation_weights[4]

array([0.        , 0.        , 0.34661247, 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.33116531, 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.32222222])

Taking a look at the original dml_obj, one can see that this combines the following estimates (only show month):

• \(\widehat{ATT}(04,03,04)\)

• \(\widehat{ATT}(05,04,05)\)

• \(\widehat{ATT}(06,05,06)\)

print(dml_obj.summary["coef"])

ATT(2025-04,2025-01,2025-02)    0.204730
ATT(2025-04,2025-02,2025-03)   -0.177923
ATT(2025-04,2025-03,2025-04)    1.145894
ATT(2025-04,2025-03,2025-05)    1.905966
ATT(2025-04,2025-03,2025-06)    3.250239
ATT(2025-05,2025-01,2025-02)    0.193202
ATT(2025-05,2025-02,2025-03)   -0.056627
ATT(2025-05,2025-03,2025-04)    0.068216
ATT(2025-05,2025-04,2025-05)    0.800469
ATT(2025-05,2025-04,2025-06)    2.111392
ATT(2025-06,2025-01,2025-02)    0.043119
ATT(2025-06,2025-02,2025-03)   -0.126025
ATT(2025-06,2025-03,2025-04)    0.011008
ATT(2025-06,2025-04,2025-05)   -0.112267
ATT(2025-06,2025-05,2025-06)    1.089835
Name: coef, dtype: float64

Anticipation

As described in the Model Guide, one can include anticipation periods \(\delta>0\) by setting the anticipation_periods parameter.

Data with Anticipation

The DGP allows to include anticipation periods via the anticipation_periods parameter. In this case the observations will be “shifted” such that units anticipate the effect earlier and the exposure effect is increased by the number of periods where the effect is anticipated.

n_obs = 4000
n_periods = 6

df_anticipation = make_did_CS2021(n_obs, dgp_type=4, n_periods=n_periods, n_pre_treat_periods=3, time_type="datetime", anticipation_periods=1)

print(df_anticipation.shape)
df_anticipation.head()

(19068, 10)

	id	y	y0	y1	d	t	Z1	Z2	Z3	Z4
1	0	190.616369	190.616369	190.105414	2025-05-01	2025-01-01	-1.446744	1.447195	-0.398511	-0.386465
2	0	179.947855	179.947855	178.818314	2025-05-01	2025-02-01	-1.446744	1.447195	-0.398511	-0.386465
3	0	169.658986	169.658986	168.386431	2025-05-01	2025-03-01	-1.446744	1.447195	-0.398511	-0.386465
4	0	159.113416	160.145743	159.113416	2025-05-01	2025-04-01	-1.446744	1.447195	-0.398511	-0.386465
5	0	150.938473	148.569647	150.938473	2025-05-01	2025-05-01	-1.446744	1.447195	-0.398511	-0.386465

To visualize the anticipation, we will again plot the “oracle” values

df_anticipation["ite"] = df_anticipation["y1"] - df_anticipation["y0"]
df_anticipation["First Treated"] = df_anticipation["d"].dt.strftime("%Y-%m").fillna("Never Treated")
agg_df_anticipation = df_anticipation.groupby(["t", "First Treated"]).agg(**agg_dictionary).reset_index()
agg_df_anticipation.head()

	t	First Treated	y_mean	y_lower_quantile	y_upper_quantile	ite_mean	ite_lower_quantile	ite_upper_quantile
0	2025-01-01	2025-04	208.513750	191.079873	225.985937	0.026761	-2.384698	2.295447
1	2025-01-01	2025-05	211.412458	194.894738	228.821351	0.003658	-2.296333	2.351912
2	2025-01-01	2025-06	213.013327	195.822673	230.140195	0.060842	-2.393095	2.292370
3	2025-01-01	Never Treated	217.179294	198.872760	233.519933	-0.047482	-2.310247	2.412780
4	2025-02-01	2025-04	208.359647	181.668285	233.830646	-0.010302	-2.345604	2.333478

One can see that the effect is already anticipated one period before the actual treatment assignment.

plot_data(agg_df_anticipation, col_name='ite')

Initialize a corresponding DoubleMLPanelData object.

dml_data_anticipation = DoubleMLPanelData(
    data=df_anticipation,
    y_col="y",
    d_cols="d",
    id_col="id",
    t_col="t",
    x_cols=["Z1", "Z2", "Z3", "Z4"],
    datetime_unit="M"
)

ATT Estimation

Let us take a look at the estimation without anticipation.

dml_obj_anticipation = DoubleMLDIDMulti(dml_data_anticipation, **default_args)
dml_obj_anticipation.fit()
dml_obj_anticipation.bootstrap(n_rep_boot=5000)
dml_obj_anticipation.plot_effects()

/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/doubleml/double_ml.py:1479: UserWarning: The estimated nu2 for d is not positive. Re-estimation based on riesz representer (non-orthogonal).
  warnings.warn(msg, UserWarning)
/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 2025-04'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-05'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-06'}, xlabel='Evaluation Period', ylabel='Effect'>])

The effects are obviously biased. To include anticipation periods, one can adjust the anticipation_periods parameter. Correspondingly, the outcome regression (and not yet treated units) are adjusted.

dml_obj_anticipation = DoubleMLDIDMulti(dml_data_anticipation, **(default_args| {"anticipation_periods": 1}))
dml_obj_anticipation.fit()
dml_obj_anticipation.bootstrap(n_rep_boot=5000)
dml_obj_anticipation.plot_effects()

/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)
/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)
/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)
/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)
/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 2025-04'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-05'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-06'}, xlabel='Evaluation Period', ylabel='Effect'>])

Group-Time Combinations

The default option gt_combinations="standard" includes all group time values with the specific choice of \(t_\text{pre} = \min(\mathrm{g}, t_\text{eval}) - 1\) (without anticipation) which is the weakest possible parallel trend assumption.

Other options are possible or only specific combinations of \((\mathrm{g},t_\text{pre},t_\text{eval})\).

All Combinations

The option gt_combinations="all" includes all relevant group time values with \(t_\text{pre} < \min(\mathrm{g}, t_\text{eval})\), including longer parallel trend assumptions. This can result in multiple estimates for the same \(ATT(\mathrm{g},t)\), which have slightly different assumptions (length of parallel trends).

dml_obj_all = DoubleMLDIDMulti(dml_data, **(default_args| {"gt_combinations": "all"}))
dml_obj_all.fit()
dml_obj_all.bootstrap(n_rep_boot=5000)
dml_obj_all.plot_effects()

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 2025-04'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-05'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-06'}, xlabel='Evaluation Period', ylabel='Effect'>])

Universal Base Period

The option gt_combinations="universal" set \(t_\text{pre} = \mathrm{g} - \delta - 1\), corresponding to a universal/constant comparison or base period.

Remark that this implies \(t_\text{pre} > t_\text{eval}\) for all pre-treatment periods (accounting for anticipation). Therefore these effects do not have the same straightforward interpretation as ATT’s.

dml_obj_universal = DoubleMLDIDMulti(dml_data, **(default_args| {"gt_combinations": "universal"}))
dml_obj_universal.fit()
dml_obj_universal.bootstrap(n_rep_boot=5000)
dml_obj_universal.plot_effects()

/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/matplotlib/cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 2025-04'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-05'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 2025-06'}, xlabel='Evaluation Period', ylabel='Effect'>])

Selected Combinations

Instead it is also possible to just submit a list of tuples containing \((\mathrm{g}, t_\text{pre}, t_\text{eval})\) combinations. E.g. only two combinations

gt_dict = {
    "gt_combinations": [
        (np.datetime64('2025-04'),
         np.datetime64('2025-01'),
         np.datetime64('2025-02')),
        (np.datetime64('2025-04'),
         np.datetime64('2025-02'),
         np.datetime64('2025-03')),
    ]
}

dml_obj_all = DoubleMLDIDMulti(dml_data, **(default_args| gt_dict))
dml_obj_all.fit()
dml_obj_all.bootstrap(n_rep_boot=5000)
dml_obj_all.plot_effects()

(<Figure size 1200x800 with 2 Axes>,
 [<Axes: title={'center': 'First Treated: 2025-04'}, xlabel='Evaluation Period', ylabel='Effect'>])