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/. Deterministic (no randomness; random_seed = null)."""4import json, os5import numpy as np, pandas as pd6from scipy import stats7import statsmodels.api as sm8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt1112OUT = "out"13os.makedirs(OUT, exist_ok=True)14FIRST_Q, LAST_Q = pd.Period("1960Q1", "Q"), pd.Period("2019Q4", "Q")15LEVEL_START = FIRST_Q - 1 # 1959Q4 levels needed for first differences1617un_raw = pd.read_csv("data/UNRATE.csv", keep_default_na=False, dtype=str)18gdp_raw = pd.read_csv("data/GDPC1.csv", keep_default_na=False, dtype=str)19n_input_rows = {"UNRATE.csv": int(len(un_raw)), "GDPC1.csv": int(len(gdp_raw))}2021data_issues, excl = [], []22def log_excl(step, rule, nb, na, why):23 excl.append(dict(step=step, rule=rule, n_before=nb, n_after=na, n_removed=nb - na, rationale=why))2425# ---- unemployment: monthly -> quarterly average26un = un_raw.copy()27un["date"] = pd.to_datetime(un["observation_date"])28un["v"] = pd.to_numeric(un["UNRATE"].replace("", np.nan), errors="coerce")29miss = un[un["v"].isna()]30for d in miss["date"]:31 q = d.to_period("Q")32 inwin = LEVEL_START <= q <= LAST_Q33 data_issues.append({"issue": f"UNRATE missing (empty field) for {d:%Y-%m}",34 "handling": ("inside required window -> quarter dropped" if inwin else35 "outside 1959Q4-2019Q4 level window; not used, no effect on estimates")})36un["q"] = un["date"].dt.to_period("Q")37assert not un["date"].duplicated().any()38g = un.groupby("q")["v"]39uq = pd.DataFrame({"u": g.mean(), "n_months": g.count()})40# ---- GDP41gd = gdp_raw.copy()42gd["q"] = pd.to_datetime(gd["observation_date"]).dt.to_period("Q")43gd["Y"] = pd.to_numeric(gd["GDPC1"].replace("", np.nan), errors="coerce")44assert not gd["q"].duplicated().any() and gd["Y"].notna().all()45gq = gd.set_index("q")[["Y"]]4647# ---- quarterly panel restricted to level window 1959Q4-2019Q448panel = uq.join(gq, how="outer").sort_index()49n_q_all = len(panel)50lev = panel.loc[LEVEL_START:LAST_Q].copy()51log_excl("1", "keep quarters 1959Q4-2019Q4 (levels window; information rule)", n_q_all, len(lev),52 "task fixes regression window 1960Q1-2019Q4; 1959Q4 levels needed only for first differences")53assert (lev["n_months"] == 3).all() and lev["Y"].notna().all(), "incomplete quarter in window"54assert len(pd.period_range(LEVEL_START, LAST_Q, freq="Q")) == len(lev)55lev["du"] = lev["u"].diff()56lev["g"] = 100 * np.log(lev["Y"] / lev["Y"].shift(1))57reg = lev.loc[FIRST_Q:LAST_Q].dropna(subset=["du", "g"])58log_excl("2", "drop 1959Q4 (no prior level for differencing; levels-only quarter)", len(lev), len(reg),59 "1959Q4 is used only as the base level for 1960Q1 differences")60n = len(reg)61assert n == 2406263# ---- OLS64y, x = reg["du"].values, reg["g"].values65X = sm.add_constant(x)66ols = sm.OLS(y, X).fit()67alpha_hat, beta_hat = map(float, ols.params)68se_iid = float(ols.bse[1])6970# Newey-West bandwidth: classic rule floor(4*(T/100)^(2/9)) (Newey & West 1994);71# cross-check with floor(0.75*T^(1/3)) (Stock & Watson textbook rule). Both = 4 for T=240.72L_nw94 = int(np.floor(4 * (n / 100) ** (2 / 9)))73L_sw = int(np.floor(0.75 * n ** (1 / 3)))74hac_lags = L_nw947576def nw_se(L, correction=False):77 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett",78 "use_correction": correction})79 return float(r.bse[1])8081# manual Bartlett-kernel check (no small-sample df correction)82def nw_manual(L):83 e = ols.resid84 u = X * e[:, None]85 S = u.T @ u86 for j in range(1, L + 1):87 w = 1 - j / (L + 1)88 G = u[j:].T @ u[:-j]89 S += w * (G + G.T)90 XtXi = np.linalg.inv(X.T @ X)91 return float(np.sqrt((XtXi @ S @ XtXi)[1, 1]))9293se_hac = nw_se(hac_lags)94assert abs(se_hac - nw_manual(hac_lags)) < 1e-1095z = stats.norm.ppf(0.975)96ci_lo, ci_hi = beta_hat - z * se_hac, beta_hat + z * se_hac97t_hac = beta_hat / se_hac98p_hac = float(2 * stats.norm.sf(abs(t_hac)))99okun_neg = bool(beta_hat < 0 and p_hac < 0.05)100101# ---- sensitivity table (written to CSV, not results.json)102resid_ac = [float(pd.Series(ols.resid).autocorr(k)) for k in range(1, 5)]103sens = pd.DataFrame([{"lags": L, "se_hac": nw_se(L), "se_hac_dfcorr": nw_se(L, True)} for L in range(0, 13)])104sens["ci_low"] = beta_hat - z * sens["se_hac"]105sens["ci_high"] = beta_hat + z * sens["se_hac"]106sens["p_value"] = 2 * stats.norm.sf(np.abs(beta_hat / sens["se_hac"]))107sens.to_csv(f"{OUT}/hac_sensitivity.csv", index=False)108pd.DataFrame(excl).to_csv(f"{OUT}/exclusions.csv", index=False)109110# ---- figure111fig, ax = plt.subplots(figsize=(6.5, 4.5))112ax.scatter(x, y, s=12, alpha=0.6, label="quarters 1960Q1–2019Q4")113xx = np.linspace(x.min(), x.max(), 100)114ax.plot(xx, alpha_hat + beta_hat * xx, color="C3",115 label=f"OLS: Δu = {alpha_hat:.3f} {beta_hat:+.3f}·g")116ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)117ax.set_xlabel("real GDP growth g_t (% q/q, not annualised)")118ax.set_ylabel("change in unemployment rate Δu_t (pp)")119ax.set_title("Okun's law, U.S. 1960Q1–2019Q4 (n=%d)" % n)120ax.legend(fontsize=8)121fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=150); plt.close(fig)122123results = {124 "task_id": "F2-conf-04",125 "beta_hat": beta_hat, "alpha_hat": alpha_hat, "r_squared": float(ols.rsquared),126 "se_iid": se_iid, "se_hac": se_hac, "hac_lags": hac_lags,127 "beta_ci_low": float(ci_lo), "beta_ci_high": float(ci_hi),128 "okun_negative_5pct": okun_neg,129 "mean_growth": float(np.mean(x)),130 "first_quarter": str(reg.index[0]), "last_quarter": str(reg.index[-1]),131 "n_analysed": int(n),132 "feasibility": "complete", "feasibility_reason": None,133 "data_issues": data_issues, "n_input_rows": n_input_rows, "random_seed": None,134 # supplementary keys135 "t_stat_hac": float(t_hac), "p_value_hac": p_hac,136 "baseline_ci_iid_low": float(beta_hat - z * se_iid), "baseline_ci_iid_high": float(beta_hat + z * se_iid),137 "hac_lags_rule_nw1994": L_nw94, "hac_lags_rule_sw": L_sw,138 "resid_autocorr_lag1_4": resid_ac,139 "se_hac_max_over_lags_0_12": float(sens["se_hac"].max()),140 "p_value_hac_max_over_lags_0_12": float(sens["p_value"].max()),141}142with open(f"{OUT}/results.json", "w") as f:143 json.dump(results, f, indent=2)144print(json.dumps({k: results[k] for k in ["beta_hat", "alpha_hat", "r_squared", "se_iid", "se_hac",145 "hac_lags", "beta_ci_low", "beta_ci_high", "okun_negative_5pct",146 "mean_growth", "n_analysed"]}, indent=1))147