1"""Okun's law (difference form), US 1960Q1-2019Q4. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads only data/, writes out/results.json, out/table1.csv, out/figure1.png, out/exclusions.csv.4No randomness is used (random_seed = null)."""5import json, os6import numpy as np, pandas as pd7import statsmodels.api as sm8from scipy import stats9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213OUT = "out"14os.makedirs(OUT, exist_ok=True)15gdp = pd.read_csv("data/GDPC1.csv", na_values=["."])16un = pd.read_csv("data/UNRATE.csv", na_values=["."])17n_input_rows = {"GDPC1.csv": int(len(gdp)), "UNRATE.csv": int(len(un))}18gdp["observation_date"] = pd.to_datetime(gdp["observation_date"])19un["observation_date"] = pd.to_datetime(un["observation_date"])2021data_issues = []22miss = un[un["UNRATE"].isna()]["observation_date"].dt.strftime("%Y-%m-%d").tolist()23if miss:24 data_issues.append({"issue": f"UNRATE missing value(s) at {miss} (blank in raw FRED CSV)",25 "handling": "outside the 1959Q4-2019Q4 window; dropped with all post-2019Q4 observations, no effect on estimates"})26data_issues.append({"issue": "Both series extend beyond 2019Q4 (to 2026) and are the 2026 vintage (revised, not real-time)",27 "handling": "all observations after 2019Q4 excluded before computation; vintage noted as limitation"})2829# Quarterly unemployment = mean of 3 months; require all 3 months present30un["q"] = un["observation_date"].dt.to_period("Q")31uq = un.groupby("q")["UNRATE"].agg(["mean", "count"])32gdp["q"] = gdp["observation_date"].dt.to_period("Q")33df = gdp.set_index("q")[["GDPC1"]].join(uq, how="left")3435excl = []36def log(step, rule, before, after, why):37 excl.append(dict(step=step, rule=rule, n_before=before, n_after=after, n_removed=before - after, rationale=why))3839n0 = len(df)40win = df.loc[pd.Period("1959Q4"):pd.Period("2019Q4")].copy()41log(1, "keep quarters 1959Q4-2019Q4 (GDPC1 quarters as unit)", n0, len(win),42 "task window; 1959Q4 needed only for lagged levels; post-2019Q4 forbidden")43assert (win["count"] == 3).all() and win["GDPC1"].notna().all(), "incomplete data inside window"44win["u"] = win["mean"]45win["du"] = win["u"].diff()46win["g"] = 100 * np.log(win["GDPC1"]).diff()47reg = win.dropna(subset=["du", "g"]).copy()48log(2, "drop 1959Q4 (no lagged level, differences undefined)", len(win), len(reg),49 "first-difference construction")50pd.DataFrame(excl).to_csv(f"{OUT}/exclusions.csv", index=False)51assert len(reg) == 240 and str(reg.index[0]) == "1960Q1" and str(reg.index[-1]) == "2019Q4"5253y, X = reg["du"].values, sm.add_constant(reg["g"].values)54n = len(y)55# Baseline: mean-only prediction56baseline_rmse = float(np.sqrt(np.mean((y - y.mean()) ** 2)))57# Full sample OLS with Newey-West HAC (Bartlett) SE58HAC_LAGS = int(np.floor(4 * (n / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb -> 459ols = sm.OLS(y, X).fit()60hac = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": HAC_LAGS, "use_correction": True})61alpha, beta = float(ols.params[0]), float(ols.params[1])62se = float(hac.bse[1])63z = stats.norm.ppf(0.975)64ci = (beta - z * se, beta + z * se)65model_rmse = float(np.sqrt(np.mean(ols.resid ** 2)))6667# Subperiods68s1 = reg.index <= pd.Period("1983Q4")69s2 = ~s170fits = {}71for name, m in [("sub1", s1), ("sub2", s2)]:72 Xs = sm.add_constant(reg["g"].values[m])73 f = sm.OLS(y[m], Xs).fit()74 fh = sm.OLS(y[m], Xs).fit(cov_type="HAC", cov_kwds={"maxlags": HAC_LAGS, "use_correction": True})75 fits[name] = (f, fh, int(m.sum()))76ssr_p = ols.ssr77ssr_u = fits["sub1"][0].ssr + fits["sub2"][0].ssr78k = 279chow_f = float(((ssr_p - ssr_u) / k) / (ssr_u / (n - 2 * k)))80chow_p = float(stats.f.sf(chow_f, k, n - 2 * k))8182# Supplementary: HAC test of slope-only difference (interaction model) and lag sensitivity83d2 = s2.astype(float)84Xi = np.column_stack([np.ones(n), reg["g"].values, d2, d2 * reg["g"].values])85inter = sm.OLS(y, Xi).fit(cov_type="HAC", cov_kwds={"maxlags": HAC_LAGS, "use_correction": True})86slope_diff_p = float(inter.pvalues[3]); int_diff_p = float(inter.pvalues[2])87R = np.array([[0, 0, 1, 0], [0, 0, 0, 1]], dtype=float)88wald = inter.wald_test(R, use_f=True, scalar=True)89hac_joint_F, hac_joint_p = float(wald.statistic), float(wald.pvalue)90lag_sens = {str(L): float(sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": True}).bse[1])91 for L in (0, 2, 8, 12)}92resid_var = {"sub1": float(fits["sub1"][0].ssr / (fits["sub1"][2] - 2)), "sub2": float(fits["sub2"][0].ssr / (fits["sub2"][2] - 2))}9394conclusion = "positive_effect" if ci[0] > 0 else ("negative_effect" if ci[1] < 0 else "no_detectable_effect")95results = {96 "task_id": "F6-conf-03",97 "feasibility": "complete",98 "feasibility_reason": None,99 "data_issues": data_issues,100 "n_input_rows": n_input_rows,101 "n_analysed": int(n),102 "random_seed": None,103 "sample_start": reg.index[0].start_time.strftime("%Y-%m-%d"),104 "sample_end": reg.index[-1].start_time.strftime("%Y-%m-%d"),105 "n_sub1": fits["sub1"][2],106 "n_sub2": fits["sub2"][2],107 "baseline_metric": baseline_rmse,108 "model_metric": model_rmse,109 "effect_estimate": beta,110 "intercept": alpha,111 "effect_se": se,112 "hac_lags": HAC_LAGS,113 "effect_ci_low": float(ci[0]),114 "effect_ci_high": float(ci[1]),115 "beta_sub1": float(fits["sub1"][0].params[1]),116 "beta_sub2": float(fits["sub2"][0].params[1]),117 "chow_f_stat": chow_f,118 "chow_p_value": chow_p,119 "stability_conclusion": "stable" if chow_p >= 0.05 else "unstable",120 "conclusion": conclusion,121 # supplementary (not required)122 "ols_classical_se_full": float(ols.bse[1]),123 "r_squared_full": float(ols.rsquared),124 "intercept_sub1": float(fits["sub1"][0].params[0]),125 "intercept_sub2": float(fits["sub2"][0].params[0]),126 "beta_sub1_hac_se": float(fits["sub1"][1].bse[1]),127 "beta_sub2_hac_se": float(fits["sub2"][1].bse[1]),128 "implied_zero_du_growth_full": float(-alpha / beta),129 "hac_slope_diff_p_value": slope_diff_p,130 "hac_intercept_diff_p_value": int_diff_p,131 "hac_joint_break_f_stat": hac_joint_F,132 "hac_joint_break_p_value": hac_joint_p,133 "hac_se_lag_sensitivity": lag_sens,134 "resid_variance_sub1": resid_var["sub1"],135 "resid_variance_sub2": resid_var["sub2"],136}137with open(f"{OUT}/results.json", "w") as fh_:138 json.dump(results, fh_, indent=2)139140rows = [("full_1960Q1_2019Q4", n, alpha, beta, se)]141for lab, key in [("sub1_1960Q1_1983Q4", "sub1"), ("sub2_1984Q1_2019Q4", "sub2")]:142 f, fh, nn = fits[key]143 rows.append((lab, nn, float(f.params[0]), float(f.params[1]), float(fh.bse[1])))144pd.DataFrame(rows, columns=["period", "n", "intercept", "slope", "slope_se"]).to_csv(f"{OUT}/table1.csv", index=False)145146fig, ax = plt.subplots(figsize=(7, 5))147g = reg["g"].values148ax.scatter(g[s1], y[s1], s=14, alpha=0.6, label="1960Q1-1983Q4", color="tab:blue")149ax.scatter(g[s2], y[s2], s=14, alpha=0.6, label="1984Q1-2019Q4", color="tab:orange")150xx = np.linspace(g.min(), g.max(), 50)151for key, c in [("sub1", "tab:blue"), ("sub2", "tab:orange")]:152 a, b = fits[key][0].params153 ax.plot(xx, a + b * xx, color=c, lw=2, label=f"fit {key}: du = {a:.2f} {b:+.3f} g")154ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)155ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)")156ax.set_ylabel("Change in unemployment rate du_t (pp)")157ax.set_title("Okun's law, US 1960Q1-2019Q4")158ax.legend(fontsize=8)159fig.tight_layout()160fig.savefig(f"{OUT}/figure1.png", dpi=120)161print(json.dumps({k: results[k] for k in ["effect_estimate", "effect_se", "effect_ci_low", "effect_ci_high",162 "beta_sub1", "beta_sub2", "chow_f_stat", "chow_p_value", "baseline_metric", "model_metric"]}, indent=1))163