1"""F2-conf-03: Monte Carlo coverage of OLS slope CIs under fixed heteroskedastic design.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3"""4import json, os5import numpy as np6import pandas as pd7from scipy import stats89ALPHA, BETA, N, R, LEVEL, SEED = 1.0, 0.5, 30, 2000, 0.95, 202610301011d = pd.read_csv("data/design.csv")12n_rows = len(d)13issues = []14assert n_rows == N and d[["x", "sigma"]].notna().all().all() and (d.sigma > 0).all()15x = d["x"].to_numpy(float); sig = d["sigma"].to_numpy(float)1617xc = x - x.mean(); Sxx = np.sum(xc**2)18w = xc / Sxx # beta_hat = beta + sum w_i sigma_i eps_i19sd_beta_exact = float(np.sqrt(np.sum(w**2 * sig**2)))20# exact E[s^2] (not E[se]): trace(M Sigma)/(n-2), reported as a diagnostic21X = np.column_stack([np.ones(N), x])22H = X @ np.linalg.solve(X.T @ X, X.T); M = np.eye(N) - H23se_conv_sqrt_Es2 = float(np.sqrt(np.trace(M * sig**2) / (N - 2) / Sxx))2425rng = np.random.default_rng(SEED)26eps = rng.standard_normal((R, N)) # only eps redrawn; x, sigma fixed27Y = ALPHA + BETA * x + sig * eps # R x N28yc = Y - Y.mean(axis=1, keepdims=True)29b = yc @ xc / Sxx30a = Y.mean(axis=1) - b * x.mean()31e = Y - a[:, None] - b[:, None] * x32se_conv = np.sqrt(np.sum(e**2, axis=1) / (N - 2) / Sxx)33hc0 = np.sum((xc**2) * e**2, axis=1) / Sxx**2 # slope element of HC0 sandwich34se_hc1 = np.sqrt(hc0 * N / (N - 2))35tcrit = stats.t.ppf(0.5 + LEVEL / 2, N - 2)36cov_conv = float(np.mean(np.abs(b - BETA) <= tcrit * se_conv))37cov_hc1 = float(np.mean(np.abs(b - BETA) <= tcrit * se_hc1))3839mcse_nominal = float(np.sqrt(LEVEL * (1 - LEVEL) / R))40mcse_conv = float(np.sqrt(cov_conv * (1 - cov_conv) / R))41mcse_hc1 = float(np.sqrt(cov_hc1 * (1 - cov_hc1) / R))42undercovers = bool(cov_conv < LEVEL - 3 * mcse_nominal)4344# cross-check one replication with statsmodels45try:46 import statsmodels.api as sm47 f = sm.OLS(Y[0], X).fit()48 assert abs(f.bse[1] - se_conv[0]) < 1e-1049 assert abs(f.get_robustcov_results("HC1").bse[1] - se_hc1[0]) < 1e-1050 sm_check = True51except ImportError:52 sm_check = None5354res = {55 "task_id": "F2-conf-03",56 "coverage_conventional": cov_conv,57 "coverage_hc1": cov_hc1,58 "sd_beta_exact": sd_beta_exact,59 "sd_beta_mc": float(np.std(b, ddof=1)),60 "mean_se_conventional": float(np.mean(se_conv)),61 "conventional_ci_undercovers": undercovers,62 "n_replications": R,63 "n_obs_per_replication": N,64 "feasibility": "complete",65 "feasibility_reason": None,66 "data_issues": issues,67 "n_input_rows": {"design.csv": n_rows},68 "random_seed": SEED,69 # supplementary70 "t_critical": float(tcrit),71 "mcse_at_nominal": mcse_nominal,72 "mcse_coverage_conventional": mcse_conv,73 "mcse_coverage_hc1": mcse_hc1,74 "undercover_threshold": LEVEL - 3 * mcse_nominal,75 "conventional_ci_undercovers_using_estimated_mcse": bool(cov_conv < LEVEL - 3 * mcse_conv),76 "hc1_ci_undercovers": bool(cov_hc1 < LEVEL - 3 * mcse_nominal),77 "mean_se_hc1": float(np.mean(se_hc1)),78 "mean_beta_hat": float(np.mean(b)),79 "mcse_sd_beta_mc_approx": float(np.std(b, ddof=1) / np.sqrt(2 * (R - 1))),80 "sqrt_expected_conv_variance_exact": se_conv_sqrt_Es2,81 "ratio_mean_se_conv_to_sd_exact": float(np.mean(se_conv) / sd_beta_exact),82 "statsmodels_crosscheck_passed": sm_check,83}84os.makedirs("out", exist_ok=True)85with open("out/results.json", "w") as fh:86 json.dump(res, fh, indent=2)87print(json.dumps(res, indent=2))88