1"""F2-conf-03: Monte Carlo coverage of OLS CIs under heteroskedasticity.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, os4import numpy as np, pandas as pd5from scipy import stats67SEED = 202610308ALPHA, BETA, R, LEVEL = 1.0, 0.5, 2000, 0.95910d = pd.read_csv("data/design.csv")11issues = []12if d.isna().any().any():13 issues.append({"issue": "missing values", "handling": "none expected; would abort"})14x, sig = d["x"].to_numpy(float), d["sigma"].to_numpy(float)15n = len(x)16assert n == 30 and (sig > 0).all()1718X = np.column_stack([np.ones(n), x])19XtX_inv = np.linalg.inv(X.T @ X)20xc = x - x.mean(); Sxx = xc @ xc21w = xc / Sxx # beta_hat = beta + sum w_i sigma_i eps_i22sd_exact = float(np.sqrt(np.sum(w**2 * sig**2)))23tcrit = stats.t.ppf(0.5 + LEVEL / 2, n - 2)2425rng = np.random.default_rng(SEED)26eps = rng.standard_normal((R, n))27Y = ALPHA + BETA * x + sig * eps # R x n28B = Y @ X @ XtX_inv # (XtX)^-1 X'y for each row (XtX_inv symmetric)29E = Y - B @ X.T # residuals30s2 = (E**2).sum(1) / (n - 2)31se_conv = np.sqrt(s2 * XtX_inv[1, 1])32# HC0 sandwich slope variance: sum (row-2 of (XtX)^-1 X')_i^2 e_i^2 ; HC1 = HC0 * n/(n-2)33a = (XtX_inv @ X.T)[1] # equals w34se_hc1 = np.sqrt((a**2 * E**2).sum(1) * n / (n - 2))35bh = B[:, 1]36cov_conv = float(np.mean(np.abs(bh - BETA) <= tcrit * se_conv))37cov_hc1 = float(np.mean(np.abs(bh - BETA) <= tcrit * se_hc1))38mcse_nom = float(np.sqrt(LEVEL * (1 - LEVEL) / R))39undercovers = bool(cov_conv < LEVEL - 3 * mcse_nom)4041res = {42 "task_id": "F2-conf-03",43 "coverage_conventional": cov_conv,44 "coverage_hc1": cov_hc1,45 "sd_beta_exact": sd_exact,46 "sd_beta_mc": float(np.std(bh, ddof=1)),47 "mean_se_conventional": float(se_conv.mean()),48 "conventional_ci_undercovers": undercovers,49 "n_replications": R,50 "n_obs_per_replication": n,51 "feasibility": "complete",52 "feasibility_reason": None,53 "data_issues": issues,54 "n_input_rows": {"design.csv": int(len(d))},55 "random_seed": SEED,56 # supplementary keys57 "t_crit_28": float(tcrit),58 "mcse_binomial_at_nominal": mcse_nom,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 "undercover_threshold": float(LEVEL - 3 * mcse_nom),62 "mean_se_hc1": float(se_hc1.mean()),63 "sd_beta_conventional_formula_homoskedastic": float(np.sqrt(np.mean(sig**2) * XtX_inv[1, 1])),64 "mean_beta_hat_mc": float(bh.mean()),65}66os.makedirs("out", exist_ok=True)67with open("out/results.json", "w") as f:68 json.dump(res, f, indent=2)69print(json.dumps(res, indent=2))70