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; deterministic (no randomness)."""4import json, math5import numpy as np, pandas as pd6import statsmodels.api as sm7from scipy import stats89d = pd.read_csv("data/series.csv")10n_raw = len(d)11issues = []12# audit checks13assert d[["month", "x", "y"]].notna().all().all(), "missing values"14per = pd.PeriodIndex(d["month"], freq="M")15assert per.is_monotonic_increasing and per.is_unique16gaps = int(((per[1:] - per[:-1]).map(lambda z: z.n) != 1).sum())17assert gaps == 01819X = sm.add_constant(d["x"].to_numpy())20y = d["y"].to_numpy()21T = len(y)22ols = sm.OLS(y, X).fit()23alpha_hat, beta_hat = map(float, ols.params)24se_iid = float(ols.bse[1])25u = ols.resid26rho = float(np.sum(u[1:] * u[:-1]) / np.sum(u ** 2))2728# Bandwidth: Andrews (1991) AR(1) plug-in for the Bartlett kernel, applied to29# the slope score v_t = (x_t - xbar) u_t. Chosen over the fixed NW rule-of-thumb30# floor(4(T/100)^(2/9)) because errors (rho~0.7) and x are both persistent,31# which makes the score persistent and the rule of thumb too short.32v = (d["x"].to_numpy() - d["x"].mean()) * u33rv = float(np.sum(v[1:] * v[:-1]) / np.sum(v ** 2))34a1 = 4 * rv ** 2 / ((1 - rv) ** 2 * (1 + rv) ** 2)35S_andrews = 1.1447 * (a1 * T) ** (1 / 3)36hac_lags = int(math.ceil(S_andrews))37nw_rule = int(math.floor(4 * (T / 100) ** (2 / 9)))3839def hac(L):40 f = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": True})41 return float(f.bse[1])4243se_hac = hac(hac_lags)44t_hac = beta_hat / se_hac45p_hac = float(2 * stats.norm.sf(abs(t_hac)))46z = float(stats.norm.ppf(0.975))47lo, hi = beta_hat - z * se_hac, beta_hat + z * se_hac4849sens = {}50for L in sorted({0, nw_rule, 8, 12, hac_lags, 24, 36}):51 s = hac(L)52 sens[str(L)] = {"se": s, "t": beta_hat / s, "p": float(2 * stats.norm.sf(abs(beta_hat / s)))}5354res = {55 "task_id": "F2-conf-02",56 "beta_hat": beta_hat,57 "alpha_hat": alpha_hat,58 "se_iid": se_iid,59 "t_iid": beta_hat / se_iid,60 "p_iid": float(2 * stats.t.sf(abs(beta_hat / se_iid), T - 2)),61 "rho_resid": rho,62 "rho_score": rv,63 "andrews_bandwidth_raw": S_andrews,64 "nw_rule_of_thumb_lags": nw_rule,65 "se_hac": se_hac,66 "hac_lags": hac_lags,67 "t_hac": t_hac,68 "p_hac": p_hac,69 "beta_ci_low": lo,70 "beta_ci_high": hi,71 "hac_sensitivity": sens,72 "memo_claim_supported": bool(p_hac < 0.05),73 "n_analysed": int(T),74 "feasibility": "complete",75 "feasibility_reason": None,76 "data_issues": issues,77 "n_input_rows": {"series.csv": int(n_raw)},78 "random_seed": None,79}80with open("out/results.json", "w") as fh:81 json.dump(res, fh, indent=2)82print(json.dumps({k: res[k] for k in ["beta_hat", "se_iid", "rho_resid", "hac_lags", "se_hac", "t_hac", "p_hac", "memo_claim_supported"]}, indent=2))83