Difference-in-differences, and the trend you have to trust

Lumen rolled a loyalty program out to some regions before others. You cannot rerun history to see what the early-adopter regions would have earned without loyalty - the counterfactual is missing, as always. But a plain before/after in the treated regions is confounded by everything else that changed over that time (seasonality, a price move, macro). Difference-in-differences fixes this by subtracting the control group's change from the treated group's change: whatever common shock hit both cancels out, leaving the treatment effect. The whole method rests on one assumption - parallel trends - and this session teaches you to build DiD and to interrogate that assumption in code.

With one treated group, one control group, and two periods (pre and post), DiD is a difference of two differences. Take the treated group's post-minus-pre change, subtract the control group's post-minus-pre change, and what remains is the effect - because any shock common to both groups sits in both differences and cancels:

DiD = (ȲT,post - ȲT,pre) - (ȲC,post - ȲC,pre)

DiD is only causal if parallel trends holds: absent the treatment, the treated and control groups would have moved by the same amount over time. This is untestable directly - it is a claim about the missing counterfactual - but you can build confidence by checking that the groups moved in parallel before treatment, using an event study. Scaling from 2x2 to many units and periods, you add unit and time fixed effects: two-way fixed effects (TWFE).

Start with the four-means version so the estimate is never a black box. Take Lumen's early-adopter regions as treated and the late-adopter regions as control, split weekly revenue into pre- and post-rollout, and compute the double difference.

import pandas as pd, numpy as np

df = pd.read_parquet("lumen_loyalty_panel.parquet")   # region, week, revenue
early = ["US-CA", "US-NY", "US-TX"]                    # got loyalty first
ROLLOUT = pd.Timestamp("2026-04-06")

df["Treat"] = df.region.isin(early).astype(int)
df["Post"]  = (df.week >= ROLLOUT).astype(int)

m = df.groupby(["Treat", "Post"])["revenue"].mean()
T_pre,  T_post = m[(1, 0)], m[(1, 1)]
C_pre,  C_post = m[(0, 0)], m[(0, 1)]

did = (T_post - T_pre) - (C_post - C_pre)
print(f"treated change  = {T_post - T_pre:,.0f}")
print(f"control change  = {C_post - C_pre:,.0f}")
print(f"DiD estimate    = {did:,.0f}")   # the effect, common shocks cancelled

Reproduce the hand number as a regression coefficient, then scale to a full panel with region and week fixed effects. The Treat×Post interaction is your estimate; TWFE absorbs every fixed region gap and every common weekly shock.

import statsmodels.formula.api as smf

ols = smf.ols("revenue ~ Treat * Post", data=df).fit(
        cov_type="cluster", cov_kwds={"groups": df["region"]})
print(ols.params["Treat:Post"])          # == the by-hand DiD
print(ols.summary())                     # SE + CI on the interaction term

from linearmodels.panel import PanelOLS

panel = df.set_index(["region", "week"])
panel["TreatPost"] = panel["Treat"] * panel["Post"]

twfe = PanelOLS.from_formula(
    "revenue ~ TreatPost + EntityEffects + TimeEffects",
    data=panel).fit(cov_type="clustered", cluster_entity=True)
print(twfe.params["TreatPost"])          # DiD with region + week fixed effects

Regress on event-time leads and lags relative to each region's rollout week. Pre-period (lead) coefficients near zero support parallel trends; post-period (lag) coefficients trace the dynamic effect. Then plot it with confidence bands.

# weeks relative to each region's own rollout; controls never treated -> large event time
df["evt"] = ((df.week - df.rollout_week) / np.timedelta64(1, "W")).round().astype(int)
df["evt"] = df["evt"].clip(-4, 6)                 # bin the tails

es = smf.ols("revenue ~ C(evt, Treatment(reference=-1)) + C(region) + C(week)",
             data=df).fit(cov_type="cluster", cov_kwds={"groups": df["region"]})

import matplotlib.pyplot as plt
coef = es.params.filter(like="C(evt")
ci   = es.conf_int().loc[coef.index]
k    = [int(s.split("[T.")[1].rstrip("]")) for s in coef.index]

plt.errorbar(k, coef.values,
             yerr=[coef.values - ci[0], ci[1] - coef.values], fmt="o")
plt.axhline(0, ls="--"); plt.axvline(-0.5, color="red", ls="--")
plt.xlabel("weeks since rollout"); plt.ylabel("coefficient")
# leads (k < 0) near zero -> parallel trends holds; lags (k >= 0) -> the effect

• Break parallel trends on purpose. Add a pre-existing upward trend to the treated regions only, rerun DiD, and watch it report a "loyalty effect" that is really just the trend. Confirm the event-study leads now slope - the tell.
• Reconcile the three numbers. Show the by-hand 2x2, the interaction β₃, and TWFE agree on a clean balanced panel, and explain what changes when the panel is unbalanced.
• Try pyfixest. Re-fit the TWFE and event study in pyfixest and compare speed and syntax to linearmodels - it is becoming the default for large panels.
• Optional: read the Goodman-Bacon decomposition of your TWFE estimate and identify which 2x2 comparisons carry negative weight under the staggered rollout.

We build the working core of DiD - the estimator, the assumption, and the honest test - and flag the modern staggered-adoption literature. The full econometric treatment lives in the sources below.