1"""F2-conf-02: OLS slope of y on x with Newey-West HAC SE. Run from workspace root."""2import json, math3import numpy as np, pandas as pd, statsmodels.api as sm4from scipy import stats56d = pd.read_csv("data/series.csv")7n_rows = len(d)8issues = []9# data checks10months = pd.to_datetime(d["month"], format="%Y-%m")11assert d["month"].is_unique and d[["x", "y"]].notna().all().all()12gaps = (months.dt.year * 12 + months.dt.month).diff().dropna()13assert (gaps == 1).all(), "non-consecutive months"14d = d.sort_values("month").reset_index(drop=True)1516X = sm.add_constant(d["x"])17ols = sm.OLS(d["y"], X).fit()18alpha_hat, beta_hat = float(ols.params["const"]), float(ols.params["x"])19se_iid = float(ols.bse["x"])20e = ols.resid.values21rho_resid = float(np.sum(e[1:] * e[:-1]) / np.sum(e ** 2)) # lag-1 autocorrelation (mean-zero residuals)2223# Bandwidth: Andrews (1991) AR(1) plug-in for the Bartlett kernel, applied to the24# slope score v_t = (x_t - xbar) * u_t, which is what the long-run variance concerns.25T = len(d)26v = (d["x"].values - d["x"].values.mean()) * e27rho_v = float(np.sum(v[1:] * v[:-1]) / np.sum(v[:-1] ** 2))28a1 = 4 * rho_v ** 2 / ((1 - rho_v) ** 2 * (1 + rho_v) ** 2)29hac_lags = int(math.floor(1.1447 * (a1 * T) ** (1 / 3)))30nw_rule = int(math.floor(4 * (T / 100) ** (2 / 9))) # Newey-West (1994) rule of thumb, for sensitivity3132def hac(L):33 r = ols.get_robustcov_results(cov_type="HAC", maxlags=L, use_correction=True)34 return float(r.bse[1])3536se_hac = hac(hac_lags)37t_hac = beta_hat / se_hac38p_hac = float(2 * stats.norm.sf(abs(t_hac))) # asymptotic normal reference39z = stats.norm.ppf(0.975)40ci = (beta_hat - z * se_hac, beta_hat + z * se_hac)4142sens = {}43for L in sorted(set([0, nw_rule, hac_lags, 12, 24, 36])):44 s = hac(L)45 sens[str(L)] = {"se": s, "t": beta_hat / s, "p_normal": float(2 * stats.norm.sf(abs(beta_hat / s))),46 "p_t": float(2 * stats.t.sf(abs(beta_hat / s), T - 2))}4748res = {49 "task_id": "F2-conf-02",50 "beta_hat": beta_hat, "alpha_hat": alpha_hat, "se_iid": se_iid,51 "t_iid": beta_hat / se_iid, "p_iid": float(ols.pvalues["x"]),52 "rho_resid": rho_resid, "rho_score": rho_v,53 "se_hac": se_hac, "hac_lags": hac_lags, "t_hac": t_hac, "p_hac": p_hac,54 "p_hac_t_dist": float(2 * stats.t.sf(abs(t_hac), T - 2)),55 "beta_ci_low": ci[0], "beta_ci_high": ci[1],56 "memo_claim_supported": bool(p_hac < 0.05),57 "n_analysed": int(T),58 "feasibility": "complete", "feasibility_reason": None,59 "data_issues": issues,60 "n_input_rows": {"series.csv": n_rows},61 "random_seed": None,62 "hac_sensitivity_by_lags": sens, "nw_rule_of_thumb_lags": nw_rule,63}64with open("out/results.json", "w") as f:65 json.dump(res, f, indent=2)66print(json.dumps(res, indent=2))67