1"""F2-conf-03: Monte Carlo coverage of OLS slope CIs under fixed heteroskedastic design.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 stats8import statsmodels.api as sm910SEED = 2026103011R = 200012ALPHA, BETA = 1.0, 0.513LEVEL = 0.951415d = pd.read_csv("data/design.csv")16n_rows = int(len(d))17data_issues = []18assert n_rows == 30 and d[["x", "sigma"]].notna().all().all() and (d.sigma > 0).all()19x = d["x"].to_numpy(float)20sig = d["sigma"].to_numpy(float)21n = len(x)22df = n - 223tcrit = stats.t.ppf(0.5 + LEVEL / 2, df)2425# (c) exact SD of betahat: betahat - beta = sum w_i sigma_i eps_i, w_i = (x_i - xbar)/Sxx26xc = x - x.mean()27Sxx = np.sum(xc ** 2)28w = xc / Sxx29sd_exact = float(np.sqrt(np.sum(w ** 2 * sig ** 2)))3031# Monte Carlo (only eps redrawn; x, sigma fixed)32rng = np.random.default_rng(SEED)33eps = rng.standard_normal((R, n))34Y = ALPHA + BETA * x[None, :] + sig[None, :] * eps35X = np.column_stack([np.ones(n), x])36XtX_inv = np.linalg.inv(X.T @ X)37P = XtX_inv @ X.T # 2 x n38B = Y @ P.T # R x 2 coefficients39resid = Y - B @ X.T # R x n40s2 = np.sum(resid ** 2, axis=1) / df41se_conv = np.sqrt(s2 * XtX_inv[1, 1])42a = P[1] # slope row of (X'X)^-1 X'43var_hc0 = np.sum((a[None, :] ** 2) * resid ** 2, axis=1)44se_hc1 = np.sqrt(var_hc0 * n / df)45bhat = B[:, 1]4647cov_conv = float(np.mean(np.abs(bhat - BETA) <= tcrit * se_conv))48cov_hc1 = float(np.mean(np.abs(bhat - BETA) <= tcrit * se_hc1))49mcse_nominal = float(np.sqrt(LEVEL * (1 - LEVEL) / R))50undercovers = bool(cov_conv < LEVEL - 3 * mcse_nominal)5152# Cross-check against statsmodels on replication 053m0 = sm.OLS(Y[0], X).fit()54assert np.isclose(m0.bse[1], se_conv[0])55assert np.isclose(m0.get_robustcov_results("HC1").bse[1], se_hc1[0])56assert np.isclose(np.sqrt(np.sum(a ** 2 * sig ** 2)), sd_exact)5758H = np.diag(X @ XtX_inv @ X.T)59results = {60 "task_id": "F2-conf-03",61 "coverage_conventional": cov_conv,62 "coverage_hc1": cov_hc1,63 "sd_beta_exact": sd_exact,64 "sd_beta_mc": float(np.std(bhat, ddof=1)),65 "mean_se_conventional": float(np.mean(se_conv)),66 "conventional_ci_undercovers": undercovers,67 "n_replications": R,68 "n_obs_per_replication": n,69 "feasibility": "complete",70 "feasibility_reason": None,71 "data_issues": data_issues,72 "n_input_rows": {"design.csv": n_rows},73 "random_seed": SEED,74 # supplementary diagnostics75 "t_critical_28": float(tcrit),76 "mcse_coverage_conventional": float(np.sqrt(cov_conv * (1 - cov_conv) / R)),77 "mcse_coverage_hc1": float(np.sqrt(cov_hc1 * (1 - cov_hc1) / R)),78 "mcse_at_nominal_095": mcse_nominal,79 "undercoverage_threshold": float(LEVEL - 3 * mcse_nominal),80 "mean_se_hc1": float(np.mean(se_hc1)),81 "mean_beta_hat": float(np.mean(bhat)),82 "sd_se_conventional_mc": float(np.std(se_conv, ddof=1)),83 "sqrt_expected_se_conv_sq_exact": float(np.sqrt(XtX_inv[1, 1] * np.sum(sig ** 2 * (1 - H)) / df)),84}85os.makedirs("out", exist_ok=True)86with open("out/results.json", "w") as f:87 json.dump(results, f, indent=2)88print(json.dumps(results, indent=2))89