1"""F2-conf-03: Monte Carlo coverage of OLS CIs under heteroskedasticity.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3"""4import json, os5import numpy as np6import pandas as pd7from scipy import stats89SEED = 2026103010ALPHA, BETA, N, R, LEVEL = 1.0, 0.5, 30, 2000, 0.951112d = pd.read_csv("data/design.csv")13n_rows = len(d)14issues = []15assert n_rows == N and d[["x", "sigma"]].notna().all().all() and (d.sigma > 0).all()16x = d["x"].to_numpy(float); sig = d["sigma"].to_numpy(float)1718# exact SD of beta_hat: beta_hat - beta = sum w_i sigma_i eps_i19xc = x - x.mean(); Sxx = np.sum(xc**2); w = xc / Sxx20sd_exact = float(np.sqrt(np.sum(w**2 * sig**2)))2122rng = np.random.default_rng(SEED)23eps = rng.standard_normal((R, N))24Y = ALPHA + BETA * x[None, :] + sig[None, :] * eps25ybar = Y.mean(axis=1, keepdims=True)26b = (Y - ybar) @ xc / Sxx27a = ybar[:, 0] - b * x.mean()28E = Y - a[:, None] - b[:, None] * x[None, :]29s2 = np.sum(E**2, axis=1) / (N - 2)30se_conv = np.sqrt(s2 / Sxx)31se_hc0 = np.sqrt(np.sum((xc[None, :] ** 2) * E**2, axis=1)) / Sxx32se_hc1 = se_hc0 * np.sqrt(N / (N - 2))33tcrit = stats.t.ppf(0.5 + LEVEL / 2, N - 2)34cov_conv = float(np.mean(np.abs(b - BETA) <= tcrit * se_conv))35cov_hc1 = float(np.mean(np.abs(b - BETA) <= tcrit * se_hc1))36mcse0 = np.sqrt(LEVEL * (1 - LEVEL) / R)37undercovers = bool(cov_conv < LEVEL - 3 * mcse0)3839# cross-check one replication against statsmodels40import statsmodels.api as sm41m = sm.OLS(Y[0], sm.add_constant(x)).fit()42assert np.isclose(m.bse[1], se_conv[0]) and np.isclose(m.HC1_se[1], se_hc1[0])4344res = {45 "task_id": "F2-conf-03",46 "coverage_conventional": cov_conv,47 "coverage_hc1": cov_hc1,48 "sd_beta_exact": sd_exact,49 "sd_beta_mc": float(np.std(b, ddof=1)),50 "mean_se_conventional": float(se_conv.mean()),51 "conventional_ci_undercovers": undercovers,52 "n_replications": R,53 "n_obs_per_replication": N,54 "feasibility": "complete",55 "feasibility_reason": None,56 "data_issues": issues,57 "n_input_rows": {"design.csv": int(n_rows)},58 "random_seed": SEED,59 # supplementary60 "t_critical": float(tcrit),61 "mcse_coverage_conventional": float(np.sqrt(cov_conv * (1 - cov_conv) / R)),62 "mcse_coverage_hc1": float(np.sqrt(cov_hc1 * (1 - cov_hc1) / R)),63 "mcse_at_nominal": float(mcse0),64 "undercoverage_threshold": float(LEVEL - 3 * mcse0),65 "mean_se_hc1": float(se_hc1.mean()),66 "sqrt_E_s2_over_Sxx_analytic": float(np.sqrt((np.sum(sig**2) - np.sum((1/N + xc**2/Sxx) * sig**2)) / (N - 2) / Sxx)),67 "baseline_coverage_conventional": cov_conv,68}69os.makedirs("out", exist_ok=True)70with open("out/results.json", "w") as f:71 json.dump(res, f, indent=2)72print(json.dumps(res, indent=2))73