1"""F2-conf-02: OLS slope of y on x with HAC (Newey-West, Bartlett) standard errors.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, math4import numpy as np, pandas as pd5import statsmodels.api as sm6from scipy import stats78d = pd.read_csv("data/series.csv")9n_rows = len(d)10issues = []11# data checks12per = pd.PeriodIndex(d["month"], freq="M")13if d[["x", "y"]].isna().any().any():14 issues.append({"issue": "missing values", "handling": "rows dropped"})15 d = d.dropna(subset=["x", "y"])16if per.duplicated().any():17 issues.append({"issue": "duplicate months", "handling": "flagged"})18gaps = int(((per[1:] - per[:-1]).map(lambda o: o.n) != 1).sum())19if gaps:20 issues.append({"issue": f"{gaps} non-consecutive month steps", "handling": "flagged; HAC assumes equal spacing"})21d = d.sort_values("month").reset_index(drop=True)2223X = sm.add_constant(d["x"].values)24y = d["y"].values25T = len(y)26ols = sm.OLS(y, X).fit()27alpha_hat, beta_hat = ols.params28se_iid = ols.bse[1]29u = ols.resid30rho_resid = float(np.sum(u[1:] * u[:-1]) / np.sum(u ** 2))3132# Bandwidth: Andrews (1991) AR(1) plug-in for the Bartlett kernel, applied to the33# slope score v_t = (x_t - xbar) u_t. Data-dependent, adapts to the strong persistence.34v = (d["x"].values - d["x"].values.mean()) * u35rv = float(np.sum(v[1:] * v[:-1]) / np.sum(v[:-1] ** 2))36a1 = 4 * rv ** 2 / ((1 - rv) ** 2 * (1 + rv) ** 2)37andrews_bw = 1.1447 * (a1 * T) ** (1 / 3)38hac_lags = int(math.ceil(andrews_bw)) # Bartlett weights 1 - j/(L+1), j<=L39nw_rule = int(math.floor(4 * (T / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb4041def hac(L):42 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": True})43 return float(r.bse[1])4445se_hac = hac(hac_lags)46t_hac = beta_hat / se_hac47p_hac = 2 * stats.norm.sf(abs(t_hac))48z = stats.norm.ppf(0.975)49ci = (beta_hat - z * se_hac, beta_hat + z * se_hac)5051sens = {}52for L in sorted({0, nw_rule, 8, 12, hac_lags, 18, 24, int(math.floor(0.75 * T ** (1 / 3)))}):53 s = hac(L)54 sens[L] = {"se": s, "t": beta_hat / s, "p_normal": 2 * stats.norm.sf(abs(beta_hat / s)),55 "p_t": 2 * stats.t.sf(abs(beta_hat / s), T - 2)}5657res = {58 "task_id": "F2-conf-02",59 "beta_hat": float(beta_hat), "alpha_hat": float(alpha_hat),60 "se_iid": float(se_iid), "rho_resid": rho_resid,61 "se_hac": se_hac, "hac_lags": hac_lags,62 "t_hac": float(t_hac), "p_hac": float(p_hac),63 "beta_ci_low": float(ci[0]), "beta_ci_high": float(ci[1]),64 "memo_claim_supported": bool(p_hac < 0.05),65 "n_analysed": int(T), "feasibility": "complete", "feasibility_reason": None,66 "data_issues": issues, "n_input_rows": {"series.csv": int(n_rows)},67 "random_seed": None,68}69with open("out/results.json", "w") as f:70 json.dump(res, f, indent=2)71diag = {"t_iid": float(beta_hat / se_iid), "p_iid": float(ols.pvalues[1]),72 "rho_score": rv, "andrews_bw": andrews_bw, "nw_rule_lags": nw_rule,73 "sensitivity": {str(k): v for k, v in sens.items()}}74with open("out/diagnostics.json", "w") as f:75 json.dump(diag, f, indent=2)76print(json.dumps(res, indent=2)); print(json.dumps(diag, indent=2))77