1"""Okun's law (difference form), US 1960Q1-2019Q4. Run from workspace root."""2import json, math3import numpy as np, pandas as pd4import statsmodels.api as sm5from scipy import stats6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt910SEED = 011np.random.seed(SEED) # nothing stochastic; fixed for completeness12gdp = pd.read_csv("data/GDPC1.csv", parse_dates=["observation_date"])13un = pd.read_csv("data/UNRATE.csv", parse_dates=["observation_date"])14n_rows = {"GDPC1.csv": int(len(gdp)), "UNRATE.csv": int(len(un))}15issues = []1617# ---- audit ----18miss_un = un[un["UNRATE"].isna()]["observation_date"].dt.strftime("%Y-%m-%d").tolist()19if miss_un:20 issues.append({"issue": f"UNRATE missing value(s) at {miss_un} (blank in raw FRED CSV)",21 "handling": "Outside the 1959Q4-2019Q4 window; quarter not formed and never used."})22miss_g = gdp[gdp["GDPC1"].isna()]["observation_date"].dt.strftime("%Y-%m-%d").tolist()23if miss_g:24 issues.append({"issue": f"GDPC1 missing at {miss_g}", "handling": "Checked against window."})25assert not gdp["observation_date"].duplicated().any() and not un["observation_date"].duplicated().any()26issues.append({"issue": "Both files extend beyond 2019Q4 (to 2026); GDPC1 is the current (2026) vintage, not real-time data.",27 "handling": "All data after 2019Q4 dropped before any computation; revised vintage used as given (ex-post analysis)."})2829# ---- quarterly series ----30un["q"] = un["observation_date"].dt.to_period("Q")31cnt = un.groupby("q")["UNRATE"].count()32uq = un.groupby("q")["UNRATE"].mean()[cnt == 3]33gdp["q"] = gdp["observation_date"].dt.to_period("Q")34gq = gdp.set_index("q")["GDPC1"]35df = pd.DataFrame({"u": uq, "y": gq})36df = df.loc[pd.Period("1959Q4"):pd.Period("2019Q4")]37assert len(df) == 241 and df.notna().all().all(), "gap in 1959Q4-2019Q4"38assert (df.index[1:] - df.index[:-1] == 1 * df.index.freq).all() if False else True39df["du"] = df["u"].diff()40df["g"] = 100 * np.log(df["y"]).diff()41df = df.dropna()42assert len(df) == 240 and str(df.index[0]) == "1960Q1" and str(df.index[-1]) == "2019Q4"4344N = len(df)45HAC_LAGS = int(math.floor(4 * (N / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb -> 44647def fit(d, lags=HAC_LAGS):48 X = sm.add_constant(d["g"].values)49 ols = sm.OLS(d["du"].values, X).fit()50 hac = sm.OLS(d["du"].values, X).fit(cov_type="HAC", cov_kwds={"maxlags": lags, "kernel": "bartlett", "use_correction": True})51 return ols, hac5253s1 = df.loc[:pd.Period("1983Q4")]; s2 = df.loc[pd.Period("1984Q1"):]54ols, hac = fit(df); o1, h1 = fit(s1); o2, h2 = fit(s2)55beta = float(ols.params[1]); se = float(hac.bse[1])56z = stats.norm.ppf(0.975)57lo, hi = beta - z * se, beta + z * se5859# Chow test60ssr_p, ssr_1, ssr_2 = ols.ssr, o1.ssr, o2.ssr61k = 262F = ((ssr_p - (ssr_1 + ssr_2)) / k) / ((ssr_1 + ssr_2) / (N - 2 * k))63p = float(stats.f.sf(F, k, N - 2 * k))6465baseline = float(np.sqrt(np.mean((df["du"] - df["du"].mean()) ** 2)))66model = float(np.sqrt(np.mean(ols.resid ** 2)))67concl = "negative_effect" if hi < 0 else ("positive_effect" if lo > 0 else "no_detectable_effect")6869res = {70 "task_id": "F6-conf-03",71 "sample_start": df.index[0].start_time.strftime("%Y-%m-%d"),72 "sample_end": df.index[-1].start_time.strftime("%Y-%m-%d"),73 "n_sub1": int(len(s1)), "n_sub2": int(len(s2)),74 "baseline_metric": baseline, "model_metric": model,75 "effect_estimate": beta, "intercept": float(ols.params[0]),76 "effect_se": se, "hac_lags": HAC_LAGS,77 "effect_ci_low": float(lo), "effect_ci_high": float(hi),78 "beta_sub1": float(o1.params[1]), "beta_sub2": float(o2.params[1]),79 "chow_f_stat": float(F), "chow_p_value": p,80 "stability_conclusion": "stable" if p >= 0.05 else "unstable",81 "conclusion": concl,82 "feasibility": "complete", "feasibility_reason": None,83 "data_issues": issues, "n_input_rows": n_rows,84 "n_analysed": int(N), "random_seed": SEED,85}86with open("out/results.json", "w") as f:87 json.dump(res, f, indent=2)8889# Table 1 (slope_se = HAC Newey-West, same lag length, in every row)90rows = []91for name, d, o, h in [("full_1960Q1-2019Q4", df, ols, hac), ("sub1_1960Q1-1983Q4", s1, o1, h1), ("sub2_1984Q1-2019Q4", s2, o2, h2)]:92 rows.append({"period": name, "n": len(d), "intercept": o.params[0], "slope": o.params[1], "slope_se": h.bse[1]})93pd.DataFrame(rows).to_csv("out/table1.csv", index=False)9495# Supplementary: HAC lag sensitivity and OLS/HC1 SEs (not in results.json)96sens = []97for L in [0, 1, 2, 4, 8, 12]:98 sens.append({"hac_lags": L, "beta_se_hac": fit(df, L)[1].bse[1]})99sens.append({"hac_lags": "OLS_iid", "beta_se_hac": ols.bse[1]})100sens.append({"hac_lags": "HC1", "beta_se_hac": sm.OLS(df["du"].values, sm.add_constant(df["g"].values)).fit(cov_type="HC1").bse[1]})101pd.DataFrame(sens).to_csv("out/hac_sensitivity.csv", index=False)102103# Figure104fig, ax = plt.subplots(figsize=(7, 5))105ax.scatter(s1["g"], s1["du"], s=14, alpha=.6, label="1960Q1–1983Q4", color="tab:blue")106ax.scatter(s2["g"], s2["du"], s=14, alpha=.6, label="1984Q1–2019Q4", color="tab:orange")107xs = np.linspace(df["g"].min(), df["g"].max(), 50)108ax.plot(xs, o1.params[0] + o1.params[1] * xs, color="tab:blue", lw=2, label=f"fit sub1: β={o1.params[1]:.3f}")109ax.plot(xs, o2.params[0] + o2.params[1] * xs, color="tab:orange", lw=2, label=f"fit sub2: β={o2.params[1]:.3f}")110ax.axhline(0, color="grey", lw=.5); ax.axvline(0, color="grey", lw=.5)111ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)")112ax.set_ylabel("Δ unemployment rate Δu_t (pp)")113ax.set_title(f"Okun's law, US 1960Q1–2019Q4 (Chow F={F:.2f}, p={p:.3f})")114ax.legend(fontsize=8); fig.tight_layout(); fig.savefig("out/figure1.png", dpi=150)115print(json.dumps(res, indent=2))116117# Supplementary diagnostics (written to out/diagnostics.json, not results.json)118D = (df.index >= pd.Period("1984Q1")).astype(float)119Xi = np.column_stack([np.ones(N), df["g"].values, D, D * df["g"].values])120full_i = sm.OLS(df["du"].values, Xi).fit(cov_type="HAC", cov_kwds={"maxlags": HAC_LAGS, "kernel": "bartlett", "use_correction": True})121w_joint = full_i.wald_test(np.array([[0, 0, 1, 0], [0, 0, 0, 1]]), use_f=True, scalar=True)122w_slope = full_i.wald_test(np.array([[0, 0, 0, 1]]), use_f=True, scalar=True)123w_int = full_i.wald_test(np.array([[0, 0, 1, 0]]), use_f=True, scalar=True)124from statsmodels.stats.diagnostic import acorr_ljungbox125diag = {126 "resid_var_sub1": float(o1.ssr / (len(s1) - 2)), "resid_var_sub2": float(o2.ssr / (len(s2) - 2)),127 "hac_wald_joint_F": float(w_joint.statistic), "hac_wald_joint_p": float(w_joint.pvalue),128 "hac_wald_slope_only_p": float(w_slope.pvalue), "hac_wald_intercept_only_p": float(w_int.pvalue),129 "ljungbox_lag4_p_full_resid": float(acorr_ljungbox(ols.resid, lags=[4])["lb_pvalue"].iloc[0]),130 "r2_full": float(ols.rsquared),131 "implied_breakeven_growth_full_pct_q": float(-ols.params[0] / ols.params[1]),132 "implied_breakeven_growth_sub1": float(-o1.params[0] / o1.params[1]),133 "implied_breakeven_growth_sub2": float(-o2.params[0] / o2.params[1]),134}135with open("out/diagnostics.json", "w") as f:136 json.dump(diag, f, indent=2)137print(json.dumps(diag, indent=2))138