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/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")11n_raw = len(d)12issues = []13excl = []14# --- audit / exclusions15n0 = len(d); d = d.dropna(subset=["month", "x", "y"])16excl.append(("1", "drop rows with missing month/x/y", n0, len(d), n0 - len(d), "regression requires complete cases"))17n0 = len(d); d = d.drop_duplicates(subset=["month"])18excl.append(("2", "drop duplicated months", n0, len(d), n0 - len(d), "one observation per month"))19per = pd.PeriodIndex(d["month"], freq="M")20d = d.iloc[np.argsort(per.asi8)].reset_index(drop=True)21per = pd.PeriodIndex(d["month"], freq="M")22gaps = int((np.diff(per.asi8) != 1).sum())23if gaps:24 issues.append({"issue": f"{gaps} gaps in monthly sequence", "handling": "kept; HAC lags computed on row order"})25pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv("out/exclusions.csv", index=False)2627X = sm.add_constant(d["x"].to_numpy()); y = d["y"].to_numpy(); T = len(y)2829# --- baseline: conventional OLS (the memo's measure)30ols = sm.OLS(y, X).fit()31alpha_hat, beta_hat = map(float, ols.params)32se_iid = float(ols.bse[1])33u = ols.resid34rho_resid = float(np.sum(u[1:] * u[:-1]) / np.sum(u * u))3536# --- bandwidth: Andrews (1991) AR(1) plug-in for the Bartlett kernel, applied to the37# slope score v_t = (x_t - xbar) u_t (the series whose long-run variance enters the slope SE).38v = (d["x"].to_numpy() - d["x"].mean()) * u39rho_v = float(np.sum(v[1:] * v[:-1]) / np.sum(v[:-1] ** 2))40a1 = 4 * rho_v**2 / ((1 - rho_v) ** 2 * (1 + rho_v) ** 2)41andrews_bw = 1.1447 * (a1 * T) ** (1 / 3)42hac_lags = int(math.ceil(andrews_bw))43nw94_rule = int(math.floor(4 * (T / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb, for comparison4445def hac(L):46 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": False})47 return float(r.bse[1])4849se_hac = hac(hac_lags)50t_hac = beta_hat / se_hac51p_hac = float(2 * stats.norm.sf(abs(t_hac)))52z = stats.norm.ppf(0.975)53lo, hi = beta_hat - z * se_hac, beta_hat + z * se_hac5455# sensitivity over bandwidths (reported, not used for the verdict)56rows = []57for L in sorted(set([0, 1, 2, nw94_rule, 6, 8, 10, 12, hac_lags, 16, 20, 24, 30, 36])):58 s = hac(L); t = beta_hat / s59 rows.append({"lags": L, "se": s, "t": t, "p_normal": 2 * stats.norm.sf(abs(t))})60sens = pd.DataFrame(rows); sens.to_csv("out/sensitivity.csv", index=False)6162res = {63 "task_id": "F2-conf-02",64 "beta_hat": beta_hat, "alpha_hat": alpha_hat, "se_iid": se_iid,65 "t_iid": beta_hat / se_iid, "p_iid": float(2 * stats.t.sf(abs(beta_hat / se_iid), T - 2)),66 "rho_resid": rho_resid, "rho_score": rho_v,67 "se_hac": se_hac, "hac_lags": hac_lags,68 "andrews_bandwidth_raw": andrews_bw, "nw1994_rule_lags": nw94_rule,69 "se_hac_nw1994_rule": hac(nw94_rule),70 "t_hac": t_hac, "p_hac": p_hac,71 "beta_ci_low": float(lo), "beta_ci_high": float(hi),72 "memo_claim_supported": bool(p_hac < 0.05),73 "sensitivity_max_p": float(sens.loc[sens.lags >= nw94_rule, "p_normal"].max()),74 "sensitivity_min_p": float(sens.loc[sens.lags >= nw94_rule, "p_normal"].min()),75 "n_analysed": int(T),76 "feasibility": "complete", "feasibility_reason": None,77 "data_issues": issues,78 "n_input_rows": {"series.csv": int(n_raw)},79 "random_seed": None,80}81with open("out/results.json", "w") as f:82 json.dump(res, f, indent=2)83print(json.dumps(res, indent=2)); print(sens.to_string(index=False))84