fix(inference): three standard-error parity fixes vs canonical R

CHANGELOG.md

@@ -40,6 +40,44 @@ 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
+- **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
+  any non-finite `df`. `MultiPeriodDiD(vcov_type="hc2_bm")` likewise no longer falls back to
+  the shared residual `df` when a coefficient's Bell-McCaffrey Satterthwaite DOF is non-finite
+  (guard-suppressed) — such a coefficient's inference is now NaN rather than silently computed
+  from a different `df`, preserving the joint-NaN inference contract. Only affects coefficients
+  whose BM DOF was declared unreliable (high-leverage / collinear); well-conditioned
+  treatment / event-study / average contrasts are unchanged.
+- **Unweighted clustered CR2 / Bell-McCaffrey per-coefficient Satterthwaite DOF no
+  longer returns non-physical values** for high-leverage FE-dummy / collinear nuisance
+  columns. The simple `(tr B)²/tr(B²)` form could produce a garbage DOF there (float64
+  noise in `trace_B2` → up to ~1e61, or a finite-but-inflated value above the cluster
+  count). The unweighted path now carries the same noise-floor guard as the weighted
+  path plus a `DOF ≤ G` (cluster-count) physical bound, NaN-ing those columns with a
+  warning. The treatment / event-study / compound-average contrasts that estimators
+  actually consume are unchanged and match R clubSandwich; this only affects direct
+  `LinearRegression(vcov_type="hc2_bm", cluster_ids=...)` callers reading full
+  per-coefficient DOF vectors.
+- **`ChaisemartinDHaultfoeuille` phase-1 placebo (`DID_M^pl`) point estimate sign
+  corrected.** The single-horizon (`L_max=None`) `placebo_effect` used the opposite
+  (forward) difference order from the multi-horizon placebo path and R
+  `did_multiplegt_dyn`, so it was sign-flipped on pure-direction panels (magnitude was
+  bit-identical). It now uses the backward-difference × switch-direction convention,
+  matching R to working precision on `joiners_only` / `leavers_only`. Mixed-direction
+  panels retain the separately-documented period-vs-cohort control-set deviation. SEs
+  (NaN by design for the single-period placebo) and the multi-horizon placebo path are
+  unaffected.
+- **`HonestDiD(method="smoothness")` now returns a finite FLCI when the estimated
+  identified set is empty**, instead of silently NaN-propagating all inference. When the
+  observed pre-trend's curvature exceeds `M`, the `δ_pre = β_pre`-pinned identified-set LP
+  is infeasible (`lb`/`ub` = NaN) — but the optimal FLCI does not depend on that LP (its
+  worst-case bias is taken over `δ ∈ Δ^SD(M)` treating β as random), so it is well-defined
+  given `(Σ, M)`. R's `HonestDiD::createSensitivityResults` returns the FLCI in exactly this
+  case; previously `fit()` returned all-NaN, so smoothness sensitivity analysis yielded no
+  inference on any event study with non-trivial pre-trend curvature. The FLCI now matches R
+  to ~1e-3 at `M=0`; a known optimizer/center divergence at intermediate `M` (up to ~9% on
+  wide pre/post windows, CI *width* unaffected) is documented in `REGISTRY.md` and deferred.
 - **`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:

TODO.md

@@ -36,6 +36,7 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo
 | `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 |
+| `HonestDiD` Δ^SD optimal-FLCI affine-estimator optimizer (Nelder-Mead over slope weights) diverges from R `.FLCI.computeFLCI` at intermediate small M (~`[0.005, 0.02]`): the CI **center** shifts up to ~9% on wide pre/post windows while the **width** matches to ~1e-3. Not a local-minimum artifact (multi-start does not move it); the two implementations land on different affine estimators of near-equal length. Coverage is unaffected. Reconcile our FLCI optimization with R's exact algorithm. Exposed when the identified-set NaN gate was removed (finite CIs now surface at all M). | `honest_did.py::_compute_optimal_flci` | SE-audit | Mid | Low |
 
 ### Performance
 

diff_diff/chaisemartin_dhaultfoeuille.py

@@ -4665,7 +4665,12 @@ def _compute_placebo(
         else:
             joiner_avg = float((y_prev[joiner_mask] - y_pre_prev[joiner_mask]).mean())
             stable0_avg = float((y_prev[stable0_mask] - y_pre_prev[stable0_mask]).mean())
-            placebo_plus_t = joiner_avg - stable0_avg
+            # Backward-difference convention (pre-period minus reference,
+            # Y_{t-2} - Y_{t-1}), matching the multi-horizon placebo path
+            # `_compute_multi_horizon_placebos` (`switcher_change - ctrl_avg`
+            # with `Y_bwd - Y_ref`, times switch_direction=+1 for joiners) and
+            # R `did_multiplegt_dyn`. Equivalently `stable0_avg - joiner_avg`.
+            placebo_plus_t = stable0_avg - joiner_avg
 
         # Leavers side: symmetric A11 distinction
         if n_01 == 0:
@@ -4678,7 +4683,10 @@ def _compute_placebo(
         else:
             stable1_avg = float((y_prev[stable1_mask] - y_pre_prev[stable1_mask]).mean())
             leaver_avg = float((y_prev[leaver_mask] - y_pre_prev[leaver_mask]).mean())
-            placebo_minus_t = stable1_avg - leaver_avg
+            # Backward-difference convention times switch_direction=-1 for
+            # leavers (see the joiners-side note above): the multi-horizon path
+            # contributes `leaver_avg - stable1_avg` here.
+            placebo_minus_t = leaver_avg - stable1_avg
 
         placebo_plus_per_t.append(placebo_plus_t)
         placebo_minus_per_t.append(placebo_minus_t)

diff_diff/estimators.py

@@ -2132,11 +2132,16 @@ def _refit_mp_absorb(w_r):
             idx = interaction_indices[period]
             effect = coefficients[idx]
             se = np.sqrt(vcov[idx, idx])
-            # Prefer per-coefficient BM DOF when available (hc2_bm path);
-            # otherwise fall back to the shared analytical df.
-            period_df = df
-            if _bm_dof_per_coef is not None and np.isfinite(_bm_dof_per_coef[idx]):
+            # Use the per-coefficient BM DOF when available (hc2_bm path);
+            # otherwise fall back to the shared analytical df. On the hc2_bm path
+            # the coefficient's OWN Bell-McCaffrey DOF governs its inference: a
+            # non-finite (guard-suppressed, unreliable) BM DOF makes t/p/CI
+            # undefined, so pass it through to `safe_inference` (which returns
+            # all-NaN) rather than silently falling back to the residual df.
+            if _bm_dof_per_coef is not None:
                 period_df = float(_bm_dof_per_coef[idx])
+            else:
+                period_df = df
             t_stat, p_value, conf_int = safe_inference(effect, se, alpha=self.alpha, df=period_df)
 
             period_effects[period] = PeriodEffect(

diff_diff/honest_did.py

@@ -2466,14 +2466,20 @@ def _compute_smoothness_bounds(
         # Construct constraints
         A_ineq, b_ineq = _construct_constraints_sd(num_pre, num_post, M)
 
-        # Solve for identified set bounds with delta_pre = beta_pre pinned
+        # Solve for the identified set bounds with delta_pre = beta_pre pinned.
+        # When the observed pre-trend's own curvature exceeds M this LP is
+        # infeasible and the ESTIMATED identified set is empty (lb/ub = NaN) - the
+        # point estimate rejects Delta^SD(M). That does NOT invalidate the FLCI: the
+        # optimal FLCI is an affine estimator whose worst-case bias is taken over
+        # delta in Delta^SD(M) treating beta as random, so it is well-defined given
+        # (sigma, M) regardless of whether the realized beta_pre lies in Delta. R's
+        # HonestDiD::createSensitivityResults returns the FLCI in exactly this case,
+        # so we compute and return it (leaving lb/ub = NaN to flag the empty
+        # estimated id-set) rather than NaN-propagating the whole result.
         lb, ub = _solve_bounds_lp(beta_pre, beta_post, l_vec, A_ineq, b_ineq, num_pre)
 
-        # Propagate infeasibility: if bounds are NaN, CI is NaN too
-        if np.isnan(lb) or np.isnan(ub):
-            return np.nan, np.nan, np.nan, np.nan
-
-        # Compute optimal FLCI (Rambachan & Roth Section 4.1)
+        # Compute optimal FLCI (Rambachan & Roth Section 4.1) - independent of the
+        # identified-set LP above.
         if sigma_full.shape[0] == num_pre + num_post:
             ci_lb, ci_ub = _compute_optimal_flci(
                 beta_pre,
@@ -2487,7 +2493,12 @@ def _compute_smoothness_bounds(
                 df=df,
             )
         else:
-            # Fallback to naive FLCI when full sigma unavailable
+            # The naive fallback FLCI extends the identified set by z*se, so it
+            # genuinely needs finite id-set bounds; when the LP was infeasible the
+            # naive CI is undefined (this branch is only reachable without the full
+            # covariance matrix).
+            if np.isnan(lb) or np.isnan(ub):
+                return lb, ub, np.nan, np.nan
             se = np.sqrt(l_vec @ sigma_post @ l_vec)
             ci_lb, ci_ub = _compute_flci(lb, ub, se, self.alpha, df=df)
 

diff_diff/linalg.py

@@ -2055,10 +2055,14 @@ def _cr2_bm_dof_inner(
     omega_all = {g: A_g_Xbi[g] @ contrasts for g in unique_clusters}
 
     dof_vec = np.empty(m)
+    # Retain max|B_{g,h}| per contrast so we can NaN-guard noise-floor
+    # degeneracies in a second pass (mirrors `_cr2_bm_dof_inner_weighted`).
+    max_abs_B_arr = np.zeros(m)
     for j in range(m):
         q = Q[:, j]
         trace_B = float(np.sum(q * q))
         trace_B2 = 0.0
+        max_abs_B = 0.0
         omega_cache = {g: omega_all[g][:, j] for g in unique_clusters}
         for g in unique_clusters:
             idx_g = cluster_idx[g]
@@ -2069,8 +2073,85 @@ def _cr2_bm_dof_inner(
                 M_gh = M[np.ix_(idx_g, idx_h)]
                 val = float(omega_g @ M_gh @ omega_h)
                 trace_B2 += val * val
+                if abs(val) > max_abs_B:
+                    max_abs_B = abs(val)
+        max_abs_B_arr[j] = max_abs_B
         dof_vec[j] = (trace_B * trace_B) / trace_B2 if trace_B2 > 0 else np.nan
 
+    # Noise-floor NaN-guard (unweighted analogue of the guard in
+    # `_cr2_bm_dof_inner_weighted`). For a high-leverage FE-dummy / collinear
+    # nuisance column, `trace_B2 = sum_{g,h} B_{g,h}²` collapses to float64
+    # accumulation noise while `trace_B` stays O(1), so `(trace_B)²/trace_B2`
+    # blows up to a non-physical DOF (observed up to ~1e61 on the absorbed-FE
+    # golden). `trace_B2 > 0` alone does not catch this because the roundoff is
+    # positive. The Satterthwaite DOF is scale-invariant, so a contrast whose
+    # `max|B_{g,h}|` sits at the accumulation floor is unreliable and
+    # BLAS-order-dependent; we return NaN with a warning rather than ship it.
+    # Two union-wise criteria:
+    #   1. Batch-relative: a contrast's max|B| is 1e-10× below the largest
+    #      contrast's, i.e. it sits at the accumulation floor relative to a
+    #      well-conditioned column (per-coefficient sweeps, `contrasts=eye(k)`).
+    #      `B_{g,h}` scales as ‖c‖² while the Satterthwaite DOF is scale-invariant,
+    #      so the comparison is done on `max|B| / ‖c‖²` - otherwise the same
+    #      contrast passed at two scales in one batch would spuriously flag the
+    #      smaller copy (‖c‖² differs by the scale²). For `contrasts=eye(k)` every
+    #      `‖c‖²=1`, so this is a no-op on the per-coefficient path.
+    #   2. Absolute (single-contrast safe): max|B| < (EPS·n·k·bread_scale)².
+    #      Calibrated for the O(1)-scale contrasts estimators pass (per-coef unit
+    #      vectors, the averaging-weight compound contrast).
+    # The treatment / event-study / compound-average contrasts that estimators
+    # actually consume are well-conditioned and unaffected.
+    _EPS = np.finfo(np.float64).eps
+    n_obs, k_X = X.shape
+    bread_scale = float(np.max(np.abs(bread_inv))) if bread_inv.size else 1.0
+    abs_noise_floor = (_EPS * n_obs * k_X * max(bread_scale, 1.0)) ** 2
+    abs_degenerate = max_abs_B_arr < abs_noise_floor
+    if m > 1:
+        contrast_sq_norm = np.einsum("ij,ij->j", contrasts, contrasts)
+        contrast_sq_norm = np.where(contrast_sq_norm > 0, contrast_sq_norm, 1.0)
+        scaled_max_B = max_abs_B_arr / contrast_sq_norm
+        max_scaled_overall = float(np.max(scaled_max_B))
+        rel_degenerate = (
+            scaled_max_B < 1e-10 * max_scaled_overall
+            if max_scaled_overall > 0
+            else np.zeros(m, dtype=bool)
+        )
+    else:
+        rel_degenerate = np.zeros(m, dtype=bool)
+    at_noise_floor = abs_degenerate | rel_degenerate
+    # Physical-bound guard. The Bell-McCaffrey Satterthwaite DOF is
+    # `(tr B)²/tr(B²)` with `B` PSD and cluster-structured, so it is bounded by
+    # `rank(B) <= G` (the number of clusters) - it can approach G when clusters
+    # are small (few obs each), so the bound is G, NOT the residual dof `n-k`
+    # (which can be smaller than G and would wrongly flag legitimate near-G
+    # DOFs on short panels). A value above G is non-physical. The simple
+    # unweighted `(tr B)²/tr(B²)` form is numerically less faithful than
+    # clubSandwich's P-array form (used on the weighted path) for high-leverage
+    # FE-dummy columns and can return a finite-but-inflated DOF there (observed
+    # ~32.7 and ~16.3 vs R's 6 and 3 on the absorbed-FE golden, G=8) that is NOT
+    # at the noise floor. Rather than ship an impossible DOF we NaN it; exact
+    # clubSandwich reproduction of these non-user-facing nuisance DOFs would
+    # require porting the P-array form and is out of scope (estimators consume
+    # only the well-conditioned treatment / event-study / compound-average
+    # contrasts, which are unaffected).
+    n_clusters = len(unique_clusters)
+    non_physical = np.isfinite(dof_vec) & (dof_vec > n_clusters + 1e-6)
+    degenerate = at_noise_floor | non_physical
+    n_degenerate = int(np.sum(degenerate))
+    if n_degenerate > 0:
+        dof_vec[degenerate] = np.nan
+        warnings.warn(
+            f"Satterthwaite DOF for {n_degenerate} of {m} contrast(s) is "
+            f"unreliable (at the float64 noise floor, or above the cluster-count "
+            f"bound G={n_clusters}); reporting NaN. This affects high-leverage "
+            f"FE-dummy / collinear nuisance coefficients; the resulting DOF is "
+            f"BLAS-order-dependent / non-physical. The coefficient SEs remain "
+            f"valid — only the Satterthwaite DOF (and any t-test / CI depending "
+            f"on it) is suppressed.",
+            UserWarning,
+            stacklevel=3,
+        )
+
     return dof_vec
 
 

diff_diff/utils.py

@@ -352,8 +352,13 @@ def safe_inference(effect, se, alpha=0.05, df=None):
     """
     if not (np.isfinite(se) and se > 0):
         return np.nan, np.nan, (np.nan, np.nan)
-    if df is not None and df <= 0:
-        # Undefined degrees of freedom (e.g., rank-deficient replicate design)
+    if df is not None and not (np.isfinite(df) and df > 0):
+        # Undefined degrees of freedom: df <= 0 (e.g., rank-deficient replicate
+        # design) OR non-finite (a guard-suppressed / non-physical Bell-McCaffrey
+        # Satterthwaite DOF, which `_cr2_bm_dof_inner` returns as NaN for
+        # high-leverage nuisance contrasts). All inference fields are NaN so a
+        # coefficient whose DOF was declared unreliable never reports finite
+        # t/p/CI - preserving the joint-NaN inference contract.
         return np.nan, np.nan, (np.nan, np.nan)
     t_stat = effect / se
     p_value = compute_p_value(t_stat, df=df)