1"""F2-conf-04: Okun's law, U.S. quarterly 1960Q1-2019Q4, OLS with Newey-West HAC SE.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/; writes out/results.json, out/exclusions.csv, out/figure1.png. No randomness."""4import json, math, os5import numpy as np6import pandas as pd7import statsmodels.api as sm8from scipy import stats9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213OUT = "out"14os.makedirs(OUT, exist_ok=True)15FIRST, LAST = pd.Period("1960Q1", "Q"), pd.Period("2019Q4", "Q")16LEVEL_START = FIRST - 1 # 1959Q4 levels needed for first differences1718un = pd.read_csv("data/UNRATE.csv", keep_default_na=False, dtype=str)19gd = pd.read_csv("data/GDPC1.csv", keep_default_na=False, dtype=str)20n_input_rows = {"UNRATE.csv": int(len(un)), "GDPC1.csv": int(len(gd))}21data_issues = []22excl = []2324def add_excl(step, rule, nb, na, why):25 excl.append(dict(step=step, rule=rule, n_before=nb, n_after=na, n_removed=nb - na, rationale=why))2627# --- unemployment: monthly -> quarterly average28un["date"] = pd.to_datetime(un["observation_date"])29un["v"] = pd.to_numeric(un["UNRATE"].replace("", np.nan), errors="coerce")30miss = un[un["v"].isna()]31if len(miss):32 inwin = [d for d in miss["date"] if LEVEL_START <= d.to_period("Q") <= LAST]33 data_issues.append({34 "issue": f"UNRATE missing (empty) for {', '.join(miss['observation_date'])}",35 "handling": ("outside the 1959Q4-2019Q4 window; does not affect the analysis" if not inwin36 else "INSIDE window: affected quarter averaged over available months")})37assert not un["date"].duplicated().any() and not gd["observation_date"].duplicated().any()38un["q"] = un["date"].dt.to_period("Q")39m = un[(un["q"] >= LEVEL_START) & (un["q"] <= LAST)]40assert m["v"].notna().all()41cnt = m.groupby("q")["v"].count()42assert (cnt == 3).all(), "incomplete quarter in window"43u = m.groupby("q")["v"].mean()4445# --- GDP46gd["q"] = pd.to_datetime(gd["observation_date"]).dt.to_period("Q")47gd["v"] = pd.to_numeric(gd["GDPC1"], errors="coerce")48assert gd["v"].notna().all()49y = gd.set_index("q")["v"]50y = y[(y.index >= LEVEL_START) & (y.index <= LAST)]5152df = pd.DataFrame({"u": u, "y": y})53assert df.notna().all().all() and len(df) == 24154df["du"] = df["u"].diff()55df["g"] = 100 * np.log(df["y"] / df["y"].shift(1))5657# exclusion ledger (in units of quarters, starting from the full GDP quarterly range)58nq_all = len(pd.period_range(gd["q"].min(), gd["q"].max(), freq="Q"))59add_excl(1, "keep quarters 1959Q4-2019Q4 (levels window)", nq_all, len(df),60 "task window 1960Q1-2019Q4 plus 1959Q4 levels for first differences; window not extended")61d = df.dropna().loc[FIRST:LAST]62add_excl(2, "drop 1959Q4 (no first difference; levels only)", len(df), len(d),63 "1959Q4 used only as lag level for 1960Q1 differences")64n = len(d)65assert n == 2406667# --- OLS: baseline iid SE and Newey-West HAC SE68X = sm.add_constant(d["g"].values)69Y = d["du"].values70ols = sm.OLS(Y, X).fit()71T = n72hac_lags = int(math.floor(4 * (T / 100) ** (2 / 9))) # Newey-West (1994) plug-in rule -> 473hac = sm.OLS(Y, X).fit(cov_type="HAC", cov_kwds={"maxlags": hac_lags, "kernel": "bartlett", "use_correction": False})74beta, alpha = float(ols.params[1]), float(ols.params[0])75se_iid, se_hac = float(ols.bse[1]), float(hac.bse[1])76z = stats.norm.ppf(0.975)77lo, hi = beta - z * se_hac, beta + z * se_hac78p_hac = float(2 * stats.norm.sf(abs(beta / se_hac)))7980# sensitivity of HAC SE to bandwidth (reported, not used for headline)81sens = {}82for L in [0, 1, 2, 4, 6, 8, 12]:83 r = sm.OLS(Y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": False})84 sens[str(L)] = float(r.bse[1])85hac_dfc = float(sm.OLS(Y, X).fit(cov_type="HAC", cov_kwds={"maxlags": hac_lags, "use_correction": True}).bse[1])86resid_ac1 = float(pd.Series(ols.resid).autocorr(1))8788res = {89 "task_id": "F2-conf-04",90 "beta_hat": beta, "alpha_hat": alpha, "r_squared": float(ols.rsquared),91 "se_iid": se_iid, "se_hac": se_hac, "hac_lags": hac_lags,92 "beta_ci_low": float(lo), "beta_ci_high": float(hi),93 "okun_negative_5pct": bool(beta < 0 and p_hac < 0.05),94 "mean_growth": float(d["g"].mean()),95 "first_quarter": str(d.index[0]), "last_quarter": str(d.index[-1]),96 "n_analysed": int(n),97 "feasibility": "complete", "feasibility_reason": None,98 "data_issues": data_issues, "n_input_rows": n_input_rows, "random_seed": None,99 # supplementary keys100 "t_hac": float(beta / se_hac), "p_value_hac": p_hac,101 "se_hac_by_lags": sens, "se_hac_lag4_dfcorrected": hac_dfc,102 "residual_autocorr_lag1": resid_ac1,103 "mean_du": float(d["du"].mean()),104 "baseline_ci_iid_low": float(beta - z * se_iid), "baseline_ci_iid_high": float(beta + z * se_iid),105}106with open(os.path.join(OUT, "results.json"), "w") as f:107 json.dump(res, f, indent=2)108pd.DataFrame(excl).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)109110fig, ax = plt.subplots(figsize=(6.5, 4.5))111ax.scatter(d["g"], d["du"], s=12, alpha=0.6)112gx = np.linspace(d["g"].min(), d["g"].max(), 100)113ax.plot(gx, alpha + beta * gx, color="red", label=f"OLS: Δu = {alpha:.3f} {beta:+.3f}·g")114ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)115ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)")116ax.set_ylabel("Δ unemployment rate Δu_t (pp)")117ax.set_title("Okun's law, U.S. 1960Q1–2019Q4 (n=240)")118ax.legend(); fig.tight_layout()119fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)120print(json.dumps({k: res[k] for k in ["beta_hat", "alpha_hat", "se_iid", "se_hac", "hac_lags",121 "beta_ci_low", "beta_ci_high", "okun_negative_5pct", "n_analysed"]}, indent=1))122