1"""Okun's law, US 1960Q1-2019Q4 (task F6-conf-03). Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads only data/, writes out/results.json, out/table1.csv, out/figure1.png, out/exclusions.csv.4No randomness is used (random_seed = null)."""5import json, os6import numpy as np7import pandas as pd8import statsmodels.api as sm9from scipy import stats10import matplotlib11matplotlib.use("Agg")12import matplotlib.pyplot as plt1314DATA, OUT = "data", "out"15os.makedirs(OUT, exist_ok=True)16FIRST_Q, LAST_Q, SPLIT_Q = pd.Period("1960Q1"), pd.Period("2019Q4"), pd.Period("1984Q1")1718gdp_raw = pd.read_csv(f"{DATA}/GDPC1.csv", parse_dates=["observation_date"])19ur_raw = pd.read_csv(f"{DATA}/UNRATE.csv", parse_dates=["observation_date"])20n_input_rows = {"GDPC1.csv": int(len(gdp_raw)), "UNRATE.csv": int(len(ur_raw))}2122data_issues = []23ur = ur_raw.copy()24ur["UNRATE"] = pd.to_numeric(ur["UNRATE"], errors="coerce")25ur["q"] = ur["observation_date"].dt.to_period("Q")26miss = ur[ur["UNRATE"].isna()]27if len(miss):28 inside = bool(((miss["q"] >= FIRST_Q - 1) & (miss["q"] <= LAST_Q)).any())29 data_issues.append({30 "issue": f"UNRATE has {len(miss)} missing monthly value(s): "31 + ", ".join(d.strftime('%Y-%m') for d in miss["observation_date"]),32 "handling": ("inside the estimation window: affected quarter dropped (3 valid months required)" if inside else33 "outside the 1959Q4-2019Q4 window needed for estimation, so no effect; quarters with "34 "missing months are never averaged (3 valid months required); no imputation")})35data_issues.append({"issue": "Both files extend beyond 2019Q4 (to 2026); post-2019 quarters are outside the question",36 "handling": "truncated to 1959Q4-2019Q4 before any computation (1959Q4 used only as lag level)"})37data_issues.append({"issue": "GDPC1 is the 2026 FRED vintage (chained 2017 dollars), not real-time data",38 "handling": "used as provided; results are conditional on this vintage (noted as limitation)"})3940qagg = ur.groupby("q")["UNRATE"].agg(["mean", "count"])41qagg = qagg[qagg["count"] == 3] # complete quarters only42u = qagg["mean"].rename("u")43gdp = gdp_raw.assign(q=gdp_raw["observation_date"].dt.to_period("Q")).set_index("q")["GDPC1"].astype(float)4445df = pd.concat([u, gdp], axis=1, join="inner").sort_index()46df = df[(df.index >= FIRST_Q - 1) & (df.index <= LAST_Q)]47assert len(df) == 241 and (np.diff(df.index.asi8) == 1).all(), "gap in quarterly series"48df["du"] = df["u"].diff()49df["g"] = 100 * np.log(df["GDPC1"]).diff()5051n0 = len(df)52d = df.dropna()53pd.DataFrame([("1", "drop 1959Q4 (lag level only; difference undefined)", n0, len(d), n0 - len(d),54 "1959Q4 level used only to form 1960Q1 differences; post-2019Q4 rows were never in the window")],55 columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(56 f"{OUT}/exclusions.csv", index=False)57assert len(d) == 240 and d.index[0] == FIRST_Q and d.index[-1] == LAST_Q5859def ols(sub):60 return sm.OLS(sub["du"].values, sm.add_constant(sub["g"].values)).fit()6162def nw_lags(m): # Newey-West (1994) rule of thumb63 return int(np.floor(4 * (m / 100) ** (2 / 9)))6465def hac_se(fit, lags):66 return float(fit.get_robustcov_results(cov_type="HAC", maxlags=lags, kernel="bartlett",67 use_correction=False).bse[1])6869n = len(d)70hac_lags = nw_lags(n)71full = ols(d)72beta, alpha = float(full.params[1]), float(full.params[0])73se = hac_se(full, hac_lags)74z = stats.norm.ppf(0.975)75ci_lo, ci_hi = beta - z * se, beta + z * se7677baseline_rmse = float(np.sqrt(np.mean((d["du"] - d["du"].mean()) ** 2)))78model_rmse = float(np.sqrt(np.mean(full.resid ** 2)))7980s1, s2 = d[d.index < SPLIT_Q], d[d.index >= SPLIT_Q]81f1, f2 = ols(s1), ols(s2)82k = 283ssr_u = f1.ssr + f2.ssr84chow_f = float(((full.ssr - ssr_u) / k) / (ssr_u / (n - 2 * k)))85chow_p = float(stats.f.sf(chow_f, k, n - 2 * k))8687table = pd.DataFrame([88 {"period": "full_1960Q1-2019Q4", "n": n, "intercept": alpha, "slope": beta, "slope_se": se},89 {"period": "sub1_1960Q1-1983Q4", "n": len(s1), "intercept": float(f1.params[0]),90 "slope": float(f1.params[1]), "slope_se": hac_se(f1, nw_lags(len(s1)))},91 {"period": "sub2_1984Q1-2019Q4", "n": len(s2), "intercept": float(f2.params[0]),92 "slope": float(f2.params[1]), "slope_se": hac_se(f2, nw_lags(len(s2)))},93])94table.to_csv(f"{OUT}/table1.csv", index=False)9596# robustness (supplementary): HAC Wald tests in the interaction model du = a + b g + c D + d D*g97D = (d.index >= SPLIT_Q).astype(float)98Xi = np.column_stack([np.ones(n), d["g"].values, D, D * d["g"].values])99fi = sm.OLS(d["du"].values, Xi).fit().get_robustcov_results(cov_type="HAC", maxlags=hac_lags,100 kernel="bartlett", use_correction=False)101w_joint = fi.wald_test(np.array([[0, 0, 1, 0], [0, 0, 0, 1]]), use_f=False, scalar=True)102w_slope = fi.wald_test(np.array([[0, 0, 0, 1]]), use_f=False, scalar=True)103w_int = fi.wald_test(np.array([[0, 0, 1, 0]]), use_f=False, scalar=True)104105sens = {str(L): hac_se(full, L) for L in (0, 2, 4, 8, 12)}106concl = "positive_effect" if ci_lo > 0 else ("negative_effect" if ci_hi < 0 else "no_detectable_effect")107108res = {109 "task_id": "F6-conf-03",110 "feasibility": "complete",111 "feasibility_reason": None,112 "data_issues": data_issues,113 "n_input_rows": n_input_rows,114 "n_analysed": int(n),115 "random_seed": None,116 "sample_start": d.index[0].start_time.strftime("%Y-%m-%d"),117 "sample_end": d.index[-1].start_time.strftime("%Y-%m-%d"),118 "n_sub1": int(len(s1)), "n_sub2": int(len(s2)),119 "baseline_metric": baseline_rmse,120 "model_metric": model_rmse,121 "effect_estimate": beta,122 "intercept": alpha,123 "effect_se": se,124 "hac_lags": hac_lags,125 "effect_ci_low": float(ci_lo), "effect_ci_high": float(ci_hi),126 "beta_sub1": float(f1.params[1]), "beta_sub2": float(f2.params[1]),127 "chow_f_stat": chow_f, "chow_p_value": chow_p,128 "stability_conclusion": "stable" if chow_p >= 0.05 else "unstable",129 "conclusion": concl,130 # supplementary keys131 "ols_classical_se": float(full.bse[1]),132 "r_squared": float(full.rsquared),133 "hac_se_by_lag": sens,134 "intercept_sub1": float(f1.params[0]), "intercept_sub2": float(f2.params[0]),135 "beta_sub1_hac_se": float(table.loc[1, "slope_se"]), "beta_sub2_hac_se": float(table.loc[2, "slope_se"]),136 "hac_lags_sub1": nw_lags(len(s1)), "hac_lags_sub2": nw_lags(len(s2)),137 "chow_df": [k, int(n - 2 * k)],138 "hac_wald_joint_p": float(w_joint.pvalue), "hac_wald_slope_diff_p": float(w_slope.pvalue),139 "hac_wald_intercept_diff_p": float(w_int.pvalue),140 "okun_growth_sub1": float(-f1.params[0] / f1.params[1]), "okun_growth_sub2": float(-f2.params[0] / f2.params[1]),141 "n_window_quarters": int(n0),142 "mean_du": float(d["du"].mean()), "mean_g": float(d["g"].mean()),143}144with open(f"{OUT}/results.json", "w") as fh:145 json.dump(res, fh, indent=2)146147fig, ax = plt.subplots(figsize=(7, 5))148ax.scatter(s1["g"], s1["du"], s=14, alpha=0.7, label="1960Q1–1983Q4")149ax.scatter(s2["g"], s2["du"], s=14, alpha=0.7, marker="^", label="1984Q1–2019Q4")150for sub, fit, ls in ((s1, f1, "-"), (s2, f2, "--")):151 xs = np.linspace(sub["g"].min(), sub["g"].max(), 50)152 ax.plot(xs, fit.params[0] + fit.params[1] * xs, ls, color="k",153 label=f"OLS {sub.index[0]}–{sub.index[-1]}: β={fit.params[1]:.3f}")154ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)155ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)")156ax.set_ylabel("Change in unemployment rate Δu_t (pp)")157ax.set_title("Okun's law, US 1960Q1–2019Q4")158ax.legend(fontsize=8)159fig.tight_layout()160fig.savefig(f"{OUT}/figure1.png", dpi=150)161print(json.dumps({k: res[k] for k in ("effect_estimate", "effect_se", "hac_lags", "effect_ci_low",162 "effect_ci_high", "beta_sub1", "beta_sub2", "chow_f_stat",163 "chow_p_value", "baseline_metric", "model_metric", "conclusion",164 "stability_conclusion")}, indent=1))165