fix(staggered,inference): CS positional base period + absorbed-FE full-K SE parity vs R

Two standard-error parity fixes vs canonical R (SE-audit D3, D4).

D3 - CallawaySantAnna base period is now selected positionally (nearest observed
period) via shared `_select_base_period` / `_valid_periods_for_group` helpers used
across every estimation path (panel fast/vectorized/covariate-reg, repeated
cross-sections, and the universal filters), matching R `did::att_gt` on gapped
(non-consecutive) grids. The pre/post split is corrected to `t < g` (independent of
anticipation; only the post/universal base uses anticipation), verified against a
deparse of `did` 2.5.1 `compute.att_gt`. Consecutive grids are byte-identical to the
old rule; gapped panels now reproduce every R cell (fewer_periods {1,3,4,6},
non_consecutive {1,2,5,7}) instead of NaN-ing cells whose calendar base was
unobserved, and universal no longer materializes the positional base as a fake zero
cell. Aggregated event-study / group SEs on gapped panels are corrected accordingly.

D4 - absorbed-FE (`absorb=`) non-clustered classical/hc1 standard errors now use the
full-K finite-sample scale (`linalg._absorbed_fe_vcov_scale`, a single scalar rescale
of the k_visible vcov; fail-closed to NaN inference when the full-K residual dof is
non-positive), matching `fixest::feols(vcov="iid"/"hetero")` and the explicit
full-dummy oracle. Applies to TwoWayFixedEffects(classical),
DifferenceInDifferences(absorb=), and MultiPeriodDiD(absorb=). Clustered SEs are
unchanged (gated on cluster_ids is None): the k_visible CR1 scale matches fixest for
FE nested in the cluster; non-nested absorbed FE (e.g. time FE with unit clustering)
are a small documented deviation tracked in TODO.md. hc2/hc2_bm (Satterthwaite DOF),
survey vcov, full-dummy fits, and SunAbraham hc1 (auto-clusters) are unaffected.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0178r8ZK6VRiXbaBiknrWjQB

Author: igerber

CHANGELOG.md

@@ -40,6 +40,34 @@ 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` base period is now selected positionally (sorted-index), matching R
+  `did::att_gt` on gapped panels.** Previously the base period used literal calendar arithmetic
+  (`t-1` pre-treatment, `g-1-anticipation` post/universal), so on non-consecutive period grids
+  (e.g. biennial surveys, skipped years) cells whose calendar base was unobserved were NaN'd as
+  `missing_period` even though R estimates them from the nearest observed period — which also
+  corrupted the aggregated event-study / group SEs. The base is now the nearest observed period
+  (largest observed `p < t` pre-treatment; largest observed `p` with `p + anticipation < g`
+  post/universal), and the pre/post split is on `t < g` (independent of anticipation, matching a
+  deparse of `did` 2.5.1 `compute.att_gt`) — resolving the internal inconsistency with the
+  library's own dCDH estimator, which already used positional neighbors. On consecutive grids
+  this is byte-identical to the old rule (all existing goldens unchanged); on gapped panels we now
+  reproduce every R cell (e.g. `fewer_periods` {1,3,4,6} 7/15 previously-NaN cells, `reg` to
+  ~1e-11 / `dr` to ~1e-4). A single shared `_select_base_period` helper is used across all
+  estimation paths (panel fast / vectorized / covariate-reg, and repeated cross-sections).
+- **Absorbed-FE (`absorb=`) non-clustered classical/hetero standard errors now use the full-K
+  finite-sample scale, matching `fixest`.** The within-transform variance scale (`sse/(n-k)` for
+  `classical`, `n/(n-k)` for `hc1`) counted only the *visible* regressors `k`, excluding the
+  absorbed fixed effects, so `TwoWayFixedEffects(vcov_type="classical")`,
+  `DifferenceInDifferences(absorb=..., vcov_type in {classical,hc1})`, and
+  `MultiPeriodDiD(absorb=..., vcov_type in {classical,hc1})` non-clustered SEs sat ~6.5% below
+  `fixest::feols(vcov="iid"/"hetero")` — even though the reported t-`df` already counted the
+  absorbed FE (an internal inconsistency). A single scalar rescale
+  (`(n-k)/(n-k-df_adjustment)`, fail-closed when the full-K residual dof is non-positive) now
+  aligns the SE's `k` with `K_full`, so the absorb path equals the explicit full-dummy
+  (`fixed_effects=`) path and fixest. **Clustered** SEs are unchanged (fixest's nested-FE `ssc`
+  convention already matches with `k_visible`); `hc2`/`hc2_bm` (leverage / Satterthwaite DOF) and
+  survey vcov are unaffected; `SunAbraham` `hc1` auto-clusters at unit and so keeps its documented
+  deviation.
 - **Non-finite degrees of freedom now fail closed to all-NaN inference.**
   `safe_inference()` previously rejected `df <= 0` but let a non-finite `df` (NaN) through,
   producing an inconsistent tuple (finite t-stat, NaN p-value/CI). It now returns all-NaN for

TODO.md

@@ -110,6 +110,7 @@ exists but parity can't be verified without a local toolchain.
 | Extend `WooldridgeDiD` `method ∈ {logit, poisson}` with `vcov_type ∈ {classical, hc2, hc2_bm}`: composing HC2 leverage + Bell-McCaffrey DOF with the QMLE pseudo-residual sandwich needs derivation + R parity vs `clubSandwich::vcovCR(glm, type="CR2")`. Rejected at `__init__`. | `wooldridge.py` | follow-up | Medium |
 | `PreTrendsPower` CS/SA `anticipation=1` R-parity fixture: R `pretrends` has no anticipation parameter, so the Python `_extract_pre_period_params` anticipation filter isn't R-parity-locked. Build a synthetic CS/SA result with `anticipation=1` and assert γ_p matches R's `slope_for_power()`. (Mechanism already covered by MC + full-VCV tests.) | `tests/test_methodology_pretrends.py`, `generate_pretrends_golden.R` | PR-C | Low |
 | Harmonize SunAbraham's HC1 within-transform finite-sample correction with `fixest::sunab()` — SA applies `n/(n-k_dm)`, fixest applies `n/(n-k_total)` (counts absorbed FE); ~1-2% SE difference, documented as a "Deviation from R" and pinned at `atol=5e-3`. Either thread `df_adjustment` or keep as an intentional, R-verified difference. | `sun_abraham.py`, `linalg.py` | follow-up | Low |
+| Absorbed-FE **clustered** CR1 with *non-nested* FE: for `absorb=[FE1,FE2], cluster=FE1` (e.g. `absorb=["unit","time"], cluster="unit"`), `fixest` counts the non-nested FE (time) in the CR1 `(n-1)/(n-k)` finite-sample denominator, but the clustered path uses only `k_visible`. D4 harmonized the *non-clustered* classical/hc1 full-K scale (`_absorbed_fe_vcov_scale`) and left the clustered path unchanged — correct for FE nested in the cluster, a small deviation for non-nested FE (documented in REGISTRY within-transform note). Thread a non-nested `df_adjustment` into the clustered CR1 factor; verify vs `fixest::feols(..., cluster=)`. | `linalg.py`, `estimators.py` | SE-audit D4 | Low |
 | Rust multiplier-bootstrap weight RNG (`generate_bootstrap_weights_batch`) seeds `Xoshiro256PlusPlus::seed_from_u64(seed+i)` per row; audit Python callers (`sdid.py`, `efficient_did_bootstrap.py`, `bootstrap_utils.py`) for parity-test gaps and, where a numpy-canonical equivalent exists, pre-generate in Python and pass through PyO3 (same fix shape as TROP RNG parity #354). | `rust/src/bootstrap.rs`, `bootstrap_utils.py` | follow-up | Medium |
 | `SyntheticDiD` bootstrap cross-language parity anchor vs R `synthdid::vcov(method="bootstrap")` or Julia `Synthdid.jl` (refit-native). Same-library validation is in place; Julia is the cleanest target. Tolerance ~1e-6 (BLAS+RNG paths preclude 1e-10). | `benchmarks/R/`, `benchmarks/julia/`, `tests/` | follow-up | Low |
 | CS R helpers hard-code `xformla = ~1`; no covariate-adjusted R benchmark for the IRLS path. | `tests/test_methodology_callaway.py` | #202 | Low |

diff_diff/estimators.py

@@ -21,6 +21,7 @@
 
 from diff_diff.linalg import (
     LinearRegression,
+    _absorbed_fe_vcov_scale,
     _expand_vcov_with_nan,
     compute_r_squared,
     solve_ols,
@@ -1991,6 +1992,29 @@ def _refit_mp_absorb(w_r):
         if survey_weights is not None and survey_weight_type == "fweight":
             n_eff_df = int(round(np.sum(survey_weights)))
         df = n_eff_df - k_effective - n_absorbed_effects
+
+        # Absorbed-FE variance scale (fixest full-K convention): the within-
+        # transform solve_ols above scales the non-clustered classical/hc1 vcov
+        # by k_visible, but the correct finite-sample count is
+        # K_full = k_effective + n_absorbed_effects (matching `df` just above and
+        # fixest feols(vcov="iid"/"hetero")). Rescale so the SE's k agrees with
+        # the t-df's. Gated exactly as LinearRegression.fit: clustered SEs keep
+        # k_visible (fixest ssc nested-FE convention), hc2/hc2_bm use
+        # leverage/Satterthwaite DOF, survey has its own df. When the full-K
+        # residual dof is non-positive the helper returns NaN and we void the
+        # vcov -> NaN inference (fail-closed, per the non-finite-df contract).
+        if (
+            n_absorbed_effects > 0
+            and effective_cluster_ids is None
+            and not _use_survey_vcov
+            and _fit_vcov_type in ("classical", "hc1")
+        ):
+            _fe_scale_mp = _absorbed_fe_vcov_scale(n_eff_df, k_effective, n_absorbed_effects)
+            if np.isnan(_fe_scale_mp):
+                vcov = np.full_like(vcov, np.nan)
+            elif _fe_scale_mp != 1.0:
+                vcov = vcov * _fe_scale_mp
+
         if resolved_survey is not None and resolved_survey.df_survey is not None:
             df = resolved_survey.df_survey
         # Replicate df: rank-deficient → NaN inference; dropped replicates → n_valid-1

diff_diff/linalg.py

@@ -267,6 +267,47 @@ def _expand_coefficients_with_nan(
     return coef_full
 
 
+def _absorbed_fe_vcov_scale(n_eff: float, k_eff: int, df_adjustment: int) -> float:
+    """Finite-sample vcov rescale for absorbed fixed effects (fixest full-K).
+
+    Within-transform (``absorb=``) fits residualize the FE out of the design, so
+    the classical ``sse/(n-k)`` and HC1 ``n/(n-k)`` variance factors use only the
+    *visible* regressor count ``k_eff = k_visible``. fixest (and the reported
+    t-``df``) count the absorbed FE too, i.e. ``K_full = k_visible +
+    df_adjustment``. This returns the scalar that maps the ``k_visible`` vcov to
+    the ``K_full`` one::
+
+        vcov_full = vcov_visible * (n_eff - k_eff) / (n_eff - k_eff - df_adjustment)
+
+    which is algebraically exact for both the classical and HC1 factors (both are
+    pure scalars in ``1/(n-k)``; the single-coefficient robust variance is
+    FWL-invariant between the demeaned and full-dummy designs, so only the
+    ``1/(n-k)`` scalar differs).
+
+    Returns:
+    - ``1.0`` when ``df_adjustment <= 0`` (no absorbed FE -- a no-op).
+    - the finite scale when the full-K residual dof
+      ``(n_eff - k_eff - df_adjustment)`` is positive.
+    - ``nan`` (fail-closed) when that full-K residual dof (or the visible dof) is
+      non-positive: the full-K variance is undefined for such a saturated
+      within-transform design, so callers void the vcov to NaN -> NaN inference
+      (per the non-finite-df fail-closed contract) rather than leaving a
+      misleading ``k_visible`` SE in place.
+
+    Callers must gate on non-clustered ``classical``/``hc1``: clustered SEs
+    follow fixest's ``ssc`` nested-FE convention (FE nested in the cluster are
+    not counted, so ``k_visible`` already matches for the nested case) and
+    ``hc2``/``hc2_bm`` use leverage / Satterthwaite DOF -- none must be rescaled.
+    """
+    denom_visible = n_eff - k_eff
+    denom_full = n_eff - k_eff - df_adjustment
+    if df_adjustment <= 0:
+        return 1.0
+    if denom_full <= 0 or denom_visible <= 0:
+        return float("nan")
+    return denom_visible / denom_full
+
+
 def _expand_vcov_with_nan(
     vcov_reduced: np.ndarray,
     k_full: int,
@@ -3897,6 +3938,29 @@ def fit(
                 n_eff_df = int(round(np.sum(_fit_weights)))
             elif np.any(_fit_weights == 0):
                 n_eff_df = int(np.count_nonzero(_fit_weights > 0))
+        # Absorbed-FE variance scale (fixest full-K convention): the classical
+        # sse/(n-k) and HC1 n/(n-k) factors computed above use k_visible, but
+        # with absorbed FE the correct finite-sample count is
+        # K_full = n_params_effective_ + df_adjustment (the t-df below already
+        # uses it). Rescale the NON-CLUSTERED iid/hetero vcov so the SE's k
+        # agrees with the t-df's and with fixest feols(vcov="iid"/"hetero").
+        # Clustered SEs keep k_visible (fixest ssc nested-FE convention already
+        # matches); hc2/hc2_bm use leverage/Satterthwaite DOF; survey has its
+        # own df; full-dummy fits carry df_adjustment == 0. When the full-K
+        # residual dof is non-positive the helper returns NaN and we void the
+        # vcov -> NaN inference (fail-closed, per the non-finite-df contract).
+        if (
+            df_adjustment > 0
+            and effective_cluster_ids is None
+            and not _use_survey_vcov
+            and _fit_vcov_type in ("classical", "hc1")
+        ):
+            _fe_scale = _absorbed_fe_vcov_scale(n_eff_df, self.n_params_effective_, df_adjustment)
+            if np.isnan(_fe_scale):
+                self.vcov_ = np.full_like(self.vcov_, np.nan)
+            elif _fe_scale != 1.0:
+                self.vcov_ = self.vcov_ * _fe_scale
+
         self.df_ = n_eff_df - self.n_params_effective_ - df_adjustment
 
         # Survey degrees of freedom: n_PSU - n_strata (overrides standard df)