CS positional base period + absorbed-FE full-K SE parity vs R

CHANGELOG.md

@@ -40,6 +40,41 @@ 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). For
+  `base_period="universal"`, each cohort's positional base is now materialized as a zero reference
+  cell (`att=0`, `se=NaN`) in `group_time_effects` / `to_dataframe("group_time")` at its positional
+  base event time (`e = base - g`, which can be `-2`, `-3`, … on gapped grids), matching R's
+  `att_gt` table and `aggte(type="dynamic")` — including the overlapping-reference case where the
+  zero base dilutes another cohort's estimated pre-trend at the same event time (analytical AND
+  multiplier-bootstrap paths verified vs `did` 2.5.1 to ~1e-5). The `group` / `simple` aggregations
+  use post-treatment cells only and exclude the reference.
+- **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

@@ -111,6 +111,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 |
@@ -128,6 +129,8 @@ Doable in principle, but no current caller and/or explicitly out of paper scope.
 
 | Issue | Location | PR | Priority |
 |-------|----------|----|----------|
+| CallawaySantAnna **unbalanced-panel event-study weighting**: a dynamic horizon pooling >1 (g,t) cell weights each cell by its per-cell aggregation weight (valid `n_treated` / `agg_weight`), but R `did::aggte` weights by the fixed cohort probability `pg = n_g/N` (or survey mass) from the group column. Balanced panels coincide exactly (verified ~1e-5, incl. universal zero-reference dilution); unbalanced panels can differ at multi-cell horizons (e.g. a two-real-cell horizon `-0.065` vs R `-0.136`). Pre-existing + general (independent of universal reference cells; documented Deviation from R in REGISTRY). Fix = weight every event-study cell by fixed cohort mass (touches `_aggregate_event_study` + the bootstrap bucket + WIF `pg` consistency); gate on balanced byte-identity + fresh unbalanced R parity. | `staggered_aggregation.py`, `staggered_bootstrap.py` | SE-audit D3 | Low |
+| CallawaySantAnna event-study bucket/weight construction is duplicated between the analytical aggregator (`staggered_aggregation.py::_aggregate_event_study`) and the multiplier bootstrap (`staggered_bootstrap.py`): both group (g,t) by `e = t - g`, apply the finite/NaN/reference masks, and read cohort weights. Both already consume the same source-materialized universal reference cells (so they agree), but the bucket logic is copy-pasted. Extract one shared helper returning per-event-time buckets (finite cells, NaN cells, reference flags, cohort weights, combined-IF inputs) used by both. Pure refactor; gate on byte-identical analytical + bootstrap output. | `staggered_aggregation.py`, `staggered_bootstrap.py` | SE-audit D3 | Low |
 | `StackedDiD` survey re-resolution intra-file dedup (raw-weight extraction ×3, compose-normalize ×3, resolve-on-stacked ×2). The cross-estimator ContinuousDiD/EfficientDiD panel-to-unit collapse consolidation LANDED (#226 shared helpers `ResolvedSurveyDesign.subset_to_units_by_row_idx` / `build_unit_first_row_index`); StackedDiD is deliberately NOT on that path (control units are duplicated across sub-experiments, so it re-resolves at stacked granularity rather than collapsing to one row per unit). The residual is stacked-specific, low value, and touches the numerically-sensitive composed-weight renormalization. Post-filter re-resolution / metadata-recompute unification across the three estimators was assessed and is not warranted — they use genuinely different mechanisms and already delegate to shared `_resolve_survey_for_fit` / `compute_survey_metadata`. | `stacked_did.py` | #226 | Low |
 | `solve_ols` residual memory floor is inherent to SVD-based lstsq after the PR-E marshalling slim (rust-side allocator high-water 15.32 -> 7.81 GB on the 2.4M x 130 clustered solve; remaining = one fused equilibrated input copy + thin-SVD U transient + vcov scores block + LAPACK gelsd internals on the python path). Further reduction requires the tall-skinny-QR algorithm change - same trade-off as the parked QR-reuse row below (library-wide last-digit perturbation + loses the second independent collinearity detector); the two land together if ever reopened. Diagnosis + measurements in docs/performance-plan.md. | `rust/src/linalg.rs::solve_ols`, `diff_diff/linalg.py` | PR-E | Low |
 | SunAbraham-family fits are SOLVER-bound after the v3.6.x demeaning speedups: on the county-class shape (3.1k units x 60 months, ~100 interaction columns) `solve_ols` = ~60% of fit time (faer SVD 1.24s + pivoted-QR rank detection 0.70s of 3.27s), and the share grows further under the Rust demean kernel; the pivoted QR alone is ~20% of the whole fit (measured attribution in docs/performance-plan.md). Candidate fix is QR-reuse (one factorization for rank detection + solve). **Parked per the 2026-07 correctness-first decision**: the change perturbs every coefficient in the library at the last digits (golden/parity/bit-identity recapture) and removes the SVD's independent singular-value truncation - a second, different near-collinearity detector that is deliberate defense-in-depth for the no-silent-failures contract (the FE-span guard episode showed borderline cases are real). Re-open only on practitioner demand for wide event studies; any implementation gates on fixest/R end-to-end parity, not just internal consistency. | `linalg.py::solve_ols`, `linalg.py::_detect_rank_deficiency` | PR-C attribution | Medium |

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)