1"""F2-conf-03: Monte Carlo coverage of OLS slope 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, R, LEVEL = 1.0, 0.5, 2000, 0.951112d = pd.read_csv("data/design.csv")13n_rows = len(d)14issues = []15assert n_rows == 30 and d[["x", "sigma"]].notna().all().all() and (d.sigma > 0).all()16x = d["x"].to_numpy(float); sig = d["sigma"].to_numpy(float); n = len(x)17X = np.column_stack([np.ones(n), x])18XtX_inv = np.linalg.inv(X.T @ X)19H = XtX_inv @ X.T # 2 x n, betahat = H y20tcrit = stats.t.ppf(0.5 + LEVEL / 2, n - 2)2122# (c) exact sd of betahat: Var = sum w_i^2 sigma_i^2, w = H[1]23w = H[1]24sd_exact = float(np.sqrt(np.sum(w**2 * sig**2)))2526# Monte Carlo (vectorised): Y is R x n27rng = np.random.default_rng(SEED)28eps = rng.standard_normal((R, n))29Y = ALPHA + BETA * x[None, :] + sig[None, :] * eps30B = Y @ H.T # R x 231resid = Y - B @ X.T32s2 = (resid**2).sum(1) / (n - 2)33se_conv = np.sqrt(s2 * XtX_inv[1, 1])34# HC0 slope variance = sum_i w_i^2 e_i^2 ; HC1 = n/(n-2) * HC035se_hc1 = np.sqrt(n / (n - 2) * (resid**2 @ (w**2)))36bhat = B[:, 1]37cov_conv = float(np.mean(np.abs(bhat - BETA) <= tcrit * se_conv))38cov_hc1 = float(np.mean(np.abs(bhat - BETA) <= tcrit * se_hc1))39mcse_nominal = float(np.sqrt(LEVEL * (1 - LEVEL) / R))40undercovers = bool(cov_conv < LEVEL - 3 * mcse_nominal)4142res = {43 "task_id": "F2-conf-03",44 "coverage_conventional": cov_conv,45 "coverage_hc1": cov_hc1,46 "sd_beta_exact": sd_exact,47 "sd_beta_mc": float(np.std(bhat, ddof=1)),48 "mean_se_conventional": float(se_conv.mean()),49 "conventional_ci_undercovers": undercovers,50 "n_replications": R,51 "n_obs_per_replication": n,52 "feasibility": "complete",53 "feasibility_reason": None,54 "data_issues": issues,55 "n_input_rows": {"design.csv": n_rows},56 "random_seed": SEED,57 # supplementary diagnostics58 "t_critical_28": float(tcrit),59 "mcse_coverage_conventional": float(np.sqrt(cov_conv * (1 - cov_conv) / R)),60 "mcse_coverage_hc1": float(np.sqrt(cov_hc1 * (1 - cov_hc1) / R)),61 "mcse_at_nominal_095": mcse_nominal,62 "undercoverage_threshold": LEVEL - 3 * mcse_nominal,63 "mean_se_hc1": float(se_hc1.mean()),64 "mean_beta_hat": float(bhat.mean()),65 # sqrt(E[s^2] * (X'X)^-1_11), E[s^2] = sum sigma_i^2 (1 - h_ii) / (n-2); >= E[se_conv] by Jensen66 "rms_se_conventional_exact": float(np.sqrt(XtX_inv[1, 1] * np.sum(sig**2 * (1 - np.einsum("ij,ji->i", X, H))) / (n - 2))),67}68os.makedirs("out", exist_ok=True)69with open("out/results.json", "w") as f:70 json.dump(res, f, indent=2)71print(json.dumps(res, indent=2))72