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 plt910np.random.seed(0)11OUT = "out"12gdp = pd.read_csv("data/GDPC1.csv", parse_dates=["observation_date"])13un = pd.read_csv("data/UNRATE.csv", parse_dates=["observation_date"])14n_rows = {"GDPC1.csv": len(gdp), "UNRATE.csv": len(un)}15issues = []1617# ---- audit18for name, df, col, freq in [("GDPC1.csv", gdp, "GDPC1", "QS"), ("UNRATE.csv", un, "UNRATE", "MS")]:19 d = df["observation_date"]20 if d.duplicated().any():21 issues.append({"issue": f"{name}: duplicated dates", "handling": "would abort"}); raise SystemExit22 full = pd.date_range(d.min(), d.max(), freq=freq)23 if len(full) != len(d):24 issues.append({"issue": f"{name}: {len(full)-len(d)} missing dates", "handling": "checked against sample window"})25 miss = df[df[col].isna()]26 for dt in miss["observation_date"]:27 inside = pd.Timestamp("1959-10-01") <= dt <= pd.Timestamp("2019-12-31")28 issues.append({"issue": f"{name}: missing value at {dt.date()}",29 "handling": ("inside sample window -- quarter dropped" if inside else30 "outside the 1959Q4-2019Q4 window used; no effect on estimates")})31 if (df[col].dropna() <= 0).any():32 issues.append({"issue": f"{name}: non-positive values", "handling": "inspected"})33issues.append({"issue": "Both series extend beyond 2019Q4 (to 2026) and are current vintages (revised data, GDP in chained 2017 dollars), not real-time",34 "handling": "observations after 2019Q4 excluded before any computation; revision/vintage effects acknowledged as limitation"})3536# ---- build quarterly data (only through 2019Q4; 1959Q4 needed for first difference)37un = un[(un.observation_date >= "1959-10-01") & (un.observation_date <= "2019-12-01")].copy()38un["q"] = un.observation_date.dt.to_period("Q")39cnt = un.groupby("q")["UNRATE"].count()40assert (cnt == 3).all(), "incomplete quarter in window"41uq = un.groupby("q")["UNRATE"].mean()42gdp = gdp[(gdp.observation_date >= "1959-10-01") & (gdp.observation_date <= "2019-10-01")].copy()43gdp["q"] = gdp.observation_date.dt.to_period("Q")44gq = gdp.set_index("q")["GDPC1"]45df = pd.DataFrame({"u": uq, "y": gq})46df["du"] = df.u.diff()47df["g"] = 100 * np.log(df.y).diff()48df = df.loc[pd.Period("1960Q1"):pd.Period("2019Q4")].dropna()49assert len(df) == 2405051def ols(d, hac=None):52 X = sm.add_constant(d["g"].values)53 m = sm.OLS(d["du"].values, X)54 return m.fit(cov_type="HAC", cov_kwds={"maxlags": hac, "kernel": "bartlett"}) if hac is not None else m.fit()5556T = len(df)57L = int(np.floor(4 * (T / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb -> 458full = ols(df, L)59s1 = df.loc[:pd.Period("1983Q4")]; s2 = df.loc[pd.Period("1984Q1"):]60r1, r2 = ols(s1, L), ols(s2, L)61# HAC lag sensitivity62sens = {lag: float(ols(df, lag).bse[1]) for lag in [0, 2, 4, 8, 12]}63print("HAC SE sensitivity by lag:", sens)64print("OLS (iid) SE:", float(ols(df).bse[1]))6566# Chow67ssr_p = float(ols(df).ssr); ssr_u = float(ols(s1).ssr + ols(s2).ssr); k = 268F = ((ssr_p - ssr_u) / k) / (ssr_u / (T - 2 * k))69p = float(stats.f.sf(F, k, T - 2 * k))7071# Robustness: HAC Wald test of break in (alpha,beta) jointly and in beta only (dummy-interaction regression)72D = (df.index >= pd.Period("1984Q1")).astype(float)73Xi = np.column_stack([np.ones(T), df.g.values, D, D * df.g.values])74mi = sm.OLS(df.du.values, Xi).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett"})75w_joint = mi.wald_test(np.array([[0, 0, 1, 0], [0, 0, 0, 1]]), use_f=False, scalar=True)76w_slope = mi.wald_test(np.array([[0, 0, 0, 1]]), use_f=False, scalar=True)77w_int = mi.wald_test(np.array([[0, 0, 1, 0]]), use_f=False, scalar=True)78rob = {"hac_se_by_lag": {str(k_): v for k_, v in sens.items()}, "ols_iid_se_beta": float(ols(df).bse[1]),79 "hac_wald_joint_break_chi2": float(w_joint.statistic), "hac_wald_joint_break_p": float(w_joint.pvalue),80 "hac_wald_slope_break_chi2": float(w_slope.statistic), "hac_wald_slope_break_p": float(w_slope.pvalue),81 "hac_wald_intercept_break_chi2": float(w_int.statistic), "hac_wald_intercept_break_p": float(w_int.pvalue),82 "slope_diff_sub2_minus_sub1": float(mi.params[3]), "slope_diff_hac_se": float(mi.bse[3]),83 "intercept_diff_sub2_minus_sub1": float(mi.params[2]), "intercept_diff_hac_se": float(mi.bse[2]),84 "breakeven_growth_sub1_pct_qq": float(-r1.params[0] / r1.params[1]),85 "breakeven_growth_sub2_pct_qq": float(-r2.params[0] / r2.params[1]),86 "mean_g_sub1": float(s1.g.mean()), "mean_g_sub2": float(s2.g.mean()),87 "durbin_watson_full_resid": float(sm.stats.stattools.durbin_watson(ols(df).resid))}88print("Robustness:", json.dumps(rob, indent=2))8990b, se = float(full.params[1]), float(full.bse[1])91z = stats.norm.ppf(0.975)92lo, hi = b - z * se, b + z * se93concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")9495res = {96 "task_id": "F6-conf-03",97 "sample_start": "1960-01-01", "sample_end": "2019-10-01",98 "n_sub1": int(len(s1)), "n_sub2": int(len(s2)),99 "baseline_metric": float(np.sqrt(np.mean((df.du - df.du.mean()) ** 2))),100 "model_metric": float(np.sqrt(np.mean(full.resid ** 2))),101 "effect_estimate": b, "intercept": float(full.params[0]),102 "effect_se": se, "hac_lags": L,103 "effect_ci_low": lo, "effect_ci_high": hi,104 "beta_sub1": float(r1.params[1]), "beta_sub2": float(r2.params[1]),105 "chow_f_stat": float(F), "chow_p_value": p,106 "stability_conclusion": "stable" if p >= 0.05 else "unstable",107 "conclusion": concl,108 "feasibility": "complete", "feasibility_reason": None,109 "data_issues": issues, "n_input_rows": n_rows,110 "n_analysed": int(T), "random_seed": 0,111}112os.makedirs(OUT, exist_ok=True)113json.dump(res, open(f"{OUT}/results.json", "w"), indent=2)114json.dump(rob, open(f"{OUT}/robustness.json", "w"), indent=2)115116pd.DataFrame([117 {"period": "full_1960Q1-2019Q4", "n": T, "intercept": full.params[0], "slope": full.params[1], "slope_se": full.bse[1]},118 {"period": "sub1_1960Q1-1983Q4", "n": len(s1), "intercept": r1.params[0], "slope": r1.params[1], "slope_se": r1.bse[1]},119 {"period": "sub2_1984Q1-2019Q4", "n": len(s2), "intercept": r2.params[0], "slope": r2.params[1], "slope_se": r2.bse[1]},120]).to_csv(f"{OUT}/table1.csv", index=False)121122fig, ax = plt.subplots(figsize=(7, 5))123ax.scatter(s1.g, s1.du, s=12, alpha=.6, label="1960Q1–1983Q4")124ax.scatter(s2.g, s2.du, s=12, alpha=.6, label="1984Q1–2019Q4")125xs = np.linspace(df.g.min(), df.g.max(), 50)126ax.plot(xs, r1.params[0] + r1.params[1] * xs, label=f"fit 1: β={r1.params[1]:.3f}")127ax.plot(xs, r2.params[0] + r2.params[1] * xs, label=f"fit 2: β={r2.params[1]:.3f}")128ax.axhline(0, c="grey", lw=.5); ax.axvline(0, c="grey", lw=.5)129ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)")130ax.set_ylabel("Δ unemployment rate Δu_t (pp)")131ax.set_title(f"Okun's law, US 1960–2019 (Chow F={F:.2f}, p={p:.3f})")132ax.legend(fontsize=8); fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=150)133print(json.dumps(res, indent=2))134