test(se-audit): G2 hetero + cluster machine-precision fixest locks on unbalanced heteroskedastic DGP

CHANGELOG.md

@@ -49,6 +49,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   `TestHonestFLCIParityR`); the M=0 result and all existing behaviour are unchanged.
 
 ### Testing
+- **fixest hetero + cluster SE machine-precision locks on an unbalanced, heteroskedastic
+  DGP (SE-audit G2 completion).** The committed `fixest_did_twfe_golden.json` gains two
+  appended scenarios (error sd varying by arm/period, ~15% rows dropped; the original
+  balanced scenarios' RNG draws precede them and reproduce value-identically): on the
+  plain-OLS DiD path, `hetero` (HC1) no longer collapses to iid and is locked against
+  `fixest` at machine precision — and the cluster-robust CR1 SE turns out to match fixest
+  **exactly** on plain OLS (balanced and unbalanced), so the former ~0.5% DiD cluster
+  band-pin is tightened to a machine-precision lock: the documented ~0.25% fixest-CR1
+  DOF-convention deviation is an absorbed-FE (within-transform) phenomenon only. The TWFE
+  cluster band-pin is retained and re-scoped to that documented non-nested-FE ssc
+  deviation (~0.3% unbalanced); TWFE `hetero` has no public unclustered surface
+  (auto-cluster-at-unit convention), so its scenario locks iid — which also pins the D4
+  full-K rescale on an UNBALANCED panel for the first time.
 - **`CallawaySantAnna` ipw R-parity yardsticks folded into the golden fixture + no-covariate
   ipw structural-parity decision recorded.** `csdid_golden_values.json` regenerated (R 4.5.2,
   did 2.5.1, DRDID 1.3.0): all pre-existing data and result blocks reproduced byte-identically;

TODO.md

@@ -62,7 +62,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 
 | Issue | Location | Origin | Effort | Priority |
 |-------|----------|--------|--------|----------|
-| SE-audit CI-lock — remaining "fiddly bits" after the second coverage batch landed (that batch pinned C2 `dof_hc2_bm`/`dof_per_coef` via CI-inversion, C3 LOO `df`, C4 estimatr HC1/CR1 intercept SE, C5 Yatchew `p`/`sigma2_lin`/`sigma2_diff`, and the G2 fixest cluster-SE band). Still deferred, each needing a golden regeneration, new computation, or a documented-deviation call: **(a) G2 machine-precision hetero/cluster SE** — needs an unbalanced/heteroskedastic-DGP regen of `fixest_did_twfe_golden.json` (hetero collapses to iid on the current balanced design; the cluster SE is the documented ~0.25% fixest-CR1 DOF-convention deviation, currently band-pinned only); **(b) PlaceboTests `boundary_gap`** — a permutation randomization-inference margin NOT computed anywhere in code (a new feature + result field, not a coverage lock); **(c) StackedDiD intercept SEs** (`se_cr1/cr2_intercept`, C1) — MEASURED to diverge ~0.3% from R: a nuisance-parameter reference-cell/parameterization gap, NOT machine-precision lockable (the event-study interaction SEs already match ~2e-13; surfacing it would add an unasserted, R-divergent public field); **(d) estimatr `classical` intercept SE** — same documented `O(1/n)` projection/DOF deviation as the slope (reference-only, excluded from parity). Plus the tolerance-tightenings C6-C8. | `benchmarks/R/generate_fixest_did_twfe_golden.R`, `tests/test_fixest_did_twfe_parity.py`, `tests/test_methodology_stacked_did.py`, `tests/test_methodology_placebo.py` | SE-audit | Mid | Low |
+| SE-audit CI-lock — remaining "fiddly bits" after the second coverage batch landed (that batch pinned C2 `dof_hc2_bm`/`dof_per_coef` via CI-inversion, C3 LOO `df`, C4 estimatr HC1/CR1 intercept SE, C5 Yatchew `p`/`sigma2_lin`/`sigma2_diff`, and the G2 fixest cluster-SE band). Still deferred, each needing a golden regeneration, new computation, or a documented-deviation call — **(a) G2 hetero/cluster is DONE (2026-07-07)**: the unbalanced/heteroskedastic-DGP regen landed, the DiD-path hetero AND cluster CR1 SEs are machine-precision-locked (the plain-OLS CR1 matches fixest exactly — the DOF-convention deviation is absorbed-FE-only), and the TWFE cluster band-pin is retained for the documented non-nested-FE ssc deviation (tracked under "Needs external reference"); TWFE has no public unclustered-hetero surface (auto-cluster convention). Remaining: **(b) PlaceboTests `boundary_gap`** — a permutation randomization-inference margin NOT computed anywhere in code (a new feature + result field, not a coverage lock); **(c) StackedDiD intercept SEs** (`se_cr1/cr2_intercept`, C1) — MEASURED to diverge ~0.3% from R: a nuisance-parameter reference-cell/parameterization gap, NOT machine-precision lockable (the event-study interaction SEs already match ~2e-13; surfacing it would add an unasserted, R-divergent public field); **(d) estimatr `classical` intercept SE** — same documented `O(1/n)` projection/DOF deviation as the slope (reference-only, excluded from parity). Plus the tolerance-tightenings C6-C8. | `benchmarks/R/generate_fixest_did_twfe_golden.R`, `tests/test_fixest_did_twfe_parity.py`, `tests/test_methodology_stacked_did.py`, `tests/test_methodology_placebo.py` | SE-audit | Mid | Low |
 | Render `docs/methodology/REPORTING.md` and `REGISTRY.md` as in-site Sphinx pages so cross-refs can use `:doc:` instead of off-site `blob/main` URLs (stable-docs readers can otherwise land on a different revision than their package version). Two paths: (a) add `myst-parser` to `conf.py` + docs extras and link with `:doc:`, or (b) convert both to `.rst`. **Note:** REGISTRY.md is ~4.5k lines of LaTeX-heavy markdown — high risk under the `-W` (warnings-as-errors) Sphinx build; budget multiple rounds. | `docs/conf.py`, `docs/api/business_report.rst`, `docs/api/diagnostic_report.rst`, tutorials 18 & 19 | follow-up | Mid | Low |
 | `ImputationDiD` covariate-path variance lacks a dedicated parity anchor — only the no-covariate staggered panel is R-parity'd, though the covariate path shares the same validated projection code. Add a small dense-design **hand-calc** for the covariate projection (no external tooling), or a covariate (time-varying X) R `didimputation` golden asserting overall/ES SE parity (the golden variant needs local R). | `tests/test_methodology_imputation.py`, `benchmarks/R/generate_didimputation_golden.R` | imputation-validation | Mid | Low |
 | Add true half-sample BRR replicate-weight regressions per estimator family (current tests use Fay-like 0.5/1.5 perturbations; `test_survey_phase6.py` covers true BRR at the helper level). | `tests/test_replicate_weight_expansion.py` | #253 | Mid | Low |

benchmarks/R/generate_fixest_did_twfe_golden.R

@@ -7,12 +7,16 @@
 # panels and fixest's feols() ATT + SE so tests can assert machine-precision SE
 # parity WITHOUT R at test time.
 #
-# Scope (this golden): the classical / iid SE, which Python matches to machine
-# precision on both the 2x2 DiD path and the within-transform TWFE path (the
-# latter also locks the SE-audit D4 full-K rescale). The cluster-robust ATT is
-# stored too; its SE carries the documented CR1 small-sample DOF-convention
-# difference vs fixest and is left to a follow-up. (`hetero`/HC1 collapses to iid
-# on these balanced 2-group designs, so it is not a distinct target here.)
+# Scope (this golden): scenarios 1-2 are the original balanced designs — the
+# classical / iid SE locks (the TWFE one also pins the SE-audit D4 full-K
+# rescale) plus the cluster blocks. Scenarios 3-4 (G2 completion, 2026-07) are
+# heteroskedastic + unbalanced so `hetero` (HC1) is a distinct target: the
+# plain-OLS DiD path locks hetero AND cluster CR1 at machine precision (the
+# CR1 DOF-convention difference vs fixest is absorbed-FE-only); the TWFE
+# cluster SE stays band-pinned for that documented non-nested-FE ssc
+# deviation, and TWFE hetero has no public unclustered Python surface
+# (auto-cluster-at-unit convention), so scenario 4 locks iid on an
+# UNBALANCED panel.
 #
 # Regenerate:  Rscript benchmarks/R/generate_fixest_did_twfe_golden.R
 # Output:      benchmarks/data/fixest_did_twfe_golden.json
@@ -97,6 +101,84 @@ twfe_golden <- list(
   cluster_unit = fit_att(twfe_m, ~unit)
 )
 
+# ---------------------------------------------------------------------------
+# Scenario 3: heteroskedastic + unbalanced 2x2 DiD (SE-audit G2 hetero lock).
+# Error sd depends on treatment arm and period (so HC1 'hetero' does NOT
+# collapse to iid) and ~15% of rows are dropped deterministically-by-draw
+# (unbalanced groups). Appended AFTER scenarios 1-2 so their RNG draws (and
+# the committed scenario 1-2 golden values) are unchanged on regeneration.
+# ---------------------------------------------------------------------------
+n_units_h <- 120
+did_h_rows <- list()
+i <- 1
+for (unit in 0:(n_units_h - 1)) {
+  is_treated <- as.integer(unit < 45)   # unequal arms: 45 treated / 75 control
+  for (period in c(0, 1)) {
+    sd_it <- 0.5 + 1.5 * is_treated + 0.8 * period   # heteroskedastic
+    y <- 10.0 + period * 2.0
+    if (is_treated == 1 && period == 1) y <- y + 3.0
+    y <- y + rnorm(1, 0, sd_it)
+    keep <- runif(1) > 0.15                          # unbalanced: drop ~15%
+    if (keep) {
+      did_h_rows[[i]] <- data.frame(unit = unit, outcome = y, treated = is_treated, post = period)
+      i <- i + 1
+    }
+  }
+}
+did_h <- do.call(rbind, did_h_rows)
+did_h_m <- feols(outcome ~ treated * post, data = did_h)
+
+did_hetero_golden <- list(
+  data = list(unit = did_h$unit, outcome = did_h$outcome, treated = did_h$treated, post = did_h$post),
+  n_obs = unbox(nrow(did_h)),
+  iid = fit_att(did_h_m, "iid"),
+  hetero = fit_att(did_h_m, "hetero"),
+  cluster_unit = fit_att(did_h_m, ~unit)
+)
+
+# ---------------------------------------------------------------------------
+# Scenario 4: heteroskedastic + unbalanced TWFE. The fixest `hetero` block is
+# stored for reference only — Python's TwoWayFixedEffects auto-clusters at
+# unit on hc1 (no public unclustered-hetero surface), so the public locks are
+# the unbalanced iid/full-K rescale plus the clustered ATT (exact) and SE
+# (band, documented non-nested-FE ssc deviation).
+# ---------------------------------------------------------------------------
+n_units_th <- 40
+n_periods_th <- 5
+twfe_h_rows <- list()
+i <- 1
+for (unit in 0:(n_units_th - 1)) {
+  is_treated <- as.integer(unit < 15)   # unequal arms
+  unit_effect <- unit * 0.2
+  for (period in 0:(n_periods_th - 1)) {
+    post <- as.integer(period >= 3)
+    sd_it <- 0.4 + 1.2 * is_treated + 0.3 * post
+    y <- 5.0 + unit_effect + period * 1.5
+    if (is_treated == 1 && post == 1) y <- y + 2.5
+    y <- y + rnorm(1, 0, sd_it)
+    keep <- runif(1) > 0.12
+    if (keep) {
+      twfe_h_rows[[i]] <- data.frame(
+        unit = unit, period = period, outcome = y, treated = is_treated, post = post
+      )
+      i <- i + 1
+    }
+  }
+}
+twfe_h <- do.call(rbind, twfe_h_rows)
+twfe_h_m <- feols(outcome ~ treated:post | unit + post, data = twfe_h)
+
+twfe_hetero_golden <- list(
+  data = list(
+    unit = twfe_h$unit, period = twfe_h$period, outcome = twfe_h$outcome,
+    treated = twfe_h$treated, post = twfe_h$post
+  ),
+  n_obs = unbox(nrow(twfe_h)),
+  iid = fit_att(twfe_h_m, "iid"),
+  hetero = fit_att(twfe_h_m, "hetero"),
+  cluster_unit = fit_att(twfe_h_m, ~unit)
+)
+
 # ---------------------------------------------------------------------------
 golden <- list(
   meta = list(
@@ -109,7 +191,9 @@ golden <- list(
     ))
   ),
   did = did_golden,
-  twfe = twfe_golden
+  twfe = twfe_golden,
+  did_hetero = did_hetero_golden,
+  twfe_hetero = twfe_hetero_golden
 )
 
 out <- "benchmarks/data/fixest_did_twfe_golden.json"