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, N, R = 1.0, 0.5, 30, 2000910d = pd.read_csv("data/design.csv")11issues = []12n_rows = len(d)13if d[["x", "sigma"]].isna().any().any():14 issues.append({"issue": "missing values", "handling": "none dropped; would fail"})15if (d["sigma"] <= 0).any():16 issues.append({"issue": "non-positive sigma", "handling": "flagged"})17if n_rows != N:18 issues.append({"issue": f"expected {N} rows, got {n_rows}", "handling": "flagged"})19x = d["x"].to_numpy(float); sig = d["sigma"].to_numpy(float)20n = len(x)2122X = np.column_stack([np.ones(n), x])23XtXi = np.linalg.inv(X.T @ X)24A = XtXi @ X.T # beta_hat = A y25a = A[1] # weights for slope26sd_exact = float(np.sqrt(np.sum(a**2 * sig**2))) # Var = sum w_i^2 sigma_i^227sd_naive_formula = float(np.sqrt(np.mean(sig**2) * XtXi[1, 1])) # for reference only2829tcrit = stats.t.ppf(0.975, n - 2)30rng = np.random.default_rng(SEED)31eps = rng.standard_normal((R, n))32Y = ALPHA + BETA * x + sig * eps # R x n33B = Y @ A.T # R x 234resid = Y - B @ X.T35s2 = (resid**2).sum(1) / (n - 2)36se_conv = np.sqrt(s2 * XtXi[1, 1])37hc0 = (resid**2) @ (a**2) # slope element of sandwich38se_hc1 = np.sqrt(hc0 * n / (n - 2))39bh = B[:, 1]40cov_conv = float(np.mean(np.abs(bh - BETA) <= tcrit * se_conv))41cov_hc1 = float(np.mean(np.abs(bh - BETA) <= tcrit * se_hc1))42mcse_conv = float(np.sqrt(cov_conv * (1 - cov_conv) / R))43mcse_hc1 = float(np.sqrt(cov_hc1 * (1 - cov_hc1) / R))44mcse_nominal = float(np.sqrt(0.95 * 0.05 / R))45# Decision rule: binomial MC SE at R (evaluated at the estimated coverage)46undercovers = bool(cov_conv < 0.95 - 3 * mcse_conv)4748res = {49 "task_id": "F2-conf-03",50 "coverage_conventional": cov_conv,51 "coverage_hc1": cov_hc1,52 "sd_beta_exact": sd_exact,53 "sd_beta_mc": float(np.std(bh, ddof=1)),54 "mean_se_conventional": float(se_conv.mean()),55 "conventional_ci_undercovers": undercovers,56 "n_replications": R,57 "n_obs_per_replication": n,58 "feasibility": "complete",59 "feasibility_reason": None,60 "data_issues": issues,61 "n_input_rows": {"design.csv": n_rows},62 "random_seed": SEED,63}64os.makedirs("out", exist_ok=True)65with open("out/results.json", "w") as f:66 json.dump(res, f, indent=2)67extra = {"mcse_coverage_conventional": mcse_conv, "mcse_coverage_hc1": mcse_hc1,68 "mcse_at_nominal_0.95": mcse_nominal,69 "undercovers_using_nominal_mcse": bool(cov_conv < 0.95 - 3 * mcse_nominal),70 "t_crit": float(tcrit), "mean_se_hc1": float(se_hc1.mean()),71 "homoskedastic_formula_sd_with_mean_sigma2": sd_naive_formula,72 "mean_beta_hat": float(bh.mean())}73with open("out/diagnostics.json", "w") as f:74 json.dump(extra, f, indent=2)75print(json.dumps(res, indent=2)); print(json.dumps(extra, indent=2))76