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, os, csv5import numpy as np6import pandas as pd7from scipy import stats89ALPHA, BETA, N, R, LEVEL = 1.0, 0.5, 30, 2000, 0.9510SEED = 202610301112d = pd.read_csv("data/design.csv")13n_raw = len(d)14issues = []15# audit checks (no filtering needed if all pass)16assert list(d.columns) == ["i", "x", "sigma"]17assert not d.isna().any().any(), "missing values in design"18assert n_raw == N, f"expected {N} rows, got {n_raw}"19assert (d["sigma"] > 0).all()20assert d["i"].tolist() == list(range(1, N + 1))21assert d["x"].std() > 022x = d["x"].to_numpy(float); sig = d["sigma"].to_numpy(float)2324os.makedirs("out", exist_ok=True)25with open("out/exclusions.csv", "w", newline="") as f:26 w = csv.writer(f)27 w.writerow(["step", "rule", "n_before", "n_after", "n_removed", "rationale"])28 w.writerow(["1", "drop rows with missing/invalid x or sigma", n_raw, n_raw, 0,29 "audit found no missing, duplicate or non-positive values; nothing removed"])3031X = np.column_stack([np.ones(N), x])32XtX_inv = np.linalg.inv(X.T @ X)33xc = x - x.mean(); Sxx = xc @ xc34wts = xc / Sxx # beta_hat - beta = sum w_i sigma_i eps_i35sd_beta_exact = float(np.sqrt(np.sum(wts**2 * sig**2)))36tcrit = stats.t.ppf(0.5 + LEVEL / 2, N - 2)3738rng = np.random.default_rng(SEED)39eps = rng.standard_normal((R, N)) # one row per replication40Y = ALPHA + BETA * x[None, :] + sig[None, :] * eps41A = X @ XtX_inv # (N,2): beta_hat = Y @ A42B = Y @ A43resid = Y - B @ X.T44s2 = (resid**2).sum(1) / (N - 2)45se_conv = np.sqrt(s2 * XtX_inv[1, 1])46# HC0 slope variance: sum_i a_i1^2 e_i^2 (= [(X'X)^-1 X' diag(e^2) X (X'X)^-1]_{11})47hc0 = (resid**2) @ (A[:, 1]**2)48se_hc1 = np.sqrt(hc0 * N / (N - 2))49bhat = B[:, 1]50cov_conv = float(np.mean(np.abs(bhat - BETA) <= tcrit * se_conv))51cov_hc1 = float(np.mean(np.abs(bhat - BETA) <= tcrit * se_hc1))52mcse_nominal = float(np.sqrt(LEVEL * (1 - LEVEL) / R))53undercovers = bool(cov_conv < LEVEL - 3 * mcse_nominal)5455res = {56 "task_id": "F2-conf-03",57 "coverage_conventional": cov_conv,58 "coverage_hc1": cov_hc1,59 "sd_beta_exact": sd_beta_exact,60 "sd_beta_mc": float(np.std(bhat, ddof=1)),61 "mean_se_conventional": float(se_conv.mean()),62 "conventional_ci_undercovers": undercovers,63 "n_replications": R,64 "n_obs_per_replication": N,65 "feasibility": "complete",66 "feasibility_reason": None,67 "data_issues": issues,68 "n_input_rows": {"design.csv": int(n_raw)},69 "random_seed": SEED,70 # supplementary keys71 "t_critical_28": float(tcrit),72 "mcse_coverage_at_nominal": mcse_nominal,73 "undercoverage_threshold": float(LEVEL - 3 * mcse_nominal),74 "mcse_coverage_conventional": float(np.sqrt(cov_conv * (1 - cov_conv) / R)),75 "mcse_coverage_hc1": float(np.sqrt(cov_hc1 * (1 - cov_hc1) / R)),76 "sd_beta_homoskedastic_formula": float(np.sqrt(np.mean(sig**2) / Sxx)),77 "mean_se_hc1": float(se_hc1.mean()),78 "mean_beta_hat": float(bhat.mean()),79 "mcse_mean_beta_hat": float(np.std(bhat, ddof=1) / np.sqrt(R)),80 "mcse_mean_se_conventional": float(np.std(se_conv, ddof=1) / np.sqrt(R)),81}82with open("out/results.json", "w") as f:83 json.dump(res, f, indent=2)84print(json.dumps(res, indent=2))85