1"""F2-conf-03: Monte Carlo coverage of OLS slope CIs under heteroskedasticity.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/design.csv; writes out/results.json deterministically."""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")13issues = []14assert d.shape[0] == 30 and d[["x", "sigma"]].notna().all().all() and (d.sigma > 0).all()15x = d["x"].to_numpy(float); s = d["sigma"].to_numpy(float); n = len(x); df = n - 216tcrit = stats.t.ppf(0.5 + LEVEL / 2, df)1718xc = x - x.mean(); Sxx = np.sum(xc ** 2); w = xc / Sxx19sd_exact = float(np.sqrt(np.sum(w ** 2 * s ** 2))) # (c) analytic2021rng = np.random.default_rng(SEED)22eps = rng.standard_normal((R, n))23Y = ALPHA + BETA * x + s * eps # R x n24bhat = (Y - Y.mean(axis=1, keepdims=True)) @ xc / Sxx25ahat = Y.mean(axis=1) - bhat * x.mean()26e = Y - ahat[:, None] - bhat[:, None] * x27se_conv = np.sqrt(np.sum(e ** 2, axis=1) / df / Sxx)28hc0 = np.sum((xc ** 2) * e ** 2, axis=1) / Sxx ** 229se_hc1 = np.sqrt(hc0 * n / df)3031cov_conv = float(np.mean(np.abs(bhat - BETA) <= tcrit * se_conv))32cov_hc1 = float(np.mean(np.abs(bhat - BETA) <= tcrit * se_hc1))33mcse_nominal = float(np.sqrt(LEVEL * (1 - LEVEL) / R))34undercovers = bool((LEVEL - cov_conv) > 3 * mcse_nominal)3536res = {37 "task_id": "F2-conf-03",38 "coverage_conventional": cov_conv,39 "coverage_hc1": cov_hc1,40 "sd_beta_exact": sd_exact,41 "sd_beta_mc": float(np.std(bhat, ddof=1)),42 "mean_se_conventional": float(np.mean(se_conv)),43 "conventional_ci_undercovers": undercovers,44 "n_replications": R,45 "n_obs_per_replication": n,46 "feasibility": "complete",47 "feasibility_reason": None,48 "data_issues": issues,49 "n_input_rows": {"design.csv": int(d.shape[0])},50 "random_seed": SEED,51 # supplementary diagnostics52 "t_crit_28": float(tcrit),53 "mcse_at_nominal": mcse_nominal,54 "undercoverage_threshold": float(LEVEL - 3 * mcse_nominal),55 "mcse_coverage_conventional": float(np.sqrt(cov_conv * (1 - cov_conv) / R)),56 "mcse_coverage_hc1": float(np.sqrt(cov_hc1 * (1 - cov_hc1) / R)),57 "mean_se_hc1": float(np.mean(se_hc1)),58 "mean_bhat_mc": float(np.mean(bhat)),59 "sd_beta_homoskedastic_formula": float(np.sqrt(np.mean(s ** 2) / Sxx)),60}61os.makedirs("out", exist_ok=True)62with open("out/results.json", "w") as f:63 json.dump(res, f, indent=2)64print(json.dumps(res, indent=2))65