diff --git a/.claude/memory.md b/.claude/memory.md index 31bd9b265..557b856ea 100644 --- a/.claude/memory.md +++ b/.claude/memory.md @@ -85,7 +85,8 @@ print(results.summary()) The canonical inference quintet is `att` / `se` / `t_stat` / `p_value` / `conf_int` (never `.pvalue` / `.ci`). Every results class exposes all five, but native STORAGE varies in 3.x (`overall_att` on the staggered family, -`avg_att` on MultiPeriodDiD, with canonical names as properties). `summary()` +`avg_att` on the MultiPeriodDiD class removed at the 4.0 cut, with canonical +names as properties). `summary()` everywhere; `to_dict()` / `to_dataframe()` on most classes. See `.claude/../docs/methodology/REGISTRY.md` per estimator. diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 7ccf78cb5..683229cd5 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -201,10 +201,9 @@ A new parameter is only complete when it is: - applied in the **bootstrap/inference** paths, not just the analytical one - reflected on the result object, so `to_dict()`/`summary()` do not misreport it - propagated to the estimators that inherit it: `TwoWayFixedEffects` defines - no `__init__` of its own, and `MultiPeriodDiD`'s 3.9 deprecation shim - forwards via `super().__init__(*args, **kwargs)` with an import-time - `__signature__` mirror - so both inherit a new - `DifferenceInDifferences` constructor parameter automatically; + no `__init__` of its own, so it inherits a new `DifferenceInDifferences` + constructor parameter automatically (the `MultiPeriodDiD` 3.9 shim that + mirrored the signature was removed at the 4.0 cut); `SyntheticDiD` defines its OWN signature (it forwards only `alpha` to `super().__init__`, pinning `cluster=None`), so a new parent parameter must be explicitly mirrored — or explicitly rejected, like its diff --git a/DEFERRED.md b/DEFERRED.md index 285c92b99..2ea98aee8 100644 --- a/DEFERRED.md +++ b/DEFERRED.md @@ -25,7 +25,6 @@ provenance and AI-review deviation-documentation: a row here (or in | LWDiD cohort-relative sensitivity exclusions: `robustness_pre_periods` / `sensitivity_no_anticipation` currently reject multi-cohort staggered inputs because their exclusion windows are defined relative to the EARLIEST adoption (later cohorts' own pre periods fall inside the global post window and survive every restriction, mislabeling the specification). Supporting staggered inputs needs per-cohort window semantics (exclude the last k periods of each cohort's own `t < g` window before its transformation), which the current row-subset design cannot express — a per-cohort masking derivation + its aggregation contract | `diff_diff/lwdid_sensitivity.py` | #588 | Low | | LWDiD estimation-time anticipation offset (CS-parity gap): no `anticipation=` parameter exists — no-anticipation is hard-assumed and the only surface is the single-cohort `sensitivity_no_anticipation` refit diagnostic. The LW device for anticipation is the anchor-window shift — drop the last k pre-treatment periods from the transformation/anchor window (`Ybar_{i,S0}` with `S0 < S-1`, per the REGISTRY LWDiD Anticipation edge case; NOT a cohort recode g → g − k) — and exposing an `anticipation=` parameter needs that derivation checked against LW 2025 NAC/CNAS before implementation | `diff_diff/lwdid.py` | LWDiD release audit | Low | | HonestDiD rejects LWDiD event studies: `LWDiDResults.aggregate('event_study')` emits the unified EventStudyResults surface but the `honest_did.py` allowlist admits CS/Stacked sources only. Admission needs a methodology decision, not an allowlist edit: LWDiD pre-period cells are transformation-based placebos on the Appendix D control pools, and Rambachan-Roth requires those coefficients plus a joint pre/post vcov to stand as valid pre-trend evidence | `diff_diff/honest_did.py`, `diff_diff/lwdid_results.py` | LWDiD release audit | Low | -| HonestDiD non-chronological declared partitions (native `MultiPeriodDiDResults` route): the Rambachan-Roth restriction matrices are built POSITIONALLY over the concatenated declared pre/post lists assuming one chronological boundary, but the native route accepts non-suffix `post_periods` / non-last-pre references and returns bounds whose restriction system does not match the Registry equations (pre-existing; surfaced by the Phase 3(a) calendar-route review, which fails closed instead). Fix = transform the declared partition into boundary form where a valid mapping exists, else reject on the native route too - needs the restriction-geometry derivation. REGISTRY HonestDiD Note records the limitation. | `diff_diff/honest_did.py` | 3(a) | Medium | | `PlaceboTests` `boundary_gap` — a permutation randomization-inference margin (SE-audit item (b)); NOT computed anywhere in code today, so this is a new feature + result field, not a coverage lock. **User-locked 2026-07-09: defer until a derivation/paper source exists** — do not design or implement from scratch. | `tests/test_methodology_placebo.py`, `diff_diff/diagnostics.py` | SE-audit | Low | | TwoStageDiD honest/pretrends container admission DEFERRED (decision revised from "widen" during the 2(b) PR-3b plan review): analytical fits carry the joint Gardner-GMM event-study covariance (M-092), but the pre-period coefficients are stage-1 residual MEANS — the reference horizon is dropped from the no-intercept Stage-2 design and the zero anchor row appended mechanically — not contrasts against a reference period, while HonestDiD's Δ^RM/Δ^SD arithmetic hard-codes the `delta_0 = 0` normalization into its boundary/bridge constraints. Admission needs either a Stage-2 re-estimation with the reference horizon in the design or a derived residual-to-reference normalization mapping (+ its variance transform). Both consumers' TypeErrors state the deferral; see the REGISTRY TwoStageDiD Note (d). | `diff_diff/honest_did.py`, `diff_diff/pretrends.py`, `diff_diff/two_stage_aggregation.py` | 2(b) PR-3b | Low | | RDDensityTest binomial windows test (R rddensity's `bino=` block: exact binomial tests of side counts in shrinking windows around the cutoff): methodology is sourced from Cattaneo, Frandsen & Titiunik (2015, JCI) and Cattaneo, Frandsen & Vazquez-Bare (2017), NOT the reviewed CJM 2020 paper, so implementing it needs its own paper review(s) first (the 2-PR methodology convention). | `diff_diff/rddensity.py` | rddensity PR-B | Low | @@ -121,7 +120,6 @@ For survey-specific limitations (`NotImplementedError` paths), see the | `SpilloverDiD` estimator-level end-to-end vcov reconstruction tests (`bread @ meat @ bread` against `res.vcov`): requires exposing the estimator's internal `X_2_kept` design arrays; the surface is currently pinned from different angles (uniform-weight bit-identity, drift goldens, manual lincom reconstruction at rtol=1e-6). | `spillover.py`, `tests/test_spillover.py` | follow-up | Low | | `SpilloverDiD` Wave E.3 `finite_mask + design-subset` hygiene not yet adopted for TwoStageDiD's analogous pattern (`two_stage.py:567-601`) — separate parity follow-up noted in the `docs/api/spillover.rst` Restrictions block. | `two_stage.py` | Wave E.3 | Low | | `HeterogeneousAdoptionDiD` `covariates=` (Theorem 6 multivariate-covariate extension) not implemented — `fit(covariates=...)` raises `NotImplementedError` via the shipped future-work trap (locked by the `test_had.py` / `test_methodology_had.py` L73 tests); the deferred work is the Theorem 6 extension itself. | `had.py` | Phase 2a | Low | -| MultiPeriodDiD wild bootstrap not supported (falls back to analytical, n_bootstrap-independent) — user-facing edge-case limitation; the 4.0 removal replaces the fallback with a raise (v4-design §4.1), and the merged TWFE event-study mode already RAISES since 3.9 (Phase 3(a)) — this row now governs only the deprecated class itself. | `estimators.py` (MultiPeriodDiD.fit wild-fallback block) | — | Low | | `predict()` raises `NotImplementedError` — rarely needed; user-facing limitation. | `estimators.py:890-911` | — | Low | ## Version-gated (v4) @@ -147,7 +145,6 @@ decisions (refactor waivers, perf trade-offs, test-infrastructure calls) are rec | **Fixed-B Stata SE-golden comparisons in the LWDiD suite run unmarked in default CI** (no `ci_params.bootstrap()` scaling, no `slow` marker): the SE tolerance derives from BOTH fixed rep counts (ours B=999, Stata's R from the golden meta) so scaling either side would invalidate the committed-golden comparison — the ci_params convention governs convergence-style tests, not fixed-B golden comparisons — and these tests are the PR #588 acceptance bar, which must run in default CI; module-scoped fit memoization bounds the cost (~1-2 s per fit locally) | `tests/test_methodology_lwdid.py` | 2026-08-16 | | **The 4.0 migration guide's code blocks are not snippet-executed.** `tests/test_doc_snippets.py` discovers a hardcoded list of `.rst` files and only `.. code-block:: python` / RST `::` bodies, so `docs/migration-4.0.md` gets no coverage. Deliberate: the guide is a MIXED document - most "after" examples (the renames, `results.att`) run on the current release, but the `field-flip` and `df-convention-flip` examples describe 4.0 behaviour that cannot run until 4.0, so a blanket execution lane would fail by construction. Closing the gap fully means a markdown-fence extractor plus a skip-marker convention for the future-API blocks - a harness change, out of scope for a docs PR. What IS gated: the appendix's ledger parity (`test_migration_guide_*`), which pins the row set and every mechanically checkable cell; and, since the first local review found all three merge examples carrying invalid keywords, `test_migration_guide_examples_bind_to_real_signatures`, which ast-parses the guide's python blocks and asserts every constructor/`fit()` keyword exists on the target signature. That is signature binding, NOT execution - it deliberately skips calls whose owner it cannot resolve (e.g. `results.aggregate(...)`), and it cannot catch a wrong VALUE or a wrong sequence of calls. The hand-written `Fix` prose remains unverifiable by any available means. | `docs/migration-4.0.md`, `tests/test_doc_snippets.py` | Phase 4 / 2026-08-09 | | **MMM launch carousel carries scoped marketing claims, not exporter-contract documentation.** The deck (a LinkedIn marketing artifact, not a docs surface) states capabilities in scoped-but-punchy form: guardrail copy says "the easy mistakes fail loudly" / "you own the design, it owns the math" (the caller-owned estimand/population/window/outcome-scale alignment lives in REGISTRY.md's MMM section and the exporter docstrings, not on slides), and slide 8 shows tutorial 29's worked staggered-boost lift row WITHOUT an on-slide linearity qualifier - the compression's linear-channel scoping is the tutorial's job, which the CTA points to. Absolutes ("no silent mis-calibration", "any estimate exports", "anything with an estimate + SE") are ban-listed by `tests/test_mmm_carousel_claims.py`; local review R1/R2 pressed for compliance-style scoping language on the slides and the marketing-appropriate scoped copy was chosen instead. | `carousel/generate_mmm_carousel.py`, `tests/test_mmm_carousel_claims.py` | mmm-carousel / 2026-08-20 | -| **MultiPeriodDiD deprecation shim loses static constructor-arg checking (3.9 window).** The M-010 shim is `__init__(*args, **kwargs)` + an import-time `__signature__` mirror of DiD's constructor: runtime introspection (get_params/set_params, `inspect.signature`) and eager validation are fully preserved, but static type checkers / IDEs cannot check constructor arguments for the deprecated class until the ledger-tracked M-010 removal. Accepted: the alternative (hand-mirroring ~20 parameters) is a drift magnet on a class with one minor version of remaining life. | `diff_diff/estimators.py` | 3(a) / 2026-08-07 | | **DCDH `sklearn.base.clone` param-identity failure won't-fix.** `ChaisemartinDHaultfoeuille._validate_paths_of_interest` unconditionally canonicalizes `paths_of_interest` into a fresh `List[Tuple[int, ...]]`, so sklearn `clone()`'s post-construction `param1 is param2` identity check fails for configured instances - a pre-existing normalization the BaseEstimator mixin PR documented rather than changed (get_params/set_params signatures are clone-compatible; the dependency-free `cls(**est.get_params())` config-equality contract is the enforced one, `tests/test_base_estimator.py`). Fixing would mean returning the caller's raw object from a validator whose job is canonicalization. | `chaisemartin_dhaultfoeuille.py` | mixin PR / 2026-08-01 | | **scikit-learn is a required dev/CI dependency, optional at runtime; clone-identity tests stay `importorskip`-guarded.** Supersedes the 2026-08-01 decision that kept scikit-learn out of the dev extra: #816 added it (unpinned) to the `dev` extra and to both `rust-test.yml` install steps, so the sklearn `clone()` round-trip tests (`test_base_estimator.py`, had/rdd/cic suites) and the DML cross-fit learner tests always run in CI. The library's runtime dependency posture is unchanged (numpy/pandas/scipy only) and the `importorskip` guards remain, so the suite still passes in an environment without scikit-learn; the always-running contract is still the dependency-free re-instantiation config-equality test. | `pyproject.toml`, `tests/test_base_estimator.py` | #816 | | **Plan-review hash gate threat model: accident prevention, NOT adversarial defense.** The ExitPlanMode content-hash gate (hook + `plan_snapshot.py`) exists to stop accidents — stale approvals, concurrent-worktree cross-talk, plans edited mid-review — all of which it closes by construction (snapshot identity + invocation-unique state tokens, 30+ behavioral tests). It does NOT and cannot defend against a malicious local process: nothing verifies review AUTHORSHIP, and such an actor can simply write a matching review file directly — no userland hook can prevent that short of signed reviews, which is out of scope for a personal workflow aid. Review findings that presuppose a hostile local actor against this gate are waived by this decision (2026-07-20, after 7 local AI-review rounds converged on ever-deeper adversarial-model refinements with no reachable fixpoint). Genuine accident vectors remain in scope and are fixed as found. | `.claude/hooks/check-plan-review.py`, `.claude/scripts/plan_snapshot.py` | 2026-07-20 | diff --git a/METHODOLOGY_REVIEW.md b/METHODOLOGY_REVIEW.md index 25cf891cf..4acd1ade0 100644 --- a/METHODOLOGY_REVIEW.md +++ b/METHODOLOGY_REVIEW.md @@ -37,7 +37,7 @@ The catalog grew incrementally over several quarters, so formats vary across the | Estimator | Module | R / Stata Reference | Status | Last Review | |-----------|--------|---------------------|--------|-------------| | DifferenceInDifferences | `estimators.py` | `fixest::feols()` | **Complete** | 2026-01-24 | -| MultiPeriodDiD | `estimators.py` | `fixest::feols()` | **Complete** | 2026-02-02 | +| TwoWayFixedEffects (event-study mode; the former MultiPeriodDiD class, removed at 4.0) | `estimators.py` (`_fit_event_study_core`) + `twfe.py` | `fixest::feols()` | **Complete** | 2026-02-02 | | TwoWayFixedEffects | `twfe.py` | `fixest::feols()` | **Complete** | 2026-02-08 | ### Staggered Treatment Estimators @@ -155,11 +155,11 @@ documented tolerances; cluster-robust and absorbed-FE behavior verified.) --- -#### MultiPeriodDiD +#### TwoWayFixedEffects event-study mode (the former MultiPeriodDiD class) | Field | Value | |-------|-------| -| Module | `estimators.py` | +| Module | `estimators.py` (`DifferenceInDifferences._fit_event_study_core`, the shared core) + `twfe.py` (the `event_study=True` surface; the `MultiPeriodDiD` class, `MultiPeriodDiDResults`, `PeriodEffect` and the `EventStudy` alias were removed at 4.0, PR 4c-ii) | | Primary Reference | Freyaldenhoven et al. (2021), Wooldridge (2010), Angrist & Pischke (2009) | | R Reference | `fixest::feols()` | | Status | **Complete** | @@ -170,23 +170,22 @@ documented tolerances; cluster-robust and absorbed-FE behavior verified.) - [x] Reference period coefficient is zero (normalized by omission from design matrix) - [x] Default reference period is last pre-period (e=-1 convention, matches fixest/did) - [x] Pre-period coefficients available for parallel trends assessment -- [x] Average ATT computed from post-treatment effects only, with covariance-aware SE -- [x] Returns PeriodEffect objects with confidence intervals for all periods +- [x] Average ATT computed from post-treatment effects only, with covariance-aware SE (removed class only; the surviving surface has no post-period-average scalar - `tests/helpers/event_study_average.py` reproduces the arithmetic) +- [x] Per-period rows with confidence intervals for all non-reference periods on `EventStudyResults` (the removed `PeriodEffect` records survive only as the M-016 FutureWarning view, removed 5.0) - [x] Supports balanced and unbalanced panels - [x] NaN inference: t_stat/p_value/CI use NaN when SE is non-finite or zero -- [x] R-style NA propagation: avg_att is NaN if any post-period effect is unidentified +- [x] R-style NA propagation (removed class): `avg_att` was NaN if any post-period effect was unidentified - [x] Rank-deficient design matrix: warns and sets NaN for dropped coefficients (R-style) - [x] Staggered adoption detection warning (via `unit` parameter) - [x] Treatment reversal detection warning - [x] Time-varying D_it detection warning (advises creating ever-treated indicator) - [x] Single pre-period warning (ATT valid but pre-trends assessment unavailable) -- [x] Post-period reference_period raises ValueError (would bias avg_att) +- [x] Post-period reference_period raises ValueError (would exclude a post-period from estimation; biased the removed `avg_att`) - [x] HonestDiD/PreTrendsPower integration uses interaction sub-VCV (not full regression VCV) - [x] All REGISTRY.md edge cases tested **Test Coverage:** -- 50 tests across `TestMultiPeriodDiD` and `TestMultiPeriodDiDEventStudy` in `tests/test_estimators.py` -- 18 new event-study specification tests added in PR #125 +- The `TestMultiPeriodDiD` / `TestMultiPeriodDiDEventStudy` suites (50 tests) and PR #125's 18 event-study specification tests were deleted with the class at 4.0 (PR 4c-ii); the surviving surface is pinned field for field against the committed pre-refactor oracle `tests/fixtures/v4_twfe_event_study_oracles.json` by `tests/test_twfe_event_study_decoupling.py`, exercised by `TestTwoWayFixedEffectsEventStudy` in `tests/test_estimators.py`, and the removed names' absence is pinned by `tests/test_v4_merge_mpd.py` **Corrections Made:** - **PR #125 (2026-02-02)**: Transformed from post-period-only estimator into full event-study @@ -1274,7 +1273,7 @@ CI and extending covariate-adjusted R parity are tracked follow-ups in `DEFERRED - [x] Non-bootstrap CS adapter consumes full `event_study_vcov` sub-block (not diag) - [x] Non-bootstrap SA adapter consumes full `event_study_vcov` sub-block (W-matrix construction `event_study_vcov = W @ vcov_cohort @ W.T` added to `SunAbrahamResults`) - [x] Bootstrap CS/SA and replicate-weight survey paths fall through to `diag(ses^2)` (analytical VCV cleared to prevent mixing with bootstrap/replicate SE overrides) -- [x] `_get_violation_weights('linear')` honors actual pre-period relative-time labels via `fit()` threading → reported MDV is in Roth's γ units on irregular and anticipation-shifted grids. For `MultiPeriodDiDResults`, supported label types are numeric (`int` / `float` / `np.int64`) and `pandas.Period` / `pandas.Timestamp` / `np.datetime64`; **genuinely non-numeric labels** (string period IDs, unranked categoricals) emit an explicit `UserWarning` and fall through to the legacy count-based normalized direction (MDV is NOT in γ units in that case — re-fit with numeric labels) +- [x] `_get_violation_weights('linear')` honors actual pre-period relative-time labels via `fit()` threading → reported MDV is in Roth's γ units on irregular and anticipation-shifted grids. On the TWFE event-study calendar container route (the removed `MultiPeriodDiDResults`'s successor), supported label types are numeric (`int` / `float` / `np.int64`) and `pandas.Period` / `pandas.Timestamp` / `np.datetime64`; **genuinely non-numeric labels** (string period IDs, unranked categoricals) emit an explicit `UserWarning` and fall through to the legacy count-based normalized direction (MDV is NOT in γ units in that case — re-fit with numeric labels) - [x] `PreTrendsPowerResults` persists fitted `violation_weights` + `pretest_form` + `nis_box_probability`; `power_at(M)` works for all four violation types on fresh fits - [x] Helper API (`compute_pretrends_power`, `compute_mdv`) accepts `violation_weights` and `pretest_form`; closes the PR-A R18 helper/class API gap - [x] Summary, `to_dict`, `to_dataframe` dispatch on `pretest_form` (NIS prints box probability; Wald prints noncentrality) @@ -1410,7 +1409,7 @@ Goldens at `benchmarks/data/r_conleyreg_conley_golden.json`; generator at `bench **Outstanding Concerns:** - **Spillover-conley dependency: resolved.** SpilloverDiD ships Conley + survey via Wave E.1/E.2/E.3 (PR #468 / #474 / #482, stratified-Conley sandwich on PSU totals with within-PSU serial Bartlett HAC for `lag_cutoff > 0`); TwoStageDiD ships Wave E.3 parity (PR #485, `fdf2cebc`). -- **Generic LinearRegression / DiD / MPD / TWFE + Conley + survey_design: deferred.** Two fail-closed contracts assert the unsupported combination: the estimator-level gate (DiD / MPD / TWFE) lives in `diff_diff/conley.py::_validate_conley_estimator_inputs` (`TestConleyDeferrals::test_did_mpd_twfe_survey_design_not_implemented`); the generic LinearRegression gate lives in `diff_diff/linalg.py::LinearRegression.fit` (`TestConleyDeferrals::test_linear_regression_survey_design_not_implemented`). The `survey_design=` surface is `LinearRegression` only; `compute_robust_vcov` does not accept `survey_design=`. +- **Generic LinearRegression / DiD / TWFE (both modes) + Conley + survey_design: deferred.** Two fail-closed contracts assert the unsupported combination: the estimator-level gate (DiD / TWFE, both modes) lives in `diff_diff/conley.py::_validate_conley_estimator_inputs` (`TestConleyDeferrals::test_did_mpd_twfe_survey_design_not_implemented`); the generic LinearRegression gate lives in `diff_diff/linalg.py::LinearRegression.fit` (`TestConleyDeferrals::test_linear_regression_survey_design_not_implemented`). The `survey_design=` surface is `LinearRegression` only; `compute_robust_vcov` does not accept `survey_design=`. - **Generic LinearRegression / compute_robust_vcov + Conley + `weights=`: deferred for any `weight_type` (`pweight` / `aweight` / `fweight`).** Weighted Conley is not implemented on the generic linalg surface. The `pweight` / `survey_design` subset additionally reflects an **open methodological question** — no canonical extension of Conley (1999) exists for weighted spatial-HAC under probability sampling. `TestConleyDeferrals::test_conley_plus_weights_not_implemented`. - **SyntheticDiD + Conley**: raises `TypeError`. SyntheticDiD uses bootstrap / jackknife / placebo variance, not the analytical sandwich. `TestConleyDeferrals::test_synthetic_did_conley_typeerror`. - **Wild bootstrap + Conley**: `NotImplementedError`. Wild bootstrap is a separate inference path that does not consume the analytical Conley sandwich (REGISTRY L3547). `TestConleyDeferrals::test_wild_bootstrap_conley_not_implemented`. diff --git a/TODO.md b/TODO.md index bc71854f0..b88812607 100644 --- a/TODO.md +++ b/TODO.md @@ -30,12 +30,12 @@ Related tracking surfaces: | Staggered-DDD power support: `simulate_power`/`simulate_mde`/`simulate_sample_size` now REJECT a staggered-configured `TripleDifference` (both registered DDD generators emit 2x2x2 data and fit with `(group, partition, post)`, so a staggered config would be simulated under the wrong design). Support needs a staggered DDD DGP profile plus fit-kwargs builder, and a decision on whether the mode is selected by profile or by the estimator's own config | `diff_diff/power.py` | 3(b) | Mid | Low | | Bootstrap-`seed` provenance on multiplier-bootstrap results containers: neither `StaggeredTripleDiffResults` nor `CallawaySantAnnaResults` carries the `seed` that generated its bootstrap SEs / p-values / sup-t bands, so a serialized result cannot report the random configuration behind its inference. NOT a 3(b) regression - `seed` reaches the engine and `get_params()` correctly (same seed reproduces the SE bit-exactly, a different seed moves it), the gap is results-object observability only, it predates the merge, and both containers inherit it from the shared `CallawaySantAnnaBootstrapMixin`. Add `seed` (and consider `n_bootstrap`/`bootstrap_weights`/`cband`) to BOTH containers plus `to_dict()`, with seeded and unseeded pins; sequence it with the M-014 container unification rather than schema-changing one container mid-merge. Precedent for exposing it: `ContinuousDiDResults`, `EfficientDiDResults`, `SyntheticDiDResults` already do | `diff_diff/staggered_triple_diff_results.py`, `diff_diff/staggered_results.py` | 3(b) | Quick | Low | | Staggered-mode cluster-robust ANALYTICAL SEs: `cluster=` raises in `TripleDifference`'s staggered mode, so clustered inference there is bootstrap-only. Implementing a clustered analytical path for the GMM-combined influence function would let the raise become a real lane | `diff_diff/_staggered_triple_diff_engine.py` | 3(b) | Heavy | Low | -| diagnostic_report admission for `EventStudyResults` surfaces (the TWFE event-study mode + `aggregate('event_study')` containers): DiagnosticReport/BusinessReport now REJECT the surface explicitly (Phase 3(a); previously a silent zero-check report / all-null headline) and practitioner_next_steps serves the generic fall-through - admission needs source-aware routing (the type-name-keyed `_APPLICABILITY`/`_HANDLERS` registries cannot discriminate the unified container's producers) and a scalar-vs-per-period headline design; MPD-native results received {parallel_trends, pretrends_power, sensitivity, bacon, design_effect} | `diff_diff/diagnostic_report.py`, `diff_diff/business_report.py`, `diff_diff/practitioner.py` | 3(a) | Mid | Medium | +| diagnostic_report admission for `EventStudyResults` surfaces (the TWFE event-study mode + `aggregate('event_study')` containers): DiagnosticReport/BusinessReport now REJECT the surface explicitly (Phase 3(a); previously a silent zero-check report / all-null headline) and practitioner_next_steps serves the generic fall-through - admission needs source-aware routing (the type-name-keyed `_APPLICABILITY`/`_HANDLERS` registries cannot discriminate the unified container's producers) and a scalar-vs-per-period headline design; MPD-native results received {parallel_trends, pretrends_power, sensitivity, bacon, design_effect}. PR 4c-ii (2026-10-25) DELETED the MPD-native branches that logic lived in (DiagnosticReport's `interaction_indices`+`vcov` joint-Wald covariance source, the `post_period_effects` heterogeneity source, the `pre_period_effects` Bonferroni/NaN-filter collector branch, the `avg_att` alias row, practitioner's `_handle_multi_period`) - recover them from git history (`git show 5a10b847^:diff_diff/diagnostic_report.py`) when admission lands; event-study results reach `_handle_generic` meanwhile | `diff_diff/diagnostic_report.py`, `diff_diff/business_report.py`, `diff_diff/practitioner.py` | 3(a) | Mid | Medium | | DiagnosticReport public skip bookkeeping omits RUNNER-level skips: `applicable_checks` reflects only gate outcomes, so a check whose runner returns `status="skipped"` (heterogeneity's empty/non-finite branches, `_pt_event_study`'s empty-coefs branch - a pre-existing pattern, now also reachable via a failed post-fit event-study derivation on bootstrapped / kit-less ImputationDiD/TwoStageDiD/ContinuousDiD fits) stays listed as applicable while `skipped_checks`/`schema["skipped"]` omit it, so automation reading the public fields can misclassify the check as completed; reconcile runner-returned skipped sections into the public bookkeeping (all checks, one convention) or resolve those availabilities at the gate | `diff_diff/diagnostic_report.py` | derived-ES review | Mid | Low | -| Align the NATIVE pretrends pre_periods= contract with the container routes' fail-closed validation: `_extract_pre_period_params`'s MPD branch silently filters an explicit `pre_periods=` selection (unknown labels, the reference, unusable-inference rows dropped without error, caller order preserved) while both container routes validate every requested label and enforce calendar chronology (the relative route since M-024; the calendar route since 3(a) R8) - the relative route also silently collapses DUPLICATE requested labels (the calendar route rejects them), and both filter on SE only while the calendar route additionally requires a finite EFFECT (3(a) R9); one contract across all three routes, with pinned rejection messages | `diff_diff/pretrends.py` | 3(a) R8 | Quick | Low | +| Align the two pretrends container routes' `pre_periods=` contracts (the native MPD branch died in PR 4c-ii, so only the relative and calendar routes remain): the relative route silently collapses DUPLICATE requested labels (the calendar route rejects them, 3(a) R8), and it filters on SE only while the calendar route additionally requires a finite EFFECT (3(a) R9); one contract across both routes, with pinned rejection messages | `diff_diff/pretrends.py` | 3(a) R8 | Quick | Low | | `EventStudyResults` inference-provenance fields: the container records no `vcov_type`/`cluster_name`/`n_clusters`/`df_convention`/Conley metadata, so a serialized surface cannot distinguish unit auto-clustering from explicit clustering, survey, Conley, or the one-way carve-out (3(a) R9 review). Adding them is a cross-producer M-092 schema amendment (six builders, to_dict/summary rendering, surface-suite pins) - follow the pre-cut amendment convention (optional fields appended last, ledger note same-diff) rather than bolting onto one producer | `diff_diff/results_base.py` | 3(a) R9 | Mid | Low | | Opt-in singleton-group pruning for TwoWayFixedEffects (static + event-study mode; reghdfe parity): singleton units/periods are currently RETAINED class-wide - the within-demeaned row is zero so points are unchanged, but N/G/residual-df count it and CR1/finite-sample SEs shift (~0.41019 -> 0.40962 measured; REGISTRY "Deviation from R" Note, R5 review) - reghdfe iteratively drops singletons by default while fixest retains them (diff-diff matches fixest); an opt-in knob needs iterative unit+period pruning with consistent cluster/survey/replicate/Conley array subsetting and a default-flip decision protocol (moves published SEs) | `diff_diff/twfe.py`, `diff_diff/estimators.py`, `diff_diff/utils.py` | 3(a) R5 | Mid | Low | -| Cohort-timing validation input for the simultaneous-adoption event-study family (TWFE `event_study=True` + MultiPeriodDiD through 3.9): an optional `first_treat=`/`cohort=` column so simultaneous adoption becomes checkable under the contract-valid time-invariant `D_i` indicator - today the staggered-adoption advisory derives timing from within-unit 0->1 transitions, so it can only fire on off-contract time-varying `D_it` input, and with valid `D_i` adoption timing is not observable in the inputs at all (REGISTRY "staggered-adoption detection limit" Notes, both sections); design questions: validate-only vs steering error, and interplay with the M-011 removal | `diff_diff/twfe.py`, `diff_diff/estimators.py` | 3(a) R2 | Mid | Medium | +| Cohort-timing validation input for the simultaneous-adoption event-study family (TWFE `event_study=True`; MultiPeriodDiD through 3.9): an optional `first_treat=`/`cohort=` column so simultaneous adoption becomes checkable under the contract-valid time-invariant `D_i` indicator - today the staggered-adoption advisory derives timing from within-unit 0->1 transitions, so it can only fire on off-contract time-varying `D_it` input, and with valid `D_i` adoption timing is not observable in the inputs at all (REGISTRY "staggered-adoption detection limit" Notes, both sections); design questions: validate-only vs steering error (the M-011 removal shipped in PR 4c-ii, so only the TWFE surface is in play) | `diff_diff/twfe.py`, `diff_diff/estimators.py` | 3(a) R2 | Mid | Medium | | CS parity with DMLDiD's replicate conventions: per-cell replicate SEs for `CallawaySantAnna` (CS leaves per-cell SEs on the weighted sqrt-sum under replicate designs; DMLDiD now computes them via `compute_replicate_if_variance`) AND flip the replace-style effective-df relays to min-cap (`min(df_survey, n_valid - 1)`) at all five sites: CS overall (`staggered.py` overall relay), the shared event-study/group aggregation sites (`staggered_aggregation.py` `df_survey_val = min(non_none_dfs)` x2), and the staggered-DDD engine twins (`_staggered_triple_diff_engine.py` overall + event-study) — the replace convention can RAISE df above the QR-rank design df (anti-conservative); latent, fires only when `n_valid < R` (all-zero/non-finite replicate column). DMLDiD's min-cap + the pinned inherited ES/group behavior are the templates (REGISTRY DMLDiD replicate Note) | `diff_diff/staggered.py`, `diff_diff/staggered_aggregation.py`, `diff_diff/_staggered_triple_diff_engine.py` | DML replicate PR | Mid | Low | | `balance_e` replay on bootstrapped ImputationDiD/TwoStageDiD fits: the unrestricted ES/group surface is relayed since 4.0 PR 3c (computed at fit, `balance_e=None`; M-021/M-022), but `aggregate('event_study', balance_e=k)` still raises there - a balanced-window replay needs retained per-target psi (ImputationDiD: the panel-backed kit retains everything the psi precompute reads, so a seeded replay is tractable) / retained per-level GMM scores (TwoStageDiD: function-locals today) | `diff_diff/imputation_results.py`, `diff_diff/two_stage_results.py`, `diff_diff/aggregation.py` | 2(b) PR-3b | Mid | Low | | ContinuousDiD post-fit bootstrap-ES REPLAY: since 4.0 PR 3c a bootstrapped fit relays the event-study surface computed at fit (M-025), so the gap is a genuine post-fit replay (fresh draws, no fit-side cost, a route for pre-4.0 pickles); a seeded replay is tractable - the multiplier draws are seeded (`np.random.default_rng(self.seed)`) - but needs the FULL per-cell `_bootstrap_info` (bread/ee_treated/Psi_eval/dPsi_*/beta_pred) the pruned kit deliberately drops, so shipping it means a kit-payload change with its own memory contract | `diff_diff/continuous_did_aggregation.py`, `diff_diff/continuous_did_results.py` | 2(b) PR-3c | Mid | Low | @@ -94,11 +94,11 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Geo-experiments tutorial: the SDiD fits on the `generate_factor_data(...)` panel (conversions ~1,500, sigma-hat ~310) hit the Frank-Wolfe iteration cap - the placebo fit warns once and the bootstrap cross-check warns on 100/100 draws - yet the narrative presents the bootstrap SE as a paper-faithful robustness check. Pre-existing on `main` under BOTH backends (surfaced when the notebook outputs were refreshed for the shape-only fit diagnostic); the committed April outputs predate the aggregated non-convergence warning. Either rescale / re-tune the DGP (or raise `max_iter` / loosen `min_decrease`) so the refits converge and re-execute, or label the bootstrap result unreliable and drop it from the "Validity evidence" claim. | `docs/tutorials/18_geo_experiments.ipynb` | local codex review of the SDiD shape-only fit-RMSE PR | Mid | Low | | Replicate Chang (2020) §4.2.1 ML-design RCS DGP — requires a penalized (Lasso-class) propensity learner or a maintainer fixture via the duck-typed `propensity_learner=` object route; native logit+linear verified pure noise at (N=500, p=100) across ~20 review seeds (att ~ −23..+16, SE ~4-10, EPV ~1.4-1.6 vs threshold 10, fitted out-of-fold clipping ~5-29%, mean ~14%); at (N=200, p=100) `outcome_learner="linear"` fails closed on control-fold rank deficiency even under an oracle propensity (use ridge/sieve there); §4.2.2 is replicated | `tests/test_methodology_dml_did.py`, `diff_diff/_learners.py` | DML PR-B2 | Mid | Low | | Optional scheduled end-to-end execution gate for the MMM tutorials (29/30): a cron-only workflow (or extension of `mmm-interop.yml`) that executes both notebooks in isolated exact-pin environments, so a stale/invalid committed posterior cannot stay green indefinitely - today the hybrid posture (deliberate: notebooks execute locally with committed outputs; CI smoke-tests the exporters without sampling; drift tests pin source + committed-output needles) leaves the MCMC claims un-re-executed in CI | `.github/workflows/mmm-interop.yml`, `docs/tutorials/29_mmm_calibration_pymc.ipynb`, `docs/tutorials/30_mmm_calibration_meridian.ipynb` | mmm-interop | Mid | Low | -| Committed `fixest::feols` event-study golden for TWFE `event_study=True` (within + pooled specs, unbalanced + covariate panels, matched CR1 cluster convention, per-period effects + vcov block) - the in-suite gates are shared-core cross-checks (TWFE-within == MPD-absorb, pooled == MPD bit-exact), so a defect common to the shared core would pass; the live-R harness (`benchmarks/R/benchmark_multiperiod.R`, `feols(y ~ treated * time_f \| unit)`) validated the within design in `docs/benchmarks.rst` but is not a committed regression test - follow the `fixest_did_twfe_golden.json` committed-golden pattern (pytest.skip when absent) | `tests/test_fixest_did_twfe_parity.py`, `benchmarks/R/` | 3(a) R2 | Mid | Medium | +| Committed `fixest::feols` event-study golden for TWFE `event_study=True` (within + pooled specs, unbalanced + covariate panels, matched CR1 cluster convention, per-period effects + vcov block) - the in-suite gate is the committed pre-refactor capture (`tests/fixtures/v4_twfe_event_study_oracles.json`, a self-capture of the shared core at 68e383bd/d45d5316), so a defect that predates the capture would pass; the live-R harness (`benchmarks/R/benchmark_multiperiod.R`, `feols(y ~ treated * time_f \| unit)`) validated the within design in `docs/benchmarks.rst` but is not a committed regression test - follow the `fixest_did_twfe_golden.json` committed-golden pattern (pytest.skip when absent) | `tests/test_fixest_did_twfe_parity.py`, `benchmarks/R/` | 3(a) R2 | Mid | Medium | | Evaluate adding the `BaseEstimator` param surface (get_params/set_params) to the exported classes that never had it - `PowerAnalysis`, `LinearRegression`, `BusinessReport`, `DiagnosticReport`, `TWFEWeightsResult` (a NEW public surface, deliberately out of the 2(c)-i pure-refactor scope; `LinearRegression` is the one `fit`-bearing class excluded from the contract suite's roster-completeness test). | `diff_diff/linalg.py`, `diff_diff/power.py` | mixin PR | Mid | Low | | Tighten the mypy suppressions that back the enforced-zero posture: burn down `prep_dgp`'s per-module `[index]` override (needs a None-vs-array restructure that preserves the seeded RNG stream), and evaluate re-enabling the globally disabled codes (`arg-type`, `return-value`, `var-annotated`, `assignment`) one at a time — `assignment` alone hid several real annotation drifts found during the 2026-07 triage. | `pyproject.toml` `[tool.mypy]`, `diff_diff/prep_dgp.py` | lint-CI | Mid | Low | | Route LWDiD into the practitioner code surfaces: `practitioner_next_steps` (`diff_diff/practitioner.py`), `agent_workflow.py`, `DiagnosticReport`, and `BusinessReport` never dispatch to or mention LWDiD — the AI guides now cover it, but the code surfaces need contract decisions (when to recommend it, what diagnostics to run) plus tests | `diff_diff/practitioner.py`, `diff_diff/agent_workflow.py`, `diff_diff/diagnostic_report.py`, `diff_diff/business_report.py` | LWDiD release audit | Mid | Low | -| df-provenance completion follow-up (PR C descoping): a Wooldridge public per-row df channel (a per-key dict that can also carry the hc2_bm BM contrast dofs the scalar `_df_analytic_fallback` cannot represent — the unified surface's Wooldridge df column stays NaN until then); ImputationDiD NON-survey per-lead df (the survey-fit scalar `event_study_df` + its bootstrap clearing shipped with the M-092 completion; what remains is surfacing the `resolve_tail_df`-derived per-lead dfs on non-survey fits, which needs the `_compute_lead_coefficients` return-path change); and the `inference_df` (DiD/MPD/SA/StackedDiD) vs `df_inference` (CS/dCDH) adapter-naming unification. NOTE (M-024): StackedDiD's `df_survey=None` on analytical fits is DELIBERATE - the container resolver's `df_inference` fallback not matching `inference_df` is load-bearing (normal-theory FLCI; see the REGISTRY StackedDiD M-024 Note and the pin in `tests/test_event_study_consumers.py`) - the unification must preserve that outcome or renegotiate it explicitly. | `diff_diff/wooldridge_results.py`, `diff_diff/imputation.py`, `diff_diff/results_base.py` | #variance-PR-C | Mid | Low | +| df-provenance completion follow-up (PR C descoping): a Wooldridge public per-row df channel (a per-key dict that can also carry the hc2_bm BM contrast dofs the scalar `_df_analytic_fallback` cannot represent — the unified surface's Wooldridge df column stays NaN until then); ImputationDiD NON-survey per-lead df (the survey-fit scalar `event_study_df` + its bootstrap clearing shipped with the M-092 completion; what remains is surfacing the `resolve_tail_df`-derived per-lead dfs on non-survey fits, which needs the `_compute_lead_coefficients` return-path change); and the `inference_df` (DiD/SA/StackedDiD) vs `df_inference` (CS/dCDH) adapter-naming unification. NOTE (M-024): StackedDiD's `df_survey=None` on analytical fits is DELIBERATE - the container resolver's `df_inference` fallback not matching `inference_df` is load-bearing (normal-theory FLCI; see the REGISTRY StackedDiD M-024 Note and the pin in `tests/test_event_study_consumers.py`) - the unification must preserve that outcome or renegotiate it explicitly. | `diff_diff/wooldridge_results.py`, `diff_diff/imputation.py`, `diff_diff/results_base.py` | #variance-PR-C | Mid | Low | | `worktree-rm` safety via a tested argv helper: the prose rewrite (ask-before-remove confirmation gate, detached-HEAD reachability/rescue, tip-identity force-delete guard) was reverted to the main version because editing the prose repeatedly reintroduced shell-injection (last: sourcing a state file built from a git-derived branch name). Restore those guards in a `worktree_rm.py` that takes the name via file ingress, invokes git through argv arrays, and has a metacharacter-branch/path injection regression test — the pattern that worked for `pr_prepare.py`/`premerge_scan.py`. | `.claude/commands/worktree-rm.md`, `.claude/scripts/` | skill-audit | Heavy | Medium | | `premerge_scan.py` should scan the staged blob (`git show :path`) for staged methodology files, not the working-tree copy — a stage-then-revert-working-copy edit currently reads the safe working version and misses the staged violation. Union staged-index findings with unstaged/untracked filesystem findings. | `.claude/scripts/premerge_scan.py` | skill-audit | Mid | Low | | Carousel PDF drift guard: the committed deck PDFs are build artifacts of their `generate_*.py` scripts (README rule: regenerate in the same commit), but nothing verifies a committed PDF matches its generator — a source edit without regeneration, or a stale raster chart, passes every claims test. Applies to ALL decks, not just the MMM one. fpdf embeds creation-date metadata, so a byte-compare needs fixed metadata or a per-page rendered-text manifest. Origin: mmm-carousel local review R1. | `carousel/`, `tests/` | mmm-carousel | Mid | Low | diff --git a/benchmarks/R/benchmark_multiperiod.R b/benchmarks/R/benchmark_multiperiod.R index 9dbb62a8e..4965ad3b4 100644 --- a/benchmarks/R/benchmark_multiperiod.R +++ b/benchmarks/R/benchmark_multiperiod.R @@ -1,5 +1,6 @@ #!/usr/bin/env Rscript -# Benchmark: MultiPeriodDiD event study (R `fixest` package) +# Benchmark: the TwoWayFixedEffects event-study arm (within spec; the 3.x +# MultiPeriodDiD arm) - R `fixest` package side # # Usage: # Rscript benchmark_multiperiod.R --data path/to/data.csv --output path/to/results.json \ @@ -98,7 +99,7 @@ if (config$warmup) { } # Run benchmark -message("Running MultiPeriodDiD estimation (fixest::feols)...") +message("Running the TwoWayFixedEffects event-study arm (within spec; the 3.x MultiPeriodDiD arm) estimation (fixest::feols)...") start_time <- Sys.time() # Regression: outcome ~ treated * time_f | unit, clustered SEs diff --git a/benchmarks/R/generate_clubsandwich_golden.R b/benchmarks/R/generate_clubsandwich_golden.R index 7c466c360..045be8a21 100644 --- a/benchmarks/R/generate_clubsandwich_golden.R +++ b/benchmarks/R/generate_clubsandwich_golden.R @@ -123,9 +123,12 @@ output$absorbed_fe_did <- list( dof_cr2 = as.numeric(ct_did_cr2$df_Satt) ) -# --- Absorbed-FE MultiPeriodDiD event-study scenario (gate lift PR) ---------- +# --- Absorbed-FE event-study scenario (gate lift PR) ------------------------- +# (Originally the 3.x MultiPeriodDiD absorb= design, removed in the 4.0 cut; the +# `mpd_*` golden keys are kept verbatim - tests/test_linalg_hc2_bm.py consumes +# them as linalg-helper parity inputs.) # Mirrors MPD(fixed_effects=["unit"]) destination of the absorb auto-route on -# MultiPeriodDiD. MPD parameterization: const + treated + period_f (non-ref) +# the former MultiPeriodDiD. MPD parameterization: const + treated + period_f (non-ref) # + treated:period_X (non-ref) + factor(unit). Build the interaction columns # explicitly so the R fit's coefficient names match MPD's `treated:period_X`. @@ -175,7 +178,7 @@ output$mpd_absorbed_fe_did <- list( # --- MPD clustered avg_att DOF scenario (Gate 6 lift PR) --------------------- # Pins clubSandwich's compound-contrast Satterthwaite DOF for the post-period- -# average ATT under cluster-robust CR2. Mirrors MultiPeriodDiD(cluster=unit, +# average ATT under cluster-robust CR2. Mirrors the former MultiPeriodDiD(cluster=unit, # vcov_type='hc2_bm', fixed_effects=['unit']) parameterization. Per-coefficient # DOFs use coef_test()$df_Satt (the canonical Satterthwaite per-coef API); # the compound contrast DOF uses Wald_test(constraints=matrix(c_avg, 1), @@ -642,7 +645,7 @@ output$weighted_did_absorbed_fe <- list( # "HTZ")$df_denom` against `_compute_cr2_bm_contrast_dof(weights=)` on an # MPD-shaped design. NOTE: like `weighted_did_absorbed_fe`, this fixture # targets the analytical surface only — estimator-level -# `MultiPeriodDiD(survey_design=...)` paths use TSL. +# the former `MultiPeriodDiD(survey_design=...)` paths used TSL. set.seed(50600) d_mpd_w <- make_mpd_panel(n_total = 15, units_per_cohort = 5, n_periods = 4, seed = 50601) diff --git a/benchmarks/python/benchmark_honest.py b/benchmarks/python/benchmark_honest.py index 8a5983889..aeeda8635 100644 --- a/benchmarks/python/benchmark_honest.py +++ b/benchmarks/python/benchmark_honest.py @@ -63,15 +63,16 @@ def __init__(self, beta, sigma, num_pre, num_post): self.vcov = sigma self.n_pre_periods = num_pre self.n_post_periods = num_post - # Create period effects for compatibility + # Plain per-period records for compatibility (not an estimator + # container; HonestDiD reads beta/vcov above). self.period_effects = [] for i, b in enumerate(beta): event_time = i - num_pre + 1 if i < num_pre else i - num_pre + 1 - self.period_effects.append(type('PeriodEffect', (), { + self.period_effects.append({ 'period': event_time, 'coefficient': b, - 'se': np.sqrt(sigma[i, i]) if i < len(sigma) else 0.1 - })()) + 'se': np.sqrt(sigma[i, i]) if i < len(sigma) else 0.1, + }) mock_results = MockEventStudyResults(beta, sigma, args.num_pre, args.num_post) diff --git a/benchmarks/python/benchmark_multiperiod.py b/benchmarks/python/benchmark_multiperiod.py index 6080f1a1e..4035d582b 100644 --- a/benchmarks/python/benchmark_multiperiod.py +++ b/benchmarks/python/benchmark_multiperiod.py @@ -2,6 +2,16 @@ """ Benchmark: MultiPeriodDiD event study (diff-diff MultiPeriodDiD). +FROZEN 3.x ARM - NOT RUNNABLE ON A 4.0 TREE. ``MultiPeriodDiD`` was removed in +the 4.0 cut (PR 4c-ii); this script is kept verbatim because the +``benchmarks/refresh_2026_07/`` harness (``run_refresh.py`` HEADLINE_PIN 3.7.0, +``run_version_story.py`` 3.5.3/3.7.0, via ``refresh_common.run_python_rep``) +executes it inside PINNED diff-diff venvs whose wheels predate +``event_study=True``, and the committed refresh results carry its labels. +The live arm consumed by ``run_benchmarks.py`` is +``benchmarks/python/benchmark_twfe_event_study.py`` (same CLI and output keys; +``TwoWayFixedEffects(...).fit(..., event_study=True, spec="within")``). + Usage: python benchmark_multiperiod.py --data path/to/data.csv --output path/to/results.json \ --n-pre 4 --n-post 4 diff --git a/benchmarks/python/benchmark_twfe_event_study.py b/benchmarks/python/benchmark_twfe_event_study.py new file mode 100644 index 000000000..c4f45434b --- /dev/null +++ b/benchmarks/python/benchmark_twfe_event_study.py @@ -0,0 +1,284 @@ +#!/usr/bin/env python3 +""" +Benchmark: TwoWayFixedEffects event study (``event_study=True, spec="within"``). + +Live successor of the frozen 3.x ``benchmark_multiperiod.py`` arm: the same +unit-FE event-study design (``outcome ~ treated x period | unit``, clustered +by unit), fit through ``TwoWayFixedEffects(cluster=...).fit(..., +event_study=True, spec="within")``, which returns the unified +``EventStudyResults`` surface instead of the removed ``MultiPeriodDiDResults`` +container. The CLI and the output JSON keys are identical to the frozen arm so +``run_benchmarks.py`` and ``compare_results.py`` consume it unchanged, and the +R side stays ``benchmarks/R/benchmark_multiperiod.R`` (the model is the same). + +Two conversions the frozen arm also carried: + +* The surface's ``event_time`` holds CALENDAR period labels; the R script and + ``compare_results.compare_event_study`` join on reference-RELATIVE + ``event_time``, so each row is emitted as ``{"period": label, "event_time": + label - reference_period, ...}``. +* The surface has no post-period average. The headline ``att``/``se`` are the + former post-average arithmetic, inlined below exactly as + ``tests/helpers/event_study_average.py`` (benchmarks cannot import + ``tests/``): equal-weight mean of the post-period ``att`` rows and + ``sqrt(sum(post vcov block) / n_post**2)`` with the Cartesian ``np.ix_`` + block (plain ``vcov[cols, cols]`` would keep only the diagonal and drop the + cross-period covariances). Student-t p-value / 95% CI use the surface's + shared post-row ``df``; R uses ``nobs - k`` without the absorbed-FE + subtraction - a pre-existing difference, and ``compare_estimates`` gates + ``att``/``se`` only. + +Usage: + python benchmark_twfe_event_study.py --data path/to/data.csv \ + --output path/to/results.json --n-pre 4 --n-post 4 +""" + +import argparse +import json +import os +import sys +from pathlib import Path +from typing import Any, Dict + + +# IMPORTANT: Parse --backend and set environment variable BEFORE importing diff_diff +# This ensures the backend configuration is respected by all modules +def _get_backend_from_args(): + """Parse --backend argument without importing diff_diff.""" + parser = argparse.ArgumentParser(add_help=False) + parser.add_argument("--backend", default="auto", choices=["auto", "python", "rust"]) + args, _ = parser.parse_known_args() + return args.backend + + +_requested_backend = _get_backend_from_args() +if _requested_backend in ("python", "rust"): + os.environ["DIFF_DIFF_BACKEND"] = _requested_backend + +# NOW import diff_diff and other dependencies (will see the env var) +import numpy as np # noqa: E402 +import pandas as pd # noqa: E402 + +# Add repo root to path for benchmarks.python.utils. When +# DIFF_DIFF_BENCH_USE_INSTALLED=1 (isolated-venv refresh runs), APPEND so the +# venv's installed diff-diff wheel wins over the dev tree; otherwise PREPEND +# (historical behavior: benchmark the working tree). +_REPO_ROOT = str(Path(__file__).parent.parent.parent) +if os.environ.get("DIFF_DIFF_BENCH_USE_INSTALLED") == "1": + sys.path.append(_REPO_ROOT) +else: + sys.path.insert(0, _REPO_ROOT) + +from benchmarks.python.utils import Timer, collect_provenance # noqa: E402 +from diff_diff import HAS_RUST_BACKEND, TwoWayFixedEffects # noqa: E402 +from diff_diff.utils import safe_inference # noqa: E402 + + +def parse_args(): + parser = argparse.ArgumentParser( + description="Benchmark TwoWayFixedEffects event-study mode (spec='within')" + ) + parser.add_argument("--data", required=True, help="Path to input CSV data") + parser.add_argument("--output", required=True, help="Path to output JSON results") + parser.add_argument("--cluster", default="unit", help="Column to cluster standard errors on") + parser.add_argument("--n-pre", type=int, required=True, help="Number of pre-treatment periods") + parser.add_argument( + "--n-post", type=int, required=True, help="Number of post-treatment periods" + ) + parser.add_argument( + "--reference-period", + type=int, + default=None, + help="Reference period (default: last pre-period = n_pre)", + ) + parser.add_argument( + "--backend", + default="auto", + choices=["auto", "python", "rust"], + help="Backend to use: auto (default), python (pure Python), rust (Rust backend)", + ) + parser.add_argument( + "--warmup", + action="store_true", + help="Run one untimed fit before the timed fit (JIT/cache warm-up)", + ) + return parser.parse_args() + + +def get_actual_backend() -> str: + """Return the actual backend being used based on HAS_RUST_BACKEND.""" + return "rust" if HAS_RUST_BACKEND else "python" + + +def post_period_average(es) -> Dict[str, float]: + """Headline post-period average of an ``EventStudyResults`` surface. + + Reproduces the former post-average arithmetic and NaN contract exactly as + ``tests/helpers/event_study_average.py``: ``att = mean(post att rows)``; + ``se = sqrt(sum(post vcov block) / n_post**2)`` over the Cartesian + ``np.ix_`` block; any NaN post ``att`` voids all five fields; a non-finite + or negative summed variance voids ``se`` (and hence p-value / CI); a zero + variance gives ``se = 0.0`` with an undefined (NaN) p-value / CI. The + p-value and CI come from the library's own ``safe_inference`` (all-NaN on + a non-positive or non-finite SE - never a finite p-value for undefined + inference) with the surface's shared post-row ``df`` (Student t; normal + when the df is undefined), matching the R script's + ``is.finite(avg_se) && avg_se > 0`` guard. + """ + nan5 = { + "att": float("nan"), + "se": float("nan"), + "pvalue": float("nan"), + "ci_lower": float("nan"), + "ci_upper": float("nan"), + } + if es.post_periods is None: + raise ValueError("post_period_average requires a surface with post_periods provenance") + post_rows = np.isin(es.event_time, list(es.post_periods)) + post_att = np.asarray(es.att, dtype=float)[post_rows] + if post_att.size == 0: + raise ValueError("post_period_average: no post-period rows on the surface") + if np.any(np.isnan(post_att)): + return nan5 + att_bar = float(np.mean(post_att)) + out = dict(nan5) + out["att"] = att_bar + if es.vcov is None or es.vcov_index is None: + return out + post_cols = np.isin(es.vcov_index, list(es.post_periods)) + n_post = int(post_cols.sum()) + if n_post == 0: + return out + avg_var = float(np.sum(np.asarray(es.vcov, dtype=float)[np.ix_(post_cols, post_cols)])) / ( + n_post**2 + ) + if not np.isfinite(avg_var) or avg_var < 0: + return out + se_bar = float(np.sqrt(avg_var)) + out["se"] = se_bar + # Shared post-row df of the surface (the fit's resolved inference df); + # undefined df -> None (normal reference inside safe_inference). + df = float(np.asarray(es.df, dtype=float)[post_rows][0]) + df_arg = df if np.isfinite(df) and df > 0 else None + _t, p_value, (ci_lo, ci_hi) = safe_inference(att_bar, se_bar, alpha=0.05, df=df_arg) + out["pvalue"] = float(p_value) + out["ci_lower"] = float(ci_lo) + out["ci_upper"] = float(ci_hi) + return out + + +def main(): + args = parse_args() + + # Get actual backend (already configured via env var before imports) + actual_backend = get_actual_backend() + print(f"Using backend: {actual_backend}") + + # Load data + print(f"Loading data from: {args.data}") + data = pd.read_csv(args.data) + + # Compute post_periods and reference_period from args + all_periods = sorted(data["time"].unique()) + n_pre = args.n_pre + post_periods = [p for p in all_periods if p > n_pre] + ref_period = args.reference_period if args.reference_period is not None else n_pre + + print(f"All periods: {all_periods}") + print(f"Post periods: {post_periods}") + print(f"Reference period: {ref_period}") + + fit_kwargs: Dict[str, Any] = dict( + outcome="outcome", + treatment="treated", + unit="unit", + time="time", + event_study=True, + spec="within", + post_periods=post_periods, + reference_period=ref_period, + ) + + # Run benchmark + print("Running TwoWayFixedEffects event-study estimation (spec='within')...") + + if args.warmup: + print("Warm-up fit (untimed)...") + TwoWayFixedEffects(cluster=args.cluster).fit(data, **fit_kwargs) + + twfe = TwoWayFixedEffects(cluster=args.cluster) + + with Timer() as timer: + results = twfe.fit(data, **fit_kwargs) + + total_time = timer.elapsed + + # Extract period effects (excluding reference period). The surface's + # event_time holds CALENDAR labels; R and compare_event_study join on + # reference-relative event_time, so convert here (as the frozen arm did). + period_effects = [] + for label, att, se, is_ref in zip( + results.event_time, results.att, results.se, results.is_reference + ): + if is_ref: + continue + period_effects.append( + { + "period": int(label), + "event_time": int(label - ref_period), + "att": float(att), + "se": float(se), + } + ) + period_effects.sort(key=lambda row: row["period"]) + + headline = post_period_average(results) + + # Build output + output = { + "estimator": "diff_diff.TwoWayFixedEffects[event_study, spec=within]", + "backend": actual_backend, + "cluster": args.cluster, + # Average treatment effect (across post-periods) + "att": headline["att"], + "se": headline["se"], + "pvalue": headline["pvalue"], + "ci_lower": headline["ci_lower"], + "ci_upper": headline["ci_upper"], + # Reference period + "reference_period": int(ref_period), + # Period-level effects + "period_effects": period_effects, + # Timing + "timing": { + "estimation_seconds": total_time, + "total_seconds": total_time, + }, + # Metadata + "metadata": { + "n_units": int(data["unit"].nunique()), + "n_periods": int(data["time"].nunique()), + "n_obs": len(data), + "n_pre": n_pre, + "n_post": len(post_periods), + "warmup": args.warmup, + }, + # Wheel/backend provenance (refresh runs hard-fail on mismatch) + "provenance": collect_provenance(), + } + + # Write output + print(f"Writing results to: {args.output}") + output_path = Path(args.output) + output_path.parent.mkdir(parents=True, exist_ok=True) + with open(output_path, "w") as f: + json.dump(output, f, indent=2) + + print(f"ATT: {headline['att']:.6f}") + print(f"SE: {headline['se']:.6f}") + print(f"Completed in {total_time:.3f} seconds") + return output + + +if __name__ == "__main__": + main() diff --git a/benchmarks/run_benchmarks.py b/benchmarks/run_benchmarks.py index 758946167..b354e669d 100644 --- a/benchmarks/run_benchmarks.py +++ b/benchmarks/run_benchmarks.py @@ -896,7 +896,13 @@ def run_multiperiod_benchmark( n_replications: int = 1, backends: Optional[List[str]] = None, ) -> Dict[str, Any]: - """Run MultiPeriodDiD event study benchmarks (Python and R) with replications.""" + """Run the TwoWayFixedEffects event-study benchmarks (Python and R) with replications. + + Live successor of the 3.x MultiPeriodDiD arm (removed in the 4.0 cut): the same + unit-FE event-study design, fit via ``TwoWayFixedEffects(...).fit(..., event_study=True, + spec="within")``. The dataset / result names keep the ``multiperiod`` key and the R + script is unchanged (the model is identical). + """ print(f"\n{'='*60}") print(f"MULTIPERIOD DID BENCHMARK ({scale})") print(f"{'='*60}") @@ -920,7 +926,8 @@ def run_multiperiod_benchmark( for backend in backends: backend_label = f"python_{'pure' if backend == 'python' else backend}" print( - f"\nRunning Python (diff_diff.MultiPeriodDiD, backend={backend}) - {n_replications} replications..." + f"\nRunning Python (diff_diff.TwoWayFixedEffects[event_study, spec=within], " + f"backend={backend}) - {n_replications} replications..." ) py_output = RESULTS_DIR / "accuracy" / f"{backend_label}_{name}_{scale}.json" py_output.parent.mkdir(parents=True, exist_ok=True) @@ -930,7 +937,7 @@ def run_multiperiod_benchmark( for rep in range(n_replications): try: py_result = run_python_benchmark( - "benchmark_multiperiod.py", + "benchmark_twfe_event_study.py", data_path, py_output, extra_args=extra_args, @@ -996,7 +1003,7 @@ def run_multiperiod_benchmark( comparison = compare_estimates( results["python"], results["r"], - "MultiPeriodDiD", + "TWFE[event_study]", scale=scale, se_rtol=0.01, python_pure_results=results.get("python_pure"), diff --git a/benchmarks/speed_review/bench_brand_awareness_survey.py b/benchmarks/speed_review/bench_brand_awareness_survey.py index f4178f0ea..71c00c578 100644 --- a/benchmarks/speed_review/bench_brand_awareness_survey.py +++ b/benchmarks/speed_review/bench_brand_awareness_survey.py @@ -21,8 +21,8 @@ from diff_diff import ( DifferenceInDifferences, - MultiPeriodDiD, SurveyDesign, + TwoWayFixedEffects, check_parallel_trends, compute_honest_did, ) @@ -125,10 +125,17 @@ def placebo_refit(): ) def honest_did_grid(): - es = MultiPeriodDiD() + # Pooled event study (repeated cross-section form, no unit FE): the + # 3.x MultiPeriodDiD call this replaces passed unit= but only the + # Conley/cluster paths read it, so the measured design was the + # pooled one. Cohort onset is period 7 (cohort_periods=[7] above), + # so post_periods = {7..12} with the last pre-period as reference. + es = TwoWayFixedEffects() es_result = es.fit( data, outcome="outcome", treatment="treat_unit", - time="period", unit="unit", reference_period=6, + time="period", event_study=True, spec="pooled", + reference_period=6, + post_periods=[p for p in sorted(data["period"].unique()) if p >= 7], survey_design=sd_tsl, ) results["event_study"] = es_result diff --git a/changelog.d/20261025-v4-pr4c-ii-b-docs-benchmarks.md b/changelog.d/20261025-v4-pr4c-ii-b-docs-benchmarks.md new file mode 100644 index 000000000..23a630097 --- /dev/null +++ b/changelog.d/20261025-v4-pr4c-ii-b-docs-benchmarks.md @@ -0,0 +1,12 @@ +### Internal +- PR 4c-ii-b (docs/benchmarks half of the MultiPeriodDiD removal): REGISTRY / REPORTING / + survey-theory / METHODOLOGY_REVIEW and the remaining user-doc prose stop describing the + removed class as live; TODO rows 33/35/38/95/99 and three DEFERRED rows that governed the + class are reworded or retired. The oracle fixture `tests/fixtures/v4_twfe_event_study_oracles. + json` drops its `mpd_twin` records and the lane metadata loses its MPD twin fields (the TWFE + records are byte-identical); the capture script documents the per-lane capture SHAs. +- Benchmarks: the live event-study arm is the new `benchmarks/python/benchmark_twfe_event_study.py` + (`spec="within"`, `cluster="unit"`; parity against R `fixest::feols(... | unit)` verified + locally - ATT to 2e-12, SE to 0.0%, per-period effects correlation 1.000); + `benchmark_multiperiod.py` and `benchmarks/refresh_2026_07/` stay frozen because that harness + runs inside pinned 3.5.3/3.7.0 venvs that predate `event_study=True`. diff --git a/docs/benchmarks.rst b/docs/benchmarks.rst index 8d10aa7bf..9bae7550d 100644 --- a/docs/benchmarks.rst +++ b/docs/benchmarks.rst @@ -37,9 +37,9 @@ where no runnable R reference exists: * - ``CallawaySantAnna`` - ``did::att_gt`` - Callaway & Sant'Anna (2021) - * - ``MultiPeriodDiD`` + * - ``TwoWayFixedEffects`` event-study mode (``spec="within"``) - ``fixest::feols`` - - Event study with treatment × period interactions + - Event study with treatment × period interactions (unit FE absorbed) * - ``SyntheticDiD`` - ``synthdid::synthdid_estimate`` - Arkhangelsky et al. (2021) @@ -155,6 +155,11 @@ registry note (benchmark SE gate is Monte Carlo-bounded). .. refresh-table-end: summary +The ``MultiPeriodDiD`` row above is the committed 2026-07 snapshot, run on the +3.7.0 wheel through the ``MultiPeriodDiD`` class removed at the 4.0 cut; the +live arm is the ``TwoWayFixedEffects`` event-study mode (``spec="within"``), +which runs the same shared core. + Basic DiD Results ~~~~~~~~~~~~~~~~~ @@ -191,11 +196,16 @@ Basic DiD Results **Validation**: PASS - Results are numerically identical across all implementations. -MultiPeriodDiD Results -~~~~~~~~~~~~~~~~~~~~~~ +TWFE Event-Study Results +~~~~~~~~~~~~~~~~~~~~~~~~ **Data**: 200 units, 8 periods (4 pre, 4 post), true ATT = 3.0 (small scale) +The table below is the committed 2026-07 snapshot, run on the 3.7.0 wheel +through the ``MultiPeriodDiD`` class removed at the 4.0 cut; the live arm is +``TwoWayFixedEffects().fit(..., event_study=True, spec="within")``, which +runs the same shared core. + .. refresh-table-start: accuracy_multiperiod .. list-table:: @@ -578,7 +588,7 @@ Dataset Sizes * - Scale - BasicDiD - - MultiPeriodDiD + - TWFE event study - CallawaySantAnna - SyntheticDiD - Observations @@ -1114,9 +1124,12 @@ When to Trust Results - **BasicDiD/TWFE**: Results are identical to R. Use with confidence. -- **MultiPeriodDiD**: Results are identical to R's ``fixest::feols`` with - ``treated * time_f | unit`` interaction syntax (unit FE absorbed). Both average - ATT and all period-level effects match to machine precision. Use with confidence. +- **TWFE event-study mode** (``spec="within"``): Results are identical to R's + ``fixest::feols`` with ``treated * time_f | unit`` interaction syntax (unit FE + absorbed). Both the post-period average and all period-level effects match to + machine precision (the committed snapshot ran through the removed + ``MultiPeriodDiD`` class; the live arm runs the same shared core). Use with + confidence. - **SyntheticDiD**: Point estimates are numerically identical (< 1e-10 diff) and standard errors agree within placebo Monte Carlo dispersion. Both diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 3f6eeb0a2..56905b852 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -2,7 +2,7 @@ This document provides the academic foundations and key implementation requirements for each estimator in diff-diff. It serves as a reference for contributors and users who want to understand the theoretical basis of the methods. -**Result-class field naming.** Headline scalar inference fields appear under one of four native naming patterns: flat `att` / `se` / `conf_int` / `p_value` / `t_stat` (`DiDResults`, `SyntheticDiDResults`, `TROPResults`, `TripleDifferenceResults`, `HeterogeneousAdoptionDiDResults`, `ChangesInChangesResults`); `overall_*` (`CallawaySantAnnaResults` and the rest of the staggered family); `overall_att_*` (`ContinuousDiDResults`, where `att` and `acrt` are parallel response curves); and `avg_*` (`MultiPeriodDiDResults`). Result classes in the prefixed `overall_*` / `overall_att_*` / `avg_*` families additionally expose the flat `att` / `se` / `conf_int` / `p_value` / `t_stat` names as read-only `@property` aliases over their canonical fields, for adapter / external-consumer compatibility (see PR for v3.3.3, motivated by `balance.interop.diff_diff`). The flat-native classes (`DiDResults`, `SyntheticDiDResults`, `TROPResults`, `TripleDifferenceResults`, `HeterogeneousAdoptionDiDResults`, `ChangesInChangesResults`) already carry these names as native dataclass fields and are unchanged by this contract. `ContinuousDiDResults` further exposes `overall_*` aliases pointing at the ATT side (so `result.overall_se` reads `result.overall_att_se`, etc.). The native field is canonical for documentation, semantics, and computation — aliases are pure read-throughs and inherit the `safe_inference()` joint-NaN consistency contract automatically. Because aliases are `@property` descriptors (not dataclass fields), they do NOT appear in `dataclasses.fields()` or `dataclasses.asdict()` output, and assignment to an alias raises `AttributeError`; serializers and field-walkers continue to see only the canonical field set. +**Result-class field naming.** Headline scalar inference fields appear under one of three native naming patterns: flat `att` / `se` / `conf_int` / `p_value` / `t_stat` (`DiDResults`, `SyntheticDiDResults`, `TROPResults`, `TripleDifferenceResults`, `HeterogeneousAdoptionDiDResults`, `ChangesInChangesResults`); `overall_*` (`CallawaySantAnnaResults` and the rest of the staggered family); and `overall_att_*` (`ContinuousDiDResults`, where `att` and `acrt` are parallel response curves) - the fourth pattern, `avg_*`, left with the removed `MultiPeriodDiDResults` at 4.0 (the TWFE event-study surface carries no post-period-average scalar). Result classes in the prefixed `overall_*` / `overall_att_*` families additionally expose the flat `att` / `se` / `conf_int` / `p_value` / `t_stat` names as read-only `@property` aliases over their canonical fields, for adapter / external-consumer compatibility (see PR for v3.3.3, motivated by `balance.interop.diff_diff`). The flat-native classes (`DiDResults`, `SyntheticDiDResults`, `TROPResults`, `TripleDifferenceResults`, `HeterogeneousAdoptionDiDResults`, `ChangesInChangesResults`) already carry these names as native dataclass fields and are unchanged by this contract. `ContinuousDiDResults` further exposes `overall_*` aliases pointing at the ATT side (so `result.overall_se` reads `result.overall_att_se`, etc.). The native field is canonical for documentation, semantics, and computation — aliases are pure read-throughs and inherit the `safe_inference()` joint-NaN consistency contract automatically. Because aliases are `@property` descriptors (not dataclass fields), they do NOT appear in `dataclasses.fields()` or `dataclasses.asdict()` output, and assignment to an alias raises `AttributeError`; serializers and field-walkers continue to see only the canonical field set. ## Table of Contents @@ -95,9 +95,9 @@ where τ is the ATT. - Rank-deficient design matrix (collinearity): warns and sets NA for dropped coefficients (R-style, matches `lm()`) - Tolerance: `1e-07` (matches R's `qr()` default), relative to largest diagonal element of R in QR decomposition - Controllable via `rank_deficient_action` parameter: "warn" (default), "error", or "silent" - - **Note (scale invariance):** the shared `diff_diff/linalg.py` rank check re-detects on column-equilibrated (unit 2-norm) columns when the raw pass reports a deficiency (adopting the higher equilibrated rank and its scale-corrected pivot selection only when a large-scale column inflated the threshold; genuine collinearity keeps the raw selection), and the least-squares solve equilibrates columns then unscales the coefficients — so rank detection and the fit are invariant to per-column scaling; for a well-scaled collinear design the dropped column is unchanged, while a scale-induced under-count adopts the scale-corrected equilibrated selection (which may differ from the raw choice but is guaranteed to retain an identified subset). This covers every covariate outcome-regression fit routed through `solve_ols` — DiD, TwoWayFixedEffects, MultiPeriodDiD, ImputationDiD, TwoStageDiD, TripleDifference, and (as of the OR scale-equilibration change) CallawaySantAnna's `_compute_all_att_gt_covariate_reg` / `_doubly_robust` and the staggered DDD engine's `_compute_or` (`_staggered_triple_diff_engine.py`, TripleDifference staggered mode) — see the CallawaySantAnna rank-deficiency Note's scope statement. **Certification fast path (2026-07):** the common full-rank case is short-circuited by a Gram/eigenvalue certification on the equilibrated Gram at the documented `_rank_guarded_inv` 1e-10 Gram threshold — two orders stricter than the QR full-rank boundary, so it never reports full rank where the pivoted-QR path would report a deficiency; deficient or uncertifiable designs (including all n < k and non-finite inputs) fall through to the pivoted-QR path verbatim, leaving every drop decision, pivot selection, and rank count unchanged. Certification never SELECTS columns (it only confirms none need dropping), so the `_rank_guarded_inv` note about equilibrated-Gram column selection being reserved for the IF inverse still holds. + - **Note (scale invariance):** the shared `diff_diff/linalg.py` rank check re-detects on column-equilibrated (unit 2-norm) columns when the raw pass reports a deficiency (adopting the higher equilibrated rank and its scale-corrected pivot selection only when a large-scale column inflated the threshold; genuine collinearity keeps the raw selection), and the least-squares solve equilibrates columns then unscales the coefficients — so rank detection and the fit are invariant to per-column scaling; for a well-scaled collinear design the dropped column is unchanged, while a scale-induced under-count adopts the scale-corrected equilibrated selection (which may differ from the raw choice but is guaranteed to retain an identified subset). This covers every covariate outcome-regression fit routed through `solve_ols` — DiD, TwoWayFixedEffects (both modes), ImputationDiD, TwoStageDiD, TripleDifference, and (as of the OR scale-equilibration change) CallawaySantAnna's `_compute_all_att_gt_covariate_reg` / `_doubly_robust` and the staggered DDD engine's `_compute_or` (`_staggered_triple_diff_engine.py`, TripleDifference staggered mode) — see the CallawaySantAnna rank-deficiency Note's scope statement. **Certification fast path (2026-07):** the common full-rank case is short-circuited by a Gram/eigenvalue certification on the equilibrated Gram at the documented `_rank_guarded_inv` 1e-10 Gram threshold — two orders stricter than the QR full-rank boundary, so it never reports full rank where the pivoted-QR path would report a deficiency; deficient or uncertifiable designs (including all n < k and non-finite inputs) fall through to the pivoted-QR path verbatim, leaving every drop decision, pivot selection, and rank count unchanged. Certification never SELECTS columns (it only confirms none need dropping), so the `_rank_guarded_inv` note about equilibrated-Gram column selection being reserved for the IF inverse still holds. - **Note (run-to-run determinism of the clustered vcov, 2026-07):** the Rust backend's clustered `compute_robust_vcov` / `solve_ols(return_vcov=True)` previously aggregated cluster scores in `HashMap` iteration order — SipHash-randomized per call — so identical inputs produced ~1e-14-level distinct vcov values across calls (the Python backend was always bit-stable). Fixed in #653: rows accumulate in first-appearance order (ascending for the factorized 0..G-1 ids the dispatcher passes, matching NumPy's groupby order); repeated identical calls are bit-identical (a 50-call probe during the fix; the shipped regression asserts 20 repeats per id-layout), NumPy parity unchanged (~1e-14 cross-backend GEMM tolerance). Regression: `tests/test_rust_backend.py::TestClusterVcovDeterminism` (contiguous + non-contiguous unsorted ids, plus a raw-kernel path). -- **Note (covariate-name collision guard):** `DifferenceInDifferences`, `MultiPeriodDiD`, and `TwoWayFixedEffects` build their `coefficients` dict by zipping a variable-name list (structural terms + user covariate column names, appended verbatim) with the coefficient vector. A covariate whose name equals a reserved structural name — `const`, the treatment/time column names, the `{treatment}:{time}` interaction (DiD), the `period_{p}` dummies / `{treatment}:period_{p}` interactions (MultiPeriodDiD), `ATT` (TwoWayFixedEffects), any fixed-effect / unit / time dummy name, or an internal `_`-prefixed working column (`_treat_time`, `_did_treatment`, `_treatment_post`) — would silently overwrite that structural coefficient (dict last-write-wins) with no error. `fit()` now calls `utils.validate_covariate_names()` before building the design, raising `ValueError` on such a collision (case-sensitive, so `Const` is allowed) and on duplicate covariate names. Reserved fixed-effect/unit/time dummy names are taken from the same `pd.get_dummies(..., drop_first=True)` call used to build them (exact match, including for `Categorical` columns). The TwoWayFixedEffects guard fires on both variance paths — the within-transform path exposes only `{"ATT": att}`, but a `_treatment_post`-named covariate would still clobber the internal interaction column. The influence-function estimators (CallawaySantAnna, etc.) key results by `(g, t)` and are unaffected. +- **Note (covariate-name collision guard):** `DifferenceInDifferences` and `TwoWayFixedEffects` (static and event-study modes) build their `coefficients` dict by zipping a variable-name list (structural terms + user covariate column names, appended verbatim) with the coefficient vector. A covariate whose name equals a reserved structural name — `const`, the treatment/time column names, the `{treatment}:{time}` interaction (DiD), the `period_{p}` dummies / `{treatment}:period_{p}` interactions (TwoWayFixedEffects event-study mode, the removed MultiPeriodDiD design), `ATT` (TwoWayFixedEffects static mode), any fixed-effect / unit / time dummy name, or an internal `_`-prefixed working column (`_treat_time`, `_did_treatment`, `_treatment_post`) — would silently overwrite that structural coefficient (dict last-write-wins) with no error. `fit()` now calls `utils.validate_covariate_names()` before building the design, raising `ValueError` on such a collision (case-sensitive, so `Const` is allowed) and on duplicate covariate names. Reserved fixed-effect/unit/time dummy names are taken from the same `pd.get_dummies(..., drop_first=True)` call used to build them (exact match, including for `Categorical` columns). The TwoWayFixedEffects guard fires on both variance paths — the within-transform path exposes only `{"ATT": att}`, but a `_treatment_post`-named covariate would still clobber the internal interaction column. The influence-function estimators (CallawaySantAnna, etc.) key results by `(g, t)` and are unaffected. **Reference implementation(s):** - R: `fixest::feols()` with interaction term @@ -115,8 +115,9 @@ where τ is the ATT. bootstrap of Cameron, Gelbach & Miller (2008), matching the defaults of R's `fwildclusterboot::boottest` (Roodman, MacKinnon, Nielsen & Webb 2019). Implemented in `diff_diff.utils.wild_bootstrap_se`; used by `DifferenceInDifferences` and (inherited) -`TwoWayFixedEffects`. (`MultiPeriodDiD` does **not** support it — it falls back to analytical -inference and the inherited `p_val_type` is inert there.) +`TwoWayFixedEffects` (static mode; the event-study mode raises `ValueError` - see the +TwoWayFixedEffects event-study section. The removed MultiPeriodDiD fell back to analytical +inference through 3.x.) - **Note:** Since 3.9 the selector fails closed (rows M-081/M-096): `inference=` accepts exactly `{"analytical", "wild_bootstrap"}` (string-typed) at construction @@ -127,7 +128,7 @@ inference and the inherited `p_val_type` is inert there.) (0 or 1 replications deterministically produce the all-NaN degenerate tuple). `TwoWayFixedEffects`' unit auto-cluster satisfies the cluster prerequisite; the survey and Conley `NotImplementedError` rejections keep precedence. - `MultiPeriodDiD` keeps its warn-and-analytical-fallback until its 4.0 removal. + The removed MultiPeriodDiD's warn-and-analytical-fallback went with the class at 4.0. *Algorithm (test of H₀: τ = r, default r = 0):* 1. **Impose the null** by dropping the interaction column and re-fitting the reduced model; @@ -336,7 +337,9 @@ where V is the VCV sub-matrix for post-treatment δ_e coefficients. 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 + `tests/test_variance_conventions.py` through 3.x; since 4.0 the live rows + are `twfe_es_within_hc1_cluster_unit` / `twfe_es_pooled_hc1_cluster_time` - + the `fixed_effects=` spelling died with the class). 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. @@ -399,7 +402,7 @@ where V is the VCV sub-matrix for post-treatment δ_e coefficients. - [ ] Supports unit and time FE via absorption - [ ] Endpoint binning for distant event times - [x] Average ATT correctly accounts for covariance between period effects -- [x] Returns PeriodEffect objects with confidence intervals +- [x] Returned `PeriodEffect` records with confidence intervals (through 3.x; since 4.0 the per-period rows live on `EventStudyResults`, with `period_effects` surviving as the M-016 FutureWarning view) - [x] Supports both balanced and unbalanced panels --- @@ -556,7 +559,8 @@ 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** CR1 SEs converge + and the TWFE event-study `spec="within"` lane with `vcov_type in {classical,hc1}` (the removed + `MultiPeriodDiD(absorb=...)`). **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 @@ -600,7 +604,7 @@ This matches the behavior of R's `fixest::feols()` with absorbed FE. convergence above (the historical ~0.25% non-nested-FE band is closed). **Three-value knob (2026-07, widened 3.9 / M-127):** `df_convention ∈ {"residual", "cluster", "normal"}` on the library-wide knob roster — `DifferenceInDifferences` / `TwoWayFixedEffects` / - `MultiPeriodDiD` / `LinearRegression` plus (since 3.9) `SunAbraham`, `WooldridgeDiD` (OLS + `LinearRegression` (the TWFE event-study mode included) plus (since 3.9) `SunAbraham`, `WooldridgeDiD` (OLS arms), `StackedDiD`, `ImputationDiD` (pretrends lead regression), and `LPDiD` (default `"cluster"`, its Stata-anchored convention). `"cluster"` switches clustered analytical t/p/CI to the Stata/fixest `G − 1` convention and is inert on unclustered and conley fits; @@ -705,7 +709,7 @@ distinction governs which `vcov_type` values an estimator can accept: regression and derive variance via `solve_ols(..., vcov_type=...)`, returning a sandwich `(X'X)^{-1} M (X'X)^{-1}` whose meat `M` is parameterized by `vcov_type ∈ {classical, hc1, hc2, hc2_bm, hc3}` (plus `conley` for spatial-HAC). `hc3` entered the shared vocabulary with the LWDiD contribution; estimators that do not implement it reject it explicitly (structural roster guard in `tests/test_estimators_vcov_type.py`). -Examples: `DifferenceInDifferences`, `MultiPeriodDiD`, `TwoWayFixedEffects`, +Examples: `DifferenceInDifferences`, `TwoWayFixedEffects` (both modes), `SunAbraham`, `StackedDiD`, `WooldridgeDiD`, `LinearRegression`. The full `vcov_type` contract is methodologically applicable because every family has a defined interpretation on the hat-matrix-bearing design (HC2 leverage @@ -1184,7 +1188,7 @@ The multiplier bootstrap uses random weights w_i with E[w]=0 and Var(w)=1: - Detection: Pivoted QR decomposition with tolerance `1e-07` (R's `qr()` default), with a **column-equilibration re-check** (unit 2-norm) that makes the rank *count* invariant to per-column scaling; the dropped-column selection is unchanged for well-scaled collinear designs (a scale-induced under-count instead adopts the scale-corrected equilibrated selection) - Handling: Warns and drops linearly dependent columns, sets NA for dropped coefficients (R-style, matches `lm()`) - Parameter: `rank_deficient_action` controls behavior: "warn" (default), "error", or "silent" - - **Note:** Rank detection and the least-squares solve are invariant to per-column scaling in BOTH the Python and Rust backends. The rank threshold is anchored to the largest pivot/singular value, so a column on a large raw scale (e.g. an unstandardized covariate in raw population or currency units) previously inflated the threshold and false-dropped the intercept/treatment/interaction to NaN on an otherwise full-rank design — or truncated the small-scale direction in the solve, returning finite-but-wrong coefficients. Detection now runs a raw pivoted QR first and only re-checks on column-equilibrated (unit 2-norm) columns when the raw pass reports a deficiency. If the equilibrated rank is higher (a scale-induced false-drop), the equilibrated rank AND its pivot selection are adopted — its first `rank` columns are independent under the scale-corrected criterion, so the retained design is guaranteed identified; otherwise (genuine collinearity, no scale disparity) the raw rank and pivot selection are kept unchanged. The solve equilibrates columns and unscales the coefficients. The Python (LAPACK `gelsd`) and Rust (faer thin-SVD) backends implement this same contract and agree at the parity-suite tolerances (coefficients ~1e-8, vcov ~1e-5), never bit-identically. This repairs the scale bug while leaving everything else unchanged: it is a no-op for full-rank well-conditioned designs (R-parity unaffected) and does **not** change which column is dropped in a *well-scaled* collinear design (the established raw pivot selection is preserved). A scale-induced under-count instead adopts the scale-corrected equilibrated selection — which may differ from the raw choice in a mixed scale+collinearity design but is guaranteed to retain an identified (full-rank) subset. This shared `diff_diff/linalg.py` behavior covers every covariate outcome-regression fit routed through `solve_ols` — DiD, TwoWayFixedEffects, MultiPeriodDiD, ImputationDiD, TwoStageDiD, and TripleDifference. **Scope (now scale-equilibrated):** CallawaySantAnna's covariate outcome-regression (`_compute_all_att_gt_covariate_reg`) and doubly-robust (`_doubly_robust`) nuisance solves — and the staggered DDD engine's per-cohort OR solve (`_compute_or`, TripleDifference staggered mode) — now route through `solve_ols` (column-equilibrated SVD/gelsd), matching TripleDifference's already-`solve_ols`-routed OR fit (`triple_diff.py:1438/1444`) and R's `lm()`/QR. The prior estimator-local `cho_solve(X'X)` / `scipy.linalg.lstsq(cond=1e-7)` fast paths were not scale-equilibrated: a covariate **correlated with another regressor at a very large scale** (e.g. a large constant offset, near-collinear with the intercept) could perturb the point-estimate ATT — and the IF SE that follows it — because the normal-equations Cholesky squares the condition number (pure orthogonal ill-scaling was already safe). The equilibrated SVD is offset-invariant to ~1e-11 where the prior solve drifted; scale-invariance is now pinned by tests (covariate + 1e6 offset → ATT(g,t) unchanged). The change is not bit-identical (cho/normal-equations → SVD) but well-scaled designs move only ~1e-12. The separate DR/OR influence-function SE rank-guard — which previously returned garbage SEs (~1e13) when these local Gram matrices were near-singular — is also **implemented** for CS / TripleDifference (both modes) via `_rank_guarded_inv` (see the "rank-guarded IF standard errors" Note above). + - **Note:** Rank detection and the least-squares solve are invariant to per-column scaling in BOTH the Python and Rust backends. The rank threshold is anchored to the largest pivot/singular value, so a column on a large raw scale (e.g. an unstandardized covariate in raw population or currency units) previously inflated the threshold and false-dropped the intercept/treatment/interaction to NaN on an otherwise full-rank design — or truncated the small-scale direction in the solve, returning finite-but-wrong coefficients. Detection now runs a raw pivoted QR first and only re-checks on column-equilibrated (unit 2-norm) columns when the raw pass reports a deficiency. If the equilibrated rank is higher (a scale-induced false-drop), the equilibrated rank AND its pivot selection are adopted — its first `rank` columns are independent under the scale-corrected criterion, so the retained design is guaranteed identified; otherwise (genuine collinearity, no scale disparity) the raw rank and pivot selection are kept unchanged. The solve equilibrates columns and unscales the coefficients. The Python (LAPACK `gelsd`) and Rust (faer thin-SVD) backends implement this same contract and agree at the parity-suite tolerances (coefficients ~1e-8, vcov ~1e-5), never bit-identically. This repairs the scale bug while leaving everything else unchanged: it is a no-op for full-rank well-conditioned designs (R-parity unaffected) and does **not** change which column is dropped in a *well-scaled* collinear design (the established raw pivot selection is preserved). A scale-induced under-count instead adopts the scale-corrected equilibrated selection — which may differ from the raw choice in a mixed scale+collinearity design but is guaranteed to retain an identified (full-rank) subset. This shared `diff_diff/linalg.py` behavior covers every covariate outcome-regression fit routed through `solve_ols` — DiD, TwoWayFixedEffects (both modes), ImputationDiD, TwoStageDiD, and TripleDifference. **Scope (now scale-equilibrated):** CallawaySantAnna's covariate outcome-regression (`_compute_all_att_gt_covariate_reg`) and doubly-robust (`_doubly_robust`) nuisance solves — and the staggered DDD engine's per-cohort OR solve (`_compute_or`, TripleDifference staggered mode) — now route through `solve_ols` (column-equilibrated SVD/gelsd), matching TripleDifference's already-`solve_ols`-routed OR fit (`triple_diff.py:1438/1444`) and R's `lm()`/QR. The prior estimator-local `cho_solve(X'X)` / `scipy.linalg.lstsq(cond=1e-7)` fast paths were not scale-equilibrated: a covariate **correlated with another regressor at a very large scale** (e.g. a large constant offset, near-collinear with the intercept) could perturb the point-estimate ATT — and the IF SE that follows it — because the normal-equations Cholesky squares the condition number (pure orthogonal ill-scaling was already safe). The equilibrated SVD is offset-invariant to ~1e-11 where the prior solve drifted; scale-invariance is now pinned by tests (covariate + 1e6 offset → ATT(g,t) unchanged). The change is not bit-identical (cho/normal-equations → SVD) but well-scaled designs move only ~1e-12. The separate DR/OR influence-function SE rank-guard — which previously returned garbage SEs (~1e13) when these local Gram matrices were near-singular — is also **implemented** for CS / TripleDifference (both modes) via `_rank_guarded_inv` (see the "rank-guarded IF standard errors" Note above). - **Note (opt-in solve_ols Cholesky fast path, 2026-07):** setting `DIFF_DIFF_SOLVE_OLS_FASTPATH=1` (positive integer; read per call) routes full-rank `solve_ols` fits through an equilibrated normal-equations Cholesky solve, certification-gated at reciprocal condition 1e-6 on the equilibrated Gram (Python: LAPACK `dpocon`, the same guard shape as the IRLS inner-solver Note under Propensity score estimation; Rust: the exact 1-norm rcond, whose inverse byproduct is reused as the sandwich bread). Any certification decline falls back VERBATIM to the SVD/gelsd solve, so the knob-off default path - and every rank-deficiency decision - stays byte-identical to the legacy behavior on both backends. The fast path never runs on a design the QR detector would call deficient: the 1e-6 solve guard is strictly inside the 1e-10 stage-0 rank certification, which is itself two orders stricter than the 1e-7 QR boundary. Within the opt-in path, parity vs the default solve is tol-bounded, not bit-level (fitted values ~1e-8 abs, SEs ~1e-6 rel; certified forward-error budget ~eps*cond(G_eq) <= 2e-10 in the equilibrated basis). The default path keeps the SVD's independent singular-value truncation as defense-in-depth; the knob is opt-in per the 2026-07 correctness-first decision. - Non-finite inference values: - Analytic SE: Returns NaN to signal invalid inference (not biased via zeroing) @@ -1326,7 +1330,7 @@ The multiplier bootstrap uses random weights w_i with E[w]=0 and Var(w)=1: - [x] Repeated cross-sections (`panel=False`) for non-panel surveys (Phase 7b) - **Note:** `anticipation` is validated at construction (non-negative integer; `bool` rejected) via the shared `utils.validate_anticipation`, and re-checked on the fit path (direct-mutation defense) — see the family-wide adoption note (ledger row [M-144]) in the TripleDifference staggered-mode section. -- **Note (`summary(alpha=)` never recomputes stored inference, ledger row M-146):** `CallawaySantAnnaResults.summary`/`print_summary` — and the sibling staggered-family results classes (StaggeredTripleDiffResults, ChaisemartinDHaultfoeuilleResults, ImputationDiDResults, EfficientDiDResults, TwoStageDiDResults, StackedDiDResults, SunAbrahamResults) — reject an `alpha` different from the fit-time value via the shared `results_base._require_fit_alpha` guard (the DMLDiDResults/EventStudyResults precedent). Previously the requested alpha relabeled the confidence-interval header while the FIT-TIME stored intervals were printed — silent coverage mislabeling (bootstrap percentile intervals cannot be reconstructed from the SE at all); `alpha=0.0` was additionally swallowed by a falsy-`or` default and now raises. Re-fit at the desired alpha instead. The guard is now FAMILY-WIDE: the non-staggered summaries adopted it too (DiDResults incl. SpilloverDiDResults by inheritance, MultiPeriodDiDResults, SyntheticDiDResults, TripleDifferenceResults, TROPResults, ContinuousDiDResults, SyntheticControlResults), closing the audit row this sentence used to track. +- **Note (`summary(alpha=)` never recomputes stored inference, ledger row M-146):** `CallawaySantAnnaResults.summary`/`print_summary` — and the sibling staggered-family results classes (StaggeredTripleDiffResults, ChaisemartinDHaultfoeuilleResults, ImputationDiDResults, EfficientDiDResults, TwoStageDiDResults, StackedDiDResults, SunAbrahamResults) — reject an `alpha` different from the fit-time value via the shared `results_base._require_fit_alpha` guard (the DMLDiDResults/EventStudyResults precedent). Previously the requested alpha relabeled the confidence-interval header while the FIT-TIME stored intervals were printed — silent coverage mislabeling (bootstrap percentile intervals cannot be reconstructed from the SE at all); `alpha=0.0` was additionally swallowed by a falsy-`or` default and now raises. Re-fit at the desired alpha instead. The guard is now FAMILY-WIDE: the non-staggered summaries adopted it too (DiDResults incl. SpilloverDiDResults by inheritance, the since-removed MultiPeriodDiDResults, SyntheticDiDResults, TripleDifferenceResults, TROPResults, ContinuousDiDResults, SyntheticControlResults), closing the audit row this sentence used to track. - **Note (SyntheticControl summary-alpha posture):** `SyntheticControlResults.summary(alpha=)` applies the same guard although the class prints no alpha-based Wald interval (`se`/`t_stat`/`p_value`/`conf_int` NaN by design; the displayed Firpo-Possebom confidence set is keyed on its own stored `gamma`): the previously dead `alpha` parameter now raises on a non-fit value, with a class-accurate message, instead of silently doing nothing. - **Note (TROP summary-alpha posture):** `TROPResults.summary(alpha=)` applies the guard purely by the family-wide never-recompute contract - its stored interval is a t interval at the fit alpha around a bootstrap SE and WOULD be arithmetically reconstructible; the tailored error message states the uniform-contract rationale rather than a non-reconstructibility claim. --- @@ -1978,7 +1982,7 @@ where weights ŵ_{g,e} = n_{g,e} / Σ_g n_{g,e} (sample share of cohort g at eve Per-event-time contrast `c_e[full_idx(g,e)] = w_{g,e}` (IW weight) and overall ATT contrast `c_overall[full_idx(g,e)] = w_e × w_{g,e}` are passed to `_compute_cr2_bm_contrast_dof` (the helper PR #465 - added for MultiPeriodDiD's post-period-average DOF). The resulting + added for the removed MultiPeriodDiD's post-period-average DOF). The resulting per-contrast DOF threads into `safe_inference(..., df=)`. Matches `clubSandwich::Wald_test(constraints=matrix(c, 1), test="HTZ")$df_denom` at atol=1e-10 (pinned in @@ -2582,7 +2586,7 @@ where `g(·)` is the link inverse (logistic or exp), `η_i` is the individual li - `method ∈ {"logit","poisson"}` + `vcov_type != "hc1"` — REJECTED at `__init__`. GLM QMLE sandwich with HC2 leverage on canonical-link pseudo-residuals (`w = p(1-p)` for logit, `w = μ_i` for Poisson) needs CR2-BM-on-GLM derivation + R parity against `clubSandwich::vcovCR(glm(...))`. Tracked in DEFERRED.md (WooldridgeDiD logit/poisson follow-up row). - `survey_design=` + `vcov_type != "hc1"` — REJECTED at `fit()` with `NotImplementedError`. Survey TSL/replicate-refit overrides analytical sandwich. Use `vcov_type="hc1"` (default) for survey designs. - `n_bootstrap > 0` + `vcov_type ∈ {"hc2","classical"}` — REJECTED at `fit()` regardless of `self.cluster` setting. The multiplier bootstrap is intrinsically clustered, but one-way vcov_type does not compose with `cluster_ids`: with `cluster=None` the auto-cluster is dropped (bootstrap has no cluster to draw weights at); with `cluster=X` the linalg validator rejects one-way + cluster_ids downstream with a less-informative error. User must drop bootstrap (`n_bootstrap=0`) or pick a cluster-compatible `vcov_type` (`hc1` or `hc2_bm`). -- **Note:** This routing is a documented synthesis of two existing methodology ingredients: the full-dummy auto-route from the Phase 1b PR 1/8 SunAbraham pattern (PR #472, which itself reused the Phase 1a Gate 1 TWFE lift from PR #469), and the clubSandwich WLS-CR2 algebra from the Phase 1a port (PR #475). The BM contrast DOF threading reuses `_compute_cr2_bm_contrast_dof` from PR #465 (MPD). No new methodology choice is introduced — the change is purely surface: extending the existing pattern from SA-OLS to WooldridgeDiD-OLS. +- **Note:** This routing is a documented synthesis of two existing methodology ingredients: the full-dummy auto-route from the Phase 1b PR 1/8 SunAbraham pattern (PR #472, which itself reused the Phase 1a Gate 1 TWFE lift from PR #469), and the clubSandwich WLS-CR2 algebra from the Phase 1a port (PR #475). The BM contrast DOF threading reuses `_compute_cr2_bm_contrast_dof` from PR #465 (the removed MPD's post-average DOF helper). No new methodology choice is introduced — the change is purely surface: extending the existing pattern from SA-OLS to WooldridgeDiD-OLS. - **Note:** Bootstrap is supported only with `vcov_type ∈ {"hc1","hc2_bm"}` (one-way `classical`/`hc2` + bootstrap is rejected at `fit()` per the previous bullet). On the supported paths, the bootstrap clusters at `self.cluster if self.cluster else unit` — i.e., it matches the user's explicit cluster column if set, falling back to unit otherwise (the panel's natural unit of variation). The bootstrap SE overrides the analytical SE for `overall_*` on `n_bootstrap > 0` paths; per-cell `(g, t)` SEs still come from the analytical vcov. - **Note (tail df, 3.9 / M-127):** every non-BM OLS arm's per-cell and aggregated t/p/CI now consume ONE `df_convention`-resolved fallback df — @@ -5056,12 +5060,12 @@ Shipped in `diff_diff/had_pretests.py` as `stute_joint_pretest()` (residuals-in **Requirements checklist (tracks implementation phase completion):** - [x] Phase 1a: Epanechnikov / triangular / uniform kernels with closed-form `κ_k` constants (`diff_diff/local_linear.py`). - [x] Phase 1a: Univariate local-linear regression at a boundary (`local_linear_fit` in `diff_diff/local_linear.py`). -- [x] Phase 1a: HC2 + Bell-McCaffrey DOF correction in `diff_diff/linalg.py` via `vcov_type="hc2_bm"` enum (both one-way and CR2 cluster-robust with Imbens-Kolesar / Pustejovsky-Tipton Satterthwaite DOF). **Weighted Bell-McCaffrey (`hc2_bm + weights`, both one-way and cluster) is now supported on the analytical linalg surface (`compute_robust_vcov` / `solve_ols` / `LinearRegression`) via the clubSandwich WLS-CR2 port.** The port matches clubSandwich's `pweight` (sampling-weight) convention only — `weight_type="aweight"` and `weight_type="fweight"` raise `NotImplementedError` on `hc2_bm + weights` (separate methodology task; CR1 / `vcov_type="hc1"` continues to support all three weight types). Estimator-level `survey_design=` paths continue to use the Taylor-series linearization (TSL) survey variance, which takes precedence over the analytical sandwich (no change in this PR). **Known precision limit on high-leverage FE-dummy coefficients**: the Satterthwaite DOF formula `(tr P)² / sum(P²)` is at the float64 noise floor for contrasts whose projection onto the design is essentially zero (typically FE-dummy coefficients where the dummy fires in a single cluster). In this regime, BLAS reduction-order differences between NumPy and R's BLAS produce 15-30% DOF discrepancies despite vcov matching at machine precision. `_cr2_bm_dof_inner_weighted` detects this regime via two criteria applied union-wise — **(a) batch-relative**: per-contrast max|P| below `1e-10 ×` the largest contrast's max|P| (catches noise-floor coefficients in a per-coefficient sweep where a non-noise reference is available); **(b) absolute single-contrast safe**: per-contrast max|P| below `(EPS × n × k × max(bread_inv_scale, 1))²` (catches single-contrast calls like MPD avg_att where no batch reference exists) — and returns NaN with a `UserWarning` rather than ship arbitrarily-different DOF on the user-facing surface. The coefficient SEs remain valid; only the DOF (and any t-test / CI that depends on it) is suppressed. The diff-diff implementation matches clubSandwich's specific algebra (R source: `CR-adjustments.R::CR2`, `clubSandwich.R::vcov_CR`, `coef_test.R::Satterthwaite_df`), **not** a textbook Pustejovsky-Tipton (2018) §3.3 transform-once derivation. clubSandwich uses W (not √W) in the hat matrix `H_gg = X_g M_U X_g' W_g`, W² in the bias-correction term `S_W = Σ_g X_g' W_g² X_g`, and unweighted residuals in the score construction `s_g = X_g' W_g A_g e_g`. The Satterthwaite DOF uses the full H1/H2/H3 array construction (`get_arrays.R::get_GH`), not the simplified `(tr B)² / tr(B²)` form (which works for the unweighted case but diverges from clubSandwich by 0.5-30% on weighted designs — see `feedback_wls_cr2_clubsandwich_parity`). Pinned at atol=1e-10 against `clubSandwich::vcovCR(..., type="CR2") + coef_test(test="Satterthwaite")$df_Satt` and `Wald_test(test="HTZ")$df_denom` on six weighted scenarios in `benchmarks/data/clubsandwich_cr2_golden.json` (vcov + non-noise-floor per-coefficient DOF + compound-contrast DOF; high-leverage FE-dummy coefficients are suppressed to NaN per the precision-limit note below) (`tests/test_methodology_wls_cr2.py`). Unweighted CR2-BM behavior is algebraically identical to prior (parity within floating-point tolerance after the 2026-07 scores-based evaluation change; regression-safe via `TestUnweightedRegressionStillBitEqual` + `TestDOFFormulaDualPathEquivalence`). The one-way weighted HC2-BM path uses clubSandwich's singleton-cluster CR2 reduction (each obs is its own cluster), routed via `_compute_bm_dof_from_contrasts` when `weights is not None`. clubSandwich pin: `≥ 0.7.0`. +- [x] Phase 1a: HC2 + Bell-McCaffrey DOF correction in `diff_diff/linalg.py` via `vcov_type="hc2_bm"` enum (both one-way and CR2 cluster-robust with Imbens-Kolesar / Pustejovsky-Tipton Satterthwaite DOF). **Weighted Bell-McCaffrey (`hc2_bm + weights`, both one-way and cluster) is now supported on the analytical linalg surface (`compute_robust_vcov` / `solve_ols` / `LinearRegression`) via the clubSandwich WLS-CR2 port.** The port matches clubSandwich's `pweight` (sampling-weight) convention only — `weight_type="aweight"` and `weight_type="fweight"` raise `NotImplementedError` on `hc2_bm + weights` (separate methodology task; CR1 / `vcov_type="hc1"` continues to support all three weight types). Estimator-level `survey_design=` paths continue to use the Taylor-series linearization (TSL) survey variance, which takes precedence over the analytical sandwich (no change in this PR). **Known precision limit on high-leverage FE-dummy coefficients**: the Satterthwaite DOF formula `(tr P)² / sum(P²)` is at the float64 noise floor for contrasts whose projection onto the design is essentially zero (typically FE-dummy coefficients where the dummy fires in a single cluster). In this regime, BLAS reduction-order differences between NumPy and R's BLAS produce 15-30% DOF discrepancies despite vcov matching at machine precision. `_cr2_bm_dof_inner_weighted` detects this regime via two criteria applied union-wise — **(a) batch-relative**: per-contrast max|P| below `1e-10 ×` the largest contrast's max|P| (catches noise-floor coefficients in a per-coefficient sweep where a non-noise reference is available); **(b) absolute single-contrast safe**: per-contrast max|P| below `(EPS × n × k × max(bread_inv_scale, 1))²` (catches single-contrast calls - the removed MPD's post-average `avg_att` was the motivating case - where no batch reference exists) — and returns NaN with a `UserWarning` rather than ship arbitrarily-different DOF on the user-facing surface. The coefficient SEs remain valid; only the DOF (and any t-test / CI that depends on it) is suppressed. The diff-diff implementation matches clubSandwich's specific algebra (R source: `CR-adjustments.R::CR2`, `clubSandwich.R::vcov_CR`, `coef_test.R::Satterthwaite_df`), **not** a textbook Pustejovsky-Tipton (2018) §3.3 transform-once derivation. clubSandwich uses W (not √W) in the hat matrix `H_gg = X_g M_U X_g' W_g`, W² in the bias-correction term `S_W = Σ_g X_g' W_g² X_g`, and unweighted residuals in the score construction `s_g = X_g' W_g A_g e_g`. The Satterthwaite DOF uses the full H1/H2/H3 array construction (`get_arrays.R::get_GH`), not the simplified `(tr B)² / tr(B²)` form (which works for the unweighted case but diverges from clubSandwich by 0.5-30% on weighted designs — see `feedback_wls_cr2_clubsandwich_parity`). Pinned at atol=1e-10 against `clubSandwich::vcovCR(..., type="CR2") + coef_test(test="Satterthwaite")$df_Satt` and `Wald_test(test="HTZ")$df_denom` on six weighted scenarios in `benchmarks/data/clubsandwich_cr2_golden.json` (vcov + non-noise-floor per-coefficient DOF + compound-contrast DOF; high-leverage FE-dummy coefficients are suppressed to NaN per the precision-limit note below) (`tests/test_methodology_wls_cr2.py`). Unweighted CR2-BM behavior is algebraically identical to prior (parity within floating-point tolerance after the 2026-07 scores-based evaluation change; regression-safe via `TestUnweightedRegressionStillBitEqual` + `TestDOFFormulaDualPathEquivalence`). The one-way weighted HC2-BM path uses clubSandwich's singleton-cluster CR2 reduction (each obs is its own cluster), routed via `_compute_bm_dof_from_contrasts` when `weights is not None`. clubSandwich pin: `≥ 0.7.0`. - **Note (scope limitation on absorbed FE):** HC2 and HC2 + Bell-McCaffrey on within-transformed designs still depend on the FULL FE hat matrix because FWL preserves coefficients and residuals but NOT the hat matrix: `h_ii = x_i' (X'X)^{-1} x_i` on the reduced design is not the diagonal of the full FE projection, and CR2's block adjustment `A_g = (I - H_gg)^{-1/2}` likewise depends on the full cluster-block hat matrix. The status across the three estimators that previously rejected this combination: - **`DifferenceInDifferences(absorb=..., vcov_type in {"hc2","hc2_bm"})` — SUPPORTED (auto-route).** When the user pairs `absorb=` with HC2 / HC2-BM, `DiD.fit()` internally promotes the absorb columns to `fixed_effects=` so the existing full-dummy code path computes the algebraically correct vcov from the full FE projection. Verified at ~1e-10 vs `lm() + sandwich::vcovHC(type="HC2")` and `lm() + clubSandwich::vcovCR(cluster=..., type="CR2")` (singleton-cluster CR2 trick for one-way HC2-BM Satterthwaite DOF; PT2018 §3.3 unweighted CR2 algebra). **User-visible surface change**: under the auto-route, the entire `DiDResults` (coefficients, vcov, residuals, fitted_values, r_squared) reflect the full-dummy fit rather than the within-transformed fit — the FE-dummy entries are included in `result.coefficients` / `result.vcov`, `r_squared` is computed on the un-demeaned outcome, and `residuals` / `fitted_values` are on the original scale. `result.att` is unaffected (FWL-equivalent). HC1/CR1 paths on `absorb=` are unchanged (no leverage term). **Survey-design scope**: when `survey_design=` is supplied, the existing survey variance path (Taylor-series linearization / replicate weights) takes precedence over the analytical HC2/HC2-BM sandwich; the auto-route only changes the FE handling (removing the prior reject) and does not redirect to the analytical small-sample sandwich on survey fits. - **`TwoWayFixedEffects(vcov_type in {"hc2","hc2_bm"})` — SUPPORTED (inline full-dummy build).** TWFE has no `absorb=` / `fixed_effects=` parameter (the unit + time FE are baked into the estimator's identity), so the same parameter-swap auto-route used for DiD-absorb / the event-study within spec is not directly applicable to the static mode. Instead, `TwoWayFixedEffects.fit()` bypasses the within-transform when `vcov_type in {"hc2","hc2_bm"}` and builds the full-dummy design `[intercept, treated×post, covariates, factor(unit), factor(time)]` explicitly, then runs OLS through the standard `solve_ols` path so the leverage correction and BM DOF compute on the full FE projection. Verified at atol=1e-10 vs `lm(y ~ treat_post + factor(unit) + factor(post)) + sandwich::vcovHC(type="HC2")` for HC2 and vs `clubSandwich::vcovCR(cluster=seq_len(n), type="CR2")` for the singleton-cluster one-way HC2-BM Satterthwaite DOF; vs `vcovCR(cluster=unit, type="CR2")` for the auto-cluster CR2-BM path. **Auto-cluster default (non-survey analytical path):** TWFE's unit auto-cluster is preserved on `hc2_bm` (routes to CR2-BM at unit) and on `hc2 + wild_bootstrap` (the bootstrap consumes the cluster structure for resampling regardless of the analytical sandwich choice); dropped on explicit `hc2 + analytical` to match the one-way contract (the linalg validator rejects `hc2 + cluster_ids`). `hc2_bm + analytical` with no explicit cluster yields the auto-cluster CR2-BM path. **Survey-design exception:** under `survey_design=` with no explicit `cluster=`, TWFE intentionally keeps the documented implicit-PSU path (the auto-cluster is NOT injected into the survey PSU structure); users who want unit-level PSU injection under a survey design must pass explicit `cluster="unit"` or set `survey_design.psu` directly. **User-visible surface change** (matches the DiD-absorb disclosure above): under HC2 / HC2-BM, `result.coefficients`, `result.vcov`, `result.residuals`, `result.fitted_values`, and `result.r_squared` reflect the full-dummy fit rather than the within-transformed reduced fit (FE-dummy entries are included, `r_squared` is computed on the un-demeaned outcome, residuals/fitted are on the original scale). `result.att`, its SE, and analytical inference are unchanged (FWL-equivalent). HC1 / CR1 / Conley / classical paths remain on the within-transform (no leverage term in those vcov families). **Survey-design scope** (mirrors the DiD-absorb auto-route contract above): when `survey_design=` is supplied, the existing survey variance path (Taylor-series linearization for analytical-weight designs, or replicate-weight variance for BRR/Fay/JK1/JKn/SDR) takes precedence over the analytical HC2/HC2-BM sandwich; the full-dummy build only changes the FE handling (removing the prior reject) and does not redirect to the analytical small-sample sandwich on survey fits. **Replicate-weight survey designs** are blocked at the estimator level: `vcov_type in {"hc2","hc2_bm"}` + replicate weights raises `NotImplementedError` because the replicate refit path re-demeans per replicate, which doesn't compose with the full-dummy build (would require per-replicate full-dummy refit); workaround: use `vcov_type="hc1"` for replicate-weight CR1. `hc2_bm + weights` (Gates 4-5) is now lifted via the clubSandwich WLS-CR2 port; the survey-design rejection here is a separate estimator-level gate (analytical sandwich vs survey TSL precedence), independent of the linalg validator. - **TWFE event-study mode, `spec="within"` + `vcov_type in {"hc2","hc2_bm"}` (the removed `MultiPeriodDiD(absorb=...)` lane) — SUPPORTED (auto-route).** Same auto-route pattern as `DifferenceInDifferences`: the shared core (`DifferenceInDifferences._fit_event_study_core`) promotes the within spec's unit FE from the absorb path to the full-dummy `fixed_effects=` path for HC2 / HC2-BM callers, so the full-dummy code path computes the algebraically correct vcov from the full FE projection on the event-study design (`treated + period_X dummies + treated:period_X interactions + factor(unit)`). Live parity: the within + `hc2` lane is pinned vs `lm() + sandwich::vcovHC(type="HC2")` on the `mpd_absorbed_fe_did` golden (`benchmarks/data/clubsandwich_cr2_golden.json`, the 5-cohort × 5-period event-study fixture the removed class's absorb tests used; `TestTWFEEventStudyAbsorbedFERParity` in `tests/test_estimators_vcov_type.py`), and the clustered `hc2_bm` per-period CR2-BM DOF is pinned vs clubSandwich by CI inversion on `mpd_clustered_avg_att_dof` (`test_twfe_event_study_cluster_hc2_bm_per_period_df_matches_clubsandwich`, same file). The one-way singleton-cluster `clubSandwich::vcovCR(cluster=1:n, type="CR2")` HC2-BM lane, the multi-absorb lane, and the survey-design / replicate-weight auto-route lanes were verified through 3.x on the removed class only (within + `hc2_bm` auto-clusters at unit, so the one-way lane has no within twin); the within spec exposes only the unit FE (multi-column `absorb=` and `fixed_effects=` are not expressible on the surviving surface), so those lanes carry no live ~1e-10 pin. The parity target is a per-period interaction `treated:period_X` because the event-study mode requires the `treated` column to be a time-invariant ever-treated indicator, which lies in the span of the intercept and the post-auto-route unit FE dummies (under `pd.get_dummies(..., drop_first=True)` the dropped reference unit is implicit in the intercept, so the exact alias relation depends on the omitted FE category — it is NOT simply "the sum of treated-cohort unit dummies"). `solve_ols` drops one column from the collinear set under R-style rank-deficiency handling; in the golden fixture (4 ever-treated cohorts of 5 units + 1 never-treated cohort of 5 units) it drops a unit dummy from the never-treated cohort (`unit_25`) and the `treated` main effect remains finite, but the specific column that gets NaN'd is pivot-order and dummy-coding dependent. Either way, the slope coefficients (`treated:period_X`) are identified and invariant to which column was dropped (so was the removed class's post-period-average `avg_att`). The surviving `EventStudyResults` surface exposes only the per-period rows and their interaction-block `vcov` / `vcov_index`, all FWL-invariant, so the removed class's full-dummy surface disclosure (`coefficients` / `vcov` / `residuals` / `fitted_values` / `r_squared` reflecting the full-dummy fit) has no surviving reader. HC1/CR1 paths on the within spec are unchanged (no leverage term). Same survey-design scope as DiD: replicate-weight variance routes through the standard `compute_replicate_vcov` path on the fixed full-dummy design rather than the per-replicate refit branch (which targets the demeaning path); since the auto-routed design does not depend on replicate weights, no refit is needed. **Redundant time-FE skip (history):** when the routed (or directly-supplied) `fixed_effects` list contained the `time` column, the removed class silently skipped emitting `