# ===== hdt_ps2_analysis_20260903.py =====
#!/usr/bin/env python3
"""hdt_ps2_analysis.py — implements HDT_PS2_PREREGISTRATION_20260902.md (sha256 2ce6a578…3f1608) exactly,
with §12 amendment 4 (3 Sep 2026) applied: see the GB_CODE_CROSSWALK block below. Otherwise byte-identical to the
2 Sep 2026 run (hdt_ps2_analysis_20260902.py), which reproduced its frozen outputs exactly before this amendment was made.
Inputs: hdt_2025_clean.csv · PS2_open_data_202603.csv · GB_2024-25.ods · lpa_metrics.csv (PlanningLens run 3).
Outputs: hdt_ps2_lpa.csv · hdt_ps2_summary.md. No narrative."""
import csv, io, math, random, sys
import numpy as np, pandas as pd

SEED = 20260902; B = 2000
rng = np.random.default_rng(SEED)
QUARTERS = ["2022 Q2","2022 Q3","2022 Q4","2023 Q1","2023 Q2","2023 Q3","2023 Q4","2024 Q1","2024 Q2","2024 Q3","2024 Q4","2025 Q1"]
MAJ = ("Total decisions; major dwellings (all)", "Total granted; major dwellings (all)", "Total refused; major dwellings (all)")
MIN = ("Total decisions; minor dwellings (all)", "Total granted; minor dwellings (all)", "Total refused; minor dwellings (all)")

# ------------------------------------------------------------------ load
hdt = pd.read_csv("/home/claude/hdt/hdt_2025_clean.csv")
hdt = hdt[hdt["hdt_pct"].notna()].copy()
ps2 = pd.read_csv("/home/claude/ps2/PS2_open_data_202603.csv", skiprows=2, encoding="cp1252", low_memory=False,
                  usecols=["Region","LPANM","LPACD","Quarter",*MAJ,*MIN])
ps2 = ps2[ps2["Quarter"].isin(QUARTERS)].copy()
for c in (*MAJ, *MIN): ps2[c] = pd.to_numeric(ps2[c], errors="coerce")
g = ps2.groupby("LPACD")
agg = pd.DataFrame({
    "lpanm": g["LPANM"].first(), "region": g["Region"].first(),
    "quarters_present": g[MAJ[0]].apply(lambda s: s.notna().sum()),
    "maj_dec": g[MAJ[0]].sum(min_count=1), "maj_gr": g[MAJ[1]].sum(min_count=1), "maj_ref": g[MAJ[2]].sum(min_count=1),
    "min_dec": g[MIN[0]].sum(min_count=1), "min_gr": g[MIN[1]].sum(min_count=1), "min_ref": g[MIN[2]].sum(min_count=1),
}).reset_index().rename(columns={"LPACD":"ons_code"})
df = hdt.merge(agg, on="ons_code", how="left")
df["matched_ps2"] = df["lpanm"].notna()
df["london"] = df["ons_code"].str.startswith("E09") | (df["region"] == "London")
df["dev_corp_or_nondistrict"] = ~df["ons_code"].str[:3].isin(["E06","E07","E08","E09"])
df["maj_appr"] = 100 * df["maj_gr"] / df["maj_dec"]; df["maj_refr"] = 100 * df["maj_ref"] / df["maj_dec"]
df["min_appr"] = 100 * df["min_gr"] / df["min_dec"]
df["maj_discrep"] = df["maj_dec"] - (df["maj_gr"] + df["maj_ref"])
df["maj_discrep_flag"] = (df["maj_discrep"].abs() > 0.02 * df["maj_dec"])
# Green Belt
gb = pd.read_excel("/home/claude/ps2/GB_2024-25.ods", engine="odf", sheet_name="Area_by_LA", header=None, skiprows=3)
gb.columns = ["gb_name","ons_code","gb_ha","land_ha","gb_pct"][:gb.shape[1]]
gb = gb[gb["ons_code"].astype(str).str.match(r"^E0[6-9]\d{6}$", na=False)].copy()
for c in ("gb_ha","land_ha","gb_pct"): gb[c] = pd.to_numeric(gb[c], errors="coerce")
# --- §12 AMENDMENT 4 (3 Sep 2026): ONS-code crosswalk. MHCLG's HDT 2025 and PS2 files carry Sheffield as E08000039 and Barnsley as
# E08000038; MHCLG's Green Belt live table carries the same two authorities as E08000019 and E08000016. Without this mapping neither
# Green Belt row joins and both authorities are treated as having no Green Belt (which is what the 2 Sep 2026 frozen run did). The
# mismatch was found after the frozen results had been read; the correction is an identifier fix, not an analytical choice, and is
# logged in the pre-registration as amendment 4. No other line of this script differs from the 2 Sep run.
GB_CODE_CROSSWALK = {"E08000019": "E08000039", "E08000016": "E08000038"}
gb["ons_code"] = gb["ons_code"].replace(GB_CODE_CROSSWALK)
df = df.merge(gb[["ons_code","gb_ha","land_ha","gb_pct"]], on="ons_code", how="left")
df["gb_any"] = df["gb_ha"].fillna(0) > 0
df["gb_25"] = df["gb_pct"].fillna(0) >= 25
# PlanningLens layer (run 3)
pl = pd.read_csv("/home/claude/hdt_out/hdt_feasibility_20260902_1830/lpa_metrics.csv")
pl["zero_read"] = (pl["major_new_decided"].fillna(0) == 0) & (pl["minor_new_decided"].fillna(0) == 0) & (pl["procedural_excluded"].fillna(0) == 0)
df = df.merge(pl[["ons_code","major_new_decided","major_new_approval_pct","zero_read"]].rename(
    columns={"major_new_decided":"pl_major_dec","major_new_approval_pct":"pl_major_appr"}), on="ons_code", how="left")
df["pl_ratio"] = df["pl_major_dec"] / df["maj_dec"]
df["pl_layer_ok"] = (df["pl_ratio"] >= 0.40) & (df["zero_read"] == False)
# bands
def band_primary(h):
    return "<75 presumption" if h < 75 else "75-<85 buffer" if h < 85 else "85-<95 action plan" if h < 95 else ">=95 none"
def band_secondary(h):
    return "<75" if h < 75 else "75-<95" if h < 95 else ">=95"
df["band_primary"] = df["hdt_pct"].apply(band_primary); df["band_secondary"] = df["hdt_pct"].apply(band_secondary)
df["elig20"] = df["matched_ps2"] & (df["quarters_present"] >= 10) & (df["maj_dec"] >= 20)
df["elig40"] = df["elig20"] & (df["maj_dec"] >= 40)
df["elig_min60"] = df["matched_ps2"] & (df["quarters_present"] >= 10) & (df["min_dec"] >= 60)
df.to_csv("/home/claude/ps2/hdt_ps2_lpa.csv", index=False)

# ------------------------------------------------------------------ stats
def spearman(x, y):
    xr = pd.Series(x).rank().values; yr = pd.Series(y).rank().values
    if len(xr) < 5 or np.std(xr) == 0 or np.std(yr) == 0: return np.nan
    return float(np.corrcoef(xr, yr)[0, 1])
def spearman_ci(x, y):
    x = np.asarray(x, float); y = np.asarray(y, float); n = len(x)
    rho = spearman(x, y)
    idx = rng.integers(0, n, size=(B, n))
    boots = np.array([spearman(x[i], y[i]) for i in idx]); boots = boots[~np.isnan(boots)]
    lo, hi = np.percentile(boots, [2.5, 97.5])
    z = np.arctanh(np.clip(rho, -0.999, 0.999)); se = 1.06 / math.sqrt(n - 3) if n > 3 else np.nan
    fz = (math.tanh(z - 1.96 * se), math.tanh(z + 1.96 * se)) if n > 3 else (np.nan, np.nan)
    return rho, lo, hi, fz, n
def interpret(rho, lo, hi):
    if np.isnan(rho): return "n/a"
    within30 = lo > -0.30 and hi < 0.30; within15 = lo > -0.15 and hi < 0.15
    if within15 and abs(rho) < 0.15: return "NO MEANINGFUL RELATIONSHIP (strong null: CI within ±0.15)"
    if within30 and abs(rho) < 0.15: return "NO MEANINGFUL RELATIONSHIP (CI within ±0.30; point |rho|<0.15)"
    if abs(rho) < 0.15: return "point estimate in null region but CI extends beyond ±0.30 — cannot claim no relationship"
    if abs(rho) < 0.30: return f"WEAK {'positive' if rho>0 else 'negative'} association"
    return f"MODERATE-OR-STRONGER {'positive' if rho>0 else 'negative'} association"
def med_ci(v):
    v = np.asarray(v, float); v = v[~np.isnan(v)]
    if len(v) == 0: return (np.nan, np.nan, np.nan, np.nan, np.nan, 0)
    idx = rng.integers(0, len(v), size=(B, len(v))); boots = np.median(v[idx], axis=1)
    q1, q3 = np.percentile(v, [25, 75])
    return (float(np.median(v)), float(np.percentile(boots, 2.5)), float(np.percentile(boots, 97.5)), q1, q3, len(v))
def diff_med_ci(a, b):
    a = np.asarray(a, float); b = np.asarray(b, float)
    if len(a) < 3 or len(b) < 3: return (np.nan, np.nan, np.nan)
    ia = rng.integers(0, len(a), size=(B, len(a))); ib = rng.integers(0, len(b), size=(B, len(b)))
    d = np.median(a[ia], axis=1) - np.median(b[ib], axis=1)
    return (float(np.median(a) - np.median(b)), float(np.percentile(d, 2.5)), float(np.percentile(d, 97.5)))
def f(v, d=1): return "n/a" if v is None or (isinstance(v, float) and np.isnan(v)) else f"{v:.{d}f}"

S = io.StringIO()
def say(s=""): S.write(s + "\n")

def block(name, sub, appr="maj_appr", dec="maj_dec", gr="maj_gr", bands="band_primary"):
    say(f"\n### {name}  (n = {len(sub)} authorities)")
    if len(sub) < 5: say("  too few authorities"); return
    rho, lo, hi, fz, n = spearman_ci(sub["hdt_pct"], sub[appr])
    say(f"- Spearman rho (HDT % vs {appr}) = **{f(rho,3)}**  bootstrap 95% CI [{f(lo,3)}, {f(hi,3)}]  Fisher-z CI [{f(fz[0],3)}, {f(fz[1],3)}]  n={n}")
    say(f"- Interpretation (pre-registered rule): **{interpret(rho, lo, hi)}**")
    say(f"\n| HDT band | n LPAs | median approval % | bootstrap 95% CI | IQR | decisions-weighted approval % | pooled decisions |")
    say("|---|---|---|---|---|---|---|")
    order = ["<75 presumption","75-<85 buffer","85-<95 action plan",">=95 none"] if bands == "band_primary" else ["<75","75-<95",">=95"]
    for b in order:
        s = sub[sub[bands] == b]
        m, l, h, q1, q3, k = med_ci(s[appr])
        pooled = 100 * s[gr].sum() / s[dec].sum() if s[dec].sum() else np.nan
        say(f"| {b} | {k} | {f(m)} | [{f(l)}, {f(h)}] | {f(q1)}–{f(q3)} | {f(pooled)} | {int(s[dec].sum()) if k else 0} |")
    if bands == "band_primary":
        s100 = sub[sub["hdt_pct"] >= 100]; m, l, h, q1, q3, k = med_ci(s100[appr])
        say(f"| (>=100 subset) | {k} | {f(m)} | [{f(l)}, {f(h)}] | {f(q1)}–{f(q3)} | {f(100*s100[gr].sum()/s100[dec].sum() if s100[dec].sum() else np.nan)} | {int(s100[dec].sum()) if k else 0} |")
        lowb = sub[sub["hdt_pct"] < 75][appr]; highb = sub[sub["hdt_pct"] >= 95][appr]
        d, dl, dh = diff_med_ci(lowb, highb)
        say(f"- Band contrast, median approval <75 minus >=95: **{f(d)} pp**  bootstrap 95% CI [{f(dl)}, {f(dh)}]  (n {len(lowb)} vs {len(highb)})")

say("# HDT 2025 × PS2 major-dwellings approval — RESULTS (numbers only, no narrative)")
say(f"Pre-registration: HDT_PS2_PREREGISTRATION_20260902.md, sha256 2ce6a5780b73f32f403a704e9504d592e8f5320d086b80fc16606033c53f1608. Seed {SEED}, {B} bootstrap resamples.")
say("Amendment 4 applied (3 Sep 2026): Green Belt ONS-code crosswalk for Sheffield (E08000019→E08000039) and Barnsley (E08000016→E08000038); see pre-registration §12.")
say("\n## Matching and exclusions")
say(f"- HDT numeric authorities: {len(df)}; matched to PS2 on ONS code: {int(df['matched_ps2'].sum())}; unmatched: {sorted(df.loc[~df['matched_ps2'],'lpa'].tolist())}")
say(f"- Authorities with <10 of 12 quarters present: {int(((df['quarters_present']<10)&df['matched_ps2']).sum())} → {sorted(df.loc[(df['quarters_present']<10)&df['matched_ps2'],'lpa'].tolist())}")
say(f"- Eligible at floor 20 major decisions: **{int(df['elig20'].sum())}**; at floor 40: {int(df['elig40'].sum())}; minor floor 60: {int(df['elig_min60'].sum())}")
say(f"- Below the 20 floor (matched, quarters ok): {sorted(df.loc[df['matched_ps2']&(df['quarters_present']>=10)&(df['maj_dec']<20),'lpa'].tolist())}")
say(f"- Consistency check decisions ≠ granted+refused by >2%: {int(df['maj_discrep_flag'].sum())} authorities: {sorted(df.loc[df['maj_discrep_flag'],'lpa'].tolist())}")
e = df[df["elig20"]]
say(f"- Pooled over eligible authorities: major decisions {int(e['maj_dec'].sum())}, granted {int(e['maj_gr'].sum())}, refused {int(e['maj_ref'].sum())} → pooled approval {f(100*e['maj_gr'].sum()/e['maj_dec'].sum())}%")
say(f"- Major decisions per eligible authority per year: median {f(e['maj_dec'].median()/3)}, IQR {f(e['maj_dec'].quantile(.25)/3)}–{f(e['maj_dec'].quantile(.75)/3)}")
say(f"- Green Belt: {int(e['gb_any'].sum())} eligible authorities with any Green Belt, {int((~e['gb_any']).sum())} without; {int(e['gb_25'].sum())} with Green Belt ≥25% of land area")

say("\n## PRIMARY: major dwellings (all), floor ≥20 decisions, HDT consequence bands")
block("Cut 1 — England, all eligible", e)
exl = e[~e["london"]]; lon = e[e["london"]]
block("Cut 2 — England excluding London", exl)
block("Cut 3 — London only", lon)
say("\n## Sensitivities (major dwellings)")
block("Cut 4a — any Green Belt (England, all)", e[e["gb_any"]]); block("Cut 4a — no Green Belt (England, all)", e[~e["gb_any"]])
block("Cut 4a′ — any Green Belt, excluding London", exl[exl["gb_any"]]); block("Cut 4a′ — no Green Belt, excluding London", exl[~exl["gb_any"]])
block("Cut 4b — Green Belt ≥25% of land (England, all)", e[e["gb_25"]]); block("Cut 4b — Green Belt <25% of land (England, all)", e[~e["gb_25"]])
block("Cut 5 — excluding HDT >200% (England, all)", e[e["hdt_pct"] <= 200]); block("Cut 5′ — excluding HDT >200%, excluding London", exl[exl["hdt_pct"] <= 200])
block("Cut 6 — excluding development corporations / non-district (England, all)", e[~e["dev_corp_or_nondistrict"]])
e40 = df[df["elig40"]]
block("Cut 7 — floor ≥40 decisions, England all", e40); block("Cut 7 — floor ≥40, excluding London", e40[~e40["london"]]); block("Cut 7 — floor ≥40, London", e40[e40["london"]])
say("\n## Secondary bands (<75 / 75–<95 / ≥95), major dwellings")
block("England excluding London — secondary bands", exl, bands="band_secondary")
say("\n## SECONDARY OUTCOME: minor dwellings (all), floor ≥60 decisions")
em = df[df["elig_min60"]]
block("Minor — England, all", em, appr="min_appr", dec="min_dec", gr="min_gr")
block("Minor — England excluding London", em[~em["london"]], appr="min_appr", dec="min_dec", gr="min_gr")
block("Minor — London", em[em["london"]], appr="min_appr", dec="min_dec", gr="min_gr")

say("\n## Anomalies (rule-based, §9) — England excluding London, primary cut")
q75 = exl["maj_appr"].quantile(.75); q25 = exl["maj_appr"].quantile(.25)
say(f"Top-quartile threshold {f(q75)}%, bottom-quartile threshold {f(q25)}%")
def row(r):
    plt = f"PL {int(r['pl_major_dec'])} / PS2 {int(r['maj_dec'])} = {f(r['pl_ratio']*100,0)}% {'OK' if r['pl_layer_ok'] else 'NOT eligible'}" if not np.isnan(r.get("pl_ratio", np.nan)) else "PL: no data"
    return f"| {r['lpa']} | {f(r['hdt_pct'],0)}% | {r['consequence']} | {int(r['maj_dec'])} | {f(r['maj_appr'])}% | {int(r['req_total'])} | {'GB' if r['gb_any'] else '-'} | {plt} |"
hdr = "| Authority | HDT | Consequence | Major decisions | Approval | Homes required | Green Belt | PlanningLens layer |\n|---|---|---|---|---|---|---|---|"
say("\n### A. HDT <75% and major approval in top quartile"); say(hdr)
for _, r in exl[(exl["hdt_pct"] < 75) & (exl["maj_appr"] >= q75)].sort_values("maj_appr", ascending=False).iterrows(): say(row(r))
say("\n### B. HDT ≥100% and major approval in bottom quartile"); say(hdr)
for _, r in exl[(exl["hdt_pct"] >= 100) & (exl["maj_appr"] <= q25)].sort_values("maj_appr").iterrows(): say(row(r))
say("\n### C. 20 largest housing requirements (eligible, excluding London)"); say(hdr)
for _, r in exl.sort_values("req_total", ascending=False).head(20).iterrows(): say(row(r))
say("\n### London panel — all eligible boroughs"); say(hdr)
for _, r in lon.sort_values("hdt_pct").iterrows(): say(row(r))
say("\n### Full eligible list (England excluding London), sorted by HDT"); say(hdr)
for _, r in exl.sort_values("hdt_pct").iterrows(): say(row(r))
say(f"\n## PlanningLens layer eligibility\n- Authorities meeting the ≥40%-of-PS2 rule and readable: {int(df['pl_layer_ok'].sum())} of {int(df['pl_ratio'].notna().sum())} with PlanningLens data; median PL/PS2 ratio {f(df['pl_ratio'].median()*100,0)}%")
open("/home/claude/ps2/hdt_ps2_summary.md", "w").write(S.getvalue())
print(S.getvalue())


# ===== hdt_ps2_exploratory_gb_20260902.py =====
import pandas as pd, numpy as np, io
rng=np.random.default_rng(20260902); B=2000
df=pd.read_csv('hdt_ps2_lpa.csv'); ex=df[df['elig20'] & ~df['london']].copy(); ex['gb_pct']=ex['gb_pct'].fillna(0)
def rho(x,y):
    x=pd.Series(x).rank().values; y=pd.Series(y).rank().values
    return float(np.corrcoef(x,y)[0,1]) if len(x)>=5 and x.std()>0 and y.std()>0 else np.nan
def rho_ci(s):
    x=s['hdt_pct'].values; y=s['maj_appr'].values; n=len(x); r=rho(x,y)
    idx=rng.integers(0,n,size=(B,n)); b=np.array([rho(x[i],y[i]) for i in idx]); b=b[~np.isnan(b)]
    return r,np.percentile(b,2.5),np.percentile(b,97.5),n
def med_ci(v):
    v=np.asarray(v,float); v=v[~np.isnan(v)]
    if len(v)<3: return np.nan,np.nan,np.nan,len(v)
    b=np.median(v[rng.integers(0,len(v),size=(B,len(v)))],axis=1); return float(np.median(v)),np.percentile(b,2.5),np.percentile(b,97.5),len(v)
def gap_ci(a,b_):
    a=np.asarray(a,float); b_=np.asarray(b_,float)
    if len(a)<3 or len(b_)<3: return np.nan,np.nan,np.nan
    d=np.median(a[rng.integers(0,len(a),size=(B,len(a)))],axis=1)-np.median(b_[rng.integers(0,len(b_),size=(B,len(b_)))],axis=1)
    return float(np.median(a)-np.median(b_)),np.percentile(d,2.5),np.percentile(d,97.5)
S=io.StringIO(); say=lambda s="": S.write(s+"\n")
say("\n\n---\n# POST-PRE-REGISTRATION EXPLORATORY ANALYSIS — Green Belt interaction (2 Sep 2026, ~19:40)")
say("**Not pre-registered.** Motivated by the pre-registered Green Belt sensitivity (cut 4b). England excluding London throughout. Seed 20260902, 2,000 resamples.")
say("\n## E1. Green Belt ≥25% vs <25% of land area, at both decision floors")
say("| Floor | Group | n | Spearman rho [95% CI] | <75 median [CI] | >=95 median [CI] | gap <75 minus >=95 [CI] |\n|---|---|---|---|---|---|---|")
for floor in (20,40):
    e=ex[ex['maj_dec']>=floor]
    for label,m in (("GB >=25%",e['gb_pct']>=25),("GB <25%",e['gb_pct']<25)):
        s=e[m]; r,lo,hi,n=rho_ci(s); a=med_ci(s[s.hdt_pct<75]['maj_appr']); b_=med_ci(s[s.hdt_pct>=95]['maj_appr']); g=gap_ci(s[s.hdt_pct<75]['maj_appr'],s[s.hdt_pct>=95]['maj_appr'])
        say(f"| >={floor} | {label} | {n} | {r:+.3f} [{lo:+.3f}, {hi:+.3f}] | {a[0]:.1f} [{a[1]:.1f}, {a[2]:.1f}] (n={a[3]}) | {b_[0]:.1f} [{b_[1]:.1f}, {b_[2]:.1f}] (n={b_[3]}) | {g[0]:+.1f} [{g[1]:+.1f}, {g[2]:+.1f}] |")
say("\n## E2. Threshold sweep — is 25% doing special work? (floor >=20)")
say("| GB share threshold | heavy n | heavy rho [CI] | light n | light rho [CI] | difference in rho |\n|---|---|---|---|---|---|")
for t in (10,15,20,25,30,40,50):
    h=ex[ex['gb_pct']>=t]; l=ex[ex['gb_pct']<t]; rh=rho_ci(h); rl=rho_ci(l)
    say(f"| >={t}% | {rh[3]} | {rh[0]:+.3f} [{rh[1]:+.3f}, {rh[2]:+.3f}] | {rl[3]} | {rl[0]:+.3f} [{rl[1]:+.3f}, {rl[2]:+.3f}] | {rh[0]-rl[0]:+.3f} |")
say("\n## E3. Permutation test — does the HDT–approval association differ between GB >=25% and <25%? (floor >=20)")
for floor in (20,40):
    e=ex[ex['maj_dec']>=floor]; heavy=(e['gb_pct']>=25).values; x=e['hdt_pct'].values; y=e['maj_appr'].values
    obs=rho(x[heavy],y[heavy])-rho(x[~heavy],y[~heavy]); nperm=5000; cnt=0
    for _ in range(nperm):
        p=rng.permutation(heavy); d=rho(x[p],y[p])-rho(x[~p],y[~p])
        if abs(d)>=abs(obs): cnt+=1
    say(f"- floor >={floor}: observed rho(heavy) − rho(light) = **{obs:+.3f}**; two-sided permutation p = **{(cnt+1)/(nperm+1):.4f}** ({nperm} label shuffles, group sizes held fixed)")
say("\n## E4. Continuous moderator — rank regression: rank(approval) ~ rank(HDT) + rank(GB share) + rank(HDT)×rank(GB share)  (floor >=20, ex-London)")
e=ex.copy(); n=len(e)
rh=e['hdt_pct'].rank()/n; rg=e['gb_pct'].rank()/n; ry=e['maj_appr'].rank()/n
X=np.column_stack([np.ones(n),rh,rg,rh*rg]); beta,res,rk,sv=np.linalg.lstsq(X,ry,rcond=None)
resid=ry-X@beta; sigma2=(resid@resid)/(n-4); cov=sigma2*np.linalg.inv(X.T@X); se=np.sqrt(np.diag(cov)); t=beta/se
# bootstrap the interaction coefficient
bi=[]
for _ in range(B):
    i=rng.integers(0,n,size=n); Xi=X[i]; yi=ry.values[i]
    try: bi.append(np.linalg.lstsq(Xi,yi,rcond=None)[0][3])
    except Exception: pass
bi=np.array(bi)
say(f"- n = {n}. Coefficients (ranks scaled 0–1): HDT {beta[1]:+.3f} (t={t[1]:+.2f}), GB share {beta[2]:+.3f} (t={t[2]:+.2f}), **HDT × GB share {beta[3]:+.3f} (t={t[3]:+.2f}), bootstrap 95% CI [{np.percentile(bi,2.5):+.3f}, {np.percentile(bi,97.5):+.3f}]**")
say(f"- Reading: a positive interaction means the HDT–approval slope rises with Green Belt share. Implied slope of rank(approval) on rank(HDT): at GB share rank 0 → {beta[1]:+.3f}; at rank 0.5 → {beta[1]+0.5*beta[3]:+.3f}; at rank 1 → {beta[1]+beta[3]:+.3f}.")
say("\n## E5. Diagnostic — PS2 vs PlanningLens approval gap across scheme-layer-eligible authorities (diagnostic only, not a correction)")
d=df[df['pl_layer_ok']==True].copy(); d['gap']=d['maj_appr']-d['pl_major_appr']
say(f"- {len(d)} authorities pass the ≥40%-recall rule (all England). PS2 approval minus PlanningLens (RM/s73-free) approval: median {d['gap'].median():+.1f} pp, IQR {d['gap'].quantile(.25):+.1f} to {d['gap'].quantile(.75):+.1f}.")
say(f"- Spearman rho of the gap with HDT %: {rho(d['hdt_pct'],d['gap']):+.3f} (n={len(d)}); with Green Belt share: {rho(d['gb_pct'].fillna(0),d['gap']):+.3f}. Median gap, GB >=25%: {d[d['gb_pct'].fillna(0)>=25]['gap'].median():+.1f} pp (n={int((d['gb_pct'].fillna(0)>=25).sum())}); GB <25%: {d[d['gb_pct'].fillna(0)<25]['gap'].median():+.1f} pp (n={int((d['gb_pct'].fillna(0)<25).sum())}).")
say("- Both measures are imperfect (PlanningLens recall in these authorities is 40–68%), so this shows only that the RM/s73 inflation is not *obviously* patterned by HDT or Green Belt in the authorities where it can be observed. It does not establish uniformity.")
open('hdt_ps2_summary.md','a').write(S.getvalue()); print(S.getvalue())