1"""F2-conf-02: OLS slope of y on x with Newey-West (Bartlett) HAC inference.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/series.csv; writes out/results.json and out/exclusions.csv. No randomness."""4import json, math, os5import numpy as np, pandas as pd6from scipy import stats7import statsmodels.api as sm89DATA = "data/series.csv"10os.makedirs("out", exist_ok=True)1112raw = pd.read_csv(DATA, dtype={"month": str})13n_raw = len(raw)14issues = []15excl = []1617# Audit: missing values, duplicates, gaps18d = raw.copy()19nb = len(d); d = d.dropna(subset=["x", "y"])20excl.append(("1", "drop rows with missing x or y", nb, len(d), nb - len(d), "OLS requires complete cases"))21nb = len(d); d = d.drop_duplicates(subset=["month"])22excl.append(("2", "drop duplicated months", nb, len(d), nb - len(d), "one observation per month"))23per = pd.PeriodIndex(d["month"], freq="M")24gaps = int(((per[1:] - per[:-1]).map(lambda z: z.n) != 1).sum())25if gaps:26 issues.append({"issue": f"{gaps} gaps/non-consecutive months", "handling": "kept; HAC lags computed on row order"})27if n_raw - len(d):28 issues.append({"issue": f"{n_raw-len(d)} rows removed (missing/duplicate)", "handling": "excluded, see out/exclusions.csv"})29d = d.sort_values("month").reset_index(drop=True)3031y = d["y"].to_numpy(float); x = d["x"].to_numpy(float); T = len(y)32X = np.column_stack([np.ones(T), x])33XtX_inv = np.linalg.inv(X.T @ X)34b = XtX_inv @ X.T @ y35u = y - X @ b36k = X.shape[1]3738# Naive iid SE (memo's measure)39s2 = u @ u / (T - k)40se_iid = math.sqrt(s2 * XtX_inv[1, 1])4142# First-order residual autocorrelation (OLS of u_t on u_{t-1}, no intercept; residuals have mean 0)43rho_resid = float(np.sum(u[1:] * u[:-1]) / np.sum(u[:-1] ** 2))4445def nw_se(L):46 """Newey-West Bartlett-kernel HAC SE of the slope, no small-sample df correction."""47 g = X * u[:, None]48 S = g.T @ g49 for j in range(1, L + 1):50 w = 1 - j / (L + 1)51 G = g[j:].T @ g[:-j]52 S += w * (G + G.T)53 V = XtX_inv @ S @ XtX_inv54 return math.sqrt(V[1, 1])5556# Bandwidth: Andrews (1991) AR(1) plug-in for the Bartlett kernel, applied to the slope score57# v_t = (x_t - xbar) u_t: S_T = 1.1447 * (alpha(1) * T)^(1/3), alpha(1) = 4 r^2 / ((1-r)^2 (1+r)^2)58v = (x - x.mean()) * u59r_v = float(np.sum(v[1:] * v[:-1]) / np.sum(v[:-1] ** 2))60alpha1 = 4 * r_v ** 2 / ((1 - r_v) ** 2 * (1 + r_v) ** 2)61andrews_bw = 1.1447 * (alpha1 * T) ** (1 / 3)62hac_lags = int(math.floor(andrews_bw))63nw_rule = int(math.floor(4 * (T / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb, for sensitivity6465se_hac = nw_se(hac_lags)66t_hac = float(b[1] / se_hac)67p_hac = float(2 * stats.norm.sf(abs(t_hac)))68z = stats.norm.ppf(0.975)69ci_low, ci_high = float(b[1] - z * se_hac), float(b[1] + z * se_hac)7071# Cross-check against statsmodels (use_correction=False -> no df correction)72sm_fit = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": hac_lags, "use_correction": False})73assert abs(sm_fit.bse[1] - se_hac) < 1e-10, (sm_fit.bse[1], se_hac)7475sens = {}76for L in sorted(set([0, nw_rule, hac_lags, 12, 24])):77 s = nw_se(L); t = b[1] / s78 sens[f"lags_{L}"] = {"se": s, "t": float(t), "p_normal": float(2 * stats.norm.sf(abs(t))),79 "p_t_df": float(2 * stats.t.sf(abs(t), T - k))}80p_iid = float(2 * stats.t.sf(abs(b[1] / se_iid), T - k))8182memo_claim_supported = bool(p_hac < 0.05)8384res = {85 "task_id": "F2-conf-02",86 "beta_hat": float(b[1]), "alpha_hat": float(b[0]),87 "se_iid": se_iid, "t_iid": float(b[1] / se_iid), "p_iid": p_iid,88 "rho_resid": rho_resid,89 "se_hac": se_hac, "hac_lags": hac_lags, "t_hac": t_hac, "p_hac": p_hac,90 "beta_ci_low": ci_low, "beta_ci_high": ci_high,91 "memo_claim_supported": memo_claim_supported,92 "n_analysed": int(T),93 "feasibility": "complete", "feasibility_reason": None,94 "data_issues": issues,95 "n_input_rows": {"series.csv": int(n_raw)},96 "random_seed": None,97 "bandwidth_rule": "Andrews (1991) AR(1) plug-in, Bartlett kernel, on slope score (x-xbar)*u; floor",98 "andrews_bandwidth_raw": float(andrews_bw), "score_rho": r_v, "nw_rule_of_thumb_lags": nw_rule,99 "hac_sensitivity": sens,100}101with open("out/results.json", "w") as f:102 json.dump(res, f, indent=2)103pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv("out/exclusions.csv", index=False)104print(f"beta={b[1]:.4f} se_iid={se_iid:.4f} se_hac({hac_lags})={se_hac:.4f} t={t_hac:.3f} p={p_hac:.4f} supported={memo_claim_supported}")105