import numpy as np, pandas as pd, matplotlib.pyplot as plt, statsmodels.formula.api as smf, warnings
warnings.filterwarnings("ignore")
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
d=pd.read_csv("emplUK.csv").sort_values(["firm","year"]).reset_index(drop=True)
for c in ["emp","wage","capital","output"]: d["l"+c]=np.log(d[c])
d["lemp_l1"]=d.groupby("firm")["lemp"].shift(1)
print(f"Arellano-Bond employment panel: {d.firm.nunique()} UK firms, {d.year.min()}-{d.year.max()}, {len(d)} firm-years (unbalanced, {d.groupby('firm').size().min()}-{d.groupby('firm').size().max()} yrs each)")
dd=d.dropna(subset=["lemp_l1"]).copy()
ols=smf.ols("lemp~lemp_l1+lwage+lcapital+loutput", dd).fit()
for c in ["lemp","lemp_l1","lwage","lcapital","loutput"]: dd[c+"_w"]=dd[c]-dd.groupby("firm")[c].transform("mean")
fe=smf.ols("lemp_w~lemp_l1_w+lwage_w+lcapital_w+loutput_w-1", dd).fit()
rho_ols=ols.params['lemp_l1']; rho_fe=fe.params['lemp_l1_w']
print(f"\n  persistence rho (coef on lagged log-employment):")
print(f"    pooled OLS = {rho_ols:.3f}   biased UP  (firm effect left in error correlates with the lag)")
print(f"    FE within  = {rho_fe:.3f}   biased DOWN (Nickell bias, O(1/T); here T is only ~7-9)")
print(f"    => Bond bracket: a consistent estimate should lie in [{rho_fe:.2f}, {rho_ols:.2f}]")
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].bar(["pooled OLS\n(up)","FE within\n(down)"],[rho_ols,rho_fe],color=[RED,ORANGE])
ax[0].axhspan(rho_fe,rho_ols,color=GREEN,alpha=.12); ax[0].set_ylabel("persistence rho"); ax[0].set_title("The two workhorses bracket the truth")
for i,v in enumerate([rho_ols,rho_fe]): ax[0].text(i,v+0.02,f"{v:.2f}",ha="center")
ex=d[d.firm.isin(d.firm.unique()[:6])]
for f,g in ex.groupby("firm"): ax[1].plot(g.year,g.lemp,marker="o",ms=3,alpha=.7)
ax[1].set_xlabel("year"); ax[1].set_ylabel("log employment"); ax[1].set_title("Employment is persistent within firms (6 example firms)")
plt.tight_layout(); plt.show()
print("Persistence is exactly what makes the lag matter -- and what will later weaken the difference-GMM instruments.")

Arellano-Bond employment panel: 140 UK firms, 1976-1984, 1031 firm-years (unbalanced, 7-9 yrs each)

  persistence rho (coef on lagged log-employment):
    pooled OLS = 0.932   biased UP  (firm effect left in error correlates with the lag)
    FE within  = 0.514   biased DOWN (Nickell bias, O(1/T); here T is only ~7-9)
    => Bond bracket: a consistent estimate should lie in [0.51, 0.93]

Persistence is exactly what makes the lag matter -- and what will later weaken the difference-GMM instruments.

firms=[g for _,g in d.groupby("firm")]
def fa(g):
    g=g.sort_values("year"); return g["lemp"].values,g["lwage"].values,g["lcapital"].values,g["loutput"].values
maxT=max(len(fa(g)[0]) for g in firms)

def anderson_hsiao(firms):                     # 2SLS: dY_t ~ [dY_{t-1},dX], instrument dY_{t-1} by level y_{t-2}
    Z,X,Y=[],[],[]
    for g in firms:
        n,w,k,ys=fa(g); T=len(n)
        for t in range(2,T):
            dw,dk,dys=w[t]-w[t-1],k[t]-k[t-1],ys[t]-ys[t-1]
            Z.append([n[t-2],dw,dk,dys]); X.append([n[t-1]-n[t-2],dw,dk,dys]); Y.append(n[t]-n[t-1])
    Z,X,Y=np.array(Z),np.array(X),np.array(Y); Pz=Z@np.linalg.pinv(Z.T@Z)@Z.T
    return (np.linalg.pinv(X.T@Pz@X)@(X.T@Pz@Y))[0]

def arellano_bond(firms,maxT):                 # one-step difference GMM, all lagged levels as instruments
    per=list(range(2,maxT)); m=len(per); lag_cols=sum(t-1 for t in per); ncol=lag_cols+3
    H=2*np.eye(m)-np.eye(m,k=1)-np.eye(m,k=-1); cs={}; c=0
    for t in per: cs[t]=c; c+=t-1
    SZZ=np.zeros((ncol,ncol)); SZX=np.zeros((ncol,4)); SZY=np.zeros(ncol)
    for g in firms:
        n,w,k,ys=fa(g); T=len(n); Zi=np.zeros((m,ncol)); Xi=np.zeros((m,4)); Yi=np.zeros(m)
        for r,t in enumerate(per):
            if t>=T: continue
            lags=n[0:t-1]; Zi[r,cs[t]:cs[t]+len(lags)]=lags
            dw,dk,dys=w[t]-w[t-1],k[t]-k[t-1],ys[t]-ys[t-1]; Zi[r,lag_cols:]=[dw,dk,dys]
            Xi[r]=[n[t-1]-n[t-2],dw,dk,dys]; Yi[r]=n[t]-n[t-1]
        SZZ+=Zi.T@H@Zi; SZX+=Zi.T@Xi; SZY+=Zi.T@Yi
    W=np.linalg.pinv(SZZ); return (np.linalg.solve(SZX.T@W@SZX, SZX.T@W@SZY))[0]

rho_ah=anderson_hsiao(firms); rho_ab=arellano_bond(firms,maxT)
print(f"persistence rho, four ways:")
print(f"  pooled OLS             = {rho_ols:.3f}  (up)")
print(f"  FE within              = {rho_fe:.3f}  (down, Nickell)")
print(f"  Anderson-Hsiao IV      = {rho_ah:.3f}  (consistent, single lag instrument -- inside the bracket)")
print(f"  Arellano-Bond diff-GMM = {rho_ab:.3f}  (all lags; falls BELOW FE -- a warning sign)")
fig,ax=plt.subplots(figsize=(8.7,4.2))
nm=["pooled\nOLS","FE\n(Nickell)","Anderson-\nHsiao","Arellano-Bond\ndiff-GMM"]; vals=[rho_ols,rho_fe,rho_ah,rho_ab]
ax.bar(nm,vals,color=[RED,ORANGE,BLUE,PURP]); ax.axhspan(rho_fe,rho_ols,color=GREEN,alpha=.12,label="Bond bracket")
for i,v in enumerate(vals): ax.text(i,v+0.015,f"{v:.2f}",ha="center")
ax.set_ylabel("persistence rho"); ax.set_title("Anderson-Hsiao lands in the bracket; diff-GMM drops below FE"); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()
print("Difference GMM sitting below the downward-biased FE estimate is the classic symptom of WEAK instruments on a persistent")
print("series: lagged levels barely move the differenced lag. The fix is system GMM.")

persistence rho, four ways:
  pooled OLS             = 0.932  (up)
  FE within              = 0.514  (down, Nickell)
  Anderson-Hsiao IV      = 0.584  (consistent, single lag instrument -- inside the bracket)
  Arellano-Bond diff-GMM = 0.336  (all lags; falls BELOW FE -- a warning sign)

Difference GMM sitting below the downward-biased FE estimate is the classic symptom of WEAK instruments on a persistent
series: lagged levels barely move the differenced lag. The fix is system GMM.

def system_gmm(firms,maxT):
    dper=list(range(2,maxT)); lper=list(range(2,maxT))
    dcols=sum(t-1 for t in dper)+3; lcols=len(lper)+3+1; ncol=dcols+lcols
    globals()["N_INSTR"]=ncol                              # exposed so the diagnostics below can report it
    md_=len(dper); ml=len(lper); m=md_+ml
    H=np.zeros((m,m)); H[:md_,:md_]=2*np.eye(md_)-np.eye(md_,k=1)-np.eye(md_,k=-1); H[md_:,md_:]=np.eye(ml)
    cs={}; c=0
    for t in dper: cs[t]=c; c+=t-1
    P=5; SZZ=np.zeros((ncol,ncol)); SZX=np.zeros((ncol,P)); SZY=np.zeros(ncol)
    for g in firms:
        n,w,k,ys=fa(g); T=len(n); Zi=np.zeros((m,ncol)); Xi=np.zeros((m,P)); Yi=np.zeros(m)
        for r,t in enumerate(dper):                                  # differenced equations
            if t>=T: continue
            lags=n[0:t-1]; Zi[r,cs[t]:cs[t]+len(lags)]=lags
            dw,dk,dys=w[t]-w[t-1],k[t]-k[t-1],ys[t]-ys[t-1]; Zi[r,dcols-3:dcols]=[dw,dk,dys]
            Xi[r]=[n[t-1]-n[t-2],dw,dk,dys,0.0]; Yi[r]=n[t]-n[t-1]
        for r,t in enumerate(lper):                                  # levels equations (BB instrument = dy_{t-1})
            row=md_+r
            if t>=T: continue
            Zi[row,dcols+r]=n[t-1]-n[t-2]; Zi[row,dcols+ml:dcols+ml+3]=[w[t],k[t],ys[t]]; Zi[row,ncol-1]=1.0
            Xi[row]=[n[t-1],w[t],k[t],ys[t],1.0]; Yi[row]=n[t]
        SZZ+=Zi.T@H@Zi; SZX+=Zi.T@Xi; SZY+=Zi.T@Yi
    W=np.linalg.pinv(SZZ); return np.linalg.solve(SZX.T@W@SZX, SZX.T@W@SZY)
theta=system_gmm(firms,maxT); rho_sys=theta[0]
# Arellano-Bond AR(2) test: second-order serial correlation in the differenced residuals
res={}
for g in firms:
    n,w,k,ys=fa(g); T=len(n)
    for t in range(2,T): res[(id(g),t)]=(n[t]-n[t-1]) - rho_ab*(n[t-1]-n[t-2])
pairs=[]
for g in firms:
    n,w,k,ys=fa(g); T=len(n)
    for t in range(4,T):
        a=res.get((id(g),t)); b=res.get((id(g),t-2))
        if a is not None and b is not None: pairs.append((a,b))
pa=np.array(pairs); r2=np.corrcoef(pa[:,0],pa[:,1])[0,1]; z2=r2*np.sqrt(len(pa)); from scipy.stats import norm
print(f"persistence rho:")
print(f"  Arellano-Bond diff-GMM     = {rho_ab:.3f}  (downward biased by persistence)")
print(f"  system GMM (Blundell-Bond) = {rho_sys:.3f}  <-- restored INTO the Bond bracket [{rho_fe:.2f},{rho_ols:.2f}]")
print(f"  AR(2) test on differenced residuals: z={z2:.2f}, p={2*(1-norm.cdf(abs(z2))):.3f}  (no 2nd-order corr -> lagged-level instruments valid)")
fig,ax=plt.subplots(figsize=(8.7,4.2))
nm=["pooled\nOLS","FE","Anderson-\nHsiao","diff-GMM","system\nGMM"]; vals=[rho_ols,rho_fe,rho_ah,rho_ab,rho_sys]
ax.bar(nm,vals,color=[RED,ORANGE,BLUE,PURP,GREEN]); ax.axhspan(rho_fe,rho_ols,color=GREEN,alpha=.12,label="Bond bracket")
for i,v in enumerate(vals): ax.text(i,v+0.015,f"{v:.2f}",ha="center")
ax.set_ylabel("persistence rho"); ax.set_title("System GMM restores a sensible, in-bracket persistence estimate"); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"  instruments used by the system-GMM moment set: {N_INSTR}   (against {d.firm.nunique()} firms)")
print()
print("Now the diagnostic this notebook told you to always report, and which it is easy to skip because")
print("it does not cooperate. The R companion runs the same one-step system GMM through plm::pgmm and")
print("gets rho = 0.658 against 0.64 here, an AR(2) p of 0.670 -- and a SARGAN OVER-IDENTIFICATION TEST")
print("THAT REJECTS AT p = 0.003.")
print()
print("Two readings, and the honest position is that neither rescues the other. Taken at face value the")
print("rejection says the moment conditions are not all valid, which would undercut the estimate this")
print("section just restored. But Sargan is known to over-reject badly when instruments proliferate, and")
print(f"the R specification asks for every available lag (2:99) while this one uses {N_INSTR} moment columns --")
print("so the test is simultaneously failing AND unreliable, which is the least useful combination.")
print()
print("What that leaves is a weaker claim than 'system GMM fixes it': the estimate is back inside the Bond")
print("bracket and the AR(2) test is clean, so the lagged-level instruments are not obviously invalid --")
print("but the over-identifying restrictions are not confirmed, and with this many instruments they")
print("cannot be. The standard remedy is to collapse the instrument set or cap the lag depth and check")
print("that rho survives; that is the honest next step rather than a footnote.")
print()
print("System GMM's lagged-difference instruments stay informative under persistence, pulling the estimate back into [FE, OLS].")

persistence rho:
  Arellano-Bond diff-GMM     = 0.336  (downward biased by persistence)
  system GMM (Blundell-Bond) = 0.635  <-- restored INTO the Bond bracket [0.51,0.93]
  AR(2) test on differenced residuals: z=-0.49, p=0.623  (no 2nd-order corr -> lagged-level instruments valid)

  instruments used by the system-GMM moment set: 42   (against 140 firms)

Now the diagnostic this notebook told you to always report, and which it is easy to skip because
it does not cooperate. The R companion runs the same one-step system GMM through plm::pgmm and
gets rho = 0.658 against 0.64 here, an AR(2) p of 0.670 -- and a SARGAN OVER-IDENTIFICATION TEST
THAT REJECTS AT p = 0.003.

Two readings, and the honest position is that neither rescues the other. Taken at face value the
rejection says the moment conditions are not all valid, which would undercut the estimate this
section just restored. But Sargan is known to over-reject badly when instruments proliferate, and
the R specification asks for every available lag (2:99) while this one uses 42 moment columns --
so the test is simultaneously failing AND unreliable, which is the least useful combination.

What that leaves is a weaker claim than 'system GMM fixes it': the estimate is back inside the Bond
bracket and the AR(2) test is clean, so the lagged-level instruments are not obviously invalid --
but the over-identifying restrictions are not confirmed, and with this many instruments they
cannot be. The standard remedy is to collapse the instrument set or cap the lag depth and check
that rho survives; that is the honest next step rather than a footnote.

System GMM's lagged-difference instruments stay informative under persistence, pulling the estimate back into [FE, OLS].

Causal Inference V(b) — Dynamic Panels¶

Lagged dependent variables, Nickell bias, and GMM — on Arellano & Bond's own employment data¶

1. The lagged dependent variable and the OLS–FE bracket¶

2. First-difference the fixed effect away, then instrument with lagged levels¶

3. Persistence, weak instruments, and system GMM¶

4. Summary¶