1"""Okun's law, U.S. 1960Q1-2019Q4, OLS with Newey-West HAC SE. Run from workspace root."""2import json, os3import numpy as np, pandas as pd4import statsmodels.api as sm5from scipy import stats6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt910OUT = "out"11FIRST, LAST, BASE = pd.Period("1960Q1"), pd.Period("2019Q4"), pd.Period("1959Q4")12issues = []1314un = pd.read_csv("data/UNRATE.csv", dtype={"observation_date": str})15gd = pd.read_csv("data/GDPC1.csv", dtype={"observation_date": str})16n_rows = {"UNRATE.csv": int(len(un)), "GDPC1.csv": int(len(gd))}17un["UNRATE"] = pd.to_numeric(un["UNRATE"], errors="coerce")18gd["GDPC1"] = pd.to_numeric(gd["GDPC1"], errors="coerce")19for d in (un, gd):20 d["date"] = pd.to_datetime(d["observation_date"])21 assert not d["date"].duplicated().any()22miss = un.loc[un["UNRATE"].isna(), "observation_date"].tolist()23if miss:24 issues.append({"issue": f"UNRATE missing value(s) at {', '.join(miss)}",25 "handling": "outside the 1959Q4-2019Q4 window; not used. Quarterly averages require all 3 months (asserted for every quarter used)."})2627un["q"] = un["date"].dt.to_period("Q")28uq = un.groupby("q")["UNRATE"].agg(["mean", "count"])29gd["q"] = gd["date"].dt.to_period("Q")30gq = gd.set_index("q")["GDPC1"]3132n0 = len(gq)33panel = pd.DataFrame({"gdp": gq}).join(uq, how="left")34panel = panel.loc[(panel.index >= BASE) & (panel.index <= LAST)]35assert (panel["count"] == 3).all() and panel["gdp"].notna().all() and (panel["gdp"] > 0).all()36n1 = len(panel)37panel["du"] = panel["mean"].diff()38panel["g"] = 100 * np.log(panel["gdp"] / panel["gdp"].shift(1))39df = panel.loc[panel.index >= FIRST].copy()40n2 = len(df)41assert df[["du", "g"]].notna().all().all()42assert n2 == 240 and df.index[0] == FIRST and df.index[-1] == LAST4344pd.DataFrame([45 dict(step=1, rule="keep quarters 1959Q4-2019Q4 (GDPC1 quarters with complete 3-month UNRATE average)",46 n_before=n0, n_after=n1, n_removed=n0 - n1, rationale="window fixed by TASK.md; 1959Q4 kept only as difference base"),47 dict(step=2, rule="drop 1959Q4 (no lagged level for first difference)", n_before=n1, n_after=n2,48 n_removed=n1 - n2, rationale="first differences need t-1 levels; regression sample 1960Q1-2019Q4"),49]).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)5051y, X = df["du"].values, sm.add_constant(df["g"].values)52T = len(y)53ols = sm.OLS(y, X).fit()54# Bandwidth: Newey-West (1994) plug-in-free rule floor(4*(T/100)^(2/9)) -> 4 for T=24055L = int(np.floor(4 * (T / 100) ** (2 / 9)))56def hac(lags):57 # classic Newey-West: Bartlett weights 1 - j/(L+1), no small-sample df correction58 return sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": lags, "use_correction": False})59h = hac(L)60# manual cross-check of the HAC SE61e = ols.resid; Xe = X * e[:, None]; S = Xe.T @ Xe62for j in range(1, L + 1):63 G = Xe[j:].T @ Xe[:-j]; S += (1 - j / (L + 1)) * (G + G.T)64XXi = np.linalg.inv(X.T @ X); se_manual = np.sqrt((XXi @ S @ XXi)[1, 1])65beta, se_h = float(ols.params[1]), float(h.bse[1])66assert abs(se_manual - se_h) < 1e-1067z = stats.norm.ppf(0.975)68lo, hi = beta - z * se_h, beta + z * se_h69p_hac = float(2 * stats.norm.sf(abs(beta / se_h)))70resid_ac = [float(pd.Series(e).autocorr(k)) for k in (1, 2, 3, 4)]7172res = {73 "task_id": "F2-conf-04",74 "beta_hat": beta, "alpha_hat": float(ols.params[0]), "r_squared": float(ols.rsquared),75 "se_iid": float(ols.bse[1]), "se_hac": se_h, "hac_lags": L,76 "beta_ci_low": float(lo), "beta_ci_high": float(hi),77 "okun_negative_5pct": bool(beta < 0 and p_hac < 0.05),78 "mean_growth": float(df["g"].mean()),79 "first_quarter": str(df.index[0]), "last_quarter": str(df.index[-1]),80 "n_analysed": int(T), "feasibility": "complete", "feasibility_reason": None,81 "data_issues": issues, "n_input_rows": n_rows, "random_seed": None,82 # supplementary83 "t_hac": float(beta / se_h), "p_value_hac": p_hac,84 "baseline_ci_iid_low": float(beta - z * ols.bse[1]), "baseline_ci_iid_high": float(beta + z * ols.bse[1]),85 "resid_autocorr_lag1_4": resid_ac,86 "sensitivity_se_hac_by_lag": {str(k): float(hac(k).bse[1]) for k in (0, 2, 4, 5, 8, 12)},87 "mean_du": float(df["du"].mean()),88}89with open(os.path.join(OUT, "results.json"), "w") as f:90 json.dump(res, f, indent=2)9192fig, ax = plt.subplots(figsize=(6.5, 4.5))93ax.scatter(df["g"], df["du"], s=12, alpha=0.6)94gg = np.linspace(df["g"].min(), df["g"].max(), 50)95ax.plot(gg, res["alpha_hat"] + beta * gg, color="C3",96 label=f"OLS: Δu = {res['alpha_hat']:.3f} + ({beta:.3f}) g")97ax.axhline(0, lw=0.5, color="k"); ax.axvline(0, lw=0.5, color="k")98ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)")99ax.set_ylabel("Δ unemployment rate Δu_t (pp)")100ax.set_title("Okun's law, U.S. 1960Q1–2019Q4 (n=240)")101ax.legend(); fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)102print(json.dumps({k: res[k] for k in ("beta_hat", "se_iid", "se_hac", "hac_lags", "beta_ci_low", "beta_ci_high", "okun_negative_5pct")}, indent=1))103