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.py3Reads only data/; writes out/results.json, out/figure1.png, out/exclusions.csv. No randomness."""4import json, os, hashlib5import numpy as np, pandas as pd6import statsmodels.api as sm7from statsmodels.tsa.stattools import acf8from scipy import stats9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213ROOT = os.getcwd(); OUT = os.path.join(ROOT, "out"); os.makedirs(OUT, exist_ok=True)14FIRST, LAST = pd.Period("1960Q1", "Q"), pd.Period("2019Q4", "Q")15EXPECTED_SHA = {"UNRATE.csv": "621ff6fadf50910cb7a21ee20ba610c247cb043b29b4e05c4c1a72fadca2a06c",16 "GDPC1.csv": "1c6c6e9f59e383e20c849dd299175c7a307526c09abdc18f18da9754fc5b3343"}17issues, excl = [], []1819def sha(p):20 return hashlib.sha256(open(p, "rb").read()).hexdigest()2122for f, h in EXPECTED_SHA.items():23 if sha(os.path.join("data", f)) != h:24 issues.append({"issue": f"{f} SHA-256 differs from data/README.md", "handling": "used as given; flagged"})2526u = pd.read_csv("data/UNRATE.csv", dtype={"UNRATE": str}, keep_default_na=False)27g = pd.read_csv("data/GDPC1.csv", dtype={"GDPC1": str}, keep_default_na=False)28n_rows = {"UNRATE.csv": int(len(u)), "GDPC1.csv": int(len(g))}29u["date"] = pd.to_datetime(u.observation_date); g["date"] = pd.to_datetime(g.observation_date)30u["UNRATE"] = pd.to_numeric(u.UNRATE.replace("", np.nan)); g["GDPC1"] = pd.to_numeric(g.GDPC1.replace("", np.nan))31assert u.date.is_unique and g.date.is_unique3233# Step 1: drop months with missing UNRATE (monthly rows)34miss = u[u.UNRATE.isna()]35if len(miss):36 issues.append({"issue": "UNRATE missing (empty field) for month(s) " + ", ".join(miss.observation_date) +37 "; outside the 1959Q4-2019Q4 window used",38 "handling": "row dropped; affects no quarter used in the analysis"})39n0 = len(u); u = u.dropna(subset=["UNRATE"])40# (monthly row drop is documented in data_issues; the exclusion ledger below is in quarterly units41# so that it chains from the quarterly panel to n_analysed)4243# Quarterly aggregation44u["q"] = u.date.dt.to_period("Q")45uq = u.groupby("q").UNRATE.agg(["mean", "count"])46g["q"] = g.date.dt.to_period("Q"); gq = g.set_index("q").GDPC147need = pd.period_range(FIRST - 1, LAST, freq="Q")48assert (uq["count"].reindex(need) == 3).all(), "incomplete quarter of monthly UNRATE in window"49assert gq.reindex(need).notna().all(), "missing GDPC1 in window"5051# Step 2: build quarterly panel over the union of quarters, then restrict to the window52allq = pd.period_range(min(uq.index.min(), gq.index.min()), max(uq.index.max(), gq.index.max()), freq="Q")53full = pd.DataFrame({"u": uq["mean"], "Y": gq}).reindex(allq)54full["du"] = full.u.diff(); full["g"] = 100 * np.log(full.Y / full.Y.shift(1))55sample = full.loc[FIRST:LAST]56excl.append(("1", "restrict to regression window 1960Q1-2019Q4 (1959Q4 levels used only for differences)",57 len(full), len(sample), "unit: quarters (union of both series 1947Q1-2026Q3); window fixed by TASK.md; quarters 1947Q1-1959Q4 and 2020Q1-2026Q3 out of window"))58n_b = len(sample); sample = sample.dropna(subset=["du", "g"])59excl.append(("2", "drop window quarters with missing du or g", n_b, len(sample), "none expected"))60pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "rationale"]).assign(61 n_removed=lambda d: d.n_before - d.n_after)[["step", "rule", "n_before", "n_after", "n_removed", "rationale"]62 ].to_csv(os.path.join(OUT, "exclusions.csv"), index=False)6364n = len(sample); assert n == 24065X = sm.add_constant(sample.g.values); y = sample.du.values66ols = sm.OLS(y, X).fit()67# Bandwidth: Newey-West (1994) rule of thumb for the Bartlett kernel, floor(4*(T/100)^(2/9))68L = int(np.floor(4 * (n / 100) ** (2 / 9)))69def hac_se(lag):70 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": lag, "kernel": "bartlett", "use_correction": False})71 return float(r.bse[1])72b, a = float(ols.params[1]), float(ols.params[0])73se_iid, se_hac = float(ols.bse[1]), hac_se(L)74z = float(stats.norm.ppf(0.975))75lo, hi = b - z * se_hac, b + z * se_hac76p_hac = float(2 * stats.norm.sf(abs(b / se_hac)))77sens = {str(l): hac_se(l) for l in [0, 1, 2, 4, 6, 8, 12]}78res_acf = acf(ols.resid, nlags=4, fft=False)[1:]7980results = {81 "task_id": "F2-conf-04", "beta_hat": b, "alpha_hat": a, "r_squared": float(ols.rsquared),82 "se_iid": se_iid, "se_hac": se_hac, "hac_lags": L,83 "beta_ci_low": float(lo), "beta_ci_high": float(hi),84 "okun_negative_5pct": bool(b < 0 and p_hac < 0.05),85 "mean_growth": float(sample.g.mean()),86 "first_quarter": str(sample.index[0]), "last_quarter": str(sample.index[-1]),87 "n_analysed": int(n), "feasibility": "complete", "feasibility_reason": None,88 "data_issues": issues, "n_input_rows": n_rows, "random_seed": None,89 # supplementary keys (not required by TASK.md)90 "t_stat_hac": float(b / se_hac), "p_value_hac": p_hac, "z_crit_95": z,91 "baseline_ci_iid_low": float(b - z * se_iid), "baseline_ci_iid_high": float(b + z * se_iid),92 "se_hac_by_lag": sens, "resid_acf_lags1to4": [float(x) for x in res_acf],93 "mean_du": float(sample.du.mean()),94}95with open(os.path.join(OUT, "results.json"), "w") as fh:96 json.dump(results, fh, indent=2)9798fig, ax = plt.subplots(figsize=(6, 4.5))99ax.scatter(sample.g, sample.du, s=12, alpha=0.7)100xs = np.linspace(sample.g.min(), sample.g.max(), 50)101ax.plot(xs, a + b * xs, "r-", label=f"OLS: Δu = {a:.3f} + ({b:.3f}) g")102ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)103ax.set_xlabel("g_t: real GDP growth, % q/q (not annualised)"); ax.set_ylabel("Δu_t: change in unemployment rate, pp")104ax.set_title("Okun's law, US 1960Q1–2019Q4 (n=240)"); ax.legend()105fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)106print(json.dumps({k: results[k] for k in ["beta_hat", "se_iid", "se_hac", "hac_lags", "beta_ci_low", "beta_ci_high", "okun_negative_5pct"]}))107