1"""Okun's law (difference form), US 1960Q1-2019Q4. Run from workspace root."""2import json, os3import numpy as np, pandas as pd4import statsmodels.api as sm5from scipy import stats6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt910OUT = "out"; os.makedirs(OUT, exist_ok=True)11issues = []12gdp = pd.read_csv("data/GDPC1.csv", parse_dates=["observation_date"])13ur = pd.read_csv("data/UNRATE.csv", parse_dates=["observation_date"])14n_rows = {"GDPC1.csv": int(len(gdp)), "UNRATE.csv": int(len(ur))}1516# --- audit ---17for name, d, col in [("GDPC1.csv", gdp, "GDPC1"), ("UNRATE.csv", ur, "UNRATE")]:18 assert not d.observation_date.duplicated().any(), name + " duplicate dates"19 miss = d.loc[d[col].isna(), "observation_date"].dt.strftime("%Y-%m-%d").tolist()20 if miss:21 inside = [m for m in miss if "1959-10-01" <= m <= "2019-12-31"]22 issues.append({"issue": f"{name}: missing value(s) at {miss}",23 "handling": "outside 1959Q4-2019Q4 estimation window; no effect on estimates" if not inside24 else "inside window: quarter dropped"})25 if (d[col].dropna() <= 0).any():26 issues.append({"issue": f"{name}: non-positive values", "handling": "flagged"})27# monthly continuity28months = pd.date_range(ur.observation_date.min(), ur.observation_date.max(), freq="MS")29assert len(months) == len(ur), "UNRATE month gaps"30qs = pd.date_range(gdp.observation_date.min(), gdp.observation_date.max(), freq="QS")31assert len(qs) == len(gdp), "GDPC1 quarter gaps"3233# --- construct quarterly series ---34ur["q"] = ur.observation_date.dt.to_period("Q")35uq = ur.groupby("q")["UNRATE"].agg(["mean", "count"])36uq = uq[uq["count"] == 3]["mean"].rename("u")37g = gdp.set_index(gdp.observation_date.dt.to_period("Q"))["GDPC1"]38df = pd.concat([uq, g], axis=1).sort_index()39df["du"] = df["u"].diff()40df["g"] = 100 * np.log(df["GDPC1"]).diff()41df = df.loc[pd.Period("1960Q1"):pd.Period("2019Q4")]42assert len(df) == 240 and df[["du", "g"]].notna().all().all()43issues.append({"issue": "Observations after 2019Q4 present in both files (incl. 2020 COVID outliers); GDPC1 is the 2026 vintage (revised data, chained 2017 $)",44 "handling": "post-2019Q4 data excluded per task; current vintage used as provided (log growth rates invariant to base year)"})45# outlier screen on in-sample data46z = (df[["du", "g"]] - df[["du", "g"]].mean()) / df[["du", "g"]].std()47big = df.index[(z.abs() > 4).any(axis=1)].astype(str).tolist()48if big:49 issues.append({"issue": f"Large (|z|>4) in-sample observations: {big}",50 "handling": "retained (genuine recession/recovery quarters); influence discussed in report"})5152def ols(d, hac=None):53 X = sm.add_constant(d["g"])54 m = sm.OLS(d["du"], X)55 return m.fit(cov_type="HAC", cov_kwds={"maxlags": hac, "kernel": "bartlett"}) if hac is not None else m.fit()5657T = len(df)58L = int(np.floor(4 * (T / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb59full = ols(df, L); full_ols = ols(df)60s1 = df.loc[:pd.Period("1983Q4")]; s2 = df.loc[pd.Period("1984Q1"):]61r1, r2 = ols(s1, L), ols(s2, L)62ssr_p = full_ols.ssr; ssr_u = ols(s1).ssr + ols(s2).ssr; k = 263F = ((ssr_p - ssr_u) / k) / (ssr_u / (T - 2 * k))64p = float(stats.f.sf(F, k, T - 2 * k))6566beta = float(full.params["g"]); se = float(full.bse["g"])67z975 = stats.norm.ppf(0.975)68lo, hi = beta - z975 * se, beta + z975 * se69concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")7071# sensitivity (reported in report, not in results keys)72sens = {f"lags_{l}": float(ols(df, l).bse["g"]) for l in (0, 2, 4, 8, 12)}73sens["ols_iid_se"] = float(full_ols.bse["g"])74sens["beta_excl_z4"] = float(ols(df.drop(pd.PeriodIndex(big, freq="Q"))).params["g"]) if big else beta75# HAC-robust Wald test of break (interaction model)76d2 = df.copy(); d2["D"] = (d2.index >= pd.Period("1984Q1")).astype(float); d2["Dg"] = d2["D"] * d2["g"]77mw = sm.OLS(d2["du"], sm.add_constant(d2[["g", "D", "Dg"]])).fit(cov_type="HAC", cov_kwds={"maxlags": L})78wt = mw.wald_test("D = 0, Dg = 0", use_f=True, scalar=True)79sens["hac_wald_break_F"] = float(wt.statistic); sens["hac_wald_break_p"] = float(wt.pvalue)80sens["hac_slope_diff"] = float(mw.params["Dg"]); sens["hac_slope_diff_p"] = float(mw.pvalues["Dg"])81sens["hac_intercept_diff"] = float(mw.params["D"]); sens["hac_intercept_diff_p"] = float(mw.pvalues["D"])82sens["breakeven_sub1"] = float(-r1.params["const"] / r1.params["g"]); sens["breakeven_sub2"] = float(-r2.params["const"] / r2.params["g"])83sens["sub1_hac_se"] = float(r1.bse["g"]); sens["sub2_hac_se"] = float(r2.bse["g"])84sens["sub1_alpha"] = float(r1.params["const"]); sens["sub2_alpha"] = float(r2.params["const"])85sens["r2_full"] = float(full_ols.rsquared)86sens["breakeven_growth_full"] = float(-full.params["const"] / beta)87json.dump(sens, open(f"{OUT}/sensitivity.json", "w"), indent=2)8889res = {90 "task_id": "F6-conf-03",91 "sample_start": df.index[0].start_time.strftime("%Y-%m-%d"),92 "sample_end": df.index[-1].start_time.strftime("%Y-%m-%d"),93 "n_sub1": int(len(s1)), "n_sub2": int(len(s2)),94 "baseline_metric": float(np.sqrt(np.mean((df.du - df.du.mean()) ** 2))),95 "model_metric": float(np.sqrt(np.mean(full_ols.resid ** 2))),96 "effect_estimate": beta, "intercept": float(full.params["const"]),97 "effect_se": se, "hac_lags": L, "effect_ci_low": float(lo), "effect_ci_high": float(hi),98 "beta_sub1": float(r1.params["g"]), "beta_sub2": float(r2.params["g"]),99 "chow_f_stat": float(F), "chow_p_value": p,100 "stability_conclusion": "stable" if p >= 0.05 else "unstable",101 "conclusion": concl, "feasibility": "complete", "feasibility_reason": None,102 "data_issues": issues, "n_input_rows": n_rows, "n_analysed": int(T), "random_seed": None,103}104json.dump(res, open(f"{OUT}/results.json", "w"), indent=2)105106pd.DataFrame([107 {"period": "full_1960Q1-2019Q4", "n": T, "intercept": full.params["const"], "slope": beta, "slope_se": se},108 {"period": "sub1_1960Q1-1983Q4", "n": len(s1), "intercept": r1.params["const"], "slope": r1.params["g"], "slope_se": r1.bse["g"]},109 {"period": "sub2_1984Q1-2019Q4", "n": len(s2), "intercept": r2.params["const"], "slope": r2.params["g"], "slope_se": r2.bse["g"]},110]).to_csv(f"{OUT}/table1.csv", index=False)111112fig, ax = plt.subplots(figsize=(7, 5))113for d, r, c, lab in [(s1, r1, "tab:blue", "1960Q1–1983Q4"), (s2, r2, "tab:red", "1984Q1–2019Q4")]:114 ax.scatter(d.g, d.du, s=12, alpha=0.6, color=c, label=f"{lab} (n={len(d)})")115 xs = np.linspace(d.g.min(), d.g.max(), 50)116 ax.plot(xs, r.params["const"] + r.params["g"] * xs, color=c,117 label=f"fit {lab}: Δu = {r.params['const']:.2f} {r.params['g']:+.3f}·g")118ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)119ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)")120ax.set_ylabel("Δ unemployment rate Δu_t (pp)")121ax.set_title("Okun's law, US 1960Q1–2019Q4"); ax.legend(fontsize=8)122fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=150)123print(json.dumps(res, indent=2)); print(json.dumps(sens, indent=2))124