diff --git a/TODO.md b/TODO.md index 96e8d23b9..b0c7e859a 100644 --- a/TODO.md +++ b/TODO.md @@ -22,7 +22,6 @@ Related tracking surfaces: | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| | Consolidate the remaining estimator-entangled DR/logit score variants (`staggered.py::_doubly_robust` + RC twins, `triple_diff.py`, `lwdid.py`, `wooldridge.py`) onto the shared `_dr_scores.py` module, each migration with its own committed oracle capture (the ContinuousDiD lift's two-tier pattern in `tests/test_dr_scores.py`); and add a ridge vcov path to `solve_ridge` if an estimator ever needs analytical ridge inference | `diff_diff/_dr_scores.py` | dml-b0 | Mid | Low | -| hc2/hc2_bm floor `1 - h_ii` at 1e-10 in the shared leverage meat, fabricating finite (if inflated) variances for leverage-one observations - hc3 now fails closed there (LWDiD fix wave) but the pre-existing hc2 family behavior is released surface; decide fail-closed vs keep-floor for hc2/hc2_bm | `diff_diff/linalg.py` | #588 | Quick | Low | | Numeric between-period cohorts (e.g. `first_treat=4.5` with integer times) are rejected by LWDiD while CallawaySantAnna estimates them and LWDiD's own datetime/Period cohorts map to the next observed period — close the dtype asymmetry by adopting the next-observed-period mapping for numeric cohorts too (contract documented in REGISTRY cohort-encodings Note + `docs/api/lwdid.rst` Input Contract). Lands only after PR #588 merges | `diff_diff/lwdid.py` | #588 | Quick | Low | | Implement the LW 2026 eq. 7.9/7.10 unit-average cohort estimand (regress per-unit post-average transformed outcomes on `[1, D_g]` vs never-treated) as an alternative to the documented cell-mass `cohort_effects` convention (REGISTRY within-cohort aggregation Note; the two differ on unbalanced panels, where cell-mass weights units by observed post periods). Needs the 7.10 regression + its covariance on the NT path. Lands only after PR #588 merges | `diff_diff/lwdid_staggered.py` | #588 | Quick | Low | | Expose cell-mass overall ATT (Stata `Post_avg` convention; = CS-simple on balanced panels) as an aggregate extra on LWDiD results — the fit's `.att` is the paper's `tau_omega` (cohort-mean-then-treated-weight, eq. 7.18); the authors' large-N display uses cell-mass weighting instead, and both are legitimate estimands (see the REGISTRY LWDiD Aggregation note). Lands only after PR #588 merges | `diff_diff/lwdid_results.py` | #588 | Quick | Low | diff --git a/changelog.d/20260905-hc2-leverage-one.md b/changelog.d/20260905-hc2-leverage-one.md new file mode 100644 index 000000000..bfe9f7624 --- /dev/null +++ b/changelog.d/20260905-hc2-leverage-one.md @@ -0,0 +1,20 @@ +### Behavioral Changes +- **HC2 and unweighted, unclustered HC2-BM now fail closed at leverage one:** + an effective observation with hat-matrix leverage at least `1 - 1e-8` + produces a warning and entirely NaN covariance (and requested degrees of + freedom), preserving point estimates while suppressing undefined inference. + Python and Rust agree; over-one leverage no longer substitutes HC1. + Older Rust extensions without the fail-closed HC2 capability use NumPy + for HC2 while retaining their other accelerations. + Weighted and clustered HC2-BM retain their separate CR2 conventions, + including all-ones probability weights. + +### Fixed +- **Zero-weight observations do not invalidate HC2/HC3 inference:** excluded + rows contribute zero to the covariance and cannot trigger the leverage + guard. Zero-frequency rows now agree with dropping those rows or expanding + the frequency counts literally. +- **LWDiD uses the shared HC2 covariance guard:** leverage-one regressions + retain their point estimate and unavailable influence contribution while + emitting one covariance warning per regression, without a duplicate local + warning. diff --git a/diff_diff/_backend.py b/diff_diff/_backend.py index 346684601..01a355df9 100644 --- a/diff_diff/_backend.py +++ b/diff_diff/_backend.py @@ -114,13 +114,13 @@ except ImportError: _rust_batched_ridge_chol_solve = None -# HC2 (leverage-corrected) robust vcov: imported independently for the same -# mixed-version reason as demean_map (a stale extension missing only this -# newer symbol degrades HC2 to the NumPy path without disabling the older -# Rust accelerations). +# HC2 requires the v2 fail-closed leverage contract. An older extension can +# export the original symbol yet return finite covariance at unit leverage, +# so symbol presence alone is insufficient. Import v2 independently: legacy +# extensions use NumPy HC2 while retaining every other Rust acceleration. try: from diff_diff._rust_backend import ( - compute_robust_vcov_hc2 as _rust_compute_robust_vcov_hc2, + compute_robust_vcov_hc2_v2 as _rust_compute_robust_vcov_hc2, ) except ImportError: _rust_compute_robust_vcov_hc2 = None diff --git a/diff_diff/estimators.py b/diff_diff/estimators.py index 07d81da20..397cd2b71 100644 --- a/diff_diff/estimators.py +++ b/diff_diff/estimators.py @@ -109,6 +109,8 @@ class DifferenceInDifferences(BaseEstimator): (library default). With ``cluster=``, uses CR1 (Liang-Zeger). - ``"hc2"``: leverage-corrected meat (one-way only). Errors with ``cluster=``; use ``"hc2_bm"`` for clustered Bell-McCaffrey. + Effective leverage ``h_ii >= 1 - 1e-8`` produces a warning and + entirely NaN covariance/inference, retaining point estimates. - ``"hc2_bm"``: one-way HC2 + Imbens-Kolesar (2016) Satterthwaite DOF; with ``cluster=``, Pustejovsky-Tipton (2018) CR2 cluster-robust. ``MultiPeriodDiD(cluster=..., vcov_type="hc2_bm")`` is supported and @@ -116,6 +118,8 @@ class DifferenceInDifferences(BaseEstimator): post-period-average ATT (see ``_compute_cr2_bm_contrast_dof`` in ``linalg.py`` and the REGISTRY.md note). Weighted CR2-BM (``survey_design=`` paths) is a separate gate. + Unweighted, unclustered fits share HC2's leverage-one NaN guard; + clustered CR2 and design-based survey inference are separate. - ``"hc3"``: jackknife-style leverage correction, meat ``e_i^2 / (1 - h_ii)^2`` (one-way only; errors with ``cluster=``). A leverage-one observation has no defined HC3 variance and the @@ -2513,12 +2517,16 @@ class MultiPeriodDiD(DifferenceInDifferences): (library default). With ``cluster=``, uses CR1 (Liang-Zeger). - ``"hc2"``: leverage-corrected meat (one-way only). Errors with ``cluster=``; use ``"hc2_bm"`` without cluster for Bell-McCaffrey. + Effective leverage ``h_ii >= 1 - 1e-8`` produces a warning and + entirely NaN covariance/inference, retaining point estimates. - ``"hc2_bm"``: one-way HC2 + Imbens-Kolesar (2016) Satterthwaite DOF per coefficient plus a contrast-aware DOF for the post-period-average ATT. With ``cluster=``, dispatches to Pustejovsky-Tipton (2018) CR2 cluster-robust with a Bell-McCaffrey Satterthwaite contrast DOF on the post-period average (see ``cluster`` above for parity details). Weighted CR2-BM (``survey_design=``) is still gated. + Unweighted, unclustered fits share HC2's leverage-one NaN guard + for covariance and all contrast DOFs, including the average ATT. - ``"hc3"``: jackknife-style leverage correction, meat ``e_i^2 / (1 - h_ii)^2`` (one-way only; errors with ``cluster=``). A leverage-one observation has no defined HC3 variance and the diff --git a/diff_diff/linalg.py b/diff_diff/linalg.py index 8252222e1..49b6b3324 100644 --- a/diff_diff/linalg.py +++ b/diff_diff/linalg.py @@ -1170,6 +1170,8 @@ def solve_ols( (default). With ``cluster_ids``, dispatches to CR1 (Liang-Zeger). - ``"hc2"``: leverage-corrected meat. One-way only; raises with ``cluster_ids`` (use ``"hc2_bm"`` for clustered Bell-McCaffrey). + An effective observation with ``h_ii >= 1 - 1e-8`` produces a + warning and entirely NaN covariance, retaining point estimates. - ``"hc2_bm"``: HC2 + Imbens-Kolesar (2016) Satterthwaite DOF one-way; Pustejovsky-Tipton (2018) CR2 Bell-McCaffrey with ``cluster_ids``. With ``weights``, dispatches to the clubSandwich WLS-CR2 port — @@ -1177,11 +1179,15 @@ def solve_ols( ``fweight`` raise ``NotImplementedError`` (port matches the ``pweight`` convention only; aweight/fweight derivations are a separate methodology task). + Unweighted, unclustered fits share HC2's leverage-one NaN guard; + weighted/clustered CR2 retains its existing convention, including + all-ones probability weights. - ``"hc3"``: jackknife-style leverage correction, meat ``e_i^2 / (1 - h_ii)^2``. One-way only; raises with ``cluster_ids``. An observation with leverage ``h_ii ~ 1`` has no defined HC3 variance and the vcov fails closed (warning + NaN) - rather than flooring ``1 - h_ii``. + rather than flooring ``1 - h_ii``. As for HC2, zero-weight rows + are excluded; frequency weights use each expanded row's leverage. - ``"conley"``: Conley (1999) spatial-HAC sandwich. Requires ``conley_coords`` (n × 2 array) and ``conley_cutoff_km`` (positive bandwidth, no default per Conley 1999 Section 5's sensitivity-grid @@ -1919,6 +1925,7 @@ def _solve_ols_numpy( _VALID_VCOV_TYPES = frozenset({"classical", "hc1", "hc2", "hc2_bm", "hc3", "conley"}) +_HC_LEVERAGE_THRESHOLD = 1.0 - 1e-8 def _validate_vcov_args( @@ -2159,6 +2166,9 @@ def compute_robust_vcov( ``sum_i (u_i^2 / (1 - h_ii)) x_i x_i'`` where ``h_ii`` are hat-matrix diagonals. No DOF adjustment beyond ``n - k``. One-way only; errors with ``cluster_ids``. + Any effective observation with ``h_ii >= 1 - 1e-8`` yields a warning + and entirely NaN covariance/DOF, as for HC3 and unweighted one-way + HC2-BM. This includes numerically over-one leverage; no HC1 fallback. - ``"hc3"``: jackknife-style leverage correction, meat ``sum_i (u_i^2 / (1 - h_ii)^2) x_i x_i'`` (matches ``sandwich::vcovHC`` type="HC3": no DOF factor). One-way only; errors with ``cluster_ids``. @@ -2171,6 +2181,10 @@ def compute_robust_vcov( with ``G=103``. **Weighted hc2_bm** (both one-way and clustered) is supported for ``weight_type="pweight"`` only via the clubSandwich WLS-CR2 port; ``aweight`` and ``fweight`` raise ``NotImplementedError``. + Weighted/clustered CR2 keeps its own adjustment convention: even + all-ones probability weights can yield finite covariance where the + unweighted one-way leverage guard returns NaN. Weighting is not a + remedy for undefined inference. - ``"conley"``: spatial HAC sandwich (Conley 1999 Eq 4.2). Requires ``conley_coords`` (n×2 array) and ``conley_cutoff_km`` (positive bandwidth). Two operating modes: cross-sectional (default) and panel @@ -2236,7 +2250,8 @@ def compute_robust_vcov( array of per-coefficient degrees of freedom. For ``classical``, ``hc1``, ``hc2``, ``hc3``: every element is ``n_eff - k``. For ``hc2_bm`` one-way: Imbens-Kolesar (2016) Satterthwaite DOF per - contrast. + contrast. The leverage-one fail-closed guard returns an entirely + NaN vector for HC2, HC3, and unweighted one-way HC2-BM. cluster_k_adjustment : int, default 0, keyword-only Signed K_reference adjustment added to the visible column count in the CLUSTERED CR1 finite-sample factor only (absorbed FE not nested @@ -2270,7 +2285,9 @@ def compute_robust_vcov( For HC2 one-way (weighted per review MEDIUM #3): h_ii = w_i * x_i' * (X'WX)^{-1} * x_i (unweighted: w_i = 1) meat = sum_i (u_i^2 / (1 - h_ii)) x_i x_i' - Guards against h_ii > 1 - eps with a fall-back to HC1 plus warning. + For fweights, use each expanded row's leverage without the w_i + multiplier. Exclude zero-weight rows from the guard and meat. + Any effective h_ii >= 1 - 1e-8 yields warning + all-NaN vcov/DOF. For HC2 + Bell-McCaffrey one-way DOF (per Imbens-Kolesar 2016): For each coefficient j, let q_j = X (X'X)^{-1} e_j, let M = I - H. @@ -2286,11 +2303,11 @@ def compute_robust_vcov( weights = _validate_weights(weights, weight_type, X.shape[0]) # Rust HC2 (one-way, unweighted, no DOF): mirrors the NumPy hc2 branch - # exactly (leverage meat, no n/(n-k) factor). The near-singular - # hat-diagonal guard stays Python-side: the kernel returns a sentinel - # error and the documented warn-and-fall-back-to-HC1 fires here, - # identical to the NumPy branch's behavior. Imported independently - # (mixed-version safe) — None on a stale extension. + # exactly (leverage meat, no n/(n-k) factor). At leverage ~1 the + # kernel returns a sentinel error, translated to warning + NaN here, + # identical to the NumPy branch's behavior. _backend imports only the + # versioned fail-closed symbol; a legacy extension uses NumPy HC2 even + # when it exports the original, denominator-flooring HC2 kernel. if ( HAS_RUST_BACKEND and _rust_compute_robust_vcov_hc2 is not None @@ -2305,21 +2322,22 @@ def compute_robust_vcov( return _rust_compute_robust_vcov_hc2(X_c, residuals_c) except ValueError as e: error_msg = str(e) - if "Hat-matrix diagonal exceeds 1" in error_msg: + if error_msg.startswith("HC2 variance is undefined:"): warnings.warn( - f"{error_msg} Falling back to HC1.", + f"{error_msg} Returning NaN vcov.", UserWarning, stacklevel=2, ) - return _compute_robust_vcov_numpy( - X, - residuals, - cluster_ids=None, - weights=None, - weight_type=weight_type, - vcov_type="hc1", - return_dof=return_dof, + return np.full((X.shape[1], X.shape[1]), np.nan) + if "Hat-matrix diagonal exceeds 1" in error_msg: + # Older kernels can signal over-one leverage with this + # sentinel. Never relabel an HC1 fallback as HC2. + warnings.warn( + f"HC2 variance is undefined: {error_msg} Returning NaN vcov.", + UserWarning, + stacklevel=2, ) + return np.full((X.shape[1], X.shape[1]), np.nan) if "Matrix inversion failed" in error_msg: raise ValueError( "Design matrix is rank-deficient (singular X'X matrix). " @@ -2436,8 +2454,9 @@ def _compute_hat_diagonals( ``sandwich::vcovHC(..., type="HC2")`` in R and matches the per-observation effective leverage under WLS. - Returns an ``(n,)`` array. Values are clamped to ``[0, 1 - 1e-10]`` to - guard against numerical `` h_ii > 1`` from near-singular designs. + Returns an unclamped ``(n,)`` array. Covariance callers check effective + observations for leverage ``h_ii >= 1 - 1e-8`` and fail closed; keeping + raw values also exposes numerical ``h_ii > 1`` from unstable designs. """ # Compute x_i' (X'WX)^{-1} x_i via a single solve rather than per-row. # np.linalg.solve(bread, X.T) has shape (k, n); multiplying element-wise by @@ -2456,7 +2475,7 @@ def _compute_hat_diagonals( if weights is not None: h_diag = weights * h_diag # Numerical guard. Do not silently clip values materially exceeding 1 — that - # indicates a real design pathology; the caller warns and falls back. + # indicates a real design pathology; the covariance caller fails closed. return np.asarray(h_diag, dtype=np.float64) @@ -3287,9 +3306,10 @@ def _compute_bm_dof_from_contrasts( Schur-product expansion (see the inline derivation) at ``O(n k^2 + k^3)`` per contrast with NO dense ``n×n`` residual-maker — the prior form's ``O(n^2 k)`` hat build limited it to ``n < 10_000``. A noise-floor - cancellation guard NaNs extreme-leverage contrasts whose expanded - denominator collapses below float precision (mirrors the clustered - scores path's guard). + cancellation guard NaNs contrasts whose expanded denominator collapses + below float precision (mirrors the clustered scores path's guard). + Before that calculation, any ``h_ii >= 1 - 1e-8`` suppresses ALL + unweighted contrast DOFs, matching the design-wide covariance guard. **Weighted** (``weights is not None``): dispatches to the clubSandwich singleton-cluster CR2 reduction (each observation is its own cluster) @@ -3318,7 +3338,9 @@ def _compute_bm_dof_from_contrasts( ------- ndarray of shape (m,) of Satterthwaite DOF per contrast column. NaN when the denominator is non-positive or at/below the cancellation noise - floor (degenerate / extreme-leverage case; see the inline guard note). + floor. On the unweighted path, any ``h_ii >= 1 - 1e-8`` returns an + entirely NaN vector without warning, matching the design-wide covariance + guard (which owns the warning). Weighted CR2 keeps its own convention. """ n, k = X.shape if contrasts.ndim != 2 or contrasts.shape[0] != k: @@ -3336,6 +3358,11 @@ def _compute_bm_dof_from_contrasts( X, cluster_ids_singleton, bread_matrix, contrasts, weights=weights ) + if np.any(h_diag >= _HC_LEVERAGE_THRESHOLD): + # Covariance already warns on this design. Remain silent for callers + # computing additional contrasts, e.g. MultiPeriodDiD's average ATT. + return np.full(contrasts.shape[1], np.nan) + # Unweighted: keep the simple (tr B)² / tr(B²) formula — algebraically # identical backward compatibility with prior unweighted Bell-McCaffrey # output (dense prior evaluation; floating-point-tolerance parity). @@ -3350,7 +3377,7 @@ def _compute_bm_dof_from_contrasts( raise # q has shape (n, m); column j is X @ (bread_inv @ contrasts[:, j]). q = X @ bread_inv_c - one_minus_h = np.maximum(1.0 - h_diag, 1e-10) + one_minus_h = 1.0 - h_diag one_minus_2h = 1.0 - 2.0 * h_diag m = contrasts.shape[1] dof = np.empty(m) @@ -3607,24 +3634,19 @@ def _compute_robust_vcov_numpy( bread_matrix, weights=None if weight_type == "fweight" else weights, ) - # Leverage-one observations make the HC3 leave-one-out residual - # undefined (and HC2 nearly so): flooring 1 - h_ii would fabricate - # an arbitrary finite variance for a perfectly-leveraged point - # (e.g. a single treated unit under [1, D]). HC3 fails closed with - # a NaN vcov instead (LWDiD fix-wave review finding); HC2/HC2-BM - # keep their long-standing floor behavior (released surface; - # pre-existing, tracked separately). This check runs BEFORE the - # generic over-one HC1 fallback below (round-10 review: numerically - # over-one leverage previously escaped into an HC1 result still - # labeled hc3 - h >= 1 - 1e-8 covers h > 1 + 1e-6 entirely). - if vcov_type == "hc3" and np.any(h_diag >= 1.0 - 1e-8): - n_lev1 = int(np.sum(h_diag >= 1.0 - 1e-8)) + # A leverage-one observation has no residual information for HC2/HC3. + # Do not fabricate finite variance by flooring its denominator or + # substituting HC1. The threshold also catches numerical h > 1. + # Zero-weight rows represent no observations, even if their fweight + # quadratic form (which omits the weight multiplier) exceeds one. + active = np.ones(n, dtype=bool) if weights is None else weights > 0 + leverage_one = active & (h_diag >= _HC_LEVERAGE_THRESHOLD) + if np.any(leverage_one): + n_lev1 = int(np.sum(leverage_one)) + family = "HC2-BM" if vcov_type == "hc2_bm" else vcov_type.upper() warnings.warn( - f"HC3 variance is undefined: {n_lev1} observation(s) have " - f"hat-matrix leverage ~1 (a perfectly-leveraged design, " - f"e.g. a single treated unit). Returning NaN vcov; use " - f"vcov_type='classical' exact inference or add treated " - f"units.", + f"{family} variance is undefined: {n_lev1} observation(s) have " + "hat-matrix leverage ~1 (h_ii >= 1 - 1e-8). Returning NaN vcov.", UserWarning, stacklevel=3, ) @@ -3635,44 +3657,26 @@ def _compute_robust_vcov_numpy( # the fail-closed semantics through safe_inference. return nan_vcov, np.full(X.shape[1], np.nan) return nan_vcov - if np.any(h_diag > 1.0 + 1e-6): - # hc2/hc2_bm only: hc3 designs with over-one leverage are - # already caught by the fail-closed guard above. - warnings.warn( - f"Hat-matrix diagonal exceeds 1 (max={h_diag.max():.6f}); " - "the design is near-singular. Falling back to HC1.", - UserWarning, - stacklevel=3, - ) - return _compute_robust_vcov_numpy( - X, - residuals, - cluster_ids=None, - weights=weights, - weight_type=weight_type, - vcov_type="hc1", - return_dof=return_dof, - ) - one_minus_h = np.maximum(1.0 - h_diag, 1e-10) + one_minus_h = np.where(active, 1.0 - h_diag, 1.0) + # Mask before arithmetic: even an extreme residual on an excluded + # row must not overflow or produce 0/0 before its zero contribution. + score_residuals = np.where(active, residuals, 0.0) # HC2 meat: sum_i (u_i^2 / (1 - h_ii)) x_i x_i'; HC3 squares the # leverage denominator (jackknife-style, sandwich::vcovHC type="HC3"). # pweight scaling matches the HC1 convention (w_i * u_i / sqrt(denom) # as score). lev_denom = one_minus_h**2 if vcov_type == "hc3" else one_minus_h if weights is not None and weight_type == "fweight": - factor = weights * (residuals**2) / lev_denom + factor = weights * (score_residuals**2) / lev_denom meat = X.T @ (X * factor[:, np.newaxis]) elif weights is not None and weight_type == "pweight": # pweight scores carry w in the score, so meat = sum (w u / sqrt(denom))^2 x x' - scaled = weights * residuals / np.sqrt(lev_denom) + scaled = weights * score_residuals / np.sqrt(lev_denom) scores_hc2 = X * scaled[:, np.newaxis] meat = scores_hc2.T @ scores_hc2 else: # aweight / unweighted: meat = sum_i (u_i^2 / denom_i) x_i x_i' - factor = (residuals**2) / lev_denom - # Zero out zero-weight rows under aweight (subpopulation invariance) - if weights is not None and np.any(weights == 0): - factor = factor * (weights > 0) + factor = (score_residuals**2) / lev_denom meat = X.T @ (X * factor[:, np.newaxis]) # Sandwich without DOF adjustment for HC2/HC3 (matches sandwich::vcovHC @@ -4403,6 +4407,11 @@ class LinearRegression: stores per-coefficient BM Satterthwaite DOF (``self._bm_dof``) and threads it into ``get_inference``. + HC2, HC3, and unweighted unclustered HC2-BM retain coefficients but + warn and return entirely NaN covariance at effective leverage + ``h_ii >= 1 - 1e-8``; coefficient SEs, tests, and intervals are NaN. + Weighted/clustered HC2-BM keeps its separate CR2 convention. + For ``"conley"`` (Conley 1999 spatial-HAC) two operating modes are supported on the `LinearRegression` / `compute_robust_vcov` surface: cross-sectional (single-period or pooled cross-section) and panel diff --git a/diff_diff/lwdid.py b/diff_diff/lwdid.py index 1e8168104..a78eb684a 100644 --- a/diff_diff/lwdid.py +++ b/diff_diff/lwdid.py @@ -2956,10 +2956,10 @@ def _ols_treatment_influence( if self.vcov_type in ("hc2", "hc3"): raw_leverage = np.sum((X @ xtx_inv) * X, axis=1) if self.vcov_type in ("hc2", "hc3") and np.any(raw_leverage >= 1.0 - 1e-8): - # Match the shared linalg fail-closed contract (round-10 - # review: clipping fabricated a finite HC3 influence - # vector for a design whose HC3 vcov is NaN, so aggregate - # inference disagreed with the cell's own). + # Match the shared HC2/HC3 covariance guard. Keep this + # safeguard when reconstructing influence contributions so + # aggregate inference cannot reuse a finite vector from a + # design whose covariance is unavailable. return np.full_like(psi, np.nan) leverage = np.clip(raw_leverage, 0.0, 1.0 - 1e-10) if self.vcov_type == "hc2": @@ -3129,25 +3129,6 @@ def _estimate_reg( used_scales[used_scales == 0] = 1.0 X_scaled = X_used / used_scales xtx_inv = np.linalg.pinv(X_scaled.T @ X_scaled) / np.outer(used_scales, used_scales) - if self.vcov_type == "hc2" and cluster_ids is None: - # Round-21 review: the shared hc2 kernel keeps its RELEASED - # 1 - h floor (tracked separately), but the NEW LWDiD surface - # must not report a fabricated finite variance for a - # perfectly-leveraged design - fail closed HERE, mirroring - # hc3 (point retained, inference NaN). - leverage_used = np.sum((X_used @ xtx_inv) * X_used, axis=1) - if np.any(leverage_used >= 1.0 - 1e-8): - n_lev1 = int(np.sum(leverage_used >= 1.0 - 1e-8)) - warnings.warn( - f"HC2 variance is undefined for this design: {n_lev1} " - f"observation(s) have hat-matrix leverage ~1 (e.g. a " - f"single treated unit). Returning NaN inference (point " - f"retained); use vcov_type='classical' exact inference " - f"or add treated units.", - UserWarning, - stacklevel=2, - ) - se = np.nan influence = self._finalize_influence( self._ols_treatment_influence( X_used, diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index e19c8601c..a6c0b0640 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -87,6 +87,8 @@ where τ is the ATT. - With "warn" (default): emits warning, sets NaN for affected coefficients - With "error": raises ValueError - With "silent": continues silently with NaN coefficients +- **Note:** Shared HC2/HC3 leverage-one policy: HC2, HC3, and unweighted, unclustered HC2-BM warn and return an entirely NaN covariance when any effective observation has `h_ii >= 1 - 1e-8`, including numerical leverage above one. Requested covariance DOF vectors and all unweighted BM contrast DOFs are also entirely NaN; the contrast-only helper stays silent because covariance computation owns the diagnostic. This is a defensive threshold convention: flooring the undefined leverage correction or substituting HC1 could fabricate finite inference. Identified coefficients, fitted values, and residuals are retained; SEs, t-statistics, p-values, and confidence intervals are unavailable. The policy applies to all causes of unit leverage, including singleton treatment cells and singleton fixed-effect levels. Zero-weight rows are excluded from both the guard and meat; under frequency weights the leverage is each expanded row's `x_i'(X'WX)^{-1}x_i` (without a weight multiplier), so zero-frequency rows cannot invalidate HC2 or HC3 and compressed covariance/DOF matches literal frequency expansion. Probability and analytical weights retain the WLS leverage `w_i x_i'(X'WX)^{-1}x_i`. +- **Note:** HC2-BM CR2 boundary: weighted or clustered HC2-BM retains its existing CR2 adjustment and DOF conventions and does not use this observation-level fail-closed guard. Consequently, supplying all-ones probability weights can retain finite CR2 covariance where an otherwise identical unweighted one-way HC2-BM call now returns NaN. These are distinct dispatch conventions at the singular-leverage boundary; weighting is not a remedy for undefined inference. Healthy-design formulas and weighted/clustered CR2 reference goldens are unchanged. - Singleton clusters (one observation): included in variance estimation; contribute to meat matrix via u_i² X_i X_i' (same formula as larger clusters with n_g=1) - 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 @@ -298,7 +300,7 @@ where V is the VCV sub-matrix for post-treatment δ_e coefficients. `_compute_cr2_bm_contrast_dof` (DOF-only for arbitrary contrasts) are thin wrappers over that shared core, so every CR2 caller routes through one implementation. The consolidation is bit-identical to the prior two-call path (proven at atol=0). -- **Note (unweighted per-coef DOF guard):** the unweighted clustered CR2-BM per-coefficient DOF (`_cr2_bm_dof_inner`, the simple `(tr B)²/tr(B²)` form) carries the same two-part reliability guard as the weighted P-array path. (1) **Noise floor:** for a high-leverage FE-dummy / collinear nuisance column `trace_B2 = Σ B_{g,h}²` collapses to float64 accumulation noise while `trace_B` stays O(1), inflating the ratio to a non-physical DOF (observed ~1e61 on the absorbed-FE golden); a contrast whose `max|B_{g,h}|` sits below the batch-relative (`1e-10×max`, computed on the scale-normalized `max|B|/‖c‖²` since `B ∝ ‖c‖²` while the DOF is scale-invariant) or absolute (`(EPS·n·k·bread_scale)²`) floor is NaN'd. (2) **Cluster-count bound:** the Bell-McCaffrey Satterthwaite DOF is `(tr B)²/tr(B²)` with `B` PSD and cluster-structured, so it is bounded by `rank(B) ≤ G` (number of clusters); the simple unweighted form is numerically less faithful than clubSandwich's P-array form on high-leverage columns and can return a finite-but-inflated DOF above `G` (observed ~32.7 and ~16.3 vs R's 6 and 3, `G=8`), which is NaN'd as non-physical. The well-conditioned contrasts estimators consume — the treatment effect, event-study coefficients, and the compound post-period-average — are unaffected and match R clubSandwich; only the non-user-facing high-leverage nuisance DOFs are suppressed (exact P-array reproduction of those is deferred). A `UserWarning` fires per fit. Regression: `tests/test_estimators_vcov_type.py::TestDiDAbsorbedFERParity::test_unweighted_cr2_bm_per_coef_dof_no_nonphysical`. **Evaluation (2026-07):** the pairwise `B` matrix is computed via the algebraic identity `B = Ω'MΩ = diag(‖ω_g‖²) − P' M_U P` with `P = X'Ω` (the Satterthwaite DOF itself is Pustejovsky-Tipton 2018 §3.1 Eq. 11; equivalently the `q=1` reduction of the AHT Eq. 13) — instead of contracting against an explicit dense `n×n` residual-maker (`O(n k + G² k)` per contrast; peak memory = two `O(n k)` input-scale score precomputes plus working buffers capped by a byte cap subject to a one-contrast lower bound (a single contrast intrinsically needs `O(n)` + `O(G k)` buffers) — q vectors, per-cluster omegas, and product buffers contrast-chunked with every width-scaled buffer counted in the chunk denominator, and the `(G, G)` pairwise matrix row-chunked (Frobenius sum and max are row-separable), so none of the `n×n` residual-maker, `O(n m)` score arrays, or `O(G²)` pairwise entries is ever held at once; chunk-count invariant to ~1 ULP — BLAS width-dependent column accumulation; algebraically identical, ~1e-12 float agreement locked by a frozen-oracle parity test). **One-way scores evaluation (2026-07):** the non-clustered unweighted BM DOF denominator `a'(M∘M)a` (`_compute_bm_dof_from_contrasts`) is likewise evaluated without the dense `n×n` residual-maker, via the Schur-product expansion `Σ a_i²(1−2h_ii) + tr((B S_a)²)` with `S_a = X'diag(a)X` — exact algebra, `O(n k² + k³)` per contrast (was `O(n²k)`), frozen-dense-oracle parity ~1e-12, with a noise-floor cancellation guard NaN-ing extreme-leverage contrasts whose expanded denominator collapses below float precision (`TestOneWayBMScoresDOF`). **Low-rank A_g factorization (2026-07):** the unweighted per-cluster adjustment operator `A_g = (I − H_gg)^{−1/2}` is evaluated from the k×k eigenproblem of `U_g'U_g` (`U_g = X_g M_U^{1/2}`; `A_g = I + (U_g Q) diag(γ) (U_g Q)'`, `γ = ((1−λ)^{−1/2}−1)/λ` via expm1/log1p, Moore-Penrose zeroing at `1−λ ≤ 1e-10` matching `_cr2_adjustment_matrix`'s convention) and only ever applied to skinny matrices — the dense `(n_g, n_g)` `A_g` (an `O(n_g³)` eigh per cluster) is never materialized for `n_g > k`; clusters with `n_g ≤ k` keep the (tiny) dense construction since it is the smaller Gram side. Algebraically identical; vcov/DOF match a frozen dense-eigh oracle at rtol 1e-12/1e-10 incl. the leverage-1 absorbed-cluster-FE and singleton-cluster regimes (`TestCR2BMLowRankAdjustment`). The weighted clubSandwich path keeps the dense construction (its `G_g` carries the `S_W` bias term). +- **Note (unweighted per-coef DOF guard):** the unweighted clustered CR2-BM per-coefficient DOF (`_cr2_bm_dof_inner`, the simple `(tr B)²/tr(B²)` form) carries the same two-part reliability guard as the weighted P-array path. (1) **Noise floor:** for a high-leverage FE-dummy / collinear nuisance column `trace_B2 = Σ B_{g,h}²` collapses to float64 accumulation noise while `trace_B` stays O(1), inflating the ratio to a non-physical DOF (observed ~1e61 on the absorbed-FE golden); a contrast whose `max|B_{g,h}|` sits below the batch-relative (`1e-10×max`, computed on the scale-normalized `max|B|/‖c‖²` since `B ∝ ‖c‖²` while the DOF is scale-invariant) or absolute (`(EPS·n·k·bread_scale)²`) floor is NaN'd. (2) **Cluster-count bound:** the Bell-McCaffrey Satterthwaite DOF is `(tr B)²/tr(B²)` with `B` PSD and cluster-structured, so it is bounded by `rank(B) ≤ G` (number of clusters); the simple unweighted form is numerically less faithful than clubSandwich's P-array form on high-leverage columns and can return a finite-but-inflated DOF above `G` (observed ~32.7 and ~16.3 vs R's 6 and 3, `G=8`), which is NaN'd as non-physical. The well-conditioned contrasts estimators consume — the treatment effect, event-study coefficients, and the compound post-period-average — are unaffected and match R clubSandwich; only the non-user-facing high-leverage nuisance DOFs are suppressed (exact P-array reproduction of those is deferred). A `UserWarning` fires per fit. Regression: `tests/test_estimators_vcov_type.py::TestDiDAbsorbedFERParity::test_unweighted_cr2_bm_per_coef_dof_no_nonphysical`. **Evaluation (2026-07):** the pairwise `B` matrix is computed via the algebraic identity `B = Ω'MΩ = diag(‖ω_g‖²) − P' M_U P` with `P = X'Ω` (the Satterthwaite DOF itself is Pustejovsky-Tipton 2018 §3.1 Eq. 11; equivalently the `q=1` reduction of the AHT Eq. 13) — instead of contracting against an explicit dense `n×n` residual-maker (`O(n k + G² k)` per contrast; peak memory = two `O(n k)` input-scale score precomputes plus working buffers capped by a byte cap subject to a one-contrast lower bound (a single contrast intrinsically needs `O(n)` + `O(G k)` buffers) — q vectors, per-cluster omegas, and product buffers contrast-chunked with every width-scaled buffer counted in the chunk denominator, and the `(G, G)` pairwise matrix row-chunked (Frobenius sum and max are row-separable), so none of the `n×n` residual-maker, `O(n m)` score arrays, or `O(G²)` pairwise entries is ever held at once; chunk-count invariant to ~1 ULP — BLAS width-dependent column accumulation; algebraically identical, ~1e-12 float agreement locked by a frozen-oracle parity test). **One-way scores evaluation (2026-07):** the non-clustered unweighted BM DOF denominator `a'(M∘M)a` (`_compute_bm_dof_from_contrasts`) is likewise evaluated without the dense `n×n` residual-maker, via the Schur-product expansion `Σ a_i²(1−2h_ii) + tr((B S_a)²)` with `S_a = X'diag(a)X` — exact algebra, `O(n k² + k³)` per contrast (was `O(n²k)`), frozen-dense-oracle parity ~1e-12, with a noise-floor cancellation guard NaN-ing contrasts whose expanded denominator collapses below float precision. A separate design-wide guard now returns all-NaN unweighted contrast DOFs silently at any `h_ii >= 1 - 1e-8`, matching the shared covariance failure (`TestOneWayBMScoresDOF`; healthy dense-oracle expectations retained). **Low-rank A_g factorization (2026-07):** the unweighted per-cluster adjustment operator `A_g = (I − H_gg)^{−1/2}` is evaluated from the k×k eigenproblem of `U_g'U_g` (`U_g = X_g M_U^{1/2}`; `A_g = I + (U_g Q) diag(γ) (U_g Q)'`, `γ = ((1−λ)^{−1/2}−1)/λ` via expm1/log1p, Moore-Penrose zeroing at `1−λ ≤ 1e-10` matching `_cr2_adjustment_matrix`'s convention) and only ever applied to skinny matrices — the dense `(n_g, n_g)` `A_g` (an `O(n_g³)` eigh per cluster) is never materialized for `n_g > k`; clusters with `n_g ≤ k` keep the (tiny) dense construction since it is the smaller Gram side. Algebraically identical; vcov/DOF match a frozen dense-eigh oracle at rtol 1e-12/1e-10 incl. the leverage-1 absorbed-cluster-FE and singleton-cluster regimes (`TestCR2BMLowRankAdjustment`). The weighted clubSandwich path keeps the dense construction (its `G_g` carries the `S_W` bias term). - **Note:** `LinearRegression.get_se()` / `get_inference()` clamp the vcov diagonal at 0 before `sqrt`. A high-leverage / degenerate coefficient (an absorbed-FE dummy near-collinear with the treatment, whose Satterthwaite DOF already hits the noise-floor @@ -2743,7 +2745,7 @@ Event-study/placebo transformations over ALL periods (Appendix D): demeaning (D. - **Note (review round 7: replay follows the fitted design; sensitivity runs the design check):** the post-fit replay mirrors `_estimate_reg`'s LW eq. 3.3 interaction gate (`N_1 > K+1` AND `N_0 > K+1`): small-arm fits use the plain `(1, D, X)` design and their replayed RI/WCR statistic matches `.att` (pre-fix the replay always interacted, so the round-5 coherence assert made small-arm fits' post-fit inference unusable — the fail-closed backstop working as designed, now with the correct design selected). `_prevalidate_frame` in the sensitivity helpers runs the full treatment-design check (absorbing treatment, common-timing onset homogeneity, D_it/cohort consistency, with fit's encode-then-normalize ordering), so structural design violations RAISE instead of being swallowed by the per-spec ValueError handler as `not_estimable`. - **Note (review round 8: calendar partition from S, tau_omega window counts, aweight leverage-family convention):** the common-timing pre/post partition derives from the SINGLE adoption period `S = min(observed treated period)`: `pre = {t < S}`, `post = {t >= S}` (the pre-fix per-period `max(D)` partition classified a post period with no observed treated rows as PRE-treatment, contaminating the rolling pre window — execution-verified: zero-effect trend panel biased to ATT 0.75); the common-timing design check likewise validates `D_it = 1[t >= S]` over observed rows, so a unit whose `t = S` row is missing is accepted (matching the staggered branch) while genuinely heterogeneous onsets still raise. The `tau_omega` complete-case semantics are CLARIFIED, not changed (round-8 reviewer proposed full-window per-period counting; NOT adopted): a unit contributes cohort g's component iff its OBSERVED post-g rows yield a finite average — partial post windows are averaged over observed rows, symmetrically for treated and control units. This is the adjudicated WS1 design pinned byte-frozen by the acceptance suite's independent reference oracle (`_complete_case_tau_omega_reference`) and its zero-drop metadata test; changing to every-period counting would change the estimand those tests pin. CAVEAT (documented): on unbalanced panels with time trends, differential post-period availability enters the composite through the observed-window averages — complete-case drops fire only when a required window is entirely missing/non-finite. Covariates must be numeric and FINITE at the front door (Inf passed the NaN check and was silently cell-filtered); `validate_staggered_data` rejects datetime64-vs-Period mixtures and Period-frequency mismatches exactly like the encoding step. **aweight + hc2/hc3 refutation (round 8):** the reviewer's proposed `w^2` score meat contradicts the documented aweight convention (this section, Weight Type Effects: "aweights use unweighted meat ... matches Stata convention" — known-heteroskedasticity WLS leaves ~homoskedastic errors) — hc2's aweight surface is RELEASED behavior retained byte-identical from main, and hc3 follows the same family branch (unweighted meat with the WLS-hat leverage). The convention is deliberate and documented, not a defect. - **Note (review round 9: one event-time convention, onset partition propagated, finite outcomes):** event-time labels follow ONE convention across the common and staggered interfaces: NUMERIC calendars use the Registry's arithmetic `r = t - g` (validated INTEGRAL — a fractional horizon raises instead of silently merging under the integer storage keys, which previously overwrote distinct horizons' estimates and covariance entries via `int(t - g)`); datetime/Period calendars use position differences on the ordered support (they are position-encoded before the staggered machinery). Pre-fix, the common interface used positional labels for ALL dtypes, so a gapped numeric calendar got different event keys per interface ({0,1} vs {0,2} on {1,2,4,6} with onset 4). The round-8 onset partition (`pre = {t < S}`) is propagated to `get_transformation_diagnostics` and the sensitivity helpers' pre/post sets (a controls-only post period is post everywhere; sensitivity subsets retain every `t >= S` period). Outcomes must be numeric and FINITE at the front door in both timing modes (Inf previously passed the NaN check and was silently np.isfinite-filtered inside staggered cells, changing the estimation sample without warning). -- **Note (review round 10: guard ordering + sensitivity coherence):** hc3's undefined-leverage fail-closed check runs BEFORE the generic over-one HC1 fallback (numerically over-one leverage previously escaped into an HC1 result still labeled hc3), and the LWDiD hc3 influence vector fails closed to NaN under the same condition instead of clipping (aggregate inference matches the cell's NaN vcov; the NaN influence drops the cell from joint aggregation). The sensitivity helpers count treated cohorts on the NORMALIZED frame (beyond-window/inf encodings no longer masquerade as extra cohorts), their BASELINE full-frame fit propagates every fit error (a configuration/support failure such as covariate-free PSM raises instead of reporting `not_estimable`; only restricted-subset fits map failures to NaN specs), and zero-post-row units are counted by the fixed-window drop warning (previously they vanished silently in the merge). +- **Note (review round 10: guard ordering + sensitivity coherence):** the shared HC2/HC3 and unweighted one-way HC2-BM guards fail closed at effective leverage `h_ii >= 1 - 1e-8`, including numerical over-one leverage, without an HC1 fallback. The LWDiD HC2/HC3 influence vector independently fails closed to NaN under the same condition instead of clipping (aggregate inference matches the cell's NaN vcov; the NaN influence drops the cell from joint aggregation). The sensitivity helpers count treated cohorts on the NORMALIZED frame (beyond-window/inf encodings no longer masquerade as extra cohorts), their BASELINE full-frame fit propagates every fit error (a configuration/support failure such as covariate-free PSM raises instead of reporting `not_estimable`; only restricted-subset fits map failures to NaN specs), and zero-post-row units are counted by the fixed-window drop warning (previously they vanished silently in the merge). - **Note (review round 11: propensity linearization, reduced-rank propensity fits, identified-rank gate):** the IPW/DR influence functions build the logit score and Hessian from the RAW fitted probabilities (the actual MLE's estimating equation — its score is ~0 at the fit; the pre-fix code used the CLIPPED probabilities, breaking the linearization whenever `pscore_trim` fired), and the weight-derivative `dw/dgamma` is ZERO for clipped observations (a clipped weight is locally constant in gamma); the clipped probabilities remain the WEIGHTING choice for the point estimator. A rank-deficient propensity model (NaN logit coefficients from dropped collinear columns, finite probabilities) CONTINUES as an IPW/DR fit on the reduced-rank propensity (score/Hessian on the kept columns) — the pre-fix code silently substituted regression adjustment under ipw/dr provenance; only a genuinely failed solve (non-finite probabilities) falls back, with its warning. The RA interaction gate (`N_1 > K+1` and `N_0 > K+1`, eq. 3.3) counts the IDENTIFIED control dimension (matrix rank), mirrored exactly by the post-fit replay — a perfectly collinear control previously flipped the gate and changed the ATT while adding no information. `validate_staggered_data` marks duplicate `(unit, time)` cells invalid (a duplicate could mask a missing cell in the row-count balance check). Tutorial 27 re-executed against the final code: the single-treated California HC3 example now TEACHES the leverage-one fail-closed boundary (classical exact-t / RI are the small-N tools), and the CS-efficiency comparison states the paper-faithful serial-correlation trade-off instead of a dominance claim. - **Note (review round 12: rank-aware DR nuisances, identified parameter counts, survivor cohort masses):** the DR outcome WLS fits through the shared rank-aware solver and every outcome-model influence term (prediction, `S_beta`, `H_beta`, `dATT/dbeta`) uses the IDENTIFIED column mask (the pre-fix raw `inv`/`pinv` Gram was not scale-equilibrated — an exactly redundant 1e12-rescaled duplicate changed the DR SE by ~2.5x); IPW/DR report the identified propensity/outcome ranks as `n_params` (nominal counts let a redundant control shrink residual df and move p-values/CIs). On the tau_omega DROPS route, `.att` and its combined influence function weight cohorts by the SURVIVING cohort masses returned by the composite helper (the Registry complete-case rule; raw masses previously kept dropped treated units in the weights) — pinned by an independent survivor-mass oracle. - **Note (review round 13: scale-equilibrated influence bread, effective-rank guard):** the RA influence reconstruction inverts the COLUMN-EQUILIBRATED Gram and unscales (`(X'X)^{-1} = D^{-1}(Xs'Xs)^{-1}D^{-1}`) — the pre-fix raw-Gram pinv silently dropped low-scale directions at large covariate units, so cell ATT/SE (from the equilibrated `solve_ols`) were unit-invariant while every AGGREGATE SE/p/CI and the multiplier-bootstrap inputs were not (execution class: rescaling one covariate by 1e7 moved the overall SE from 0.128 to 0.028 with no warning). Aggregate-inference unit-invariance is pinned across the overall and event-study surfaces. The exact-inference small-sample guard uses the EFFECTIVE (equilibrated) design rank, so a redundant-column design with positive effective residual df fits while a genuinely saturated design still raises. docs/index.rst and the practitioner tree scope the heterogeneous-trends claim to `rolling='detrend'` and describe PSM as point-estimation-only. @@ -2755,7 +2757,7 @@ Event-study/placebo transformations over ALL periods (Appendix D): demeaning (D. - **Note (review round 18):** the common-timing time-scale contract (Period rejected for detrend/detrendq; trend/seasonal transforms require numeric/datetime/Period time) lives in one shared validator called by BOTH `fit()` and `get_transformation_diagnostics()` (diagnostics previously reached the transforms' raw float-conversion errors). `randomization_inference` validates array shapes/lengths BEFORE the non-finite-outcome filter (a mismatched length combined with a non-finite y previously raised a raw boolean-index IndexError). RI citations point at the LW 2026 small-sample paper (the 2025 Section-5 reference concerned detrending, not RI), and the api-docs no longer call RI "assumption-free" (it does not require normality, conditional on the complete-randomization assignment mechanism). - **Note (review round 19: family-consistent multiplier contributions — deliberate, externally validated):** the influence contributions feeding the event-study multiplier bootstrap are NORMALIZED TO THE REQUESTED ANALYTICAL VARIANCE FAMILY (classical: per-cell scalar rescale to the classical magnitude; hc1/CR1: the small-sample factor; hc2/hc3: leverage adjustment), not the raw Appendix E.2 contributions. Consequences: per-cell SCALAR adjustments (classical/hc1/CR1) leave the sup-t critical value INVARIANT (draws and SEs scale together and the normalized statistic cancels the factor) while the per-event bootstrap SEs report magnitudes consistent with the requested family rather than the raw asymptotic form — a deliberate coherence choice, so a fit's analytical and bootstrap surfaces answer in the same family. External validation: the RA/hc1 configuration's multiplier-bootstrap SEs are gated against the AUTHORS' Stata package's high-B multiplier bootstrap within the Monte-Carlo bound (acceptance suite, `test_walmart_eventstudy_se_vs_stata`). PSM continues under a rank-deficient (finite-probability) propensity fit exactly like ipw/dr — matching needs only the probabilities — with the regression-point fail-closed fallback reserved for genuinely non-finite propensity fits. - **Note (review round 20):** point-only PSM is EXEMPT from the common-timing exact-OLS residual-df guard (its inference is NaN by contract and `df_inference=None`; the guard previously rejected valid matching fits whose nominal propensity width exhausted an OLS df count PSM never uses — the staggered path already retained the point). `get_transformation_diagnostics` rejects an all-never-treated staggered panel ("No treated cohorts found", matching `fit_staggered`) instead of returning an empty `by_cohort` that read as success. Non-numeric, non-datetime common-timing time columns must be ORDERED CATEGORICALS (encoded to their codes before any comparison; plain object labels are rejected — lexicographic order breaks at 'Q10' < 'Q2', silently corrupting the pre/post partition), replacing the former plain-string demean acceptance; declared-order chronology pinned on a Q1..Q10 zero-effect trend panel. -- **Note (review round 21):** the NEW LWDiD surface fails closed for `vcov_type='hc2'` at leverage-one designs (warning + NaN inference, point retained), mirroring hc3 — the SHARED hc2/hc2_bm kernel keeps its released `1 - h` floor for the pre-existing estimators pending the tracked family decision (TODO row), so the fabricated-variance path is unreachable from LWDiD while released surfaces are unchanged. The common-timing HEADLINE and unit-bootstrap SEs run through the same scale-equivariant degenerate-SE guard as the staggered/event surfaces (an exactly fitted panel previously reported se ~ 1e-16 with t ~ 1e16). LWDiD plots render the FITTED interval endpoints (per-row t/df, fitted alpha, cband when present; sensitivity specs now carry `conf_int`) instead of a fabricated normal-theory `+/-1.96*SE`; inference-unavailable rows keep the omit-interval rule. PSM docstrings describe 1:`n_neighbors` matching (1:1 default). +- **Note (review round 21):** LWDiD fails closed for `vcov_type='hc2'` at leverage-one designs (warning + NaN inference, point retained), mirroring hc3. The shared covariance kernel now enforces this policy for HC2 and unweighted one-way HC2-BM as well; LWDiD relies on that diagnostic and its existing covariance-to-SE propagation, with no duplicate local covariance guard. The independent influence-function safeguard is retained so these regressions supply no usable influence contribution to joint aggregation. The common-timing HEADLINE and unit-bootstrap SEs run through the same scale-equivariant degenerate-SE guard as the staggered/event surfaces (an exactly fitted panel previously reported se ~ 1e-16 with t ~ 1e16). LWDiD plots render the FITTED interval endpoints (per-row t/df, fitted alpha, cband when present; sensitivity specs now carry `conf_int`) instead of a fabricated normal-theory `+/-1.96*SE`; inference-unavailable rows keep the omit-interval rule. PSM docstrings describe 1:`n_neighbors` matching (1:1 default). - **Note (review round 23):** staggered overall/cohort masses count treated units CONTRIBUTING to each cohort's estimable post cells on every route (a raw cohort member with no estimable cell no longer raises its cohort's weight in `.att`, its combined influence function, or `cohort_effects[g]['n_treated']` — previously only the tau_omega drops route recomputed masses); pinned by a contributing-unit oracle (4/5 vs the raw 4/8 weighting on a mostly-unobserved cohort). `plot_cohort_trends(cohort=)` is IMPLEMENTED (one trajectory per treated cohort, never-treated control line, per-cohort onset markers — previously the parameter was accepted and silently ignored). `validate_staggered_data` marks missing unit/time values as ERRORS (fit rejects the same frame; warning-only let `valid: True` disagree with fit). Input-contract docs state the unit-constant covariate/cluster rule applies to BOTH timing paths and add `_lwdid_season` to the reserved-name list. - **Note (review round 24, final):** verdict "Looks good — no unmitigated P0 or P1 findings"; the three P2 nits are resolved: `randomization_inference` requires `n_reps >= 10` up front (the reliable-inference floor made smaller values fail after the permutation loop with a misleading hint); the HC3 leverage-one fail-closed path honors the `return_dof` contract with a length-k NaN vector (was `None`); datetime cohort positions relabel to the CANONICAL observed period (two raw between-period labels mapping to the same onset previously collided with row-order-dependent survivor). - **Note (round-2 refutations, evidence-anchored):** two reviewer claims were checked and REFUTED by execution: (1) pre-treatment placebo transformations — the implementation applies one per-cohort transformation over the full `t < g` pre window with anchor exclusions and the D.3 placebo control pools, and matches the authors' Stata `lwdid` 2.4.2 at full precision (~1e-9) on every Walmart placebo cell `r in [-22, -3]` (the fail-closed label-set gate pins the surface), so the horizon-specific future-window reading is not what the reference implementation does; (2) the IPW influence function's `p_bar = n_1/n` normalization (Lunceford-Davidian linearization of the Hajek ATT) was compared against the proposed finite-sample `B_hat = sum_ctrl(w)/n` variant by Monte Carlo (400 reps, strong propensity heterogeneity): the variants are first-order equivalent and `B_hat` calibrated no better (SE/SD 0.865 vs 0.873), so the implemented convention stands. diff --git a/docs/methodology/variance-conventions.md b/docs/methodology/variance-conventions.md index ee99a8c92..a975fff8e 100644 --- a/docs/methodology/variance-conventions.md +++ b/docs/methodology/variance-conventions.md @@ -149,6 +149,8 @@ output). family applies no CR1 finite-sample factor, so it has no cell on the axis this matrix measures. + **Effective-leverage guard:** HC2, HC3, and unweighted one-way HC2-BM return warning + all-NaN covariance/DOF when any positive-weight observation has `h_ii >= 1 - 1e-8` (all rows are effective without weights); identified point estimates remain. Zero-weight rows contribute nothing and cannot trigger the guard, preserving frequency-expansion parity even with zero counts. Weighted/clustered HC2-BM retains CR2, including the all-ones-pweight boundary; its output can differ from unweighted HC2-BM at leverage one. No over-one HC1 fallback is used in the leverage-family branch. + ## Tail-df landscape (converged, 3.9 / M-127) ONE three-value knob (`df_convention ∈ {"residual", "cluster", "normal"}`) diff --git a/rust/src/lib.rs b/rust/src/lib.rs index 7ad95f9a3..6f504384a 100644 --- a/rust/src/lib.rs +++ b/rust/src/lib.rs @@ -55,6 +55,13 @@ fn _rust_backend(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(linalg::solve_ols_chol, m)?)?; m.add_function(wrap_pyfunction!(linalg::compute_robust_vcov, m)?)?; m.add_function(wrap_pyfunction!(linalg::compute_robust_vcov_hc2, m)?)?; + // Versioned capability: Python must not dispatch HC2 to older kernels + // that silently floor the leverage denominator. Keep the original name + // as an alias for callers using the extension directly. + m.add( + "compute_robust_vcov_hc2_v2", + m.getattr("compute_robust_vcov_hc2")?, + )?; // Batched ridge-regularized SPD solve (EfficientDiD per-unit weights) m.add_function(wrap_pyfunction!( diff --git a/rust/src/linalg.rs b/rust/src/linalg.rs index 660c6af73..86db4b291 100644 --- a/rust/src/linalg.rs +++ b/rust/src/linalg.rs @@ -242,12 +242,11 @@ pub fn compute_robust_vcov<'py>( /// Mirrors the NumPy `_compute_robust_vcov_numpy` unweighted `hc2` branch /// exactly (sandwich::vcovHC type="HC2" convention): /// h_i = x_i' (X'X)^{-1} x_i -/// meat = X' diag(u_i^2 / max(1 - h_i, 1e-10)) X +/// meat = X' diag(u_i^2 / (1 - h_i)) X /// vcov = (X'X)^{-1} meat (X'X)^{-1} (NO n/(n-k) factor) -/// A hat diagonal exceeding 1 + 1e-6 signals a near-singular design; this -/// returns the sentinel error "Hat-matrix diagonal exceeds 1" so the Python -/// dispatcher can reproduce the documented warn-and-fall-back-to-HC1 -/// behavior (the guard decision stays in one place, Python-side). +/// A hat diagonal >= 1 - 1e-8 makes HC2 variance unavailable. Return a +/// dedicated sentinel error so the Python dispatcher emits a UserWarning +/// and returns an all-NaN covariance, matching the NumPy guard. /// /// # Arguments /// * `x` - Design matrix (n, k) @@ -281,19 +280,20 @@ pub fn compute_robust_vcov_hc2<'py>( let x_bread = x_arr.dot(&xtx_inv); // (n, k) let h_diag: Array1 = (&x_bread * &x_arr).sum_axis(Axis(1)); - let h_max = h_diag.iter().cloned().fold(f64::NEG_INFINITY, f64::max); - if h_max > 1.0 + 1e-6 { + // Keep the cutoff aligned with Python's _HC_LEVERAGE_THRESHOLD. + let n_leverage_one = h_diag.iter().filter(|&&h| h >= 1.0 - 1e-8).count(); + if n_leverage_one > 0 { return Err(PyErr::new::(format!( - "Hat-matrix diagonal exceeds 1 (max={:.6}); the design is near-singular.", - h_max + "HC2 variance is undefined: {} observation(s) have hat-matrix leverage ~1 (h_ii >= 1 - 1e-8).", + n_leverage_one ))); } - // meat = X' diag(u^2 / max(1 - h, 1e-10)) X + // meat = X' diag(u^2 / (1 - h)) X; the guard ensures positive denominators. let factor: Array1 = residuals_arr .iter() .zip(h_diag.iter()) - .map(|(u, h)| u * u / (1.0 - h).max(1e-10)) + .map(|(u, h)| u * u / (1.0 - h)) .collect(); let factor_col = factor.insert_axis(Axis(1)); // (n, 1) let x_weighted = &x_arr * &factor_col; // (n, k) diff --git a/tests/test_estimators_vcov_type.py b/tests/test_estimators_vcov_type.py index 9ba7709d9..af2980520 100644 --- a/tests/test_estimators_vcov_type.py +++ b/tests/test_estimators_vcov_type.py @@ -39,6 +39,56 @@ def _make_did_panel(n_units: int = 30, seed: int = 20260420) -> pd.DataFrame: return pd.DataFrame(rows) +class TestLeverageOneEstimatorInference: + """Shared fail-closed covariance reaches scalar and event-study results.""" + + @pytest.mark.parametrize("vcov_type", ["hc2", "hc2_bm"]) + @pytest.mark.parametrize("event_study", [False, True]) + def test_single_treated_observation_per_period(self, vcov_type, event_study): + from tests.conftest import assert_nan_inference + + rng = np.random.default_rng(182) + data = pd.DataFrame( + [ + dict(unit=u, time=t, treated=int(u == 0), y=rng.normal() + int(u == 0) * (t + 1)) + for u in range(12) + for t in range(4 if event_study else 2) + ] + ) + kwargs = dict(outcome="y", treatment="treated") + if event_study: + cls = MultiPeriodDiD + kwargs.update(time="time", unit="unit", post_periods=[2, 3], reference_period=1) + else: + cls = DifferenceInDifferences + kwargs.update(post="time") + with warnings.catch_warnings(): + warnings.simplefilter("ignore", FutureWarning) + baseline = cls(vcov_type="hc1").fit(data, **kwargs) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + result = cls(vcov_type=vcov_type).fit(data, **kwargs) + family = "HC2-BM" if vcov_type == "hc2_bm" else "HC2" + # MPD's additional contrast-DOF calculation must not warn a second time. + assert sum(f"{family} variance is undefined" in str(w.message) for w in caught) == 1 + assert np.isnan(result.vcov).all() and np.isnan(result.se) + assert result.att == pytest.approx(baseline.att) + assert_nan_inference( + dict( + se=result.se, t_stat=result.t_stat, p_value=result.p_value, conf_int=result.conf_int + ) + ) + np.testing.assert_allclose(result.residuals, baseline.residuals, atol=1e-13) + np.testing.assert_allclose(result.fitted_values, baseline.fitted_values, atol=1e-13) + if event_study: + for period, effect in result.period_effects.items(): + if period == result.reference_period: + continue + assert effect.effect == pytest.approx(baseline.period_effects[period].effect) + assert np.isnan(effect.se) + assert_nan_inference(vars(effect)) + + # ============================================================================= # robust <-> vcov_type alias resolution # ============================================================================= diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 48888a1d8..2fb62730a 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -2769,13 +2769,12 @@ def test_matches_frozen_dense_oracle(self, kw): assert fin.all(), "oracle produced NaN on a well-conditioned design" np.testing.assert_allclose(dof[fin], oracle[fin], rtol=1e-10) - def test_noise_floor_guard_nans_leverage_one_contrast(self): - """A dummy column firing on exactly one observation gives that row - leverage 1: for the dummy's own coefficient the expanded - denominator's two terms cancel at ~1e20 scale down to the float - noise floor, so the guard must NaN it (the prior dense den > 0 - would have kept the noise and inflated the DOF) while every - ordinary contrast in the same design stays finite.""" + def test_leverage_one_guard_nans_all_contrasts_silently(self): + """A single-observation dummy suppresses all unweighted contrast DOFs. + + The covariance guard owns the warning; this helper stays silent + when computing further contrasts from the same failed design. + """ from diff_diff.linalg import _compute_bm_dof_from_contrasts rng = np.random.default_rng(9) @@ -2785,11 +2784,11 @@ def test_noise_floor_guard_nans_leverage_one_contrast(self): bread = X.T @ X h_diag = np.einsum("ij,ij->i", X @ np.linalg.pinv(bread), X) assert h_diag.max() > 1 - 1e-12 - with warnings.catch_warnings(): - warnings.simplefilter("ignore") + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") dof = _compute_bm_dof_from_contrasts(X, bread, h_diag, np.eye(k)) - assert np.isnan(dof[3]), "leverage-1 dummy coefficient must NaN" - assert np.isfinite(dof[:3]).all(), "ordinary contrasts must stay finite" + assert dof.shape == (4,) and np.isnan(dof).all() + assert not caught class TestCR2BMLowRankAdjustment: diff --git a/tests/test_linalg_hc2_bm.py b/tests/test_linalg_hc2_bm.py index 95f9b0ef8..6490ae689 100644 --- a/tests/test_linalg_hc2_bm.py +++ b/tests/test_linalg_hc2_bm.py @@ -16,24 +16,30 @@ clustered) is now supported via the clubSandwich WLS-CR2 port. Parity against ``clubSandwich::vcovCR(lm(weights=w), type="CR2") + coef_test(test= "Satterthwaite")$df_Satt`` is locked in ``tests/test_methodology_wls_cr2.py``; -this file's tests cover backward-compat (unweighted is bit-equal to prior). +this file covers healthy-design compatibility and the leverage-one NaN policy. """ from __future__ import annotations +import warnings + import numpy as np import pytest from diff_diff.linalg import ( + LinearRegression, + _compute_bm_dof_from_contrasts, _compute_bm_dof_oneway, _compute_cr2_bm, _compute_cr2_bm_contrast_dof, _compute_cr2_bm_vcov_and_dof, _compute_hat_diagonals, + _compute_robust_vcov_numpy, _cr2_adjustment_matrix, compute_robust_vcov, solve_ols, ) +from tests.conftest import assert_nan_inference # ============================================================================= # Fixtures: deterministic OLS datasets with hand-computable properties @@ -66,6 +72,193 @@ def _fit_unweighted(X, y): # ============================================================================= +class TestLeverageOneInference: + """Effective leverage-one observations suppress the whole HC covariance.""" + + @staticmethod + def _design(): + X = np.column_stack([np.ones(4), [0.0, 0.0, 0.0, 1.0]]) + return X, np.array([0.0, 1.0, 2.0, 5.0]) + + @pytest.mark.parametrize("vcov_type", ["hc2", "hc2_bm", "hc3"]) + @pytest.mark.parametrize("return_dof", [False, True]) + @pytest.mark.parametrize("compute", [compute_robust_vcov, _compute_robust_vcov_numpy]) + def test_leverage_one_nan_shapes_and_warning(self, vcov_type, return_dof, compute): + X, _ = self._design() + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + result = compute( + X, + np.array([-1.0, 0.0, 1.0, 0.0]), + vcov_type=vcov_type, + return_dof=return_dof, + ) + vcov = result[0] if return_dof else result + assert vcov.shape == (2, 2) and np.isnan(vcov).all() + if return_dof: + assert result[1].shape == (2,) and np.isnan(result[1]).all() + assert len(caught) == 1 + assert caught[0].category is UserWarning + family = "HC2-BM" if vcov_type == "hc2_bm" else vcov_type.upper() + assert f"{family} variance is undefined: 1 observation(s)" in str(caught[0].message) + assert "Returning NaN vcov" in str(caught[0].message) + + @pytest.mark.parametrize("vcov_type", ["hc2", "hc2_bm", "hc3"]) + @pytest.mark.parametrize("h", [1 - 2e-8, 1 - 1e-8, 1 - 5e-9, 1.0, 1.00001]) + def test_inclusive_threshold_and_over_one(self, monkeypatch, vcov_type, h): + # Drive the exact comparison independently of BLAS rounding at the cutoff. + import diff_diff.linalg as la + + X, _ = self._design() + monkeypatch.setattr( + la, "_compute_hat_diagonals", lambda *a, **k: np.array([h, 0.2, 0.2, 0.2]) + ) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + v = la._compute_robust_vcov_numpy( + X, np.array([1.0, -1.0, 0.0, 0.0]), vcov_type=vcov_type + ) + if h >= 1 - 1e-8: + assert np.isnan(v).all() and len(caught) == 1 + assert "variance is undefined" in str(caught[0].message) + else: + assert np.isfinite(v).all() and not caught + + @pytest.mark.parametrize("vcov_type", ["hc2", "hc2_bm"]) + @pytest.mark.parametrize("rank_reduced", [False, True]) + def test_solve_and_regression_preserve_points(self, vcov_type, rank_reduced): + X, y = self._design() + if rank_reduced: + X = np.column_stack([X, X[:, 1]]) + baseline = solve_ols( + X, y, return_fitted=True, return_vcov=False, rank_deficient_action="silent" + ) + with pytest.warns(UserWarning, match="variance is undefined"): + coef, resid, fitted, vcov = solve_ols( + X, y, vcov_type=vcov_type, return_fitted=True, rank_deficient_action="silent" + ) + for actual, expected in zip((coef, resid, fitted), baseline[:3]): + np.testing.assert_allclose(actual, expected, atol=1e-14, equal_nan=True) + assert np.isnan(vcov).all() + with pytest.warns(UserWarning, match="variance is undefined"): + reg = LinearRegression( + vcov_type=vcov_type, include_intercept=False, rank_deficient_action="silent" + ).fit(X, y) + np.testing.assert_allclose(reg.coefficients_, baseline[0], equal_nan=True) + for j in np.flatnonzero(np.isfinite(coef)): + inference = reg.get_inference(j) + assert np.isfinite(inference.coefficient) and np.isnan(inference.se) + assert_nan_inference(vars(inference)) + + @pytest.mark.parametrize("return_dof", [False, True]) + def test_hc2_bm_all_ones_weights_keep_cr2_boundary(self, return_dof): + X, _ = self._design() + resid = np.array([-1.0, 0.0, 1.0, 0.0]) + with pytest.warns(UserWarning, match="HC2-BM variance is undefined") as caught: + plain = compute_robust_vcov(X, resid, vcov_type="hc2_bm", return_dof=return_dof) + assert len(caught) == 1 + assert np.isnan(plain[0] if return_dof else plain).all() + if return_dof: + assert np.isnan(plain[1]).all() + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + weighted = compute_robust_vcov( + X, + resid, + weights=np.ones(4), + weight_type="pweight", + vcov_type="hc2_bm", + return_dof=return_dof, + ) + dof = _compute_bm_dof_from_contrasts( + X, X.T @ X, np.array([1 / 3, 1 / 3, 1 / 3, 1.0]), np.eye(2), weights=np.ones(4) + ) + assert not caught + expected = np.array([[1.0, -1.0], [-1.0, 1.0]]) / 3 + np.testing.assert_allclose(weighted[0] if return_dof else weighted, expected, atol=1e-12) + np.testing.assert_allclose(dof, [2.0, 2.0], atol=1e-12) + if return_dof: + np.testing.assert_allclose(weighted[1], [2.0, 2.0], atol=1e-12) + + @pytest.mark.parametrize("weight_type", ["pweight", "aweight"]) + @pytest.mark.parametrize("delta", [0.0, 5e-9, 2e-8]) + def test_weighted_hc2_near_one(self, weight_type, delta): + # Weighted bread = 1; the first row's WLS leverage is 1 - delta. + w = np.array([2.0, 3.0]) + X = np.sqrt(np.array([(1 - delta) / w[0], delta / w[1]]))[:, None] + resid = np.array([0.1, -0.2]) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + vcov, dof = compute_robust_vcov( + X, resid, weights=w, weight_type=weight_type, vcov_type="hc2", return_dof=True + ) + if delta < 1e-8: + assert np.isnan(vcov).all() and np.isnan(dof).all() + assert len(caught) == 1 and "HC2 variance is undefined" in str(caught[0].message) + else: + assert not caught and np.isfinite(vcov).all() + np.testing.assert_array_equal(dof, [1.0]) + + @pytest.mark.parametrize("weight_type", ["pweight", "aweight", "fweight"]) + @pytest.mark.parametrize("vcov_type, expected", [("hc2", 1 / 3), ("hc3", 4 / 9)]) + @pytest.mark.parametrize("z", [2.0, 10.0]) + def test_zero_weight_rows_have_no_contribution(self, weight_type, vcov_type, expected, z): + X, y, w = ( + np.array([[1.0], [1.0], [z]]), + np.array([0.0, 2.0, 7.0]), + np.array([2.0, 2.0, 0.0]), + ) + # Includes exact unit leverage and over-one fweight quadratic forms. + resid = y - X[:, 0] + with ( + warnings.catch_warnings(record=True) as caught, + np.errstate(divide="raise", invalid="raise", over="raise"), + ): + warnings.simplefilter("always") + v, df = compute_robust_vcov( + X, resid, weights=w, weight_type=weight_type, vcov_type=vcov_type, return_dof=True + ) + dropped, df_dropped = compute_robust_vcov( + X[:2], + resid[:2], + weights=w[:2], + weight_type=weight_type, + vcov_type=vcov_type, + return_dof=True, + ) + coef, _, v_fit = solve_ols( + X, y, weights=w, weight_type=weight_type, vcov_type=vcov_type + ) + assert not caught and np.isfinite(v).all() + np.testing.assert_allclose(v, dropped, atol=1e-14) + np.testing.assert_allclose(v_fit, dropped, atol=1e-14) + np.testing.assert_allclose(coef, [1.0], atol=1e-14) + np.testing.assert_array_equal(df, df_dropped) + if weight_type == "fweight": + Xe, ye = np.repeat(X, w.astype(int), axis=0), np.repeat(y, w.astype(int)) + _, re, ve = solve_ols(Xe, ye, vcov_type=vcov_type) + _, df_e = compute_robust_vcov(Xe, re, vcov_type=vcov_type, return_dof=True) + np.testing.assert_allclose(v, [[expected]], atol=1e-14) + np.testing.assert_allclose(v, ve, atol=1e-14) + np.testing.assert_array_equal(df, [3.0]) + np.testing.assert_array_equal(df, df_e) + + @pytest.mark.parametrize("vcov_type", ["hc2", "hc3"]) + @pytest.mark.parametrize("count", [1, 2]) + def test_fweight_singleton_vs_repeated(self, vcov_type, count): + X, y = self._design() + w = np.array([2, 2, 2, count]) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + _, _, vc = solve_ols(X, y, weights=w, weight_type="fweight", vcov_type=vcov_type) + _, _, ve = solve_ols(np.repeat(X, w, axis=0), np.repeat(y, w), vcov_type=vcov_type) + if count == 1: + assert np.isnan(vc).all() and np.isnan(ve).all() and len(caught) == 2 + else: + assert not caught and np.isfinite(vc).all() + np.testing.assert_allclose(vc, ve, atol=1e-12) + + class TestClassicalVcov: def test_matches_sigma_squared_inverse_XtX(self, small_ols_dataset): """V = sigma^2 * (X'X)^{-1}.""" diff --git a/tests/test_lwdid.py b/tests/test_lwdid.py index 7bef8a227..9b6fbd051 100644 --- a/tests/test_lwdid.py +++ b/tests/test_lwdid.py @@ -4309,6 +4309,29 @@ class TestReviewRound21Guards: KW = dict(outcome="y", unit="unit", time="time", treatment="treat") + @pytest.mark.parametrize("vcov_type", ["hc2", "hc3"]) + def test_regression_leverage_one_warns_once_and_discards_influence(self, vcov_type): + est = LWDiD(rolling="demean", vcov_type=vcov_type) + treatment = np.array([0.0, 0.0, 0.0, 1.0]) + y = np.array([0.0, 1.0, 2.0, 5.0]) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + att, se, coefs, vcov, n_params, influence = est._estimate_reg( + y, treatment, None, None, 4 + ) + np.testing.assert_allclose(coefs, [1.0, 4.0], atol=1e-14) + assert att == pytest.approx(4.0) and n_params == 2 + assert np.isnan(se) and np.isnan(vcov).all() and influence is None + assert len(caught) == 1 + assert f"{vcov_type.upper()} variance is undefined" in str(caught[0].message) + + # Keep the independent influence safeguard even after removing the + # duplicate covariance guard: aggregations must not reuse finite IFs. + X = np.column_stack([np.ones(4), treatment]) + psi = est._ols_treatment_influence(X, np.linalg.inv(X.T @ X), y - X @ coefs, 4, 2, None) + assert np.isnan(psi).all() + assert est._finalize_influence(psi, se) is None + def test_hc2_leverage_one_fails_closed_on_lwdid(self): rng = np.random.default_rng(0) rows = [] @@ -4317,8 +4340,12 @@ def test_hc2_leverage_one_fails_closed_on_lwdid(self): d = 1 if (u < 1 and t >= 4) else 0 # single treated unit rows.append(dict(unit=u, time=t, treat=d, y=rng.normal() + d)) df = pd.DataFrame(rows) - with pytest.warns(UserWarning, match="HC2 variance is undefined"): + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") res = LWDiD(rolling="demean", vcov_type="hc2").fit(df, **self.KW) + # Three post-period regressions plus the headline regression; each + # emits its own shared covariance diagnostic, with no local duplicate. + assert sum("HC2 variance is undefined" in str(w.message) for w in caught) == 4 assert np.isfinite(res.att) from tests.conftest import assert_nan_inference diff --git a/tests/test_rust_backend.py b/tests/test_rust_backend.py index c8dae731f..60fd10568 100644 --- a/tests/test_rust_backend.py +++ b/tests/test_rust_backend.py @@ -9,6 +9,8 @@ Tests are skipped if the Rust backend is not available. """ +import warnings + import numpy as np import pandas as pd import pytest @@ -3538,14 +3540,77 @@ def test_shape_validation_errors(self): _batched_chol_symbol(np.zeros((2, 3, 3)), np.zeros(5)) +class TestHC2BackendCompatibility: + """Older successful HC2 kernels must not bypass the leverage-one policy.""" + + @pytest.mark.parametrize("backend_mode", ["auto", "rust"]) + @pytest.mark.parametrize("return_dof", [False, True]) + @pytest.mark.parametrize("leverage_one", [False, True]) + def test_legacy_symbol_uses_numpy(self, monkeypatch, backend_mode, return_dof, leverage_one): + import runpy + + import diff_diff._backend as backend + import diff_diff.linalg as la + + native = pytest.importorskip("diff_diff._rust_backend") + X = np.column_stack([np.ones(4), [0.0, 0.0, 0.0 if leverage_one else 1.0, 1.0]]) + residuals = np.array([-1.0, 0.0, 1.0, 0.0]) + calls = [] + + def legacy_hc2(X, residuals): + # Supplied base kernel: unit leverage returns successfully because + # its denominator floor turns 0/0 into a finite contribution. + calls.append(True) + bread_inv = np.linalg.inv(X.T @ X) + h = np.sum((X @ bread_inv) * X, axis=1) + factor = residuals**2 / np.maximum(1.0 - h, 1e-10) + return bread_inv @ (X.T @ (X * factor[:, None])) @ bread_inv + + if leverage_one: + np.testing.assert_allclose( + legacy_hc2(X, residuals), np.array([[1.0, -1.0], [-1.0, 1.0]]) / 3 + ) + calls.clear() + monkeypatch.setattr(native, "compute_robust_vcov_hc2", legacy_hc2) + monkeypatch.delattr(native, "compute_robust_vcov_hc2_v2", raising=False) + monkeypatch.setenv("DIFF_DIFF_BACKEND", backend_mode) + # Execute import selection in an isolated namespace, avoiding reloads + # that change class identities in other modules/tests. + selected = runpy.run_path(backend.__file__) + monkeypatch.setattr(la, "HAS_RUST_BACKEND", selected["HAS_RUST_BACKEND"]) + monkeypatch.setattr( + la, "_rust_compute_robust_vcov_hc2", selected["_rust_compute_robust_vcov_hc2"] + ) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + result = la.compute_robust_vcov(X, residuals, vcov_type="hc2", return_dof=return_dof) + vcov = result[0] if return_dof else result + assert vcov.shape == (2, 2) + if leverage_one: + assert np.isnan(vcov).all() + assert len(caught) == 1 and "HC2 variance is undefined" in str(caught[0].message) + if return_dof: + assert result[1].shape == (2,) and np.isnan(result[1]).all() + else: + assert not caught + np.testing.assert_allclose(vcov, [[0.5, -0.5], [-0.5, 1.0]], atol=1e-14) + if return_dof: + np.testing.assert_array_equal(result[1], [2.0, 2.0]) + assert not calls, "legacy HC2 must never be used, even if its symbol exists" + assert selected["_rust_compute_robust_vcov_hc2"] is None + assert selected["HAS_RUST_BACKEND"] + assert selected["_rust_solve_ols"] is native.solve_ols + assert selected["_rust_compute_robust_vcov"] is native.compute_robust_vcov + + @pytest.mark.skipif(not HAS_RUST_BACKEND, reason="Rust backend not available") class TestRustHC2Vcov: """Rust HC2 (leverage-corrected) vcov parity with the NumPy hc2 branch. The kernel mirrors `_compute_robust_vcov_numpy`'s unweighted `hc2` path - exactly (hat diagonals off the same bread, `u^2 / max(1 - h, 1e-10)` meat, - NO n/(n-k) factor); the near-singular hat-diagonal guard stays Python-side - (sentinel error -> documented warn-and-fall-back-to-HC1).""" + exactly (hat diagonals off the same bread, `u^2 / (1 - h)` meat, + NO n/(n-k) factor). At leverage ~1 the native sentinel becomes a + Python warning and all-NaN covariance.""" @staticmethod def _design(n=400, k=5, seed=0): @@ -3575,10 +3640,8 @@ def test_hc2_kernel_direct(self): np.testing.assert_allclose(v, v.T, rtol=0, atol=1e-12) assert np.all(np.diag(v) > 0) - def test_exact_unit_leverage_clamp_parity(self): - """h_ii == 1 exactly (a one-obs dummy is its own perfect predictor) - does NOT trip the > 1 + 1e-6 guard in either backend — both take the - max(1 - h, 1e-10) clamp path and must agree.""" + def test_exact_unit_leverage_fails_closed_in_both_backends(self): + """A one-observation dummy produces NaN covariance in both backends.""" from diff_diff.linalg import _compute_robust_vcov_numpy, compute_robust_vcov n = 60 @@ -3589,14 +3652,14 @@ def test_exact_unit_leverage_clamp_parity(self): y = X @ np.array([1.0, 2.0, 0.5]) + rng.normal(size=n) resid = y - X @ np.linalg.lstsq(X, y, rcond=None)[0] - v_rust_path = compute_robust_vcov(X, resid, vcov_type="hc2") - v_numpy_path = _compute_robust_vcov_numpy(X, resid, None, vcov_type="hc2") - np.testing.assert_allclose(v_rust_path, v_numpy_path, rtol=1e-9, atol=1e-12) + for compute in (compute_robust_vcov, _compute_robust_vcov_numpy): + with pytest.warns(UserWarning, match="HC2 variance is undefined") as caught: + v = compute(X, resid, vcov_type="hc2") + assert len(caught) == 1 + assert v.shape == (3, 3) and np.isnan(v).all() - def test_sentinel_error_falls_back_to_hc1_with_warning(self, monkeypatch): - """The kernel's near-singular sentinel error must reproduce the NumPy - branch's warn-and-fall-back-to-HC1 through the dispatcher (the guard - decision is Python-side; the kernel only signals).""" + def test_legacy_sentinel_fails_closed_with_warning(self, monkeypatch): + """A legacy over-one sentinel must never return mislabeled HC1.""" import diff_diff.linalg as la X, resid = self._design() @@ -3607,10 +3670,70 @@ def _sentinel(*a, **k): ) monkeypatch.setattr(la, "_rust_compute_robust_vcov_hc2", _sentinel) - with pytest.warns(UserWarning, match="Falling back to HC1"): + with pytest.warns(UserWarning, match="HC2 variance is undefined") as caught: + v = la.compute_robust_vcov(X, resid, vcov_type="hc2") + assert len(caught) == 1 and np.isnan(v).all() + assert "Falling back to HC1" not in str(caught[0].message) + + @pytest.mark.parametrize("symbol", ["compute_robust_vcov_hc2", "compute_robust_vcov_hc2_v2"]) + @pytest.mark.parametrize("delta", [0.0, 5e-9, 2e-8]) + def test_native_threshold_sentinel_and_public_parity(self, delta, symbol): + import diff_diff._rust_backend as rust_backend + + from diff_diff.linalg import _compute_robust_vcov_numpy, compute_robust_vcov + + compute_robust_vcov_hc2 = getattr(rust_backend, symbol) + X = np.sqrt(np.array([[1 - delta], [delta]])) + resid = np.array([0.1, -0.2]) + if delta < 1e-8: + with pytest.raises(ValueError, match="HC2 variance is undefined: 1 observation"): + compute_robust_vcov_hc2(X, resid) + for compute in (compute_robust_vcov, _compute_robust_vcov_numpy): + with pytest.warns(UserWarning, match="HC2 variance is undefined") as caught: + v = compute(X, resid, vcov_type="hc2") + assert len(caught) == 1 and np.isnan(v).all() + else: + native = compute_robust_vcov_hc2(X, resid) + numpy = _compute_robust_vcov_numpy(X, resid, vcov_type="hc2") + assert np.isfinite(native).all() + np.testing.assert_allclose(native, numpy, rtol=1e-8) + + def test_public_dispatch_uses_versioned_kernel(self, monkeypatch): + import diff_diff._rust_backend as native + + import diff_diff.linalg as la + + assert la._rust_compute_robust_vcov_hc2 is native.compute_robust_vcov_hc2_v2 + calls = [] + + def tracked(X, residuals): + calls.append(True) + return native.compute_robust_vcov_hc2_v2(X, residuals) + + monkeypatch.setattr(la, "_rust_compute_robust_vcov_hc2", tracked) + X, residuals = self._design() + vcov = la.compute_robust_vcov(X, residuals, vcov_type="hc2") + assert calls == [True] and np.isfinite(vcov).all() + + @pytest.mark.parametrize("fallback", ["missing", "unstable"]) + def test_leverage_one_fails_closed_on_numpy_fallback(self, monkeypatch, fallback): + import diff_diff.linalg as la + + X = np.column_stack([np.ones(4), [0.0, 0.0, 0.0, 1.0]]) + resid = np.array([-1.0, 0.0, 1.0, 0.0]) + + def unstable(*args, **kwargs): + raise ValueError("Matrix inversion numerically unstable (residual check failed)") + + monkeypatch.setattr( + la, "_rust_compute_robust_vcov_hc2", None if fallback == "missing" else unstable + ) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") v = la.compute_robust_vcov(X, resid, vcov_type="hc2") - v_hc1 = la._compute_robust_vcov_numpy(X, resid, None, vcov_type="hc1") - np.testing.assert_allclose(v, v_hc1, rtol=1e-12, atol=1e-15) + assert np.isnan(v).all() + assert sum("HC2 variance is undefined" in str(w.message) for w in caught) == 1 + assert len(caught) == (1 if fallback == "missing" else 2) def test_dispatch_declined_for_dof_and_weights(self): """return_dof / weights / cluster requests stay on the NumPy path