1"""F2-conf-02: OLS slope of y on x with Newey-West HAC inference. Run from workspace root."""2import json, math, os3import numpy as np, pandas as pd4import statsmodels.api as sm5from scipy import stats67d = pd.read_csv("data/series.csv")8n_raw = len(d)9d = d.dropna(subset=["x", "y"])10months = pd.PeriodIndex(d["month"], freq="M")11assert (np.diff(months.asi8) == 1).all(), "month gaps"12n = len(d)13X = sm.add_constant(d["x"].to_numpy()); y = d["y"].to_numpy()1415# P2 baseline: conventional iid OLS (memo's measure)16ols = sm.OLS(y, X).fit()17alpha_hat, beta_hat = map(float, ols.params)18se_iid = float(ols.bse[1]); t_iid = beta_hat / se_iid19p_iid = float(2 * stats.t.sf(abs(t_iid), n - 2))20u = ols.resid21rho_resid = float(np.sum(u[1:] * u[:-1]) / np.sum(u * u))2223# Bandwidth: Andrews (1991) AR(1) plug-in for the Bartlett kernel applied to the24# slope score v_t = (x_t - xbar) * u_t: S_T = 1.1447 * (a1 * T)^(1/3),25# a1 = 4 rho^2 / ((1-rho)^2 (1+rho)^2). Rounded up to an integer.26v = (d["x"].to_numpy() - d["x"].mean()) * u27rho_v = float(np.sum(v[1:] * v[:-1]) / np.sum(v * v))28a1 = 4 * rho_v**2 / ((1 - rho_v) ** 2 * (1 + rho_v) ** 2)29andrews_bw = 1.1447 * (a1 * n) ** (1 / 3)30hac_lags = int(math.ceil(andrews_bw))31nw_rule = int(math.floor(4 * (n / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb3233def hac(L):34 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": True})35 return float(r.bse[1])3637se_hac = hac(hac_lags)38t_hac = beta_hat / se_hac39p_hac = float(2 * stats.norm.sf(abs(t_hac)))40z = stats.norm.ppf(0.975)41ci = (beta_hat - z * se_hac, beta_hat + z * se_hac)4243sens = {}44for L in sorted({nw_rule, 6, 8, 12, 18, 24, hac_lags}):45 s = hac(L); sens[str(L)] = {"se_hac": s, "t_hac": beta_hat / s, "p_hac": float(2 * stats.norm.sf(abs(beta_hat / s)))}4647res = {48 "task_id": "F2-conf-02",49 "beta_hat": beta_hat, "alpha_hat": alpha_hat, "se_iid": se_iid,50 "baseline_t_iid": t_iid, "baseline_p_iid": p_iid,51 "rho_resid": rho_resid, "rho_score": rho_v, "andrews_bandwidth_raw": andrews_bw,52 "nw_rule_of_thumb_lags": nw_rule,53 "se_hac": se_hac, "hac_lags": hac_lags, "t_hac": t_hac, "p_hac": p_hac,54 "beta_ci_low": ci[0], "beta_ci_high": ci[1],55 "se_ratio_hac_to_iid": se_hac / se_iid,56 "memo_claim_supported": bool(p_hac < 0.05),57 "hac_sensitivity": sens,58 "n_analysed": int(n), "feasibility": "complete", "feasibility_reason": None,59 "data_issues": [], "n_input_rows": {"series.csv": int(n_raw)}, "random_seed": None,60}61os.makedirs("out", exist_ok=True)62with open("out/results.json", "w") as f:63 json.dump(res, f, indent=2)64print(json.dumps(res, indent=2))65