1"""Okun's law, US 1960Q1-2019Q4 (task F6-conf-03). Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads data/ only; writes out/results.json, out/table1.csv, out/figure1.png, out/exclusions.csv."""4import json, os5import numpy as np, pandas as pd6import statsmodels.api as sm7from scipy import stats8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt1112OUT = "out"13os.makedirs(OUT, exist_ok=True)14gdp = pd.read_csv("data/GDPC1.csv", parse_dates=["observation_date"])15un = pd.read_csv("data/UNRATE.csv", parse_dates=["observation_date"])16n_input_rows = {"GDPC1.csv": int(len(gdp)), "UNRATE.csv": int(len(un))}17data_issues = []18miss = un[un["UNRATE"].isna()]19if len(miss):20 data_issues.append({"issue": f"UNRATE missing for {', '.join(miss.observation_date.dt.strftime('%Y-%m-%d'))}",21 "handling": "Outside the 1959Q4-2019Q4 window; removed by the 2019Q4 cutoff, no imputation needed."})22data_issues.append({"issue": "Both files extend beyond 2019Q4 (to 2026), outside the question.",23 "handling": "All rows after 2019-12-31 dropped before any computation."})2425# --- quarterly unemployment: simple average of the 3 months; require complete quarters in window26un["q"] = un.observation_date.dt.to_period("Q")27uq = un.groupby("q")["UNRATE"].agg(["mean", "count"])28gdp["q"] = gdp.observation_date.dt.to_period("Q")29df = gdp.set_index("q")[["GDPC1"]].join(uq, how="left")3031LAG0, START, END, SPLIT = pd.Period("1959Q4"), pd.Period("1960Q1"), pd.Period("2019Q4"), pd.Period("1984Q1")32excl = []33def step(d, mask, name, rationale):34 nb = len(d); d2 = d[mask]35 excl.append(dict(step=len(excl) + 1, rule=name, n_before=nb, n_after=len(d2), n_removed=nb - len(d2), rationale=rationale))36 return d237df = step(df, df.index >= LAG0, "drop quarters before 1959Q4", "Not needed; 1959Q4 kept only as lag level for 1960Q1 differences")38df = step(df, df.index <= END, "drop quarters after 2019Q4", "Information rule: post-2019Q4 quarters outside the question")39assert (df["count"] == 3).all() and df["mean"].notna().all() and df["GDPC1"].notna().all(), "incomplete data in window"40df = df.assign(u=df["mean"], du=df["mean"].diff(), g=100 * np.log(df["GDPC1"]).diff())41df = step(df, df.index >= START, "drop 1959Q4 (lag-only row)", "First differences undefined for the 1959Q4 row in-sample")42pd.DataFrame(excl).to_csv(f"{OUT}/exclusions.csv", index=False)43assert len(df) == 240 and df[["du", "g"]].notna().all().all()4445def nw_lags(n): # Newey-West (1994) rule of thumb46 return int(np.floor(4 * (n / 100) ** (2 / 9)))4748def fit(d):49 X = sm.add_constant(d["g"].values)50 L = nw_lags(len(d))51 r = sm.OLS(d["du"].values, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": True})52 return r, L5354full, L = fit(df)55s1, s2 = df[df.index < SPLIT], df[df.index >= SPLIT]56r1, L1 = fit(s1); r2, L2 = fit(s2)5758n = len(df); k = 259ssr_p, ssr_1, ssr_2 = full.ssr, r1.ssr, r2.ssr60F = ((ssr_p - (ssr_1 + ssr_2)) / k) / ((ssr_1 + ssr_2) / (n - 2 * k))61p = float(stats.f.sf(F, k, n - 2 * k))6263# supplementary robustness: HAC Wald test of break in (alpha, beta) and in beta alone64D = (df.index >= SPLIT).astype(float)65Xb = np.column_stack([np.ones(n), df.g.values, D, D * df.g.values])66rb = sm.OLS(df.du.values, Xb).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": True})67w_joint = rb.wald_test(np.array([[0, 0, 1, 0], [0, 0, 0, 1]]), scalar=True, use_f=False)68w_slope = rb.wald_test(np.array([[0, 0, 0, 1]]), scalar=True, use_f=False)6970beta, se = float(full.params[1]), float(full.bse[1])71z = stats.norm.ppf(0.975)72lo, hi = beta - z * se, beta + z * se73concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")74baseline = float(np.sqrt(np.mean((df.du - df.du.mean()) ** 2)))75model = float(np.sqrt(np.mean(full.resid ** 2)))7677tab = pd.DataFrame([78 dict(period="full 1960Q1-2019Q4", n=n, intercept=full.params[0], slope=full.params[1], slope_se=full.bse[1]),79 dict(period="sub1 1960Q1-1983Q4", n=len(s1), intercept=r1.params[0], slope=r1.params[1], slope_se=r1.bse[1]),80 dict(period="sub2 1984Q1-2019Q4", n=len(s2), intercept=r2.params[0], slope=r2.params[1], slope_se=r2.bse[1]),81])82tab.to_csv(f"{OUT}/table1.csv", index=False, float_format="%.6f")8384fig, ax = plt.subplots(figsize=(7, 5))85ax.scatter(s1.g, s1.du, s=14, alpha=.6, label="1960Q1-1983Q4", color="tab:blue")86ax.scatter(s2.g, s2.du, s=14, alpha=.6, label="1984Q1-2019Q4", color="tab:orange", marker="s")87for d, r, c in [(s1, r1, "tab:blue"), (s2, r2, "tab:orange")]:88 xs = np.linspace(d.g.min(), d.g.max(), 50)89 ax.plot(xs, r.params[0] + r.params[1] * xs, color=c, lw=2, label=f"fit: du = {r.params[0]:.3f} + {r.params[1]:.3f} g")90ax.axhline(0, color="grey", lw=.5); ax.axvline(0, color="grey", lw=.5)91ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)"); ax.set_ylabel("Change in unemployment rate Δu_t (pp)")92ax.set_title(f"Okun's law, US 1960Q1-2019Q4 (Chow F={F:.2f}, p={p:.3f})"); ax.legend(fontsize=8)93fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)9495res = dict(96 task_id="F6-conf-03", feasibility="complete", feasibility_reason=None, data_issues=data_issues,97 n_input_rows=n_input_rows, n_analysed=int(full.nobs), random_seed=None,98 sample_start=str(df.index[0].start_time.date()), sample_end=str(df.index[-1].start_time.date()),99 n_sub1=int(len(s1)), n_sub2=int(len(s2)),100 baseline_metric=baseline, model_metric=model,101 effect_estimate=beta, intercept=float(full.params[0]), effect_se=se, hac_lags=int(L),102 effect_ci_low=float(lo), effect_ci_high=float(hi),103 beta_sub1=float(r1.params[1]), beta_sub2=float(r2.params[1]),104 chow_f_stat=float(F), chow_p_value=p,105 stability_conclusion="stable" if p >= 0.05 else "unstable", conclusion=concl,106 # supplementary (not required)107 ols_se_beta_iid=float(sm.OLS(df.du.values, sm.add_constant(df.g.values)).fit().bse[1]),108 hac_lags_sub1=int(L1), hac_lags_sub2=int(L2), r_squared=float(full.rsquared),109 ssr_pooled=float(ssr_p), ssr_sub1=float(ssr_1), ssr_sub2=float(ssr_2),110 intercept_sub1=float(r1.params[0]), intercept_sub2=float(r2.params[0]),111 hac_wald_break_joint_chi2=float(w_joint.statistic), hac_wald_break_joint_p=float(w_joint.pvalue),112 hac_wald_break_slope_chi2=float(w_slope.statistic), hac_wald_break_slope_p=float(w_slope.pvalue),113 sigma_resid_sub1=float(np.sqrt(ssr_1 / len(s1))), sigma_resid_sub2=float(np.sqrt(ssr_2 / len(s2))),114 se_beta_sub1=float(r1.bse[1]), se_beta_sub2=float(r2.bse[1]),115)116with open(f"{OUT}/results.json", "w") as f:117 json.dump(res, f, indent=2)118print(json.dumps({k: res[k] for k in ["effect_estimate", "effect_se", "effect_ci_low", "effect_ci_high", "beta_sub1", "beta_sub2", "chow_f_stat", "chow_p_value", "baseline_metric", "model_metric"]}, indent=1))119