1"""F2-conf-03: Monte Carlo coverage of OLS slope CIs under heteroskedasticity (fixed design).2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/design.csv; writes out/results.json deterministically."""4import json5from pathlib import Path6import numpy as np7import pandas as pd8from scipy import stats910ALPHA, BETA, N, R, LEVEL = 1.0, 0.5, 30, 2000, 0.9511SEED = 20261030 # fixed seed for numpy.random.default_rng (PCG64)1213d = pd.read_csv("data/design.csv")14n_rows = int(len(d))15issues = []16assert n_rows == N and d[["x", "sigma"]].notna().all().all() and (d["sigma"] > 0).all()17x = d["x"].to_numpy(float)18sig = d["sigma"].to_numpy(float)1920xc = x - x.mean()21Sxx = float(xc @ xc)22w = xc / Sxx # beta_hat = beta + sum(w_i sigma_i eps_i)23sd_beta_exact = float(np.sqrt(np.sum(w**2 * sig**2)))24sd_beta_conv_formula = float(np.sqrt(np.mean(sig**2) / Sxx)) # naive homoskedastic SD (diagnostic only)2526tcrit = float(stats.t.ppf(0.5 + LEVEL / 2, N - 2))2728rng = np.random.default_rng(SEED)29eps = rng.standard_normal((R, N))30Y = ALPHA + BETA * x[None, :] + sig[None, :] * eps31X = np.column_stack([np.ones(N), x])32XtX_inv = np.linalg.inv(X.T @ X)33B = Y @ X @ XtX_inv # (R,2) OLS coefficients (XtX_inv symmetric)34bhat = B[:, 1]35E = Y - B @ X.T # residuals36s2 = np.sum(E**2, axis=1) / (N - 2)37se_conv = np.sqrt(s2 * XtX_inv[1, 1])38# HC0 sandwich for slope: (X'X)^-1 X' diag(e^2) X (X'X)^-1, element [1,1]; HC1 = HC0 * n/(n-2)39A = X @ XtX_inv # row i: a_i' ; slope weight = A[:,1] (= w)40hc0 = (E**2) @ (A[:, 1] ** 2)41se_hc1 = np.sqrt(hc0 * N / (N - 2))4243cov_conv = float(np.mean(np.abs(bhat - BETA) <= tcrit * se_conv))44cov_hc1 = float(np.mean(np.abs(bhat - BETA) <= tcrit * se_hc1))45mcse_nominal = float(np.sqrt(LEVEL * (1 - LEVEL) / R))46mcse_conv = float(np.sqrt(cov_conv * (1 - cov_conv) / R))47mcse_hc1 = float(np.sqrt(cov_hc1 * (1 - cov_hc1) / R))48undercovers = bool(cov_conv < LEVEL - 3 * mcse_nominal)4950res = {51 "task_id": "F2-conf-03",52 "coverage_conventional": cov_conv,53 "coverage_hc1": cov_hc1,54 "sd_beta_exact": sd_beta_exact,55 "sd_beta_mc": float(np.std(bhat, ddof=1)),56 "mean_se_conventional": float(np.mean(se_conv)),57 "conventional_ci_undercovers": undercovers,58 "n_replications": R,59 "n_obs_per_replication": N,60 "feasibility": "complete",61 "feasibility_reason": None,62 "data_issues": issues,63 "n_input_rows": {"design.csv": n_rows},64 "random_seed": SEED,65 # supplementary (not required) diagnostics66 "t_critical_28": tcrit,67 "mcse_binomial_nominal": mcse_nominal,68 "mcse_coverage_conventional": mcse_conv,69 "mcse_coverage_hc1": mcse_hc1,70 "undercoverage_threshold": LEVEL - 3 * mcse_nominal,71 "mean_se_hc1": float(np.mean(se_hc1)),72 "sd_beta_homoskedastic_formula": sd_beta_conv_formula,73 "mean_beta_hat_mc": float(np.mean(bhat)),74}75Path("out").mkdir(exist_ok=True)76Path("out/results.json").write_text(json.dumps(res, indent=2))77for k, v in res.items():78 print(f"{k}: {v}")79