1"""F2-conf-02: OLS slope of y on x with HAC (Newey-West, Bartlett) SE. Run from workspace root."""2import json, numpy as np, pandas as pd, statsmodels.api as sm3from scipy import stats45d = pd.read_csv("data/series.csv")6n_rows = len(d)7issues = []8per = pd.PeriodIndex(d["month"], freq="M")9assert d[["x", "y"]].notna().all().all(), "missing values"10assert not per.duplicated().any()11assert ((per[1:] - per[:-1]).map(lambda o: o.n) == 1).all(), "gaps"12d = d.sort_values("month").reset_index(drop=True)13T = len(d)14X = sm.add_constant(d["x"].values); y = d["y"].values15ols = sm.OLS(y, X).fit()16a, b = ols.params; se_iid = ols.bse[1]17u = ols.resid18rho_resid = float(np.corrcoef(u[1:], u[:-1])[0, 1])1920# Bandwidth: Andrews (1991) AR(1) plug-in for the Bartlett kernel, applied to the21# slope score v_t = (x_t - xbar) * u_t. This adapts to the stated persistence (rho~0.7)22# unlike the fixed Newey-West (1994) rule-of-thumb floor(4(T/100)^(2/9)).23v = (d["x"].values - d["x"].mean()) * u24rv = float(np.sum(v[1:] * v[:-1]) / np.sum(v[:-1] ** 2))25alpha1 = 4 * rv**2 / ((1 - rv) ** 2 * (1 + rv) ** 2)26andrews_bw = 1.1447 * (alpha1 * T) ** (1 / 3)27hac_lags = int(np.ceil(andrews_bw))28nw_rule = int(np.floor(4 * (T / 100) ** (2 / 9)))2930def hac(L):31 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": True, "kernel": "bartlett"})32 return float(r.bse[1])3334se_hac = hac(hac_lags)35t_hac = b / se_hac36p_hac = 2 * stats.norm.sf(abs(t_hac))37z = stats.norm.ppf(0.975)38sens = {}39for L in sorted({0, nw_rule, 8, 12, hac_lags, 18, 24}):40 s = hac(L); sens[str(L)] = {"se": s, "t": b / s, "p": float(2 * stats.norm.sf(abs(b / s)))}4142res = {43 "task_id": "F2-conf-02",44 "beta_hat": float(b), "alpha_hat": float(a), "se_iid": float(se_iid),45 "t_iid": float(ols.tvalues[1]), "p_iid": float(ols.pvalues[1]),46 "rho_resid": rho_resid, "se_hac": se_hac, "hac_lags": hac_lags,47 "t_hac": float(t_hac), "p_hac": float(p_hac),48 "beta_ci_low": float(b - z * se_hac), "beta_ci_high": float(b + z * se_hac),49 "memo_claim_supported": bool(p_hac < 0.05),50 "n_analysed": int(T), "feasibility": "complete", "feasibility_reason": None,51 "data_issues": issues, "n_input_rows": {"series.csv": int(n_rows)}, "random_seed": None,52 "bandwidth_details": {"rule": "Andrews (1991) AR(1) plug-in, Bartlett kernel, on slope score",53 "score_ar1": rv, "andrews_bw_raw": andrews_bw, "nw1994_rule_of_thumb": nw_rule},54 "hac_lag_sensitivity": sens,55}56json.dump(res, open("out/results.json", "w"), indent=2)57print(json.dumps(res, indent=2))58