1"""Okun's law, US 1960Q1-2019Q4 (task F6-conf-03). 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 np7import pandas as pd8import statsmodels.api as sm9from scipy import stats10import matplotlib11matplotlib.use("Agg")12import matplotlib.pyplot as plt1314OUT = "out"15os.makedirs(OUT, exist_ok=True)16S0, S1 = pd.Period("1960Q1", "Q"), pd.Period("2019Q4", "Q")17SPLIT = pd.Period("1984Q1", "Q")1819gdp_raw = pd.read_csv("data/GDPC1.csv")20un_raw = pd.read_csv("data/UNRATE.csv")21n_input_rows = {"GDPC1.csv": int(len(gdp_raw)), "UNRATE.csv": int(len(un_raw))}22data_issues = []2324# --- unemployment: quarterly average of monthly values (require all 3 months) ---25un = un_raw.copy()26un["date"] = pd.to_datetime(un["observation_date"])27un["UNRATE"] = pd.to_numeric(un["UNRATE"], errors="coerce")28nan_months = un.loc[un["UNRATE"].isna(), "observation_date"].tolist()29if nan_months:30 data_issues.append({"issue": f"UNRATE missing value(s) for month(s) {nan_months}",31 "handling": "Quarter containing a missing month is not averaged (set missing); "32 "all such months lie outside the 1959Q4-2019Q4 window, so no effect on estimates."})33un["q"] = un["date"].dt.to_period("Q")34uq = un.groupby("q")["UNRATE"].agg(["mean", "count"])35uq.loc[uq["count"] < 3, "mean"] = np.nan # incomplete quarters36u = uq["mean"].rename("u")3738gdp = gdp_raw.copy()39gdp["q"] = pd.to_datetime(gdp["observation_date"]).dt.to_period("Q")40assert not gdp["q"].duplicated().any()41g_level = gdp.set_index("q")["GDPC1"].astype(float)4243# --- merge and exclusion ledger (unit: quarters of GDPC1) ---44ledger = []45def log(step, rule, before, after, why):46 ledger.append(dict(step=step, rule=rule, n_before=before, n_after=after,47 n_removed=before - after, rationale=why))4849df = pd.DataFrame({"gdp": g_level})50n0 = len(df)51df = df.join(u, how="left")52n1 = int(df["u"].notna().sum())53df_ok = df[df["u"].notna()].copy()54log(1, "drop GDPC1 quarters without a complete 3-month UNRATE quarterly average", n0, n1,55 "UNRATE starts 1948-01 (1947 quarters have no unemployment data); 2025Q4 has a missing month")56# differences require consecutive quarters57full = df.reindex(pd.period_range(df.index.min(), df.index.max(), freq="Q"))58full["du"] = full["u"].diff()59full["g"] = 100 * np.log(full["gdp"]).diff()60df_d = full.loc[df_ok.index].dropna(subset=["du", "g"])61log(2, "drop quarters whose first difference cannot be formed (previous quarter missing)", n1, len(df_d),62 "first available quarter (1948Q1) and the quarter after the incomplete 2025Q4")63d = df_d.loc[(df_d.index >= S0) & (df_d.index <= S1)].copy()64log(3, "restrict to regression window 1960Q1-2019Q4", len(df_d), len(d),65 "TASK.md information rule: quarters outside 1960Q1-2019Q4 must not enter any estimate")66pd.DataFrame(ledger).to_csv(f"{OUT}/exclusions.csv", index=False)6768expected = pd.period_range(S0, S1, freq="Q")69assert len(d) == 240 and (d.index == expected).all(), "sample is not the 240 consecutive quarters"70# 1959Q4 levels used for 1960Q1 differences: check present71assert pd.notna(full.loc[pd.Period("1959Q4", "Q"), "u"]) and pd.notna(full.loc[pd.Period("1959Q4", "Q"), "gdp"])7273# --- baseline ---74du, g = d["du"].values, d["g"].values75n = len(d)76baseline_rmse = float(np.sqrt(np.mean((du - du.mean()) ** 2)))7778# --- full-sample OLS with Newey-West HAC ---79HAC_LAGS = int(np.floor(4 * (n / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb -> 4 for n=24080X = sm.add_constant(g)81ols = sm.OLS(du, X).fit()82hac = sm.OLS(du, X).fit(cov_type="HAC", cov_kwds={"maxlags": HAC_LAGS, "use_correction": True})83alpha, beta = float(ols.params[0]), float(ols.params[1])84se_b = float(hac.bse[1])85z = stats.norm.ppf(0.975)86ci = (beta - z * se_b, beta + z * se_b)87model_rmse = float(np.sqrt(np.mean(ols.resid ** 2)))88# lag-length sensitivity (reported, not used for the headline)89hac_sens = {}90for L in [0, 2, 4, 8, 12]:91 r = sm.OLS(du, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": True})92 hac_sens[str(L)] = float(r.bse[1])93ols_se_classical = float(ols.bse[1])9495# --- subperiods and Chow test ---96m1 = d.index < SPLIT97def fit(mask):98 y, x = du[mask], sm.add_constant(g[mask])99 r = sm.OLS(y, x).fit()100 rh = sm.OLS(y, x).fit(cov_type="HAC", cov_kwds={"maxlags": HAC_LAGS, "use_correction": True})101 return r, rh102r1, r1h = fit(m1)103r2, r2h = fit(~m1)104ssr_p = float(ols.ssr); ssr_u = float(r1.ssr + r2.ssr)105k = 2106df2 = n - 2 * k107F = ((ssr_p - ssr_u) / k) / (ssr_u / df2)108p = float(stats.f.sf(F, k, df2))109110# robustness (supplementary): interacted model, HAC Wald tests of equality111D = (~m1).astype(float)112Xi = np.column_stack([np.ones(n), g, D, D * g])113ri = sm.OLS(du, Xi).fit(cov_type="HAC", cov_kwds={"maxlags": HAC_LAGS, "use_correction": True})114w_joint = ri.wald_test(np.array([[0, 0, 1, 0], [0, 0, 0, 1]]), use_f=False, scalar=True)115w_slope = ri.wald_test(np.array([[0, 0, 0, 1]]), use_f=False, scalar=True)116w_int = ri.wald_test(np.array([[0, 0, 1, 0]]), use_f=False, scalar=True)117resid_var_sub1 = float(r1.ssr / (m1.sum() - 2)); resid_var_sub2 = float(r2.ssr / ((~m1).sum() - 2))118119if ci[0] > 0: concl = "positive_effect"120elif ci[1] < 0: concl = "negative_effect"121else: concl = "no_detectable_effect"122123data_issues.append({"issue": "GDPC1 is the current (2026) vintage and extends to 2026Q3; UNRATE extends to 2026-09",124 "handling": "Observations after 2019Q4 dropped before estimation (exclusion step 3); "125 "revised (not real-time) GDP is used, consistent with an ex-post descriptive estimand."})126127res = {128 "task_id": "F6-conf-03",129 "feasibility": "complete",130 "feasibility_reason": None,131 "data_issues": data_issues,132 "n_input_rows": n_input_rows,133 "n_analysed": int(n),134 "random_seed": None,135 "sample_start": str(d.index[0].start_time.date()),136 "sample_end": str(d.index[-1].start_time.date()),137 "n_sub1": int(m1.sum()),138 "n_sub2": int((~m1).sum()),139 "baseline_metric": baseline_rmse,140 "model_metric": model_rmse,141 "effect_estimate": beta,142 "intercept": alpha,143 "effect_se": se_b,144 "hac_lags": HAC_LAGS,145 "effect_ci_low": float(ci[0]),146 "effect_ci_high": float(ci[1]),147 "beta_sub1": float(r1.params[1]),148 "beta_sub2": float(r2.params[1]),149 "chow_f_stat": float(F),150 "chow_p_value": p,151 "stability_conclusion": "stable" if p >= 0.05 else "unstable",152 "conclusion": concl,153 # supplementary (not required by TASK.md)154 "intercept_sub1": float(r1.params[0]),155 "intercept_sub2": float(r2.params[0]),156 "beta_sub1_hac_se": float(r1h.bse[1]),157 "beta_sub2_hac_se": float(r2h.bse[1]),158 "effect_se_ols_classical": ols_se_classical,159 "effect_se_hac_by_lag": hac_sens,160 "r_squared": float(ols.rsquared),161 "chow_df": [k, int(df2)],162 "ssr_pooled": ssr_p,163 "ssr_unrestricted": ssr_u,164 "hac_wald_joint_equality_p": float(w_joint.pvalue),165 "hac_wald_slope_equality_p": float(w_slope.pvalue),166 "hac_wald_intercept_equality_p": float(w_int.pvalue),167 "resid_var_sub1": resid_var_sub1,168 "resid_var_sub2": resid_var_sub2,169 "growth_for_stable_u_sub1": float(-r1.params[0] / r1.params[1]),170 "growth_for_stable_u_sub2": float(-r2.params[0] / r2.params[1]),171}172with open(f"{OUT}/results.json", "w") as f:173 json.dump(res, f, indent=2)174175pd.DataFrame([176 dict(period="full 1960Q1-2019Q4", n=n, intercept=alpha, slope=beta, slope_se=se_b),177 dict(period="sub1 1960Q1-1983Q4", n=int(m1.sum()), intercept=float(r1.params[0]), slope=float(r1.params[1]), slope_se=float(r1h.bse[1])),178 dict(period="sub2 1984Q1-2019Q4", n=int((~m1).sum()), intercept=float(r2.params[0]), slope=float(r2.params[1]), slope_se=float(r2h.bse[1])),179]).to_csv(f"{OUT}/table1.csv", index=False)180181fig, ax = plt.subplots(figsize=(7, 5))182ax.scatter(g[m1], du[m1], s=14, alpha=0.6, label="1960Q1-1983Q4")183ax.scatter(g[~m1], du[~m1], s=14, alpha=0.6, marker="^", label="1984Q1-2019Q4")184for r, mask, lab in [(r1, m1, "fit 1960-1983"), (r2, ~m1, "fit 1984-2019")]:185 xs = np.linspace(g[mask].min(), g[mask].max(), 50)186 ax.plot(xs, r.params[0] + r.params[1] * xs,187 label=f"{lab}: du = {r.params[0]:.3f} {r.params[1]:+.3f} g")188ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="grey", lw=0.5)189ax.set_xlabel("Real GDP growth g_t (% q/q, not annualised)")190ax.set_ylabel("Change in unemployment rate du_t (pp)")191ax.set_title("Okun's law, US 1960Q1-2019Q4")192ax.legend(fontsize=8)193fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120)194195print(f"n={n} beta={beta:.4f} HAC se={se_b:.4f} (L={HAC_LAGS}) CI=[{ci[0]:.4f},{ci[1]:.4f}] "196 f"beta1={r1.params[1]:.4f} beta2={r2.params[1]:.4f} Chow F={F:.3f} p={p:.4g}")197