1"""F2-conf-04: Okun's law, US quarterly 1960Q1-2019Q4, OLS with Newey-West HAC SE.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, os, math4import numpy as np, pandas as pd5import statsmodels.api as sm6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt910OUT = "out"11FIRST, LAST = pd.Period("1960Q1", "Q"), pd.Period("2019Q4", "Q")12LEVEL_START = FIRST - 1 # 1959Q4 levels used for first differences1314un_raw = pd.read_csv("data/UNRATE.csv")15gdp_raw = pd.read_csv("data/GDPC1.csv")16n_input_rows = {"UNRATE.csv": int(len(un_raw)), "GDPC1.csv": int(len(gdp_raw))}1718data_issues = []19un = un_raw.copy()20un["UNRATE"] = pd.to_numeric(un["UNRATE"], errors="coerce")21un["q"] = pd.to_datetime(un["observation_date"]).dt.to_period("Q")22miss = un[un["UNRATE"].isna()]23for _, r in miss.iterrows():24 inwin = LEVEL_START <= r["q"] <= LAST25 data_issues.append({"issue": f"UNRATE missing for {r['observation_date']}",26 "handling": ("inside window: quarter would be dropped" if inwin else27 "outside the 1959Q4-2019Q4 window; no effect on the analysis")})28gdp = gdp_raw.copy()29gdp["GDPC1"] = pd.to_numeric(gdp["GDPC1"], errors="coerce")30gdp["q"] = pd.to_datetime(gdp["observation_date"]).dt.to_period("Q")31for _, r in gdp[gdp["GDPC1"].isna()].iterrows():32 data_issues.append({"issue": f"GDPC1 missing for {r['observation_date']}", "handling": "quarter dropped"})3334# quarterly average of the three monthly observations; require all 3 present35uq = un.groupby("q")["UNRATE"].agg(["mean", "count"])36uq = uq[uq["count"] == 3]["mean"].rename("u")37df = gdp.set_index("q")[["GDPC1"]].join(uq, how="left").sort_index()38excl = []39n0 = len(df) # n_raw = GDPC1 quarters40df = df.dropna(subset=["u", "GDPC1"])41excl.append(("0", "drop quarters lacking GDPC1 or a complete 3-month UNRATE average", n0, len(df),42 "1947Q1-Q4 precede UNRATE start; 2025Q4 lacks Oct-2025 UNRATE; all outside window"))43n0 = len(df)44df = df[(df.index >= LEVEL_START) & (df.index <= LAST)]45excl.append(("1", "keep quarters 1959Q4-2019Q4 (levels)", n0, len(df),46 "task window; 1959Q4 kept only for differencing"))47assert len(df) == 241 and (df.index[1:].astype(int) - df.index[:-1].astype(int) == 1).all()48df["du"] = df["u"].diff()49df["g"] = 100 * np.log(df["GDPC1"] / df["GDPC1"].shift(1))50n1 = len(df)51df = df.dropna(subset=["du", "g"])52excl.append(("2", "drop 1959Q4 (no prior level for differences)", n1, len(df),53 "first-difference loses one observation"))54pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "rationale"]).assign(55 n_removed=lambda d: d.n_before - d.n_after)[["step", "rule", "n_before", "n_after", "n_removed", "rationale"]56].to_csv(os.path.join(OUT, "exclusions.csv"), index=False)5758T = len(df)59X = sm.add_constant(df["g"].values)60y = df["du"].values61ols = sm.OLS(y, X).fit()62# Newey-West (1994) plug-in style rule of thumb: floor(4*(T/100)^(2/9))63hac_lags = int(math.floor(4 * (T / 100) ** (2 / 9)))64hac = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": hac_lags, "kernel": "bartlett",65 "use_correction": False})66beta, alpha = float(ols.params[1]), float(ols.params[0])67se_iid, se_hac = float(ols.bse[1]), float(hac.bse[1])68z = 1.95996398454005469lo, hi = beta - z * se_hac, beta + z * se_hac70t_hac = beta / se_hac71sens = {}72for L in [0, 1, 2, 4, 6, 8, 12]:73 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": False})74 sens[str(L)] = round(float(r.bse[1]), 10)75resid = ols.resid76rho = [float(np.corrcoef(resid[k:], resid[:-k])[0, 1]) for k in range(1, 5)]7778results = {79 "task_id": "F2-conf-04",80 "beta_hat": beta, "alpha_hat": alpha, "r_squared": float(ols.rsquared),81 "se_iid": se_iid, "se_hac": se_hac, "hac_lags": hac_lags,82 "beta_ci_low": lo, "beta_ci_high": hi,83 "okun_negative_5pct": bool(beta < 0 and abs(t_hac) > z),84 "mean_growth": float(df["g"].mean()),85 "first_quarter": str(df.index[0]), "last_quarter": str(df.index[-1]),86 "n_analysed": int(T),87 "feasibility": "complete", "feasibility_reason": None,88 "data_issues": data_issues, "n_input_rows": n_input_rows, "random_seed": None,89 # supplementary keys90 "t_stat_hac": float(t_hac),91 "p_value_hac_two_sided": float(2 * (1 - __import__("scipy.stats").stats.norm.cdf(abs(t_hac)))),92 "se_hac_by_lags": sens,93 "resid_autocorr_lag1_4": rho,94 "mean_du": float(df["du"].mean()),95}96with open(os.path.join(OUT, "results.json"), "w") as f:97 json.dump(results, f, indent=2)9899fig, ax = plt.subplots(figsize=(6.5, 4.5))100ax.scatter(df["g"], df["du"], s=12, alpha=0.6)101xs = np.linspace(df["g"].min(), df["g"].max(), 50)102ax.plot(xs, alpha + beta * xs, color="red", label=f"OLS: Δu = {alpha:.3f} {beta:+.3f}·g")103ax.axhline(0, lw=0.5, color="k"); ax.axvline(0, lw=0.5, color="k")104ax.set_xlabel("g_t: real GDP growth, % q/q (not annualised)")105ax.set_ylabel("Δu_t: change in unemployment rate, pp")106ax.set_title("Okun's law, US 1960Q1–2019Q4 (n=%d)" % T); ax.legend()107fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)108print(json.dumps({k: results[k] for k in ["beta_hat", "se_iid", "se_hac", "hac_lags", "beta_ci_low",109 "beta_ci_high", "okun_negative_5pct", "n_analysed"]}, indent=1))110