1"""F2-conf-02: re-estimate OLS slope of y on x with a HAC (Newey-West, Bartlett) SE.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, math5import numpy as np, pandas as pd, statsmodels.api as sm6from scipy import stats78df = pd.read_csv("data/series.csv")9n_raw = len(df)10issues = []11excl = []12# --- audit / filters (none expected to remove rows) ---13n0 = len(df); df = df.dropna(subset=["x", "y"]); excl.append(("1", "drop rows with missing x or y", n0, len(df), n0 - len(df), "OLS requires complete cases"))14n0 = len(df); df = df.drop_duplicates(subset=["month"]); excl.append(("2", "drop duplicate months", n0, len(df), n0 - len(df), "one observation per month"))15per = pd.PeriodIndex(df["month"], freq="M")16gaps = int((np.diff(per.asi8) != 1).sum())17if gaps:18 issues.append({"issue": f"{gaps} non-consecutive month steps", "handling": "kept; HAC lags interpreted as observation lags"})19pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv("out/exclusions.csv", index=False)2021y = df["y"].to_numpy(float); x = df["x"].to_numpy(float); T = len(y)22X = sm.add_constant(x)23ols = sm.OLS(y, X).fit()24alpha_hat, beta_hat = float(ols.params[0]), float(ols.params[1])25se_iid = float(ols.bse[1])26u = ols.resid27rho_resid = float(np.sum(u[1:] * u[:-1]) / np.sum(u * u))2829# --- bandwidth: Andrews (1991) AR(1) plug-in for the Bartlett kernel, applied to the slope score30v = (x - x.mean()) * u31rv = float(np.sum(v[1:] * v[:-1]) / np.sum(v[:-1] ** 2))32a1 = 4 * rv**2 / ((1 - rv) ** 2 * (1 + rv) ** 2)33andrews_bw = 1.1447 * (a1 * T) ** (1 / 3)34hac_lags = int(math.floor(andrews_bw))35nw_rule_lags = int(math.floor(4 * (T / 100) ** (2 / 9)))3637def hac_fit(L):38 return sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": False})3940hac = hac_fit(hac_lags)41se_hac = float(hac.bse[1])4243# manual cross-check of the sandwich (Bartlett weights 1 - j/(L+1))44Z = X * u[:, None]45S = Z.T @ Z46for j in range(1, hac_lags + 1):47 G = Z[j:].T @ Z[:-j]48 S += (1 - j / (hac_lags + 1)) * (G + G.T)49XtXi = np.linalg.inv(X.T @ X)50se_manual = float(np.sqrt((XtXi @ S @ XtXi)[1, 1]))51assert abs(se_manual - se_hac) / se_hac < 1e-8, (se_manual, se_hac)5253t_hac = beta_hat / se_hac54p_hac = float(2 * stats.norm.sf(abs(t_hac)))55z = float(stats.norm.ppf(0.975))56ci_low, ci_high = beta_hat - z * se_hac, beta_hat + z * se_hac5758sens = {}59for L in sorted(set([0, nw_rule_lags, hac_lags, 12, 16, 24])):60 f = hac_fit(L); se = float(f.bse[1]); t = beta_hat / se61 sens[str(L)] = {"se": se, "t": t, "p": float(2 * stats.norm.sf(abs(t)))}6263res = {64 "task_id": "F2-conf-02",65 "beta_hat": beta_hat, "alpha_hat": alpha_hat,66 "se_iid": se_iid, "t_iid": float(ols.tvalues[1]), "p_iid": float(ols.pvalues[1]),67 "rho_resid": rho_resid,68 "se_hac": se_hac, "hac_lags": hac_lags, "t_hac": float(t_hac), "p_hac": p_hac,69 "beta_ci_low": float(ci_low), "beta_ci_high": float(ci_high),70 "memo_claim_supported": bool(p_hac < 0.05),71 "n_analysed": int(T),72 "feasibility": "complete", "feasibility_reason": None,73 "data_issues": issues,74 "n_input_rows": {"series.csv": int(n_raw)},75 "random_seed": None,76 "baseline_se_iid_t_p": [se_iid, float(ols.tvalues[1]), float(ols.pvalues[1])],77 "hac_bandwidth_details": {"method": "Andrews (1991) AR(1) plug-in, Bartlett kernel, on score (x-xbar)*u",78 "rho_score": rv, "andrews_bw_continuous": andrews_bw,79 "nw_rule_of_thumb_lags": nw_rule_lags, "small_sample_correction": False,80 "reference_distribution": "standard normal"},81 "hac_lag_sensitivity": sens,82 "hac_lag_sensitivity_min_p_excl_0": min(v_["p"] for k, v_ in sens.items() if k != "0"),83}84with open("out/results.json", "w") as fh:85 json.dump(res, fh, indent=2)86print(f"beta={beta_hat:.4f} se_iid={se_iid:.4f} se_hac={se_hac:.4f} (L={hac_lags}) t={t_hac:.3f} p={p_hac:.4f} "87 f"CI=[{ci_low:.4f},{ci_high:.4f}] memo_supported={p_hac < 0.05}")88