1"""F2-conf-02: OLS slope of y on x with Newey-West HAC SE. Reads data/ only; writes out/results.json."""2import json, os, numpy as np, pandas as pd3import statsmodels.api as sm4from scipy import stats56ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))7d = pd.read_csv(os.path.join(ROOT, "data", "series.csv"))8n_raw = len(d)9issues = []10# audit checks (no rows removed unless a problem is found)11per = pd.PeriodIndex(d["month"], freq="M")12assert d[["x", "y"]].notna().all().all(), "missing values"13assert not d["month"].duplicated().any(), "duplicate months"14assert (np.diff(per.asi8) == 1).all(), "gaps in months"15d = d.sort_values("month").reset_index(drop=True)16n = len(d)1718X = sm.add_constant(d["x"].to_numpy())19y = d["y"].to_numpy()20ols = sm.OLS(y, X).fit()21alpha_hat, beta_hat = map(float, ols.params)22se_iid = float(ols.bse[1])23u = ols.resid24rho_resid = float(np.sum(u[1:] * u[:-1]) / np.sum(u ** 2))2526# Bandwidth: Andrews (1991) AR(1) plug-in for the Bartlett kernel, applied to the27# score of the slope v_t = (x_t - xbar) u_t. S = 1.1447 (a1 * T)^(1/3),28# a1 = 4 r^2 / ((1-r)^2 (1+r)^2), r = AR(1) coefficient of v_t. Truncated to integer lags.29v = (d["x"].to_numpy() - d["x"].mean()) * u30r_v = float(np.sum(v[1:] * v[:-1]) / np.sum(v[:-1] ** 2))31a1 = 4 * r_v ** 2 / ((1 - r_v) ** 2 * (1 + r_v) ** 2)32S_andrews = 1.1447 * (a1 * n) ** (1 / 3)33hac_lags = int(np.floor(S_andrews))34nw_rule = int(np.floor(4 * (n / 100) ** (2 / 9))) # Newey-West (1994) rule-of-thumb, for sensitivity3536def hac(L):37 f = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": True})38 se = float(f.bse[1]); t = beta_hat / se39 return se, t, float(2 * stats.norm.sf(abs(t)))4041se_hac, t_hac, p_hac = hac(hac_lags)42z = stats.norm.ppf(0.975)43ci_low, ci_high = beta_hat - z * se_hac, beta_hat + z * se_hac44sens = {str(L): dict(zip(["se", "t", "p"], hac(L))) for L in sorted({nw_rule, 8, 12, 16, 24, hac_lags})}45memo_supported = bool(p_hac < 0.05)4647res = {48 "task_id": "F2-conf-02",49 "beta_hat": beta_hat, "alpha_hat": alpha_hat, "se_iid": se_iid,50 "t_iid": float(ols.tvalues[1]), "p_iid": float(ols.pvalues[1]),51 "rho_resid": rho_resid, "rho_score": r_v, "andrews_bandwidth": float(S_andrews),52 "se_hac": se_hac, "hac_lags": hac_lags, "t_hac": t_hac, "p_hac": p_hac,53 "beta_ci_low": float(ci_low), "beta_ci_high": float(ci_high),54 "se_ratio_hac_to_iid": se_hac / se_iid,55 "hac_sensitivity": sens, "nw_rule_of_thumb_lags": nw_rule,56 "sensitivity_all_insignificant": bool(all(s["p"] >= 0.05 for s in sens.values())),57 "memo_claim_supported": memo_supported,58 "n_analysed": int(n), "feasibility": "complete", "feasibility_reason": None,59 "data_issues": issues, "n_input_rows": {"series.csv": int(n_raw)}, "random_seed": None,60}61os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)62with open(os.path.join(ROOT, "out", "results.json"), "w") as fh:63 json.dump(res, fh, indent=2)64pd.DataFrame([{"step": "S0", "rule": "none (no missing, duplicate or gap rows)", "n_before": n_raw,65 "n_after": n, "n_removed": n_raw - n, "rationale": "data audit found no problems"}]66 ).to_csv(os.path.join(ROOT, "out", "exclusions.csv"), index=False)67print(f"beta={beta_hat:.4f} se_iid={se_iid:.4f} se_hac={se_hac:.4f} (L={hac_lags}) t={t_hac:.2f} p={p_hac:.4f} memo_supported={memo_supported}")68