1"""Okun's law, US 1960Q1-2019Q4, OLS with Newey-West HAC SE. Run from workspace root."""2import json, math3import numpy as np, pandas as pd, statsmodels.api as sm4from scipy import stats5import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt67u = pd.read_csv("data/UNRATE.csv"); g = pd.read_csv("data/GDPC1.csv")8n_rows = {"UNRATE.csv": int(len(u)), "GDPC1.csv": int(len(g))}9issues = []10u["date"] = pd.to_datetime(u.observation_date); g["date"] = pd.to_datetime(g.observation_date)11miss = u[u.UNRATE.isna()]12if len(miss):13 issues.append({"issue": f"UNRATE has {len(miss)} missing monthly value(s): " + ", ".join(miss.observation_date),14 "handling": "Outside the 1959Q4-2019Q4 levels window; not used. Quarterly averages inside the window are checked to have all 3 months."})15u["q"] = u.date.dt.to_period("Q"); g["q"] = g.date.dt.to_period("Q")16uq = u.groupby("q").UNRATE.agg(["mean", "count"])17gq = g.set_index("q").GDPC118win = pd.period_range("1959Q4", "2019Q4", freq="Q")19assert (uq.loc[win, "count"] == 3).all(), "incomplete quarter in window"20assert gq.loc[win].notna().all() and len(gq.loc[win]) == len(win)21df = pd.DataFrame({"u": uq.loc[win, "mean"], "y": gq.loc[win]})22df["du"] = df.u.diff(); df["g"] = 100 * np.log(df.y / df.y.shift(1))23df = df.loc[pd.Period("1960Q1"):].dropna()24T = len(df); assert T == 2402526X = sm.add_constant(df.g.values); y = df.du.values27ols = sm.OLS(y, X).fit()28# Bandwidth: Newey-West (1994) plug-in rule of thumb floor(4*(T/100)^(2/9))29L = int(math.floor(4 * (T / 100) ** (2 / 9)))30hac = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": False})31b, a = ols.params[1], ols.params[0]32se_hac = float(hac.bse[1]); z = stats.norm.ppf(0.975)33lo, hi = b - z * se_hac, b + z * se_hac34p_hac = 2 * stats.norm.sf(abs(b / se_hac))3536# sensitivity (reported in report, not a required key)37sens = {}38for lag in [0, 2, 4, 6, 8, 12]:39 for corr in [False, True]:40 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": lag, "use_correction": corr})41 sens[f"L{lag}_{'dfcorr' if corr else 'nocorr'}"] = float(r.bse[1])42resid_ac = [float(pd.Series(ols.resid).autocorr(k)) for k in range(1, 5)]4344res = {45 "task_id": "F2-conf-04", "beta_hat": float(b), "alpha_hat": float(a), "r_squared": float(ols.rsquared),46 "se_iid": float(ols.bse[1]), "se_hac": se_hac, "hac_lags": L,47 "beta_ci_low": float(lo), "beta_ci_high": float(hi),48 "okun_negative_5pct": bool(b < 0 and p_hac < 0.05),49 "mean_growth": float(df.g.mean()), "first_quarter": str(df.index[0]), "last_quarter": str(df.index[-1]),50 "n_analysed": int(T), "feasibility": "complete", "feasibility_reason": None,51 "data_issues": issues, "n_input_rows": n_rows, "random_seed": None,52}53with open("out/results.json", "w") as f: json.dump(res, f, indent=2)54with open("out/diagnostics.json", "w") as f:55 json.dump({"p_hac": float(p_hac), "t_hac": float(b / se_hac), "se_sensitivity": sens,56 "resid_autocorr_1to4": resid_ac, "alpha_se_hac": float(hac.bse[0])}, f, indent=2)5758fig, ax = plt.subplots(figsize=(6, 4.5))59ax.scatter(df.g, df.du, s=12, alpha=0.6)60xs = np.linspace(df.g.min(), df.g.max(), 50); ax.plot(xs, a + b * xs, "r-",61 label=f"OLS: Δu = {a:.3f} {b:+.3f}·g")62ax.set_xlabel("Real GDP growth g_t (%, q/q, not annualised)"); ax.set_ylabel("Δu_t (pp)")63ax.set_title("Okun's law, US 1960Q1–2019Q4 (T=240)"); ax.legend(); fig.tight_layout()64fig.savefig("out/figure1.png", dpi=150)65print(json.dumps(res, indent=2))66