1"""F6-conf-03 Okun's law replication. 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"; os.makedirs(OUT, exist_ok=True)13SEED = None # no randomness used14START, END = pd.Timestamp("1960-01-01"), pd.Timestamp("2019-10-01")15SPLIT = pd.Timestamp("1984-01-01")1617gdp = pd.read_csv("data/GDPC1.csv")18un = pd.read_csv("data/UNRATE.csv")19n_input = {"GDPC1.csv": int(len(gdp)), "UNRATE.csv": int(len(un))}20gdp["date"] = pd.to_datetime(gdp["observation_date"]); un["date"] = pd.to_datetime(un["observation_date"])21gdp["GDPC1"] = pd.to_numeric(gdp["GDPC1"], errors="coerce")22un["UNRATE"] = pd.to_numeric(un["UNRATE"], errors="coerce")2324data_issues = []25miss_un = un.loc[un["UNRATE"].isna(), "observation_date"].tolist()26if miss_un:27 data_issues.append({"issue": f"UNRATE missing (empty value) for month(s) {miss_un}",28 "handling": "Outside 1959Q4-2019Q4 window; affected quarter set to missing and in any case excluded by the window restriction, no imputation."})29miss_g = gdp.loc[gdp["GDPC1"].isna(), "observation_date"].tolist()30if miss_g:31 data_issues.append({"issue": f"GDPC1 missing for {miss_g}", "handling": "dropped if outside window; fatal if inside"})32data_issues.append({"issue": "Files extend beyond 2019Q4 (GDPC1 to 2026Q3, UNRATE to 2026-09), including COVID-19 quarters",33 "handling": "Excluded by the information rule: no quarter after 2019Q4 enters any estimate."})34data_issues.append({"issue": "GDPC1 is the current (2026) vintage, revised data, not real-time",35 "handling": "Used as provided; noted as limitation (estimates reflect revised GDP)."})3637# Quarterly unemployment: mean of 3 months, require complete quarter38un["q"] = un["date"].dt.to_period("Q")39uq = un.groupby("q")["UNRATE"].agg(["mean", "count"])40uq.loc[uq["count"] < 3, "mean"] = np.nan # incomplete quarter -> missing41g = gdp.set_index(gdp["date"].dt.to_period("Q"))["GDPC1"]42df = pd.DataFrame({"u": uq["mean"], "y": g}).sort_index()43df["du"] = df["u"].diff()44df["g"] = 100 * np.log(df["y"]).diff()45df["date"] = df.index.to_timestamp()4647excl = []48n0 = len(df)49step1 = df[(df["date"] >= START) & (df["date"] <= END)]50excl.append(["1", "keep quarters 1960Q1-2019Q4 (union of GDPC1/UNRATE quarterly index)", n0, len(step1), n0 - len(step1),51 "information rule: regression sample; pre-1960 used only as lag for 1960Q1 difference; post-2019 excluded"])52step2 = step1.dropna(subset=["du", "g"])53excl.append(["2", "drop quarters with missing du or g", len(step1), len(step2), len(step1) - len(step2),54 "no missing values expected inside the window"])55pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(f"{OUT}/exclusions.csv", index=False)56d = step2.copy()57assert len(d) == 240, len(d)58assert d["date"].max() <= END5960def nw_lags(n): # Newey-West (1994) plug-in rule of thumb61 return int(np.floor(4 * (n / 100) ** (2 / 9)))6263def fit(sub):64 X = sm.add_constant(sub["g"].values); y = sub["du"].values65 L = nw_lags(len(sub))66 ols = sm.OLS(y, X).fit()67 hac = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": True})68 return ols, hac, L6970full_ols, full_hac, L = fit(d)71s1 = d[d["date"] < SPLIT]; s2 = d[d["date"] >= SPLIT]72o1, h1, L1 = fit(s1); o2, h2, L2 = fit(s2)7374# baseline / model metrics75baseline_rmse = float(np.sqrt(np.mean((d["du"] - d["du"].mean()) ** 2)))76model_rmse = float(np.sqrt(np.mean(full_ols.resid ** 2)))77beta = float(full_ols.params[1]); alpha = float(full_ols.params[0])78se = float(full_hac.bse[1])79z = stats.norm.ppf(0.975)80lo, hi = beta - z * se, beta + z * se8182# Chow test83n, k = len(d), 284ssr_p, ssr_1, ssr_2 = full_ols.ssr, o1.ssr, o2.ssr85F = ((ssr_p - (ssr_1 + ssr_2)) / k) / ((ssr_1 + ssr_2) / (n - 2 * k))86p = float(stats.f.sf(F, k, n - 2 * k))878889# Supplementary robustness for the stability question (not the required Chow test):90# HAC Wald tests in a pooled regression with subperiod-2 dummy and interaction.91D = (d["date"] >= SPLIT).astype(float).values92Xi = np.column_stack([np.ones(n), d["g"].values, D, D * d["g"].values])93hi_fit = sm.OLS(d["du"].values, Xi).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": True})94w_joint = hi_fit.wald_test(np.array([[0, 0, 1, 0], [0, 0, 0, 1]]), use_f=True, scalar=True)95w_slope = hi_fit.wald_test(np.array([[0, 0, 0, 1]]), use_f=True, scalar=True)96var_ratio = (ssr_1 / (len(s1) - 2)) / (ssr_2 / (len(s2) - 2))9798conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")99res = {100 "task_id": "F6-conf-03",101 "sample_start": d["date"].min().strftime("%Y-%m-%d"),102 "sample_end": d["date"].max().strftime("%Y-%m-%d"),103 "n_sub1": int(len(s1)), "n_sub2": int(len(s2)),104 "baseline_metric": baseline_rmse, "model_metric": model_rmse,105 "effect_estimate": beta, "intercept": alpha, "effect_se": se, "hac_lags": L,106 "effect_ci_low": float(lo), "effect_ci_high": float(hi),107 "beta_sub1": float(o1.params[1]), "beta_sub2": float(o2.params[1]),108 "chow_f_stat": float(F), "chow_p_value": p,109 "stability_conclusion": "stable" if p >= 0.05 else "unstable",110 "conclusion": conclusion,111 "feasibility": "complete", "feasibility_reason": None,112 "data_issues": data_issues, "n_input_rows": n_input,113 "n_analysed": int(n), "random_seed": SEED,114 # supplementary (not required)115 "ols_classical_se": float(full_ols.bse[1]),116 "r_squared": float(full_ols.rsquared),117 "se_sub1_hac": float(h1.bse[1]), "se_sub2_hac": float(h2.bse[1]),118 "hac_lags_sub1": L1, "hac_lags_sub2": L2,119 "intercept_sub1": float(o1.params[0]), "intercept_sub2": float(o2.params[0]),120 "ssr_pooled": float(ssr_p), "ssr_sub1": float(ssr_1), "ssr_sub2": float(ssr_2),121 "supp_hac_wald_joint_F": float(w_joint.statistic), "supp_hac_wald_joint_p": float(w_joint.pvalue),122 "supp_hac_wald_slope_F": float(w_slope.statistic), "supp_hac_wald_slope_p": float(w_slope.pvalue),123 "supp_resid_var_ratio_sub1_over_sub2": float(var_ratio),124 "supp_var_ratio_F_p_two_sided": float(2 * min(stats.f.sf(var_ratio, len(s1) - 2, len(s2) - 2), stats.f.cdf(var_ratio, len(s1) - 2, len(s2) - 2))),125 "sensitivity_effect_se_by_lag": {str(l): float(sm.OLS(d["du"].values, sm.add_constant(d["g"].values)).fit(126 cov_type="HAC", cov_kwds={"maxlags": l, "kernel": "bartlett", "use_correction": True}).bse[1]) for l in (0, 2, 4, 8, 12)},127}128with open(f"{OUT}/results.json", "w") as f:129 json.dump(res, f, indent=2)130131tab = pd.DataFrame([132 ["full_1960Q1-2019Q4", n, alpha, beta, se],133 ["sub1_1960Q1-1983Q4", len(s1), o1.params[0], o1.params[1], h1.bse[1]],134 ["sub2_1984Q1-2019Q4", len(s2), o2.params[0], o2.params[1], h2.bse[1]],135], columns=["period", "n", "intercept", "slope", "slope_se"])136tab.to_csv(f"{OUT}/table1.csv", index=False)137138fig, ax = plt.subplots(figsize=(7, 5))139ax.scatter(s1["g"], s1["du"], s=14, alpha=.6, label="1960Q1–1983Q4")140ax.scatter(s2["g"], s2["du"], s=14, alpha=.6, marker="^", label="1984Q1–2019Q4")141xs = np.linspace(d["g"].min(), d["g"].max(), 50)142ax.plot(xs, o1.params[0] + o1.params[1] * xs, lw=2, label=f"fit 1: β={o1.params[1]:.3f}")143ax.plot(xs, o2.params[0] + o2.params[1] * xs, lw=2, ls="--", label=f"fit 2: β={o2.params[1]:.3f}")144ax.axhline(0, c="grey", lw=.5); ax.axvline(0, c="grey", lw=.5)145ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)")146ax.set_ylabel("Δ unemployment rate Δu_t (pp)")147ax.set_title(f"Okun's law, US 1960Q1–2019Q4 (Chow F={F:.2f}, p={p:.3f})")148ax.legend(fontsize=8); fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120)149150print(json.dumps({k: res[k] for k in ["effect_estimate", "effect_se", "hac_lags", "effect_ci_low", "effect_ci_high",151 "beta_sub1", "beta_sub2", "chow_f_stat", "chow_p_value", "stability_conclusion",152 "conclusion", "baseline_metric", "model_metric"]}, indent=1))153