1"""F2-conf-02: OLS slope of y on x with Newey-West HAC inference.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/series.csv; writes out/results.json. No randomness used."""4import json, math5from pathlib import Path6import numpy as np7import pandas as pd8from scipy import stats9import statsmodels.api as sm1011ROOT = Path(__file__).resolve().parent.parent12OUT = ROOT / "out"1314df = pd.read_csv(ROOT / "data" / "series.csv", dtype={"month": str})15n_raw = len(df)16issues = []17# data checks18per = pd.PeriodIndex(df["month"], freq="M")19assert per.is_monotonic_increasing and not per.duplicated().any()20assert ((per[1:] - per[:-1]).map(lambda d: d.n) == 1).all(), "gaps in months"21d = df.dropna(subset=["x", "y"])22assert len(d) == n_raw23T = len(d)2425X = sm.add_constant(d["x"].to_numpy())26y = d["y"].to_numpy()27ols = sm.OLS(y, X).fit()28alpha_hat, beta_hat = map(float, ols.params)29se_iid = float(ols.bse[1])30t_iid = float(ols.tvalues[1]); p_iid = float(ols.pvalues[1])31u = ols.resid32rho_resid = float(np.sum(u[1:] * u[:-1]) / np.sum(u ** 2))3334# Bandwidth: Andrews (1991) AR(1) plug-in for the Bartlett kernel, applied to the35# slope score v_t = (x_t - xbar) * u_t: S = 1.1447 * (a1 * T)^(1/3),36# a1 = 4 rho_v^2 / ((1-rho_v)^2 (1+rho_v)^2). Rounded up to an integer lag count.37v = (d["x"].to_numpy() - d["x"].mean()) * u38rho_v = float(np.sum(v[1:] * v[:-1]) / np.sum(v[:-1] ** 2))39a1 = 4 * rho_v ** 2 / ((1 - rho_v) ** 2 * (1 + rho_v) ** 2)40andrews_bw = 1.1447 * (a1 * T) ** (1 / 3)41nw_rule = int(math.floor(4 * (T / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb42hac_lags = int(max(math.ceil(andrews_bw), nw_rule))4344def hac(L):45 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": True})46 return float(r.bse[1])4748se_hac = hac(hac_lags)49t_hac = beta_hat / se_hac50p_hac = float(2 * stats.norm.sf(abs(t_hac)))51z = float(stats.norm.ppf(0.975))52ci_lo, ci_hi = beta_hat - z * se_hac, beta_hat + z * se_hac5354# sensitivity of the verdict to the bandwidth (reported, not used for the primary key)55sens = {}56for L in sorted(set([nw_rule, hac_lags, 8, 12, 18, 24, int(math.floor(0.75 * T ** (1 / 3)))])):57 s = hac(L); t = beta_hat / s58 sens[str(L)] = {"se_hac": s, "t_hac": t, "p_hac": float(2 * stats.norm.sf(abs(t)))}5960res = {61 "task_id": "F2-conf-02",62 "beta_hat": beta_hat,63 "alpha_hat": alpha_hat,64 "se_iid": se_iid,65 "rho_resid": rho_resid,66 "se_hac": se_hac,67 "hac_lags": hac_lags,68 "t_hac": t_hac,69 "p_hac": p_hac,70 "beta_ci_low": ci_lo,71 "beta_ci_high": ci_hi,72 "memo_claim_supported": bool(p_hac < 0.05),73 "n_analysed": T,74 "feasibility": "complete",75 "feasibility_reason": None,76 "data_issues": issues,77 "n_input_rows": {"series.csv": n_raw},78 "random_seed": None,79 # supplementary keys80 "baseline_t_iid": t_iid,81 "baseline_p_iid": p_iid,82 "rho_score": rho_v,83 "andrews_bandwidth_raw": andrews_bw,84 "nw1994_rule_lags": nw_rule,85 "se_ratio_hac_over_iid": se_hac / se_iid,86 "hac_sensitivity_by_lag": sens,87}88OUT.mkdir(exist_ok=True)89(OUT / "results.json").write_text(json.dumps(res, indent=2))90print(f"beta={beta_hat:.4f} alpha={alpha_hat:.4f} se_iid={se_iid:.4f} (t={t_iid:.2f}, p={p_iid:.4f})")91print(f"rho_resid={rho_resid:.3f} rho_score={rho_v:.3f} lags={hac_lags} se_hac={se_hac:.4f} t={t_hac:.2f} p={p_hac:.4f}")92print(f"CI=[{ci_lo:.4f},{ci_hi:.4f}] memo_claim_supported={res['memo_claim_supported']}")93for L, s in sens.items(): print(L, {k: round(v, 4) for k, v in s.items()})94