1"""Okun's law, U.S. 1960Q1-2019Q4, OLS with Newey-West HAC SE. Reads only data/, writes out/."""2import json, hashlib, os3import numpy as np, pandas as pd4import matplotlib5matplotlib.use("Agg")6import matplotlib.pyplot as plt7from scipy import stats89ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))10D, O = os.path.join(ROOT, "data"), os.path.join(ROOT, "out")11FIRST, LAST = pd.Period("1960Q1", "Q"), pd.Period("2019Q4", "Q")1213def sha(p):14 return hashlib.sha256(open(p, "rb").read()).hexdigest()1516un = pd.read_csv(os.path.join(D, "UNRATE.csv"))17gd = pd.read_csv(os.path.join(D, "GDPC1.csv"))18n_rows = {"UNRATE.csv": int(len(un)), "GDPC1.csv": int(len(gd))}19data_issues = []20miss = un[un["UNRATE"].isna()]21for d in miss["observation_date"]:22 data_issues.append({"issue": f"UNRATE missing (empty) value at {d}",23 "handling": "Quarter containing it has <3 monthly values and is dropped; it lies outside the 1959Q4-2019Q4 levels window, so results are unaffected."})2425# quarterly average unemployment, only quarters with all 3 months26un["q"] = pd.to_datetime(un["observation_date"]).dt.to_period("Q")27g = un.groupby("q")["UNRATE"]28uq = g.mean()[g.count() == 3].rename("u")29n_un_q = int(un["q"].nunique())30gd["q"] = pd.to_datetime(gd["observation_date"]).dt.to_period("Q")31gq = gd.set_index("q")["GDPC1"].astype(float)3233ledger = []34def log(step, rule, nb, na, why):35 ledger.append(dict(step=step, rule=rule, n_before=nb, n_after=na, n_removed=nb - na, rationale=why))3637df = pd.concat([gq, uq], axis=1, join="inner").sort_index()38log(1, "keep GDPC1 quarters that have a complete (3-month) UNRATE quarterly average", len(gq), len(df),39 "UNRATE starts 1948Q1 (GDPC1 starts 1947Q1); 2025Q4 lacks October 2025 UNRATE")40# differences only between consecutive calendar quarters41idx = df.index42consec = pd.Series([False] + [idx[i].ordinal - idx[i - 1].ordinal == 1 for i in range(1, len(idx))], index=idx)43df["du"] = df["u"].diff()44df["g"] = 100 * np.log(df["GDPC1"] / df["GDPC1"].shift(1))45df.loc[~consec, ["du", "g"]] = np.nan46nb = len(df); df = df.dropna(subset=["du", "g"])47log(2, "drop quarters whose previous calendar quarter is unavailable (first differences undefined)", nb, len(df),48 "first available quarter and quarters following a dropped quarter")49nb = len(df); df = df[(df.index >= FIRST) & (df.index <= LAST)]50log(3, "restrict to regression window 1960Q1-2019Q4", nb, len(df), "task-defined sample; no extension of window")51assert len(df) == 240 and df.index[0] == FIRST and df.index[-1] == LAST52pd.DataFrame(ledger).to_csv(os.path.join(O, "exclusions.csv"), index=False)5354# OLS55y = df["du"].to_numpy(); x = df["g"].to_numpy(); T = len(y)56X = np.column_stack([np.ones(T), x])57XtXi = np.linalg.inv(X.T @ X)58b = XtXi @ X.T @ y59e = y - X @ b60k = X.shape[1]61r2 = 1 - (e @ e) / np.sum((y - y.mean()) ** 2)62se_iid = float(np.sqrt((e @ e) / (T - k) * XtXi[1, 1]))6364def nw_se(L):65 """Newey-West (1987) Bartlett-kernel HAC, no small-sample df correction."""66 u = X * e[:, None]67 S = u.T @ u68 for l in range(1, L + 1):69 w = 1 - l / (L + 1)70 G = u[l:].T @ u[:-l]71 S += w * (G + G.T)72 V = XtXi @ S @ XtXi73 return float(np.sqrt(V[1, 1]))7475# Bandwidth: Newey-West (1994) / Stock-Watson rule of thumb floor(4*(T/100)^(2/9))76L = int(np.floor(4 * (T / 100) ** (2 / 9)))77se_hac = nw_se(L)78z = stats.norm.ppf(0.975)79lo, hi = float(b[1] - z * se_hac), float(b[1] + z * se_hac)80t_hac = float(b[1] / se_hac)81p_hac = float(2 * stats.norm.sf(abs(t_hac)))82neg = bool(b[1] < 0 and p_hac < 0.05)8384# cross-check with statsmodels (no df correction)85import statsmodels.api as sm86m = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": False})87assert abs(m.bse[1] - se_hac) < 1e-10 and abs(m.params[1] - b[1]) < 1e-1288m0 = sm.OLS(y, X).fit()89assert abs(m0.bse[1] - se_iid) < 1e-109091sens = {str(l): nw_se(l) for l in [0, 1, 2, 4, 6, 8, 12]}92rho1 = float(np.corrcoef(e[1:], e[:-1])[0, 1])9394res = {95 "task_id": "F2-conf-04",96 "beta_hat": float(b[1]), "alpha_hat": float(b[0]), "r_squared": float(r2),97 "se_iid": se_iid, "se_hac": se_hac, "hac_lags": L,98 "beta_ci_low": lo, "beta_ci_high": hi,99 "okun_negative_5pct": neg,100 "mean_growth": float(x.mean()),101 "first_quarter": str(df.index[0]), "last_quarter": str(df.index[-1]),102 "n_analysed": int(T),103 "feasibility": "complete", "feasibility_reason": None,104 "data_issues": data_issues, "n_input_rows": n_rows, "random_seed": None,105 # supplementary106 "baseline_beta_ci_low_iid": float(b[1] - z * se_iid), "baseline_beta_ci_high_iid": float(b[1] + z * se_iid),107 "t_stat_hac": t_hac, "p_value_hac": p_hac, "critical_value": float(z),108 "se_hac_by_lags": sens, "resid_lag1_autocorr": rho1,109 "mean_du": float(y.mean()),110 "input_sha256": {f: sha(os.path.join(D, f)) for f in ["UNRATE.csv", "GDPC1.csv"]},111}112with open(os.path.join(O, "results.json"), "w") as f:113 json.dump(res, f, indent=2)114115fig, ax = plt.subplots(figsize=(6.5, 4.5))116ax.scatter(x, y, s=12, alpha=0.6, label="quarters 1960Q1-2019Q4")117xs = np.linspace(x.min(), x.max(), 100)118ax.plot(xs, b[0] + b[1] * xs, "r-", label=f"OLS: Δu = {b[0]:.3f} {b[1]:+.3f} g")119ax.set_xlabel("Real GDP growth g_t (%, q/q, not annualised)")120ax.set_ylabel("Change in unemployment rate Δu_t (pp)")121ax.set_title("Okun's law, U.S. 1960Q1-2019Q4 (n=240)")122ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)123ax.legend(fontsize=8); fig.tight_layout()124fig.savefig(os.path.join(O, "figure1.png"), dpi=120)125print(f"beta={b[1]:.4f} se_iid={se_iid:.4f} se_hac(L={L})={se_hac:.4f} CI=[{lo:.4f},{hi:.4f}] neg5%={neg}")126