import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
d=pd.read_csv("card_krueger.csv")
nj,pa=d[d.state==1],d[d.state==0]
print(f"Card-Krueger: {len(d)} fast-food stores ({len(nj)} NJ treated, {len(pa)} PA control), FTE employment before/after")
print(f"\nFirst stage (did the policy bite?): starting wage")
print(f"   NJ: ${nj.wage.mean():.2f} -> ${nj.wage2.mean():.2f}   (+${nj.wage2.mean()-nj.wage.mean():.2f})")
print(f"   PA: ${pa.wage.mean():.2f} -> ${pa.wage2.mean():.2f}   ({pa.wage2.mean()-pa.wage.mean():+.2f})")
njb,nja=nj.fte.mean(),nj.fte2.mean(); pab,paa=pa.fte.mean(),pa.fte2.mean()
did=(nja-njb)-(paa-pab)
print(f"\n2x2 FTE employment:")
print(f"                before   after   change")
print(f"   NJ (treated) {njb:6.2f}  {nja:6.2f}  {nja-njb:+6.2f}")
print(f"   PA (control) {pab:6.2f}  {paa:6.2f}  {paa-pab:+6.2f}")
print(f"   DiD = ({nja-njb:+.2f}) - ({paa-pab:+.2f}) = {did:+.2f}   <-- employment ROSE in NJ relative to PA")
fig,ax=plt.subplots(1,2,figsize=(13,4.6))
ax[0].bar(["NJ before","NJ after","PA before","PA after"],[nj.wage.mean(),nj.wage2.mean(),pa.wage.mean(),pa.wage2.mean()],color=[RED,RED,BLUE,BLUE],alpha=.8)
ax[0].axhline(4.25,color=GREY,ls=":",label="old min wage $4.25"); ax[0].axhline(5.05,color=GREEN,ls=":",label="new NJ min $5.05")
ax[0].set_ylabel("starting wage ($/hr)"); ax[0].set_ylim(4,5.3); ax[0].set_title("First stage: NJ wages jump, PA flat"); ax[0].legend(fontsize=8)
ax[1].plot([0,1],[njb,nja],"o-",color=RED,lw=2.5,ms=9,label="NJ (treated)")
ax[1].plot([0,1],[pab,paa],"s-",color=BLUE,lw=2.5,ms=9,label="PA (control)")
cf=njb+(paa-pab)                                  # NJ counterfactual under parallel trends
ax[1].plot([0,1],[njb,cf],"--",color=GREY,lw=2,label="NJ counterfactual (parallel)")
ax[1].annotate("",xy=(1,nja),xytext=(1,cf),arrowprops=dict(arrowstyle="<->",color=GREEN,lw=2))
ax[1].text(1.02,(nja+cf)/2,f"DiD\n{did:+.2f}",color=GREEN,fontsize=10)
ax[1].set_xticks([0,1]); ax[1].set_xticklabels(["before","after"]); ax[1].set_ylabel("FTE employment"); ax[1].set_title("Difference-in-differences: NJ vs its parallel-trends counterfactual"); ax[1].legend(fontsize=8,loc="lower left")
plt.tight_layout(); plt.show()
print(f"\nThe minimum wage rose in NJ, wages followed, and employment did NOT fall relative to PA (DiD {did:+.2f}).")
print("The second difference (PA) nets out the common 1992 downturn; the gap that remains is the causal estimate.")

Card-Krueger: 410 fast-food stores (331 NJ treated, 79 PA control), FTE employment before/after

First stage (did the policy bite?): starting wage
   NJ: $4.61 -> $5.08   (+$0.47)
   PA: $4.63 -> $4.62   (-0.01)

2x2 FTE employment:
                before   after   change
   NJ (treated)  20.44   21.03   +0.59
   PA (control)  23.33   21.17   -2.17
   DiD = (+0.59) - (-2.17) = +2.75   <-- employment ROSE in NJ relative to PA

The minimum wage rose in NJ, wages followed, and employment did NOT fall relative to PA (DiD +2.75).
The second difference (PA) nets out the common 1992 downturn; the gap that remains is the causal estimate.

# reshape to long and run the DiD regression with the interaction
long=pd.concat([
    pd.DataFrame({"store":d.sheet,"y":d.fte,"nj":d.state,"after":0}),
    pd.DataFrame({"store":d.sheet,"y":d.fte2,"nj":d.state,"after":1})]).dropna(subset=["y"]).reset_index(drop=True)
long["did"]=long.nj*long.after
X=np.column_stack([np.ones(len(long)),long.nj,long.after,long.did]); Y=long.y.values
b=np.linalg.lstsq(X,Y,rcond=None)[0]
# cluster-robust SE by store
res=Y-X@b; XtXinv=np.linalg.inv(X.T@X); meat=np.zeros((4,4))
for s,g in long.groupby("store"):
    ix=long.index.get_indexer(g.index); Xg=X[ix]; ug=res[ix]; meat+=Xg.T@np.outer(ug,ug)@Xg
V=XtXinv@meat@XtXinv; se=np.sqrt(np.diag(V))
print("DiD regression  Y = b0 + b1*NJ + b2*after + b3*(NJ x after):")
for nm,est,s in zip(["intercept","NJ","after","NJ x after (DiD)"],b,se):
    print(f"   {nm:20s} {est:+7.3f}  (cluster SE {s:.3f})")
print(f"\nThe interaction b3 = {b[3]:+.2f} reproduces the 2x2 DiD exactly, now with a clustered standard error.")
print("DiD = twoway fixed effects with one treatment dummy: group FE (NJ) + time FE (after) + treatment interaction.")

DiD regression  Y = b0 + b1*NJ + b2*after + b3*(NJ x after):
   intercept            +23.331  (cluster SE 1.342)
   NJ                    -2.892  (cluster SE 1.432)
   after                 -2.166  (cluster SE 1.214)
   NJ x after (DiD)      +2.754  (cluster SE 1.302)

The interaction b3 = +2.75 reproduces the 2x2 DiD exactly, now with a clustered standard error.
DiD = twoway fixed effects with one treatment dummy: group FE (NJ) + time FE (after) + treatment interaction.

sim=pd.read_csv("staggered_sim.csv")
trueATT=sim.loc[sim.treat==1,"eff"].mean()
# naive twoway FE (double-demean)
g2=sim.copy()
for c in ["y","treat"]:
    g2[c]=g2[c]-g2.groupby("id")[c].transform("mean"); g2[c]=g2[c]-g2.groupby("t")[c].transform("mean")
twfe=np.linalg.lstsq(g2[["treat"]].values,g2["y"].values,rcond=None)[0][0]
# Goodman-Bacon: enumerate 2x2 DiDs between timing groups, classify + weight
def bacon(df):
    groups=sorted(df.g.unique()); Y=df.pivot(index="id",columns="t",values="y"); gi=df.groupby("id")["g"].first()
    comps=[]
    treatable=[g for g in groups if g!=999]
    for k in treatable:                                    # k = treated group in this 2x2
        for l in groups:                                   # l = comparison group
            if l==k: continue
            if l==999:  kind,window="vs never-treated",(df.t.min(),df.t.max())          # clean
            elif l>k:   kind,window="vs later-treated (clean, pre)",(df.t.min(),l-1)     # clean (before l treated)
            else:       kind,window="vs earlier-treated (FORBIDDEN)",(l, df.t.max())     # forbidden: window opens when l was treated, so the control is ALREADY treated across the pre-period
            t0,t1=window
            if t1<=t0 or k<t0 or k>t1: continue
            tr=gi[gi==k].index; co=gi[gi==l].index
            pre=[t for t in Y.columns if t0<=t<k]; post=[t for t in Y.columns if k<=t<=t1]
            if not pre or not post: continue
            dd=(Y.loc[tr,post].mean().mean()-Y.loc[tr,pre].mean().mean())-(Y.loc[co,post].mean().mean()-Y.loc[co,pre].mean().mean())
            comps.append((k,l,kind,dd,len(tr)*len(co)))
    return pd.DataFrame(comps,columns=["treated_g","control_g","type","did2x2","weight"])
bc=bacon(sim); bc["w"]=bc.weight/bc.weight.sum()
print(f"true average ATT        = {trueATT:.3f}")
print(f"naive twoway FE         = {twfe:.3f}   <- biased DOWN by forbidden comparisons")
print(f"\nGoodman-Bacon decomposition (weight on each comparison type):")
print(bc.groupby("type").apply(lambda x:pd.Series({"avg 2x2 DiD":np.average(x.did2x2,weights=x.weight),"total weight":x.w.sum()})).round(3).to_string())
fig,ax=plt.subplots(figsize=(8.5,4))
cols={"vs never-treated":GREEN,"vs later-treated (clean, pre)":BLUE,"vs earlier-treated (FORBIDDEN)":RED}
for ty,gp in bc.groupby("type"): ax.scatter(gp.w,gp.did2x2,s=80,color=cols[ty],label=ty,alpha=.8,edgecolor="k")
ax.axhline(trueATT,color="k",ls="--",label=f"true ATT {trueATT:.1f}"); ax.set_xlabel("Goodman-Bacon weight"); ax.set_ylabel("2x2 DiD estimate")
ax.set_title("The forbidden comparisons (red) drag the TWFE estimate below the truth"); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()
print()
print(f"  weighted average over all 2x2 comparisons = {np.average(bc.did2x2,weights=bc.weight):.3f}")
print(f"  naive twoway FE                           = {twfe:.3f}")
print(f"  true ATT                                  = {trueATT:.3f}")
print()
print("Read the FORBIDDEN row first, because it carries the whole mechanism. Those comparisons use an")
print("ALREADY-TREATED group as the control, across a window in which that control's own effect is still")
print("growing. Their average 2x2 estimate sits far below the truth, and they carry real weight -- which")
print("is how a pooled regression built from individually sensible comparisons returns a biased answer.")
print()
print("TWFE is a weighted average of ALL of these. The decomposition here uses simple cell-count weights")
print("rather than Goodman-Bacon's variance weights, so it approximates the pooled estimate rather than")
print("reproducing it exactly; the two land close, and the direction and source of the bias are the point.")
print("This is the staggered-adoption trap.")

true average ATT        = 3.580
naive twoway FE         = 2.677   <- biased DOWN by forbidden comparisons

Goodman-Bacon decomposition (weight on each comparison type):

                                avg 2x2 DiD  total weight
type                                                     
vs earlier-treated (FORBIDDEN)        0.843         0.271
vs later-treated (clean, pre)         2.756         0.271
vs never-treated                      3.407         0.458

  weighted average over all 2x2 comparisons = 2.536
  naive twoway FE                           = 2.677
  true ATT                                  = 3.580

Read the FORBIDDEN row first, because it carries the whole mechanism. Those comparisons use an
ALREADY-TREATED group as the control, across a window in which that control's own effect is still
growing. Their average 2x2 estimate sits far below the truth, and they carry real weight -- which
is how a pooled regression built from individually sensible comparisons returns a biased answer.

TWFE is a weighted average of ALL of these. The decomposition here uses simple cell-count weights
rather than Goodman-Bacon's variance weights, so it approximates the pooled estimate rather than
reproducing it exactly; the two land close, and the direction and source of the bias are the point.
This is the staggered-adoption trap.

Y=sim.pivot(index="id",columns="t",values="y"); gi=sim.groupby("id")["g"].first()
groups=sorted([g for g in gi.unique() if g!=999]); times=sorted(sim.t.unique())
att=[]
for g in groups:
    for t in times:
        base=g-1
        if base not in Y.columns or t not in Y.columns or t==base: continue
        tr=gi[gi==g].index; co=gi[gi>t].index                # not-yet-treated (incl never) as clean controls
        if len(tr)==0 or len(co)==0: continue
        dd=(Y.loc[tr,t]-Y.loc[tr,base]).mean()-(Y.loc[co,t]-Y.loc[co,base]).mean()
        att.append((g,t,t-g,dd,len(tr)))
A=pd.DataFrame(att,columns=["g","t","e","att","n"])
overall=np.average(A.loc[A.e>=0,"att"],weights=A.loc[A.e>=0,"n"])
# event-study: average ATT by event-time e (relative to treatment)
es=A.groupby("e").apply(lambda x:np.average(x.att,weights=x.n))
print(f"true ATT                       = {trueATT:.3f}")
print(f"naive twoway FE                = {twfe:.3f}   (biased)")
print(f"Callaway-SantAnna (from scratch)= {overall:.3f}   <- recovers the truth")
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
names=["true ATT","naive TWFE","Callaway-\nSant'Anna"]; vals=[trueATT,twfe,overall]
ax[0].bar(names,vals,color=[GREY,RED,GREEN]); ax[0].axhline(trueATT,color="k",ls="--")
for i,v in enumerate(vals): ax[0].text(i,v+0.05,f"{v:.2f}",ha="center")
ax[0].set_ylabel("estimated ATT"); ax[0].set_title("CS recovers the truth; TWFE is biased")
ax[1].axhline(0,color="k",lw=.6); ax[1].axvline(-0.5,color=GREY,ls=":")
ax[1].plot(es.index,es.values,"o-",color=PURP,lw=2)
ax[1].scatter([e for e in es.index if e<0],[es[e] for e in es.index if e<0],color=BLUE,s=60,zorder=5,label="pre (should be ~0)")
ax[1].scatter([e for e in es.index if e>=0],[es[e] for e in es.index if e>=0],color=GREEN,s=60,zorder=5,label="post (dynamic effect)")
ax[1].set_xlabel("event time (periods since treatment)"); ax[1].set_ylabel("ATT"); ax[1].set_title("Event study: flat pre-trends, growing dynamic effect"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print("Pre-treatment event-study coefficients are ~0 (parallel-trends check passes); the post-treatment effect GROWS with")
print("time since adoption -- exactly the dynamic heterogeneity that biased TWFE. CS handles it by using only clean controls.")

true ATT                       = 3.580
naive twoway FE                = 2.677   (biased)
Callaway-SantAnna (from scratch)= 3.553   <- recovers the truth

Pre-treatment event-study coefficients are ~0 (parallel-trends check passes); the post-treatment effect GROWS with
time since adoption -- exactly the dynamic heterogeneity that biased TWFE. CS handles it by using only clean controls.

Causal Inference VI — Difference-in-Differences¶

The 2×2, parallel trends, and why staggered adoption broke (and fixed) applied econometrics¶

1. The canonical 2×2 — Card & Krueger (1994)¶

2. DiD as a regression, and the parallel-trends assumption¶

3. Staggered adoption — why naive twoway FE is biased¶

4. The fix — Callaway & Sant'Anna group-time ATTs, and the event study¶

5. Summary¶