1"""F2-conf-02: OLS slope of y on x with HAC (Newey-West, Bartlett) SE.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/series.csv; writes out/results.json, out/exclusions.csv, out/hac_sensitivity.csv.4No randomness is used (random_seed = null)."""5import json, math6import numpy as np, pandas as pd7import statsmodels.api as sm8from scipy import stats910d = pd.read_csv("data/series.csv", dtype={"month": str})11n_raw = len(d)12issues = []13# --- audit / exclusions ---14excl = []15def step(name, rule, mask, why):16 global d17 nb = len(d); d = d[mask].copy(); excl.append(dict(step=name, rule=rule, n_before=nb, n_after=len(d), n_removed=nb-len(d), rationale=why))18step("1", "drop rows with missing month/x/y", d[["month","x","y"]].notna().all(axis=1), "OLS requires complete cases")19step("2", "drop duplicated months", ~d.duplicated("month"), "one observation per month")20per = pd.PeriodIndex(d["month"], freq="M")21gaps = int((np.diff(per.asi8) != 1).sum())22if gaps: issues.append({"issue": f"{gaps} gaps/out-of-order months", "handling": "kept; HAC uses row order"})23pd.DataFrame(excl).to_csv("out/exclusions.csv", index=False)24n = len(d)2526y = d["y"].to_numpy(); x = d["x"].to_numpy()27X = sm.add_constant(x)28ols = sm.OLS(y, X).fit()29alpha_hat, beta_hat = ols.params30se_iid = ols.bse[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_T = 1.1447 * (alpha1 * T)^(1/3),36# alpha1 = 4 rho^2 / ((1-rho)^2 (1+rho)^2). Rounded up to an integer.37v = (x - x.mean()) * u38rho_v = float(np.sum(v[1:]*v[:-1]) / np.sum(v**2))39alpha1 = 4*rho_v**2 / ((1-rho_v)**2 * (1+rho_v)**2)40andrews_bw = 1.1447 * (alpha1 * n) ** (1/3)41hac_lags = int(math.ceil(andrews_bw))42nw_rule = int(math.floor(4 * (n/100) ** (2/9))) # Newey-West (1994) rule of thumb, for reference4344def hac(L):45 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": True})46 return r47r = hac(hac_lags)48se_hac = float(r.bse[1])49t_hac = float(beta_hat / se_hac)50p_hac = float(2 * stats.norm.sf(abs(t_hac)))51z = stats.norm.ppf(0.975)52ci = (float(beta_hat - z*se_hac), float(beta_hat + z*se_hac))5354sens = []55for L in sorted(set([0, 1, 2, 3, nw_rule, 6, 8, hac_lags, 12, 16, 20, 24])):56 rr = hac(L); se = float(rr.bse[1]); t = beta_hat/se57 sens.append(dict(lags=L, se_hac=se, t=t, p_normal=float(2*stats.norm.sf(abs(t)))))58sens = pd.DataFrame(sens); sens.to_csv("out/hac_sensitivity.csv", index=False)5960res = {61 "task_id": "F2-conf-02",62 "beta_hat": float(beta_hat), "alpha_hat": float(alpha_hat),63 "se_iid": float(se_iid), "t_iid": float(ols.tvalues[1]), "p_iid": float(ols.pvalues[1]),64 "rho_resid": rho_resid,65 "se_hac": se_hac, "hac_lags": hac_lags, "t_hac": t_hac, "p_hac": p_hac,66 "beta_ci_low": ci[0], "beta_ci_high": ci[1],67 "memo_claim_supported": bool(p_hac < 0.05),68 "n_analysed": int(n),69 "feasibility": "complete", "feasibility_reason": None,70 "data_issues": issues,71 "n_input_rows": {"series.csv": int(n_raw)},72 "random_seed": None,73 "bandwidth_details": {"rule": "Andrews (1991) AR(1) plug-in, Bartlett kernel, on slope score (x-xbar)*u, rounded up",74 "rho_score": rho_v, "andrews_bw_raw": andrews_bw, "nw1994_rule_of_thumb_lags": nw_rule,75 "small_sample_correction": "statsmodels use_correction=True (n/(n-k))",76 "p_value_reference": "standard normal"},77 "hac_sensitivity": sens.to_dict(orient="records"),78 "sensitivity_all_p_above_005_for_lags_ge_1": bool((sens.loc[sens.lags >= 1, "p_normal"] > 0.05).all()),79}80with open("out/results.json", "w") as f:81 json.dump(res, f, indent=2)82print(json.dumps({k: res[k] for k in ["beta_hat","alpha_hat","se_iid","rho_resid","se_hac","hac_lags","t_hac","p_hac","beta_ci_low","beta_ci_high","memo_claim_supported"]}, indent=1))83print(sens.to_string(index=False))84