test(se-audit): CI-lock flagship + previously-unasserted SE parity (G2 + C1-C3/C5)

CHANGELOG.md

@@ -7,6 +7,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Testing
+- **CI-locked standard-error parity for flagship and previously-unasserted paths (SE-audit
+  coverage batch).** These surfaces computed SEs matching R but had no CI assertion pinning them
+  (the latent-risk pattern that once hid the CallawaySantAnna reg-method gap):
+  - **New `fixest` DiD/TWFE golden** (`benchmarks/data/fixest_did_twfe_golden.json`, generated by
+    `benchmarks/R/generate_fixest_did_twfe_golden.R`): the flagship 2×2 DiD and TWFE classical/iid
+    ATT **and SE** now match `fixest::feols` to ~1e-16 in CI without R at test time (the prior
+    live-`Rscript` tests only asserted `att` at rtol=1e-3 and skipped in CI). The TWFE assertion
+    also locks the D4 full-K within-transform rescale.
+  - **CR2 backbone**: `TwoWayFixedEffects(vcov_type="hc2_bm")`'s default CR2-clustered-at-unit SE
+    and the one-way HC2-BM SE (previously only their Satterthwaite DOF was asserted); the
+    `MultiPeriodDiD` average-effect CR2 SE (previously only checked finite).
+  - **StackedDiD** event-study coefficient point estimates; **PlaceboTests** leave-one-out
+    `t_stat`/`p_value`; **DIDHAD** QUG order-statistic test (`t_stat`/`p_value`, previously
+    unasserted against the golden). All at ~1e-10.
+
 ### Added
 - **`ContinuousDiD` discrete-treatment saturated regression** (`treatment_type="discrete"`) for
   multi-valued / discrete dose (CGBS 2024 Eq. 4.1). Each distinct dose level gets its own effect

TODO.md

@@ -66,6 +66,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 
 | Issue | Location | Origin | Effort | Priority |
 |-------|----------|--------|--------|----------|
+| SE-audit CI-lock follow-up ("fiddly bits" deferred from the G2+C1-C3/C5 coverage batch). The `fixest_did_twfe_golden.json` stores `cluster_unit` ATT+SE for DiD and TWFE but only the ATT is asserted (the cluster SE carries the documented CR1 DOF-convention difference vs fixest; the hetero SE collapses to iid on these balanced 2-group designs). Also deferred, each needing plumbing or convention work: surface the 2SLS intercept SE for estimatr `se_intercept` parity (C4); surface StackedDiD intercept SEs (`se_cr1/cr2_intercept`, C1); add a `df` field + `boundary_gap` computation for PlaceboTests (C3); assert Yatchew `p`/`sigma2_lin/diff` through the N/(N-1) convention shift + event-map (C5); assert `dof_hc2_bm`/`dof_per_coef` via CI-inversion (C2). Plus the tolerance-tightenings C6-C8. | `tests/test_fixest_did_twfe_parity.py`, `tests/test_estimatr_iv_robust_parity.py`, `tests/test_methodology_stacked_did.py`, `tests/test_methodology_placebo.py`, `tests/test_did_had_parity.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

@@ -0,0 +1,117 @@
+# Generate a committed fixest golden for the flagship 2x2 DiD and TWFE
+# standard-error parity paths (SE-audit item G2).
+#
+# The core DiD / TWFE fixest SE parity was previously only checked by
+# skip-guarded live-Rscript tests that never run in CI -- and those tests only
+# assert `att` (rtol=1e-3), never the SE. This golden materializes deterministic
+# 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.)
+#
+# Regenerate:  Rscript benchmarks/R/generate_fixest_did_twfe_golden.R
+# Output:      benchmarks/data/fixest_did_twfe_golden.json
+# Environment: R 4.5.2, fixest 0.14.2.
+
+suppressMessages(library(fixest))
+suppressMessages(library(jsonlite))
+
+set.seed(12345)
+
+fit_att <- function(model, vcov_spec) {
+  ct <- coeftable(model, vcov = vcov_spec)
+  idx <- which(rownames(ct) == "treated:post")
+  ci <- confint(model, vcov = vcov_spec)[idx, ]
+  list(
+    att = unbox(ct[idx, "Estimate"]),
+    se = unbox(ct[idx, "Std. Error"]),
+    t_stat = unbox(ct[idx, "t value"]),
+    p_value = unbox(ct[idx, "Pr(>|t|)"]),
+    ci_lower = unbox(ci[[1]]),
+    ci_upper = unbox(ci[[2]])
+  )
+}
+
+# ---------------------------------------------------------------------------
+# Scenario 1: 2x2 DiD (200 units x 2 periods), outcome ~ treated * post
+# ---------------------------------------------------------------------------
+n_units_did <- 200
+did_rows <- list()
+i <- 1
+for (unit in 0:(n_units_did - 1)) {
+  is_treated <- as.integer(unit < n_units_did %/% 2)
+  for (period in c(0, 1)) {
+    y <- 10.0 + period * 2.0
+    if (is_treated == 1 && period == 1) y <- y + 3.0
+    y <- y + rnorm(1, 0, 1)
+    did_rows[[i]] <- data.frame(unit = unit, outcome = y, treated = is_treated, post = period)
+    i <- i + 1
+  }
+}
+did <- do.call(rbind, did_rows)
+did_m <- feols(outcome ~ treated * post, data = did)
+
+did_golden <- list(
+  data = list(unit = did$unit, outcome = did$outcome, treated = did$treated, post = did$post),
+  n_obs = unbox(nrow(did)),
+  iid = fit_att(did_m, "iid"),
+  cluster_unit = fit_att(did_m, ~unit)
+)
+
+# ---------------------------------------------------------------------------
+# Scenario 2: TWFE (50 units x 4 periods), outcome ~ treated:post | unit + post
+# ---------------------------------------------------------------------------
+n_units_twfe <- 50
+n_periods_twfe <- 4
+twfe_rows <- list()
+i <- 1
+for (unit in 0:(n_units_twfe - 1)) {
+  is_treated <- as.integer(unit < n_units_twfe %/% 2)
+  unit_effect <- unit * 0.2
+  for (period in 0:(n_periods_twfe - 1)) {
+    post <- as.integer(period >= n_periods_twfe %/% 2)
+    y <- 5.0 + unit_effect + period * 1.5
+    if (is_treated == 1 && post == 1) y <- y + 2.5
+    y <- y + rnorm(1, 0, 1)
+    twfe_rows[[i]] <- data.frame(
+      unit = unit, period = period, outcome = y, treated = is_treated, post = post
+    )
+    i <- i + 1
+  }
+}
+twfe <- do.call(rbind, twfe_rows)
+twfe_m <- feols(outcome ~ treated:post | unit + post, data = twfe)
+
+twfe_golden <- list(
+  data = list(
+    unit = twfe$unit, period = twfe$period, outcome = twfe$outcome,
+    treated = twfe$treated, post = twfe$post
+  ),
+  n_obs = unbox(nrow(twfe)),
+  iid = fit_att(twfe_m, "iid"),
+  cluster_unit = fit_att(twfe_m, ~unit)
+)
+
+# ---------------------------------------------------------------------------
+golden <- list(
+  meta = list(
+    generator = unbox("benchmarks/R/generate_fixest_did_twfe_golden.R"),
+    r_version = unbox(paste(R.version$major, R.version$minor, sep = ".")),
+    fixest_version = unbox(as.character(packageVersion("fixest"))),
+    description = unbox(paste(
+      "fixest feols() ATT + SE golden for the flagship 2x2 DiD",
+      "(outcome ~ treated*post) and TWFE (| unit + post) paths, SE-audit G2."
+    ))
+  ),
+  did = did_golden,
+  twfe = twfe_golden
+)
+
+out <- "benchmarks/data/fixest_did_twfe_golden.json"
+writeLines(toJSON(golden, pretty = TRUE, digits = 16, auto_unbox = FALSE), out)
+cat("Wrote", out, "\n")