Note

https://docs.doubleml.org/stable/examples/did/py_rep_cs.ipynb

Python: Repeated Cross-Sectional 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_cs_CS2021

Data#

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

We will observe approximately n_obs units over n_periods. The parameter lambda_t determines the probability of observing a unit i in time period t. The parameter lambda_t is set to 0.5 for all time periods, which means that each unit has a 50% chance of being observed in each time period.

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 = 10000
n_periods = 6

df = make_did_cs_CS2021(n_obs, dgp_type=4, include_never_treated=True, n_periods=n_periods, n_pre_treat_periods=3,
                        lambda_t=0.5, time_type="float")
df["ite"] = df["y1"] - df["y0"]

print(df.shape)
df.head()

(59997, 11)

[2]:

	id	y	y0	y1	d	t	Z1	Z2	Z3	Z4	ite
0	0	225.811296	225.811296	224.451202	5.0	0	2.420755	0.180447	0.567191	1.517293	-1.360094
8	1	222.552154	222.552154	223.610274	inf	2	1.067125	-0.673612	0.462135	-1.004924	1.058120
10	1	229.451191	229.451191	231.191981	inf	4	1.067125	-0.673612	0.462135	-1.004924	1.740790
11	1	234.775399	234.775399	233.427297	inf	5	1.067125	-0.673612	0.462135	-1.004924	-1.348102
13	2	235.625437	235.625437	234.333672	inf	1	1.178469	0.381408	0.681317	1.028896	-1.291765

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`).

The data set with repeated cross-sectional data is generated on the basis of a panel data set with the following data generating process (DGP). To obtain repeated cross-sectional data, the number of generated individuals is increased to \(\frac{n_{obs}}{\lambda_t}\), where \(\lambda_t\) denotes the probability to observe a unit at each time period (time constant).

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.

[3]:

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

[3]:

t
0    10045
1    10043
2    10029
3    10003
4    10019
5     9858
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
3.0    15162
4.0    14445
5.0    15068
inf    15322
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
3.0    15162
4.0    14445
5.0    15068
inf    15322
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

# 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", "d"]).agg(**agg_dictionary).reset_index()
agg_df.head()

[6]:

	t	d	y_mean	y_lower_quantile	y_upper_quantile	ite_mean	ite_lower_quantile	ite_upper_quantile
0	0	3.0	209.056308	198.508749	219.117143	-0.025680	-2.342504	2.335146
1	0	4.0	210.437395	200.276142	220.659240	0.002298	-2.277681	2.320602
2	0	5.0	212.360509	202.534866	222.366200	-0.023649	-2.384957	2.323247
3	0	inf	214.142785	204.051782	224.058405	-0.019320	-2.303132	2.245458
4	1	3.0	208.311441	187.639125	228.170627	-0.021297	-2.345081	2.373114

[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["d"].nunique()
    color_palette = sns.color_palette("colorblind", n_colors=n_colors)

    sns.lineplot(
        data=df,
        x='t',
        y=f'{col_name}_mean',
        hue='d',
        style='d',
        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='d', 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')

/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/matplotlib/colors.py:2560: RuntimeWarning: invalid value encountered in divide
  resdat /= (vmax - vmin)

../../_images/examples_did_py_rep_cs_16_1.png

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\]

[9]:

plot_data(agg_df, col_name='ite')

/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/matplotlib/colors.py:2560: RuntimeWarning: invalid value encountered in divide
  resdat /= (vmax - vmin)

../../_images/examples_did_py_rep_cs_18_1.png

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

[10]:

dml_data = DoubleMLPanelData(
    data=df,
    y_col="y",
    d_cols="d",
    id_col="id",
    t_col="t",
    x_cols=["Z1", "Z2", "Z3", "Z4"],
)
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
Static panel data: False
No. Unique Ids: 19697
No. Observations: 59997
------------------ DataFrame info    ------------------
<class 'pandas.DataFrame'>
RangeIndex: 59997 entries, 0 to 59996
Columns: 11 entries, id to ite
dtypes: float64(9), int64(2)
memory usage: 5.0 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.

For repeated cross-sectional data, we additionally specify the argument

panel=False

[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",
    "anticipation_periods": 0,
    "n_folds": 5,
    "n_rep": 1,
    "panel": False,
}

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
Static panel data: False
No. Unique Ids: 19697
No. Observations: 59997
------------------ 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_g_d0_t0 RMSE: [[1.77579806 2.50375945 3.29986152 3.29626809 3.28445381 1.79397993
  2.53000541 3.21751864 4.20608296 4.23743218 1.79752168 2.48530973
  3.26835588 4.283126   5.2131769 ]]
Learner ml_g_d0_t1 RMSE: [[2.51257668 3.28340222 4.2904653  5.28968936 6.22924973 2.5096209
  3.27316897 4.26912562 5.32767024 6.21020534 2.5179741  3.29874381
  4.3303128  5.38143205 6.2888387 ]]
Learner ml_g_d1_t0 RMSE: [[1.81828037 2.34589552 3.22469841 3.27245126 3.15798811 1.81249813
  2.45209869 3.35622929 4.2539725  4.35351645 1.76503842 2.54474312
  3.46558247 4.39139169 5.22256304]]
Learner ml_g_d1_t1 RMSE: [[2.38252719 3.19400688 4.17726451 4.8673076  6.04120741 2.47143505
  3.33715962 4.31700617 5.31910936 6.4321066  2.5088117  3.45505077
  4.29592703 5.25812449 6.31272156]]
Classification:
Learner ml_m Log Loss: [[0.62754317 0.62820192 0.62920744 0.62910079 0.62650023 0.65125757
  0.6543574  0.6571436  0.65450699 0.65596116 0.66864343 0.66910646
  0.66850173 0.66714126 0.66555039]]

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

------------------ Fit summary       ------------------
                  coef   std err          t         P>|t|     2.5 %    97.5 %
ATT(3.0,0,1) -0.043606  0.104116  -0.418827  6.753425e-01 -0.247670  0.160457
ATT(3.0,1,2)  0.168154  0.132632   1.267817  2.048635e-01 -0.091801  0.428108
ATT(3.0,2,3)  0.698672  0.193707   3.606847  3.099407e-04  0.319013  1.078331
ATT(3.0,2,4)  1.680938  0.197985   8.490220  0.000000e+00  1.292894  2.068982
ATT(3.0,2,5)  2.644529  0.236352  11.188954  0.000000e+00  2.181288  3.107770
ATT(4.0,0,1) -0.049775  0.096484  -0.515888  6.059330e-01 -0.238881  0.139331
ATT(4.0,1,2) -0.008764  0.128945  -0.067964  9.458145e-01 -0.261492  0.243964
ATT(4.0,2,3) -0.158548  0.163196  -0.971517  3.312908e-01 -0.478407  0.161311
ATT(4.0,3,4)  1.123582  0.214559   5.236695  1.634776e-07  0.703053  1.544111
ATT(4.0,3,5)  1.969160  0.253273   7.774841  7.549517e-15  1.472754  2.465567
ATT(5.0,0,1)  0.043189  0.090061   0.479551  6.315469e-01 -0.133327  0.219705
ATT(5.0,1,2)  0.028818  0.122828   0.234617  8.145064e-01 -0.211921  0.269556
ATT(5.0,2,3) -0.013297  0.157399  -0.084479  9.326758e-01 -0.321793  0.295199
ATT(5.0,3,4) -0.118488  0.207658  -0.570591  5.682770e-01 -0.525491  0.288515
ATT(5.0,4,5)  0.898225  0.236700   3.794780  1.477742e-04  0.434301  1.362149

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=3, t_{\text{pre}}=0, t_{\text{eval}}=1)\) which estimates the pre-trend from time period \(0\) to \(1\) even if the actual treatment occured in time period \(3\).

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(3.0,0,1)	-0.346921	0.259708
ATT(3.0,1,2)	-0.218236	0.554543
ATT(3.0,2,3)	0.134357	1.262987
ATT(3.0,2,4)	1.104160	2.257717
ATT(3.0,2,5)	1.955980	3.333078
ATT(4.0,0,1)	-0.330857	0.231307
ATT(4.0,1,2)	-0.384412	0.366885
ATT(4.0,2,3)	-0.633978	0.316882
ATT(4.0,3,4)	0.498519	1.748645
ATT(4.0,3,5)	1.231314	2.707006
ATT(5.0,0,1)	-0.219180	0.305557
ATT(5.0,1,2)	-0.329010	0.386646
ATT(5.0,2,3)	-0.471837	0.445243
ATT(5.0,3,4)	-0.723446	0.486471
ATT(5.0,4,5)	0.208661	1.587790

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()

[14]:

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 3.0'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 4.0'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 5.0'}, xlabel='Evaluation Period', ylabel='Effect'>])

../../_images/examples_did_py_rep_cs_30_1.png

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  CI upper
ATT(3.0,0,1) -0.518764    -0.346434 -0.043606     0.259221  0.430225
ATT(3.0,1,2) -0.460032    -0.241805  0.168154     0.578112  0.797115
ATT(3.0,2,3) -0.150474     0.169095  0.698672     1.228248  1.547137
ATT(3.0,2,4)  0.750094     1.076510  1.680938     2.285367  2.611635
ATT(3.0,2,5)  1.556490     1.945446  2.644529     3.343611  3.733957
ATT(4.0,0,1) -0.501544    -0.342700 -0.049775     0.243150  0.402147
ATT(4.0,1,2) -0.612357    -0.400319 -0.008764     0.382792  0.595605
ATT(4.0,2,3) -0.931519    -0.663052 -0.158548     0.345957  0.615258
ATT(4.0,3,4)  0.131558     0.484254  1.123582     1.762910  2.117247
ATT(4.0,3,5)  0.843258     1.261746  1.969160     2.676574  3.092661
ATT(5.0,0,1) -0.385454    -0.237183  0.043189     0.323561  0.471936
ATT(5.0,1,2) -0.554128    -0.351569  0.028818     0.409205  0.611312
ATT(5.0,2,3) -0.772644    -0.513396 -0.013297     0.486802  0.746199
ATT(5.0,3,4) -1.086428    -0.744677 -0.118488     0.507701  0.850143
ATT(5.0,4,5) -0.245297     0.143918  0.898225     1.652532  2.043338

------------------ Robustness Values ------------------
              H_0     RV (%)   RVa (%)
ATT(3.0,0,1)  0.0   0.437498  0.000620
ATT(3.0,1,2)  0.0   1.241747  0.000406
ATT(3.0,2,3)  0.0   3.938757  2.159247
ATT(3.0,2,4)  0.0   8.119944  6.583931
ATT(3.0,2,5)  0.0  10.878060  9.330447
ATT(4.0,0,1)  0.0   0.516422  0.000625
ATT(4.0,1,2)  0.0   0.068009  0.000434
ATT(4.0,2,3)  0.0   0.952834  0.000463
ATT(4.0,3,4)  0.0   5.211897  3.605618
ATT(4.0,3,5)  0.0   8.127121  6.441712
ATT(5.0,0,1)  0.0   0.467949  0.000617
ATT(5.0,1,2)  0.0   0.230345  0.000538
ATT(5.0,2,3)  0.0   0.080814  0.000484
ATT(5.0,3,4)  0.0   0.574877  0.000648
ATT(5.0,4,5)  0.0   3.562075  2.035034

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()

[16]:

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 3.0'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 4.0'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 5.0'}, xlabel='Evaluation Period', ylabel='Effect'>])

../../_images/examples_did_py_rep_cs_35_1.png

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.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(
/opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/sklearn/linear_model/_logistic.py:451: OptimizeWarning: Unknown solver options: iprint
  opt_res = optimize.minimize(

[17]:

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 3.0'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 4.0'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 5.0'}, xlabel='Evaluation Period', ylabel='Effect'>])

../../_images/examples_did_py_rep_cs_37_2.png

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
3.0    1.906640
4.0    1.505688
5.0    1.025252
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.371321 0.119222 11.502225    0.0 1.13765 1.604993
------------------ Aggregated Effects         ------------------
         coef   std err          t         P>|t|     2.5 %    97.5 %
3.0  1.674713  0.147021  11.390952  0.000000e+00  1.386556  1.962870
4.0  1.546371  0.186516   8.290819  2.220446e-16  1.180806  1.911936
5.0  0.898225  0.236700   3.794780  1.477742e-04  0.434301  1.362149
------------------ Additional Information     ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0

/opt/hostedtoolcache/Python/3.12.13/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(

../../_images/examples_did_py_rep_cs_42_2.png

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
3    0.908554
4    1.499078
5    1.960356
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.314948 0.104874 12.538317    0.0 1.109398 1.520498
------------------ Aggregated Effects         ------------------
       coef   std err          t    P>|t|     2.5 %    97.5 %
3  0.698672  0.193707   3.606847  0.00031  0.319013  1.078331
4  1.409009  0.167348   8.419649  0.00000  1.081013  1.737004
5  1.837164  0.179978  10.207717  0.00000  1.484414  2.189915
------------------ Additional Information     ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0

/opt/hostedtoolcache/Python/3.12.13/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(

../../_images/examples_did_py_rep_cs_47_2.png

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_treated = df[df["d"] != np.inf].copy()
df_treated["e"] = df_treated["t"] - df_treated["d"]
df_treated.groupby("e")["ite"].mean().iloc[1:]

[22]:

e
-4.0    0.017538
-3.0   -0.013980
-2.0   -0.022272
-1.0    0.011879
 0.0    0.984695
 1.0    1.974179
 2.0    2.873069
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.789818 0.131107 13.651632    0.0 1.532854 2.046782
------------------ Aggregated Effects         ------------------
          coef   std err          t     P>|t|     2.5 %    97.5 %
-4.0  0.043189  0.090061   0.479551  0.631547 -0.133327  0.219705
-3.0 -0.009649  0.069499  -0.138840  0.889576 -0.145864  0.126566
-2.0 -0.022118  0.064403  -0.343425  0.731279 -0.148346  0.104110
-1.0 -0.034159  0.082357  -0.414767  0.678312 -0.195576  0.127258
0.0   0.903366  0.103107   8.761459  0.000000  0.701280  1.105451
1.0   1.821559  0.159842  11.396004  0.000000  1.508275  2.134844
2.0   2.644529  0.236352  11.188954  0.000000  2.181288  3.107770
------------------ Additional Information     ------------------
Score function: observational
Control group: never_treated
Anticipation periods: 0

/opt/hostedtoolcache/Python/3.12.13/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(

[23]:

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

../../_images/examples_did_py_rep_cs_51_3.png

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

[24]:

print(aggregated_eventstudy)

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

------------------ Overall Aggregated Effects ------------------
    coef  std err         t  P>|t|    2.5 %   97.5 %
1.789818 0.131107 13.651632    0.0 1.532854 2.046782
------------------ Aggregated Effects         ------------------
          coef   std err          t     P>|t|     2.5 %    97.5 %
-4.0  0.043189  0.090061   0.479551  0.631547 -0.133327  0.219705
-3.0 -0.009649  0.069499  -0.138840  0.889576 -0.145864  0.126566
-2.0 -0.022118  0.064403  -0.343425  0.731279 -0.148346  0.104110
-1.0 -0.034159  0.082357  -0.414767  0.678312 -0.195576  0.127258
0.0   0.903366  0.103107   8.761459  0.000000  0.701280  1.105451
1.0   1.821559  0.159842  11.396004  0.000000  1.508275  2.134844
2.0   2.644529  0.236352  11.188954  0.000000  2.181288  3.107770
------------------ 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.33938444, 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.3233352 , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.33728036])

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(3.0,0,1)   -0.043606
ATT(3.0,1,2)    0.168154
ATT(3.0,2,3)    0.698672
ATT(3.0,2,4)    1.680938
ATT(3.0,2,5)    2.644529
ATT(4.0,0,1)   -0.049775
ATT(4.0,1,2)   -0.008764
ATT(4.0,2,3)   -0.158548
ATT(4.0,3,4)    1.123582
ATT(4.0,3,5)    1.969160
ATT(5.0,0,1)    0.043189
ATT(5.0,1,2)    0.028818
ATT(5.0,2,3)   -0.013297
ATT(5.0,3,4)   -0.118488
ATT(5.0,4,5)    0.898225
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_cs_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()

(19137, 10)

[28]:

	id	y	y0	y1	d	t	Z1	Z2	Z3	Z4
2	0	227.691951	227.691951	226.453598	2025-06-01	2025-02-01	0.333050	-0.700270	0.507965	0.056144
3	0	231.231007	231.231007	232.442887	2025-06-01	2025-03-01	0.333050	-0.700270	0.507965	0.056144
4	0	238.174454	238.174454	237.948797	2025-06-01	2025-04-01	0.333050	-0.700270	0.507965	0.056144
5	0	245.418160	245.488350	245.418160	2025-06-01	2025-05-01	0.333050	-0.700270	0.507965	0.056144
9	1	198.076349	198.076349	198.233518	2025-05-01	2025-02-01	-1.038059	0.512974	-0.112859	0.090253

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

[29]:

df_anticipation["ite"] = df_anticipation["y1"] - df_anticipation["y0"]
agg_df_anticipation = df_anticipation.groupby(["t", "d"]).agg(**agg_dictionary).reset_index()
agg_df_anticipation.head()

[29]:

	t	d	y_mean	y_lower_quantile	y_upper_quantile	ite_mean	ite_lower_quantile	ite_upper_quantile
0	2025-01-01	2025-04-01	209.880281	193.322150	227.050682	0.022818	-2.188473	2.493069
1	2025-01-01	2025-05-01	210.721953	194.098509	227.148429	0.038161	-2.211419	2.353118
2	2025-01-01	2025-06-01	213.198377	197.476641	229.428866	0.034106	-2.254008	2.396582
3	2025-02-01	2025-04-01	209.711220	184.142451	235.766397	0.040827	-2.290011	2.397375
4	2025-02-01	2025-05-01	210.732985	185.330234	236.352557	0.002229	-2.292104	2.211288

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')

../../_images/examples_did_py_rep_cs_66_0.png

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()

[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'>])

../../_images/examples_did_py_rep_cs_70_1.png

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()

[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'>])

../../_images/examples_did_py_rep_cs_72_1.png

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()

[34]:

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 3.0'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 4.0'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 5.0'}, xlabel='Evaluation Period', ylabel='Effect'>])

../../_images/examples_did_py_rep_cs_75_1.png

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.

[35]:

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()

[35]:

(<Figure size 1200x800 with 4 Axes>,
 [<Axes: title={'center': 'First Treated: 3.0'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 4.0'}, ylabel='Effect'>,
  <Axes: title={'center': 'First Treated: 5.0'}, xlabel='Evaluation Period', ylabel='Effect'>])

../../_images/examples_did_py_rep_cs_77_1.png

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

[36]:

gt_dict = {
    "gt_combinations": [
        (4.0, 1, 2),
        (4.0, 1, 3),
        ]
}

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()

[36]:

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

../../_images/examples_did_py_rep_cs_79_1.png