# F2-conf-03: Monte Carlo coverage of OLS confidence intervals with heteroskedastic errors (n = 30)

## Question
The model has a fixed design: y_i = 1 + 0.5 x_i + σ_i ε_i, with ε_i iid N(0,1), and x_i and σ_i taken from `data/design.csv` (n = 30). The question is whether the nominal 95% interval β̂ ± t_{0.975}(28)·se for the OLS slope reaches 95% coverage. Two standard errors are compared: (i) the conventional OLS standard error and (ii) the HC1 robust standard error. Two further quantities are reported: the exact SD of β̂ and E[se_conv].

## Data and exclusions
The only input is `data/design.csv`. It has 30 rows (`n_input_rows`), SHA-256 b1a17a52…74f0, and matches the bundle manifest. It has no missing values, no duplicates and all σ_i > 0. No rows were excluded (`out/exclusions.csv`), and `data_issues` is empty. σ_i rises toward both ends of the x range, so the largest error variances sit at the high-leverage points. This is the pattern under which the conventional SE is biased downward.

## Method
- **Exact SD.** β̂ − β = Σ w_i σ_i ε_i, where w_i = (x_i − x̄)/S_xx. So SD(β̂) = sqrt(Σ w_i² σ_i²), computed analytically (`sd_beta_exact`).
- **Monte Carlo.** R = 2000 replications (`n_replications`) with numpy `default_rng` seed 20261030 (`random_seed`). Only ε is redrawn. In each replication:
  - fit OLS;
  - compute se_conv = sqrt(s²·[(X'X)⁻¹]₂₂), with s² = RSS/28;
  - compute se_HC1 = sqrt((n/(n−2))·Σ a_i² e_i²), where a = row 2 of (X'X)⁻¹X';
  - check whether β = 0.5 lies inside β̂ ± 2.0484·se (`t_crit_28`).
- **Undercoverage rule.** The conventional interval counts as undercovering if its coverage is below 0.95 − 3·sqrt(0.95·0.05/2000) = 0.9354 (`undercover_threshold`). The binomial SE here is evaluated at the nominal p (`mcse_binomial_at_nominal`).
- **Validation.** One replication was cross-checked against statsmodels OLS/HC1 (exact match).

## Results
- **Conventional interval coverage:** 0.844 (`coverage_conventional`). This is the naïve baseline.
- **HC1 interval coverage:** 0.902 (`coverage_hc1`).
- **Exact SD(β̂):** 0.2680 (`sd_beta_exact`). The Monte Carlo SD is 0.2696 (`sd_beta_mc`), which agrees.
- **Mean se_conv:** 0.1903 (`mean_se_conventional`). That is about 71% of the true SD. The formula that assumes constant variance gives 0.1962 (`sd_beta_conventional_formula_homoskedastic`).
- **Mean se_HC1:** 0.2408 (`mean_se_hc1`).
- **Conclusion:** the conventional 95% interval does **not** reach its nominal coverage (`conventional_ci_undercovers` = true). It undercovers by about 10.6 percentage points. HC1 does better but still falls short of 95% at n = 30.

## Uncertainty
- **Monte Carlo SE of each coverage estimate** (binomial, evaluated at p̂):
  - conventional: 0.0081 (`mcse_coverage_conventional`), giving an approximate 95% MC interval of [0.828, 0.860];
  - HC1: 0.0066 (`mcse_coverage_hc1`), giving approximately [0.889, 0.915].
- **At nominal p = 0.95,** the binomial SE is 0.0049 (`mcse_binomial_at_nominal`).
- **Size of the shortfall:** both coverages lie more than 10 nominal SEs below 0.95. The conclusion does not depend on which SE is used.
- **Other Monte Carlo checks:**
  - the MC SD of β̂ is within one SE (≈0.004) of the exact value;
  - mean β̂ = 0.4962 (`mean_beta_hat_mc`), which is within one MC SE (≈0.006) of 0.5.

## Limitations and open issues
- The results hold conditionally on this particular fixed design and on normal errors. They are not general statements about OLS.
- With R = 2000, coverage is known only to about ±0.016 (95% MC interval). The conclusions above are well clear of that margin.
- The task specifies the HC1 interval with t(28) critical values. Small-sample alternatives (HC2/HC3, or Satterthwaite/Bell–McCaffrey degrees of freedom) would probably come closer to 95%, but they were not evaluated because the task did not ask for them.
- No part of the task was infeasible.