1"""F2-conf-04: Okun's law, US quarterly 1960Q1-2019Q4, OLS + Newey-West HAC SE.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, os4import numpy as np, pandas as pd5import statsmodels.api as sm6from scipy import stats7import matplotlib8matplotlib.use("Agg")9import matplotlib.pyplot as plt1011OUT = "out"12un = pd.read_csv("data/UNRATE.csv")13gd = pd.read_csv("data/GDPC1.csv")14n_rows = {"UNRATE.csv": int(len(un)), "GDPC1.csv": int(len(gd))}15issues = []1617un["date"] = pd.to_datetime(un["observation_date"])18gd["date"] = pd.to_datetime(gd["observation_date"])19un["UNRATE"] = pd.to_numeric(un["UNRATE"], errors="coerce")20gd["GDPC1"] = pd.to_numeric(gd["GDPC1"], errors="coerce")2122miss_u = un.loc[un["UNRATE"].isna(), "observation_date"].tolist()23if miss_u:24 issues.append({"issue": f"UNRATE missing value(s) at {miss_u} (empty field in the FRED file)",25 "handling": "Outside the 1959Q4-2019Q4 window needed for estimation; no effect. "26 "Quarterly averages require all 3 months present, so affected quarters would be NaN."})27miss_g = gd.loc[gd["GDPC1"].isna(), "observation_date"].tolist()28if miss_g:29 issues.append({"issue": f"GDPC1 missing value(s) at {miss_g}", "handling": "Quarter set to NaN."})3031un["q"] = un["date"].dt.to_period("Q")32uq = un.groupby("q")["UNRATE"].agg(["mean", "count"])33uq.loc[uq["count"] < 3, "mean"] = np.nan # require all 3 months34gd["q"] = gd["date"].dt.to_period("Q")35assert not gd["q"].duplicated().any()36df = pd.DataFrame({"u": uq["mean"]}).join(gd.set_index("q")["GDPC1"], how="outer").sort_index()37# ensure consecutive quarters38full = pd.period_range(df.index.min(), df.index.max(), freq="Q")39df = df.reindex(full)40df["du"] = df["u"].diff()41df["g"] = 100 * np.log(df["GDPC1"]).diff()42s = df.loc[pd.Period("1960Q1"):pd.Period("2019Q4")]43assert len(s) == 240 and s[["du", "g"]].notna().all().all()4445y = s["du"].values46X = sm.add_constant(s["g"].values)47T = len(y)48ols = sm.OLS(y, X).fit()49# Newey-West (1994) rule-of-thumb bandwidth for the Bartlett kernel: floor(4*(T/100)^(2/9))50L = int(np.floor(4 * (T / 100) ** (2 / 9)))51hac = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": True})52beta, alpha = float(ols.params[1]), float(ols.params[0])53se_iid = float(ols.bse[1]); se_hac = float(hac.bse[1])54z = stats.norm.ppf(0.975)55ci = (beta - z * se_hac, beta + z * se_hac)56p_hac = float(2 * stats.norm.sf(abs(beta / se_hac)))5758# sensitivity of HAC SE to bandwidth (reported in report, not a required key)59sens = {}60for lag in [0, 1, 2, 4, 6, 8, 12, 16]:61 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": lag, "kernel": "bartlett", "use_correction": True})62 sens[lag] = float(r.bse[1])63acf_res = [float(np.corrcoef(ols.resid[k:], ols.resid[:-k])[0, 1]) for k in range(1, 5)]6465res = {66 "task_id": "F2-conf-04",67 "beta_hat": beta, "alpha_hat": alpha, "r_squared": float(ols.rsquared),68 "se_iid": se_iid, "se_hac": se_hac, "hac_lags": L,69 "beta_ci_low": float(ci[0]), "beta_ci_high": float(ci[1]),70 "okun_negative_5pct": bool(beta < 0 and p_hac < 0.05),71 "mean_growth": float(s["g"].mean()),72 "first_quarter": str(s.index[0]), "last_quarter": str(s.index[-1]),73 "n_analysed": int(T),74 "feasibility": "complete", "feasibility_reason": None,75 "data_issues": issues, "n_input_rows": n_rows, "random_seed": None,76}77with open(os.path.join(OUT, "results.json"), "w") as f:78 json.dump(res, f, indent=2)79with open(os.path.join(OUT, "diagnostics.json"), "w") as f:80 json.dump({"p_value_hac_normal": p_hac, "t_hac": beta / se_hac,81 "hac_se_by_lag": sens, "resid_acf_1to4": acf_res,82 "dw": float(sm.stats.durbin_watson(ols.resid))}, f, indent=2)8384fig, ax = plt.subplots(figsize=(6.5, 4.5))85ax.scatter(s["g"], s["du"], s=12, alpha=0.6)86gx = np.linspace(s["g"].min(), s["g"].max(), 50)87ax.plot(gx, alpha + beta * gx, "r-", label=f"OLS: du = {alpha:.3f} {beta:+.3f} g")88ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)")89ax.set_ylabel("Change in unemployment rate du_t (pp)")90ax.set_title("Okun's law, US 1960Q1-2019Q4 (n=240)")91ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)92ax.legend(); fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=130)93print(json.dumps(res, indent=2))94