Difference-in-differences in Python: didregress and parallel trends

A policy starts in some states and not in others. Comparing treated and untreated states after the policy confuses its effect with everything else that differs between them. Comparing the treated states before and after confuses it with everything else that changed over time. Difference in differences (DiD) combines the two comparisons: the change in the treated groups minus the change in the control groups. Fixed differences between groups and shocks common to all groups cancel out.

What difference in differences estimates

In the simplest case there are two groups and two periods. The estimate is a difference of four means:

δ̂ = (ȳ_treated,after − ȳ_treated,before) − (ȳ_control,after − ȳ_control,before)

With many groups and periods the same idea becomes a two-way fixed-effects (TWFE) regression, y_gt = a_g + c_t + δ D_gt + x_gt′b + e_gt. The group effects a_g absorb permanent differences between groups, the time effects c_t absorb shocks that hit every group, and D_gt is 1 for a treated group in a treated period. The coefficient δ is the average treatment effect on the treated (ATET).

The key assumption is parallel trends: without the policy, the treated groups' average outcome would have moved in parallel with the controls'. Levels may differ; trends may not. A second assumption is no anticipation: outcomes do not react before the policy starts. Neither can be verified. Pre-treatment data can reveal diverging trends, but a test that does not reject is not proof that trends are parallel. Roth, Sant'Anna, Bilinski and Poe (2023) review the assumptions and the recent extensions.

Inference needs care. Policies are usually assigned to whole groups, and outcomes within a group are correlated over time. Bertrand, Duflo and Mullainathan (2004) show that standard errors that ignore this serial correlation are far too small and reject true null hypotheses much too often. Clustering by group is the standard remedy, and it needs enough groups.

The single TWFE coefficient is reliable when every treated group starts at the same date. With staggered adoption and effects that differ across cohorts or change over time, δ is a weighted average of many two-group comparisons. Some of them use already-treated groups as controls, and the resulting weights on group-time effects can be negative (Goodman-Bacon 2021). The staggered DiD page covers estimators built for that case.

In OpenEconometrics

oe.didregress fits the TWFE model. The outcome goes in y, the 0/1 treated-and-post indicator in treatment, the group and period identifiers in group and time, and time-varying controls in x. Both sets of fixed effects are absorbed by alternating projections, so no dummy matrix is built. The same call handles a panel or repeated cross-sections with several observations per group and period. The example below corresponds to the Stata command didregress (y income) (policy), group(state) time(year); xtdidregress is the panel form of the same workflow.

• covariance="robust" (the default) clusters at the group level, as didregress does by default. "cluster" with cluster= clusters on another column, such as a region that contains several states. "HC1" and "nonrobust" do not cluster.
• Clustered fits use the small-sample factor G/(G − 1) · (N − 1)/(N − K) and t inference with G − 1 degrees of freedom.
• result.tests["parallel_trends"] corresponds to estat ptrends: an F test of a linear trend specific to the treated groups before treatment. result.tests["granger"] corresponds to estat granger: a joint F test of the treatment leads, that is, of anticipation effects.
• weights= with weight_type="aweight", "fweight" or "pweight" weights the regression.

The adoption pattern (common, staggered or not absorbing), the treated and control group counts and the first treated periods are in result.extra. The treatment-effects family guide gives the formulas, degrees-of-freedom rules and test definitions.

Supported options

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 50 states observed from 2008 to 2017. Twenty states adopt a policy in 2013 and keep it. Treated states have a higher outcome level to begin with, which the state effects absorb. All states share a trend of 0.3 a year, a control income enters with coefficient 0.5, and the errors follow an AR(1) process with coefficient 0.7 within each state. The true ATET is 1.5. The script then refits the model with unclustered HC1 standard errors to show what ignoring the serial correlation does.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(n_states=50, n_treated=20, seed=11):
    """State-year panel, 2008-2017; 20 states adopt the policy in 2013. True ATET = 1.5."""
    rng = np.random.default_rng(seed)
    rows = []
    for s in range(n_states):
        treated = s < n_treated
        state_effect = rng.normal(scale=2.0) + (1.0 if treated else 0.0)  # treated states differ in levels
        e = 0.0
        for t in range(2008, 2018):
            e = 0.7 * e + rng.normal()                  # AR(1) errors within each state
            income = 0.1 * (t - 2008) + rng.normal()
            d = int(treated and t >= 2013)
            y = state_effect + 0.3 * (t - 2008) + 0.5 * income + 1.5 * d + e
            rows.append((s, t, d, income, y))
    return pd.DataFrame(rows, columns=["state", "year", "policy", "income", "y"])


if __name__ == "__main__":
    df = simulate()
    did = oe.didregress(data=df, y="y", treatment="policy", group="state", time="year", x=["income"])
    print(did.summary())
    hc1 = oe.didregress(data=df, y="y", treatment="policy", group="state", time="year", x=["income"],
                        covariance="HC1").coefficients[0]
    print(f"Unclustered HC1: std. error {hc1.std_error:.4f}, 95% CI [{hc1.ci_low:.3f}, {hc1.ci_high:.3f}]")

With openecon 0.3.20a1 from PyPI the script prints:

Difference-in-differences regression — y
Observations: 500  |  Covariance: robust  |  Confidence: 95%

Term               Estimate  Std. error        t     P>|stat|  CI lower  CI upper
-----------------  --------  ----------  -------  -----------  --------  --------
[ATET]
ATET:r1vs0.policy   1.08076    0.332948  3.24602    0.0021132  0.411672   1.74984
[Controls]
income             0.500826   0.0527432  9.49555  1.09053e-12  0.394834  0.606817

r_squared: 0.919216  |  r_squared_within: 0.21879  |  rmse: 1.02114  |  df_resid: 488
n_groups: 50  |  n_periods: 10  |  n_clusters: 50
Parallel-trends test: treated-group linear trend before treatment = 0 (estat ptrends): F(1, 49) = 0.882372, p = 0.3522
Anticipation (Granger) test: treatment leads = 0 (estat granger): F(4, 49) = 0.563438, p = 0.6903
Unclustered HC1: std. error 0.1977, 95% CI [0.692, 1.469]

Reading the output

• The ATET is imprecise but covers the truth. The estimate is 1.081 against a true 1.5, with a 95% confidence interval from 0.412 to 1.750. With only 20 treated states and persistent state-level shocks, a gap of this size is ordinary sampling error.
• Clustering matters. The group-clustered standard error is 0.333; the unclustered HC1 one is 0.198, about 40% smaller. In this sample the narrower HC1 interval, 0.692 to 1.469, misses the true value. That is one draw, not a proof, but it is the pattern Bertrand, Duflo and Mullainathan (2004) document: with serially correlated errors, unclustered intervals cover the truth far less often than their nominal level.
• Degrees of freedom come from the clusters. The tests use 49 denominator degrees of freedom, one fewer than the 50 states, not the 488 residual degrees of freedom.
• The control is estimated well. The income coefficient is 0.5008 against a true 0.5.
• The pre-treatment checks do not reject. The parallel-trends test gives F(1, 49) = 0.88 with p = 0.35, and the joint test of four treatment leads gives F(4, 49) = 0.56 with p = 0.69. Trends hold by construction here. In real data such results only say that the pre-period shows no clear divergence; the tests may have little power.

Limitations

• The common oe.predict and oe.margins adapters do not cover oe.didregress; use its own result tables and family helpers.
• Cluster-robust covariance for oe.didregress allows at most 2 clustering dimensions.
• The parallel-trends and Granger tests are computed only when every treated group starts in the same period, treatment is absorbing, never-treated controls exist and there are at least two pre-treatment periods. Otherwise they are omitted and extra["tests_note"] says why; use oe.eventstudy for pre-trends.
• With more than one adoption cohort the fit adds a warning: the single TWFE coefficient can be a negatively weighted average of heterogeneous effects. Use oe.csdid in that case.
• Stata's vce(hc2) with Bell–McCaffrey degrees of freedom and the wild cluster bootstrap for didregress are not implemented.
• If the treatment indicator has no variation left after the group and time effects are absorbed, for example because every group is treated at the same date, the fit raises no_within_variation.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Bertrand, M., Duflo, E. and Mullainathan, S. (2004). How much should we trust differences-in-differences estimates? Quarterly Journal of Economics, 119(1), 249–275. doi:10.1162/003355304772839588
• Goodman-Bacon, A. (2021). Difference-in-differences with variation in treatment timing. Journal of Econometrics, 225(2), 254–277. doi:10.1016/j.jeconom.2021.03.014
• Roth, J., Sant'Anna, P. H. C., Bilinski, A. and Poe, J. (2023). What's trending in difference-in-differences? A synthesis of the recent econometrics literature. Journal of Econometrics, 235(2), 2218–2244. doi:10.1016/j.jeconom.2023.03.008