feat(3.9): converge clustered CR1 on K_reference - D1+D2 fix (M-126)

CHANGELOG.md

@@ -8,6 +8,42 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Fixed
+- **Clustered CR1 SEs converge on `K_reference` — the D1+D2 fix of the 3.9
+  variance-consolidation program** ([M-126]). The clustered `hc1` CR1
+  finite-sample factor `(G/(G−1))·((n−1)/(n−k))` now uses
+  `K_reference = explicit columns + (1 for the absorbed constant when the
+  design carries no intercept column) + rank(absorbed FE not nested in the
+  cluster | nested ones)` instead of the visible column count, threaded as a
+  new keyword-only `cluster_k_adjustment` on `solve_ols` /
+  `compute_robust_vcov` / `LinearRegression.fit` / `wild_bootstrap_se`
+  (additive, default 0). Two user-visible defects close at once. **D2
+  (anti-conservative):** absorbed FE not nested in the cluster were never
+  counted — clustered SEs were understated by 0.10% (n=2500), 1.30% (n=200),
+  5.51% (n=60); `WooldridgeDiD`'s SEs now match Stata `jwdid`/reghdfe at
+  machine precision on every committed arm AND on the new G≈20..500 subsample
+  ladder golden (historical ratios 1.0280→1.0010, all now 1.0 at ~1e-15..1e-14),
+  `SunAbraham` matches `fixest::sunab` (measured 3.1e-14 relative, was ~1-2%),
+  and the TWFE fixest cluster gates tightened from a 0.5% band to
+  exact / rtol 1e-9. **D1 (the 10.35% idiom split):** `absorb=[unit, time]`
+  and `fixed_effects=[unit, time]` on the same clustered model returned SEs
+  differing by exactly `sqrt((360−2)/(360−66))`; the full-dummy lane now
+  subtracts the cluster-nested dummy rank, so the two documented-equivalent
+  spellings return the identical SE. Applies to DiD/TWFE/MPD (including MPD's
+  built-in period dummies as its time-FE block), SunAbraham, WooldridgeDiD's
+  OLS within path, ImputationDiD's pretrends lead regression, LPDiD's nested
+  dummy blocks, and the wild-cluster-bootstrap reported SE (the correction
+  constant cancels in the studentized statistic, so bootstrap p-values are
+  invariant). StackedDiD (clubSandwich CR1S by construction) and Wooldridge
+  `cohort_trends` full-dummy keep their documented conventions, bit-identical.
+  The non-nested term is the exact rank given the nested set — anchored on a
+  disconnected panel against BOTH reghdfe and fixest
+  `ssc(K.fixef="full", K.exact=TRUE)` at ~1e-16, with the one-df deviation of
+  the nested composition from both references (neither implements
+  nested-drop + exact-remainder) pinned exactly in two new goldens. Tail df
+  conventions are deliberately unchanged (PR C of the program). The
+  `docs/benchmarks.rst` TWFE "SE Rel Diff 0.1%" cell is generated by the
+  R-dependent benchmark refresh and still shows the pre-fix value; the
+  refresh re-run is tracked in TODO.md (expected movement 0.1% → 0.0%).
 - **Absorbed-FE degrees of freedom over-counted on disconnected and hierarchical
   panels.** The df adjustment for absorbed fixed effects (`TwoWayFixedEffects`,
   `SunAbraham`, `DifferenceInDifferences(absorb=)`, `MultiPeriodDiD(absorb=)`)
@@ -25,12 +61,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - **New CI-enforced variance-convention audit matrix**
   (`tests/test_variance_conventions.py` +
   `docs/methodology/variance-conventions.md`): pins, per estimator surface, the
-  visible `k` reaching the shared clustered CR1 denominator and the tail-df
+  `k` reaching the shared clustered CR1 denominator and the tail-df
   convention passed to inference, with each cell classified as a documented
-  defect (scheduled for the 3.9 variance-consolidation program) or a declared
-  legitimate exception with its reason. This is the map for the follow-up PRs
-  that converge the clustered CR1 `k` on the reghdfe/fixest nested convention
-  and the tail df on `t(G−1)`.
+  defect or a declared legitimate exception with its reason. The clustered
+  CR1 `k` cells now pin the converged `K_reference` accounting (see the
+  [M-126] entry above); the remaining defect rows are the tail-df items
+  scheduled for PR C of the 3.9 variance-consolidation program.
 - **`WooldridgeDiD` silently dropped genuine post-treatment effects
   (issue #724).** With `control_group="never_treated"`, the ETWFE design emitted
   every cohort×time indicator *including* the reference period — which is

DEFERRED.md

@@ -50,8 +50,6 @@ exists but parity can't be verified without a local toolchain.
 | Extend `WooldridgeDiD` `method ∈ {logit, poisson}` with `vcov_type ∈ {classical, hc2, hc2_bm}`: composing HC2 leverage + Bell-McCaffrey DOF with the QMLE pseudo-residual sandwich needs derivation + R parity vs `clubSandwich::vcovCR(glm, type="CR2")`. Rejected at `__init__`. | `wooldridge.py` | follow-up | Medium |
 | Multi-constraint CR2 parallel-trends test (AHT/HTZ) for `hc2_bm` fits: DiagnosticReport's PT check routes `vcov_type="hc2_bm"` sources to Bonferroni over the BM-adjusted per-row p-values because the generic chi-square joint Wald would discard the CR2 small-sample correction (see REPORTING.md "hc2_bm parallel-trends policy"). The proper joint test is the AHT/HTZ Wald with a Satterthwaite-style denominator df over the pre-period contrast block; needs derivation for the stacked/pooled WLS-CR2 layout + parity vs `clubSandwich::Wald_test(..., test="HTZ")`. | `diagnostic_report.py`, `linalg.py` | vcov/df round-trip PR | Low |
 | `PreTrendsPower` CS/SA `anticipation=1` R-parity fixture: R `pretrends` has no anticipation parameter, so the Python `_extract_pre_period_params` anticipation filter isn't R-parity-locked. Build a synthetic CS/SA result with `anticipation=1` and assert γ_p matches R's `slope_for_power()`. (Mechanism already covered by MC + full-VCV tests.) | `tests/test_methodology_pretrends.py`, `generate_pretrends_golden.R` | PR-C | Low |
-| Harmonize SunAbraham's HC1 within-transform finite-sample correction with `fixest::sunab()` — SA applies `n/(n-k_dm)`, fixest applies `n/(n-k_total)` (counts absorbed FE); ~1-2% SE difference, documented as a "Deviation from R" and pinned at `atol=5e-3`. Either thread `df_adjustment` or keep as an intentional, R-verified difference. Mechanism + measured inventory now in `docs/methodology/variance-conventions.md` (defect D2 family; 3.9 consolidation PR B). | `sun_abraham.py`, `linalg.py` | follow-up | Low |
-| Absorbed-FE **clustered** CR1 with *non-nested* FE: for `absorb=[FE1,FE2], cluster=FE1` (e.g. `absorb=["unit","time"], cluster="unit"`), `fixest` counts the non-nested FE (time) in the CR1 `(n-1)/(n-k)` finite-sample denominator, but the clustered path uses only `k_visible`. D4 harmonized the *non-clustered* classical/hc1 full-K scale (`_absorbed_fe_vcov_scale`) and left the clustered path unchanged — correct for FE nested in the cluster, a small deviation for non-nested FE (documented in REGISTRY within-transform note). Thread a non-nested `df_adjustment` into the clustered CR1 factor; verify vs `fixest::feols(..., cluster=)`. Closed form now derived and externally verified (Stata ~1e-15 / fixest ~1e-12) — see `docs/methodology/variance-conventions.md` defect D2; scheduled as PR B of the 3.9 consolidation program. | `linalg.py`, `estimators.py` | SE-audit D4 | Low |
 | Rust multiplier-bootstrap weight RNG (`generate_bootstrap_weights_batch`) seeds `Xoshiro256PlusPlus::seed_from_u64(seed+i)` per row; audit Python callers (`sdid.py`, `efficient_did_bootstrap.py`, `bootstrap_utils.py`) for parity-test gaps and, where a numpy-canonical equivalent exists, pre-generate in Python and pass through PyO3 (same fix shape as TROP RNG parity #354). | `rust/src/bootstrap.rs`, `bootstrap_utils.py` | follow-up | Medium |
 | `SyntheticDiD` bootstrap cross-language parity anchor vs R `synthdid::vcov(method="bootstrap")` or Julia `Synthdid.jl` (refit-native). Same-library validation is in place; Julia is the cleanest target. Tolerance ~1e-6 (BLAS+RNG paths preclude 1e-10). | `benchmarks/R/`, `benchmarks/julia/`, `tests/` | follow-up | Low |
 | CS R helpers hard-code `xformla = ~1`; no covariate-adjusted R benchmark for the IRLS path. | `tests/test_methodology_callaway.py` | #202 | Low |

METHODOLOGY_REVIEW.md

@@ -650,7 +650,7 @@ and covariate-adjusted specifications.)
 **Deviations from the paper / from R / library extensions:** See REGISTRY.md `## WooldridgeDiD (ETWFE)` → `### Deviations from the paper / from R / library extensions` block for the consolidated list (HC1 finite-sample factor, QMLE sandwich `(n-1)/(n-k)` term, nonlinear-vs-fixest direct QMLE, logit cohort+time additive dummies, anticipation + aggregation, cell-count default with opt-in cohort-share).
 
 **Outstanding Concerns:**
-- **Stata `jwdid` golden values**: RESOLVED (#729, extended here). `benchmarks/data/etwfe_cs_stata_golden.json` ships four arms -- `csdid`, `jwdid`, `jwdid_never` (the issue #724 reference-period anchor) and `jwdid_alltreated` (the W2025 Section 5.4 anchor) -- generated with local StataSE 19 and pinned by `tests/test_etwfe_cs_stata_parity.py`. Point estimates match to ~1e-15. What remains open is the `hc1` SE gap (library SEs uniformly SMALLER than jwdid's by a cluster-count-dependent factor: 1.0280@G=20, 1.0132@G=40, 1.00264@G=191, 1.0010@G=500), tracked in TODO.md; each arm pins its own measured ratio because the constant does not transfer between cluster counts. No QMLE (logit/Poisson) golden exists yet -- every arm is linear `jwdid` -- so the QMLE cluster-SE comparison in DEFERRED.md stays open. R `etwfe` side covered in PR-B Stage D.
+- **Stata `jwdid` golden values**: RESOLVED (#729, extended here). `benchmarks/data/etwfe_cs_stata_golden.json` ships four arms -- `csdid`, `jwdid`, `jwdid_never` (the issue #724 reference-period anchor) and `jwdid_alltreated` (the W2025 Section 5.4 anchor) -- plus the G≈20..500 subsample `ladder` block, generated with local StataSE 19 and pinned by `tests/test_etwfe_cs_stata_parity.py`. Point estimates match to ~1e-15. The historical `hc1` SE gap (library SEs uniformly SMALLER by a cluster-count-dependent factor: 1.0280@G=20 down to 1.0010@G=500) was defect D2 of the 3.9 variance program and is CLOSED by the K_reference convergence: every arm and every ladder rung now sits at SE ratio 1.0 (spreads ~1e-15..1e-14; `docs/methodology/variance-conventions.md`). No QMLE (logit/Poisson) golden exists yet -- every arm is linear `jwdid` -- so the QMLE cluster-SE comparison in DEFERRED.md stays open. R `etwfe` side covered in PR-B Stage D.
 - **Response-scale APE / log-link bridge for Poisson + logit R parity** (DEFERRED.md WooldridgeDiD follow-up cluster row, added in PR-B): direct cell-level numerical parity between diff-diff's response-scale ATT and R `etwfe` log-link coefficients requires either `emfx()`-based APE extraction on the R side or link-function inversion with baseline-mean adjustment.
 - **QMLE sandwich Stata-parity `qmle` weight type** (DEFERRED.md WooldridgeDiD follow-up cluster row): diff-diff's `(G/(G-1)) × ((n-1)/(n-k))` is conservative vs Stata's `G/(G-1)` only; awaiting Stata golden values to confirm material difference.
 - **Repeated cross-sections** (paper p. 2581 → Deb et al. 2024): not in 2025 paper's main body; future PR.