1"""F2-conf-04: Okun's law, U.S. quarterly 1960Q1-2019Q4, OLS + 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, math5from pathlib import Path6import numpy as np, pandas as pd7import statsmodels.api as sm8from scipy import stats9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213ROOT = Path(__file__).resolve().parent.parent14DATA, OUT = ROOT / "data", ROOT / "out"15OUT.mkdir(exist_ok=True)16FIRST, LAST = pd.Period("1960Q1", "Q"), pd.Period("2019Q4", "Q")17LEVEL_START = FIRST - 1 # 1959Q41819un = pd.read_csv(DATA / "UNRATE.csv")20gd = pd.read_csv(DATA / "GDPC1.csv")21n_input_rows = {"UNRATE.csv": int(len(un)), "GDPC1.csv": int(len(gd))}22data_issues = []23excl = []2425# --- unemployment: monthly -> quarterly average (require all 3 months) ---26un["date"] = pd.to_datetime(un["observation_date"])27assert not un["date"].duplicated().any()28miss = un[un["UNRATE"].isna()]29for d in miss["observation_date"]:30 p = pd.Period(d, "M").asfreq("Q")31 inside = LEVEL_START <= p <= LAST32 data_issues.append({"issue": f"UNRATE missing value at {d} (empty field)",33 "handling": ("inside the needed 1959Q4-2019Q4 range; affected quarter would be dropped"34 if inside else35 "outside the 1959Q4-2019Q4 range needed for the regression; dropped, no effect on results")})36# monthly UNRATE is an input to the quarterly panel, not the unit of analysis; the exclusion37# ledger follows the quarterly panel (unit of analysis). Monthly rows: window filter then dropna.38un["q"] = un["date"].dt.to_period("Q")39un = un[(un["q"] >= LEVEL_START) & (un["q"] <= LAST)].dropna(subset=["UNRATE"])40qu = un.groupby("q")["UNRATE"].agg(["mean", "count"])41assert (qu["count"] == 3).all(), "incomplete quarter in window"42u = qu["mean"]4344# --- GDP ---45gd["q"] = pd.to_datetime(gd["observation_date"]).dt.to_period("Q")46assert not gd["q"].duplicated().any() and gd["GDPC1"].notna().all() and (gd["GDPC1"] > 0).all()47n0 = len(gd); gd = gd[(gd["q"] >= LEVEL_START) & (gd["q"] <= LAST)]48excl.append(("G1", "keep GDPC1 quarters in 1959Q4-2019Q4", n0, len(gd), "window rule: outside 1960Q1-2019Q4 only 1959Q4 levels are used"))49y = gd.set_index("q")["GDPC1"]5051n0 = len(y)52df = pd.concat({"u": u, "gdp": y}, axis=1, join="inner").sort_index()53excl.append(("M1", "inner-join with quarters having 3 non-missing UNRATE months", n0, len(df), "u_t requires complete quarter (none missing in window)"))54assert df.index.is_monotonic_increasing and len(df) == 241 and df.notna().all().all()55assert (df.index[1:] - df.index[:-1] == 1 * df.index.freq).all() if df.index.freq else True56df["du"] = df["u"].diff()57df["g"] = 100 * np.log(df["gdp"] / df["gdp"].shift(1))58n0 = len(df); reg = df.loc[FIRST:LAST].dropna()59excl.append(("Q1", "drop 1959Q4 (levels-only quarter used for first differences)", n0, len(reg), "differencing"))60assert len(reg) == 2406162pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "rationale"]).assign(63 n_removed=lambda d: d.n_before - d.n_after)[["step", "rule", "n_before", "n_after", "n_removed", "rationale"]64].to_csv(OUT / "exclusions.csv", index=False)6566# --- OLS ---67T = len(reg)68X = sm.add_constant(reg["g"].values)69Y = reg["du"].values70ols = sm.OLS(Y, X).fit()71# Newey-West (1994) plug-in rule of thumb for Bartlett kernel: floor(4 (T/100)^(2/9))72hac_lags = int(math.floor(4 * (T / 100) ** (2 / 9)))73hac = 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)77ci = (beta - z * se_hac, beta + z * se_hac)78t_hac = beta / se_hac79p_hac = 2 * stats.norm.sf(abs(t_hac))8081# sensitivity over bandwidths (not required keys; reported as extras)82sens = {}83for L in [0, 1, 2, 4, 8, 12, int(math.floor(0.75 * T ** (1 / 3)))]:84 r = sm.OLS(Y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": False})85 sens[str(L)] = float(r.bse[1])86hac_dfc = float(sm.OLS(Y, X).fit(cov_type="HAC", cov_kwds={"maxlags": hac_lags, "use_correction": True}).bse[1])87resid_ac1 = float(np.corrcoef(ols.resid[1:], ols.resid[:-1])[0, 1])8889fq = lambda p: f"{p.year}Q{p.quarter}"90res = {91 "task_id": "F2-conf-04",92 "beta_hat": beta, "alpha_hat": alpha, "r_squared": float(ols.rsquared),93 "se_iid": se_iid, "se_hac": se_hac, "hac_lags": hac_lags,94 "beta_ci_low": float(ci[0]), "beta_ci_high": float(ci[1]),95 "okun_negative_5pct": bool(beta < 0 and p_hac < 0.05),96 "mean_growth": float(reg["g"].mean()),97 "first_quarter": fq(reg.index[0]), "last_quarter": fq(reg.index[-1]),98 "n_analysed": int(T),99 "feasibility": "complete", "feasibility_reason": None,100 "data_issues": data_issues, "n_input_rows": n_input_rows, "random_seed": None,101 # extras102 "baseline_ci_iid_low": float(beta - z * se_iid), "baseline_ci_iid_high": float(beta + z * se_iid),103 "t_stat_hac": float(t_hac), "p_value_hac": float(p_hac), "z_crit": float(z),104 "se_hac_small_sample_corrected": hac_dfc,105 "se_hac_by_lags": sens, "resid_lag1_autocorr": resid_ac1,106 "mean_du": float(reg["du"].mean()),107 "hac_lag_rule": "floor(4*(T/100)^(2/9)) (Newey-West 1994 rule of thumb, Bartlett kernel)",108}109(OUT / "results.json").write_text(json.dumps(res, indent=2) + "\n")110111fig, ax = plt.subplots(figsize=(6.5, 4.5))112ax.scatter(reg["g"], reg["du"], s=12, alpha=0.7)113gg = np.linspace(reg["g"].min(), reg["g"].max(), 50)114ax.plot(gg, alpha + beta * gg, color="C3", label=f"OLS: Δu = {alpha:.3f} {beta:+.3f}·g")115ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)116ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)")117ax.set_ylabel("Change in unemployment rate Δu_t (pp)")118ax.set_title("Okun's law, U.S. 1960Q1–2019Q4 (n=240)")119ax.legend(); fig.tight_layout(); fig.savefig(OUT / "figure1.png", dpi=150)120121print(f"beta={beta:.4f} alpha={alpha:.4f} R2={ols.rsquared:.3f} se_iid={se_iid:.4f} "122 f"se_hac(L={hac_lags})={se_hac:.4f} CI=[{ci[0]:.4f},{ci[1]:.4f}] neg5%={res['okun_negative_5pct']}")123