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:
[1]:
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.
[2]:
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)
[2]:
id | y | y0 | y1 | d | t | Z1 | Z2 | Z3 | Z4 | ite | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 206.875950 | 206.875950 | 210.244681 | NaT | 2025-01-01 | -0.636072 | -0.306817 | -0.008854 | -0.713575 | 3.368731 |
1 | 0 | 205.583966 | 205.583966 | 203.155250 | NaT | 2025-02-01 | -0.636072 | -0.306817 | -0.008854 | -0.713575 | -2.428716 |
2 | 0 | 200.146297 | 200.146297 | 199.230898 | NaT | 2025-03-01 | -0.636072 | -0.306817 | -0.008854 | -0.713575 | -0.915399 |
3 | 0 | 195.186428 | 195.186428 | 196.007859 | NaT | 2025-04-01 | -0.636072 | -0.306817 | -0.008854 | -0.713575 | 0.821431 |
4 | 0 | 191.320321 | 191.320321 | 193.731495 | NaT | 2025-05-01 | -0.636072 | -0.306817 | -0.008854 | -0.713575 | 2.411174 |
Data Details#
Here, we slightly abuse the definition of the potential outcomes. :math:`Y_{i,t}(1)` corresponds to the (potential) outcome if unit :math:`i` would have received treatment at time period :math:`mathrm{g}` (where the group :math:`mathrm{g}` is drawn with probabilities based on :math:`Z`).
More specifically
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
such that
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.
[3]:
df.groupby("t").size()
[3]:
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
.
[4]:
df.groupby("d", dropna=False).size()
[4]:
d
2025-04-01 7908
2025-05-01 7248
2025-06-01 7230
NaT 7614
dtype: int64
Here, the group indicates the first treated period and NaT
units are never treated. To simplify plotting and pands
[5]:
df.groupby("d", dropna=False).size()
[5]:
d
2025-04-01 7908
2025-05-01 7248
2025-06-01 7230
NaT 7614
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
.
[6]:
# 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()
[6]:
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.812223 | 199.096920 | 218.885511 | -0.083445 | -2.313984 | 2.212231 |
1 | 2025-01-01 | 2025-05 | 210.245070 | 200.760505 | 219.621769 | -0.023824 | -2.292686 | 2.214080 |
2 | 2025-01-01 | 2025-06 | 212.473879 | 202.076274 | 222.517507 | 0.060412 | -2.272092 | 2.360931 |
3 | 2025-01-01 | Never Treated | 214.227001 | 204.229352 | 224.505554 | 0.003376 | -2.357390 | 2.408004 |
4 | 2025-02-01 | 2025-04 | 208.569308 | 189.361738 | 228.191098 | 0.031492 | -2.277407 | 2.405074 |
[7]:
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
[8]:
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):
[9]:
plot_data(agg_df, col_name='ite')

DoubleMLPanelData#
Finally, we can construct our DoubleMLPanelData
, specifying
y_col
: the outcomed_cols
: the group variable indicating the first treated period for each unitid_col
: the unique identification column for each unitt_col
: the time columnx_cols
: the additional pre-treatment controlsdatetime_unit
: unit required fordatetime
columns and plotting
[10]:
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. Observations: 5000
------------------ 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 learnerml_m
propensity Score learnercontrol_group
the control group for the parallel trend assumptiongt_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.
[11]:
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.
[12]:
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. Observations: 5000
------------------ 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.94836173 1.91295948 1.88367726 2.84355234 3.91845868 1.92203792
1.94989465 1.90964501 1.96929681 2.94232565 1.96761991 1.92005734
1.90752777 1.96513263 1.92848882]]
Learner ml_g1 RMSE: [[1.98970561 1.94187627 1.90828267 2.77429828 3.95426764 1.92668876
1.92162779 1.95755607 1.96497636 2.98936692 1.94272935 1.93162389
1.87872611 2.000074 1.90798703]]
Classification:
Learner ml_m Log Loss: [[0.67644408 0.67995509 0.67906511 0.67482084 0.67668477 0.70107629
0.70858079 0.70217582 0.71041569 0.70018404 0.73005905 0.72121977
0.72502083 0.72985873 0.73199646]]
------------------ 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.024565 0.136528 -0.179928 8.572093e-01
ATT(2025-04,2025-02,2025-03) -0.114599 0.152428 -0.751820 4.521595e-01
ATT(2025-04,2025-03,2025-04) 0.787095 0.136314 5.774119 7.735665e-09
ATT(2025-04,2025-03,2025-05) 1.920021 0.224070 8.568832 0.000000e+00
ATT(2025-04,2025-03,2025-06) 2.766289 0.297954 9.284286 0.000000e+00
ATT(2025-05,2025-01,2025-02) -0.009058 0.112638 -0.080417 9.359055e-01
ATT(2025-05,2025-02,2025-03) -0.159359 0.110496 -1.442211 1.492429e-01
ATT(2025-05,2025-03,2025-04) 0.035412 0.108987 0.324921 7.452406e-01
ATT(2025-05,2025-04,2025-05) 0.943754 0.103478 9.120325 0.000000e+00
ATT(2025-05,2025-04,2025-06) 1.814022 0.160667 11.290580 0.000000e+00
ATT(2025-06,2025-01,2025-02) -0.024554 0.089369 -0.274755 7.835048e-01
ATT(2025-06,2025-02,2025-03) 0.039897 0.084496 0.472175 6.368022e-01
ATT(2025-06,2025-03,2025-04) -0.058614 0.084248 -0.695731 4.865973e-01
ATT(2025-06,2025-04,2025-05) 0.100747 0.092205 1.092638 2.745527e-01
ATT(2025-06,2025-05,2025-06) 0.921136 0.098476 9.353918 0.000000e+00
2.5 % 97.5 %
ATT(2025-04,2025-01,2025-02) -0.292156 0.243025
ATT(2025-04,2025-02,2025-03) -0.413352 0.184155
ATT(2025-04,2025-03,2025-04) 0.519924 1.054266
ATT(2025-04,2025-03,2025-05) 1.480851 2.359191
ATT(2025-04,2025-03,2025-06) 2.182310 3.350267
ATT(2025-05,2025-01,2025-02) -0.229824 0.211708
ATT(2025-05,2025-02,2025-03) -0.375928 0.057210
ATT(2025-05,2025-03,2025-04) -0.178198 0.249022
ATT(2025-05,2025-04,2025-05) 0.740941 1.146567
ATT(2025-05,2025-04,2025-06) 1.499121 2.128923
ATT(2025-06,2025-01,2025-02) -0.199714 0.150605
ATT(2025-06,2025-02,2025-03) -0.125712 0.205505
ATT(2025-06,2025-03,2025-04) -0.223736 0.106509
ATT(2025-06,2025-04,2025-05) -0.079972 0.281466
ATT(2025-06,2025-05,2025-06) 0.728127 1.114145
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.
[13]:
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
[13]:
2.5 % | 97.5 % | |
---|---|---|
ATT(2025-04,2025-01,2025-02) | -0.417241 | 0.368110 |
ATT(2025-04,2025-02,2025-03) | -0.553005 | 0.323808 |
ATT(2025-04,2025-03,2025-04) | 0.395035 | 1.179155 |
ATT(2025-04,2025-03,2025-05) | 1.275561 | 2.564480 |
ATT(2025-04,2025-03,2025-06) | 1.909329 | 3.623248 |
ATT(2025-05,2025-01,2025-02) | -0.333021 | 0.314905 |
ATT(2025-05,2025-02,2025-03) | -0.477164 | 0.158445 |
ATT(2025-05,2025-03,2025-04) | -0.278050 | 0.348874 |
ATT(2025-05,2025-04,2025-05) | 0.646136 | 1.241372 |
ATT(2025-05,2025-04,2025-06) | 1.351920 | 2.276123 |
ATT(2025-06,2025-01,2025-02) | -0.281592 | 0.232483 |
ATT(2025-06,2025-02,2025-03) | -0.203125 | 0.282919 |
ATT(2025-06,2025-03,2025-04) | -0.300923 | 0.183695 |
ATT(2025-06,2025-04,2025-05) | -0.164449 | 0.365943 |
ATT(2025-06,2025-05,2025-06) | 0.637905 | 1.204368 |
A visualization of the effects can be obtained via the plot_effects()
method.
Remark that the plot used joint confidence intervals per default.
[14]:
dml_obj.plot_effects()
/opt/hostedtoolcache/Python/3.12.10/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)
[14]:
(<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.
[15]:
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.340187 -0.124104 -0.024565 0.074973
ATT(2025-04,2025-02,2025-03) -0.465306 -0.217896 -0.114599 -0.011301
ATT(2025-04,2025-03,2025-04) 0.457959 0.689886 0.787095 0.884304
ATT(2025-04,2025-03,2025-05) 1.410744 1.762022 1.920021 2.078020
ATT(2025-04,2025-03,2025-06) 2.059645 2.549535 2.766289 2.983042
ATT(2025-05,2025-01,2025-02) -0.285830 -0.107138 -0.009058 0.089022
ATT(2025-05,2025-02,2025-03) -0.436663 -0.259404 -0.159359 -0.059314
ATT(2025-05,2025-03,2025-04) -0.258781 -0.080728 0.035412 0.151553
ATT(2025-05,2025-04,2025-05) 0.668328 0.835039 0.943754 1.052469
ATT(2025-05,2025-04,2025-06) 1.372394 1.635046 1.814022 1.992997
ATT(2025-06,2025-01,2025-02) -0.272968 -0.125462 -0.024554 0.076353
ATT(2025-06,2025-02,2025-03) -0.209257 -0.070581 0.039897 0.150374
ATT(2025-06,2025-03,2025-04) -0.301825 -0.163334 -0.058614 0.046107
ATT(2025-06,2025-04,2025-05) -0.156178 -0.005306 0.100747 0.206800
ATT(2025-06,2025-05,2025-06) 0.663870 0.825694 0.921136 1.016578
CI upper
ATT(2025-04,2025-01,2025-02) 0.312036
ATT(2025-04,2025-02,2025-03) 0.246146
ATT(2025-04,2025-03,2025-04) 1.104487
ATT(2025-04,2025-03,2025-05) 2.466576
ATT(2025-04,2025-03,2025-06) 3.478889
ATT(2025-05,2025-01,2025-02) 0.283278
ATT(2025-05,2025-02,2025-03) 0.128602
ATT(2025-05,2025-03,2025-04) 0.332497
ATT(2025-05,2025-04,2025-05) 1.227260
ATT(2025-05,2025-04,2025-06) 2.259770
ATT(2025-06,2025-01,2025-02) 0.223608
ATT(2025-06,2025-02,2025-03) 0.290097
ATT(2025-06,2025-03,2025-04) 0.185177
ATT(2025-06,2025-04,2025-05) 0.359726
ATT(2025-06,2025-05,2025-06) 1.179324
------------------ Robustness Values ------------------
H_0 RV (%) RVa (%)
ATT(2025-04,2025-01,2025-02) 0.0 0.749067 0.000625
ATT(2025-04,2025-02,2025-03) 0.0 3.322724 0.000421
ATT(2025-04,2025-03,2025-04) 0.0 21.808972 14.322950
ATT(2025-04,2025-03,2025-05) 0.0 30.793505 25.865657
ATT(2025-04,2025-03,2025-06) 0.0 32.045896 25.312285
ATT(2025-05,2025-01,2025-02) 0.0 0.280761 0.000592
ATT(2025-05,2025-02,2025-03) 0.0 4.735707 0.000573
ATT(2025-05,2025-03,2025-04) 0.0 0.924600 0.000615
ATT(2025-05,2025-04,2025-05) 0.0 23.176854 19.425637
ATT(2025-05,2025-04,2025-06) 0.0 26.473291 22.932973
ATT(2025-06,2025-01,2025-02) 0.0 0.738625 0.000575
ATT(2025-06,2025-02,2025-03) 0.0 1.094119 0.000519
ATT(2025-06,2025-03,2025-04) 0.0 1.690558 0.000522
ATT(2025-06,2025-04,2025-05) 0.0 2.852153 0.000665
ATT(2025-06,2025-05,2025-06) 0.0 25.392824 21.109726
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).
[16]:
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.10/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)
[16]:
(<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
[17]:
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.10/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)
[17]:
(<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}\).
[18]:
df_post_treatment = df[df["t"] >= df["d"]]
df_post_treatment.groupby("d")["ite"].mean()
[18]:
d
2025-04-01 1.996048
2025-05-01 1.477051
2025-06-01 1.005521
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.
[19]:
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.388452 0.111107 12.496551 0.0 1.170687 1.606217
------------------ Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
2025-04 1.824468 0.200263 9.110374 0.0 1.431960 2.216976
2025-05 1.378888 0.119882 11.502023 0.0 1.143923 1.613853
2025-06 0.921136 0.098476 9.353918 0.0 0.728127 1.114145
------------------ Additional Information ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0
/home/runner/work/doubleml-docs/doubleml-docs/doubleml-for-py/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
[20]:
df_post_treatment.groupby("t")["ite"].mean()
[20]:
t
2025-04-01 1.006032
2025-05-01 1.512341
2025-06-01 2.017277
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"
.
[21]:
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.367427 0.126067 10.846837 0.0 1.12034 1.614513
------------------ Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
2025-04 0.787095 0.136314 5.774119 7.735665e-09 0.519924 1.054266
2025-05 1.453144 0.142385 10.205737 0.000000e+00 1.174075 1.732214
2025-06 1.862041 0.155875 11.945767 0.000000e+00 1.556533 2.167550
------------------ Additional Information ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0
/home/runner/work/doubleml-docs/doubleml-docs/doubleml-for-py/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
[22]:
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:]
[22]:
e
-122 days -0.016965
-92 days -0.015412
-61 days 0.049267
-31 days 0.016838
0 days 0.991768
31 days 2.004498
59 days 2.965808
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).
[23]:
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 %
1.838909 0.168093 10.939805 0.0 1.509452 2.168366
------------------ Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
-4 months -0.024554 0.089369 -0.274755 0.783505 -0.199714 0.150605
-3 months 0.015389 0.073491 0.209398 0.834137 -0.128651 0.159429
-2 months -0.079205 0.075151 -1.053941 0.291910 -0.226498 0.068089
-1 months 0.003521 0.085479 0.041192 0.967143 -0.164014 0.171056
0 months 0.881108 0.074469 11.831955 0.000000 0.735153 1.027064
1 months 1.869329 0.165298 11.308817 0.000000 1.545350 2.193308
2 months 2.766289 0.297954 9.284286 0.000000 2.182310 3.350267
------------------ Additional Information ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0
/home/runner/work/doubleml-docs/doubleml-docs/doubleml-for-py/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(
[23]:
(<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 aggregationaggregation_weights
: weights for the aggregation
To clarify, e.g. for the eventstudy aggregation
[24]:
print(aggregated_eventstudy)
================== DoubleMLDIDAggregation Object ==================
Event Study Aggregation
------------------ Overall Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
1.838909 0.168093 10.939805 0.0 1.509452 2.168366
------------------ Aggregated Effects ------------------
coef std err t P>|t| 2.5 % 97.5 %
-4 months -0.024554 0.089369 -0.274755 0.783505 -0.199714 0.150605
-3 months 0.015389 0.073491 0.209398 0.834137 -0.128651 0.159429
-2 months -0.079205 0.075151 -1.053941 0.291910 -0.226498 0.068089
-1 months 0.003521 0.085479 0.041192 0.967143 -0.164014 0.171056
0 months 0.881108 0.074469 11.831955 0.000000 0.735153 1.027064
1 months 1.869329 0.165298 11.308817 0.000000 1.545350 2.193308
2 months 2.766289 0.297954 9.284286 0.000000 2.182310 3.350267
------------------ Additional Information ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0
Here, the overall effect aggregation aggregates each effect with positive exposure
[25]:
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
[26]:
# the weights for e=0 correspond to the fifth element of the aggregation weights
aggregated_eventstudy.aggregation_weights[4]
[26]:
array([0. , 0. , 0.3532565 , 0. , 0. ,
0. , 0. , 0. , 0.32377379, 0. ,
0. , 0. , 0. , 0. , 0.32296971])
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)\)
[27]:
print(dml_obj.summary["coef"])
ATT(2025-04,2025-01,2025-02) -0.024565
ATT(2025-04,2025-02,2025-03) -0.114599
ATT(2025-04,2025-03,2025-04) 0.787095
ATT(2025-04,2025-03,2025-05) 1.920021
ATT(2025-04,2025-03,2025-06) 2.766289
ATT(2025-05,2025-01,2025-02) -0.009058
ATT(2025-05,2025-02,2025-03) -0.159359
ATT(2025-05,2025-03,2025-04) 0.035412
ATT(2025-05,2025-04,2025-05) 0.943754
ATT(2025-05,2025-04,2025-06) 1.814022
ATT(2025-06,2025-01,2025-02) -0.024554
ATT(2025-06,2025-02,2025-03) 0.039897
ATT(2025-06,2025-03,2025-04) -0.058614
ATT(2025-06,2025-04,2025-05) 0.100747
ATT(2025-06,2025-05,2025-06) 0.921136
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.
[28]:
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()
(19326, 10)
[28]:
id | y | y0 | y1 | d | t | Z1 | Z2 | Z3 | Z4 | |
---|---|---|---|---|---|---|---|---|---|---|
1 | 0 | 219.562183 | 219.562183 | 218.415852 | 2025-05-01 | 2025-01-01 | 0.320162 | 1.418377 | 0.282289 | -0.008114 |
2 | 0 | 222.821309 | 222.821309 | 223.645736 | 2025-05-01 | 2025-02-01 | 0.320162 | 1.418377 | 0.282289 | -0.008114 |
3 | 0 | 228.535420 | 228.535420 | 226.788844 | 2025-05-01 | 2025-03-01 | 0.320162 | 1.418377 | 0.282289 | -0.008114 |
4 | 0 | 234.835139 | 232.698282 | 234.835139 | 2025-05-01 | 2025-04-01 | 0.320162 | 1.418377 | 0.282289 | -0.008114 |
5 | 0 | 236.435587 | 236.734813 | 236.435587 | 2025-05-01 | 2025-05-01 | 0.320162 | 1.418377 | 0.282289 | -0.008114 |
To visualize the anticipation, we will again plot the “oracle” values
[29]:
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()
[29]:
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.538922 | 192.116895 | 225.549934 | -0.018434 | -2.227808 | 2.312673 |
1 | 2025-01-01 | 2025-05 | 210.496990 | 193.994532 | 226.893456 | -0.066004 | -2.366184 | 2.395901 |
2 | 2025-01-01 | 2025-06 | 213.620343 | 195.707889 | 230.905193 | -0.041665 | -2.387054 | 2.436900 |
3 | 2025-01-01 | Never Treated | 217.498480 | 200.147312 | 235.095746 | -0.016816 | -2.457736 | 2.211449 |
4 | 2025-02-01 | 2025-04 | 208.322516 | 183.404305 | 234.093308 | -0.026305 | -2.479926 | 2.315162 |
One can see that the effect is already anticipated one period before the actual treatment assignment.
[30]:
plot_data(agg_df_anticipation, col_name='ite')

Initialize a corresponding DoubleMLPanelData
object.
[31]:
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.
[32]:
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.10/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)
[32]:
(<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.
[33]:
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.10/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.10/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)
[33]:
(<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).
[34]:
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()
/home/runner/work/doubleml-docs/doubleml-docs/doubleml-for-py/doubleml/double_ml.py:1470: UserWarning: The estimated nu2 for d is not positive. Re-estimation based on riesz representer (non-orthogonal).
warnings.warn(msg, UserWarning)
[34]:
(<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
[35]:
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()
[35]:
(<Figure size 1200x800 with 2 Axes>,
[<Axes: title={'center': 'First Treated: 2025-04'}, xlabel='Evaluation Period', ylabel='Effect'>])
