1"""Okun's law, US 1960Q1-2019Q4, OLS with Newey-West HAC SE. Run from workspace root."""2import json, os3import numpy as np, pandas as pd4import statsmodels.api as sm5import matplotlib; matplotlib.use("Agg")6import matplotlib.pyplot as plt7from scipy import stats89D, O = "data", "out"10os.makedirs(O, exist_ok=True)11FIRST, LAST = pd.Period("1960Q1", "Q"), pd.Period("2019Q4", "Q")1213un = pd.read_csv(f"{D}/UNRATE.csv")14gd = pd.read_csv(f"{D}/GDPC1.csv")15n_rows = {"UNRATE.csv": int(len(un)), "GDPC1.csv": int(len(gd))}16issues = []17excl = []1819# --- audit ---20for name, df, col in [("UNRATE", un, "UNRATE"), ("GDPC1", gd, "GDPC1")]:21 df["date"] = pd.to_datetime(df["observation_date"])22 df[col] = pd.to_numeric(df[col], errors="coerce")23 assert not df["date"].duplicated().any(), f"duplicate dates in {name}"24 miss = df.loc[df[col].isna(), "observation_date"].tolist()25 if miss:26 inside = [m for m in miss if FIRST - 1 <= pd.Period(m, "M").asfreq("Q") <= LAST]27 issues.append({"issue": f"{name}: missing (empty) value(s) at {miss}",28 "handling": "outside the 1959Q4-2019Q4 levels window; ignored, no effect on estimates" if not inside29 else "INSIDE window - affected quarters cannot be formed"})3031# --- quarterly unemployment: mean of 3 months, require all 3 non-missing ---32un["q"] = un["date"].dt.to_period("Q")33uq = un.groupby("q")["UNRATE"].agg(["mean", "count"]) # count = non-missing months34gd["q"] = gd["date"].dt.to_period("Q")35gq = gd.set_index("q")["GDPC1"]3637qs = pd.period_range(FIRST - 1, LAST, freq="Q") # 1959Q4..2019Q4 levels38lev = pd.DataFrame(index=qs)39lev["u"] = uq["mean"].reindex(qs)40lev["nmon"] = uq["count"].reindex(qs)41lev["gdp"] = gq.reindex(qs)42assert (lev["nmon"] == 3).all(), "incomplete months in window"43assert lev["gdp"].notna().all() and (lev["gdp"] > 0).all()44lev["du"] = lev["u"].diff()45lev["g"] = 100 * np.log(lev["gdp"] / lev["gdp"].shift(1))4647q_all = pd.period_range(min(uq.index.min(), gq.index.min()), max(uq.index.max(), gq.index.max()), freq="Q")48excl.append(("1", "restrict to 1959Q4-2019Q4 (calendar quarters spanned by either file)", len(q_all), len(qs),49 len(q_all) - len(qs), "task window plus one lag quarter; window not extended"))50dat = lev.loc[FIRST:LAST, ["du", "g"]]51excl.append(("2", "drop 1959Q4 (lag-only level quarter)", len(qs), len(dat), len(qs) - len(dat),52 "1959Q4 only supplies levels for 1960Q1 differences"))53assert dat.notna().all().all() and len(dat) == 24054pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(f"{O}/exclusions.csv", index=False)5556y, x = dat["du"].values, dat["g"].values57X = sm.add_constant(x)58T = len(y)5960# --- baseline: OLS with iid SE ---61ols = sm.OLS(y, X).fit()62# --- HAC: Newey-West (1994) rule-of-thumb bandwidth floor(4*(T/100)^(2/9)) ---63L = int(np.floor(4 * (T / 100) ** (2 / 9)))64hac = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": False})65beta, alpha = float(ols.params[1]), float(ols.params[0])66se_hac = float(hac.bse[1])67z = stats.norm.ppf(0.975)68lo, hi = beta - z * se_hac, beta + z * se_hac69p_hac = float(2 * stats.norm.sf(abs(beta / se_hac)))7071# sensitivity over bandwidths (reported, not used for decision)72sens = {}73for l in [0, 2, 4, 6, 8, 12]:74 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": l, "kernel": "bartlett", "use_correction": False})75 sens[str(l)] = float(r.bse[1])76resid_ac1 = float(np.corrcoef(ols.resid[1:], ols.resid[:-1])[0, 1])7778res = {79 "task_id": "F2-conf-04",80 "beta_hat": beta, "alpha_hat": alpha, "r_squared": float(ols.rsquared),81 "se_iid": float(ols.bse[1]), "se_hac": se_hac, "hac_lags": L,82 "beta_ci_low": float(lo), "beta_ci_high": float(hi),83 "okun_negative_5pct": bool(beta < 0 and p_hac < 0.05),84 "mean_growth": float(np.mean(x)),85 "first_quarter": str(dat.index[0]), "last_quarter": str(dat.index[-1]),86 "n_analysed": int(T),87 "feasibility": "complete", "feasibility_reason": None,88 "data_issues": issues, "n_input_rows": n_rows, "random_seed": None,89 "baseline_ci_iid_low": float(beta - z * ols.bse[1]), "baseline_ci_iid_high": float(beta + z * ols.bse[1]),90 "t_hac": float(beta / se_hac), "p_value_hac": p_hac,91 "resid_lag1_autocorr": resid_ac1,92 "se_hac_sensitivity_by_lags": sens,93 "mean_du": float(np.mean(y)),94}95with open(f"{O}/results.json", "w") as f:96 json.dump(res, f, indent=2)9798fig, ax = plt.subplots(figsize=(6.5, 4.5))99ax.scatter(x, y, s=12, alpha=0.6)100xx = np.linspace(x.min(), x.max(), 100)101ax.plot(xx, alpha + beta * xx, "r-", label=f"OLS: du = {alpha:.3f} + ({beta:.3f}) g")102ax.axhline(0, c="grey", lw=0.5); ax.axvline(0, c="grey", lw=0.5)103ax.set_xlabel("g_t: real GDP growth, % q/q (not annualised)")104ax.set_ylabel("change in quarterly unemployment rate, pp")105ax.set_title("Okun's law, US 1960Q1-2019Q4 (n=240)"); ax.legend()106fig.tight_layout(); fig.savefig(f"{O}/figure1.png", dpi=120)107print(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"]}))108