1"""F2-conf-02: OLS slope of y on x with Newey-West HAC inference.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/series.csv (memo numbers are checked against, not used). No randomness."""4import json, math5import numpy as np, pandas as pd, statsmodels.api as sm6from scipy import stats78d = pd.read_csv("data/series.csv")9n_raw = len(d)10issues = []11# --- data checks (no exclusions needed if all pass) ---12per = pd.PeriodIndex(d["month"], freq="M")13gaps = int(((per[1:].astype("int64") - per[:-1].astype("int64")) != 1).sum())14n_dup = int(d["month"].duplicated().sum())15n_na = int(d[["x", "y"]].isna().any(axis=1).sum())16if gaps: issues.append({"issue": f"{gaps} non-consecutive month steps", "handling": "kept; HAC assumes equal spacing"})17if n_dup: issues.append({"issue": f"{n_dup} duplicate months", "handling": "reported"})18if n_na:19 issues.append({"issue": f"{n_na} rows with missing x/y", "handling": "dropped"})20 d = d.dropna(subset=["x", "y"])21d = d.sort_values("month").reset_index(drop=True)22T = len(d)2324X = sm.add_constant(d["x"].values)25y = d["y"].values26ols = sm.OLS(y, X).fit() # conventional iid SE (memo baseline)27alpha_hat, beta_hat = map(float, ols.params)28se_iid = float(ols.bse[1])29u = ols.resid30rho_resid = float(np.sum(u[1:] * u[:-1]) / np.sum(u * u)) # standard first-order autocorrelation coefficient3132# --- bandwidth: Andrews (1991) AR(1) plug-in for Bartlett kernel, applied to the slope score v_t=(x_t-xbar)u_t ---33v = (d["x"].values - d["x"].values.mean()) * u34rho_v = float(np.sum(v[1:] * v[:-1]) / np.sum(v * v))35a1 = 4 * rho_v**2 / ((1 - rho_v) ** 2 * (1 + rho_v) ** 2)36S_andrews = 1.1447 * (a1 * T) ** (1 / 3)37nw_rule = int(math.floor(4 * (T / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb, reported as sensitivity38hac_lags = max(int(math.floor(S_andrews)), nw_rule)3940def hac(L):41 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": True})42 se = float(r.bse[1]); t = beta_hat / se43 return se, t, float(2 * stats.norm.sf(abs(t)))4445se_hac, t_hac, p_hac = hac(hac_lags)46z = stats.norm.ppf(0.975)47ci_lo, ci_hi = beta_hat - z * se_hac, beta_hat + z * se_hac48sens = {str(L): dict(zip(["se", "t", "p"], hac(L))) for L in sorted({nw_rule, 4, 8, 12, 16, 24, hac_lags})}4950res = {51 "task_id": "F2-conf-02",52 "beta_hat": beta_hat, "alpha_hat": alpha_hat, "se_iid": se_iid,53 "t_iid": beta_hat / se_iid, "p_iid": float(ols.pvalues[1]),54 "rho_resid": rho_resid, "se_hac": se_hac, "hac_lags": hac_lags,55 "t_hac": t_hac, "p_hac": p_hac, "beta_ci_low": ci_lo, "beta_ci_high": ci_hi,56 "memo_claim_supported": bool(p_hac < 0.05),57 "n_analysed": int(T), "feasibility": "complete", "feasibility_reason": None,58 "data_issues": issues, "n_input_rows": {"series.csv": int(n_raw)}, "random_seed": None,59 "bandwidth_details": {"rho_score": rho_v, "andrews_alpha1": a1, "andrews_S": S_andrews,60 "nw_rule_of_thumb_lags": nw_rule, "rule": "max(floor(Andrews AR(1) plug-in), floor(4(T/100)^(2/9)))"},61 "hac_sensitivity": sens,62 "se_ratio_hac_over_iid": se_hac / se_iid,63 "memo_reproduced": {"alpha": round(alpha_hat, 4) == 2.0925, "beta": round(beta_hat, 4) == 0.1017,64 "se": round(se_iid, 4) == 0.0448, "t": round(beta_hat / se_iid, 2) == 2.27,65 "p": round(float(ols.pvalues[1]), 4) == 0.0238},66 "any_sensitivity_p_lt_005": bool(min(s["p"] for s in sens.values()) < 0.05),67}68with open("out/results.json", "w") as f:69 json.dump(res, f, indent=2)70print(f"beta={beta_hat:.4f} se_iid={se_iid:.4f} se_hac={se_hac:.4f} (L={hac_lags}) t_hac={t_hac:.3f} p_hac={p_hac:.4f} "71 f"CI=[{ci_lo:.4f},{ci_hi:.4f}] memo_supported={res['memo_claim_supported']}")72