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 plt910FIRST, LAST = pd.Period("1960Q1", "Q"), pd.Period("2019Q4", "Q")11issues, excl = [], []1213u_raw = pd.read_csv("data/UNRATE.csv")14g_raw = pd.read_csv("data/GDPC1.csv")15n_rows = {"UNRATE.csv": int(len(u_raw)), "GDPC1.csv": int(len(g_raw))}1617# --- unemployment: monthly -> quarterly average18u = u_raw.copy()19u["date"] = pd.to_datetime(u["observation_date"])20u["UNRATE"] = pd.to_numeric(u["UNRATE"], errors="coerce")21miss = u.loc[u["UNRATE"].isna(), "observation_date"].tolist()22if miss:23 issues.append({"issue": f"UNRATE missing value(s) at {', '.join(miss)}",24 "handling": "outside the 1959Q4-2019Q4 span needed; dropped by window filter, no effect on estimates"})25u["q"] = u["date"].dt.to_period("Q")26# quarterly average only for quarters with all 3 monthly values present27cnt_all = u.groupby("q")["UNRATE"].count()28uq_all = u.groupby("q")["UNRATE"].mean()[cnt_all == 3]2930# --- GDP31g = g_raw.copy()32g["q"] = pd.to_datetime(g["observation_date"]).dt.to_period("Q")33g["GDPC1"] = pd.to_numeric(g["GDPC1"], errors="coerce")34assert g["GDPC1"].notna().all() and not g["q"].duplicated().any()35gq = g.set_index("q")["GDPC1"]3637# Exclusion ledger in quarterly units, starting from the 319 GDPC1 quarters38panel = pd.DataFrame({"y": gq})39n0 = len(panel)40panel["u"] = uq_all.reindex(panel.index)41panel = panel.dropna()42excl.append(("merge_unrate", "keep GDPC1 quarters with a complete (3-month, non-missing) UNRATE quarterly average",43 n0, len(panel), n0 - len(panel),44 "UNRATE (945 monthly rows) averaged to quarters; 1947 has no UNRATE, 2025Q4 lacks Oct-2025 value"))45n0 = len(panel)46panel = panel.loc[FIRST - 1:LAST]47excl.append(("window", "keep 1959Q4-2019Q4 (1959Q4 only for differencing)", n0, len(panel), n0 - len(panel),48 "information rule: fixed sample window, no extension"))49assert len(panel) == 241 and (cnt_all.loc[FIRST - 1:LAST] == 3).all()50df = panel.sort_index()51df["du"] = df["u"].diff()52df["g"] = 100 * np.log(df["y"] / df["y"].shift(1))53n0 = len(df)54df = df.loc[FIRST:LAST].dropna()55excl.append(("differencing", "drop 1959Q4 (levels used only to form 1960Q1 differences)", n0, len(df), n0 - len(df),56 "first difference unavailable for base quarter"))57T = len(df)58assert T == 2405960X = sm.add_constant(df["g"].values)61y = df["du"].values62ols = sm.OLS(y, X).fit()63# Newey-West (1994) plug-in rule of thumb: floor(4*(T/100)^(2/9))64L = int(math.floor(4 * (T / 100) ** (2 / 9)))65hac = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": False})6667# manual check of HAC SE68e = ols.resid; Xe = X * e[:, None]69S = Xe.T @ Xe70for j in range(1, L + 1):71 w = 1 - j / (L + 1); G = Xe[j:].T @ Xe[:-j]; S += w * (G + G.T)72XtXi = np.linalg.inv(X.T @ X)73se_manual = float(np.sqrt((XtXi @ S @ XtXi)[1, 1]))74assert abs(se_manual - hac.bse[1]) < 1e-107576beta, alpha = float(ols.params[1]), float(ols.params[0])77se_hac = float(hac.bse[1]); z = stats.norm.ppf(0.975)78lo, hi = beta - z * se_hac, beta + z * se_hac79p_hac = float(2 * stats.norm.sf(abs(beta / se_hac)))8081sens = {}82for l in [0, 1, 2, 4, 6, 8, 12]:83 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": l, "kernel": "bartlett", "use_correction": False})84 sens[str(l)] = float(r.bse[1])8586res = {87 "task_id": "F2-conf-04",88 "beta_hat": beta, "alpha_hat": alpha, "r_squared": float(ols.rsquared),89 "se_iid": float(ols.bse[1]), "se_hac": se_hac, "hac_lags": L,90 "beta_ci_low": float(lo), "beta_ci_high": float(hi),91 "okun_negative_5pct": bool(beta < 0 and p_hac < 0.05),92 "mean_growth": float(df["g"].mean()),93 "first_quarter": str(df.index[0]), "last_quarter": str(df.index[-1]),94 "n_analysed": int(T),95 "feasibility": "complete", "feasibility_reason": None,96 "data_issues": issues, "n_input_rows": n_rows, "random_seed": None,97 # supplementary98 "t_stat_hac": float(beta / se_hac), "p_value_hac": p_hac,99 "residual_acf_lag1": float(pd.Series(e).autocorr(1)),100 "se_hac_by_lags": sens,101 "beta_ci_iid_low": float(beta - z * ols.bse[1]), "beta_ci_iid_high": float(beta + z * ols.bse[1]),102}103with open("out/results.json", "w") as f:104 json.dump(res, f, indent=2)105pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv("out/exclusions.csv", index=False)106107fig, ax = plt.subplots(figsize=(6.5, 4.5))108ax.scatter(df["g"], df["du"], s=12, alpha=0.6)109xs = np.linspace(df["g"].min(), df["g"].max(), 100)110ax.plot(xs, alpha + beta * xs, color="red", label=f"OLS: du = {alpha:.3f} {beta:+.3f} g")111ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)112ax.set_xlabel("Real GDP growth g_t (%, q/q, not annualised)")113ax.set_ylabel("Change in unemployment rate du_t (pp)")114ax.set_title("Okun's law, U.S. 1960Q1-2019Q4 (n=240)"); ax.legend()115fig.tight_layout(); fig.savefig("out/figure1.png", dpi=120)116print(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))117