1"""F2-conf-03: Monte Carlo coverage of OLS CIs under heteroskedasticity (fixed design)."""2import json, os3import numpy as np, pandas as pd4from scipy import stats56ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))7df = pd.read_csv(os.path.join(ROOT, "data", "design.csv"))8issues = []9n_rows = len(df)10if df[["x", "sigma"]].isna().any().any(): issues.append({"issue": "missing values", "handling": "none expected; would abort"})11if (df["sigma"] <= 0).any(): issues.append({"issue": "non-positive sigma", "handling": "flagged"})12if not np.array_equal(df["i"].values, np.arange(1, 31)): issues.append({"issue": "index not 1..30", "handling": "used rows as given"})1314ALPHA, BETA, R, SEED = 1.0, 0.5, 2000, 2026103015x = df["x"].to_numpy(float); sig = df["sigma"].to_numpy(float); n = len(x)16X = np.column_stack([np.ones(n), x])17XtXi = np.linalg.inv(X.T @ X)18A = XtXi @ X.T # 2 x n, beta_hat = A y19tcrit = stats.t.ppf(0.975, n - 2)2021# (c) exact SD: Var(b1) = sum_i a_i^2 sigma_i^222a1 = A[1]23sd_exact = float(np.sqrt(np.sum(a1**2 * sig**2)))2425rng = np.random.default_rng(SEED)26eps = rng.standard_normal((R, n))27Y = ALPHA + BETA * x + sig * eps # R x n28B = Y @ A.T # R x 229E = Y - B @ X.T # residuals30s2 = (E**2).sum(1) / (n - 2)31se_conv = np.sqrt(s2 * XtXi[1, 1])32hc0 = (E**2) @ (a1**2) # sandwich var for slope33se_hc1 = np.sqrt(hc0 * n / (n - 2))34b1 = B[:, 1]35cov_conv = float(np.mean(np.abs(b1 - BETA) <= tcrit * se_conv))36cov_hc1 = float(np.mean(np.abs(b1 - BETA) <= tcrit * se_hc1))37mcse_conv = float(np.sqrt(0.95 * 0.05 / R)) # binomial SE at nominal p38under = bool(cov_conv < 0.95 - 3 * mcse_conv)3940res = {41 "task_id": "F2-conf-03",42 "coverage_conventional": cov_conv,43 "coverage_hc1": cov_hc1,44 "sd_beta_exact": sd_exact,45 "sd_beta_mc": float(np.std(b1, ddof=1)),46 "mean_se_conventional": float(se_conv.mean()),47 "conventional_ci_undercovers": under,48 "n_replications": R,49 "n_obs_per_replication": n,50 "feasibility": "complete",51 "feasibility_reason": None,52 "data_issues": issues,53 "n_input_rows": {"design.csv": n_rows},54 "random_seed": SEED,55 # supplementary56 "mcse_binomial_at_0.95": mcse_conv,57 "mcse_coverage_conventional_at_estimate": float(np.sqrt(cov_conv * (1 - cov_conv) / R)),58 "mcse_coverage_hc1_at_estimate": float(np.sqrt(cov_hc1 * (1 - cov_hc1) / R)),59 "mean_se_hc1": float(se_hc1.mean()),60 "expected_se_conv_sq_exact": float(XtXi[1, 1] * np.sum(sig**2 * (1 - np.diag(X @ A))) / (n - 2)),61 "t_crit": float(tcrit),62}63with open(os.path.join(ROOT, "out", "results.json"), "w") as f:64 json.dump(res, f, indent=2)65print(json.dumps(res, indent=2))66