From 9e8daa8afd73cbaedd3deca13651462231d8c2de Mon Sep 17 00:00:00 2001 From: igerber Date: Thu, 30 Jul 2026 21:49:19 -0400 Subject: [PATCH 1/3] feat(3.9): converge clustered CR1 on K_reference - D1+D2 fix (M-126) PR B of the 3.9 variance-consolidation program: the clustered hc1 CR1 finite-sample factor (G/(G-1))*((n_eff-1)/(n_eff-k)) now uses K_reference = explicit columns + (1 for the absorbed constant when X carries no intercept col) + rank(absorbed FE not nested in the cluster | nested ones) threaded as a NEW keyword-only signed `cluster_k_adjustment: int = 0` on solve_ols / compute_robust_vcov / LinearRegression.fit / wild_bootstrap_se, computed by the new utils helpers (absorbed_fe_cr1_k_increment, cluster_nested_fe_dims, shared factorizations with absorbed_fe_rank). Fixes two defects from the measured variance-conventions inventory: - D2 (anti-conservative): absorbed FE not nested in the cluster were never counted - clustered SEs understated by 0.10% (n=2500) / 1.30% (n=200) / 5.51% (n=60). WooldridgeDiD now matches Stata jwdid/reghdfe at machine precision on all committed mpdta arms AND the new G~20..500 subsample ladder golden (historical ratios 1.0280@G=20 down to 1.0010@G=500, all now 1.0 at spreads ~1e-15..1e-14; the ladder doubles as the K-accounting probe: reghdfe df_a == increment - 1 per rung). SunAbraham matches fixest::sunab at 3.1e-14 rel (was ~1-2%; closes the DEFERRED row). TWFE fixest cluster gates tightened from a 0.5% band to exact / rtol 1e-9. - D1 (the 10.35% idiom split): absorb=[unit,time] vs 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 both documented-equivalent spellings return the identical externally-anchored SE (holds under collinearity drops and on disconnected panels - kernel k is the design rank). Wired surfaces: DiD/MPD/TWFE (absorb increment; fixed_effects negative nested rank; MPD's built-in period dummies are its time-FE block and follow the nested convention), SunAbraham, WooldridgeDiD OLS within, ImputationDiD pretrends lead regression, LPDiD nested dummy blocks (with an InvalidClusterKAdjustment re-raise so the broad except cannot swallow the contract), and the WCB corr constant (cancels in |t*| vs |t0|: bootstrap p-values invariant; reported SE aligns with the analytical convention). StackedDiD (clubSandwich CR1S by construction) and Wooldridge cohort_trends full-dummy keep their documented conventions bit-identically. Rust lanes apply the correction as an exact None-preserving scalar rescale of the finished vcov in the Python wrappers (zero Rust changes, <= 1 ulp; zero-adjustment surfaces bit-identical). Contracts: front-door InvalidClusterKAdjustment(ValueError) validation on every route (type half checked unconditionally at entry - classical/survey fits and degenerate WCB returns reject non-ints too); fail-closed all-NaN vcov on all three saturation sides (n_eff-k <= 0, n_eff-k_inf <= 0, k_inf <= 0) on both backends, pinned at the exact boundaries with one-step finite recovery; survey fits are the documented inert VALUE exception. External anchors beyond Stata/fixest parity: the exact non-nested RANK term is anchored on a DISCONNECTED two-way panel against BOTH reghdfe (2.8e-17) and fixest ssc(K.fixef="full", K.exact=TRUE) (5.6e-16) via two new goldens (fixest_cr1_nonnested_golden.json, reghdfe_kref_golden.json - deterministic RNG-free DGP); no external reference implements nested-drop + exact-remainder (fixest's K.exact composes incoherently with its nested drop; reghdfe's pairwise correction skips nested-dropped dims), so that one-df deviation is pinned exactly as sqrt((n-10)/(n-11)) against both and documented as a labeled REGISTRY deviation. N>=3 absorbed dims inherit the documented D3 sum(levels)-N+1 approximation (limitation pinned; TODO row). Tail df (df_, dof_vec, inference_df, df_convention) deliberately unchanged - PR C scope. Ledger: docs/v4-deprecations.yaml M-126 (+ v4-matrix snapshot 105); REGISTRY/variance-conventions/CHANGELOG/benchmark READMEs rewritten to the converged convention; DEFERRED rows 53+54 closed; TODO ladder row closed into the benchmarks.rst refresh row; perf note re-measured (increment 3.2 ms/call at 186k rows, tracked). --- CHANGELOG.md | 46 +- DEFERRED.md | 2 - METHODOLOGY_REVIEW.md | 2 +- TODO.md | 4 +- .../R/generate_fixest_cr1_nonnested_golden.R | 164 +++++ .../R/generate_fixest_did_twfe_golden.R | 14 +- benchmarks/README.md | 16 +- benchmarks/data/etwfe_cs_stata_golden.json | 153 +++++ .../data/fixest_cr1_nonnested_golden.json | 61 ++ benchmarks/data/reghdfe_kref_golden.json | 35 + benchmarks/stata/README.md | 58 +- benchmarks/stata/generate_etwfe_cs_golden.do | 105 ++- .../stata/generate_reghdfe_kref_golden.do | 207 ++++++ diff_diff/estimators.py | 123 +++- diff_diff/imputation.py | 15 + diff_diff/linalg.py | 227 ++++++- diff_diff/lpdid.py | 30 +- diff_diff/sun_abraham.py | 19 +- diff_diff/twfe.py | 54 +- diff_diff/utils.py | 199 +++++- diff_diff/wooldridge.py | 25 + docs/methodology/REGISTRY.md | 124 ++-- docs/methodology/variance-conventions.md | 90 ++- docs/performance-plan.md | 24 +- docs/v4-deprecations.yaml | 13 + docs/v4-design.md | 6 +- tests/test_estimators_vcov_type.py | 25 +- tests/test_etwfe_cs_stata_parity.py | 215 ++++-- tests/test_fixest_did_twfe_parity.py | 22 +- tests/test_linalg.py | 302 +++++++++ tests/test_methodology_sun_abraham.py | 21 +- tests/test_methodology_twfe.py | 54 +- tests/test_methodology_wooldridge.py | 32 +- tests/test_v4_matrix.py | 19 +- tests/test_variance_conventions.py | 634 +++++++++++++++++- tests/test_wooldridge.py | 16 +- 36 files changed, 2824 insertions(+), 332 deletions(-) create mode 100644 benchmarks/R/generate_fixest_cr1_nonnested_golden.R create mode 100644 benchmarks/data/fixest_cr1_nonnested_golden.json create mode 100644 benchmarks/data/reghdfe_kref_golden.json create mode 100644 benchmarks/stata/generate_reghdfe_kref_golden.do diff --git a/CHANGELOG.md b/CHANGELOG.md index 4fc45db91..7929b2746 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/DEFERRED.md b/DEFERRED.md index 10d99e00c..066f7cd12 100644 --- a/DEFERRED.md +++ b/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 | diff --git a/METHODOLOGY_REVIEW.md b/METHODOLOGY_REVIEW.md index e5147f8e8..e17fca4d1 100644 --- a/METHODOLOGY_REVIEW.md +++ b/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. diff --git a/TODO.md b/TODO.md index 77421c89f..57c9754eb 100644 --- a/TODO.md +++ b/TODO.md @@ -24,7 +24,7 @@ Related tracking surfaces: | `absorbed_fe_rank` N>=3 general rank: the helper keeps `sum(levels-1)` for 3+ absorbed dims, exact for independent connected dims but an over-count for duplicated/nested triples (measured `a(5),b(4),c==b(4)`: true 7 vs formula 10) and for disconnected N-way graphs. Two-way is component-exact. Deriving general N-way FE rank is a hypergraph problem; do it with a reference (fixest's `fixef.rm`/reghdfe df_a) rather than a guess. See `docs/methodology/variance-conventions.md` D3. | `diff_diff/utils.py` | #variance-inventory | Mid | Low | | `SyntheticControl` conformal (CWZ 2021) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies — the residual-permutation shortcut is only valid for time-permutation-invariant proxies (SC/Lasso/DiD); an AR proxy needs innovation permutation. | `diff_diff/conformal.py`, `diff_diff/synthetic_control_results.py` | CWZ-2021 | Heavy | Low | | Make the post-fit `results.aggregate("event_study")` container consumable downstream. `EventStudyResults` is rejected by all THREE consumers that read a CS event study — `compute_honest_did` (`honest_did.py`, dispatches on `CallawaySantAnnaResults` and raises `TypeError`), `compute_pretrends_power` (`pretrends.py`, same), and `plot_event_study` (`visualization`, same) — so `fit(aggregate="event_study")` is still the only route for them and their error messages say so explicitly. Needs an `EventStudyResults` branch in each extraction path (consuming `event_time` / `is_reference` / `vcov` / `vcov_index` / per-row `df`) PLUS `base_period` and `anticipation` provenance, which the unified container does not carry and HonestDiD needs for its universal-base-period warning and pre-period classification. Gate with end-to-end tests: `compute_honest_did(res.aggregate("event_study"))` at `base_period="universal"`, and `compute_pretrends_power(...)` at `anticipation=1`. | `diff_diff/honest_did.py`, `diff_diff/pretrends.py`, `diff_diff/results_base.py` | #726 | Mid | Medium | -| Derive and fix the `WooldridgeDiD` `hc1` SE gap vs Stata `jwdid`. **MECHANISM NOW DERIVED** (see `docs/methodology/variance-conventions.md`, defect D2): the clustered CR1 denominator uses `k_visible` and never counts absorbed FE not nested in the cluster; the closed form `K_reference = explicit cols + (1 if no intercept col) + rank(non-nested FE | nested)` reproduces jwdid/reghdfe to ~1e-15 on three arms and fixest to ~1e-12 on two, retrodicting the G=20/G=40 rungs. Remaining work (PR B of the 3.9 consolidation program) is the threading + Rust port + Stata ladder, not derivation. Measured: every SE is uniformly SMALLER than `jwdid`'s, by 1.0280 at G=20 / 1.0132 at G=40 / 1.00264 at G=191 / 1.0010 at G=500 (ATT(g,t) points match exactly). The gap tracks `sqrt(G/(G-1))` but sits consistently above it, and `solve_ols` already applies the full CR1 `(G/(G-1))*((n-1)/(n-k))`, so a missing cluster factor is ruled out -- the likely source is the within-transform `k` accounting vs `hdfe`/`reghdfe`'s. Derive the exact factor FIRST; a `sqrt(G/(G-1))` patch would match only approximately and would force a loose tolerance on the parity test, defeating its purpose. Move the REGISTRY note and the pinned ratio in `tests/test_etwfe_cs_stata_parity.py` together with the fix. **Required artifact:** a committed subsample LADDER in the Stata golden (rosters = first N units per `first_treat` by sorted `countyreal`; rungs spanning G≈20..500, storing Stata `G`/`n`/`df_a`/`rank`/`df_r` and per-cell `att`/`se`), with parameterized ratio assertions. Only G=500 (full panel) and G=191 (all-eventually-treated arm) are pinned today, so the few-cluster behavior — where the gap is materially largest (~2.8%) — is ungated, and the ladder is also the instrument for comparing Stata's `df_a`/`rank` against the library's within-transform `k`, which is the leading hypothesis for the factor. | `diff_diff/wooldridge.py`, `diff_diff/linalg.py` | #723-followup | Mid | Medium | +| Re-run the R-dependent benchmark refresh so `docs/benchmarks.rst`'s TWFE "SE Rel Diff 0.1%" cell reflects the 3.9 K_reference convergence (expected 0.1% -> 0.0%; the table is generated, never hand-edited — the movement is noted in the CHANGELOG entry). | `docs/benchmarks.rst`, `benchmarks/R/` | #variance-inventory | Quick | Low | | `SunAbraham`: a cohort not observed at its own reference relative period (`e = -1 - anticipation`) makes that cohort's block collinear, so QR drops an unnamed column (`dropping 1 of 12 columns (column 9)`) and `overall_att` comes back **NaN**. Found by auditing the sibling estimator while fixing the ETWFE analogue (#724); PRE-EXISTING, not introduced there. Lower severity than #724 — that returned a silently WRONG finite number, this returns NaN with a rank warning — but the event-study surface still looks complete, so a user may not notice the loss. SA already omits its reference explicitly and tracks `_reference_observed`, so the fix is per-cohort support for that flag rather than the ETWFE-style redesign. | `diff_diff/sun_abraham.py` | #724-audit | Mid | Low | | Define `N_g` (W2025 Eqs. 7.4/7.6) for UNBALANCED panels where comparison-support filtering removes every observation of some units in an estimated cohort, then replace the fail-closed guard with the defined behavior. `_n_g_per_cohort` is read off the final sample, so those units vanish from the cohort-share weights; measured on a cohort supplied with 100 units of which 90 appear only at a dropped period, `aggregate(weights="cohort_share")` moves 1.8078 -> 3.8157. The paper assumes a balanced panel and does not say whether `N_g` counts the supplied cohort or the surviving units, and the two disagree materially, so `aggregate` currently raises naming the cohorts and counts ([M-125]); `weights="cell"` is unaffected and balanced panels never trip it. Settle the estimand (likely: count the supplied cohort, since ATT(g,t) is a cohort-level quantity, but that weights units with no retained observation) and gate with a test computing Eq. 7.4 by hand on unequal cohort sizes. | `diff_diff/wooldridge_results.py`, `diff_diff/wooldridge.py` | #729-followup | Mid | Medium | | `WooldridgeDiD` + `survey_design=` does not support DOMAIN ESTIMATION, so BOTH row-deleting paths are currently REFUSED (`NotImplementedError`, all three methods) rather than performed: unidentified-cohort exclusion ([M-123]) and comparison-support period filtering ([M-125]). One fix unblocks both. Implementing it properly means zero-padding the excluded rows' weights while retaining strata/PSU/FPC, per REGISTRY *Subpopulation Analysis (Phase 6)* / Lumley (2004) 3.4, so TSL variance and `df_survey = n_PSU - n_strata` use the full design (naive deletion measured 22 -> 14 on a two-stratum panel). `SurveyDesign.subpopulation()` already implements the contract and SpilloverDiD Wave E.3 is the in-repo precedent; the blocker is that the weighted within-transform rejects zero-weight units, shared machinery behind 7 estimators. Landing it would turn both refusals back into supported fits. Gate with a `SurveyDesign.subpopulation()` parity test on ATT, TSL SE and survey df where the excluded cohort exhausts a PSU. | `diff_diff/wooldridge.py`, `diff_diff/utils.py` | #724-codex-R4/R5 | Heavy | Medium | @@ -49,7 +49,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| -| Reuse the demeaner's factorized codes in `absorbed_fe_rank` instead of re-factorizing: the helper adds ~1.9 ms per absorbed fit at 186k rows (7.7% of the fastest Rust-served TWFE fit; see `docs/performance-plan.md` "Component-aware absorbed-FE rank"), and both the helper and `demean_by_groups` factorize the same group columns. Threading the codes through the four call sites halves the factorize work; the `connected_components` call itself is ~1.1 ms. Deliberately not done in the correctness PR. | `diff_diff/utils.py` | #variance-inventory | Quick | Low | +| Reuse the demeaner's factorized codes in `absorbed_fe_rank`/`absorbed_fe_cr1_k_increment` instead of re-factorizing: at 186k rows the rank helper adds ~1.9 ms per absorbed fit (7.7% of the fastest Rust-served TWFE fit) and the K_reference increment ~3.2 ms per clustered-hc1 absorbed fit (~13%; see `docs/performance-plan.md` "Component-aware absorbed-FE rank"), and the helpers and `demean_by_groups` factorize the same group columns. Threading the codes through the call sites halves the factorize work; the `connected_components` call itself is ~1.1 ms. Deliberately not done in the correctness PRs. | `diff_diff/utils.py` | #variance-inventory | Quick | Low | | `EfficientDiD` conditional path: the largest remaining O(n) stage is the sieve/nuisance construction outside the tiled pass (~9s at 10k). (The `_ridge_solve_weights` Python-prep shave landed 2026-07-07 — the `omega_stack[rest]` fancy-index copy and tail scatter are skipped when no row is zero-masked, byte-identical outputs; the `zero_mask` abs scan itself remains, needed for correctness.) | `efficient_did_covariates.py` | CS-scaling | Mid | Low | | `_rq_fit` LP assembly is dense (`A_eq = [X, I, -I]` with dense identity blocks, rebuilt per cell fit): a `scipy.sparse` construction would cut memory and likely HiGHS time for large cells / bootstrap-heavy covariate CiC/QDiD fits. CAVEAT before doing it: a different matrix representation can change HiGHS's vertex selection at degenerate/tied QR optima - end-to-end covariate goldens are tie-selection-gated (fine), but the `qr_cases` tight coefficient matches may shift to the equal-loss branch; re-run the parity suite and re-calibrate if needed. | `diff_diff/changes_in_changes.py::_rq_fit` | covariates PR | Quick | Low | | Evaluate flipping `DIFF_DIFF_SOLVE_OLS_FASTPATH` default-ON after an opt-in soak (the 2026-07 certified normal-equations Cholesky fast path, both backends). A flip needs: golden/parity-suite recapture at the tol-bounded posture (fitted ~1e-8 abs / SE ~1e-6 rel — the default today is byte-pinned in several benchmark conventions), certification-rate telemetry across real workloads (any decline is silent-correct but forfeits the speedup), and the staged default-flip protocol used for `df_convention` (v4-class change). Lifecycle tracked in docs/v4-deprecations.yaml (M-008). | `diff_diff/linalg.py::_resolve_solve_ols_fastpath`, `rust/src/linalg.rs::solve_ols_chol` | CS-scaling | Mid | Low | diff --git a/benchmarks/R/generate_fixest_cr1_nonnested_golden.R b/benchmarks/R/generate_fixest_cr1_nonnested_golden.R new file mode 100644 index 000000000..ca3295aa5 --- /dev/null +++ b/benchmarks/R/generate_fixest_cr1_nonnested_golden.R @@ -0,0 +1,164 @@ +# Golden: fixest clustered CR1 with a NON-NESTED absorbed-FE rank term. +# +# The committed clustered fixest arms (fixest_did_twfe_golden.json) exercise +# only connected two-way panels, where the non-nested rank term is the trivial +# delta = T - 1. What no golden covered is the RANK of the non-nested absorbed +# FE on an irregular design -- the exact quantity the K_reference convention +# adds to the clustered CR1 k (variance-conventions.md D2). +# +# DGP: a DISCONNECTED two-way [unit, time] panel -- two blocks with disjoint +# period ranges (C = 2 bipartite components) -- plus a fully-crossed connected +# control panel. Three arms: +# +# 1. disconnected$crossed_cluster (THE parity anchor): clustered by +# c5 = (unit + time) %% 5, which crosses both dims, so NOTHING is nested +# and the K_reference increment is the exact full span rank +# U + T - C = 28 (K = 29 with x). fixest ssc(K.fixef = "full", +# K.exact = TRUE) computes exactly that (df.K = 29, agreeing with iid +# K.exact and with reghdfe's df_a = 28 + x); the default approximate +# count gives df.K = 30, and K.exact under the DEFAULT +# K.fixef = "nested" gives df.K = 28 -- on the clustered path fixest's +# nested handling removes 1 df even when nothing is nested (recorded as +# k_exact_nested for documentation). All three are recorded; diff-diff +# must match the full+K.exact side at 1e-12 and differ from the other +# two (the same deviation-from-R-default as the kexact sibling golden). +# This arm also pins the zero-nested-dim max(..., 1) floor end-to-end. +# +# 2. disconnected$unit_cluster (DOCUMENTATION arm, no parity target): +# clustered by unit (nested, dropped; exact remainder rank of time given +# unit is 28 - 20 = 8, so K_reference = x + constant + 8 = 10, denominator +# n - 10). NO fixest ssc configuration produces that composition: +# the default nested+approximate count gives df.K = 11 (time counted +# T - 1 = 9), and ssc(K.exact = TRUE) composes INCOHERENTLY with the +# nested drop (df.K = 28 -- it removes only 1 df for the fully-nested +# unit FE). Both are recorded so the consuming test can pin that the +# library matches NEITHER: its exact-remainder K = 10 sits one df below +# fixest's default (Stata reghdfe lands on the same 11 via its own +# approximate remainder -- see reghdfe_kref_golden.json). NO external +# reference implements nested-drop + exact-remainder; the library's +# convention is the consistent extension of the exact-rank principle, +# and the test pins the deviation exactly: +# se_default / se_library == sqrt((n - 10) / (n - 11)). +# +# 3. connected (control): fully-crossed two-way panel clustered by unit, +# default ssc -- at C = 1 with independent dims the two conventions +# coincide (rank term T - 1), pinning the regular case. +# +# Regenerate: Rscript benchmarks/R/generate_fixest_cr1_nonnested_golden.R +suppressMessages(library(fixest)) +suppressMessages(library(jsonlite)) + +set.seed(11) + +# --- Disconnected two-way panel (C = 2) -------------------------------------- +d1 <- rbind( + expand.grid(unit = 0:9, time = 0:4), + expand.grid(unit = 10:19, time = 5:9) +) +d1$x <- rnorm(nrow(d1)) +d1$out <- rnorm(nrow(d1)) + 0.5 * d1$x + 0.2 * (d1$unit %% 3) + 0.1 * d1$time +d1$c5 <- (d1$unit + d1$time) %% 5 + +m1_unit_default <- feols(out ~ x | unit + time, data = d1, cluster = ~unit) +m1_unit_kexact <- feols(out ~ x | unit + time, data = d1, cluster = ~unit, + ssc = ssc(K.exact = TRUE)) +m1_cross_kexact <- feols(out ~ x | unit + time, data = d1, cluster = ~c5, + ssc = ssc(K.fixef = "full", K.exact = TRUE)) +m1_cross_kexact_nested <- feols(out ~ x | unit + time, data = d1, cluster = ~c5, + ssc = ssc(K.exact = TRUE)) +m1_cross_default <- feols(out ~ x | unit + time, data = d1, cluster = ~c5) + +# --- Connected control (C = 1), clustered by unit, default ssc --------------- +d2 <- expand.grid(unit = 0:19, time = 0:9) +d2$x <- rnorm(nrow(d2)) +d2$out <- rnorm(nrow(d2)) + 0.5 * d2$x + 0.2 * (d2$unit %% 3) + 0.1 * d2$time + +m2 <- feols(out ~ x | unit + time, data = d2, cluster = ~unit) + +golden <- list( + meta = list( + generator = "benchmarks/R/generate_fixest_cr1_nonnested_golden.R", + r_version = paste(R.version$major, R.version$minor, sep = "."), + fixest_version = as.character(packageVersion("fixest")), + description = paste( + "Clustered CR1 with a non-nested absorbed-FE RANK term on a", + "disconnected two-way panel (C=2). crossed_cluster: nothing nested,", + "K.exact=TRUE is the parity anchor (exact span rank 28, df.K=29) vs", + "the default approximate count (df.K=30). unit_cluster: documentation", + "arm -- no fixest ssc reproduces nested-drop + exact-remainder", + "(library K=10 vs default df.K=11; reghdfe agrees with fixest's 11", + "via its own approximate remainder, see reghdfe_kref_golden.json;", + "the library deviation is pinned as se ratio sqrt((n-10)/(n-11))).", + "connected: control arm where the conventions agree." + ) + ), + disconnected = list( + data = list( + unit = d1$unit, + time = d1$time, + x = d1$x, + out = d1$out, + c5 = d1$c5 + ), + n_obs = nrow(d1), + unit_cluster = list( + n_clusters = length(unique(d1$unit)), + coef = unname(coef(m1_unit_default)[["x"]]), + default = list( + se = unname(se(m1_unit_default)[["x"]]), + df_k = degrees_freedom(m1_unit_default, "k") + ), + nested_k_exact = list( + se = unname(se(m1_unit_kexact)[["x"]]), + df_k = degrees_freedom(m1_unit_kexact, "k") + ) + ), + crossed_cluster = list( + n_clusters = length(unique(d1$c5)), + coef = unname(coef(m1_cross_kexact)[["x"]]), + k_exact = list( + se = unname(se(m1_cross_kexact)[["x"]]), + df_k = degrees_freedom(m1_cross_kexact, "k") + ), + k_exact_nested = list( + se = unname(se(m1_cross_kexact_nested)[["x"]]), + df_k = degrees_freedom(m1_cross_kexact_nested, "k") + ), + default = list( + se = unname(se(m1_cross_default)[["x"]]), + df_k = degrees_freedom(m1_cross_default, "k") + ) + ) + ), + connected = list( + data = list( + unit = d2$unit, + time = d2$time, + x = d2$x, + out = d2$out + ), + n_obs = nrow(d2), + n_clusters = length(unique(d2$unit)), + coef = unname(coef(m2)[["x"]]), + cluster_default = list( + se = unname(se(m2)[["x"]]), + df_k = degrees_freedom(m2, "k") + ) + ) +) + +path <- "benchmarks/data/fixest_cr1_nonnested_golden.json" +write_json(golden, path, digits = NA, auto_unbox = TRUE, pretty = TRUE) +cat("wrote", path, "\n") +dd <- golden$disconnected +cat(sprintf("unit_cluster: coef=%.15f default df.K=%d se=%.15f | nested+K.exact df.K=%d se=%.15f\n", + dd$unit_cluster$coef, + dd$unit_cluster$default$df_k, dd$unit_cluster$default$se, + dd$unit_cluster$nested_k_exact$df_k, dd$unit_cluster$nested_k_exact$se)) +cat(sprintf("crossed_cluster: coef=%.15f k_exact df.K=%d se=%.15f | default df.K=%d se=%.15f\n", + dd$crossed_cluster$coef, + dd$crossed_cluster$k_exact$df_k, dd$crossed_cluster$k_exact$se, + dd$crossed_cluster$default$df_k, dd$crossed_cluster$default$se)) +cat(sprintf("connected: coef=%.15f df.K=%d se=%.15f\n", + golden$connected$coef, + golden$connected$cluster_default$df_k, golden$connected$cluster_default$se)) diff --git a/benchmarks/R/generate_fixest_did_twfe_golden.R b/benchmarks/R/generate_fixest_did_twfe_golden.R index 712c17eca..b61ca9fc4 100644 --- a/benchmarks/R/generate_fixest_did_twfe_golden.R +++ b/benchmarks/R/generate_fixest_did_twfe_golden.R @@ -11,12 +11,14 @@ # 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. +# plain-OLS DiD path locks hetero AND cluster CR1 at machine precision, and +# the TWFE cluster SE is likewise locked (exact on the balanced arm, rtol +# 1e-9 on the hetero arm) — the historical ~0.25% non-nested-FE ssc band +# was defect D2, closed by the 3.9 K_reference convergence (the within +# CR1 factor now counts the non-nested time FE exactly as fixest's +# ssc(K.fixef="nested") does). 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 diff --git a/benchmarks/README.md b/benchmarks/README.md index dc667d554..f7604893d 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -55,14 +55,18 @@ STATA=/Applications/Stata/StataSE.app/Contents/MacOS/stata-se $STATA -b do benchmarks/stata/requirements.do # one-time SSC install $STATA -b do benchmarks/stata/generate_lpdid_ra_golden.do $STATA -b do benchmarks/stata/generate_imputation_loo_golden.do +$STATA -b do benchmarks/stata/generate_etwfe_cs_golden.do +$STATA -b do benchmarks/stata/generate_reghdfe_kref_golden.do ``` The `LPDiD` arm uses only **native** Stata commands (`teffects`), pinned by `version 19`. The `ImputationDiD` arm depends on SSC packages -(`did_imputation`/`reghdfe`/`ftools`/`require`), which `version 19` does NOT pin +(`did_imputation`/`reghdfe`/`ftools`/`require`), the ETWFE/CS arm on +`drdid`/`csdid`/`jwdid`/`hdfe`, and the reghdfe K_reference arm on `reghdfe`; +`version 19` does NOT pin SSC packages (SSC has no version history) — install them once via `requirements.do` (the -generators do not auto-install) and each golden records `ssc_versions` for drift -detection. See `stata/README.md`. +generators do not auto-install) and each golden records version/checksum +metadata for drift detection. See `stata/README.md`. ## Directory Structure @@ -81,7 +85,9 @@ benchmarks/ │ ├── README.md # Stata arm docs │ ├── requirements.do # one-time SSC install (did_imputation etc.) │ ├── generate_lpdid_ra_golden.do # LPDiD RA SE vs teffects ra -│ └── generate_imputation_loo_golden.do # ImputationDiD LOO SE vs did_imputation leaveout +│ ├── generate_imputation_loo_golden.do # ImputationDiD LOO SE vs did_imputation leaveout +│ ├── generate_etwfe_cs_golden.do # ETWFE/CS vs jwdid + csdid (+ subsample ladder) +│ └── generate_reghdfe_kref_golden.do # clustered CR1 K_reference vs reghdfe (disconnected panel) ├── python/ │ ├── utils.py # Common utilities │ ├── benchmark_callaway.py # CallawaySantAnna @@ -105,6 +111,8 @@ benchmarks/ | `DifferenceInDifferences` | `fixest::feols` | Standard DiD | ✓ Integrated | | `LPDiD` (RA SE) | Stata `teffects ra ... atet` | Dube, Girardi, Jorda & Taylor (2025) | ✓ Integrated | | `ImputationDiD` (LOO SE) | Stata `did_imputation, leaveout` | Borusyak, Jaravel & Spiess (2024) App. A.9 | ✓ Integrated | +| `WooldridgeDiD` / `CallawaySantAnna` | Stata `jwdid` / `csdid` (+ G≈20..500 SE ladder) | Wooldridge (2025) / Callaway & Sant'Anna (2021) | ✓ Integrated | +| Clustered CR1 `K_reference` | Stata `reghdfe` + R `fixest` (disconnected-panel arms) | reghdfe/fixest ssc conventions | ✓ Integrated | | `HonestDiD` | `HonestDiD::createSensitivityResults` | Rambachan & Roth (2023) | Planned | Note: HonestDiD benchmark scripts exist but are not yet integrated into the main runner. diff --git a/benchmarks/data/etwfe_cs_stata_golden.json b/benchmarks/data/etwfe_cs_stata_golden.json index 9127f9d0c..9a2a9a36f 100644 --- a/benchmarks/data/etwfe_cs_stata_golden.json +++ b/benchmarks/data/etwfe_cs_stata_golden.json @@ -78,5 +78,158 @@ "2006.first_treat#2006.year#c.__tr__": {"att": 0.018480380807537715, "se": 0.023079607844173931}, "_cons": {"att": 6.0092361449684812, "se": 0.002610873988496222} } + }, + "ladder": { + "roster_rule": "first N units per first_treat cohort by ascending countyreal", + "jwdid_cmd": "jwdid lemp, ivar(countyreal) tvar(year) gvar(first_treat)", + "rungs": { + "5": { + "n_per_cohort": 5, + "G": 20, + "n": 100, + "df_a": 4, + "df_a_initial": 25, + "df_a_nested": 20, + "df_a_redundant": 21, + "rank": 7, + "df_r": 19, + "cells": { + "2004bn.first_treat#2004bn.year#c.__tr__": {"att": -0.10832957930457381, "se": 0.076893226627397518}, + "2004bn.first_treat#2005.year#c.__tr__": {"att": -0.22965545092069178, "se": 0.11209286160589602}, + "2004bn.first_treat#2006.year#c.__tr__": {"att": -0.32585134732690207, "se": 0.13186371530957683}, + "2004bn.first_treat#2007.year#c.__tr__": {"att": -0.23682831134208129, "se": 0.2108974734725157}, + "2006.first_treat#2006.year#c.__tr__": {"att": -0.024998538945235873, "se": 0.094600926914735178}, + "2006.first_treat#2007.year#c.__tr__": {"att": 0.023014548893669263, "se": 0.19677660222842464}, + "2007.first_treat#2007.year#c.__tr__": {"att": 0.18969187074521637, "se": 0.22132822310637179}, + "_cons": {"att": 5.7042511846897668, "se": 0.037061606201004925} + } + }, + "10": { + "n_per_cohort": 10, + "G": 40, + "n": 200, + "df_a": 4, + "df_a_initial": 45, + "df_a_nested": 40, + "df_a_redundant": 41, + "rank": 7, + "df_r": 39, + "cells": { + "2004bn.first_treat#2004bn.year#c.__tr__": {"att": -0.083434385847123746, "se": 0.04250992636950051}, + "2004bn.first_treat#2005.year#c.__tr__": {"att": -0.17253213090254785, "se": 0.066186048646319828}, + "2004bn.first_treat#2006.year#c.__tr__": {"att": -0.23902552822183262, "se": 0.080763375004200832}, + "2004bn.first_treat#2007.year#c.__tr__": {"att": -0.17513270262204281, "se": 0.13158755835030761}, + "2006.first_treat#2006.year#c.__tr__": {"att": 0.031550170960287845, "se": 0.061041569681307008}, + "2006.first_treat#2007.year#c.__tr__": {"att": 0.080306527133632524, "se": 0.11810866121039022}, + "2007.first_treat#2007.year#c.__tr__": {"att": 0.10774816304839206, "se": 0.12776860159956716}, + "_cons": {"att": 6.2293502307864141, "se": 0.023843276209866905} + } + }, + "20": { + "n_per_cohort": 20, + "G": 80, + "n": 400, + "df_a": 4, + "df_a_initial": 85, + "df_a_nested": 80, + "df_a_redundant": 81, + "rank": 7, + "df_r": 79, + "cells": { + "2004bn.first_treat#2004bn.year#c.__tr__": {"att": -0.050296780093957293, "se": 0.027922218938938603}, + "2004bn.first_treat#2005.year#c.__tr__": {"att": -0.13542856343230192, "se": 0.039492659012540655}, + "2004bn.first_treat#2006.year#c.__tr__": {"att": -0.20008527380675198, "se": 0.046128481856089629}, + "2004bn.first_treat#2007.year#c.__tr__": {"att": -0.18210212485364996, "se": 0.075968261983847152}, + "2006.first_treat#2006.year#c.__tr__": {"att": 0.022421425340464273, "se": 0.038206050766612543}, + "2006.first_treat#2007.year#c.__tr__": {"att": -0.049009832224334804, "se": 0.073878283470593747}, + "2007.first_treat#2007.year#c.__tr__": {"att": -0.046680258002183714, "se": 0.074820594961499759}, + "_cons": {"att": 6.3673294513565919, "se": 0.013633266094166221} + } + }, + "40": { + "n_per_cohort": 40, + "G": 140, + "n": 700, + "df_a": 4, + "df_a_initial": 145, + "df_a_nested": 140, + "df_a_redundant": 141, + "rank": 7, + "df_r": 139, + "cells": { + "2004bn.first_treat#2004bn.year#c.__tr__": {"att": -0.037538524773723235, "se": 0.024361540847770407}, + "2004bn.first_treat#2005.year#c.__tr__": {"att": -0.10207637871169582, "se": 0.033929507546374769}, + "2004bn.first_treat#2006.year#c.__tr__": {"att": -0.1547044358526109, "se": 0.041817458754951305}, + "2004bn.first_treat#2007.year#c.__tr__": {"att": -0.16505539622965759, "se": 0.050869187898973586}, + "2006.first_treat#2006.year#c.__tr__": {"att": -0.0021379792212352033, "se": 0.02752902821801578}, + "2006.first_treat#2007.year#c.__tr__": {"att": -0.085566179995634969, "se": 0.042626013986087855}, + "2007.first_treat#2007.year#c.__tr__": {"att": -0.12886817377512769, "se": 0.043719804770829832}, + "_cons": {"att": 6.1192560211196723, "se": 0.0075272456594156758} + } + }, + "80": { + "n_per_cohort": 80, + "G": 220, + "n": 1100, + "df_a": 4, + "df_a_initial": 225, + "df_a_nested": 220, + "df_a_redundant": 221, + "rank": 7, + "df_r": 219, + "cells": { + "2004bn.first_treat#2004bn.year#c.__tr__": {"att": -0.038064162359172392, "se": 0.023310631552240963}, + "2004bn.first_treat#2005.year#c.__tr__": {"att": -0.09468617423676351, "se": 0.032205191489775339}, + "2004bn.first_treat#2006.year#c.__tr__": {"att": -0.14401918878546274, "se": 0.037373692875051802}, + "2004bn.first_treat#2007.year#c.__tr__": {"att": -0.11133522632938793, "se": 0.039615736809310853}, + "2006.first_treat#2006.year#c.__tr__": {"att": 0.0062590788827515797, "se": 0.021606895089511906}, + "2006.first_treat#2007.year#c.__tr__": {"att": -0.034134199058526972, "se": 0.030170801124917792}, + "2007.first_treat#2007.year#c.__tr__": {"att": -0.042602090415718386, "se": 0.026540554402909169}, + "_cons": {"att": 5.982348942937544, "se": 0.0039836114140290053} + } + }, + "200": { + "n_per_cohort": 200, + "G": 391, + "n": 1955, + "df_a": 4, + "df_a_initial": 396, + "df_a_nested": 391, + "df_a_redundant": 392, + "rank": 7, + "df_r": 390, + "cells": { + "2004bn.first_treat#2004bn.year#c.__tr__": {"att": -0.020982803705344198, "se": 0.022301800464451078}, + "2004bn.first_treat#2005.year#c.__tr__": {"att": -0.075557064505768559, "se": 0.030772264403912838}, + "2004bn.first_treat#2006.year#c.__tr__": {"att": -0.1212832568526639, "se": 0.036222222496439312}, + "2004bn.first_treat#2007.year#c.__tr__": {"att": -0.091115867098532824, "se": 0.034862701796343105}, + "2006.first_treat#2006.year#c.__tr__": {"att": 0.016924854687288118, "se": 0.021080448857745834}, + "2006.first_treat#2007.year#c.__tr__": {"att": -0.025984995955942378, "se": 0.02528976663783444}, + "2007.first_treat#2007.year#c.__tr__": {"att": -0.033501041359260235, "se": 0.019630192588614169}, + "_cons": {"att": 5.8311662331405465, "se": 0.002130009919481389} + } + }, + "500": { + "n_per_cohort": 500, + "G": 500, + "n": 2500, + "df_a": 4, + "df_a_initial": 505, + "df_a_nested": 500, + "df_a_redundant": 501, + "rank": 7, + "df_r": 499, + "cells": { + "2004bn.first_treat#2004bn.year#c.__tr__": {"att": -0.019372363675909551, "se": 0.022381770443295975}, + "2004bn.first_treat#2005.year#c.__tr__": {"att": -0.078319099062053255, "se": 0.030487838520761604}, + "2004bn.first_treat#2006.year#c.__tr__": {"att": -0.13607811444031598, "se": 0.035455486551175253}, + "2004bn.first_treat#2007.year#c.__tr__": {"att": -0.10470747157659592, "se": 0.033874305460323405}, + "2006.first_treat#2006.year#c.__tr__": {"att": 0.0025138619419130628, "se": 0.019932816904001841}, + "2006.first_treat#2007.year#c.__tr__": {"att": -0.039192735591725474, "se": 0.024008748273466048}, + "2007.first_treat#2007.year#c.__tr__": {"att": -0.043106032808698021, "se": 0.018431147223806212}, + "_cons": {"att": 5.7780697609148426, "se": 0.001544041086555641} + } + } + } } } diff --git a/benchmarks/data/fixest_cr1_nonnested_golden.json b/benchmarks/data/fixest_cr1_nonnested_golden.json new file mode 100644 index 000000000..484113e45 --- /dev/null +++ b/benchmarks/data/fixest_cr1_nonnested_golden.json @@ -0,0 +1,61 @@ +{ + "meta": { + "generator": "benchmarks/R/generate_fixest_cr1_nonnested_golden.R", + "r_version": "4.5.2", + "fixest_version": "0.14.2", + "description": "Clustered CR1 with a non-nested absorbed-FE RANK term on a disconnected two-way panel (C=2). crossed_cluster: nothing nested, K.exact=TRUE is the parity anchor (exact span rank 28, df.K=29) vs the default approximate count (df.K=30). unit_cluster: documentation arm -- no fixest ssc reproduces nested-drop + exact-remainder (library K=10 vs default df.K=11; reghdfe agrees with fixest's 11 via its own approximate remainder, see reghdfe_kref_golden.json; the library deviation is pinned as se ratio sqrt((n-10)/(n-11))). connected: control arm where the conventions agree." + }, + "disconnected": { + "data": { + "unit": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19], + "time": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9], + "x": [-0.591031102584369, 0.026594369016167, -1.51655309708187, -1.36265334929581, 1.17848915603162, -0.934151319673365, 1.32360564621202, 0.624917789944714, -0.0457229557060947, -1.00412057567228, -0.828433236663958, -0.348351725172981, -1.53829339708705, -0.255565245464953, -1.14994503270755, 0.0123269676794679, -0.222969540747622, 0.887771647911303, -0.592155279876746, -0.655718118580886, -0.682517622215019, -0.0158581927282744, -0.442604785355045, 0.352557499448816, 0.0731705822914032, 0.0071588003806459, -0.187600110602252, -0.76570064552539, -0.221056820744958, -0.983588587470196, -1.10428404215097, -0.938150214178364, 0.678624244182574, -1.57749786529483, -0.869938458225642, 0.484677045357101, -0.186052698644391, 1.54555470018366, -0.6113800701282, -0.347756487336345, -1.63651630553942, 0.0203814401000371, 0.891742680167841, -0.872749683252251, 0.890050829839478, -0.343874346129348, -2.186781368328, 0.880058181492043, 0.723856556146425, 0.219852680260502, 0.78987057400076, -0.229993900077565, -0.818502482627721, 0.499734157285632, 0.159192348086018, 0.542626428177669, -0.156645054182717, 0.438793320642639, 1.48787059522956, 0.0601650977954321, -0.849012867932317, 2.3396930572073, -0.121202952207389, -1.95020737105763, 0.538711490308542, 1.69351481682805, -0.790968220502282, -1.0752606000877, -0.607875122371882, 0.754401660588849, 0.453476062515215, -0.123433680446261, -0.76309682228376, 0.22827008805461, 1.11946191647044, 0.156573183627664, -0.688772134082379, 0.45294960414968, -1.06754665004997, 0.40156514874452, -0.064773690927254, 0.315496290100479, -0.60568155466456, -0.907584691581957, 2.26160897578211, -0.603226708199307, -1.29786209560729, 0.506451202763822, -0.853334262388665, -1.50603178934849, 1.20232889939251, -1.0278653857247, 0.938269957785018, -0.543154658413624, 0.513095125060458, -0.352590882519352, 1.32653328660507, -1.14015184737515, 1.41310993795933, -0.602178783855351], + "out": [-0.737184170153307, 0.456935957489848, 0.238504603011144, -0.801367388626954, -1.28049889822428, 0.518549925503594, 1.15119424340073, -0.494137399218387, 1.64679181295401, 0.621855433406279, 0.5282671448034, 1.74643378958206, 0.178519143872121, -2.32668805909625, -0.354136172135164, -0.0169009325503089, -0.42906450655316, 1.04876362843334, 0.172548305107774, -0.122734409417013, 2.32126330298461, -0.247917911036566, 1.59008679893854, 1.08744200994463, 0.233967471414275, 2.01481956671627, 1.3950072824351, 0.602622469329669, 1.76845234821618, -1.16569037986164, -0.390744758893371, -1.16427444088891, 1.36872219509093, -1.57118256743879, -0.452289557288066, 0.864185854048759, -0.14228915291334, -0.985526787478794, 0.107777175344013, 0.406722668158443, -1.97438727150224, 1.15437045270228, -0.29263646244231, -0.328834770906268, 0.933829460729813, 0.886622991203856, -0.304361612878549, 3.05557076675472, 0.838259300561206, 0.888588405221773, -0.500980091376268, 1.45607677406328, 0.516472226957554, 1.63946365822013, -0.0111229710586899, 1.61304855863769, 2.31931154142744, 2.70579685000927, 0.477541317594566, 0.448157767278439, 0.381504312023121, 1.60387803840271, 1.12470984846955, -0.0685459267833957, -0.0389932403572389, 1.34909547252119, -1.01417086638884, 0.209456977499621, 0.327208984381329, 1.14217290587027, 1.64255367957849, 1.3706675689923, 1.26156909881477, 1.3257825103295, 0.57032583824924, 1.84080598364095, -0.435050349332937, 2.47603566869651, 1.04814547630723, 2.26278209169537, 0.870095578312731, 1.26451612021621, 1.16676190012435, -0.178536818121551, 4.60692505420721, 0.232563460701728, -0.035230478995931, 1.30242949577884, -0.14091114464718, 0.590231371545861, 1.67093622810471, 0.28824315882779, 2.24062184477509, 0.902695476306297, 1.97164966634418, 3.49508664491642, 3.47749997356729, -0.504241973759709, 1.78084308309533, 1.7821214629558], + "c5": [0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3] + }, + "n_obs": 100, + "unit_cluster": { + "n_clusters": 20, + "coef": 0.541860525668977, + "default": { + "se": 0.138999961325538, + "df_k": 11 + }, + "nested_k_exact": { + "se": 0.154540899077657, + "df_k": 28 + } + }, + "crossed_cluster": { + "n_clusters": 5, + "coef": 0.541860525668977, + "k_exact": { + "se": 0.100145823751911, + "df_k": 29 + }, + "k_exact_nested": { + "se": 0.0994479349445786, + "df_k": 28 + }, + "default": { + "se": 0.10085861441085, + "df_k": 30 + } + } + }, + "connected": { + "data": { + "unit": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19], + "time": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9], + "x": [0.658280963931637, 1.6519984790595, -1.25586749806704, 0.881040288183123, 0.742550527843682, 0.00874156116443356, -0.361468251189953, -0.529600211031989, -1.0428576337799, -1.27349634060697, 0.732037330726618, 0.677262471344541, -0.845015487649319, -0.498274750112178, 1.42302871042578, -0.480506954012273, -0.771268652413693, 0.938922062221907, -1.01471312207972, -0.157114037087788, 0.136448620191636, 0.365970650786999, 0.366380810417199, -0.420786855167701, -0.899572390614362, -0.0253200951818512, 0.0248139759795345, 0.561003817282416, -1.63889606880835, -0.625944959411033, 0.969328438022189, 0.0429000719587969, -2.04013379709133, 0.285952287820425, -1.88282766956995, -0.15071111772361, -0.133243874305348, 1.37539682070027, -1.0176746005728, 0.809377489293562, 0.956915831868467, 0.638131349174454, 1.79197290706804, 1.17179272004281, 1.44011834104833, -1.3062799631651, -0.0247597816198378, 0.490431109972915, -1.68627076518603, -0.829620828933711, -0.238654125576964, -0.0257958838562224, 0.900170667043693, 1.42927873774721, 0.0594457822358079, -0.548887950270233, 0.818184409669585, -0.445460941455792, -0.360500355932039, -0.159547758190255, 0.357706152482205, -2.22376258208559, -0.115372825486446, -1.22403834511744, -1.1951542674148, 0.631337247585679, -1.69665002159048, 0.227301516675197, -0.310209536983869, 1.22453842424653, -0.658199629531311, -0.326082409258126, 2.65076184606774, 0.838799943338105, -0.802962786970869, -0.156717323599945, -1.03630108518638, -1.4153066216368, 0.908108596964723, -0.59922785350864, -0.393417116411639, -1.02156633137411, -0.164512802159358, 0.0310750763469802, 0.000848999425567642, 1.34345582777401, -1.65334228451632, 1.6324982037773, 0.481806535961803, -1.64022351077949, -1.04125034329438, 0.0803375767690832, 0.154902757379373, 0.641111910950676, 0.1274287396582, 1.57046515324945, 0.191960608321051, -0.454277536315632, -0.818878987885141, 0.30114351448628, -1.55005090160424, -0.72519599476199, -0.797741726688203, 2.06550044725063, 2.04925872291846, 0.786209454002233, 0.180870969788242, -0.679377494412951, 1.46378330659275, 0.44470118685535, -1.58363734000052, 0.120832955460222, -0.729606528263807, -0.516810356092618, 1.13520998312328, -0.318695466230325, 0.092775340607271, 0.627523394963977, -0.287241679815022, 0.276924545339165, 0.358783696955915, -0.922000345109948, -0.998950602390343, -0.277454746901607, 0.153750787820796, -1.10651395716324, -0.129230344755356, 1.29274044847617, 0.155493169075694, 0.138991329543313, -0.925818058934539, -0.637474716672803, -2.38451955993373, 0.207068446842626, 0.0879848794802258, 0.669134520790577, -1.86699096741038, -1.52312605001272, 2.21778959004011, -2.21854507482819, -1.53874401269311, 0.633034071175666, 0.411134597421937, -0.583943269741735, 0.942578304685943, -0.53446994539148, -1.4474161300306, -0.679668866316591, -0.0786438655603861, -0.156476923799872, 0.177034979712513, -0.221341777694144, 0.901519800134285, -0.0493141925062218, -0.437957469271032, 0.471701948482876, 0.714834222047201, 0.201292100781617, -0.413182587809343, 0.343263735624697, -0.878764930773286, 1.18811492619173, 1.4922595012095, 0.755229521488336, -1.40053133255127, -0.283978333895583, 0.0106066584837929, 0.54215921011528, -0.466296133932515, -1.75666299917714, 0.543870286797416, 1.36260369840394, -0.930340019048195, -2.17722027627304, 1.28352379632988, -0.356008316210986, 0.974272532073816, -1.50163132469859, -0.370988265767208, -0.221351382294098, -0.910645480229413, -0.199418781216622, 1.22222474121889, 0.565356304051217, 0.176239307145796, 0.855901753470314, 1.49613647214405, 1.1528573153111, -0.0483852988347887, 1.06778374570831, 1.54351782896409, 0.660799968278101, -0.34142841859183, 2.0116464159095, 0.441210221475438, 0.501538713657385, 0.186055522276737, -1.64169132295117, -0.100677077030551, -0.327677331246451], + "out": [1.67217393509495, -0.0138609212756668, 0.954952516409018, -2.36158665413303, 1.61187839280633, 0.62226491981596, 0.44870207674649, -1.02308662871166, -1.77679754345796, 1.43180843035406, -0.661965673869088, 1.50227201525199, 1.16233896132477, 1.14630691661376, 1.18562285408769, -0.94902781001005, -0.683584205319051, -0.0815676888813331, -0.195699785039419, 0.728774483640628, -1.39054263646554, 1.23172509501718, 0.926847417384485, 0.591060056746635, 0.47930261319081, -0.749587292514106, -0.208768884402223, -0.37193541086031, -0.393625178071474, 0.123445906148093, 0.588331203611123, 1.16020676797268, 0.033338027679953, -0.99367559456281, 0.243489452363023, 0.465672576614394, -1.27000074809735, 0.706058089171324, -0.819656105444303, 1.0956072819775, -0.0341408165241097, 1.85950216229798, 2.16042542791799, 0.838304238336711, 2.39279103702044, 0.986117429537003, 0.407820594766349, 1.94902915613975, -1.15165847204151, -0.265304766674087, -1.08649456341694, 2.35730669111378, -0.476222303161506, 2.8932195369101, 1.96794948054524, -0.0775073686003742, 0.760819273369391, -0.200243063204433, -1.52803566643601, -1.43363693632432, 1.41837421956449, -0.627805514932412, 0.504089620060319, -0.0621130924886617, -0.489743520066725, 0.419867070060005, -2.03868597236594, -0.614268878788981, 0.31431341664119, 0.138730841128633, -0.210794715499112, -1.23480491388974, 2.08585858711627, 1.13150346756659, 0.49376470795364, -0.425558160570985, 0.81229713142919, 1.70398456679263, 1.11672616119073, 0.359554597190721, 0.333631882691019, 1.69367148455035, -0.521554828359728, 0.435964010954009, 2.36795327699438, 2.60241309903433, 1.34272054859362, 1.17578540232443, 0.963330661029498, -1.61028644488004, 0.465444636151222, 0.69889743514619, 0.368067598751124, -0.926697545455255, -0.925042806285489, 1.18370060542729, 1.97994092959834, 1.48866458540702, 2.32346474662498, 2.88678476505183, 0.670114671167082, -1.00583755380142, 0.913341428774616, 2.46552511874389, 0.799286151423684, 1.9052863256192, -0.166284888273899, 0.620703168459027, 1.53658131142483, -0.851450149083804, 0.787582820300557, 1.77871499912492, 0.0597033564559818, 0.171375924527051, 1.93829015629417, 0.737358342380028, 2.11588427813494, 1.59864427081103, -0.699097160523057, 1.32615789101931, -1.14579159112779, -0.232637376589728, 3.52636847926857, -0.896495200964187, 1.64985393123085, 2.49749902244923, 0.938317323898118, 3.92091151806142, 0.308703149689999, 0.500893883047175, 0.287670348332621, 0.287746106445999, -0.81816290431473, 0.828068251045895, 1.41935664102232, 1.05138459524522, -0.545626116980175, 1.17219733557648, 0.922494150263573, 0.152556207011097, -0.791304767592319, 2.24630844997329, 0.83906563275975, 0.669036369407718, 1.31787774474273, 0.560530244091883, -0.346937266228471, 1.09437435636816, 2.50971468705926, 2.34143984515735, 0.368446435145024, 0.657976278793875, 2.37547874701935, -1.63864501278171, -0.360758068055504, 1.11411265452308, 0.350255607177908, 3.98536105257243, 2.35212232031869, -0.008756256635848, -0.721938059178798, 2.1847751758733, 1.34459423776449, 1.62161863026222, -0.201680530293451, 0.824751285784786, 1.60669821394087, 2.91532137510409, 0.918352647259831, -0.99160056670138, 1.04474087640803, 1.21068461569046, -0.501512047294814, -1.44953943426246, 2.99030015495577, 0.332983005645697, 3.13058372813184, -0.211344797647042, -0.915979097258319, -0.703987059551341, -0.081942162840058, -0.0456587935848688, 0.373302371235629, 1.6485293386234, 1.61111970591265, 1.45433488453055, 1.94405748612844, 2.31112365724083, 1.21098717594445, 2.47329571955322, 3.43164607721463, 3.54179985780053, 0.759628283831951, 3.45177706507973, 2.29158224108689, 0.965668594621043, 2.46522872769825, -0.290066939613683, 1.31016177929553, 0.320694062628146] + }, + "n_obs": 200, + "n_clusters": 20, + "coef": 0.521987502502566, + "cluster_default": { + "se": 0.052426784826336, + "df_k": 11 + } + } +} diff --git a/benchmarks/data/reghdfe_kref_golden.json b/benchmarks/data/reghdfe_kref_golden.json new file mode 100644 index 000000000..0fe9983c0 --- /dev/null +++ b/benchmarks/data/reghdfe_kref_golden.json @@ -0,0 +1,35 @@ +{ + "meta": { + "generator": "benchmarks/stata/generate_reghdfe_kref_golden.do", + "reghdfe_version": "version 6.13.1 10Jan2026", + "stata_version": 19, + "dgp": "deterministic disconnected two-way panel; see generator header", + "note": "cross_cluster = exact-span parity anchor (K=29, matches library + fixest full/K.exact). unit_cluster = documentation arm: reghdfe counts the nested-remainder approximately (implied K=11) where the library uses the exact remainder (K=10); the test pins se ratio == sqrt((N-10)/(N-11)). No external reference implements nested-drop + exact-remainder." + }, + "unit_cluster": { + "cmd": "reghdfe out x, absorb(unit time) vce(cluster unit)", + "n": 100, + "G": 20, + "df_a": 9, + "df_a_initial": 30, + "df_a_nested": 20, + "df_a_redundant": 21, + "rank": 1, + "df_r": 19, + "coef": 0.5367231638418074, + "se": 0.12545333917294271 + }, + "cross_cluster": { + "cmd": "reghdfe out x, absorb(unit time) vce(cluster c5)", + "n": 100, + "G": 5, + "df_a": 28, + "df_a_initial": 30, + "df_a_nested": 0, + "df_a_redundant": 2, + "rank": 1, + "df_r": 4, + "coef": 0.5367231638418074, + "se": 0.092431570011765093 + } +} diff --git a/benchmarks/stata/README.md b/benchmarks/stata/README.md index 860c4bb3c..50c58be11 100644 --- a/benchmarks/stata/README.md +++ b/benchmarks/stata/README.md @@ -224,18 +224,23 @@ external anchors, with the ETWFE-vs-CS gap recorded rather than asserted away. count, with ATTs matching to ~1e-15 — and warns about the reduction rather than passing it over in silence. `n` and `n_units` are serialized because the row count is the finding, not incidental. -- **ETWFE SEs** — do **not** match `jwdid`. Every cell is uniformly SMALLER - than Stata's, by a factor that shrinks as the cluster count grows: 1.0280 at - G=20, 1.0132 at G=40, 1.00264 at G=191, 1.0010 at G=500. Each arm therefore - pins its OWN measured ratio; the constant does not transfer between arms. - **The mechanism is not identified.** - The gap tracks `sqrt(G/(G-1))` but sits consistently ABOVE it, and the CR1 - factor in `linalg.py` already applies `(G/(G-1)) * ((n-1)/(n-k))` — so a - missing cluster term is ruled out. The test pins the observed ratio and its - within-fit uniformity rather than loosening a tolerance, so a change in - magnitude *or* uniformity fails loudly; it records the gap as MEASURED, not - diagnosed. See the REGISTRY `## WooldridgeDiD (ETWFE)` note and the `TODO.md` - row. +- **ETWFE SEs** — match `jwdid` (reghdfe) at machine precision under the 3.9 + `K_reference` convergence. The historical uniform gap (1.0280 at G=20, + 1.0132 at G=40, 1.00264 at G=191, 1.0010 at G=500) was defect D2: the + library's clustered CR1 factor counted only the visible treatment-cell + columns, omitting the absorbed FE not nested in the unit cluster + (`K_reference = cells + T` on this no-intercept within design). All arms now + gate the SE ratio at 1.0 (`rtol=1e-9`; measured spreads ~1e-14..1e-15). +- **The subsample `ladder` block** — jwdid re-run on deterministic rosters + (the first N units per `first.treat` cohort by ascending `countyreal`, + N ∈ {5, 10, 20, 40, 80, 200, 500} → G ∈ {20, 40, 80, 140, 220, 391, 500}), + each rung storing reghdfe's df accounting (`df_a` + its + initial/nested/redundant decomposition, `rank`, `df_r`) alongside per-cell + `att`/`se`. This gates the few-cluster behavior (where the historical gap + was largest, ~2.8% at G=20) and doubles as the K-accounting probe: the + consuming test asserts `df_a == absorbed_fe_cr1_k_increment − 1` at every + cluster count. See the REGISTRY `## WooldridgeDiD (ETWFE)` hc1 note and + `docs/methodology/variance-conventions.md`. ## Input panel @@ -247,6 +252,35 @@ failure mode that produced the false assertion. Both the generator and the Python test assert the digest, so a swapped panel cannot silently retarget the parity. +--- + +# `reghdfe` anchor for the clustered CR1 `K_reference` on a disconnected panel + +`benchmarks/stata/generate_reghdfe_kref_golden.do` produces +`benchmarks/data/reghdfe_kref_golden.json`, consumed by +`tests/test_variance_conventions.py::TestReghdfeKReferenceParity`. The DGP is +**deterministic (no RNG)** — integer formulas the Python test rebuilds verbatim +— so no data is embedded or shipped: a disconnected two-way panel (units 0-9 +observed in periods 0-4, units 10-19 in periods 5-9; C=2 components, span rank +`U + T − C = 28`). + +Two arms: + +- **`cross_cluster`** (cluster crosses both FE dims → nothing nested): + reghdfe's pairwise dof method computes the exact span rank (`df_a = 28`, + denominator `N − 29`) and the library matches at ~1e-17, agreeing with + fixest `ssc(K.fixef="full", K.exact=TRUE)` (see + `benchmarks/R/generate_fixest_cr1_nonnested_golden.R`). This is the parity + anchor for the exact non-nested RANK term. +- **`unit_cluster`** (unit FE nested in the cluster): a DOCUMENTATION arm — + reghdfe counts the nested-remainder approximately (`df_a = T − 1 = 9`, + implied K = 11; its pairwise correction skips pairs containing a + nested-dropped dim) where the library uses the exact remainder given the + nested span (K = 10). No external reference implements that composition, so + the test pins the deviation exactly: `se_reghdfe / se_library == + sqrt((N−10)/(N−11))`. On CONNECTED designs the two coincide — the jwdid + subsample ladder above pins machine-precision agreement at every G. + ## Known constraints - **Batch mode always exits 0**, even on a hard error (`r(NNN);`). Never trust the diff --git a/benchmarks/stata/generate_etwfe_cs_golden.do b/benchmarks/stata/generate_etwfe_cs_golden.do index 6e1195fbb..4b1c90e22 100644 --- a/benchmarks/stata/generate_etwfe_cs_golden.do +++ b/benchmarks/stata/generate_etwfe_cs_golden.do @@ -44,13 +44,16 @@ *! - Both commands are run WITHOUT covariates. `csdid`'s method is therefore *! immaterial (reg/ipw/dr coincide absent covariates); `method(reg)` is *! passed only to make the no-covariate path explicit. -*! - SEs are emitted alongside the point estimates. The library's ETWFE `hc1` -*! SEs are uniformly SMALLER than `jwdid`'s, by a factor that shrinks as the -*! cluster count grows. The MECHANISM is not identified: the gap tracks -*! sqrt(G/(G-1)) but sits consistently above it, and the CR1 factor at -*! linalg.py already applies (G/(G-1))*((n-1)/(n-k)), so a missing cluster -*! term is ruled out. The consuming test pins the POINT estimates tightly and -*! records the SE ratio as MEASURED, not diagnosed. +*! - SEs are emitted alongside the point estimates. The historical uniform SE +*! gap vs `jwdid` (1.0280 at G=20 shrinking to 1.0010 at G=500) was defect +*! D2 of the 3.9 variance program: the library's clustered CR1 factor used +*! only the visible treatment-cell count, omitting the absorbed FE not +*! nested in the unit cluster. Under the K_reference convergence +*! (docs/methodology/variance-conventions.md) the SEs match reghdfe's at +*! machine precision, and the `ladder` block below pins that agreement at +*! every cluster count alongside reghdfe's own df accounting +*! (df_a / rank / df_r), which the Python side cross-checks against +*! absorbed_fe_cr1_k_increment. *! - The `never` arm is the external anchor for issue #724: `jwdid ... never` *! omits the `g-1` reference cell per cohort (W2025 Eq. 6.1/6.4), which the *! library previously left to QR rank detection. @@ -213,6 +216,54 @@ local ja_N = e(N) local ja_nunits : word count `_at_units' restore +* ------------------------------------------------------------------------------ +* Subsample LADDER: jwdid on nested rosters spanning G ~ 20..500, the artifact +* required by the K_reference convergence work (TODO row "Stata subsample +* ladder"). Roster rule (deterministic, reproduced verbatim by the Python +* side): the first N units per first_treat cohort by ASCENDING countyreal, +* N in {5, 10, 20, 40, 80, 200, 500}; cohort sizes are 309/20/40/131, so the +* rungs land at G = 20/40/80/140/220/391/500. Each rung records reghdfe's own +* df accounting (df_a and its initial/nested/redundant decomposition, rank, +* df_r) so the golden pins not just the SEs but the K-accounting that +* produces them - the Python test cross-checks df_a against +* absorbed_fe_cr1_k_increment per rung. +* ------------------------------------------------------------------------------ +local ladder_rungs "5 10 20 40 80 200 500" +foreach N of local ladder_rungs { + preserve + keep countyreal first_treat + duplicates drop + bysort first_treat (countyreal): gen _idx = _n + quietly keep if _idx <= `N' + keep countyreal + tempfile roster`N' + quietly save `roster`N'' + restore, preserve + quietly merge m:1 countyreal using `roster`N'', keep(match) nogenerate + jwdid lemp, ivar(countyreal) tvar(year) gvar(first_treat) + matrix lad_b_`N' = e(b) + matrix lad_V_`N' = e(V) + local lad_names_`N' : colnames lad_b_`N' + local lad_n_`N' = e(N) + local lad_G_`N' = e(N_clust) + local lad_dfa_`N' = e(df_a) + local lad_dfai_`N' = e(df_a_initial) + local lad_dfan_`N' = e(df_a_nested) + local lad_dfared_`N' = e(df_a_redundant) + local lad_rank_`N' = e(rank) + local lad_dfr_`N' = e(df_r) + restore +} +* Pin the roster rule itself: a changed rule or panel would move these. +assert `lad_G_5' == 20 +assert `lad_G_10' == 40 +assert `lad_G_20' == 80 +assert `lad_G_40' == 140 +assert `lad_G_80' == 220 +assert `lad_G_200' == 391 +assert `lad_G_500' == 500 +assert `lad_n_500' == 2500 + * ------------------------------------------------------------------------------ * Emit the golden. * ------------------------------------------------------------------------------ @@ -323,6 +374,46 @@ forvalues i = 1/`k' { } } file write `fh' _n " }" _n +file write `fh' " }," _n + +* --- subsample ladder (K_reference gate at every cluster count) --------------- +file write `fh' `" "ladder": {"' _n +file write `fh' `" "roster_rule": "first N units per first_treat cohort by ascending countyreal","' _n +file write `fh' `" "jwdid_cmd": "jwdid lemp, ivar(countyreal) tvar(year) gvar(first_treat)","' _n +file write `fh' `" "rungs": {"' +local rsep "" +foreach N of local ladder_rungs { + file write `fh' "`rsep'" _n `" "`N'": {"' _n + file write `fh' `" "n_per_cohort": `N',"' _n + file write `fh' `" "G": `lad_G_`N'',"' _n + file write `fh' `" "n": `lad_n_`N'',"' _n + file write `fh' `" "df_a": `lad_dfa_`N'',"' _n + file write `fh' `" "df_a_initial": `lad_dfai_`N'',"' _n + file write `fh' `" "df_a_nested": `lad_dfan_`N'',"' _n + file write `fh' `" "df_a_redundant": `lad_dfared_`N'',"' _n + file write `fh' `" "rank": `lad_rank_`N'',"' _n + file write `fh' `" "df_r": `lad_dfr_`N'',"' _n + file write `fh' `" "cells": {"' + local csep "" + local k = colsof(lad_b_`N') + forvalues i = 1/`k' { + local nm : word `i' of `lad_names_`N'' + scalar bval = lad_b_`N'[1, `i'] + scalar seval = sqrt(lad_V_`N'[`i', `i']) + if !missing(bval) { + _jnum bval + local a = r(s) + _jnum seval + local s4 = r(s) + file write `fh' "`csep'" _n `" "`nm'": {"att": `a', "se": `s4'}"' + local csep "," + } + } + file write `fh' _n " }" _n + file write `fh' " }" + local rsep "," +} +file write `fh' _n " }" _n file write `fh' " }" _n file write `fh' "}" _n file close `fh' diff --git a/benchmarks/stata/generate_reghdfe_kref_golden.do b/benchmarks/stata/generate_reghdfe_kref_golden.do new file mode 100644 index 000000000..f120bf687 --- /dev/null +++ b/benchmarks/stata/generate_reghdfe_kref_golden.do @@ -0,0 +1,207 @@ +*! Golden generator: reghdfe K_reference anchor on a DISCONNECTED two-way panel. +*! +*! Purpose +*! External anchor for the K_reference clustered-CR1 accounting on an +*! IRREGULAR (disconnected) design, in two arms: +*! +*! cross_cluster -- the PARITY anchor for the exact non-nested RANK term: +*! with nothing nested in the cluster, reghdfe's pairwise dof method +*! computes the exact span rank (df_a = 28 = U + T - C, denominator +*! N - 29), agreeing with fixest ssc(K.fixef="full", K.exact=TRUE) and +*! with the library at machine precision. +*! +*! unit_cluster -- a DOCUMENTATION arm for the nested COMPOSITION: with +*! unit nested-dropped, reghdfe counts the time REMAINDER at the +*! per-dim approximation T - 1 = 9 (its pairwise correction only covers +*! pairs where both dims survive), implying K = rank + df_a + constant +*! = 11. The library's K_reference counts the remainder at its exact +*! rank GIVEN the nested span (28 - 20 = 8, so K = 10) -- the consistent +*! extension of the exact-rank principle both references apply when +*! nothing is nested. NO external reference implements that composition +*! (fixest's K.exact composes incoherently with its nested drop, see +*! benchmarks/R/generate_fixest_cr1_nonnested_golden.R), so this arm pins +*! the DEVIATION exactly: se_reghdfe / se_library == +*! sqrt((N - 10) / (N - 11)), one df, and the df_a decomposition records +*! why. On connected designs (exact == approx remainder) the two agree at +*! machine precision -- the jwdid subsample ladder +*! (etwfe_cs_stata_golden.json) pins that at every cluster count. +*! +*! Design (deterministic, NO RNG -- the Python side rebuilds the frame from +*! the same integer formulas, so no data needs to be embedded or shipped) +*! 100 rows: units 0-9 observed in periods 0-4, units 10-19 in periods 5-9 +*! (two bipartite components, C = 2; span rank of [unit FE, time FE] +*! including constants = 20 + 10 - 2 = 28). +*! unit = floor((_n-1)/5) time = mod(_n-1,5) + 5*(unit>=10) +*! x = (mod(_n*7, 13) - 6)/13 z = (mod(_n*11, 17) - 8)/17 +*! out = 0.5*x + 0.2*mod(unit,3) + 0.1*time + z +*! Every operation is exact integer arithmetic followed by one IEEE +*! division / multiply-add chain evaluated left-to-right in both +*! languages, so the frames agree bit-for-bit. +*! +*! Arms +*! unit_cluster : vce(cluster unit). Unit FE nested in the cluster -> +*! dropped (nested = 20); reghdfe counts the time +*! remainder approximately (df_a = 9 = T - 1, redundant = +*! nested + the global constant only), denominator +*! N - 11. Library: exact remainder 8 -> K_reference = +*! 10 (the documented one-df deviation above). +*! cross_cluster: vce(cluster c5), c5 = mod(unit+time, 5) crossing both +*! dims -> nothing nested, df_a = exact span = 28, +*! denominator N - 29 (agrees with the library AND fixest +*! ssc(K.fixef="full", K.exact=TRUE), df.K = 29). +*! +*! Consuming test +*! tests/test_variance_conventions.py::TestReghdfeKReferenceParity +*! +*! Outputs (checked into the repo) +*! benchmarks/data/reghdfe_kref_golden.json +*! +*! Usage (run from the repo root) +*! /Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ +*! benchmarks/stata/generate_reghdfe_kref_golden.do +*! Then confirm the log is clean: grep -E '^r\([0-9]+\);' generate_reghdfe_kref_golden.log +*! (Stata batch mode ALWAYS exits 0, even on error - never trust the exit code.) +*! +*! SSC dependencies +*! reghdfe (with ftools) + +version 19 +clear all +set more off +set type double + +* JSON number formatter (same as generate_etwfe_cs_golden.do). +capture program drop _jnum +program define _jnum, rclass + args x fmt + if "`fmt'" == "" local fmt "%21.17g" + local s = strtrim(string(`x', "`fmt'")) + if substr("`s'", 1, 1) == "." local s = "0" + "`s'" + else if substr("`s'", 1, 2) == "-." local s = "-0" + substr("`s'", 2, .) + return local s "`s'" +end + +capture which reghdfe +if _rc { + display as error "Missing SSC package reghdfe. Run benchmarks/stata/requirements.do first." + exit 111 +} + +* ------------------------------------------------------------------------------ +* Deterministic disconnected panel (see header). +* ------------------------------------------------------------------------------ +set obs 100 +gen unit = floor((_n - 1) / 5) +gen time = mod(_n - 1, 5) + 5 * (unit >= 10) +gen x = (mod(_n * 7, 13) - 6) / 13 +gen z = (mod(_n * 11, 17) - 8) / 17 +gen out = 0.5 * x + 0.2 * mod(unit, 3) + 0.1 * time + z +gen c5 = mod(unit + time, 5) + +* ------------------------------------------------------------------------------ +* Arm 1: cluster = unit (nested drop + exact remainder rank). +* ------------------------------------------------------------------------------ +reghdfe out x, absorb(unit time) vce(cluster unit) +assert e(df_a) == 9 +assert e(df_a_nested) == 20 +assert e(rank) == 1 +local u_coef = _b[x] +local u_se = _se[x] +local u_dfa = e(df_a) +local u_dfai = e(df_a_initial) +local u_dfan = e(df_a_nested) +local u_dfared = e(df_a_redundant) +local u_rank = e(rank) +local u_dfr = e(df_r) +local u_G = e(N_clust) +local u_N = e(N) + +* ------------------------------------------------------------------------------ +* Arm 2: cluster = c5 (nothing nested; exact span rank). +* ------------------------------------------------------------------------------ +reghdfe out x, absorb(unit time) vce(cluster c5) +assert e(df_a) == 28 +assert e(df_a_nested) == 0 +assert e(df_a_redundant) == 2 +local x_coef = _b[x] +local x_se = _se[x] +local x_dfa = e(df_a) +local x_dfai = e(df_a_initial) +local x_dfan = e(df_a_nested) +local x_dfared = e(df_a_redundant) +local x_rank = e(rank) +local x_dfr = e(df_r) +local x_G = e(N_clust) +local x_N = e(N) + +* ------------------------------------------------------------------------------ +* reghdfe version (drift signal; SSC has no version history). +* ------------------------------------------------------------------------------ +local rv "unknown" +capture findfile reghdfe.ado +if _rc == 0 { + tempname vh + file open `vh' using "`r(fn)'", read text + file read `vh' line + local n = 0 + while r(eof) == 0 & `n' < 15 { + if substr(`"`macval(line)'"', 1, 2) == "*!" { + local body = strtrim(substr(`"`macval(line)'"', 3, .)) + if strpos(lower(`"`body'"'), "version") > 0 { + local rv `"`body'"' + continue, break + } + } + local ++n + file read `vh' line + } + file close `vh' + local rv = subinstr(`"`rv'"', `"""', "'", .) + local rv = subinstr(`"`rv'"', "\", "/", .) +} + +* ------------------------------------------------------------------------------ +* Emit the golden. +* ------------------------------------------------------------------------------ +local sver = c(stata_version) + +tempname fh +file open `fh' using "benchmarks/data/reghdfe_kref_golden.json", write replace text +file write `fh' "{" _n +file write `fh' `" "meta": {"' _n +file write `fh' `" "generator": "benchmarks/stata/generate_reghdfe_kref_golden.do","' _n +file write `fh' `" "reghdfe_version": "`rv'","' _n +file write `fh' `" "stata_version": `sver',"' _n +file write `fh' `" "dgp": "deterministic disconnected two-way panel; see generator header","' _n +file write `fh' `" "note": "cross_cluster = exact-span parity anchor (K=29, matches library + fixest full/K.exact). unit_cluster = documentation arm: reghdfe counts the nested-remainder approximately (implied K=11) where the library uses the exact remainder (K=10); the test pins se ratio == sqrt((N-10)/(N-11)). No external reference implements nested-drop + exact-remainder.""' _n +file write `fh' " }," _n +foreach arm in u x { + if "`arm'" == "u" { + file write `fh' `" "unit_cluster": {"' _n + file write `fh' `" "cmd": "reghdfe out x, absorb(unit time) vce(cluster unit)","' _n + } + else { + file write `fh' `" "cross_cluster": {"' _n + file write `fh' `" "cmd": "reghdfe out x, absorb(unit time) vce(cluster c5)","' _n + } + file write `fh' `" "n": ``arm'_N',"' _n + file write `fh' `" "G": ``arm'_G',"' _n + file write `fh' `" "df_a": ``arm'_dfa',"' _n + file write `fh' `" "df_a_initial": ``arm'_dfai',"' _n + file write `fh' `" "df_a_nested": ``arm'_dfan',"' _n + file write `fh' `" "df_a_redundant": ``arm'_dfared',"' _n + file write `fh' `" "rank": ``arm'_rank',"' _n + file write `fh' `" "df_r": ``arm'_dfr',"' _n + _jnum ``arm'_coef' + local a = r(s) + _jnum ``arm'_se' + local s = r(s) + file write `fh' `" "coef": `a',"' _n + file write `fh' `" "se": `s'"' _n + if "`arm'" == "u" file write `fh' " }," _n + else file write `fh' " }" _n +} +file write `fh' "}" _n +file close `fh' + +display "Wrote benchmarks/data/reghdfe_kref_golden.json" diff --git a/diff_diff/estimators.py b/diff_diff/estimators.py index 272d0d8f2..65c0e9d96 100644 --- a/diff_diff/estimators.py +++ b/diff_diff/estimators.py @@ -29,8 +29,10 @@ from diff_diff.results import DiDResults, MultiPeriodDiDResults, PeriodEffect from diff_diff.utils import ( WildBootstrapResults, + absorbed_fe_cr1_k_increment, absorbed_fe_rank, build_fe_dummy_blocks, + cluster_nested_fe_dims, demean_by_groups, fe_dummy_names, pre_demean_norms, @@ -503,6 +505,10 @@ def fit( has_intercept_col=True, weights=survey_weights, ) + # Stash the raw FE columns: the clustered-CR1 K_reference + # increment needs them AFTER the effective cluster resolves, + # but the in-place demean below overwrites them with floats. + _fe_cols_raw = working_data[list(absorb)].copy() # Method of alternating projections: for N > 1 absorbed dimensions a # single sequential sweep is only exact on balanced (orthogonal-FE) # panels; demean_by_groups iterates to the exact (W)LS-FWL residual. @@ -656,6 +662,44 @@ def fit( _conley_time_arr = None _conley_unit_arr = None + # Clustered-CR1 K_reference adjustment (variance-conventions.md D2/D1): + # absorbed FE not nested in the cluster ADD their conditional rank; + # cluster-nested explicit FE dummies SUBTRACT theirs. Computed only + # for the effective-hc1 clustered analytical lane — under + # wild_bootstrap the analytical fit is deliberately unclustered + # (adjustment travels through the WCB wiring instead), and under a + # survey design the survey variance replaces the CR1 sandwich + # wholesale (moot by design). + _cr1_k_adj = 0 + if ( + _fit_vcov_type == "hc1" + and self.inference != "wild_bootstrap" + and effective_cluster_ids is not None + and resolved_survey is None + ): + if absorbed_vars: + _cr1_k_adj = absorbed_fe_cr1_k_increment( + _fe_cols_raw, + list(absorb), + effective_cluster_ids, + has_intercept_col=True, + weights=survey_weights, + ) + elif fixed_effects: + _nested_fe = cluster_nested_fe_dims( + working_data, + list(fixed_effects), + effective_cluster_ids, + weights=survey_weights, + ) + if _nested_fe: + _cr1_k_adj = -absorbed_fe_rank( + working_data, + _nested_fe, + has_intercept_col=True, + weights=survey_weights, + ) + # Don't forward `robust=self.robust` when the vcov_type has been # remapped; `robust=False + vcov_type="hc1"` would otherwise trip # the conflict check inside `LinearRegression.__init__`. The @@ -677,7 +721,7 @@ def fit( conley_unit=_conley_unit_arr, conley_lag_cutoff=self.conley_lag_cutoff, df_convention=self.df_convention, - ).fit(X, y, df_adjustment=n_absorbed_effects) + ).fit(X, y, df_adjustment=n_absorbed_effects, cluster_k_adjustment=_cr1_k_adj) coefficients = reg.coefficients_ residuals = reg.residuals_ @@ -762,8 +806,35 @@ def _refit_did_absorb(w_r): # test-inversion based; no reference t-distribution, so no # effective inference df). _inference_df_used = None + # K_reference adjustment for the bootstrap's own CR1 factors, + # computed against the RAW cluster ids the bootstrap partitions + # on (NOT effective_cluster_ids — the analytical-lane gate above + # deliberately passed 0 under wild_bootstrap). + _wcb_k_adj = 0 + if absorbed_vars: + _wcb_k_adj = absorbed_fe_cr1_k_increment( + _fe_cols_raw, + list(absorb), + cluster_ids, + has_intercept_col=True, + weights=survey_weights, + ) + elif fixed_effects: + _nested_wcb = cluster_nested_fe_dims( + working_data, + list(fixed_effects), + cluster_ids, + weights=survey_weights, + ) + if _nested_wcb: + _wcb_k_adj = -absorbed_fe_rank( + working_data, + _nested_wcb, + has_intercept_col=True, + weights=survey_weights, + ) se, p_value, conf_int, t_stat, vcov, _ = self._run_wild_bootstrap_inference( - X, y, residuals, cluster_ids, att_idx + X, y, residuals, cluster_ids, att_idx, cluster_k_adjustment=_wcb_k_adj ) else: # Use analytical inference from LinearRegression @@ -869,6 +940,7 @@ def _run_wild_bootstrap_inference( residuals: np.ndarray, cluster_ids: np.ndarray, coefficient_index: int, + cluster_k_adjustment: int = 0, ) -> Tuple[float, float, Tuple[float, float], float, np.ndarray, WildBootstrapResults]: """ Run wild cluster bootstrap inference. @@ -885,6 +957,12 @@ def _run_wild_bootstrap_inference( Cluster identifiers for each observation. coefficient_index : int Index of the coefficient to compute inference for. + cluster_k_adjustment : int, default 0 + Signed K_reference adjustment for the bootstrap's CR1 factors + (nestedness computed by the caller against THESE raw cluster + ids). Applied to both the analytical SE inside + ``wild_bootstrap_se`` and the stored-vcov recompute below, so + ``se == sqrt(vcov[j, j])`` stays exact. Returns ------- @@ -903,6 +981,7 @@ def _run_wild_bootstrap_inference( seed=self.seed, return_distribution=False, p_val_type=self.p_val_type, + cluster_k_adjustment=cluster_k_adjustment, ) self._bootstrap_results = bootstrap_results @@ -927,6 +1006,7 @@ def _run_wild_bootstrap_inference( X, y, cluster_ids=cluster_ids, + cluster_k_adjustment=cluster_k_adjustment, return_vcov=True, rank_deficient_action="silent", ) @@ -1770,6 +1850,10 @@ def fit( # type: ignore[override] has_intercept_col=True, weights=survey_weights, ) + # Stash the raw FE columns for the clustered-CR1 K_reference + # increment (computed after the effective cluster resolves; the + # in-place demean below overwrites these with floats). + _fe_cols_raw_mp = working_data[list(absorb)].copy() # Method of alternating projections (exact for unbalanced panels; a # single sequential sweep is exact only on balanced orthogonal-FE panels). working_data, _ = demean_by_groups( # count superseded by absorbed_fe_rank above @@ -1976,6 +2060,40 @@ def fit( # type: ignore[override] _conley_time_arr = None _conley_unit_arr = None + # Clustered-CR1 K_reference adjustment (variance-conventions.md D1/D2). + # Absorb lane: the absorbed increment (non-nested rank + no intercept + # term — MPD's X carries an intercept). Non-absorb lanes: the built-in + # period dummies ARE MPD's time-FE block (a supplied time FE is + # skipped as redundant above), so they and any placed `_mp_fes` + # dummies SUBTRACT their joint rank when nested in the cluster — + # preserving MPD's absorb/fixed_effects equivalence. Survey lanes + # pass 0 (survey vcov replaces CR1 wholesale). + _cr1_k_adj_mp = 0 + if _fit_vcov_type == "hc1" and effective_cluster_ids is not None and not _use_survey_vcov: + if absorb: + _cr1_k_adj_mp = absorbed_fe_cr1_k_increment( + _fe_cols_raw_mp, + list(absorb), + effective_cluster_ids, + has_intercept_col=True, + weights=survey_weights, + ) + else: + _mp_fe_blocks = list(_mp_fes) if fixed_effects else [] + _nested_mp = cluster_nested_fe_dims( + working_data, + _mp_fe_blocks + [time], + effective_cluster_ids, + weights=survey_weights, + ) + if _nested_mp: + _cr1_k_adj_mp = -absorbed_fe_rank( + working_data, + _nested_mp, + has_intercept_col=True, + weights=survey_weights, + ) + # Note: Wild bootstrap for multi-period effects is complex (multiple coefficients) # For now, we use analytical inference even if inference="wild_bootstrap" coefficients, residuals, fitted, vcov = solve_ols( # type: ignore[call-overload, misc] # mypy gives up on the Optional-arg union explosion ("Not all union combinations were tried") @@ -1988,6 +2106,7 @@ def fit( # type: ignore[override] rank_deficient_action=self.rank_deficient_action, weights=survey_weights, weight_type=survey_weight_type, + cluster_k_adjustment=_cr1_k_adj_mp, vcov_type=_fit_vcov_type, conley_coords=_conley_coords_arr, conley_cutoff_km=self.conley_cutoff_km, diff --git a/diff_diff/imputation.py b/diff_diff/imputation.py index 06a55c226..817fd7310 100644 --- a/diff_diff/imputation.py +++ b/diff_diff/imputation.py @@ -39,6 +39,7 @@ from diff_diff.linalg import solve_ols from diff_diff.utils import ( _iterative_fe_solve, + absorbed_fe_cr1_k_increment, demean_by_groups, pre_demean_norms, safe_inference, @@ -2414,6 +2415,19 @@ def _compute_lead_coefficients( _ols_weights = survey_weights_0 _ols_weight_type = "pweight" if survey_weights_0 is not None else None _use_survey_vcov = resolved_survey_full is not None + # Clustered-CR1 K_reference increment (variance-conventions.md D2): + # the demeaned lead design carries NO intercept column, so the + # absorbed constant contributes the +1 term. Survey fits pass 0 + # (cluster_ids is None there and the survey vcov replaces CR1). + _cr1_k_adj_imp = 0 + if not _use_survey_vcov: + _cr1_k_adj_imp = absorbed_fe_cr1_k_increment( + df_0, + [time, unit], + cluster_ids, + has_intercept_col=False, + weights=survey_weights_0, + ) try: result = solve_ols( X_dm, @@ -2421,6 +2435,7 @@ def _compute_lead_coefficients( weights=_ols_weights, weight_type=_ols_weight_type, cluster_ids=None if _use_survey_vcov else cluster_ids, + cluster_k_adjustment=_cr1_k_adj_imp, return_vcov=True, rank_deficient_action=self.rank_deficient_action, column_names=all_x_cols, diff --git a/diff_diff/linalg.py b/diff_diff/linalg.py index e541ff2e1..651d43f74 100644 --- a/diff_diff/linalg.py +++ b/diff_diff/linalg.py @@ -362,10 +362,12 @@ def _absorbed_fe_vcov_scale(n_eff: float, k_eff: int, df_adjustment: int) -> flo (per the non-finite-df fail-closed contract) rather than leaving a misleading ``k_visible`` SE in place. - Callers must gate on non-clustered ``classical``/``hc1``: clustered SEs - follow fixest's ``ssc`` nested-FE convention (FE nested in the cluster are - not counted, so ``k_visible`` already matches for the nested case) and - ``hc2``/``hc2_bm`` use leverage / Satterthwaite DOF -- none must be rescaled. + Callers must gate on non-clustered ``classical``/``hc1``: the CLUSTERED + CR1 factor has its own K_reference accounting via the + ``cluster_k_adjustment`` seam (fixest's ``ssc`` nested-FE convention drops + FE nested in the cluster and counts non-nested FE — see + ``docs/methodology/variance-conventions.md`` D2), and ``hc2``/``hc2_bm`` + use leverage / Satterthwaite DOF — neither must be rescaled here. """ denom_visible = n_eff - k_eff denom_full = n_eff - k_eff - df_adjustment @@ -376,6 +378,90 @@ def _absorbed_fe_vcov_scale(n_eff: float, k_eff: int, df_adjustment: int) -> flo return denom_visible / denom_full +class InvalidClusterKAdjustment(ValueError): + """Raised when ``cluster_k_adjustment`` is passed where it cannot apply. + + A dedicated subclass so broad ``except ValueError`` handlers around + clustered solves (e.g. LPDiD's unclustered-fallback wrapper) can re-raise + it instead of silently degrading a misuse to NaN inference. Deliberately + NOT exported from ``diff_diff.__init__`` (module-private precedent: + ``balancing.BalanceError``). + """ + + +def _validate_cluster_k_adjustment_type(cluster_k_adjustment: int) -> None: + """Type half of the K_reference contract, checked BEFORE any zero fast + path or truthiness gate: ``0.0 == 0`` and ``False == 0`` are truthy + comparisons, so a value-first check would silently accept non-int zeros + on short-circuited routes (classical fits, survey fits, degenerate + wild-bootstrap returns) and make the contract route-dependent.""" + if isinstance(cluster_k_adjustment, bool) or not isinstance( + cluster_k_adjustment, (int, np.integer) + ): + raise InvalidClusterKAdjustment( + f"cluster_k_adjustment must be an int, got {type(cluster_k_adjustment).__name__}" + ) + + +def _validate_cluster_k_adjustment( + cluster_k_adjustment: int, + cluster_ids: Optional[np.ndarray], + vcov_type: str, +) -> None: + """Front-door contract for the clustered-CR1 K_reference adjustment. + + A nonzero adjustment is meaningful only on the clustered ``hc1`` lane + (the CR1 finite-sample factor). Raises :class:`InvalidClusterKAdjustment` + on any other combination so misuse fails loudly on EVERY route — + including ``return_vcov=False`` calls and Rust dispatches that never + reach the numpy kernel. + """ + _validate_cluster_k_adjustment_type(cluster_k_adjustment) + if cluster_k_adjustment == 0: + return + if cluster_ids is None: + raise InvalidClusterKAdjustment( + "cluster_k_adjustment is nonzero but cluster_ids is None; the " + "K_reference adjustment applies only to the clustered CR1 factor." + ) + if vcov_type != "hc1": + raise InvalidClusterKAdjustment( + f"cluster_k_adjustment is nonzero but vcov_type is '{vcov_type}'; " + "the K_reference adjustment applies only to the clustered hc1 " + "(CR1) family." + ) + + +def _rescale_cr1_k_adjustment( + vcov: Optional[np.ndarray], + n_eff: float, + k: int, + cluster_k_adjustment: int, +) -> Optional[np.ndarray]: + """Exact scalar map from a ``k``-denominator clustered CR1 vcov to the + ``k + cluster_k_adjustment`` one. + + The CR1 factor ``(G/(G-1)) * ((n-1)/(n-k))`` is a pure scalar on the + finished vcov, so the corrected matrix is exactly + ``vcov * (n_eff - k) / (n_eff - k - adj)`` — used to correct the Rust + lanes without touching the Rust kernels. None-preserving + (``return_vcov=False`` puts ``vcov=None`` in the tail slot). Fails closed + to all-NaN when the corrected residual dof ``n_eff - k - adj`` or the + corrected count ``k + adj`` is non-positive; the VISIBLE-side saturation + (``n_eff <= k``) never reaches this helper — the Rust kernels already + return all-NaN there, matching the numpy kernel's kept visible-k guard. + """ + if vcov is None: + return None + if cluster_k_adjustment == 0: + return vcov + denom_old = n_eff - k + denom_new = n_eff - k - cluster_k_adjustment + if denom_new <= 0 or (k + cluster_k_adjustment) <= 0: + return np.full_like(vcov, np.nan) + return vcov * (denom_old / denom_new) + + def _expand_vcov_with_nan( vcov_reduced: np.ndarray, k_full: int, @@ -919,6 +1005,7 @@ def solve_ols( rank_deficient_action: str = ..., column_names: Optional[List[str]] = ..., skip_rank_check: bool = ..., + cluster_k_adjustment: int = ..., weights: Optional[np.ndarray] = ..., weight_type: str = ..., vcov_type: str = ..., @@ -945,6 +1032,7 @@ def solve_ols( rank_deficient_action: str = ..., column_names: Optional[List[str]] = ..., skip_rank_check: bool = ..., + cluster_k_adjustment: int = ..., weights: Optional[np.ndarray] = ..., weight_type: str = ..., vcov_type: str = ..., @@ -971,6 +1059,7 @@ def solve_ols( rank_deficient_action: str = ..., column_names: Optional[List[str]] = ..., skip_rank_check: bool = ..., + cluster_k_adjustment: int = ..., weights: Optional[np.ndarray] = ..., weight_type: str = ..., vcov_type: str = ..., @@ -999,6 +1088,7 @@ def solve_ols( rank_deficient_action: str = "warn", column_names: Optional[List[str]] = None, skip_rank_check: bool = False, + cluster_k_adjustment: int = 0, weights: Optional[np.ndarray] = None, weight_type: str = "pweight", vcov_type: str = "hc1", @@ -1051,6 +1141,18 @@ def solve_ols( rank-deficient matrices. Use only when you know the design matrix is full rank. If the matrix is actually rank-deficient, results may be incorrect (minimum-norm solution instead of R-style NA handling). + cluster_k_adjustment : int, default 0, keyword-only + Signed K_reference adjustment added to the visible column count in + the CLUSTERED CR1 finite-sample factor only: absorbed FE not nested + in the cluster add their conditional rank; cluster-nested explicit + FE dummies subtract theirs (the reghdfe / fixest + ``ssc(K.fixef="nested")`` convention; see + ``docs/methodology/variance-conventions.md`` defect D2). Validated + at this front door on EVERY route — a nonzero value with + ``cluster_ids=None`` or ``vcov_type != "hc1"`` raises + :class:`InvalidClusterKAdjustment` even on ``return_vcov=False`` + calls. Never changes coefficients, residuals, or the reported + tail df. weights : ndarray of shape (n,), optional Observation weights for Weighted Least Squares. When provided, minimizes sum(w_i * (y_i - X_i @ beta)^2). Weights should be @@ -1223,6 +1325,11 @@ def solve_ols( if cluster_ids is not None or weights is not None: _validate_vcov_args(vcov_type, cluster_ids, weights) + # Front-door K_reference-adjustment validation (same rationale as the + # conley guard above): runs BEFORE routing/backend branching so + # `return_vcov=False` calls and Rust dispatches cannot bypass it. + _validate_cluster_k_adjustment(cluster_k_adjustment, cluster_ids, vcov_type) + # WLS transformation: apply sqrt(w) scaling to X and y # This happens BEFORE routing to Rust or NumPy backends — they receive # pre-transformed X_w, y_w and solve standard OLS. @@ -1277,6 +1384,10 @@ def solve_ols( ) if result is not None and diagnostics_out is not None: diagnostics_out["solve_ols_fastpath"] = "chol_rust" + if result is not None and cluster_k_adjustment: + result = result[:-1] + ( + _rescale_cr1_k_adjustment(result[-1], float(n), k, cluster_k_adjustment), + ) if result is None and ( HAS_RUST_BACKEND and _rust_solve_ols is not None @@ -1303,6 +1414,13 @@ def solve_ols( stacklevel=2, ) result = None # Force Python fallback below + # K_reference rescale AFTER the sentinel check: a Rust all-NaN is + # the documented rank sentinel (passes through); the rescale's own + # fail-closed NaN must not be produced before that check runs. + if result is not None and cluster_k_adjustment: + result = result[:-1] + ( + _rescale_cr1_k_adjustment(result[-1], float(n), k, cluster_k_adjustment), + ) if result is None: result = _solve_ols_numpy( X, @@ -1315,6 +1433,7 @@ def solve_ols( _skip_rank_check=True, _fastpath=fastpath, _diagnostics_out=diagnostics_out, + cluster_k_adjustment=cluster_k_adjustment, vcov_type=vcov_type, conley_coords=conley_coords, conley_cutoff_km=conley_cutoff_km, @@ -1361,6 +1480,10 @@ def solve_ols( ) if result is not None and diagnostics_out is not None: diagnostics_out["solve_ols_fastpath"] = "chol_rust" + if result is not None and cluster_k_adjustment: + result = result[:-1] + ( + _rescale_cr1_k_adjustment(result[-1], float(n), k, cluster_k_adjustment), + ) if result is None and ( HAS_RUST_BACKEND and _rust_solve_ols is not None @@ -1390,6 +1513,13 @@ def solve_ols( stacklevel=2, ) result = None # Force Python fallback below + # K_reference rescale AFTER the non-finite check above — the + # rescale's own fail-closed all-NaN would otherwise be misread as + # backend instability (spurious warning + pointless numpy re-run). + if result is not None and cluster_k_adjustment: + result = result[:-1] + ( + _rescale_cr1_k_adjustment(result[-1], float(n), k, cluster_k_adjustment), + ) if result is None: result = _solve_ols_numpy( @@ -1404,6 +1534,7 @@ def solve_ols( _fastpath=fastpath, _cert_info=cert_out, _diagnostics_out=diagnostics_out, + cluster_k_adjustment=cluster_k_adjustment, vcov_type=vcov_type, conley_coords=conley_coords, conley_cutoff_km=conley_cutoff_km, @@ -1445,6 +1576,7 @@ def solve_ols( cluster_ids, weights=weights, weight_type=weight_type, + cluster_k_adjustment=cluster_k_adjustment, vcov_type=vcov_type, conley_coords=conley_coords, conley_cutoff_km=conley_cutoff_km, @@ -1464,6 +1596,7 @@ def solve_ols( cluster_ids, weights=weights, weight_type=weight_type, + cluster_k_adjustment=cluster_k_adjustment, vcov_type=vcov_type, conley_coords=conley_coords, conley_cutoff_km=conley_cutoff_km, @@ -1494,6 +1627,7 @@ def _solve_ols_numpy( column_names: Optional[List[str]] = ..., _precomputed_rank_info: Optional[Tuple[int, np.ndarray, np.ndarray]] = ..., _skip_rank_check: bool = ..., + cluster_k_adjustment: int = ..., _fastpath: bool = ..., _cert_info: Optional[dict] = ..., _diagnostics_out: Optional[dict] = ..., @@ -1520,6 +1654,7 @@ def _solve_ols_numpy( column_names: Optional[List[str]] = ..., _precomputed_rank_info: Optional[Tuple[int, np.ndarray, np.ndarray]] = ..., _skip_rank_check: bool = ..., + cluster_k_adjustment: int = ..., _fastpath: bool = ..., _cert_info: Optional[dict] = ..., _diagnostics_out: Optional[dict] = ..., @@ -1546,6 +1681,7 @@ def _solve_ols_numpy( column_names: Optional[List[str]] = ..., _precomputed_rank_info: Optional[Tuple[int, np.ndarray, np.ndarray]] = ..., _skip_rank_check: bool = ..., + cluster_k_adjustment: int = ..., _fastpath: bool = ..., _cert_info: Optional[dict] = ..., _diagnostics_out: Optional[dict] = ..., @@ -1574,6 +1710,7 @@ def _solve_ols_numpy( column_names: Optional[List[str]] = None, _precomputed_rank_info: Optional[Tuple[int, np.ndarray, np.ndarray]] = None, _skip_rank_check: bool = False, + cluster_k_adjustment: int = 0, _fastpath: bool = False, _cert_info: Optional[dict] = None, _diagnostics_out: Optional[dict] = None, @@ -1716,6 +1853,7 @@ def _solve_ols_numpy( X_reduced, residuals, cluster_ids, + cluster_k_adjustment=cluster_k_adjustment, vcov_type=vcov_type, conley_coords=conley_coords, conley_cutoff_km=conley_cutoff_km, @@ -1757,6 +1895,7 @@ def _solve_ols_numpy( X, residuals, cluster_ids, + cluster_k_adjustment=cluster_k_adjustment, vcov_type=vcov_type, conley_coords=conley_coords, conley_cutoff_km=conley_cutoff_km, @@ -1781,6 +1920,7 @@ def _validate_vcov_args( vcov_type: str, cluster_ids: Optional[np.ndarray], weights: Optional[np.ndarray], + cluster_k_adjustment: int = 0, ) -> None: """Shared validation for ``vcov_type`` / ``cluster_ids`` / ``weights`` combinations. @@ -1813,6 +1953,9 @@ def _validate_vcov_args( raise ValueError( f"vcov_type must be one of {sorted(_VALID_VCOV_TYPES)}; " f"got {vcov_type!r}" ) + # Mirrored K_reference-adjustment contract for direct compute_robust_vcov + # / kernel callers (solve_ols routes enforce it at its own front door). + _validate_cluster_k_adjustment(cluster_k_adjustment, cluster_ids, vcov_type) if vcov_type in ("classical", "hc2") and cluster_ids is not None: msg = { "classical": ( @@ -1967,6 +2110,7 @@ def compute_robust_vcov( vcov_type: str = "hc1", return_dof: bool = False, *, + cluster_k_adjustment: int = 0, conley_coords: Optional[np.ndarray] = None, conley_cutoff_km: Optional[float] = None, conley_metric: ConleyMetric = "haversine", @@ -2065,6 +2209,14 @@ def compute_robust_vcov( array of per-coefficient degrees of freedom. For ``classical``, ``hc1``, ``hc2``: every element is ``n_eff - k``. For ``hc2_bm`` one-way: Imbens-Kolesar (2016) Satterthwaite DOF per contrast. + cluster_k_adjustment : int, default 0, keyword-only + Signed K_reference adjustment added to the visible column count in + the CLUSTERED CR1 finite-sample factor only (absorbed FE not nested + in the cluster add rank; cluster-nested explicit FE dummies + subtract; see ``docs/methodology/variance-conventions.md`` D2). + Nonzero values require ``cluster_ids`` and ``vcov_type="hc1"`` + (raises :class:`InvalidClusterKAdjustment` otherwise). Never affects + the reported ``dof_vec`` (tail df). Returns ------- @@ -2083,7 +2235,9 @@ def compute_robust_vcov( For cluster-robust (CR1, Liang-Zeger): meat = sum_g (X_g' u_g)(X_g' u_g)' - adjustment = (G / (G-1)) * ((n-1) / (n-k)) + adjustment = (G / (G-1)) * ((n-1) / (n - (k + cluster_k_adjustment))) + (k = reduced design rank; the signed adjustment carries the + K_reference absorbed-FE accounting, default 0) For HC2 one-way (weighted per review MEDIUM #3): h_ii = w_i * x_i' * (X'WX)^{-1} * x_i (unweighted: w_i = 1) @@ -2097,7 +2251,7 @@ def compute_robust_vcov( The cluster-robust CR1 computation is vectorized using pandas groupby. """ - _validate_vcov_args(vcov_type, cluster_ids, weights) + _validate_vcov_args(vcov_type, cluster_ids, weights, cluster_k_adjustment) # Validate weights before dispatching to backend if weights is not None: @@ -2178,7 +2332,12 @@ def compute_robust_vcov( cluster_ids_int = pd.factorize(cluster_ids)[0].astype(np.int64) try: - return _rust_compute_robust_vcov(X, residuals, cluster_ids_int) + _rust_vcov = _rust_compute_robust_vcov(X, residuals, cluster_ids_int) + if cluster_k_adjustment: + _rust_vcov = _rescale_cr1_k_adjustment( + _rust_vcov, float(X.shape[0]), X.shape[1], cluster_k_adjustment + ) + return _rust_vcov except ValueError as e: # Translate Rust errors to consistent Python error messages or fallback error_msg = str(e) @@ -2204,6 +2363,7 @@ def compute_robust_vcov( weight_type=weight_type, vcov_type=vcov_type, return_dof=return_dof, + cluster_k_adjustment=cluster_k_adjustment, conley_coords=conley_coords, conley_cutoff_km=conley_cutoff_km, conley_metric=conley_metric, @@ -2223,6 +2383,7 @@ def compute_robust_vcov( weight_type=weight_type, vcov_type=vcov_type, return_dof=return_dof, + cluster_k_adjustment=cluster_k_adjustment, conley_coords=conley_coords, conley_cutoff_km=conley_cutoff_km, conley_metric=conley_metric, @@ -3223,6 +3384,7 @@ def _compute_robust_vcov_numpy( vcov_type: str = "hc1", return_dof: bool = False, *, + cluster_k_adjustment: int = 0, conley_coords: Optional[np.ndarray] = None, conley_cutoff_km: Optional[float] = None, conley_metric: ConleyMetric = "haversine", @@ -3248,7 +3410,7 @@ def _compute_robust_vcov_numpy( # directly and previously bypassed the raise, letting unsupported # combinations (cluster + classical, cluster + hc2, cluster + weights + # hc2_bm) silently produce wrong inference. Reviewer P0 fix. - _validate_vcov_args(vcov_type, cluster_ids, weights) + _validate_vcov_args(vcov_type, cluster_ids, weights, cluster_k_adjustment) n, k = X.shape @@ -3499,7 +3661,17 @@ def _compute_robust_vcov_numpy( # NaN vcov so downstream inference is degenerate (NaN) rather than raising # ZeroDivisionError — consistent with the library's all-or-nothing NaN # convention for undefined inference. - if n_eff - k <= 0: + # + # The clustered K_reference count k_inf = k + cluster_k_adjustment adds + # two more fail-closed sides (D2 correction): the VISIBLE-k side is KEPT + # (saturation is a property of residual dof — at n_eff == k the residuals + # are identically zero and the sandwich is 0 under ANY denominator, so a + # positive corrected denominator is not a licence to compute), plus + # n_eff - k_inf <= 0 and k_inf <= 0 (the latter covers a negative + # adjustment exceeding the reduced design width, where the factor would + # silently DEFLATE the SE). + _k_inf = k + cluster_k_adjustment + if n_eff - k <= 0 or (cluster_k_adjustment and (n_eff - _k_inf <= 0 or _k_inf <= 0)): nan_vcov = np.full((k, k), np.nan) if return_dof: return nan_vcov, np.full(k, np.nan, dtype=np.float64) @@ -3541,8 +3713,13 @@ def _compute_robust_vcov_numpy( if n_clusters < 2: raise ValueError(f"Need at least 2 clusters for cluster-robust SEs, got {n_clusters}") - # Small-sample adjustment - adjustment = (n_clusters / (n_clusters - 1)) * ((n_eff - 1) / (n_eff - k)) + # Small-sample adjustment. k_inf = k + cluster_k_adjustment is the + # K_reference count: visible columns plus absorbed FE not nested in + # the cluster, minus cluster-nested explicit FE dummies (reghdfe / + # fixest ssc(K.fixef="nested") convention; see + # docs/methodology/variance-conventions.md D2). The reported tail df + # (dof_vec below) deliberately stays on the visible k. + adjustment = (n_clusters / (n_clusters - 1)) * ((n_eff - 1) / (n_eff - _k_inf)) # Sum scores within each cluster using pandas groupby (vectorized) cluster_scores = pd.DataFrame(scores).groupby(cluster_ids).sum().values @@ -4341,6 +4518,7 @@ def fit( *, cluster_ids: Optional[np.ndarray] = None, df_adjustment: int = 0, + cluster_k_adjustment: int = 0, ) -> "LinearRegression": """ Fit OLS regression. @@ -4357,6 +4535,18 @@ def fit( df_adjustment : int, default 0 Additional degrees of freedom adjustment (e.g., for absorbed fixed effects). The effective df will be n - k - df_adjustment. + cluster_k_adjustment : int, default 0 + Signed K_reference adjustment for the CLUSTERED CR1 finite-sample + factor only (absorbed FE not nested in the cluster add rank; + cluster-nested explicit FE dummies subtract; see + ``docs/methodology/variance-conventions.md`` D2). Unlike the + non-negative ``df_adjustment``, this value is signed, never moves + the reported residual ``df_``, and is validated fail-loud: a + nonzero value on a fit that resolves to classical/unclustered + raises ``InvalidClusterKAdjustment``. Under a survey design that + computes survey variance the kwarg is INERT by design — the + survey vcov replaces the CR1 sandwich wholesale, so estimator + callers pass 0 on those lanes. Returns ------- @@ -4482,6 +4672,20 @@ def fit( _effective_survey_design, effective_cluster_ids ) + # Fail-loud K_reference contract at FIT level, BEFORE the + # classical-vs-robust branch (the classical branch below never + # reaches solve_ols's front door). The TYPE half is unconditional — + # classical and survey fits must reject a non-int too, else the + # contract is route-dependent. Survey-vcov fits are the deliberate + # VALUE exception: the survey variance replaces the CR1 sandwich + # wholesale, so a (well-typed) adjustment is inert by design there + # (documented above) and estimator callers pass 0 on those lanes. + _validate_cluster_k_adjustment_type(cluster_k_adjustment) + if cluster_k_adjustment and not _use_survey_vcov: + _validate_cluster_k_adjustment( + cluster_k_adjustment, effective_cluster_ids, _fit_vcov_type + ) + if _fit_vcov_type != "classical" or effective_cluster_ids is not None: # Use solve_ols with robust/cluster SEs. # When survey vcov will be used, skip standard vcov computation. @@ -4500,6 +4704,7 @@ def fit( rank_deficient_action=self.rank_deficient_action, weights=_fit_weights, weight_type=_fit_weight_type, + cluster_k_adjustment=(0 if _use_survey_vcov else cluster_k_adjustment), vcov_type=_fit_vcov_type, conley_coords=self.conley_coords, conley_cutoff_km=self.conley_cutoff_km, diff --git a/diff_diff/lpdid.py b/diff_diff/lpdid.py index 2bec51adc..4afd0dbe5 100644 --- a/diff_diff/lpdid.py +++ b/diff_diff/lpdid.py @@ -4,9 +4,9 @@ import numpy as np import pandas as pd -from diff_diff.linalg import _rank_guarded_inv, solve_ols +from diff_diff.linalg import InvalidClusterKAdjustment, _rank_guarded_inv, solve_ols from diff_diff.lpdid_results import LPDiDResults -from diff_diff.utils import safe_inference +from diff_diff.utils import absorbed_fe_rank, cluster_nested_fe_dims, safe_inference __all__ = ["LPDiD", "LPDiDResults"] @@ -859,17 +859,41 @@ def _estimate_sample( use_cluster_vcov = len(pd.unique(cluster_ids)) >= 2 vcov = None if use_cluster_vcov: + # Clustered-CR1 K_reference adjustment (variance-conventions.md + # D1 family): LPDiD's built-in `_event_time` dummies ARE its time + # FE (the authors' reference recipe absorbs time FE), and inline + # absorb dummies are user FE — both SUBTRACT their joint rank + # when nested in the cluster. The nested test runs on the same + # `sample` rows/cluster the solve uses. + _lp_fe_dims = [] + if include_time_fe: + _lp_fe_dims.append("_event_time") + _lp_fe_dims.extend(absorb_columns) + _cr1_k_adj_lp = 0 + if _lp_fe_dims: + _nested_lp = cluster_nested_fe_dims( + sample, _lp_fe_dims, cluster_ids, weights=weights + ) + if _nested_lp: + _cr1_k_adj_lp = -absorbed_fe_rank( + sample, _nested_lp, has_intercept_col=True, weights=weights + ) try: coef, _, vcov = solve_ols( design, response, cluster_ids=cluster_ids, + cluster_k_adjustment=_cr1_k_adj_lp, return_vcov=True, rank_deficient_action=self.rank_deficient_action, column_names=column_names, weights=weights, ) - except (ValueError, ZeroDivisionError): + except (ValueError, ZeroDivisionError) as _lp_exc: + # NEVER swallow a K-adjustment contract violation into a + # silent unclustered se=NaN refit (no-silent-failure rule). + if isinstance(_lp_exc, InvalidClusterKAdjustment): + raise coef, _, _ = solve_ols( design, response, diff --git a/diff_diff/sun_abraham.py b/diff_diff/sun_abraham.py index a199647a4..b646431c1 100644 --- a/diff_diff/sun_abraham.py +++ b/diff_diff/sun_abraham.py @@ -24,6 +24,7 @@ from diff_diff.results import _format_survey_block, _get_significance_stars from diff_diff.results_base import BaseResults from diff_diff.utils import ( + absorbed_fe_cr1_k_increment, absorbed_fe_rank, pre_demean_norms, safe_inference, @@ -1567,6 +1568,8 @@ def _fit_saturated_regression( # within-transform) at atol=1e-8 — fixest is the natural R parity # anchor for SA's HC1 default. use_full_dummy = vcov_type in ("hc2", "hc2_bm", "classical") + # Nonzero only on the clustered-hc1 within-transform branch below. + _cr1_k_adj_sa = 0 if use_full_dummy: # Full-dummy auto-route: build [intercept, interactions, @@ -1664,6 +1667,20 @@ def _fit_saturated_regression( has_intercept_col=False, weights=survey_weights, ) + # Clustered-CR1 K_reference increment (variance-conventions.md + # D2): the saturated design has NO intercept column, so the + # absorbed constant contributes the +1 term. Computed against the + # SAME cluster array the LinearRegression below uses (SA + # auto-clusters at unit by default; explicit cluster= otherwise); + # survey designs replace the CR1 sandwich, so they pass 0. + if vcov_type == "hc1" and cluster_ids is not None and resolved_survey is None: + _cr1_k_adj_sa = absorbed_fe_cr1_k_increment( + df, + [unit, time], + cluster_ids, + has_intercept_col=False, + weights=survey_weights, + ) # Interactions occupy columns 0..n_interactions-1 (no intercept) coef_offset = 0 @@ -1700,7 +1717,7 @@ def _fit_saturated_regression( conley_time=_cl_time, conley_unit=_cl_unit, conley_lag_cutoff=self.conley_lag_cutoff, - ).fit(X, y, df_adjustment=df_adj) + ).fit(X, y, df_adjustment=df_adj, cluster_k_adjustment=_cr1_k_adj_sa) vcov = reg.vcov_ diff --git a/diff_diff/twfe.py b/diff_diff/twfe.py index 84f42024b..fce0ba520 100644 --- a/diff_diff/twfe.py +++ b/diff_diff/twfe.py @@ -16,8 +16,10 @@ from diff_diff.linalg import LinearRegression from diff_diff.results import DiDResults from diff_diff.utils import ( + absorbed_fe_cr1_k_increment, absorbed_fe_rank, build_fe_dummy_blocks, + cluster_nested_fe_dims, fe_dummy_names, pre_demean_norms, snap_absorbed_regressors, @@ -529,6 +531,28 @@ def fit( # type: ignore[override] survey_cluster_ids if self.inference != "wild_bootstrap" else None ) + # Clustered-CR1 K_reference increment for the within-transform design + # (variance-conventions.md D2): absorbed unit/time FE not nested in + # the cluster add their conditional rank. Gated on the effective-hc1 + # clustered analytical lane — under wild_bootstrap the constructor + # cluster above is None (the WCB wiring carries its own adjustment), + # the full-dummy branch has no analytical CR1 (hc2/hc2_bm only), and + # survey designs replace the CR1 sandwich wholesale. + _cr1_k_adj_twfe = 0 + if ( + not use_full_dummy + and _fit_vcov_type == "hc1" + and _conley_cluster_override is not None + and resolved_survey is None + ): + _cr1_k_adj_twfe = absorbed_fe_cr1_k_increment( + data, + [unit, time], + _conley_cluster_override, + has_intercept_col=True, + weights=survey_weights, + ) + if self.rank_deficient_action == "error": reg = LinearRegression( include_intercept=False, @@ -547,7 +571,7 @@ def fit( # type: ignore[override] conley_unit=_conley_unit_arr, conley_lag_cutoff=self.conley_lag_cutoff, df_convention=self.df_convention, - ).fit(X, y, df_adjustment=df_adjustment) + ).fit(X, y, df_adjustment=df_adjustment, cluster_k_adjustment=_cr1_k_adj_twfe) else: # Suppress generic warning, TWFE provides context-specific messages below with warnings.catch_warnings(): @@ -569,7 +593,7 @@ def fit( # type: ignore[override] conley_unit=_conley_unit_arr, conley_lag_cutoff=self.conley_lag_cutoff, df_convention=self.df_convention, - ).fit(X, y, df_adjustment=df_adjustment) + ).fit(X, y, df_adjustment=df_adjustment, cluster_k_adjustment=_cr1_k_adj_twfe) coefficients = reg.coefficients_ residuals = reg.residuals_ @@ -693,8 +717,32 @@ def _refit_twfe(w_r): # test-inversion based; no reference t-distribution, so no # effective inference df). _inference_df_used = None + # K_reference adjustment for the bootstrap's own CR1 factors, + # against the RAW cluster ids it partitions on. Within design: + # the absorbed increment. Full-dummy design (hc2/hc2_bm + WCB): + # the explicit unit/time dummy blocks SUBTRACT their nested rank + # (D1 convention) — the analytical hc2 vcov has no CR1 factor, + # but the bootstrap computes its own. + _wcb_k_adj_twfe = 0 + if cluster_ids is not None: + if use_full_dummy: + _nested_wcb = cluster_nested_fe_dims( + data, [unit, time], cluster_ids, weights=survey_weights + ) + if _nested_wcb: + _wcb_k_adj_twfe = -absorbed_fe_rank( + data, _nested_wcb, has_intercept_col=True, weights=survey_weights + ) + else: + _wcb_k_adj_twfe = absorbed_fe_cr1_k_increment( + data, + [unit, time], + cluster_ids, + has_intercept_col=True, + weights=survey_weights, + ) se, p_value, conf_int, t_stat, vcov, _ = self._run_wild_bootstrap_inference( - X, y, residuals, cluster_ids, att_idx + X, y, residuals, cluster_ids, att_idx, cluster_k_adjustment=_wcb_k_adj_twfe ) else: # Use analytical inference from LinearRegression diff --git a/diff_diff/utils.py b/diff_diff/utils.py index 46ee594a8..6d1c013a1 100644 --- a/diff_diff/utils.py +++ b/diff_diff/utils.py @@ -23,6 +23,7 @@ _rust_sc_weight_fw_with_convergence, _rust_sdid_unit_weights, ) +from diff_diff.linalg import _validate_cluster_k_adjustment_type from diff_diff.linalg import compute_robust_vcov as _compute_robust_vcov_linalg from diff_diff.linalg import solve_ols as _solve_ols_linalg @@ -709,6 +710,8 @@ def wild_bootstrap_se( seed: Optional[int] = None, return_distribution: bool = False, p_val_type: str = "two-tailed", + *, + cluster_k_adjustment: int = 0, ) -> WildBootstrapResults: """ Compute wild cluster bootstrap standard errors and p-values. @@ -732,6 +735,15 @@ def wild_bootstrap_se( the original estimate — the studentized bootstrap drives the p-value and CI, not a re-scaled bootstrap dispersion. + ``cluster_k_adjustment`` (keyword-only, default 0) applies the signed + K_reference count to BOTH the analytical CR1 SE and the bootstrap ``corr`` + factor (see ``docs/methodology/variance-conventions.md`` D1/D2), keeping + ``se``, ``t_stat_original`` (= effect/se) and any returned + ``bootstrap_distribution`` on the corrected scale. The constant cancels in + the studentized statistic, so p-values are invariant and CI endpoints move + only within the test-inversion bisection tolerance (~1e-10) — never assert + exact CI equality across adjustment values. + Parameters ---------- X : np.ndarray @@ -811,7 +823,12 @@ def wild_bootstrap_se( MacKinnon, J. G., & Webb, M. D. (2018). The wild bootstrap for few (treated) clusters. The Econometrics Journal, 21(2), 114-135. """ - # Validate inputs + # Validate inputs. The adjustment TYPE is checked at entry — before the + # dropped-coefficient / saturation early returns below — so a non-int + # cannot be silently absorbed into a degenerate NaN result on those + # routes (the value/family contract is enforced by solve_ols itself at + # the clustered vcov call). + _validate_cluster_k_adjustment_type(cluster_k_adjustment) valid_weight_types = ["rademacher", "webb", "mammen"] if weight_type not in valid_weight_types: raise ValueError(f"weight_type must be one of {valid_weight_types}, got '{weight_type}'") @@ -878,11 +895,18 @@ def _degenerate() -> WildBootstrapResults: X_eff = X[:, kept] j_eff = int(np.sum(kept[:coefficient_index])) # position of the coef among kept columns k_eff = X_eff.shape[1] - if n <= k_eff: # no residual degrees of freedom -> CR1 undefined + # Fail closed on ALL THREE saturation sides, matching the kernel: the + # visible n <= k_eff (residuals identically zero), the K_reference side + # n <= k_eff + adj, and a negative adjustment exceeding the reduced width + # (k_eff + adj <= 0, where the corr factor below would silently DEFLATE). + _k_inf_eff = k_eff + cluster_k_adjustment + if n <= k_eff or n <= _k_inf_eff or _k_inf_eff <= 0: return _degenerate() # Now the cluster-robust (CR1) vcov is well-defined; it studentizes the test. - _, _, vcov_original = _solve_ols_linalg(X, y, cluster_ids=cluster_ids, return_vcov=True) + _, _, vcov_original = _solve_ols_linalg( + X, y, cluster_ids=cluster_ids, return_vcov=True, cluster_k_adjustment=cluster_k_adjustment + ) if vcov_original is None: return _degenerate() se_a = float(np.sqrt(vcov_original[coefficient_index, coefficient_index])) @@ -913,10 +937,11 @@ def _degenerate() -> WildBootstrapResults: m_y = y - fit_y_red m_xj = xj - fit_xj_red - # CR1 small-sample correction. NOTE: this constant cancels in |t*| vs |t0| - # (it scales se* and se_a identically), so it affects only the reported SE, - # not the p-value or CI. Kept for fidelity with the analytical CR1 SE. - corr = (n_clusters / (n_clusters - 1)) * ((n - 1) / (n - k_eff)) + # CR1 small-sample correction on the K_reference count (k_eff + adj). + # NOTE: this constant cancels in |t*| vs |t0| (it scales se* and se_a + # identically), so it affects only the reported SE, not the p-value or + # CI. Kept for fidelity with the analytical CR1 SE. + corr = (n_clusters / (n_clusters - 1)) * ((n - 1) / (n - _k_inf_eff)) # Cluster membership: indicator matrix C (G, n) for fast per-cluster score sums. cluster_pos = {c: i for i, c in enumerate(unique_clusters)} @@ -2978,16 +3003,37 @@ def absorbed_fe_rank( if not group_vars: return 0 + codes_list, levels, _ = _factorize_fe_codes(data, group_vars, weights) + rank = _fe_rank_from_codes(codes_list, levels) + if rank == 0: + return 0 + return rank - 1 if has_intercept_col else rank + + +def _factorize_fe_codes( + data: pd.DataFrame, + group_vars: List[str], + weights: Optional[np.ndarray], +) -> Tuple[List[np.ndarray], List[int], Optional[np.ndarray]]: + """Shared factorization for the absorbed-FE helpers. + + Applies the positive-weight row filter (REGISTRY zero-weight-padding + guarantee), validates against NaN group keys, and returns per-dim + integer codes with level counts plus the positive-row mask (None when no + filtering applied) so callers can align companion arrays (cluster ids). + """ frame = data + mask: Optional[np.ndarray] = None if weights is not None: w = np.asarray(weights, dtype=np.float64) if w.shape[0] != len(data): raise ValueError(f"weights length ({w.shape[0]}) must match data rows ({len(data)})") positive = w > 0 if not positive.any(): - return 0 + return [np.empty(0, dtype=np.intp) for _ in group_vars], [0] * len(group_vars), positive if not positive.all(): frame = data.loc[positive] + mask = positive codes_list = [] for g in group_vars: @@ -3001,11 +3047,21 @@ def absorbed_fe_rank( ) codes_list.append(codes) levels = [int(c.max()) + 1 if c.size else 0 for c in codes_list] + return codes_list, levels, mask + + +def _fe_rank_from_codes(codes_list: List[np.ndarray], levels: List[int]) -> int: + """Rank of the FE dummy span INCLUDING the shared constant. + + N == 2: ``sum(levels) - C`` with C the connected components of the + bipartite level graph (Abowd-Creecy-Kramarz); N == 1: ``levels``; + N >= 3: ``sum(levels) - N + 1`` (the documented approximate form). + """ total_levels = sum(levels) if total_levels == 0: return 0 - n_dims = len(group_vars) + n_dims = len(codes_list) if n_dims == 2: from scipy.sparse import coo_matrix from scipy.sparse.csgraph import connected_components @@ -3022,11 +3078,128 @@ def absorbed_fe_rank( # Weak connectivity of the directed one-way bipartite graph equals # undirected connectivity, and skips materializing A + A.T (~2.4x). n_components = int(connected_components(adjacency, directed=True, connection="weak")[0]) - rank = total_levels - n_components - else: - rank = total_levels - n_dims + 1 + return total_levels - n_components + return total_levels - n_dims + 1 - return rank - 1 if has_intercept_col else rank + +def _factorize_cluster_codes( + cluster_ids: np.ndarray, + n_rows: int, + mask: Optional[np.ndarray], +) -> np.ndarray: + """Factorize the cluster array over the same positive-weight rows.""" + cids = np.asarray(cluster_ids) + if cids.shape[0] != n_rows: + raise ValueError(f"cluster_ids length ({cids.shape[0]}) must match data rows ({n_rows})") + if mask is not None: + cids = cids[mask] + ccodes = pd.factorize(cids, sort=False)[0] + if ccodes.size and ccodes.min() < 0: + # "missing values" is the established cluster-validation contract + # phrase (tests match on it); keep it in the message. + raise ValueError( + "cluster_ids contain missing values (NaN keys); drop or impute " + "those rows before fitting." + ) + return ccodes + + +def _nested_dims_from_codes( + group_vars: List[str], + codes_list: List[np.ndarray], + levels: List[int], + ccodes: np.ndarray, +) -> List[str]: + """Dims whose every level maps to exactly one cluster id (O(n) scatter).""" + nested = [] + for g, codes, n_levels in zip(group_vars, codes_list, levels): + if n_levels == 0: + continue + lo = np.full(n_levels, np.iinfo(np.int64).max, dtype=np.int64) + hi = np.full(n_levels, -1, dtype=np.int64) + np.minimum.at(lo, codes, ccodes) + np.maximum.at(hi, codes, ccodes) + present = hi >= 0 + if np.array_equal(lo[present], hi[present]): + nested.append(g) + return nested + + +def cluster_nested_fe_dims( + data: pd.DataFrame, + group_vars: List[str], + cluster_ids: np.ndarray, + *, + weights: Optional[np.ndarray] = None, +) -> List[str]: + """Absorbed FE dims nested in the cluster partition. + + A dim ``g`` is nested iff every level of ``g`` maps to exactly one + cluster id — the reghdfe/fixest ``ssc(K.fixef = "nested")`` condition + under which the dim's parameters are dropped from the clustered CR1 + count. Evaluated over POSITIVE-weight rows only, mirroring + :func:`absorbed_fe_rank`. A coarser cluster (e.g. unit FE clustered by + state) is nested; a finer or crossing cluster is not. + """ + if not group_vars: + return [] + codes_list, levels, mask = _factorize_fe_codes(data, group_vars, weights) + if sum(levels) == 0: + return [] + ccodes = _factorize_cluster_codes(cluster_ids, len(data), mask) + return _nested_dims_from_codes(group_vars, codes_list, levels, ccodes) + + +def absorbed_fe_cr1_k_increment( + data: pd.DataFrame, + group_vars: List[str], + cluster_ids: np.ndarray, + *, + has_intercept_col: bool, + weights: Optional[np.ndarray] = None, +) -> int: + """The clustered-CR1 ``k`` increment for absorbed fixed effects. + + Implements the K_reference bracket (see ``docs/methodology/ + variance-conventions.md`` defect D2):: + + increment = (0 if has_intercept_col else 1) + + rank(all absorbed dims, incl. constant) + - max(rank(cluster-nested dims, incl. constant), 1) + + The ``+1`` term is the rank of the ABSORBED constant, so it exists only + when an FE set is absorbed; with no absorbed dims — or an EMPTY effective + support (all rows zero-weight) — the increment is 0. The ``max(..., 1)`` + floor is load-bearing for the zero-nested-dim case: with no absorbed dim + nested in the cluster the conditioning set is the shared constant alone + (rank 1). Verified against Stata reghdfe 3.2.9 (via jwdid, ~1e-15) and + R fixest 0.14.2 (~1e-12); see the REGISTRY absorbed-FE notes. + + Both rank terms inherit :func:`absorbed_fe_rank`'s dimensionality + contract: exact (component-aware) for one and two dims; for THREE or + more dims the documented D3 ``sum(levels) - N + 1`` approximation + applies — exact for independent connected dims, an over-count for + duplicated/nested triples (REGISTRY absorbed-FE note; TODO.md N-way + row), so the increment is approximate on such designs. + """ + if not group_vars: + return 0 + codes_list, levels, mask = _factorize_fe_codes(data, group_vars, weights) + if sum(levels) == 0: + # Empty effective support: nothing absorbed, increment 0 (guarded so + # the max(..., 1) floor cannot drive the formula to -1). + return 0 + ccodes = _factorize_cluster_codes(cluster_ids, len(data), mask) + nested_names = _nested_dims_from_codes(group_vars, codes_list, levels, ccodes) + full_rank = _fe_rank_from_codes(codes_list, levels) + if nested_names: + nested_idx = [i for i, g in enumerate(group_vars) if g in nested_names] + nested_rank = _fe_rank_from_codes( + [codes_list[i] for i in nested_idx], [levels[i] for i in nested_idx] + ) + else: + nested_rank = 0 + return (0 if has_intercept_col else 1) + full_rank - max(nested_rank, 1) def demean_by_groups( diff --git a/diff_diff/wooldridge.py b/diff_diff/wooldridge.py index d0d383255..c58081740 100644 --- a/diff_diff/wooldridge.py +++ b/diff_diff/wooldridge.py @@ -22,6 +22,7 @@ from diff_diff.linalg import compute_robust_vcov, solve_logit, solve_ols, solve_poisson from diff_diff.utils import ( + absorbed_fe_cr1_k_increment, pre_demean_norms, safe_inference, snap_absorbed_regressors, @@ -2244,6 +2245,29 @@ def _fit_ols( else: _cl_coords = _cl_time = _cl_unit = None + # Clustered-CR1 K_reference increment (variance-conventions.md D2 — + # the original jwdid SE-gap finding): the within-transform design has + # NO intercept column (coef_offset = 0), so the absorbed constant + # contributes the +1 term and the delta is T on the unit-clustered + # arms. Nestedness runs against the SAME resolved cluster array + # passed below. The cohort_trends full-dummy branch keeps 0 (its + # documented L1-convention opt-in), as do survey fits (TSL replaces + # CR1) and one-way/conley families. + _cr1_k_adj_w = 0 + if ( + not use_full_dummy + and self.vcov_type == "hc1" + and cluster_ids is not None + and resolved is None + ): + _cr1_k_adj_w = absorbed_fe_cr1_k_increment( + sample, + [unit, time], + cluster_ids, + has_intercept_col=False, + weights=survey_weights, + ) + coefs, resids, vcov = solve_ols( # type: ignore[call-overload] # mypy union-explosion limitation at this many Optional args X, y, @@ -2253,6 +2277,7 @@ def _fit_ols( column_names=col_names if not use_full_dummy else None, weights=survey_weights, weight_type=survey_weight_type, + cluster_k_adjustment=_cr1_k_adj_w, vcov_type=self.vcov_type, conley_coords=_cl_coords, conley_cutoff_km=self.conley_cutoff_km, diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 6676c60e1..6165f3100 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -135,8 +135,13 @@ inference and the inherited `p_val_type` is inert there.) located by outward bracketing + bisection on the (monotone, step) rejection frequency. The CI is therefore exactly consistent with the p-value (`0 ∈ CI ⟺ p ≥ alpha`) and may be asymmetric. 6. The reported `se` is `se_a` (analytical CR1); `p_val_type ∈ {"two-tailed" (default), - "equal-tailed"}`. CR1 uses the standard `(G/(G−1))((N−1)/(N−k))` correction, which cancels in - `|t*|` vs `|t₀|` so it affects only the reported SE, not the p-value or CI. + "equal-tailed"}`. CR1 uses the standard `(G/(G−1))((N−1)/(N−K_reference))` correction — + under absorbed/nested FE the `k` carries the same signed `cluster_k_adjustment` as the + analytical lane (3.9 K_reference convergence) — which cancels in `|t*|` vs `|t₀|` so it + affects only the reported SE (and `t_stat_original` / any requested + `bootstrap_distribution`, which share `se_a`'s scale), not the p-value; CI endpoints move + only within the bisection tolerance through `se_a`'s bracketing role + (`tests/test_variance_conventions.py::TestKReferenceConvergence`). *Verification — R parity:* validated against `fwildclusterboot::boottest()` defaults on a fixed few-cluster golden (`benchmarks/R/generate_wild_cluster_boot_golden.R` → @@ -289,6 +294,19 @@ where V is the VCV sub-matrix for post-treatment δ_e coefficients. its fit site — MultiPeriodDiD is a second implementation, not an alias of DifferenceInDifferences) — see the TwoWayFixedEffects section's absorbed-FE degrees-of-freedom note for the shared convention +- **Note (clustered CR1 k — built-in period dummies, 3.9 K_reference):** MPD's + built-in `period_*` dummies ARE its time-FE block — a time column supplied via + `fixed_effects=` is skipped as redundant with them, and `absorb=[..., time]` is + the absorbed spelling of the same block — so they follow the nested convention: + under a time-like cluster their rank is subtracted from the clustered CR1 + factor's `k` (negative `cluster_k_adjustment`), and under `absorb=` the + non-nested increment applies, preserving MPD's absorb/fixed_effects clustered-SE + equivalence (the D1 rule; audit rows `mpd_absorb_hc1_cluster_unit` / + `mpd_fixed_effects_hc1_cluster_unit` / `mpd_plain_hc1_cluster_time` in + `tests/test_variance_conventions.py`). The drop is an equivalence/reference + rule, not a blanket nested-dummies rule: a built-in dummy block subtracts + exactly when it is the estimator's FE block with an absorbed spelling or a + reference absorbing it. *Edge cases:* - Reference period: omitted from design matrix; coefficient is zero by construction. @@ -416,16 +434,34 @@ This matches the behavior of R's `fixest::feols()` with absorbed FE. Previously the within-transform SE used `k_visible`, sitting ~6.5% below fixest even though the t-`df` already used `K_full` (an internal inconsistency). Applies to `TwoWayFixedEffects(vcov_type="classical")`, `DifferenceInDifferences(absorb=..., vcov_type in {classical,hc1})`, - and `MultiPeriodDiD(absorb=..., vcov_type in {classical,hc1})`. **Clustered** SEs are unchanged - (this fix is gated on `cluster_ids is None`): the clustered CR1 `k_visible` scale matches fixest - for absorbed FE **nested** in the cluster (e.g. unit FE with unit clustering, per fixest's `ssc` - nested-FE convention, which does not count nested FE). **Known limitation (deviation from fixest):** - when an absorbed FE is *not* nested in the cluster (e.g. `absorb=["unit","time"]` clustered by - `unit`, where the time FE are non-nested), fixest counts the non-nested FE in the CR1 - finite-sample denominator, but the current clustered path uses only `k_visible` -- a small, - pre-existing deviation left out of this D4 (non-clustered) scope and tracked in `DEFERRED.md`. This - is distinct from the SunAbraham / Wooldridge `hc1` deviations below (whose event-study / - aggregation paths auto-cluster or use a different k-convention). `hc2`/`hc2_bm` use leverage / + and `MultiPeriodDiD(absorb=..., vcov_type in {classical,hc1})`. **Clustered** CR1 SEs converge + on the same accounting (3.9, K_reference): the finite-sample factor's `k` is + `K_reference = k_visible + (1 for the absorbed constant when X carries no intercept column) + + rank(absorbed FE not nested in the cluster | nested ones)`, computed by + `absorbed_fe_cr1_k_increment` and threaded as the signed `cluster_k_adjustment` seam through + `solve_ols`/`compute_robust_vcov`/`LinearRegression.fit`/`wild_bootstrap_se`. This reproduces + Stata reghdfe 3.2.9 (via `jwdid`) at machine precision on the committed mpdta arms and the + G≈20..500 subsample ladder (`benchmarks/data/etwfe_cs_stata_golden.json`), and fixest's + nested-FE convention (`ssc(K.fixef = "nested")`) on connected designs + (`tests/test_fixest_did_twfe_parity.py`, tightened to ~1e-9). Nested-in-cluster FE are + dropped exactly as both references do. **Note (deviation from R, exact non-nested rank):** + the non-nested term is the exact RANK given the nested set for one- and two-dimensional + absorbed sets, extending the D3/K.exact principle above; for THREE or more absorbed dims + both ranks inherit the documented D3 `sum(levels) − N + 1` approximation (exact for + independent connected dims, an over-count for duplicated/nested triples — same limitation + note as D3, tracked in `TODO.md`; pinned by + `tests/test_variance_conventions.py::TestClusterKIncrement::test_three_dim_increment_inherits_the_documented_d3_approximation`). + On connected designs the exact two-way rank equals the references' per-dim count; on + DISCONNECTED designs it does not, and NO external reference implements the exact + composition — fixest's `ssc(K.exact=TRUE)` composes incoherently with its nested drop, and + reghdfe's pairwise correction skips pairs containing a nested-dropped dim, so both land one + df higher per extra component. The deviation is pinned as an exact one-df SE ratio against + BOTH references (`benchmarks/data/fixest_cr1_nonnested_golden.json`, + `benchmarks/data/reghdfe_kref_golden.json`; + `tests/test_variance_conventions.py::TestFixestCr1NonNestedParity` / + `TestReghdfeKReferenceParity`), while the crossed-cluster (nothing-nested) arm matches both + references at machine precision. This is distinct from the SunAbraham / Wooldridge `hc1` + tail-df items below (PR C scope). `hc2`/`hc2_bm` use leverage / Satterthwaite DOF and are unaffected. - **Note (deviation from R — clustered CR1 inference df):** under unit clustering, Python's CR1 t-statistics/p-values/CIs use the **residual df** from the fitted design @@ -437,7 +473,8 @@ This matches the behavior of R's `fixest::feols()` with absorbed FE. numerically zero, which is what `tests/test_methodology_twfe.py::test_pvalue_matches_r_twfe` pins; the convention itself is locked at a distinguishable moderate |t|~1.8 by `tests/test_methodology_twfe.py::test_moderate_t_pins_residual_df_convention`, where the - t(148) and t(49) tails differ by ~5%). The SE itself carries the separate ~0.25% non-nested-FE ssc band above. **Opt-in knob (2026-07):** + t(148) and t(49) tails differ by ~5%). The SE itself matches fixest under the K_reference + convergence above (the historical ~0.25% non-nested-FE band is closed). **Opt-in knob (2026-07):** `df_convention="cluster"` on `DifferenceInDifferences` / `TwoWayFixedEffects` / `MultiPeriodDiD` (and the `LinearRegression` linalg surface) switches clustered analytical t/p/CI to the Stata/fixest `G − 1` convention — fallback-level only (survey df and per-coefficient @@ -450,21 +487,27 @@ This matches the behavior of R's `fixest::feols()` with absorbed FE. `tests/test_estimators_vcov_type.py::TestDfConvention` (G−1 tail match, precedence, no-op default). The full-dummy (`fixed_effects=`) idiom carries `df_adjustment == 0` and is unchanged — its residual t-df already matched fixest's full-K count because every FE column sits in `k_visible`. **That claim is scoped to - the residual df and the non-clustered variance families; it does NOT extend to the - clustered CR1 finite-sample factor** — see the deviation note below. + the residual df and the non-clustered variance families; the clustered CR1 + finite-sample factor instead applies the negative nested-rank adjustment** — see + the deviation note below. - **Deviation from R (clustered CR1, full-dummy `fixed_effects=` path):** with - `fixed_effects=[unit, time]` and `cluster=unit`, the CR1 factor's `k` counts every - FE dummy including the cluster-NESTED unit FE (measured `k = 66` on the audit - panel), which fixest's default `ssc(fixef.K = "nested")` drops. The full-dummy - clustered SE is therefore CONSERVATIVE relative to fixest — the opposite direction - from the `absorb=` path's anti-conservative `k_visible` (measured `k = 2` on the - same model). The two documented-equivalent idioms consequently return clustered - SEs differing by 10.35% on an identical model with an identical ATT - (`sqrt((360−2)/(360−66))`, pinned by - `tests/test_variance_conventions.py::test_d1_divergence_is_pinned`), and NEITHER - matches the reghdfe/fixest nested convention. Full measured map: - `docs/methodology/variance-conventions.md` (defect D1); convergence of both paths - on `K_reference` is scheduled as PR B of the 3.9 variance-consolidation program. + `fixed_effects=[unit, time]` and `cluster=unit`, the CR1 factor subtracts the + rank of the cluster-NESTED FE dummies from the design rank (audit panel: + `k = 66 → K_reference = 7`), so the two documented-equivalent idioms return + IDENTICAL clustered SEs + (`tests/test_variance_conventions.py::test_d1_convergence_is_pinned`; the + historical 10.35% split, `sqrt((360−2)/(360−66))`, is closed). This is a + deliberate deviation from a LITERAL R comparison in the other direction: + fixest's `ssc(K.fixef = "nested")` drops nested ABSORBED FE only, not nested + explicit dummy REGRESSORS — `feols(y ~ ... + i(unit) + i(time), cluster = ~unit)` + and `lm()` + clubSandwich CR1S both count all 66 columns where K_reference drops + the nested 59. No external reference implements K_reference on explicit dummies, + so this side is anchored by (a) the `absorb=` side's external parity (Stata + jwdid/reghdfe + fixest arms) and (b) the internal equality gate — two + documented-equivalent spellings of the same model must return the same SE (the + D1 equivalence, which also holds under collinearity drops and on disconnected + designs because the kernel `k` is the design RANK). Full map: + `docs/methodology/variance-conventions.md` (D1/D2, fixed 3.9). *Edge cases:* - Singleton units/periods are automatically dropped @@ -1570,15 +1613,15 @@ where weights ŵ_{g,e} = n_{g,e} / Σ_g n_{g,e} (sample share of cohort g at eve closed. The channel clears under bootstrap (percentile p/CIs used no df) - deliberately narrower than the `event_study_vcov` clear, which also fires under replicate refits whose rows DID use a genuine df. -- **Deviation from R (HC1 finite-sample correction):** SA's - within-transform HC1 SE differs from `fixest::sunab(cluster=~unit)` - by ~1-2% on typical panel sizes. fixest's correction counts the - absorbed unit + time FE in the effective parameter count - (`n / (n - k_total)`) whereas SA's `solve_ols` counts only the - within-transformed design columns (`n / (n - k_dm)`). The IW - aggregation step is otherwise identical. Tracked as a follow-up - (harmonizing the correction or documenting it as an intentional - difference). +- **HC1 finite-sample correction (K_reference convergence, 3.9):** SA's + clustered CR1 factor counts the absorbed unit + time FE via the + `cluster_k_adjustment` increment (`K_reference = k_dm + 1 + rank(non-nested + FE | nested)`), matching `fixest::sunab(cluster=~unit)` — measured gap + ~5e-15 on the R-parity fixture (`tests/test_estimators_vcov_type.py`, the + sunab gate). Historically `solve_ols` counted only the within-transformed + design columns (`n / (n - k_dm)`), sitting ~1-2% below fixest on typical + panel sizes; that follow-up (DEFERRED row 53) is closed. The IW + aggregation step is unchanged. - The within-transform residual df subtracts the absorbed unit + time FE via the component-aware rank (`diff_diff.utils.absorbed_fe_rank`, no-intercept form — the saturated design carries no intercept column) — see the @@ -2063,7 +2106,7 @@ where `g(·)` is the link inverse (logistic or exp), `η_i` is the individual li - **Note:** QMLE sandwich uses `weight_type="aweight"` which applies `(G/(G-1)) * ((n-1)/(n-k))` small-sample adjustment. Stata `jwdid` uses `G/(G-1)` only. The `(n-1)/(n-k)` term is conservative (inflates SEs slightly). For typical ETWFE panels where n >> k, the difference is negligible. *Variance families (`vcov_type`, OLS path only):* -- `hc1` (default) — CR1 Liang-Zeger cluster-robust on the within-transformed design. Bit-equal to prior behavior (FWL preserves the score). The natural R anchor is `fixest::feols(y ~ | unit + time, cluster=~unit)` or Stata `jwdid` (both within-transform). **Deviation from Stata `jwdid` (measured 2026-07-26, `tests/test_etwfe_cs_stata_parity.py`):** the ATT(g,t) POINT estimates match `jwdid` exactly (~3e-8 on the `mpdta` panel, i.e. Stata's log-output rounding), but every `hc1` SE is SMALLER than `jwdid`'s by a factor that is **uniform across cells** (spread < 1e-6 within a fit) and shrinks as the cluster count grows: 1.0280 at G=20, 1.0132 at G=40, 1.0086 at G=60, 1.0046 at G=110, 1.00264 at G=191, 1.0010 at G=500. **Two ratios are PINNED by CI** (`tests/test_etwfe_cs_stata_parity.py`): G=500 on the full `mpdta` panel, and G=191 on the all-eventually-treated arm — each arm measures its own, because the constant does not transfer between cluster counts. The remaining smaller-G figures were measured ad hoc on subsampled panels during the #724 investigation and no committed artifact reproduces them, so treat them as the shape of the trend rather than as regression-gated constants. Committing that subsample ladder as a golden block is part of the derivation work tracked in `TODO.md` — the ladder is the instrument the derivation needs, not a separate chore. The library is therefore systematically **anti-conservative** relative to the reference - negligibly with many clusters (~0.1% at G=500) and materially in few-cluster designs (~2.8% at G=20). **The mechanism IS now derived (2026-07, `docs/methodology/variance-conventions.md` defect D2):** the clustered CR1 factor's `k` counts only the visible treatment-cell columns, omitting the absorbed FE not nested in the cluster; the closed form `K_reference = explicit columns + (1 if no intercept column) + rank(non-nested FE | nested)` reproduces `jwdid`/reghdfe to ~1e-15 on all three committed arms and retrodicts the ad-hoc G=20/G=40 rungs (predicted 1.028016 vs recorded 1.0280; 1.013210 vs 1.0132). The correction itself is deliberately NOT yet applied — it lands with the rest of the clustered-CR1 convergence (PR B of the 3.9 variance-consolidation program), together with the committed subsample ladder, so every affected surface moves once under one ledger row rather than piecemeal. This is distinct from the `lm + clubSandwich` deviation below, whose factor is `k`-based. `CallawaySantAnna` shows no such gap - its SEs match Stata `csdid` outright - which localizes this to the ETWFE path rather than a library-wide convention. **Deviation from R `lm + clubSandwich::vcovCR(type="CR1S")`:** the full-dummy `lm` SE differs by a factor of `sqrt((n - k_within) / (n - k_total))` because clubSandwich's `(n-1)/(n-p)` finite-sample correction counts ALL columns (intercept + treatment + unit dummies + time dummies = `k_total`) while WooldridgeDiD's `solve_ols` on the within-transformed design counts only the treatment-cell columns (`k_within`). On the 240-obs / 51-column R-parity fixture this is ~11%; on typical larger panels (n >> k_total) the gap shrinks to <2%. No public WooldridgeDiD code path exposes the `lm + CR1S` (CR1 cluster-robust on the full-dummy design) finite-sample correction — `vcov_type="hc2_bm"` routes to the CR2 Bell-McCaffrey sandwich on the full-dummy design (different variance estimator entirely), not CR1S. Users who need exact `lm + clubSandwich::vcovCR(type="CR1S")` parity must call `solve_ols` directly on a full-dummy design or fit via R. Same deviation pattern as SunAbraham PR #472 (`fixest::sunab` vs `lm + clubSandwich`). +- `hc1` (default) — CR1 Liang-Zeger cluster-robust on the within-transformed design. Point estimates, residuals, and the unscaled cluster-score meat are unchanged by the within-transform (FWL preserves the score); the CR1 finite-sample scalar uses `K_reference` (3.9). The natural R anchor is `fixest::feols(y ~ | unit + time, cluster=~unit)` or Stata `jwdid` (both within-transform). **Parity with Stata `jwdid` (K_reference convergence, 3.9; `tests/test_etwfe_cs_stata_parity.py`):** the ATT(g,t) POINT estimates match `jwdid` exactly (~3e-8 on the `mpdta` panel, i.e. Stata's log-output rounding), and the `hc1` SEs now match `jwdid`/reghdfe at machine precision (ratio 1.0, spreads ~1e-15..1e-14) on all three committed arms AND at every rung of the committed G≈20..500 subsample ladder (`etwfe_cs_stata_golden.json` `ladder` block). Historically the SEs were uniformly SMALLER by a cluster-count-dependent factor (1.0280 at G=20 down to 1.0010 at G=500) — defect D2 (`docs/methodology/variance-conventions.md`): the clustered CR1 factor's `k` counted only the visible treatment-cell columns, omitting the absorbed FE not nested in the cluster. The applied correction is the closed form `K_reference = explicit columns + (1 if no intercept column) + rank(non-nested FE | nested)` (here `cells + T`), threaded as `cluster_k_adjustment`; the ladder doubles as the K-accounting probe (reghdfe `df_a == absorbed_fe_cr1_k_increment − 1` per rung). This is distinct from the `lm + clubSandwich` deviation below, whose factor is `k`-based. `CallawaySantAnna` shows no such gap - its SEs match Stata `csdid` outright - which localizes this to the ETWFE path rather than a library-wide convention. **Deviation from R `lm + clubSandwich::vcovCR(type="CR1S")`:** the full-dummy `lm` SE differs by a factor of `sqrt((n - K_reference) / (n - k_total))` because clubSandwich's `(n-1)/(n-p)` finite-sample correction counts ALL columns (intercept + treatment + unit dummies + time dummies = `k_total`) while WooldridgeDiD's within-transform CR1 counts `K_reference` (cells + constant + the non-nested time-FE rank; the cluster-NESTED unit dummies are the remaining difference). On the 240-obs / 51-column R-parity fixture (`k_total = 51`, `K_reference = 12`) this is `sqrt(228/189)` ≈ 10% (was ~11% under the historical cells-only count); on typical larger panels (n >> k_total) the gap shrinks to <2%. No public WooldridgeDiD code path exposes the `lm + CR1S` (CR1 cluster-robust on the full-dummy design) finite-sample correction — `vcov_type="hc2_bm"` routes to the CR2 Bell-McCaffrey sandwich on the full-dummy design (different variance estimator entirely), not CR1S. Users who need exact `lm + clubSandwich::vcovCR(type="CR1S")` parity must call `solve_ols` directly on a full-dummy design or fit via R. Same deviation pattern as SunAbraham PR #472 (`fixest::sunab` vs `lm + clubSandwich`). - `hc2_bm` — CR2 Bell-McCaffrey via auto-route to full-dummy design (`[intercept, X_design, unit_dummies, time_dummies]`), then `solve_ols(..., vcov_type="hc2_bm")` through the clubSandwich port (PR #475). FWL does NOT preserve the hat matrix; HC2 leverage + BM DOF require the full-projection design. Per-coefficient SE matches `clubSandwich::vcovCR(lm(...), cluster=~unit, type="CR2")` at atol=1e-10. Per-cell `(g, t)` inference fields use `coef_test()$df_Satt` Bell-McCaffrey DOF (pinned at atol=1e-6 from CI half-width inversion). Aggregated inference (overall ATT + `.aggregate("group" | "calendar" | "event")`) uses contrast-specific BM DOFs from `_compute_cr2_bm_contrast_dof` (matches R `Wald_test(constraints=matrix(w, 1), vcov=vcov_CR2, test="HTZ")$df_denom`); the overall ATT contrast DOF is computed at fit time, the other three aggregations lazily on each `.aggregate(...)` call from BM artifacts (the REDUCED kept-column `X` / `cluster_ids` / bread matrix + the reduced-space coef-index map) stored on the Results object — using the reduced design after rank-deficient drops keeps the bread non-singular and matches the subspace `solve_ols` actually estimated in. Fail-closed across all surfaces: when BM DOF is unavailable (helper raises or returns non-finite), the affected inference fields are NaN — not normal-theory fallback (per `feedback_bm_contrast_dof_fail_closed`). - `classical`, `hc2` — supported via auto-route to full-dummy AND auto-drop of the unit auto-cluster (one-way families don't compose with `cluster_ids` per the linalg validator). Set `self.cluster=None` (default) for these; explicit `cluster="state"` + one-way family raises at the linalg validator. SE matches `summary(lm(...))$coefficients` (classical) and `sandwich::vcovHC(type="HC2")` respectively. Per-cell + aggregate p-values/CIs use the residual DOF `n - rank(X)` (matches R `lm()` / `coef_test()` t-distribution under both classical OLS SE and `sandwich::vcovHC` defaults) — not normal-theory, so inference is correct under small samples. - `conley` (spatial-HAC, Conley 1999) — supported on the **OLS path** via the within-transform design (or the full-dummy design when `cohort_trends=True`, like the other full-dummy families — see the cohort-trends row below), threading the `conley_*` params through `solve_ols` / `conley.py` (`conley_lag_cutoff=0` = within-period spatial only; `>0` adds within-unit Bartlett serial — the panel-aware path, since `conley_time`/`conley_unit` are always supplied, not pooled cross-sectional). Reuses the already-`conleyreg`-validated machinery (no new variance code). The unit auto-cluster is dropped on the conley path (an explicit `cluster=` enables the spatial+cluster product kernel); `survey_design=` / `weights` / `n_bootstrap>0` are rejected, and `method ∈ {logit, poisson}` + conley remains rejected (the `method != "ols"` guard — a QMLE-on-pseudo-residuals Conley sandwich is a separate derivation). FWL-composability (the within-transform conley SE equals the full-dummy conley SE) is pinned in `tests/test_conley_vcov.py::TestConleyWooldridge::test_fwl_composability_vs_full_dummy`. @@ -2144,7 +2187,7 @@ where `g(·)` is the link inverse (logistic or exp), `η_i` is the individual li - **Identification (paper Section 8 / Eq. 8.1):** each treated cohort must have at least 2 pre-treatment periods (`t < g - anticipation`) for `dg_i · t` to be separately identified from cohort + time FE. `fit()` raises `ValueError` when the contract is violated. - **OLS-path only:** `cohort_trends=True` is rejected at `__init__` for `method ∈ {"logit", "poisson"}` per paper Section 8's OLS scope. `NotImplementedError` cites the paper section explicitly. - **Auto-routes to full-dummy mode** regardless of `vcov_type` (matching the absorb→fixed_effects auto-route pattern). Composing `dg_i · t` with the within-transformation yields `(dg_i − mean(dg_i)) · (t − mean(t))`, which is algebraically correct but non-trivial to verify on every panel shape; routing to the existing full-dummy auto-route used by `vcov_type ∈ {classical, hc2, hc2_bm}` keeps math closure verified against PR #483's R-parity goldens. UX implication: `cohort_trends=True` is silently more expensive than `cohort_trends=False` (carries N unit dummies); for very high-cardinality panels, the design-size warning at `wooldridge.py` fires. -- **`vcov_type="hc1"` finite-sample correction under `cohort_trends=True`:** the full-dummy auto-route changes the HC1 finite-sample factor from `(n-1)/(n-k_within)` (within-transform default) to `(n-1)/(n-k_total)` (full-dummy: counts intercept + treatment + unit + time + cohort-trend columns). On typical panels where `n >> k_total` the gap is small (<2%); on small panels it can reach ~10%. This is a documented opt-in deviation specific to `cohort_trends=True` — users who need the within-transform HC1 finite-sample factor with cohort trends should use `vcov_type="hc1"` + `cohort_trends=False` and supply the cohort-trend interactions through a custom design (out-of-scope for the standard library surface). +- **`vcov_type="hc1"` finite-sample correction under `cohort_trends=True`:** the full-dummy auto-route changes the HC1 finite-sample factor from `(n-1)/(n-K_reference)` (within-transform default; K_reference counts the non-nested absorbed FE, see § Variance families) to `(n-1)/(n-k_total)` (full-dummy: counts intercept + treatment + unit + time + cohort-trend columns, nested included — the L1 clubSandwich-CR1S-style opt-in convention, which the K_reference adjustment deliberately does not touch). On typical panels where `n >> k_total` the gap is small (<2%); on small panels it can reach ~10%. This is a documented opt-in deviation specific to `cohort_trends=True` — users who need the within-transform HC1 finite-sample factor with cohort trends should use `vcov_type="hc1"` + `cohort_trends=False` and supply the cohort-trend interactions through a custom design (out-of-scope for the standard library surface). - **Result attribute:** `WooldridgeDiDResults.cohort_trend_coefs: Dict[g → δ_g]` populated under `cohort_trends=True`; empty dict otherwise. - **Note:** Polynomial-trend extensions (`"quadratic"`, `"cubic"` per paper p. 2572 footnote) are NOT yet exposed — `cohort_trends` is a binary `True/False` flag for linear `dg_i · t` only. - **Note:** `cohort_trends=True` + `survey_design` is **NOT yet supported** (raises `NotImplementedError` at `fit()`). The full-dummy auto-route composed with the survey TSL variance has not been validated against R-parity goldens. Tracked in DEFERRED.md follow-up. @@ -2156,7 +2199,7 @@ where `g(·)` is the link inverse (logistic or exp), `η_i` is the individual li Consolidated list of substantive deviations from the W2025 paper and from R `etwfe`. Each is documented in the relevant section above with a labeled `**Note:**` or `**Deviation from R:**` line. AI PR reviewer recognizes these as documented (P3 informational) per the project's documented-deviation convention. 1. **Cell-count default for aggregation** (vs paper Eq. 7.4 / 7.6 cohort-share). `aggregate(weights="cell")` (default) matches Stata `jwdid_estat`. The opt-in `weights="cohort_share"` exposes the paper-Eq. 7.4 / 7.6 forms. Cohort-share is supported only for `type="simple"` and `type="event"`. See § Aggregations Note. -2. **HC1 finite-sample correction `(n-1)/(n-k_within)`** (vs R `lm + clubSandwich::vcovCR(type="CR1S")` which uses `(n-1)/(n-k_total)`). On 240-obs / 51-col fixture ~11%; on typical panels <2%. See § Variance families Deviation from R. +2. **HC1 finite-sample correction `(n-1)/(n-K_reference)`** (vs R `lm + clubSandwich::vcovCR(type="CR1S")` which uses `(n-1)/(n-k_total)`, nested dummies included). K_reference counts the non-nested absorbed FE + constant (matches Stata jwdid/reghdfe at machine precision); the remaining CR1S gap is the cluster-nested unit-dummy count. See § Variance families. 3. **QMLE sandwich `(G/(G-1)) · ((n-1)/(n-k))`** (vs Stata `jwdid` `G/(G-1)` only). Conservative; for typical panels n >> k the difference is negligible. Tracked in DEFERRED.md (WooldridgeDiD follow-up cluster row). See § Method Note. 4. **Nonlinear methods via direct QMLE** (vs R `etwfe` fixest backend). Avoids statsmodels/fixest dependency. See § Method Deviation from R. 5. **Logit cohort+time additive dummies** (not unit FE) to avoid incidental-parameters bias in short panels. Matches Stata `jwdid method(logit)`. See § Edge cases Note. @@ -2233,13 +2276,14 @@ Eq. 12 reuses the absorbing clean control and only restricts the treated set (a The paper specifies no standard-error formula (Section 1 defers to "standard, well-understood techniques"); the reference Stata `lpdid` uses `vce(cluster unit)`. The entries below document diff-diff's inference and scope choices. 1. **Note:** Standard errors are **cluster-robust at the unit level by default** - `cluster=None` auto-clusters at the unit identifier and the results record `cluster_name`/`n_clusters` - with a `t(G-1)` reference distribution (G = realized clusters in each horizon's clean-control sample). Matches Stata `lpdid` `vce(cluster unit)`; the paper prescribes no SE. -2. **Note:** The regression-adjustment (RA) covariate path (`reweight=True` with covariates/absorb) reports an **influence-function cluster variance** `sum_c (sum_{i in c} psi_i)^2 / n^2`, in the same family as `ImputationDiD`'s Theorem-3 / BJS variance (see "IF-based variance estimators vs analytical-sandwich estimators" above). Its single Gram inversion is routed through `linalg._rank_guarded_inv` (finite SE on the identified subspace under near-collinearity; NaN at rank 0). Unlike the default/weighted `solve_ols` `hc1`-cluster path - which applies the `(G/(G-1))*((n-1)/(n-k))` finite-sample factor - the RA IF variance carries **no finite-sample factor**, while both paths share the `t(G-1)` reference. **PR-B2 validated this asymmetry as faithful to the authors' own tooling**, not a defect: the no-factor RA convention matches the canonical Stata `teffects ra ... atet vce(cluster)`, while the default path matches `feols`/`reghdfe`. **Originally inferred from the authors' `lpdid_regression_adjustment.do` degrees-of-freedom comments, this is now MEASURED against a runnable `teffects`:** `benchmarks/stata/generate_lpdid_ra_golden.do` runs `teffects ra (Dy x i.time) (tdiff), atet vce(cluster unit)` on an independently reconstructed clean sample (a Stata port of the `generate_lpdid_golden.R` `prep`/`clean_h` recipe), and the library RA IF SE matches it to **~1e-16 at all 7 event-study horizons** - the SE ratio is a flat 1.0 (no `G/(G-1)`, no `(n-1)/(n-k)`), directly confirming the no-finite-sample-factor convention (golden `benchmarks/data/lpdid_ra_stata_golden.json`; test `tests/test_lpdid_ra_stata_parity.py`, which also gates the Stata clean sample's realized size and cluster count `(e(N), e(N_clust))` against the library's `(n_obs, n_clusters)`). The RA *point* estimate is triple-anchored: Python == R full-interaction (`i.dtreat##(i.time c.x)`) == Stata `teffects`, to ~1e-10..1e-13 (`tests/test_methodology_lpdid.py::test_ra_covariate_point`). Because Stata is node-locked and cannot run in CI (goldens are committed, exactly like the R arm), the CI-runnable guards remain the regression pin (`test_ra_covariate_se_regression_pin` / `RA_SE_PIN`) and the ungated Monte-Carlo coverage study `benchmarks/python/coverage_lpdid_ra.py` (~0.95 empirical coverage of the true effect at cluster counts G in {30, 100, 300}). +2. **Note:** The regression-adjustment (RA) covariate path (`reweight=True` with covariates/absorb) reports an **influence-function cluster variance** `sum_c (sum_{i in c} psi_i)^2 / n^2`, in the same family as `ImputationDiD`'s Theorem-3 / BJS variance (see "IF-based variance estimators vs analytical-sandwich estimators" above). Its single Gram inversion is routed through `linalg._rank_guarded_inv` (finite SE on the identified subspace under near-collinearity; NaN at rank 0). Unlike the default/weighted `solve_ols` `hc1`-cluster path - which applies the `(G/(G-1))*((n-1)/(n-K_reference))` finite-sample factor (see the nested-dummy Note below) - the RA IF variance carries **no finite-sample factor**, while both paths share the `t(G-1)` reference. **PR-B2 validated this asymmetry as faithful to the authors' own tooling**, not a defect: the no-factor RA convention matches the canonical Stata `teffects ra ... atet vce(cluster)`, while the default path matches `feols`/`reghdfe`. **Originally inferred from the authors' `lpdid_regression_adjustment.do` degrees-of-freedom comments, this is now MEASURED against a runnable `teffects`:** `benchmarks/stata/generate_lpdid_ra_golden.do` runs `teffects ra (Dy x i.time) (tdiff), atet vce(cluster unit)` on an independently reconstructed clean sample (a Stata port of the `generate_lpdid_golden.R` `prep`/`clean_h` recipe), and the library RA IF SE matches it to **~1e-16 at all 7 event-study horizons** - the SE ratio is a flat 1.0 (no `G/(G-1)`, no `(n-1)/(n-k)`), directly confirming the no-finite-sample-factor convention (golden `benchmarks/data/lpdid_ra_stata_golden.json`; test `tests/test_lpdid_ra_stata_parity.py`, which also gates the Stata clean sample's realized size and cluster count `(e(N), e(N_clust))` against the library's `(n_obs, n_clusters)`). The RA *point* estimate is triple-anchored: Python == R full-interaction (`i.dtreat##(i.time c.x)`) == Stata `teffects`, to ~1e-10..1e-13 (`tests/test_methodology_lpdid.py::test_ra_covariate_point`). Because Stata is node-locked and cannot run in CI (goldens are committed, exactly like the R arm), the CI-runnable guards remain the regression pin (`test_ra_covariate_se_regression_pin` / `RA_SE_PIN`) and the ungated Monte-Carlo coverage study `benchmarks/python/coverage_lpdid_ra.py` (~0.95 empirical coverage of the true effect at cluster counts G in {30, 100, 300}). 3. **Note:** Direct covariate inclusion (`reweight=False` with covariates/absorb) emits a `UserWarning`: per online Appendix B.2.2 it preserves the non-negative LP-DiD weighting result only under linear and homogeneous covariate effects, so the regression-adjustment path (`reweight=True`) is preferred. 4. **Deviation from R:** Scope - non-absorbing treatment (Section 4.2) implements the **entry-effect** estimands (`non_absorbing="first_entry"` / `"effect_stabilization"`, PR-C1). **PR-C2 R-parity-validated both modes against an INDEPENDENT `fixest::feols` reconstruction of the paper's Eq. 12 / Eq. 13 clean-sample restrictions** (point and SE match to ~1e-13/~1e-15 for the variance-weighted variants; the `effect_stabilization` reweighted point matches and its SE is pinned as a regression guard - a small weighted-cluster convention difference vs feols; `tests/test_methodology_lpdid.py::TestLPDiDNonAbsorbingParityR`). The recipe's independence was demonstrated when an earlier draft's Eq. 12 control off-by-one diverged from the already-correct library and was corrected against the paper, plus a hand-computed Python micro-check. **`alexCardazzi/lpdid`'s `nonabsorbing_lag` is NOT a faithful Eq. 13** (it clamps `treat_diff[<0]<-0`, so its clean-control window blocks only treatment turn-*ons*; it reuses a forward placebo window; and it NA-excludes pre-panel-treated rows where the library clamps pre-`min_t` to untreated): it diverges ~0.01-0.05 from Eq. 13 even on a monotone no-off-switch panel, so it is **recorded in the golden `meta` as a divergent third-party reference, not a parity gate** (the alexCardazzi-pooled precedent). The library's "no treatment change" (both directions) and backward placebo window are the more paper-faithful choices. `first_entry` (Eq. 12) has no R-package analogue (anchored on the independent feols recipe only). Appendix-C exit-event dynamics and the Stata canonical SE remain deferred follow-ups. 5. **Note:** LP-DiD's per-unit quantities (outcome lags `ylags`, first-difference lags `dylags`, integer-`pmd` premean baselines, treatment-entry detection) are **calendar** quantities (`t-1`, `t-k`), so the estimator requires integer-valued, globally consecutive `time` labels. A unit with an **interior time gap** is handled by reindexing that unit to its complete interior calendar grid `[min_t, max_t]`, computing the features on the grid, then **restricting back to the observed rows** - so a lag/first-difference spanning a gap is NaN and the observation fails closed (never the previous-*observed* row), and no synthetic gap row enters a regression. A gap-free panel skips this entirely and is bit-identical. **Entry = first OBSERVED treated period** (`min(t | D_it=1)`): an unobserved pre-onset gap cannot move a cohort earlier, the only well-defined convention when the true switch falls in an unobserved period. 6. **Note (pooled estimand):** The pooled pre/post ATT (the headline `results.att` is the pooled-post row) is the **unit-equal-weighted average of each unit-event-time's mean long difference** over the window - `mean_h(y_{i,t+h}) - baseline_{i,t}`, one observation per (unit, event-time), regressed on the treatment-switch indicator with event-time fixed effects on the **fixed-composition** sample (only units observing *every* pooled target, with clean controls required through `max(h)`). This equals the mean of the per-horizon event-study coefficients on a balanced panel. **PR-B2 validated it against the authors' runnable R reference**: the pooled estimand matches the authors' own R pooled recipe (`danielegirardi/lpdid`: a `slider` window-mean minus `y_{t-1}` on the clean-through-window-end sample) to ~1e-13 (`tests/test_methodology_lpdid.py::test_pooled`). A prior version of this note speculated the authors used a horizon-**stacked** pooled regression; the authors' R reference in fact uses this same fixed-composition mean-long-difference, so that speculation was incorrect. Unlike the event-study variants (where `alexCardazzi` is a cross-check gate), pooled is anchored to the authors' recipe **only**: `alexCardazzi`'s pooled uses a **laxer** clean-control window, so it differs and is recorded in the golden `meta` for transparency, not as a parity target. 7. **Deviation from R:** `no_composition` is intentionally more faithful to the paper's fixed-composition intent (Section 3.6) than the R packages: it fixes the realized sample across *all* post horizons (every post coefficient shares one sample, even on unbalanced panels) and excludes cohorts with `p_g > T-H`, whereas `alexCardazzi/lpdid` uses a looser per-horizon sample and a stricter `treat_date < T-H` cutoff. It therefore has **no exact R-package anchor** and is validated by the pure-Python tests in `tests/test_lpdid.py` (the R-parity golden omits it; `alexCardazzi`'s looser-semantics value is recorded in the golden `meta`). 8. **Note (survey design):** Complex-survey support (`survey_design=SurveyDesign(...)`, PR-D1) covers the **variance-weighted default path** (`reweight=False`, with or without direct-inclusion covariates): each horizon's long-difference regression is fit by WLS on the survey probability weights, and the SE is the stratified-PSU **Taylor-linearization (Binder 1983 TSL)** sandwich `meat = sum_h (1-f_h)*(n_h/(n_h-1))*sum_j (S_hj - S_h_bar)(S_hj - S_h_bar)'` with `df = n_PSU - n_strata`, reusing the shared `diff_diff/survey.py` helpers (`compute_survey_vcov` / `_compute_stratified_psu_meat`). The design is re-resolved on each realized (post-clean-control) sample so weights/strata/PSU align with the regression rows; with no explicit PSU the unit (LP-DiD's default cluster) is injected as the PSU. Supports pweight + strata + PSU + FPC + lonely-PSU handling. It **rejects** `survey_design` combined with `reweight=True` (the equally-weighted / regression-adjustment IF path has no validated survey reference - the same gap as the RA SE in Deviation #2), replicate-weight designs, and non-pweight (fweight/aweight) types, each a deferred follow-up. The non-survey path is byte-for-byte unchanged (gated on `survey_design is None`). **PR-D2 validated all three survey paths end-to-end against `survey::svyglm`** - per-horizon point/SE/df + pooled for the variance-weighted full design (strata+PSU+FPC), the weights-only unit-injected-PSU design, and the direct-covariate variant (`tests/test_methodology_lpdid.py::TestLPDiDSurveyParityR`; point ~1e-6, SE ~1e-5, df exact via the per-design `n_PSU - n_strata` / `n_PSU - 1` formula). `svyglm` is itself the reference implementation of the Binder TSL sandwich, so it anchors the variance directly (no third-party survey-package gate is needed); the clean-sample construction is independently cross-checked in the generator (the unweighted variance-weighted event study matches `alexCardazzi/lpdid` to <1e-8, and selection is weight-independent). A dedicated survey panel (`benchmarks/data/lpdid_survey_panel.csv`, own seed) keeps the absorbing / non-absorbing goldens byte-identical. +9. **Note (clustered CR1 k, nested dummy blocks — 3.9 K_reference):** LPDiD's per-horizon design carries two built-in dummy blocks — the `_event_time` (time-FE) dummies added under `include_time_fe=True` and the inline `drop_first` `absorb=` dummies. The authors' reference recipe ABSORBS the time FE (`feols(... | time, vcov = ~unit)`, `benchmarks/R/generate_lpdid_golden.R`), so both blocks follow the nested convention: when a block is nested in the raw `cluster` (e.g. a time-like cluster for `_event_time`, or an absorb dim coarser than the cluster), its rank is SUBTRACTED from the CR1 factor's `k` via a negative `cluster_k_adjustment`, converging on the same `K_reference` accounting as `feols`/`reghdfe`. Non-nested configs are no-ops (unit clustering leaves `_event_time` non-nested → adjustment 0, preserving bit-identical Rust-lane output). The solve's broad `except (ValueError, ZeroDivisionError)` re-solve re-raises `InvalidClusterKAdjustment` so a contract violation can never degrade to a silent `se=NaN` (`tests/test_variance_conventions.py::TestKReferenceConvergence::test_lpdid_never_swallows_the_contract_raise`). The `t(G-1)` tail reference (Note 1) is unchanged. ### Implementation Checklist diff --git a/docs/methodology/variance-conventions.md b/docs/methodology/variance-conventions.md index 06676b498..906fd4dbc 100644 --- a/docs/methodology/variance-conventions.md +++ b/docs/methodology/variance-conventions.md @@ -22,50 +22,78 @@ produced wrong figures three separate times while this inventory was drafted). | surface | CR1 `k` (multiset) | tail df (multiset) | status | reason | |---|---|---|---|---| -| `did_absorb_hc1_cluster_unit` | 2 | 294 | **defect** | D2: CR1 k omits absorbed FE not nested in the cluster (time) | -| `did_fixed_effects_hc1_cluster_unit` | 66 | 294 | **defect** | D1: same model as did_absorb yet k=66 vs 2 -> SEs differ 10.35%; full-dummy k also counts the cluster-nested unit FE the references drop | +| `did_absorb_hc1_cluster_unit` | 7 | 294 | **legitimate** | K_reference (D2 fixed): k = 2 visible + rank(time given unit) = 5; matches reghdfe/fixest ssc(K.fixef='nested') | +| `did_fixed_effects_hc1_cluster_unit` | 7 | 294 | **legitimate** | K_reference (D1 fixed): 66 visible minus the 59 cluster-nested unit dummies -> identical SE to did_absorb (documented deviation from a literal explicit-dummy R comparison, which counts all 66) | | `did_plain_hc1_cluster_unit` | 4 | 356 | **legitimate** | no absorbed FE: visible k is the whole design; nothing is omitted | -| `twfe_hc1_cluster_unit_time_post` | 2 | 298 | **defect** | D2 (within-transform k_visible); tail df is residual n-K_full | -| `wooldridge_hc1_within` | 9 | None, None, None, None, None, None, None, None, None, None | **defect** | D2 (k_visible=cells only) + normal-theory tail df with no df_convention knob | -| `sun_abraham_hc1` | 15 | 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, None, None, None, None, None, None, None, None | **defect** | D2 + D4: residual df per cohort-period cell but normal theory on aggregates | +| `twfe_hc1_cluster_unit_time_post` | 3 | 298 | **legitimate** | K_reference (D2 fixed): 2 visible + rank(post given unit) = 1; matches fixest cluster arm at rel 0 (committed golden) | +| `wooldridge_hc1_within` | 15 | None, None, None, None, None, None, None, None, None, None | **defect** | CR1 k converged on K_reference (D2 fixed: 9 cells + T = 15, no intercept col -> +1 term; jwdid arms at ratio 1.0); tail df is still normal theory with no df_convention knob (PR C) | +| `sun_abraham_hc1` | 21 | 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, None, None, None, None, None, None, None, None | **defect** | CR1 k converged on K_reference (D2 fixed: 15 cells + 6, no intercept col; fixest sunab parity ~5e-15); D4 remains: residual df per cohort-period cell but normal theory on aggregates (PR C) | | `stacked_did_hc1` | 6 | None | **legitimate** | L1: k_total is clubSandwich CR1S by construction (stacked_did.py pins vcovCR(type='CR1S') at atol=1e-10); normal-theory tail df is an open PR C question | | `lpdid_pre2_post2` | 4, 4, 5, 5, 5, 6 | 59, 59, 59, 59, 59, 59 | **legitimate** | L2: G-1 tail df (Stata/fixest convention) — the convergence target | +| `mpd_absorb_hc1_cluster_unit` | 11 | 290, 290, 290, 290, 290, 290 | **legitimate** | K_reference: 6 visible + rank(time given unit) = 5; equals the fixed_effects form's 70 - 59 (MPD absorb/fixed_effects equivalence) | +| `mpd_fixed_effects_hc1_cluster_unit` | 11 | 290, 290, 290, 290, 290, 290 | **legitimate** | K_reference: 70 visible (incl. built-in period dummies, MPD's time-FE block) minus the 59 cluster-nested unit dummies = 11 — identical to the absorb form | +| `mpd_plain_hc1_cluster_time` | 7 | 348, 348, 348, 348, 348, 348 | **legitimate** | the NESTED orientation of the built-in period dummies: 12 visible minus their rank 5 under a time cluster (under-subtraction is caught here; the unit-cluster rows catch over-) | +| `lpdid_absorb_nested_cluster_grp` | 4, 4, 5, 5, 5, 6 | 1, 1, 1, 1, 1, 1 | **legitimate** | LPDiD absorb dummies nested in the cluster subtract their rank (adj -1 per horizon: region == grp here); _event_time stays counted (unit-level cluster does not nest time); G-1 tail df (L2) | | `imputation_default` | — | None | **legitimate** | L3: BJS imputation variance, not the shared CR1 sandwich | -| `imputation_pretrends_event_study` | unpinned | unpinned | **defect** | pretrends lead regression runs the shared clustered CR1 with k_visible | +| `imputation_pretrends_event_study` | unpinned | unpinned | **defect** | pretrends lead regression: CR1 k converged on K_reference (D2 fixed); normal-theory tail df remains (PR C family) | | `two_stage_default` | — | None | **legitimate** | L3: Gardner two-stage variance, not the shared CR1 sandwich | | `callaway_santanna_default` | — | None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None | **legitimate** | L3: influence-function variance anchored to Stata csdid | -cr1_k is the sorted multiset of visible column counts reaching the shared -clustered CR1 denominator (linalg._compute_robust_vcov_numpy with +cr1_k is the sorted multiset of K_reference counts reaching the shared +clustered CR1 denominator — visible columns + the signed cluster_k_adjustment +(linalg._compute_robust_vcov_numpy with vcov_type="hc1"; a clustered call in any other family fails the row, so a surface cannot silently switch clustered family behind an unchanged design width); — means the surface's *contract* is that it never calls it. tail_df is the multiset of df values passed to safe_inference/safe_inference_batch (None = normal theory). unpinned marks the one contract row whose test asserts the shared CR1 IS reached but deliberately pins no exact values — -a literal there would be brittle configuration-detail (last measured: -k=3, normal-theory tail df on all 8 event-study leads/lags). Captured under the canonical Python backend; Rust -and Python agree to <= 8e-15 on every surface because both implement the same -conventions today. +a literal there would be brittle configuration-detail (its expected_adjustment +IS pinned: +6, the [time, unit] increment on df_0). Captured under the +canonical Python backend; on adjusted clustered surfaces the Rust lanes apply +the same correction as an exact scalar rescale of the finished vcov +(vcov * (n-k)/(n-k-adj), <= 1 ulp on top of the default lane's <= 8e-15 +cross-backend agreement; zero-adjustment surfaces keep bit-identical Rust +output). -## Defects (scheduled: 3.9 consolidation program) +## Fixed defects (3.9 consolidation program) -- **D1 — absorb= vs fixed_effects=: 10.35% SE split on the same model.** - Identical ATT, se 0.2414226781 (k=2) vs 0.2664071714 (k=66), ratio - 1.103489 = sqrt((360-2)/(360-66)) exactly. The absorb-side deviation from - fixest is documented (REGISTRY "Known limitation (deviation from fixest)", - DEFERRED.md); the *user-facing consequence* — two documented-equivalent kwargs - disagreeing by 10%, in opposite directions from the reference — was not, until - this inventory. Fix: PR B converges both on K_reference. -- **D2 — clustered CR1 never counts absorbed FE not nested in the cluster.** - k_visible in the denominator; _absorbed_fe_vcov_scale exists but is - gated on cluster_ids is None. Anti-conservative: SEs understated by 0.10% - (n=2500), 1.30% (n=200), 5.51% (n=60). The correction +- **D1 — absorb= vs fixed_effects=: 10.35% SE split on the same model. + FIXED (3.9, K_reference convergence).** Historically: identical ATT, se + 0.2414226781 (k=2) vs 0.2664071714 (k=66), ratio 1.103489 = + sqrt((360-2)/(360-66)) exactly, in opposite directions from the reference. + Both idioms now produce the identical K_reference = 7 on the audit panel + (absorb: 2 + 5; full-dummy: 66 - 59 nested rank) and the identical SE + 0.2414226781 * sqrt((360-2)/(360-7)), pinned by + test_d1_convergence_is_pinned. The equivalence is rank-based, so it holds + under collinearity drops and on disconnected panels (kernel k is the design + RANK; the two-nested-dim discriminator in TestKReferenceConvergence). + The full-dummy side deviates from a LITERAL explicit-dummy R comparison + (which counts all 66) — the labeled deviation note in REGISTRY + "Deviation from R (clustered CR1, full-dummy fixed_effects= path)". +- **D2 — clustered CR1 never counted absorbed FE not nested in the cluster. + FIXED (3.9, cluster_k_adjustment seam).** Historically anti-conservative: + SEs understated by 0.10% (n=2500), 1.30% (n=200), 5.51% (n=60). The + correction K_reference = explicit cols + (1 if no intercept col) + rank(non-nested FE | nested) - reproduces Stata reghdfe 3.2.9 (via jwdid) to ~1e-15 on three arms and - R fixest 0.14.2 to ~1e-12 on two, and retrodicts the two measured-but-unpinned - subsample rungs (predicted 1.028016 vs recorded 1.0280 at G=20; 1.013210 vs - 1.0132 at G=40). Fix: PR B. + (the +1 is the rank of the ABSORBED constant, so it exists only when FE are + absorbed and X carries no intercept column) is computed by + utils.absorbed_fe_cr1_k_increment and threaded keyword-only through + solve_ols / compute_robust_vcov / LinearRegression.fit / wild_bootstrap_se. + External anchors: Stata reghdfe 3.2.9 (via jwdid) at machine precision on + all three committed mpdta arms AND at every rung of the committed G≈20..500 + subsample ladder (ratio 1.0, spreads ~1e-15..1e-14 — the ladder also probes + the accounting itself: reghdfe df_a == increment - 1 per rung); R fixest + 0.14.2 cluster arms tightened to exact/1e-9. The non-nested term is the + exact RANK given the nested set: on disconnected designs this deviates by + one df per extra component from BOTH references (neither implements the + exact composition — fixest's K.exact composes incoherently with its nested + drop, reghdfe's pairwise correction skips nested-dropped dims), pinned as an + exact sqrt((n-10)/(n-11)) SE ratio in the committed + fixest_cr1_nonnested_golden.json / reghdfe_kref_golden.json arms, while the + nothing-nested crossed-cluster arms match both references at machine + precision. Fail-closed: n_eff - k <= 0 (visible saturation), n_eff - k_inf + <= 0, or k_inf <= 0 each yield the all-NaN vcov on every backend. - **D3 — absorbed rank assumed independent, connected FE dimensions.** sum(levels - 1) over-counted on disconnected panels (true rank sum(levels) - C) and on hierarchical specs @@ -96,7 +124,8 @@ conventions today. shared CR1 sandwich. CS is anchored to Stata csdid outright. **ImputationDiD is conditional**: its pretrends=True + aggregate="event_study" lead regression DOES run the shared clustered CR1 - and inherits D2 there (its own matrix row). + and carries the K_reference increment there (+6, the [time, unit] no-intercept + increment on df_0 — pinned via expected_adjustment on its matrix row). - **L4 — hc2/hc2_bm** (leverage / Satterthwaite DOF — no CR1 factor), **survey TSL** (n_PSU - n_strata over the full design), and **Wooldridge cohort_trends full-dummy** (documented opt-in landing on the L1 @@ -110,7 +139,8 @@ Three conventions are live: normal theory (Wooldridge, StackedDiD — no df_convention knob), residual n - K_full (DiD/MPD/TWFE default; df_convention="cluster" opts into G-1), and G-1 (LPDiD, hardcoded). At |t| = 2 normal theory understates the t(G-1) p-value by 24.2% at -G=20, 13.3% at G=40, 1.2% at G=500 — larger than the D2 SE gap. PR C decisions: +G=20, 13.3% at G=40, 1.2% at G=500 — larger than the D2 SE gap was before its +fix. PR C decisions: (1) extending the knob to Wooldridge/SunAbraham/StackedDiD/ImputationDiD needs NEW ledger rows (M-004/M-005/M-006 cover only DiD/TWFE/LinearRegression); (2) the two-value knob cannot express normal theory, so either a third value diff --git a/docs/performance-plan.md b/docs/performance-plan.md index 86f0732ae..dba55237d 100644 --- a/docs/performance-plan.md +++ b/docs/performance-plan.md @@ -10,15 +10,21 @@ This document outlines the strategy for improving diff-diff's performance on lar `pd.factorize` calls + a sparse weak-connectivity `connected_components` on the bipartite level graph; the directed weak-connectivity form skips materializing `A + A.T`, ~2.4x over the naive build). Measured on the county-class shape -(3,100 units x 60 periods, 186k rows): **1.9 ms/call**, ~10 ns/row. Against the -Rust-served TWFE hc1 fit on that shape (24 ms) the share is **7.7%** — above the -2% materiality gate set in the consolidation plan, but the denominator is the -library's fastest fit configuration; on fits with covariates, bootstrap, or the -pure-Python backend the share falls well under 2%, and the absolute cost is flat -O(n) with no allocation cliffs. The remaining mitigation — reusing the factorized -codes `demean_by_groups`/`within_transform` already produce instead of -re-factorizing — requires threading demeaner internals through four call sites -and is deliberately deferred (tracked alongside the D3 N-way row in TODO.md); +(3,100 units x 60 periods, 186k rows): **1.9 ms/call**, ~10 ns/row. The +K_reference convergence (3.9 PR B) added +`diff_diff.utils.absorbed_fe_cr1_k_increment` on clustered-hc1 absorbed fits: +**3.2 ms/call** on the same shape — the shared factorizations plus the +`np.minimum.at`/`np.maximum.at` nested-dims scatter, a cluster factorize, and +the second (nested-set) rank composition. Against the Rust-served TWFE hc1 fit +on that shape (24 ms) the shares are **7.7%** (rank alone, unclustered fits) +and **~13%** (the increment, clustered fits) — above the 2% materiality gate +set in the consolidation plan, but the denominator is the library's fastest fit +configuration; on fits with covariates, bootstrap, or the pure-Python backend +the share falls well under 2%, and the absolute cost is flat O(n) with no +allocation cliffs. The remaining mitigation — reusing the factorized codes +`demean_by_groups`/`within_transform` already produce instead of re-factorizing +— requires threading demeaner internals through the call sites and is +deliberately deferred (tracked alongside the D3 N-way row in TODO.md); correctness on disconnected/hierarchical panels ships first. ## Opt-in solve_ols normal-equations Cholesky fast path (v3.7, 2026-07) diff --git a/docs/v4-deprecations.yaml b/docs/v4-deprecations.yaml index 2d013e4d6..a1b025b28 100644 --- a/docs/v4-deprecations.yaml +++ b/docs/v4-deprecations.yaml @@ -1327,6 +1327,19 @@ rows: test_ref: tests/test_aggregate_contract.py code_refs: [diff_diff/staggered.py, diff_diff/staggered_results.py, diff_diff/aggregation.py] notes: "balance_e moves from fit() onto aggregate() with [M-020]. Previously tracked only as the prose 'balance_e moves to aggregate() in the same PR' inside M-020's notes, which no test could assert - the same un-rowed-obligation class the gating-completeness amendment closed. Applies to event-study aggregation ONLY (the shipped code threads it nowhere else), so aggregate(type='simple'|'group', balance_e=...) raises rather than silently ignoring it. The other three balance_e sites (ImputationDiD, TwoStageDiD, EfficientDiD) get their own rows in the PRs that migrate them." + - id: M-126 + kind: behavior + group: variance-cr1-convergence + old: "diff_diff.linalg:solve_ols" + new: null + introduced_in: "3.9" + deprecated_in: null + removed_in: null + status: done + phase: 2 + test_ref: tests/test_variance_conventions.py + code_refs: [diff_diff/linalg.py, diff_diff/utils.py, diff_diff/estimators.py, diff_diff/twfe.py, diff_diff/sun_abraham.py, diff_diff/wooldridge.py, diff_diff/imputation.py, diff_diff/lpdid.py, docs/methodology/REGISTRY.md, docs/methodology/variance-conventions.md] + notes: "Clustered CR1 K_reference convergence (D1+D2 of the 3.9 variance-consolidation program, PR B). The clustered hc1 CR1 finite-sample factor (G/(G-1))*((n_eff-1)/(n_eff-k)) now uses K_reference = explicit columns + (1 for the absorbed constant when X carries no intercept column) + rank(absorbed FE not nested in the cluster | nested ones), threaded as the NEW keyword-only signed `cluster_k_adjustment: int = 0` on solve_ols / compute_robust_vcov / LinearRegression.fit / wild_bootstrap_se (additive kwargs, default 0 - no deprecation cycle) and computed by the new utils.absorbed_fe_cr1_k_increment / cluster_nested_fe_dims helpers. Fixes D2 (absorbed FE not nested in the cluster were never counted: SEs understated 0.10% n=2500 / 1.30% n=200 / 5.51% n=60) and D1 (absorb= vs fixed_effects= returned clustered SEs differing by exactly sqrt((360-2)/(360-66)) = 10.35% on the same model; the full-dummy lane now subtracts the nested dummy RANK, so both idioms return the identical externally-anchored SE). Wired surfaces: DiD/MPD/TWFE (absorb increment; fixed_effects negative nested rank; MPD's built-in period dummies are its time-FE block and follow the nested convention), SunAbraham (+rank, closes the fixest::sunab ~1-2% gap, measured 3.1e-14 rel), WooldridgeDiD within (the original jwdid SE-gap finding: machine-precision parity on all committed arms + the new G~20..500 subsample ladder), ImputationDiD pretrends lead regression, LPDiD nested dummy blocks (with an InvalidClusterKAdjustment re-raise so the broad except cannot swallow the contract), and the wild-cluster-bootstrap corr constant (cancels in |t*| vs |t0|: p-values invariant, reported SE aligns). StackedDiD (clubSandwich CR1S by construction) and Wooldridge cohort_trends full-dummy keep their documented L1 convention (adjustment 0, bit-identical). Rust lanes apply the correction as an exact None-preserving scalar rescale of the finished vcov in the Python wrappers (zero Rust changes; <= 1 ulp; zero-adjustment surfaces bit-identical). Fail-closed on all three sides: n_eff-k <= 0 (visible saturation, residuals identically zero), n_eff-k_inf <= 0, k_inf <= 0 -> all-NaN vcov on every backend; front-door validation raises InvalidClusterKAdjustment(ValueError) on a nonzero adjustment without cluster_ids or off the hc1 family (survey fits are the documented inert exception). External anchors: Stata reghdfe 3.2.9 via jwdid at machine precision on all arms and ladder rungs; fixest cluster arms tightened to exact/1e-9; the exact non-nested RANK term anchored on a disconnected panel vs reghdfe AND fixest ssc(K.fixef='full', K.exact=TRUE) (crossed-cluster arms, ~1e-16), with the nested-composition one-df deviation from both references pinned exactly (no external reference implements nested-drop + exact-remainder; fixest_cr1_nonnested_golden.json / reghdfe_kref_golden.json). Tail df (df_, dof_vec, inference_df, df_convention) deliberately unchanged - PR C scope (M-004..M-006 family). introduced_in gates the 3.9 cut; deprecated_in stays null so the early-flip guard does not fire against the PR that ships it." - id: M-125 kind: behavior group: etwfe-reference-period diff --git a/docs/v4-design.md b/docs/v4-design.md index 327fb46c4..1c2c39151 100644 --- a/docs/v4-design.md +++ b/docs/v4-design.md @@ -743,11 +743,11 @@ forever - a removed symbol resurrecting is a test failure. class row (schema-enforced). Top-level `diff_diff:Name` class/function rows and alias rows also assert `__all__` membership consistent with their status (stale `import *` entries fail). The shipped row ids are a - committed snapshot in the enforcement test (104 as of the ETWFE - comparison-support row: Phase 1 + the diagnostic-family amendment + + committed snapshot in the enforcement test (105 as of the clustered-CR1 + K_reference convergence row: Phase 1 + the diagnostic-family amendment + the M-092/M-093 results-contract rows + the M-094..M-096 amendment rows + the M-097..M-115 completeness sweep + M-117/M-122 + the ETWFE - reference-period pair M-123/M-124 + M-125; + reference-period pair M-123/M-124 + M-125 + M-126; the snapshot extends by a new id range in the same diff that appends rows): ids are never deleted or reused, and the test fails if any snapshot id disappears. diff --git a/tests/test_estimators_vcov_type.py b/tests/test_estimators_vcov_type.py index d869cce15..d244fee6d 100644 --- a/tests/test_estimators_vcov_type.py +++ b/tests/test_estimators_vcov_type.py @@ -2629,16 +2629,15 @@ def test_twfe_classical_matches_full_dummy_oracle(self): np.testing.assert_allclose(tw.att, fe.att, rtol=1e-9) np.testing.assert_allclose(tw.se, fe.se, rtol=1e-9) - def test_absorb_cluster_not_rescaled(self): - """The absorbed-FE full-K rescale must NOT touch clustered SEs. - - The rescale is gated on ``cluster_ids is None``, so the cluster-absorb - SE stays at ``k_visible`` and differs from the full-dummy path here. - (Full fixest cluster parity is a *separate*, out-of-scope matter: fixest - counts non-nested absorbed FE in the CR1 denominator, so for - ``absorb=["unit","time"], cluster="unit"`` the non-nested time FE would - need counting -- a documented pre-existing limitation, see REGISTRY; - this test only pins that D4 does not rescale the cluster path.) + def test_absorb_cluster_converges_on_k_reference(self): + """absorb= and fixed_effects= clustered SEs are now IDENTICAL. + + Both lanes land on the same K_reference count (visible columns + the + conditional rank of absorbed FE not nested in the cluster; the + full-dummy lane subtracts its cluster-nested dummies) — the D1 + convergence from the 3.9 variance-consolidation program, anchored by + the Stata reghdfe (jwdid) and fixest cluster arms. Before the fix + this test asserted the two lanes DIFFERED (the pre-fix defect). """ df = _make_absorb_panel() ab = DifferenceInDifferences(vcov_type="hc1", cluster="unit").fit( @@ -2658,10 +2657,8 @@ def test_absorb_cluster_not_rescaled(self): covariates=["x"], ) assert np.isfinite(ab.se) - assert not np.isclose(ab.se, fe.se, rtol=1e-6), ( - "cluster-absorb SE must keep k_visible (nested-FE), not be rescaled to " - "the full-dummy K_full value" - ) + np.testing.assert_allclose(ab.att, fe.att, rtol=0, atol=1e-10) + np.testing.assert_allclose(ab.se, fe.se, rtol=1e-9) def test_absorb_hc2_bm_not_rescaled(self): """hc2_bm auto-routes absorb -> full-dummy and uses Satterthwaite DOF; diff --git a/tests/test_etwfe_cs_stata_parity.py b/tests/test_etwfe_cs_stata_parity.py index f4b61b848..3a019594c 100644 --- a/tests/test_etwfe_cs_stata_parity.py +++ b/tests/test_etwfe_cs_stata_parity.py @@ -232,29 +232,23 @@ def test_att_matches_jwdid(self, never_fit): err_msg=f"never_treated ATT{key} != Stata jwdid never", ) - def test_se_gap_matches_the_not_yet_treated_arm(self, never_fit): - """The SE gap is a control-group-INDEPENDENT finite-sample factor. - - The golden stores ``jwdid_never`` SEs, so the path this PR actually - changed gets its own SE gate rather than inheriting the not-yet-treated - one. Measured here across all 12 never-treated cells (7 post + 5 - placebo): the ratio is uniform to ~1e-14 -- two orders tighter than the - 1e-5 the sibling test allows -- and its mean, 1.00100634, is the SAME - factor the not-yet-treated arm shows. - - That agreement is the point. Whatever produces the gap does NOT depend - on the control group, on the cell being post-treatment vs placebo, or - on the cell count (12 here vs 7 there), which constrains any mechanism - later proposed for it (tracked in TODO.md). As with the sibling test, - no closed form is asserted -- only the observed factor. + def test_se_matches_jwdid_never_treated_arm(self, never_fit): + """SEs now reproduce Stata jwdid (reghdfe) at machine precision. + + The historical 1.00100634 SE ratio on this arm was the D2 defect: + the clustered CR1 factor's ``k`` omitted the absorbed FE not nested + in the unit cluster. Under the K_reference convergence + (variance-conventions.md, 3.9 program) all 12 never-treated cells + sit at ratio 1.0 with spread ~1e-14 — the same machine-precision + agreement the not-yet-treated and all-eventually-treated arms show. + Two-sided by design: a strict inequality on a quantity at + 1.0 ± float noise would be a coin flip. """ stata = _stata_jwdid_cells(_golden(), "jwdid_never") ratios = [rec["se"] / never_fit.group_time_effects[key]["se"] for key, rec in stata.items()] assert len(ratios) == 12, f"expected 12 never_treated cells, got {len(ratios)}" - assert max(ratios) - min(ratios) < 1e-5, f"SE ratio is not uniform: {ratios}" - mean_ratio = float(np.mean(ratios)) - assert mean_ratio > 1.0, "library SE should be below jwdid's, not above" - np.testing.assert_allclose(mean_ratio, 1.001006, rtol=1e-4) + assert max(ratios) - min(ratios) < 1e-9, f"SE ratio is not uniform: {ratios}" + np.testing.assert_allclose(float(np.mean(ratios)), 1.0, rtol=1e-9) class TestETWFEvsStataJwdid: @@ -275,45 +269,27 @@ def test_att_matches_jwdid(self, fits): compared += 1 assert compared == 7, f"expected 7 ETWFE cells, compared {compared}" - def test_se_is_uniformly_below_jwdid(self, fits): - """MEASURED deviation, recorded without asserting a mechanism. - - Every ``hc1`` SE is smaller than ``jwdid``'s by a factor that is - UNIFORM across cells (spread < 1e-6 within a fit), so it is a - finite-sample convention difference rather than noise or a per-cell - bug. The magnitude shrinks as the cluster count grows -- 1.0280 at - G=20, 1.0132 at G=40, 1.0010 at G=500 -- so the library is - systematically anti-conservative relative to the reference, negligibly - with many clusters and materially with few. - - Deliberately NOT asserted: a closed form. The gap tracks - ``sqrt(G/(G-1))`` closely but is consistently ABOVE it (by ~0.2% at - G=20, ~0.001% at G=500), and ``solve_ols`` already applies - ``(G/(G-1)) * ((n-1)/(n-k))`` -- so "the library omits G/(G-1)" is - NOT the explanation, and no verified formula has been derived. Pinning - a wrong mechanism here would be worse than pinning none: this test - locks the OBSERVED gap at this panel's G, and the derivation is tracked - in TODO.md. See the REGISTRY WooldridgeDiD note. - - SCOPE: this pins the FULL-PANEL ratio (G=500) only. The smaller-G - figures quoted in the REGISTRY note were measured ad hoc on subsampled - panels and are not gated here, so a regression in the few-cluster - behavior -- where the gap is materially largest -- would not fail CI. - Committing that ladder is listed as a required artifact on the - derivation row in TODO.md. + def test_se_matches_jwdid(self, fits): + """SEs reproduce Stata jwdid (reghdfe 3.2.9) at machine precision. + + The historical uniform 1.001006 gap on this G=500 arm was defect D2: + the clustered CR1 factor used only the visible treatment-cell count, + omitting the absorbed unit/time FE not nested in the unit cluster. + The mechanism was derived in closed form — + ``K_reference = cells + T`` on this no-intercept within design — and + shipped as the 3.9 K_reference convergence + (docs/methodology/variance-conventions.md); all 7 cells now sit at + ratio 1.0 with spread ~1e-14. Two-sided by design (a strict + inequality at 1.0 ± float noise is a coin flip). The subsample + LADDER (G≈20..500) gates the few-cluster behavior separately. """ etwfe, _ = fits stata = _stata_jwdid_cells(_golden()) ratios = [rec["se"] / etwfe.group_time_effects[key]["se"] for key, rec in stata.items()] assert len(ratios) == 7 - # Uniform across every cell: a finite-sample factor, not noise. - assert max(ratios) - min(ratios) < 1e-5, f"SE ratio is not uniform: {ratios}" - # Stata is LARGER (we are anti-conservative), by the amount observed on - # this G=500 panel. Regenerating against a different panel moves this. - mean_ratio = float(np.mean(ratios)) - assert mean_ratio > 1.0, "library SE should be below jwdid's, not above" - np.testing.assert_allclose(mean_ratio, 1.001006, rtol=1e-4) + assert max(ratios) - min(ratios) < 1e-9, f"SE ratio is not uniform: {ratios}" + np.testing.assert_allclose(float(np.mean(ratios)), 1.0, rtol=1e-9) class TestCallawaySantAnnaVsStataCsdid: @@ -334,9 +310,10 @@ def test_att_matches_csdid(self, fits): assert compared == 7, f"expected 7 CS cells, compared {compared}" def test_se_matches_csdid(self, fits): - """Unlike ETWFE, the CS SEs agree outright - which is what makes the - ETWFE SE gap a real finding rather than a convention difference the - whole library shares.""" + """CS SEs agree with Stata csdid outright (IF variance, no CR1 + factor). Historically this contrast is what isolated the ETWFE SE + gap as a CR1-k accounting defect rather than a library-wide + convention difference; both estimators now match their references.""" _, cs = fits stata = _stata_csdid_cells(_golden()) for key, rec in stata.items(): @@ -427,25 +404,123 @@ def test_estimation_sample_size_matches_jwdid(self, alltreated_fit): assert golden["n"] == 764 assert alltreated_fit.n_obs == golden["n"] - def test_se_gap_is_a_freshly_measured_ratio_at_this_cluster_count(self, alltreated_fit): - """The `hc1` SE gap is cluster-count dependent, so the ratio is measured - HERE rather than inherited. + def test_se_matches_jwdid_at_this_cluster_count(self, alltreated_fit): + """The G=191 arm matches jwdid at machine precision too. - The sibling arms pin 1.001006 at G=500. This arm has G=191 and shows - 1.00264201 -- copying the sibling constant would fail on first run. The - observed sequence (1.0280@G=20, 1.0132@G=40, 1.00264@G=191, - 1.0010@G=500) is monotone in G, which constrains any mechanism later - proposed for the gap (still open, TODO.md). - - As with the sibling tests, no closed form is asserted -- only the - observed factor, and that the library's SE sits BELOW jwdid's. + Pre-fix this arm showed the LARGEST pinned gap (1.00264201 at G=191 + vs 1.001006 at G=500 — the D2 defect is cluster-count dependent), + which is exactly why it gets its own gate rather than inheriting the + sibling constant. Under the K_reference convergence all 4 cells sit + at ratio 1.0 with spread ~1e-15, confirming the fix holds at a + materially smaller G than the full panel. Two-sided by design. """ stata = _stata_jwdid_cells(_golden(), "jwdid_alltreated") ratios = [ rec["se"] / alltreated_fit.group_time_effects[key]["se"] for key, rec in stata.items() ] assert len(ratios) == 4, f"expected 4 all-treated cells, got {len(ratios)}" - assert max(ratios) - min(ratios) < 1e-5, f"SE ratio is not uniform: {ratios}" - mean_ratio = float(np.mean(ratios)) - assert mean_ratio > 1.0, "library SE should be below jwdid's, not above" - np.testing.assert_allclose(mean_ratio, 1.002642, rtol=1e-4) + assert max(ratios) - min(ratios) < 1e-9, f"SE ratio is not uniform: {ratios}" + np.testing.assert_allclose(float(np.mean(ratios)), 1.0, rtol=1e-9) + + +_LADDER_RUNGS = [5, 10, 20, 40, 80, 200, 500] + + +def _ladder_subsample(df: pd.DataFrame, n_per_cohort: int) -> pd.DataFrame: + """The generator's roster rule, verbatim: the first ``n_per_cohort`` units + per ``first_treat`` cohort by ascending ``countyreal``.""" + units = df.drop_duplicates("countyreal")[["countyreal", "first_treat"]] + keep = set( + units.sort_values("countyreal").groupby("first_treat").head(n_per_cohort)["countyreal"] + ) + return df[df["countyreal"].isin(keep)].copy() + + +class TestSubsampleLadderVsStataJwdid: + """The K_reference gate at every cluster count (the ``ladder`` block). + + The D2 defect was cluster-count dependent — measured SE ratios of 1.0280 + at G=20, 1.0132 at G=40, 1.00264 at G=191, 1.0010 at G=500 — so a + full-panel gate alone could mask a partial fix that only converges at + large G. Each rung refits the library on the identical roster and pins + the ratio at 1.0 (observed spreads ~1e-15..1e-14 across G=20..500). + + The ladder doubles as the K-ACCOUNTING probe: reghdfe's own df + decomposition is recorded per rung, and ``df_a`` must equal the library's + ``absorbed_fe_cr1_k_increment`` minus the reported constant — the two + implementations agreeing on WHY the factor is what it is, not just on + the resulting number. + """ + + @pytest.fixture(scope="class") + @staticmethod + def ladder(): + # staticmethod: a class-scoped fixture defined as an instance method + # is deprecated (PytestRemovedIn10Warning). + return _golden()["ladder"] + + def test_rung_set_and_roster_rule_are_pinned(self, ladder): + """A changed roster rule or panel would move G/n; pin both.""" + assert sorted(int(k) for k in ladder["rungs"]) == _LADDER_RUNGS + expected_G = {5: 20, 10: 40, 20: 80, 40: 140, 80: 220, 200: 391, 500: 500} + df = _panel() + for n_per_cohort, rung in ((int(k), v) for k, v in ladder["rungs"].items()): + sub = _ladder_subsample(df, n_per_cohort) + assert rung["G"] == expected_G[n_per_cohort] + assert sub["countyreal"].nunique() == rung["G"] + assert len(sub) == rung["n"] + + @pytest.mark.parametrize("n_per_cohort", _LADDER_RUNGS) + def test_se_and_att_match_jwdid_at_every_cluster_count(self, ladder, n_per_cohort): + rung = ladder["rungs"][str(n_per_cohort)] + sub = _ladder_subsample(_panel(), n_per_cohort) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + fit = WooldridgeDiD(method="ols", control_group="not_yet_treated").fit( + sub, outcome="lemp", unit="countyreal", time="year", cohort="first_treat" + ) + stata = _stata_jwdid_cells({"rung": rung["cells"]}, "rung") + assert len(stata) == 7, f"expected 7 cells at G={rung['G']}, got {len(stata)}" + for key, rec in stata.items(): + np.testing.assert_allclose( + fit.group_time_effects[key]["att"], + rec["att"], + rtol=0, + atol=1e-9, + err_msg=f"ladder G={rung['G']} ATT{key} != Stata jwdid", + ) + ratios = [rec["se"] / fit.group_time_effects[key]["se"] for key, rec in stata.items()] + assert max(ratios) - min(ratios) < 1e-9, f"G={rung['G']}: SE ratio is not uniform: {ratios}" + np.testing.assert_allclose(float(np.mean(ratios)), 1.0, rtol=1e-9) + + @pytest.mark.parametrize("n_per_cohort", _LADDER_RUNGS) + def test_stata_df_accounting_matches_the_library_increment(self, ladder, n_per_cohort): + """reghdfe's df_a == absorbed_fe_cr1_k_increment - 1 at every rung. + + reghdfe reports the absorbed constant separately (``_cons`` in e(b), + ``report_constant``), so its ``df_a`` covers only the FE ranks beyond + it: df_a_initial (all FE levels) minus the cluster-nested unit FE and + the shared constant. The library folds that constant into the + increment (``has_intercept_col=False`` -> +1), hence the -1. On this + panel: increment = 1 + (G + 5 - 1) - G = 5 = T at every G, i.e. + K_reference = 7 cells + 5 = 12. + """ + from diff_diff.utils import absorbed_fe_cr1_k_increment + + rung = ladder["rungs"][str(n_per_cohort)] + sub = _ladder_subsample(_panel(), n_per_cohort) + increment = absorbed_fe_cr1_k_increment( + sub, + ["countyreal", "year"], + sub["countyreal"].to_numpy(), + has_intercept_col=False, + ) + assert rung["df_a"] == increment - 1 + # The decomposition pins WHICH FE were dropped: all G unit FE are + # cluster-nested; the redundant count adds the shared constant. + assert rung["df_a_nested"] == rung["G"] + assert rung["df_a_redundant"] == rung["G"] + 1 + assert rung["df_a_initial"] == rung["G"] + 5 + # Visible rank = the 7 treatment cells; cluster df = G-1. + assert rung["rank"] == 7 + assert rung["df_r"] == rung["G"] - 1 diff --git a/tests/test_fixest_did_twfe_parity.py b/tests/test_fixest_did_twfe_parity.py index ff81ffdd5..b7673ec98 100644 --- a/tests/test_fixest_did_twfe_parity.py +++ b/tests/test_fixest_did_twfe_parity.py @@ -95,16 +95,18 @@ def test_did_cluster_se_matches_fixest_exactly(self): np.testing.assert_allclose(res.se, exp["se"], atol=1e-10, rtol=0) def test_twfe_cluster_att_matches_fixest(self): - """The TWFE cluster-robust ATT matches fixest exactly; the SE stays - band-pinned. The residual gap is the documented fixest-CR1 ssc - convention for absorbed FE that are NOT nested in the cluster: with - ``absorb=[unit, time]`` and ``cluster=unit``, fixest counts the - non-nested time FE in the (n-1)/(n-k) denominator while the - within-transform path uses k_visible (measured ~0.25% balanced / - ~0.3% unbalanced; tracked in DEFERRED.md under "Needs external reference"). - The band pins that we never regress BEYOND the known deviation.""" + """TWFE cluster-robust ATT and SE both match fixest. + + The historical ~0.25% SE band was defect D2 (the within-transform CR1 + factor used k_visible, omitting the non-nested time FE fixest counts); + under the 3.9 K_reference convergence the ``twfe`` arm matches at + machine precision (measured 0.0 relative here, 1.9e-16 elsewhere) and + is locked exactly. The ``twfe_hetero`` arm carries a measured 5.2e-11 + relative residual (MAP-demean / BLAS-order dependent), pinned at + rtol=1e-9 (~20x headroom) rather than a machine-epsilon literal.""" golden = _load_golden() - for key in ("twfe", "twfe_hetero"): + tolerances = {"twfe": dict(atol=1e-10, rtol=0), "twfe_hetero": dict(rtol=1e-9)} + for key, tol in tolerances.items(): assert key in golden, f"required golden block {key!r} missing — regenerate the fixture" df = _build_df(golden[key]) res = TwoWayFixedEffects(vcov_type="hc1", cluster="unit").fit( @@ -112,7 +114,7 @@ def test_twfe_cluster_att_matches_fixest(self): ) exp = golden[key]["cluster_unit"] np.testing.assert_allclose(res.att, exp["att"], atol=1e-10, rtol=0) - assert res.se == pytest.approx(exp["se"], rel=0.005) + np.testing.assert_allclose(res.se, exp["se"], **tol) @_SKIP diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 8cf026ba4..eb37aeb3b 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -3362,3 +3362,305 @@ def test_conley_vcov_spot_check(self, monkeypatch): assert diag["solve_ols_fastpath"] == "chol_numpy" np.testing.assert_allclose(c_on, c_off, rtol=0, atol=1e-8) np.testing.assert_allclose(np.sqrt(np.diag(v_on)), np.sqrt(np.diag(v_off)), rtol=1e-6) + + +class TestClusterKAdjustmentSeam: + """The K_reference seam (variance-conventions.md D1/D2): front-door + validation on EVERY route, fail-closed saturation on all three sides, + weighted-lane forwarding, and the Rust-lane scalar rescale.""" + + @pytest.fixture + def clustered_data(self): + rng = np.random.default_rng(3) + n = 120 + X = np.column_stack([np.ones(n), rng.normal(size=(n, 2))]) + y = X @ np.array([1.0, 0.5, -0.2]) + rng.normal(size=n) + cl = np.repeat(np.arange(12), 10) + return X, y, cl + + # ---- front-door contracts (Verification 1b) -------------------------- + + def test_front_door_raises_without_cluster(self, clustered_data): + from diff_diff.linalg import InvalidClusterKAdjustment + + X, y, _ = clustered_data + for rv in (True, False): + with pytest.raises(InvalidClusterKAdjustment, match="cluster_ids is None"): + solve_ols(X, y, return_vcov=rv, cluster_k_adjustment=2) + + def test_front_door_raises_on_non_hc1_family(self, clustered_data): + from diff_diff.linalg import InvalidClusterKAdjustment + + X, y, cl = clustered_data + for rv in (True, False): + with pytest.raises(InvalidClusterKAdjustment, match="hc2_bm"): + solve_ols( + X, y, cluster_ids=cl, vcov_type="hc2_bm", return_vcov=rv, cluster_k_adjustment=2 + ) + + def test_front_door_requires_int(self, clustered_data): + from diff_diff.linalg import InvalidClusterKAdjustment + + X, y, cl = clustered_data + with pytest.raises(InvalidClusterKAdjustment, match="must be an int"): + solve_ols(X, y, cluster_ids=cl, cluster_k_adjustment=2.5) # type: ignore[arg-type] + with pytest.raises(InvalidClusterKAdjustment, match="must be an int"): + solve_ols(X, y, cluster_ids=cl, cluster_k_adjustment=True) # type: ignore[arg-type] + # Non-int ZEROS must fail too: `0.0 == 0` and `False == 0`, so a zero + # fast path ahead of the type check would silently accept them. + with pytest.raises(InvalidClusterKAdjustment, match="must be an int"): + solve_ols(X, y, cluster_ids=cl, cluster_k_adjustment=0.0) # type: ignore[arg-type] + with pytest.raises(InvalidClusterKAdjustment, match="must be an int"): + solve_ols(X, y, cluster_ids=cl, cluster_k_adjustment=False) # type: ignore[arg-type] + + def test_compute_robust_vcov_mirrors_the_contract(self, clustered_data): + from diff_diff.linalg import InvalidClusterKAdjustment + + X, y, cl = clustered_data + residuals = y - X @ np.linalg.lstsq(X, y, rcond=None)[0] + with pytest.raises(InvalidClusterKAdjustment): + compute_robust_vcov(X, residuals, cluster_k_adjustment=2) + with pytest.raises(InvalidClusterKAdjustment): + compute_robust_vcov(X, residuals, cl, vcov_type="hc2_bm", cluster_k_adjustment=2) + + def test_linear_regression_fit_raises_symmetrically(self, clustered_data): + from diff_diff.linalg import InvalidClusterKAdjustment + + X, y, _ = clustered_data + with pytest.raises(InvalidClusterKAdjustment): + LinearRegression(vcov_type="classical", include_intercept=False).fit( + X, y, cluster_k_adjustment=2 + ) + + def test_type_contract_is_not_route_dependent(self, clustered_data): + """Non-int adjustments are rejected on the SHORT-CIRCUITED routes too: + a classical fit (truthiness gate, never reaches solve_ols's front + door with the kwarg) and a wild_bootstrap_se call whose degenerate + early-returns fire before the clustered solve (both would previously + swallow `0.0` / `False` silently).""" + from diff_diff.linalg import InvalidClusterKAdjustment + from diff_diff.utils import wild_bootstrap_se + + X, y, cl = clustered_data + for bad in (0.0, False): + with pytest.raises(InvalidClusterKAdjustment, match="must be an int"): + LinearRegression(vcov_type="classical", include_intercept=False).fit( + X, y, cluster_k_adjustment=bad # type: ignore[arg-type] + ) + # a saturated design (n == k) hits the degenerate return before the + # clustered solve; the entry type check must fire first + n_sat = X.shape[1] + with pytest.raises(InvalidClusterKAdjustment, match="must be an int"): + wild_bootstrap_se( + X[:n_sat], + y[:n_sat], + np.zeros(n_sat), + cl[:n_sat], + 0, + cluster_k_adjustment=1.5, # type: ignore[arg-type] + ) + + def test_valid_path_coefficient_only(self, clustered_data): + """return_vcov=False + nonzero adjustment succeeds (None passthrough).""" + X, y, cl = clustered_data + coef, resid, vcov = solve_ols( + X, y, cluster_ids=cl, return_vcov=False, cluster_k_adjustment=3 + ) + assert vcov is None + assert np.all(np.isfinite(coef)) + + def test_kwarg_is_keyword_only_and_positional_slots_are_unchanged(self, clustered_data): + """Signature contract: cluster_k_adjustment is KEYWORD-ONLY on both + public carriers and no pre-existing positional slot moved — a legacy + positional call can never bind another argument to it. + compute_robust_vcov accepts positionals through return_dof exactly as + before; solve_ols takes only (X, y) positionally.""" + import inspect + + from diff_diff.linalg import compute_robust_vcov + + for fn in (solve_ols, compute_robust_vcov): + param = inspect.signature(fn).parameters["cluster_k_adjustment"] + assert param.kind is inspect.Parameter.KEYWORD_ONLY + assert param.default == 0 + positional = [ + n + for n, p in inspect.signature(compute_robust_vcov).parameters.items() + if p.kind is inspect.Parameter.POSITIONAL_OR_KEYWORD + ] + assert positional == [ + "X", + "residuals", + "cluster_ids", + "weights", + "weight_type", + "vcov_type", + "return_dof", + ] + # and a legacy full-positional call still binds correctly + X, y, cl = clustered_data + residuals = y - X @ np.linalg.lstsq(X, y, rcond=None)[0] + legacy = compute_robust_vcov(X, residuals, cl, None, "pweight", "hc1", False) + kw = compute_robust_vcov(X, residuals, cluster_ids=cl, vcov_type="hc1") + np.testing.assert_array_equal(legacy, kw) + + # ---- the factor itself ---------------------------------------------- + + def test_adjustment_is_the_exact_scalar_factor(self, clustered_data): + """vcov(adj) == vcov(0) * (n-k)/(n-k-adj) on every lane, both signs.""" + X, y, cl = clustered_data + n, k = X.shape + _, _, v0 = solve_ols(X, y, cluster_ids=cl) + for adj in (5, -1): + _, _, v = solve_ols(X, y, cluster_ids=cl, cluster_k_adjustment=adj) + expect = (n - k) / (n - k - adj) + np.testing.assert_allclose(v / v0, expect, rtol=0, atol=1e-12) + + def test_zero_adjustment_bit_identical(self, clustered_data): + X, y, cl = clustered_data + _, _, v0 = solve_ols(X, y, cluster_ids=cl) + _, _, v1 = solve_ols(X, y, cluster_ids=cl, cluster_k_adjustment=0) + assert np.array_equal(v0, v1) + + def test_tail_dof_vec_stays_on_visible_k(self, clustered_data): + """The adjustment never moves the reported dof (tail df is PR C).""" + X, y, cl = clustered_data + residuals = y - X @ np.linalg.lstsq(X, y, rcond=None)[0] + from diff_diff.linalg import _compute_robust_vcov_numpy + + _, dof0 = _compute_robust_vcov_numpy(X, residuals, cl, return_dof=True) + _, dof5 = _compute_robust_vcov_numpy( + X, residuals, cl, return_dof=True, cluster_k_adjustment=5 + ) + np.testing.assert_array_equal(dof0, dof5) + + # ---- fail-closed saturation, all three sides ------------------------- + + def test_fail_closed_visible_saturation_survives_negative_adjustment(self): + """n == k_visible -> all-NaN on BOTH backends even when the corrected + denominator would be positive (residuals are identically zero; the + adjudicated saturation semantics).""" + rng = np.random.default_rng(11) + n = 6 + X = np.column_stack([np.eye(3)[np.arange(n) % 3], rng.normal(size=(n, 3))]) + assert X.shape == (6, 6) + y = rng.normal(size=n) + cl = np.repeat([0, 1], 3) + coef, resid, vcov = solve_ols( + X, y, cluster_ids=cl, cluster_k_adjustment=-2, rank_deficient_action="silent" + ) + assert vcov is not None and np.isnan(vcov).all() + + def test_fail_closed_k_inf_sides(self, clustered_data): + from diff_diff.linalg import _compute_robust_vcov_numpy + + X, y, cl = clustered_data + n, k = X.shape + residuals = y - X @ np.linalg.lstsq(X, y, rcond=None)[0] + # n - k_inf <= 0 + v_hi = _compute_robust_vcov_numpy(X, residuals, cl, cluster_k_adjustment=n) + assert np.isnan(v_hi).all() + # k_inf <= 0 + v_lo = _compute_robust_vcov_numpy(X, residuals, cl, cluster_k_adjustment=-k) + assert np.isnan(v_lo).all() + + def test_fail_closed_k_inf_exact_boundary(self, clustered_data): + """Both k_inf guard sides pinned at their EXACT boundary, with the + one-step recovery on each: n - k_inf == 0 is NaN while == 1 is + finite, and k_inf == 0 is NaN while == 1 is finite. The deep-past- + boundary cases above cannot distinguish an off-by-one guard.""" + from diff_diff.linalg import _compute_robust_vcov_numpy + + X, y, cl = clustered_data + n, k = X.shape + residuals = y - X @ np.linalg.lstsq(X, y, rcond=None)[0] + # saturation side: k_inf == n -> NaN; k_inf == n - 1 -> finite + assert np.isnan( + _compute_robust_vcov_numpy(X, residuals, cl, cluster_k_adjustment=n - k) + ).all() + assert np.isfinite( + _compute_robust_vcov_numpy(X, residuals, cl, cluster_k_adjustment=n - k - 1) + ).all() + # deflation side: k_inf == 0 -> NaN; k_inf == 1 -> finite + assert np.isnan(_compute_robust_vcov_numpy(X, residuals, cl, cluster_k_adjustment=-k)).all() + assert np.isfinite( + _compute_robust_vcov_numpy(X, residuals, cl, cluster_k_adjustment=-(k - 1)) + ).all() + + # ---- weighted external-vcov forwarding sites ------------------------- + + def test_weighted_full_rank_forwarding(self, clustered_data): + """The weighted external-vcov lane bypasses _solve_ols_numpy; the + adjustment must reach its direct kernel call (full-rank site).""" + X, y, cl = clustered_data + n, k = X.shape + w = np.ones(n) + w[:10] = 2.0 + _, _, v0 = solve_ols(X, y, cluster_ids=cl, weights=w) + _, _, v5 = solve_ols(X, y, cluster_ids=cl, weights=w, cluster_k_adjustment=5) + expect = (n - k) / (n - k - 5) + np.testing.assert_allclose(v5 / v0, expect, rtol=0, atol=1e-12) + + def test_weighted_rank_deficient_forwarding(self, clustered_data): + """The reduced-design weighted site (collinear column dropped).""" + X, y, cl = clustered_data + n = X.shape[0] + X_bad = np.column_stack([X, X[:, 1] + X[:, 2]]) # collinear 4th col + w = np.ones(n) + w[:10] = 2.0 + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + _, _, v0 = solve_ols( + X_bad, y, cluster_ids=cl, weights=w, rank_deficient_action="silent" + ) + _, _, v5 = solve_ols( + X_bad, + y, + cluster_ids=cl, + weights=w, + rank_deficient_action="silent", + cluster_k_adjustment=5, + ) + kept = ~np.isnan(np.diag(v0)) + k_red = int(kept.sum()) + expect = (n - k_red) / (n - k_red - 5) + np.testing.assert_allclose( + v5[np.ix_(kept, kept)] / v0[np.ix_(kept, kept)], expect, rtol=0, atol=1e-12 + ) + + # ---- Rust lanes ------------------------------------------------------ + + def test_rust_instability_fallback_forwards_adjustment(self, clustered_data): + """The documented Rust-instability numpy re-run must carry the + adjustment (the :2199-class forwarding site).""" + from unittest.mock import patch + + if not HAS_RUST_BACKEND: + pytest.skip("Rust backend not available") + X, y, cl = clustered_data + n, k = X.shape + residuals = y - X @ np.linalg.lstsq(X, y, rcond=None)[0] + v0 = compute_robust_vcov(X, residuals, cl) + + def mock_rust_vcov(*args, **kwargs): + raise ValueError("Matrix inversion numerically unstable") + + with patch("diff_diff.linalg._rust_compute_robust_vcov", mock_rust_vcov): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + v5 = compute_robust_vcov(X, residuals, cl, cluster_k_adjustment=5) + expect = (n - k) / (n - k - 5) + np.testing.assert_allclose(v5 / v0, expect, rtol=0, atol=1e-10) + + def test_rust_and_numpy_agree_on_adjusted_vcov(self, clustered_data, monkeypatch): + """Cross-backend parity for the corrected clustered vcov (the scalar + rescale on the Rust lane vs the in-kernel factor on numpy).""" + if not HAS_RUST_BACKEND: + pytest.skip("Rust backend not available") + import diff_diff.linalg as lmod + + X, y, cl = clustered_data + _, _, v_rust = solve_ols(X, y, cluster_ids=cl, cluster_k_adjustment=5) + monkeypatch.setattr(lmod, "HAS_RUST_BACKEND", False) + _, _, v_py = solve_ols(X, y, cluster_ids=cl, cluster_k_adjustment=5) + np.testing.assert_allclose(v_rust, v_py, rtol=1e-12) diff --git a/tests/test_methodology_sun_abraham.py b/tests/test_methodology_sun_abraham.py index 9df29722f..66ca6b0a4 100644 --- a/tests/test_methodology_sun_abraham.py +++ b/tests/test_methodology_sun_abraham.py @@ -8,9 +8,11 @@ - classical → R lm() summary at atol=1e-10 (full-dummy) - hc2 → R sandwich::vcovHC(type="HC2") at atol=1e-10 (full-dummy) - hc2_bm → R clubSandwich::vcovCR(cluster=..., type="CR2") + coef_test$df_Satt at atol=1e-10 (full-dummy) -- hc1 event-study e=0 → R fixest::sunab + cluster=~unit at atol=5e-3 - (within-transform; documented HC1 finite-sample-correction deviation, - see REGISTRY.md SunAbraham section and DEFERRED.md row tracking the gap) +- hc1 event-study e=0 → R fixest::sunab + cluster=~unit at rtol=1e-12 + (within-transform; the historical ~5e-3 finite-sample-correction gap was + defect D2 of the 3.9 variance program — the CR1 factor omitted the absorbed + unit+time FE — closed by the K_reference convergence; measured residual + 3.1e-14 relative, pinned with ~30x headroom) """ import json @@ -255,11 +257,8 @@ def test_sa_hc1_event_study_e0_se_close_to_fixest_sunab(self, sa_panel, sa_golde r_se = sa_golden["sunab_hc1_event_study_e0_se"] if r_se is None or (isinstance(r_se, float) and np.isnan(r_se)): pytest.skip("fixest::sunab parity SE missing from golden") - # Absolute tolerance pin matching the documented deviation; harmonized - # with REGISTRY.md SunAbraham deviation note and DEFERRED.md row. - assert abs(sa_se - r_se) < 5e-3, ( - f"SA hc1 event-study e=0 SE {sa_se:.6f} diverges from fixest::sunab " - f"event-study SE {r_se:.6f} by {abs(sa_se - r_se):.2e} > 5e-3. " - "Expected ~2e-3 from finite-sample correction (see REGISTRY.md " - "deviation note); investigate if larger." - ) + # Machine-precision pin under the 3.9 K_reference convergence + # (measured residual 3.1e-14 relative; ~30x headroom). The historical + # ~5e-3 band was the D2 finite-sample-correction gap, now closed — + # see the REGISTRY.md SunAbraham HC1 note. + np.testing.assert_allclose(sa_se, r_se, rtol=1e-12) diff --git a/tests/test_methodology_twfe.py b/tests/test_methodology_twfe.py index 76d214cb7..be5fdd207 100644 --- a/tests/test_methodology_twfe.py +++ b/tests/test_methodology_twfe.py @@ -469,13 +469,13 @@ def test_att_matches_r_twfe(self, r_twfe_results, r_benchmark_panel_data): ) def test_se_matches_r_twfe(self, r_twfe_results, r_benchmark_panel_data): - """Cluster-robust SE within the documented fixest-CR1 band. + """Cluster-robust SE matches fixest under the K_reference convergence. - Measured gap is ~2.5e-3 relative — the documented ~0.25% - absorbed-FE CR1 DOF-convention deviation (fixest counts non-nested - FE in the finite-sample denominator; see SE-audit G2 + the - band-pin in tests/test_fixest_did_twfe_parity.py). rtol=0.005 - pins that we never regress BEYOND the known band.""" + Measured gap is ~1.4e-15 relative — the historical ~2.5e-3 band was + defect D2 (the CR1 factor omitted the non-nested time FE fixest + counts), closed in 3.9. Pinned at rtol=1e-9 rather than a + machine-epsilon literal because this MAP-demean lane measured up to + 5.2e-11 on other fixtures (BLAS-order dependent).""" data, _ = r_benchmark_panel_data py_results = self._run_python_twfe(data) @@ -483,7 +483,7 @@ def test_se_matches_r_twfe(self, r_twfe_results, r_benchmark_panel_data): np.testing.assert_allclose( py_results.se, r_twfe_results["se"], - rtol=0.005, + rtol=1e-9, err_msg=f"SE mismatch: Python={py_results.se:.6f}, R={r_twfe_results['se']:.6f}", ) @@ -574,14 +574,14 @@ def test_moderate_t_pins_residual_df_convention(self): ), "p-value matches the cluster-df convention — the documented deviation flipped" def test_ci_matches_r_twfe(self, r_twfe_results, r_benchmark_panel_data): - """CI bounds within the documented-deviation band (measured ~1.9e-3). - - The bounds inherit BOTH documented deviations through - att ± crit*se: the ~0.25% CR1 SE finite-sample-scale band AND the - clustered inference-df convention (Python t(residual df=148) - critical value vs fixest t(G−1=49); REGISTRY - §TwoWayFixedEffects). The ATT itself matches at machine - precision.""" + """CI bounds within the remaining df-convention band (measured ~1.7e-3). + + Post-K_reference the SE side matches fixest at machine precision, so + the bounds' residual gap through att ± crit*se is SOLELY the clustered + inference-df convention (Python t(residual df=148) critical value vs + fixest t(G−1=49); REGISTRY §TwoWayFixedEffects — PR C scope). + Tightened 0.005 → 0.003 accordingly (measured 1.7e-3/1.4e-3). The + ATT itself matches at machine precision.""" data, _ = r_benchmark_panel_data py_results = self._run_python_twfe(data) @@ -589,13 +589,13 @@ def test_ci_matches_r_twfe(self, r_twfe_results, r_benchmark_panel_data): np.testing.assert_allclose( py_results.conf_int[0], r_twfe_results["ci_lower"], - rtol=0.005, + rtol=0.003, err_msg=f"CI lower mismatch: Python={py_results.conf_int[0]:.6f}, R={r_twfe_results['ci_lower']:.6f}", ) np.testing.assert_allclose( py_results.conf_int[1], r_twfe_results["ci_upper"], - rtol=0.005, + rtol=0.003, err_msg=f"CI upper mismatch: Python={py_results.conf_int[1]:.6f}, R={r_twfe_results['ci_upper']:.6f}", ) @@ -618,8 +618,9 @@ def test_att_matches_r_with_covariate( def test_se_matches_r_with_covariate( self, r_twfe_results_with_covariate, r_benchmark_panel_data_with_covariate ): - """SE with covariate within the documented CR1 band (measured - ~2.5e-3; see test_se_matches_r_twfe).""" + """SE with covariate matches fixest under the K_reference convergence + (the historical ~2.5e-3 D2 band is closed; see + test_se_matches_r_twfe for the rtol=1e-9 headroom rationale).""" data, _ = r_benchmark_panel_data_with_covariate py_results = self._run_python_twfe(data, covariates=["x1"]) @@ -627,7 +628,7 @@ def test_se_matches_r_with_covariate( np.testing.assert_allclose( py_results.se, r_twfe_results_with_covariate["se"], - rtol=0.005, + rtol=1e-9, err_msg=f"SE w/ cov mismatch: Python={py_results.se:.6f}, R={r_twfe_results_with_covariate['se']:.6f}", ) @@ -892,12 +893,16 @@ def test_df_adjustment_for_absorbed_fe(self): df_adjustment = n_units + n_times - 2 cluster_ids = data["unit"].values + # Clustered-CR1 K_reference: the manual reference carries the same + # non-nested absorbed rank (time | unit) as the estimator — with the + # unit cluster, the unit FE is nested (dropped) and the time FE adds + # its conditional rank n_times - 1. reg = LinearRegression( include_intercept=False, robust=True, cluster_ids=cluster_ids, rank_deficient_action="silent", - ).fit(X, y, df_adjustment=df_adjustment) + ).fit(X, y, df_adjustment=df_adjustment, cluster_k_adjustment=n_times - 1) manual_se = reg.get_inference(1).se np.testing.assert_allclose( @@ -1114,13 +1119,16 @@ def test_cluster_se_differs_from_hc1_se(self): f"HC1 SE ({hc1_se:.6f}) — auto-clustering must be active" ) - # Also verify TWFE SE matches a manually computed cluster SE + # Also verify TWFE SE matches a manually computed cluster SE. The + # manual reference carries the clustered-CR1 K_reference increment: + # under the unit cluster the unit FE is nested (dropped) and the + # time FE contributes its conditional rank n_times - 1. cluster_reg = LinearRegression( include_intercept=False, robust=True, cluster_ids=data["unit"].values, rank_deficient_action="silent", - ).fit(X, y, df_adjustment=df_adjustment) + ).fit(X, y, df_adjustment=df_adjustment, cluster_k_adjustment=n_times - 1) manual_cluster_se = cluster_reg.get_inference(1).se np.testing.assert_allclose( diff --git a/tests/test_methodology_wooldridge.py b/tests/test_methodology_wooldridge.py index fdcfd9684..589deb5f3 100644 --- a/tests/test_methodology_wooldridge.py +++ b/tests/test_methodology_wooldridge.py @@ -38,8 +38,10 @@ The ``hc1`` variant is NOT pinned against R in ``TestWooldridgeParityR`` because the diff-diff within-transform finite-sample correction -``(n-1)/(n-k_dm)`` differs from ``lm + clubSandwich::vcovCR(type="CR1S")``'s -``(n-1)/(n-k_total)`` correction; see ``docs/methodology/REGISTRY.md`` +``(n-1)/(n-K_reference)`` (cells + constant + non-nested time-FE rank, 3.9 +K_reference convergence — matches Stata jwdid/reghdfe) differs from +``lm + clubSandwich::vcovCR(type="CR1S")``'s ``(n-1)/(n-k_total)`` correction +(all columns, nested unit dummies included); see ``docs/methodology/REGISTRY.md`` "Variance families" → "Deviation from R" for the algebra. The hc1 path is locked instead by ``tests/test_wooldridge.py::TestWooldridgeVcovType::test_hc1_se_bit_equal_to_pre_pr_baseline`` @@ -2141,8 +2143,10 @@ class TestW2025LibraryDeviations: Section 8 gap closes via Stage C). The 5 substantive deviations locked here: - 1. **HC1 finite-sample correction** ``(n-1)/(n-k_dm)`` (within-transform) - vs R ``lm + CR1S`` ``(n-1)/(n-k_total)``. + 1. **HC1 finite-sample correction** ``(n-1)/(n-K_reference)`` + (within-transform, K_reference = cells + constant + non-nested + time-FE rank; matches jwdid/reghdfe) vs R ``lm + CR1S`` + ``(n-1)/(n-k_total)`` (nested unit dummies included). 2. **QMLE sandwich `(G/(G-1)) × ((n-1)/(n-k))`** vs Stata ``jwdid`` ``G/(G-1)`` only (logit/Poisson paths). 3. **Nonlinear methods via direct QMLE** vs R ``etwfe`` fixest backend. @@ -2154,13 +2158,14 @@ class TestW2025LibraryDeviations: """ def test_hc1_within_transform_se_differs_from_naive_full_design(self) -> None: - """Deviation 1: ``vcov_type='hc1'`` uses within-transform ``(n-1)/(n-k_dm)``. + """Deviation 1: ``vcov_type='hc1'`` uses within-transform ``(n-1)/(n-K_reference)``. - The within-transformed design has ``k_dm`` columns (only the - treatment-cell + covariate columns); R's ``lm + CR1S`` on the - full-dummy design has ``k_total`` columns (including all unit - + time dummies). The two SEs differ by the factor - ``sqrt((n-k_total) / (n-k_dm))``. + The within-transformed design carries ``K_reference`` in the CR1 + factor (treatment cells + the absorbed constant + the non-nested + time-FE rank; 3.9 convergence); R's ``lm + CR1S`` on the full-dummy + design counts all ``k_total`` columns (nested unit dummies + included). The two SEs differ by the factor + ``sqrt((n-k_total) / (n-K_reference))``. This test verifies the library uses the within-transform factor (the documented deviation) by comparing ``vcov_type='hc1'`` SE @@ -2197,7 +2202,7 @@ def test_hc1_within_transform_se_differs_from_naive_full_design(self) -> None: assert se_hc1 != se_classical, ( f"hc1 SE = classical SE for (g={sample_key[0]}, t={sample_key[1]}). " f"This contradicts REGISTRY Deviation #4 (HC1 uses within-transform " - f"(n-1)/(n-k_dm) finite-sample factor)." + f"(n-1)/(n-K_reference) finite-sample factor)." ) def test_qmle_sandwich_inflates_se_vs_stata_jwdid_reference(self) -> None: @@ -2248,8 +2253,9 @@ def test_nonlinear_methods_use_direct_qmle_not_fixest_backend(self) -> None: R ``etwfe`` uses ``fixest`` for nonlinear paths; the library uses direct QMLE via ``compute_robust_vcov`` to avoid a statsmodels/ - fixest dependency. This results in HC1 finite-sample factor - ``(n-1)/(n-k_dm)`` rather than fixest's ``(n-1)/(n-k_total)``. + fixest dependency. On the OLS path the HC1 finite-sample factor is + ``(n-1)/(n-K_reference)`` (3.9); the QMLE paths keep their own + documented factor (deviation #2). Lock by verifying logit + Poisson runs without statsmodels/fixest installed — if either dep had crept in, an ``ImportError`` would diff --git a/tests/test_v4_matrix.py b/tests/test_v4_matrix.py index b5fcb5e79..fcac5de0f 100644 --- a/tests/test_v4_matrix.py +++ b/tests/test_v4_matrix.py @@ -115,11 +115,12 @@ # completeness sweep over public FUNCTIONS and the dCDH results mirror # (M-097..M-115) = 99, plus Phase 2b PR 1's two rows (M-117, M-122) = 101, # plus the ETWFE reference-period pair (M-123, M-124) = 103, plus the -# comparison-support row (M-125) = 104. +# comparison-support row (M-125) = 104, plus the clustered-CR1 K_reference +# convergence row (M-126) = 105. # Ids are never reused and terminal rows are never # deleted, so the ledger only grows - raise the floor when rows are added; a # lower parse count means scanner/format drift or an illegal row deletion. -ROW_COUNT_FLOOR = 104 +ROW_COUNT_FLOOR = 105 # Committed snapshot of the shipped id set ("ids are never deleted or reused" # contract - a delete-one-add-one edit keeps the count above the floor but trips @@ -133,7 +134,9 @@ # dCDH results mirror and the fourth `robust` site; (117,117)/(122,122) = # Phase 2b PR 1 (balance_e onto aggregate(), AggregationResult); (123,124) = # the ETWFE reference-period fix (unidentified-cohort exclusion, and the -# fail-closed guard on a design with no estimable post-treatment cell). +# fail-closed guard on a design with no estimable post-treatment cell); +# (126,126) = the clustered-CR1 K_reference convergence (D1+D2 of the 3.9 +# variance-consolidation program). # M-116 and # M-118..M-121 are reserved for the later 2b PRs, not deleted - ids are # never reused, so a gap here is intentional. @@ -152,6 +155,7 @@ (117, 117), (122, 124), (125, 125), + (126, 126), ] EXPECTED_INITIAL_IDS = frozenset( f"M-{n:03d}" for lo, hi in _INITIAL_ID_RANGES for n in range(lo, hi + 1) @@ -548,12 +552,13 @@ def test_initial_ids_never_deleted(): """The shipped id set is immutable: ids are never deleted or reused (spec section 11). ROW_COUNT_FLOOR alone would let a delete-one-add-one edit pass; this snapshot cannot. - Extends as rows ship (104 as of the ETWFE comparison-support row: Phase 1 + - diagnostic-family + M-092/M-093 + M-094..M-096 + the M-097..M-115 - public-function completeness sweep + M-117/M-122 + M-123/M-124 + M-125).""" + Extends as rows ship (105 as of the clustered-CR1 K_reference convergence row: + Phase 1 + diagnostic-family + M-092/M-093 + M-094..M-096 + the M-097..M-115 + public-function completeness sweep + M-117/M-122 + M-123/M-124 + M-125 + + M-126).""" missing = sorted(EXPECTED_INITIAL_IDS - set(_ROW_IDS)) assert not missing, f"ledger rows deleted (ids are permanent): {missing}" - assert len(EXPECTED_INITIAL_IDS) == 104 + assert len(EXPECTED_INITIAL_IDS) == 105 def test_version_tuple_pads_to_three_components(): diff --git a/tests/test_variance_conventions.py b/tests/test_variance_conventions.py index 849a746b6..3e5536603 100644 --- a/tests/test_variance_conventions.py +++ b/tests/test_variance_conventions.py @@ -100,6 +100,10 @@ class Capture: def __init__(self, monkeypatch): self.cr1_k: list = [] self.tail_df: list = [] + # Per-call K_reference adjustments on the clustered hc1 lane, so a + # row can pin the increment itself even where the exact k is a + # contract detail. + self.cr1_adjustments: list = [] # Clustered calls in a NON-hc1 family: (vcov_type, k). The CR1 factor # this matrix audits is the clustered hc1 denominator, so a surface # that silently switches clustered family must fail its row even when @@ -117,7 +121,12 @@ def spy_vcov(*a, **k): b.apply_defaults() if b.arguments.get("cluster_ids") is not None: if b.arguments.get("vcov_type") == "hc1": - self.cr1_k.append(int(b.arguments["X"].shape[1])) + # cr1_k records the K_REFERENCE count reaching the CR1 + # denominator: visible columns + the signed + # cluster_k_adjustment (variance-conventions.md D1/D2). + _adj = int(b.arguments.get("cluster_k_adjustment") or 0) + self.cr1_k.append(int(b.arguments["X"].shape[1]) + _adj) + self.cr1_adjustments.append(_adj) else: self.unexpected_clustered.append( (str(b.arguments.get("vcov_type")), int(b.arguments["X"].shape[1])) @@ -175,22 +184,26 @@ def snapshot(self): fit=lambda df: diff_diff.DifferenceInDifferences(cluster="unit").fit( df, outcome="y", treatment="grp", time="post", absorb=["unit", "time"] ), - cr1_k=(2,), + cr1_k=(7,), tail_df=(294.0,), - status="defect", - reason="D2: CR1 k omits absorbed FE not nested in the cluster (time)", + status="legitimate", + reason=( + "K_reference (D2 fixed): k = 2 visible + rank(time given unit) = 5; " + "matches reghdfe/fixest ssc(K.fixef='nested')" + ), ), dict( key="did_fixed_effects_hc1_cluster_unit", fit=lambda df: diff_diff.DifferenceInDifferences(cluster="unit").fit( df, outcome="y", treatment="grp", time="post", fixed_effects=["unit", "time"] ), - cr1_k=(66,), + cr1_k=(7,), tail_df=(294.0,), - status="defect", + status="legitimate", reason=( - "D1: same model as did_absorb yet k=66 vs 2 -> SEs differ 10.35%; " - "full-dummy k also counts the cluster-nested unit FE the references drop" + "K_reference (D1 fixed): 66 visible minus the 59 cluster-nested " + "unit dummies -> identical SE to did_absorb (documented deviation " + "from a literal explicit-dummy R comparison, which counts all 66)" ), ), dict( @@ -208,30 +221,41 @@ def snapshot(self): fit=lambda df: diff_diff.TwoWayFixedEffects(vcov_type="hc1", cluster="unit").fit( df, outcome="y", treatment="grp", time="post", unit="unit" ), - cr1_k=(2,), + cr1_k=(3,), tail_df=(298.0,), - status="defect", - reason="D2 (within-transform k_visible); tail df is residual n-K_full", + status="legitimate", + reason=( + "K_reference (D2 fixed): 2 visible + rank(post given unit) = 1; " + "matches fixest cluster arm at rel 0 (committed golden)" + ), ), dict( key="wooldridge_hc1_within", fit=lambda df: diff_diff.WooldridgeDiD(method="ols").fit( df, outcome="y", unit="unit", time="time", cohort="first_treat" ), - cr1_k=(9,), + cr1_k=(15,), tail_df=(None,) * 10, status="defect", - reason="D2 (k_visible=cells only) + normal-theory tail df with no df_convention knob", + reason=( + "CR1 k converged on K_reference (D2 fixed: 9 cells + T = 15, " + "no intercept col -> +1 term; jwdid arms at ratio 1.0); tail df " + "is still normal theory with no df_convention knob (PR C)" + ), ), dict( key="sun_abraham_hc1", fit=lambda df: diff_diff.SunAbraham().fit( df, outcome="y", unit="unit", time="time", first_treat="first_treat" ), - cr1_k=(15,), + cr1_k=(21,), tail_df=(280.0,) * 15 + (None,) * 8, status="defect", - reason="D2 + D4: residual df per cohort-period cell but normal theory on aggregates", + reason=( + "CR1 k converged on K_reference (D2 fixed: 15 cells + 6, no " + "intercept col; fixest sunab parity ~5e-15); D4 remains: residual " + "df per cohort-period cell but normal theory on aggregates (PR C)" + ), ), dict( key="stacked_did_hc1", @@ -257,6 +281,77 @@ def snapshot(self): status="legitimate", reason="L2: G-1 tail df (Stata/fixest convention) — the convergence target", ), + dict( + key="mpd_absorb_hc1_cluster_unit", + fit=lambda df: diff_diff.MultiPeriodDiD(cluster="unit").fit( + df, + outcome="y", + treatment="grp", + time="time", + absorb=["unit", "time"], + reference_period=1, + ), + cr1_k=(11,), + tail_df=(290.0,) * 6, + status="legitimate", + reason=( + "K_reference: 6 visible + rank(time given unit) = 5; equals the " + "fixed_effects form's 70 - 59 (MPD absorb/fixed_effects equivalence)" + ), + ), + dict( + key="mpd_fixed_effects_hc1_cluster_unit", + fit=lambda df: diff_diff.MultiPeriodDiD(cluster="unit").fit( + df, + outcome="y", + treatment="grp", + time="time", + fixed_effects=["unit", "time"], + reference_period=1, + ), + cr1_k=(11,), + tail_df=(290.0,) * 6, + status="legitimate", + reason=( + "K_reference: 70 visible (incl. built-in period dummies, MPD's " + "time-FE block) minus the 59 cluster-nested unit dummies = 11 — " + "identical to the absorb form" + ), + ), + dict( + key="mpd_plain_hc1_cluster_time", + fit=lambda df: diff_diff.MultiPeriodDiD(cluster="time").fit( + df, outcome="y", treatment="grp", time="time", reference_period=1 + ), + cr1_k=(7,), + tail_df=(348.0,) * 6, + status="legitimate", + reason=( + "the NESTED orientation of the built-in period dummies: 12 " + "visible minus their rank 5 under a time cluster (under-" + "subtraction is caught here; the unit-cluster rows catch over-)" + ), + ), + dict( + key="lpdid_absorb_nested_cluster_grp", + fit=lambda df: diff_diff.LPDiD(pre_window=2, post_window=2, cluster="grp").fit( + df.assign(region=df["unit"] % 2), + outcome="y", + unit="unit", + time="time", + treatment="treated", + absorb=["region"], + ), + cr1_k=(4, 4, 5, 5, 5, 6), + tail_df=(1.0,) * 6, + status="legitimate", + reason=( + "LPDiD absorb dummies nested in the cluster subtract their rank " + "(adj -1 per horizon: region == grp here); _event_time stays " + "counted (unit-level cluster does not nest time); G-1 tail df (L2)" + ), + expected_adjustment=-1, + ), dict( key="imputation_default", fit=lambda df: diff_diff.ImputationDiD().fit( @@ -277,10 +372,17 @@ def snapshot(self): first_treat="first_treat", aggregate="event_study", ), - cr1_k=None, # resolved at collection: must be NON-empty (D2 applies here) + cr1_k=None, # contract row: shared CR1 REACHED; exact k is config-detail tail_df=None, + # The increment IS pinned even though k is not: every clustered-hc1 + # call on this surface must carry the [time, unit] no-intercept + # increment (1 + 65 - 60 = 6 on the audit panel). + expected_adjustment=6, status="defect", - reason="pretrends lead regression runs the shared clustered CR1 with k_visible", + reason=( + "pretrends lead regression: CR1 k converged on K_reference " + "(D2 fixed); normal-theory tail df remains (PR C family)" + ), ), dict( key="two_stage_default", @@ -342,6 +444,16 @@ def _assert_row(row, monkeypatch): ) else: assert cr1_k == row["cr1_k"], f"{row['key']}: cr1_k {cr1_k} != {row['cr1_k']}" + if row.get("expected_adjustment") is not None: + # Contract rows pin the K_reference increment itself even where the + # exact k is configuration-detail: every clustered-hc1 call must + # carry exactly the expected adjustment. + assert cap.cr1_adjustments and all( + a == row["expected_adjustment"] for a in cap.cr1_adjustments + ), ( + f"{row['key']}: expected every cluster_k_adjustment == " + f"{row['expected_adjustment']}, captured {cap.cr1_adjustments}" + ) if row["tail_df"] is not None: assert tail_df == row["tail_df"], f"{row['key']}: tail_df {tail_df} != {row['tail_df']}" # Self-check: a row claiming tail-df expectations must actually capture some @@ -390,8 +502,14 @@ def test_capture_flags_non_hc1_clustered_family(monkeypatch): assert all(fam == "hc2_bm" for fam, _ in cap.unexpected_clustered) -def test_d1_divergence_is_pinned(): - """The absorb-vs-fixed_effects SE split: same model, same ATT, k=2 vs 66.""" +def test_d1_convergence_is_pinned(): + """absorb= and fixed_effects= now return the IDENTICAL clustered SE. + + Before the K_reference fix the two documented-equivalent idioms differed + by exactly sqrt((360-2)/(360-66)) = 10.35% (defect D1). Both lanes now + land on K_reference = 7, and the common SE equals the old absorb-side + value rescaled by the K change — the externally anchored prediction. + """ df = make_panel() with warnings.catch_warnings(): warnings.simplefilter("ignore") @@ -402,9 +520,11 @@ def test_d1_divergence_is_pinned(): df, outcome="y", treatment="grp", time="post", fixed_effects=["unit", "time"] ) np.testing.assert_allclose(a.att, f.att, rtol=0, atol=1e-10) - ratio = f.se / a.se - predicted = np.sqrt((360 - 2) / (360 - 66)) - np.testing.assert_allclose(ratio, predicted, rtol=1e-9) + np.testing.assert_allclose(f.se / a.se, 1.0, rtol=1e-9) + # measured pre-fix absorb SE (k=2) rescaled to K_reference = 7: + old_absorb_se = 0.2414226781 + predicted = old_absorb_se * np.sqrt((360 - 2) / (360 - 7)) + np.testing.assert_allclose(a.se, predicted, rtol=1e-8) # --------------------------------------------------------------------------- @@ -451,6 +571,122 @@ def fmt(t): # --------------------------------------------------------------------------- +class TestClusterKIncrement: + """The K_reference helpers: nested-dim classification + the CR1 k + increment (variance-conventions.md D2), on the shared audit panel.""" + + def test_nested_dim_classification(self): + from diff_diff.utils import cluster_nested_fe_dims + + df = make_panel() + u, t = df["unit"].values, df["time"].values + assert cluster_nested_fe_dims(df, ["unit", "time"], u) == ["unit"] + assert cluster_nested_fe_dims(df, ["unit", "time"], t) == ["time"] + # coarser cluster: unit is a function of grp's preimage -> nested + assert cluster_nested_fe_dims(df, ["unit", "time"], df["grp"].values) == ["unit"] + # constructed crossing cluster: varies WITHIN units and periods + assert cluster_nested_fe_dims(df, ["unit", "time"], (u + t) % 2) == [] + assert cluster_nested_fe_dims(df, [], u) == [] + + def test_increment_arithmetic_checks(self): + """The five checks that retrodict the externally-verified arms.""" + from diff_diff.utils import absorbed_fe_cr1_k_increment as incr + + df = make_panel() + u, t = df["unit"].values, df["time"].values + # audit panel absorb form: 0 + 65 - 60 = 5 -> K = 2 + 5 = 7 + assert incr(df, ["unit", "time"], u, has_intercept_col=True) == 5 + # Wooldridge within (no intercept col): 1 + 65 - 60 = 6 = T + assert incr(df, ["unit", "time"], u, has_intercept_col=False) == 6 + # cluster=year inverts nesting: 0 + 65 - 6 = 59 -> K = 61 = 66 - 5 + assert incr(df, ["unit", "time"], t, has_intercept_col=True) == 59 + # ZERO nested dims (the max(...,1) floor case): 0 + 65 - 1 = 64 + assert incr(df, ["unit", "time"], (u + t) % 2, has_intercept_col=True) == 64 + # no absorbed FE -> 0 regardless of intercept flag (the +1 is the + # ABSORBED constant's rank, absent without absorbed dims) + assert incr(df, [], u, has_intercept_col=False) == 0 + + def test_empty_effective_support_returns_zero(self): + """All rows zero-weight: increment 0, never -1 via the floor.""" + from diff_diff.utils import absorbed_fe_cr1_k_increment as incr + + df = make_panel() + u = df["unit"].values + w = np.zeros(len(df)) + assert incr(df, ["unit", "time"], u, has_intercept_col=True, weights=w) == 0 + assert incr(df, ["unit", "time"], u, has_intercept_col=False, weights=w) == 0 + + def test_positive_weight_rows_only(self): + """Zero-weight padding is inference-invariant (REGISTRY guarantee).""" + from diff_diff.utils import absorbed_fe_cr1_k_increment as incr + + df = make_panel() + u = df["unit"].values + base = incr(df, ["unit", "time"], u, has_intercept_col=True) + padded = pd.concat( + [df, pd.DataFrame({"unit": [999], "time": [0], "grp": [0]})], + ignore_index=True, + ) + w = np.r_[np.ones(len(df)), 0.0] + cl = np.r_[u, 999] + assert incr(padded, ["unit", "time"], cl, has_intercept_col=True, weights=w) == base + + def test_nan_keys_raise_actionable(self): + from diff_diff.utils import absorbed_fe_cr1_k_increment as incr + + df = make_panel().astype({"unit": "float64"}) + u = df["unit"].values.copy() + df.loc[0, "unit"] = np.nan + with pytest.raises(ValueError, match="'unit' contains NaN group keys"): + incr(df, ["unit", "time"], u, has_intercept_col=True) + df2 = make_panel() + cl = df2["unit"].values.astype("float64").copy() + cl[0] = np.nan + with pytest.raises(ValueError, match="cluster_ids contain missing values"): + incr(df2, ["unit", "time"], cl, has_intercept_col=True) + + def test_hierarchical_nested_pair_rank(self): + """Two nested dims (unit + state, cluster=state): the nested rank is + the component-aware PAIR rank, not sum(levels-1) — the discriminating + property the full-dummy negative adjustment relies on.""" + from diff_diff.utils import absorbed_fe_rank, cluster_nested_fe_dims + + df = make_panel() + df["state"] = df["unit"] % 6 + state = df["state"].values + nested = cluster_nested_fe_dims(df, ["unit", "state"], state) + assert nested == ["unit", "state"] + pair_rank = absorbed_fe_rank(df, ["unit", "state"], has_intercept_col=True) + naive = (60 - 1) + (6 - 1) + assert pair_rank == 59 # sum(levels) - C - 1 = 66 - 6 - 1, C = n_states + assert pair_rank < naive + + def test_three_dim_increment_inherits_the_documented_d3_approximation(self): + """LIMITATION PIN (documented deviation, not an endorsement): for + THREE or more absorbed dims the increment inherits absorbed_fe_rank's + D3 ``sum(levels) - N + 1`` approximation, which OVER-counts on + duplicated/nested triples. On a(5) x b(4) with c == b, cluster = a: + the true dummy-space rank is 5 + 4 - 1 = 8 (c adds nothing), so the + true increment is 8 - 5 = 3; the formula returns + (5 + 4 + 4) - 3 + 1 = 11 -> increment 6. Pinned so the day the exact + N-way rank ships (TODO.md N-way row; REGISTRY absorbed-FE note) this + test fails loudly and flips to the exact value. One- and two-dim + increments are exact (the tests above).""" + from diff_diff.utils import absorbed_fe_cr1_k_increment + + n = np.arange(80) + df = pd.DataFrame({"a": n % 5, "b": (n // 5) % 4}) + df["c"] = df["b"] # duplicated dim: adds NOTHING to the span + inc3 = absorbed_fe_cr1_k_increment( + df, ["a", "b", "c"], df["a"].to_numpy(), has_intercept_col=True + ) + inc2_true = absorbed_fe_cr1_k_increment( + df, ["a", "b"], df["a"].to_numpy(), has_intercept_col=True + ) + assert inc2_true == 3 # exact two-dim: (5 + 4 - 1) - 5 + assert inc3 == 6 # the documented D3 over-count (true value: 3) + + class TestAbsorbedFeRank: """Component-aware absorbed-FE rank: bit-identity where C=1, correctness where the old ``sum(levels - 1)`` count over-stated rank.""" @@ -810,3 +1046,357 @@ def test_rank_matches_fixest_exact_k(self): # and the naive count reproduces the R DEFAULT naive = (d["state"].nunique() - 1) + (d["state_year"].nunique() - 1) assert naive + 2 == g["iid_default"]["df_k"] + + +# --------------------------------------------------------------------------- +# R + Stata parity: the clustered CR1 non-nested RANK term (K_reference D2). +# --------------------------------------------------------------------------- + +_CR1_NONNESTED_GOLDEN = ( + __import__("pathlib").Path(__file__).parent.parent + / "benchmarks" + / "data" + / "fixest_cr1_nonnested_golden.json" +) + +_REGHDFE_KREF_GOLDEN = ( + __import__("pathlib").Path(__file__).parent.parent + / "benchmarks" + / "data" + / "reghdfe_kref_golden.json" +) + + +def _reghdfe_kref_frame() -> pd.DataFrame: + """The generator's deterministic DGP, rebuilt from the same integer + formulas (no RNG, no embedded data): a disconnected two-way panel -- + units 0-9 observed in periods 0-4, units 10-19 in periods 5-9 (C=2).""" + n = np.arange(1, 101) + unit = (n - 1) // 5 + time = (n - 1) % 5 + 5 * (unit >= 10) + x = ((n * 7) % 13 - 6) / 13 + z = ((n * 11) % 17 - 8) / 17 + out = 0.5 * x + 0.2 * (unit % 3) + 0.1 * time + z + return pd.DataFrame({"unit": unit, "time": time, "x": x, "out": out, "c5": (unit + time) % 5}) + + +def _demeaned_cluster_fit(d: pd.DataFrame, cluster_col: str): + """Absorb [unit, time], fit x on the demeaned design with the wired + K_reference adjustment; returns (coef, se, increment).""" + from diff_diff.linalg import LinearRegression + from diff_diff.utils import absorbed_fe_cr1_k_increment, demean_by_groups + + inc = absorbed_fe_cr1_k_increment( + d, ["unit", "time"], d[cluster_col].to_numpy(), has_intercept_col=False + ) + dm, _ = demean_by_groups(d.copy(), ["out", "x"], ["unit", "time"]) + reg = LinearRegression(include_intercept=False, cluster_ids=d[cluster_col].to_numpy()).fit( + dm[["x"]].values, dm["out"].values, cluster_k_adjustment=inc + ) + return float(reg.coefficients_[0]), float(np.sqrt(reg.vcov_[0, 0])), inc + + +@pytest.mark.skipif( + not _CR1_NONNESTED_GOLDEN.exists(), + reason=( + "fixest_cr1_nonnested_golden.json not present; regenerate via " + "`Rscript benchmarks/R/generate_fixest_cr1_nonnested_golden.R`." + ), +) +class TestFixestCr1NonNestedParity: + """The non-nested absorbed-FE RANK term in the clustered CR1 k, anchored + against fixest on a DISCONNECTED two-way panel (C=2, span rank + U + T - C = 28). + + ``crossed_cluster`` (nothing nested) is the parity arm: the library's + K = x + exact span = 29 matches ``ssc(K.fixef="full", K.exact=TRUE)`` at + machine precision and DIFFERS from both the default approximate count + (df.K=30) and the K.exact-under-nested-K.fixef quirk (df.K=28 -- the + clustered nested path removes 1 df even with nothing nested). It also + exercises the zero-nested-dim ``max(..., 1)`` floor end-to-end. + + ``unit_cluster`` is the documentation arm: NO fixest ssc implements + nested-drop + exact-remainder (library K = 10); the default lands on 11 + (approximate remainder T-1) and the deviation is pinned as an exact + one-df SE ratio.""" + + @staticmethod + def _load(): + import json + + with open(_CR1_NONNESTED_GOLDEN) as fh: + return json.load(fh) + + def test_crossed_cluster_matches_k_exact_full(self): + g = self._load() + d = pd.DataFrame(g["disconnected"]["data"]) + arm = g["disconnected"]["crossed_cluster"] + coef, se, inc = _demeaned_cluster_fit(d, "c5") + assert inc == 28 # exact span U + T - C, nothing nested (the floor case) + assert 1 + inc == arm["k_exact"]["df_k"] == 29 + np.testing.assert_allclose(coef, arm["coef"], rtol=0, atol=1e-12) + np.testing.assert_allclose(se, arm["k_exact"]["se"], rtol=0, atol=1e-12) + # Discriminating: the default approximate count (30) and the + # nested-K.fixef quirk (28) must NOT match. + assert arm["default"]["df_k"] == 30 + assert arm["k_exact_nested"]["df_k"] == 28 + n = g["disconnected"]["n_obs"] + for other in ("default", "k_exact_nested"): + assert abs(se - arm[other]["se"]) > 1e-4 + np.testing.assert_allclose( + arm["default"]["se"] / se, np.sqrt((n - 29) / (n - 30)), rtol=1e-12 + ) + + def test_unit_cluster_deviation_is_exactly_one_df(self): + """The nested COMPOSITION arm: library K=10 (exact remainder 8), + fixest default df.K=11 (approximate remainder T-1=9). Pinned as the + exact one-df ratio so any drift in either convention fails loudly.""" + g = self._load() + d = pd.DataFrame(g["disconnected"]["data"]) + arm = g["disconnected"]["unit_cluster"] + coef, se, inc = _demeaned_cluster_fit(d, "unit") + assert inc == 9 # 1 (constant) + 28 (span) - 20 (nested unit rank) + assert arm["default"]["df_k"] == 11 + assert arm["nested_k_exact"]["df_k"] == 28 # the incoherent composition + np.testing.assert_allclose(coef, arm["coef"], rtol=0, atol=1e-12) + n = g["disconnected"]["n_obs"] + np.testing.assert_allclose( + arm["default"]["se"] / se, np.sqrt((n - 10) / (n - 11)), rtol=1e-12 + ) + + def test_connected_control_matches_default_ssc(self): + """At C=1 with a nested unit cluster the conventions coincide + (remainder T-1 exactly): parity against plain default ssc.""" + g = self._load() + d = pd.DataFrame(g["connected"]["data"]) + coef, se, inc = _demeaned_cluster_fit(d, "unit") + assert 1 + inc == g["connected"]["cluster_default"]["df_k"] == 11 + np.testing.assert_allclose(coef, g["connected"]["coef"], rtol=0, atol=1e-12) + np.testing.assert_allclose(se, g["connected"]["cluster_default"]["se"], rtol=0, atol=1e-12) + + +@pytest.mark.skipif( + not _REGHDFE_KREF_GOLDEN.exists(), + reason=( + "reghdfe_kref_golden.json not present; regenerate via " + "`stata-se -b do benchmarks/stata/generate_reghdfe_kref_golden.do`." + ), +) +class TestReghdfeKReferenceParity: + """The same disconnected design against Stata reghdfe (dof pairwise), the + second reference. ``cross_cluster`` agrees with the library (and with + fixest full+K.exact) at machine precision; ``unit_cluster`` documents + that reghdfe TOO approximates the nested remainder (df_a = 9 = T-1, + implied K = 11), so the library's exact-remainder K = 10 matches no + external reference on disconnected nested designs -- pinned as the same + exact one-df ratio. On CONNECTED designs the compositions coincide and + the jwdid subsample ladder pins machine-precision agreement at every + cluster count.""" + + @staticmethod + def _load(): + import json + + with open(_REGHDFE_KREF_GOLDEN) as fh: + return json.load(fh) + + def test_cross_cluster_exact_span_parity(self): + g = self._load()["cross_cluster"] + d = _reghdfe_kref_frame() + coef, se, inc = _demeaned_cluster_fit(d, "c5") + # reghdfe's df_a IS the exact span here (constants included); + # its denominator N - rank - df_a == N - (1 + inc) == N - 29. + assert g["df_a"] == inc == 28 + assert g["df_a_nested"] == 0 + assert g["df_a_redundant"] == 2 # the C=2 pairwise correction + assert g["df_r"] == g["G"] - 1 + np.testing.assert_allclose(coef, g["coef"], rtol=0, atol=1e-12) + np.testing.assert_allclose(se, g["se"], rtol=0, atol=1e-12) + + def test_unit_cluster_composition_deviation_is_exactly_one_df(self): + g = self._load()["unit_cluster"] + d = _reghdfe_kref_frame() + coef, se, inc = _demeaned_cluster_fit(d, "unit") + assert inc == 9 # library: 1 + 28 - 20 (exact remainder) + # reghdfe: nested drop (20) but APPROXIMATE remainder -- redundant + # counts only the global constant, so df_a = T - 1 = 9 and the + # implied K = rank + df_a + constant = 11 vs the library's 10. + assert g["df_a"] == 9 + assert g["df_a_nested"] == 20 + assert g["df_a_redundant"] == 21 + assert g["df_r"] == g["G"] - 1 + np.testing.assert_allclose(coef, g["coef"], rtol=0, atol=1e-12) + n = g["n"] + np.testing.assert_allclose(g["se"] / se, np.sqrt((n - 10) / (n - 11)), rtol=1e-12) + + +class TestKReferenceConvergence: + """End-to-end properties of the D1/D2 K_reference convergence that the + matrix rows cannot express.""" + + def test_two_nested_dim_full_dummy_uses_rank_not_column_count(self): + """The ONLY design that discriminates rank from raw dummy count. + + With fixed_effects=["unit","state"] and cluster="state" (state = + unit % 6), BOTH dims are cluster-nested and hierarchical: the pair + rank incl. the constant is 66 - 6 - 1 = 59, while the raw drop-first + dummy count is 59 + 5 = 64. The state dummies are collinear with the + unit dummies, so the kernel keeps 62 of 66 columns — and 62 - 59 = 3 + equals the absorb form's K exactly (the collinearity-drop equivalence + two reviewers proved algebraically: k is the design RANK, so + K_ref(full) == K_ref(absorb) identically). A raw-count adjustment + (-64) would give k_inf = -2 and a NaN vcov instead. + """ + df = make_panel() + df = df.assign(state=df["unit"] % 6) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + a = diff_diff.DifferenceInDifferences(cluster="state").fit( + df, outcome="y", treatment="grp", time="post", absorb=["unit", "state"] + ) + f = diff_diff.DifferenceInDifferences(cluster="state").fit( + df, outcome="y", treatment="grp", time="post", fixed_effects=["unit", "state"] + ) + assert np.isfinite(a.se) and a.se > 0 + np.testing.assert_allclose(f.se, a.se, rtol=1e-12) + + @staticmethod + def _saturation_panel(n_units, n_periods, drop=0): + """Deterministic (RNG-free) absorbed-fit fixture; the mod-arithmetic + noise term is outside the unit+time+interaction span on the larger + grids, so above-boundary fits carry genuine residuals.""" + rows = [] + for u in range(n_units): + for t in range(n_periods): + rows.append( + dict( + unit=u, + time=t, + grp=u % 2, + post=int(t >= n_periods - 1), + y=0.5 * u + 0.25 * t + 0.3 * ((u * 13 + t * 7 + (u * t) % 11) % 17), + ) + ) + df = pd.DataFrame(rows) + return df.iloc[: len(df) - drop].reset_index(drop=True) if drop else df + + def _saturation_fit(self, df): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return diff_diff.DifferenceInDifferences(cluster="unit").fit( + df, outcome="y", treatment="grp", time="post", absorb=["unit", "time"] + ) + + def test_clustered_fail_closed_below_the_boundary(self): + """Clustered absorbed fits whose effective dof are exhausted fail + closed to the FULL all-NaN inference tuple — never a silent ~0 SE. + + Both fixtures sit at or below the composed boundary (n = 4 and n = 5 + against K_reference = 4 with absorbed rank 4 on the U=2, T=3 grid; + verified by construction via absorbed_fe_cr1_k_increment == 2 with + kernel k == 2). The exact two-sided n - k_inf boundary in isolation + is pinned at the seam (test_linalg + TestClusterKAdjustmentSeam::test_fail_closed_k_inf_exact_boundary), + where the tail-df lane cannot fire first; here the assertion is the + estimator-level composition: ALL-NaN, no mixed tuple, no finite ~0. + """ + from diff_diff.utils import absorbed_fe_cr1_k_increment + + for drop in (2, 1): + df = self._saturation_panel(2, 3, drop=drop) + inc = absorbed_fe_cr1_k_increment( + df, ["unit", "time"], df["unit"].to_numpy(), has_intercept_col=True + ) + assert inc == 2 # K_reference = kernel k (2) + 2 = 4 >= n - 1 + r = self._saturation_fit(df) + vals = (r.se, r.t_stat, r.p_value, *r.conf_int) + assert all(np.isnan(v) for v in vals), f"n={len(df)}: expected all-NaN, got {vals}" + + def test_clustered_finite_above_the_boundary(self): + """Clearly above the boundary (U=4, T=4: n=16, K_reference = 2 + 3, + absorbed rank 7 -> 7 genuine residual dof) the same configuration + returns a complete finite tuple with a NON-degenerate SE — the + recovery side, guarding against an over-eager fail-closed.""" + r = self._saturation_fit(self._saturation_panel(4, 4)) + vals = (r.se, r.t_stat, r.p_value, *r.conf_int) + assert all(np.isfinite(v) for v in vals) + assert r.se > 1e-8, f"degenerate ~0 SE on a non-saturated fit: {r.se}" + + def test_wcb_identity_and_p_invariance(self): + """WCB: se == sqrt(vcov[att, att]) exactly on adjusted fits across + the absorb, fixed_effects, and TWFE hc2 full-dummy lanes; p-values + are invariant to the corr change (the factor cancels in |t*| vs |t0|) + and CI endpoints move only within the bisection tolerance.""" + df = make_panel() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + r_ab = diff_diff.DifferenceInDifferences( + cluster="unit", inference="wild_bootstrap", n_bootstrap=99, seed=42 + ).fit(df, outcome="y", treatment="grp", time="post", absorb=["unit", "time"]) + r_fe = diff_diff.DifferenceInDifferences( + cluster="unit", inference="wild_bootstrap", n_bootstrap=99, seed=42 + ).fit(df, outcome="y", treatment="grp", time="post", fixed_effects=["unit", "time"]) + r_tw = diff_diff.TwoWayFixedEffects( + vcov_type="hc2", + cluster="unit", + inference="wild_bootstrap", + n_bootstrap=99, + seed=42, + ).fit(df, outcome="y", treatment="grp", time="post", unit="unit") + assert r_ab.se == np.sqrt(r_ab.vcov[3, 3]) + assert r_fe.se == np.sqrt(r_fe.vcov[3, 3]) + assert r_tw.se == np.sqrt(r_tw.vcov[1, 1]) + # the two idioms agree on the corrected SE exactly (cross-design + # bootstrap p-values are NOT compared: draws tie at |t0| within + # float noise, so p can differ by one bootstrap step legitimately) + np.testing.assert_allclose(r_ab.se, r_fe.se, rtol=1e-12) + + def test_wcb_p_invariant_to_the_adjustment(self): + """The corr constant cancels in |t*| vs |t0|: for a FIXED design and + seed, the bootstrap p-value is IDENTICAL under any adjustment, the + reported SE scales by exactly the K_reference factor, and CI + endpoints agree to the bisection tolerance (never exactly).""" + from diff_diff.utils import wild_bootstrap_se + + rng = np.random.default_rng(9) + n, G = 120, 12 + X = np.column_stack([np.ones(n), rng.normal(size=n), rng.normal(size=n)]) + y = X @ np.array([0.5, 0.3, 0.0]) + rng.normal(size=n) + cl = np.repeat(np.arange(G), n // G) + r0 = wild_bootstrap_se(X, y, np.zeros(n), cl, 1, n_bootstrap=199, seed=7) + r5 = wild_bootstrap_se( + X, y, np.zeros(n), cl, 1, n_bootstrap=199, seed=7, cluster_k_adjustment=5 + ) + assert r5.p_value == r0.p_value + expect = np.sqrt((n - 3) / (n - 3 - 5)) + np.testing.assert_allclose(r5.se / r0.se, expect, rtol=0, atol=1e-12) + np.testing.assert_allclose( + r5.t_stat_original * r5.se, r0.t_stat_original * r0.se, rtol=1e-12 + ) + np.testing.assert_allclose( + (r5.ci_lower, r5.ci_upper), (r0.ci_lower, r0.ci_upper), rtol=1e-6 + ) + + def test_lpdid_never_swallows_the_contract_raise(self, monkeypatch): + """LPDiD's broad except around the clustered solve must RE-RAISE + InvalidClusterKAdjustment rather than degrade to a silent unclustered + se=NaN refit (the no-silent-failure rule).""" + import diff_diff.lpdid as lp + from diff_diff.linalg import InvalidClusterKAdjustment + + df = make_panel() + orig = lp.solve_ols + + def inject_bad_adjustment(*a, **k): + if k.get("cluster_ids") is not None: + k["cluster_k_adjustment"] = 10**6 # forces n - k_inf <= 0... valid int + k["vcov_type"] = "hc2_bm" # nonzero adj + non-hc1 -> front door raises + return orig(*a, **k) + + monkeypatch.setattr(lp, "solve_ols", inject_bad_adjustment) + with pytest.raises(InvalidClusterKAdjustment): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + diff_diff.LPDiD(pre_window=2, post_window=2).fit( + df, outcome="y", unit="unit", time="time", treatment="treated" + ) diff --git a/tests/test_wooldridge.py b/tests/test_wooldridge.py index 1826c77a8..573db3c3b 100644 --- a/tests/test_wooldridge.py +++ b/tests/test_wooldridge.py @@ -1886,20 +1886,20 @@ def test_hc1_se_bit_equal_to_pre_pr_baseline(self): """HC1 within-transform path must match the frozen baseline at atol=1e-14. Baseline originally captured on the Phase 1b PR 3/8 branch with - ``_make_vcov_panel(seed=202605211230)`` (FWL preserves the CR1 - cluster-robust score, so the ``vcov_type`` branching kept HC1 bit-equal - to the prior hard-coded HC1 behavior). Recaptured under the v3.6.x - factorize-once + bincount demeaner, which moved the values by ~1e-15 - (bincount accumulation vs pandas' compensated grouped mean - see - REGISTRY "Absorbed Fixed Effects"); the 1e-14 lock semantics are - unchanged. + ``_make_vcov_panel(seed=202605211230)``; recaptured under the v3.6.x + factorize-once + bincount demeaner (~1e-15 move), and recaptured + AGAIN under the 3.9 clustered-CR1 K_reference convergence + (variance-conventions.md D2): the CR1 factor now counts the absorbed + FE not nested in the unit cluster, matching Stata reghdfe (jwdid) at + machine precision on the committed arms. The ATT is unchanged; the SE + moved from the pre-fix 0.03149488781317813 by exactly the K change. """ df = _make_vcov_panel() res = WooldridgeDiD(method="ols", vcov_type="hc1").fit( df, outcome="y", unit="unit", time="time", cohort="cohort" ) assert res.overall_att == pytest.approx(0.9178849934516233, abs=1e-14) - assert res.overall_se == pytest.approx(0.03149488781317813, abs=1e-14) + assert res.overall_se == pytest.approx(0.031906603167527435, abs=1e-14) def test_hc2_bm_finite_and_inflates_over_hc1(self): df = _make_vcov_panel() From b5817fe1b4e9a20fd5da4c0d3dff3a32948bbd56 Mon Sep 17 00:00:00 2001 From: igerber Date: Fri, 31 Jul 2026 05:48:43 -0400 Subject: [PATCH 2/3] docs: refresh stale k_visible comments to the K_reference convention (review round 1, 2x P2) - The _absorbed_fe_vcov_scale gate comments in linalg.py and the MPD twin in estimators.py still said "clustered SEs keep k_visible (fixest ssc nested-FE convention)"; clustered CR1 now carries the K_reference accounting through cluster_k_adjustment inside the kernel, and the two paths are mutually exclusive on cluster_ids - the comments now say so. - test_variance_conventions.py module docstring updated: cr1_k records K_reference counts (visible + the signed adjustment), and the remaining documented-defect rows are the tail-df conventions (PR C), not the closed anti-conservative k accounting. --- diff_diff/estimators.py | 6 ++++-- diff_diff/linalg.py | 12 ++++++++---- tests/test_variance_conventions.py | 17 ++++++++++------- 3 files changed, 22 insertions(+), 13 deletions(-) diff --git a/diff_diff/estimators.py b/diff_diff/estimators.py index 65c0e9d96..960d04251 100644 --- a/diff_diff/estimators.py +++ b/diff_diff/estimators.py @@ -2272,8 +2272,10 @@ def _refit_mp_absorb(w_r): # by k_visible, but the correct finite-sample count is # K_full = k_effective + n_absorbed_effects (matching `df` just above and # fixest feols(vcov="iid"/"hetero")). Rescale so the SE's k agrees with - # the t-df's. Gated exactly as LinearRegression.fit: clustered SEs keep - # k_visible (fixest ssc nested-FE convention), hc2/hc2_bm use + # the t-df's. Gated exactly as LinearRegression.fit: clustered CR1 + # carries the K_reference accounting through `cluster_k_adjustment` + # inside the kernel instead (never this rescale — the gates are + # mutually exclusive on cluster_ids), hc2/hc2_bm use # leverage/Satterthwaite DOF, survey has its own df. When the full-K # residual dof is non-positive the helper returns NaN and we void the # vcov -> NaN inference (fail-closed, per the non-finite-df contract). diff --git a/diff_diff/linalg.py b/diff_diff/linalg.py index 651d43f74..eb36530c2 100644 --- a/diff_diff/linalg.py +++ b/diff_diff/linalg.py @@ -4926,10 +4926,14 @@ def fit( # K_full = n_params_effective_ + df_adjustment (the t-df below already # uses it). Rescale the NON-CLUSTERED iid/hetero vcov so the SE's k # agrees with the t-df's and with fixest feols(vcov="iid"/"hetero"). - # Clustered SEs keep k_visible (fixest ssc nested-FE convention already - # matches); hc2/hc2_bm use leverage/Satterthwaite DOF; survey has its - # own df; full-dummy fits carry df_adjustment == 0. When the full-K - # residual dof is non-positive the helper returns NaN and we void the + # Clustered CR1 does NOT take this rescale: its factor carries the + # K_reference accounting through `cluster_k_adjustment` inside the + # kernel itself (nested-FE drop + non-nested rank; 3.9 convergence), + # and the two paths cannot stack — this gate requires + # cluster_ids is None while the adjustment requires the opposite. + # hc2/hc2_bm use leverage/Satterthwaite DOF; survey has its own df; + # full-dummy fits carry df_adjustment == 0. When the full-K residual + # dof is non-positive the helper returns NaN and we void the # vcov -> NaN inference (fail-closed, per the non-finite-df contract). if ( df_adjustment > 0 diff --git a/tests/test_variance_conventions.py b/tests/test_variance_conventions.py index 3e5536603..1b1199ab8 100644 --- a/tests/test_variance_conventions.py +++ b/tests/test_variance_conventions.py @@ -3,16 +3,19 @@ This is the CI-enforced inventory behind ``docs/methodology/variance-conventions.md``. Each row pins, for one (estimator, fit configuration) cell: -- ``cr1_k`` — the sorted multiset of visible column counts ``k`` reaching the - shared clustered CR1 denominator (``linalg._compute_robust_vcov_numpy``), or - ``()`` when the surface's contract is that it makes NO shared-CR1 call. +- ``cr1_k`` — the sorted multiset of K_reference counts reaching the shared + clustered CR1 denominator: visible columns + the signed + ``cluster_k_adjustment`` (``linalg._compute_robust_vcov_numpy``; 3.9 D1/D2 + convergence), or ``()`` when the surface's contract is that it makes NO + shared-CR1 call. - ``tail_df`` — the sorted multiset of ``df`` values passed to ``safe_inference`` / ``safe_inference_batch`` (``None`` = normal theory). -The point is visibility, not endorsement: several pinned values are DOCUMENTED -DEFECTS (anti-conservative k accounting, mixed tail-df conventions) scheduled to -change in the 3.9 consolidation program. Every row carries ``status`` and, for -legitimate differences, a ``reason``. When a later PR changes a convention, the +The point is visibility, not endorsement: the clustered-CR1 ``k`` cells now +pin the converged K_reference accounting, while the remaining DOCUMENTED +DEFECT rows are the mixed tail-df conventions scheduled for PR C of the 3.9 +consolidation program. Every row carries ``status`` and, for legitimate +differences, a ``reason``. When a later PR changes a convention, the expected literal changes HERE, in one reviewable table. Instrumentation notes (each guards against a failure mode that produced wrong From dd821ea6c6688bcbf877f1d72ea52fbd1c6ccc69 Mon Sep 17 00:00:00 2001 From: igerber Date: Fri, 31 Jul 2026 07:08:49 -0400 Subject: [PATCH 3/3] fix: py3.9 collection error - ladder fixture as module-level, not staticmethod @pytest.fixture stacked on @staticmethod collects fine on Python >= 3.10 but crashes pytest collection on 3.9 (staticmethod has no __name__ there), which is the library floor and exactly one CI matrix leg - the form passed on every other leg and locally (3.14). The module-level fixture form satisfies both constraints: no PytestRemovedIn10Warning (the reason the staticmethod form was introduced) and 3.9-compatible collection. --- tests/test_etwfe_cs_stata_parity.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/tests/test_etwfe_cs_stata_parity.py b/tests/test_etwfe_cs_stata_parity.py index 3a019594c..c6244b134 100644 --- a/tests/test_etwfe_cs_stata_parity.py +++ b/tests/test_etwfe_cs_stata_parity.py @@ -426,6 +426,15 @@ def test_se_matches_jwdid_at_this_cluster_count(self, alltreated_fit): _LADDER_RUNGS = [5, 10, 20, 40, 80, 200, 500] +@pytest.fixture(scope="module") +def ladder(): + """The golden's ladder block. Module-level: a class-scoped fixture + defined as an instance method is deprecated (PytestRemovedIn10Warning), + and the @staticmethod form breaks pytest COLLECTION on Python 3.9 + (staticmethod has no __name__ there) — the library floor's CI leg.""" + return _golden()["ladder"] + + def _ladder_subsample(df: pd.DataFrame, n_per_cohort: int) -> pd.DataFrame: """The generator's roster rule, verbatim: the first ``n_per_cohort`` units per ``first_treat`` cohort by ascending ``countyreal``.""" @@ -452,13 +461,6 @@ class TestSubsampleLadderVsStataJwdid: the resulting number. """ - @pytest.fixture(scope="class") - @staticmethod - def ladder(): - # staticmethod: a class-scoped fixture defined as an instance method - # is deprecated (PytestRemovedIn10Warning). - return _golden()["ladder"] - def test_rung_set_and_roster_rule_are_pinned(self, ladder): """A changed roster rule or panel would move G/n; pin both.""" assert sorted(int(k) for k in ladder["rungs"]) == _LADDER_RUNGS