1"""Okun's law, U.S. 1960Q1-2019Q4, OLS with Newey-West HAC SE. 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 plt910issues = []11u = pd.read_csv("data/UNRATE.csv", parse_dates=["observation_date"])12g = pd.read_csv("data/GDPC1.csv", parse_dates=["observation_date"])13n_rows = {"UNRATE.csv": int(len(u)), "GDPC1.csv": int(len(g))}1415for name, df, col in [("UNRATE.csv", u, "UNRATE"), ("GDPC1.csv", g, "GDPC1")]:16 if df["observation_date"].duplicated().any():17 raise SystemExit(f"duplicate dates in {name}")18 miss = df[df[col].isna()]19 if len(miss):20 dates = ", ".join(d.strftime("%Y-%m") for d in miss["observation_date"])21 issues.append({"issue": f"{name}: {len(miss)} missing value(s) at {dates}",22 "handling": "Outside the 1959Q4-2019Q4 window needed; dropped, no effect on the analysis."})2324u = u.dropna()25u["q"] = u["observation_date"].dt.to_period("Q")26cnt = u.groupby("q")["UNRATE"].count()27uq = u.groupby("q")["UNRATE"].mean()28g["q"] = g["observation_date"].dt.to_period("Q")29gq = g.dropna().set_index("q")["GDPC1"]3031win = pd.period_range("1959Q4", "2019Q4", freq="Q")32assert (cnt.reindex(win) == 3).all(), "incomplete months in window"33assert gq.reindex(win).notna().all(), "missing GDP in window"34uw, gw = uq.reindex(win), gq.reindex(win)35du = uw.diff().iloc[1:]36gr = (100 * np.log(gw)).diff().iloc[1:]37assert len(du) == 2403839X = sm.add_constant(gr.values)40y = du.values41T = len(y)42L = int(math.floor(4 * (T / 100) ** (2 / 9))) # Newey-West (1994) plug-in rule of thumb -> 443ols = sm.OLS(y, X).fit()44hac = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": False})45beta, alpha = float(ols.params[1]), float(ols.params[0])46se_hac = float(hac.bse[1])47z = stats.norm.ppf(0.975)48lo, hi = beta - z * se_hac, beta + z * se_hac49p_hac = float(2 * stats.norm.sf(abs(beta / se_hac)))5051# sensitivity (printed only)52sens = {}53for lag in [0, 1, 2, 4, 6, 8, 12, 16]:54 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": lag, "use_correction": False})55 sens[lag] = float(r.bse[1])56resid_ac = [float(pd.Series(ols.resid).autocorr(k)) for k in range(1, 5)]5758res = {59 "task_id": "F2-conf-04",60 "beta_hat": beta, "alpha_hat": alpha, "r_squared": float(ols.rsquared),61 "se_iid": float(ols.bse[1]), "se_hac": se_hac, "hac_lags": L,62 "beta_ci_low": float(lo), "beta_ci_high": float(hi),63 "okun_negative_5pct": bool(beta < 0 and p_hac < 0.05),64 "mean_growth": float(gr.mean()),65 "first_quarter": str(du.index[0]), "last_quarter": str(du.index[-1]),66 "n_analysed": int(T), "feasibility": "complete", "feasibility_reason": None,67 "data_issues": issues, "n_input_rows": n_rows, "random_seed": None,68}69with open("out/results.json", "w") as f:70 json.dump(res, f, indent=2)7172print(json.dumps(res, indent=2))73print("p_hac", p_hac, "t_hac", beta / se_hac)74print("HAC SE by lag:", sens)75print("resid autocorr 1-4:", resid_ac)76r_corr = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": True})77print("HAC SE with df correction:", float(r_corr.bse[1]))7879fig, ax = plt.subplots(figsize=(6, 4.5))80ax.scatter(gr.values, y, s=12, alpha=0.7)81xs = np.linspace(gr.min(), gr.max(), 50)82ax.plot(xs, alpha + beta * xs, color="red",83 label=f"OLS: Δu = {alpha:.3f} {beta:+.3f}·g")84ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)85ax.set_xlabel("Real GDP growth g_t (%, q/q, not annualised)")86ax.set_ylabel("Δu_t (pp)")87ax.set_title("Okun's law, U.S. 1960Q1–2019Q4 (n=240)")88ax.legend(); fig.tight_layout(); fig.savefig("out/figure1.png", dpi=150)89