1"""F2-conf-02: OLS slope of y on x with Newey-West HAC standard errors.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json4import numpy as np5import pandas as pd6import statsmodels.api as sm7from scipy import stats89d = pd.read_csv("data/series.csv")10n_rows = len(d)11issues = []12# data checks13m = pd.PeriodIndex(d["month"], freq="M")14gaps = int(((m[1:] - m[:-1]).map(lambda o: o.n) != 1).sum())15if d[["x", "y"]].isna().any().any():16 issues.append({"issue": "missing values", "handling": "rows dropped"})17 d = d.dropna(subset=["x", "y"])18if d["month"].duplicated().any():19 issues.append({"issue": "duplicate months", "handling": "first kept"})20 d = d.drop_duplicates("month")21if gaps:22 issues.append({"issue": f"{gaps} gaps in monthly index", "handling": "kept; HAC lags count observations"})2324X = sm.add_constant(d["x"].to_numpy())25y = d["y"].to_numpy()26n = len(y)27ols = sm.OLS(y, X).fit()28alpha_hat, beta_hat = ols.params29se_iid = ols.bse[1]30u = ols.resid31rho_resid = float(np.sum(u[1:] * u[:-1]) / np.sum(u ** 2))3233# Bandwidth: Andrews (1991) AR(1) plug-in for the Bartlett kernel, applied to the34# score series v_t = (x_t - xbar) u_t, which is what the HAC long-run variance targets.35v = (d["x"].to_numpy() - d["x"].mean()) * u36rv = float(np.sum((v[1:] - v.mean()) * (v[:-1] - v.mean())) / np.sum((v - v.mean()) ** 2))37a1 = 4 * rv ** 2 / ((1 - rv) ** 2 * (1 + rv) ** 2)38hac_lags = int(np.floor(1.1447 * (a1 * n) ** (1 / 3)))39nw_rule = int(np.floor(4 * (n / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb, for sensitivity4041def hac(L):42 r = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "kernel": "bartlett", "use_correction": True})43 return float(r.bse[1])4445se_hac = hac(hac_lags)46t_hac = beta_hat / se_hac47p_hac = 2 * stats.norm.sf(abs(t_hac))48z = stats.norm.ppf(0.975)49sens = {}50for L in sorted(set([nw_rule, hac_lags, 12, 18, 24])):51 s = hac(L)52 sens[str(L)] = {"se_hac": s, "t": beta_hat / s, "p": float(2 * stats.norm.sf(abs(beta_hat / s)))}5354res = {55 "task_id": "F2-conf-02",56 "beta_hat": float(beta_hat),57 "alpha_hat": float(alpha_hat),58 "se_iid": float(se_iid),59 "t_iid": float(ols.tvalues[1]),60 "p_iid": float(ols.pvalues[1]),61 "rho_resid": rho_resid,62 "rho_score": rv,63 "se_hac": se_hac,64 "hac_lags": hac_lags,65 "hac_lags_nw_rule_of_thumb": nw_rule,66 "t_hac": float(t_hac),67 "p_hac": float(p_hac),68 "beta_ci_low": float(beta_hat - z * se_hac),69 "beta_ci_high": float(beta_hat + z * se_hac),70 "memo_claim_supported": bool(p_hac < 0.05),71 "hac_sensitivity": sens,72 "n_analysed": int(n),73 "feasibility": "complete",74 "feasibility_reason": None,75 "data_issues": issues,76 "n_input_rows": {"series.csv": int(n_rows)},77 "random_seed": None,78}79with open("out/results.json", "w") as f:80 json.dump(res, f, indent=2)81print(json.dumps(res, indent=2))82