fix(staggered): CS panel reg/ipw influence functions + per-cell SEs to DRDID parity (#619)

Author: igerber

CHANGELOG.md

@@ -30,6 +30,32 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   `survey_design=SurveyDesign(psu=<cluster_col>)`. No behavior change for unclustered fits.
 
 ### Fixed
+- **`CallawaySantAnna` panel reg/ipw standard errors now match `DRDID` /
+  R `did` exactly** (point estimates unchanged — the omitted terms are mean-zero,
+  which is why ATTs always matched R at ~1e-11 while SEs drifted). Two defects:
+  (1) the covariate-reg influence function omitted `DRDID::reg_did_panel`'s OLS
+  estimation-effect term (`asy.lin.rep.ols %*% M1`) — per-cell SEs were 4-13% and
+  aggregated (simple/group/event-study) SEs 3-20% from R, and **anti-conservative**
+  on some designs (0.958x on the two-period golden fixture) despite the prior
+  documentation calling the plug-in "conservative"; (2) the unweighted covariate-ipw
+  per-cell SE used a weighted *population* variance never scaled by an effective
+  sample size — **~7x inflated** — and its influence function lacked
+  `std_ipw_did_panel`'s propensity-score estimation-effect correction
+  (`asy.lin.rep.ps %*% M2`; aggregated SEs ~2.4% off). The survey ipw branch already
+  carried both corrections (Phase 7a); the fix mirrors it method-uniformly. All panel
+  reg/ipw per-cell SEs (including no-covariate, where ipw's
+  `var_c*(1-p)/(n_c*p)` plug-in equaled R only at treated-share 0.5) are now derived
+  from the same influence function that feeds aggregation: `sqrt(sum(phi^2))`, the
+  convention dr always used. Everything downstream of the per-cell IFs inherits the
+  fix: t-stats/p-values/confidence intervals, aggregated SEs, `event_study_vcov`
+  (HonestDiD input), the bare-`cluster=` CR1 per-cell override, and
+  multiplier-bootstrap SEs. The reg estimation-effect projection is evaluated in the
+  centered basis (`1/n_c + (x-x̄_c)'G^{-1}(x̄_t-x̄_c)`), which is offset-invariant;
+  reg/ipw fits with collinear covariates now fire the same aggregate rank-guard
+  warning as dr. Post-fix parity: per-cell and aggregated reg/ipw vs the R golden
+  fixtures at ~5e-12, and vs fresh R `did` 2.5.1 ipw aggregations at ~1e-10
+  (previously-unasserted golden SE blocks are now enabled, plus R-free numpy
+  reconstructions of both DRDID influence functions).
 - **ImputationDiD/TwoStageDiD covariate fits with zero-weight replicate designs (JK1/plain
   BRR) now produce finite SEs.** Replicate weights that zero out whole PSUs reach Step 1
   unmasked; the previous per-estimator pandas demeaning loops divided 0/0 on zero-total-weight

TODO.md

@@ -32,7 +32,9 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo
 | `ContinuousDiD` CGBS-2024 extensions. (a) `covariates=` kwarg — **DONE (reg/dr)**; remaining: (b) discrete-treatment saturated regression (integer dose currently warned, not routed to per-level coefficients); (c) lowest-dose-as-control per Remark 3.1 when `P(D=0)=0`. Also deferred from the covariate work: `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate — documented `NotImplementedError`), and `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF). | `continuous_did.py` | CGBS-2024 | Heavy | Low |
 | `ImputationDiD` LOO conservative-variance refinement (BJS 2024 Supp. Appendix A.9) — a finite-sample improvement to the auxiliary-model residuals reducing overfit of `tau_tilde_g` to `epsilon`. Asymptotic Theorem-3 variance is implemented and matches R `didimputation` (which also omits LOO by default). | `imputation.py` | imputation-validation | Mid | Low |
 | `TwoWayFixedEffects(vcov_type in {hc2, hc2_bm})` with replicate-weight designs raises `NotImplementedError` (`twfe.py:~233`). The replicate path re-demeans per replicate, which doesn't compose with the full-dummy HC2/HC2-BM build — a correct impl needs per-replicate full-dummy refit. Workaround: `hc1` for replicate-weight CR1. | `twfe.py::fit` | follow-up | Heavy | Low |
-| `CallawaySantAnna` reg-method aggregated SEs sit 3-20% from the R golden fixtures (`two_period` simple/group/dynamic, `dynamic_effects` dynamic) while the ATTs match at 1e-11 and the dr-method aggregated SEs match at ≤7e-6. Pre-existing (byte-identical deltas on the pre-fast-path tree); golden aggregated-SE assertions were therefore enabled for dr scenarios only (`test_golden_simple_aggregation_se`, `test_golden_event_study_aggregation_se`). Investigate the reg-path per-cell IF vs R `DRDID::reg_did_panel`'s. | `staggered.py::_outcome_regression`, `tests/test_csdid_ported.py` | CS-scaling | Mid | Medium |
+| `CallawaySantAnna` DR no-covariate per-cell SE still uses the ddof=1 plug-in `sqrt(var_t/n_t + var_c/n_c)` rather than the IF-based `sqrt(sum(phi^2))` the reg/ipw paths switched to (O(1/n) from R's `drdid_panel`; aggregated SEs unaffected — they already consume the IF). Switch for method-uniform per-cell/aggregated consistency. | `staggered.py::_doubly_robust` (no-cov branch) | CS-scaling | Quick | Low |
+| `CallawaySantAnna` no-covariate ipw treats the propensity as unconditional (no correction term); R's `did` fits an intercept-only logit whose estimation effect is identically zero in the IF, so values match — decide whether to document-only (current REGISTRY note) or mirror R's code path for structural parity. | `staggered.py::_ipw_estimation` (no-cov branch) | CS-scaling | Quick | Low |
+| Fold the R `did` 2.5.1 ipw aggregation yardsticks (hardcoded with provenance in `test_golden_ipw_aggregation_se_vs_r_did_251`) into `csdid_golden_values.json` on the next fixture regeneration — the generator already emits the ipw `aggte` blocks; switch the test to read the JSON. | `benchmarks/R/generate_csdid_test_values.R`, `tests/test_csdid_ported.py` | CS-scaling | Quick | Low |
 | TWFE's HC2/HC2-BM inline full-dummy build (`twfe.py:280-315`) duplicates the dummy-construction logic in `DifferenceInDifferences(fixed_effects=...)` (`estimators.py:478-486`). Extract a shared helper, or delegate TWFE's HC2/HC2-BM path to DiD's `fixed_effects=` branch (with TWFE-specific cluster-default threading), to reduce drift risk on FE naming / survey behavior / result-surface conventions. Substantive refactor — touches both estimators. | `twfe.py::fit`, `estimators.py::DifferenceInDifferences.fit` | follow-up | Heavy | Low |
 
 ### Performance
@@ -65,6 +67,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 | Render `docs/methodology/REPORTING.md` and `REGISTRY.md` as in-site Sphinx pages so cross-refs can use `:doc:` instead of off-site `blob/main` URLs (stable-docs readers can otherwise land on a different revision than their package version). Two paths: (a) add `myst-parser` to `conf.py` + docs extras and link with `:doc:`, or (b) convert both to `.rst`. **Note:** REGISTRY.md is ~4.5k lines of LaTeX-heavy markdown — high risk under the `-W` (warnings-as-errors) Sphinx build; budget multiple rounds. | `docs/conf.py`, `docs/api/business_report.rst`, `docs/api/diagnostic_report.rst`, tutorials 18 & 19 | follow-up | Mid | Low |
 | `ImputationDiD` covariate-path variance lacks a dedicated parity anchor — only the no-covariate staggered panel is R-parity'd, though the covariate path shares the same validated projection code. Add a small dense-design **hand-calc** for the covariate projection (no external tooling), or a covariate (time-varying X) R `didimputation` golden asserting overall/ES SE parity (the golden variant needs local R). | `tests/test_methodology_imputation.py`, `benchmarks/R/generate_didimputation_golden.R` | imputation-validation | Mid | Low |
 | Add true half-sample BRR replicate-weight regressions per estimator family (current tests use Fay-like 0.5/1.5 perturbations; `test_survey_phase6.py` covers true BRR at the helper level). | `tests/test_replicate_weight_expansion.py` | #253 | Mid | Low |
+| `test_bootstrap_se_matches_pre_pr4_baseline` (dCDH bit-identity guard) fails under FULL-SUITE order but passes standalone/as-a-file: its `pure_python` backend detection reads a state some earlier test perturbs, so it compares against the wrong per-backend baseline arm (observed failure message says `backend=pure-python` while the run produces the other arm's value). Reproduced identically on pristine `origin/main` (2026-07-05, `b56931a6`), so it pre-dates the CS reg/ipw IF fix. Make the backend detection order-robust (e.g. re-resolve from `diff_diff._backend` at call time, or isolate via monkeypatch). | `tests/test_survey_dcdh.py::TestBootstrapCellPeriod` | CS-scaling | Quick | Low |
 | Port the CI `<notebook-prose>` extraction into the reviewer-eval harness so `docs/tutorials/*.ipynb` cases (currently guarded out of `verify-corpus`/`run`) can be reviewed with CI-equivalent context. | `tools/reviewer-eval/adapters/ci_prompt.py` | local-review | Mid | Low |
 
 ---

benchmarks/R/generate_csdid_test_values.R

@@ -108,6 +108,13 @@ data2_ipw <- build_sim_dataset(sp2_ipw)
 res2_ipw <- att_gt(yname = "Y", xformla = ~X, data = data2_ipw, tname = "period",
                    idname = "id", gname = "G", est_method = "ipw",
                    bstrap = FALSE, cband = FALSE)
+# Aggregations for the ipw scenario. Until the fixture is regenerated, these
+# values are hardcoded (from did 2.5.1 / DRDID 1.3.0, 2026-07-05) in
+# tests/test_csdid_ported.py::test_golden_ipw_aggregation_se_vs_r_did_251;
+# on the next regeneration the test can switch to reading them from the JSON.
+agg2_ipw_simple <- aggte(res2_ipw, type = "simple", bstrap = FALSE, cband = FALSE)
+agg2_ipw_dynamic <- aggte(res2_ipw, type = "dynamic", bstrap = FALSE, cband = FALSE)
+agg2_ipw_group <- aggte(res2_ipw, type = "group", bstrap = FALSE, cband = FALSE)
 
 scenarios$with_covariates_dr <- list(
   data = export_data(data2),
@@ -125,7 +132,12 @@ scenarios$with_covariates_ipw <- list(
   data = export_data(data2_ipw),
   params = list(est_method = "ipw", control_group = "nevertreated",
                 xformla = "~X", base_period = "varying"),
-  results = list(group_time = extract_gt(res2_ipw))
+  results = list(
+    group_time = extract_gt(res2_ipw),
+    simple = list(att = agg2_ipw_simple$overall.att, se = agg2_ipw_simple$overall.se),
+    group = extract_agg(agg2_ipw_group),
+    dynamic = extract_agg(agg2_ipw_dynamic)
+  )
 )
 
 # ---------------------------------------------------------------------------