From 79e096249ead30da04b191f609a620dad071494c Mon Sep 17 00:00:00 2001
From: igerber <isaac.gerber@gmail.com>
Date: Sun, 17 May 2026 21:55:01 -0400
Subject: [PATCH 1/4] linalg: add cluster-aware CR2 Bell-McCaffrey contrast
 DOF; wire MPD avg_att inference
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Closes Gate 6 of the six HC2/HC2-BM NotImplementedError gates:
MultiPeriodDiD(cluster=..., vcov_type="hc2_bm") at estimators.py:1657
previously raised NotImplementedError because _compute_cr2_bm returns
per-coefficient Satterthwaite DOF only — the post-period-average ATT
(`avg_att = (1/n_post) Sum_{t >= t_treat} beta_t`) is a compound
contrast that needed a cluster-aware contrast DOF helper.

New _compute_cr2_bm_contrast_dof in diff_diff/linalg.py generalizes the
per-coefficient loop in _compute_cr2_bm to arbitrary (k, m) contrast
matrices using the identical Pustejovsky-Tipton 2018 Section 4 algebra
(`q = X bread_inv c`, `omega_g = A_g X_g bread_inv c`,
`DOF = trace(B)^2 / trace(B^2)`). _compute_cr2_bm is refactored to
call the new helper via a private _cr2_bm_dof_inner with
`contrasts=eye(k)`; refactor regression at atol=1e-10 confirms the
per-coefficient DOFs are preserved (matmul ordering differs slightly
from the prior inline loop).

MultiPeriodDiD.fit() extends its existing avg_att DOF block (introduced
in PR #459) to branch on effective_cluster_ids: one-way
_compute_bm_dof_from_contrasts when None, cluster-aware
_compute_cr2_bm_contrast_dof otherwise. Cluster IDs are per-observation
length n and are NOT subscripted by the rank-deficient column-drop
mask `_kept` (which indexes coefficients, not observations).

R parity verified at atol=1e-10 against clubSandwich's
Wald_test(constraints=matrix(c, 1), test="HTZ")$df_denom on a new
mpd_clustered_avg_att_dof fixture in
benchmarks/data/clubsandwich_cr2_golden.json. On a 1-row constraint
matrix, HTZ reduces to a Satterthwaite t-test and its df_denom IS the
BM Satterthwaite DOF. The pre-flight smoke test against this same R
target passed at atol=1e-13 before any source edits.

Tests:
- TestCR2BMContrastDOF (4 new tests): refactor regression vs library,
  R-parity for compound contrast, shape validation, cluster-count
  validation.
- test_multi_period_cluster_plus_hc2_bm_rejected flipped to
  test_multi_period_cluster_plus_hc2_bm_produces_finite_inference
  (end-to-end MPD wire-through with finite avg_att / period_effects
  inference assertions).

After this PR, 3 of 6 HC2/HC2-BM gates are lifted (DiD-absorb #458,
MPD-absorb #459, MPD-cluster-contrast-DOF this PR). Remaining: TWFE
absorb (Gate 1), weighted HC2-BM (Gates 4-5).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 CHANGELOG.md                                 |   1 +
 benchmarks/R/generate_clubsandwich_golden.R  |  60 ++++++++
 benchmarks/data/clubsandwich_cr2_golden.json |  20 ++-
 diff_diff/estimators.py                      |  54 ++++---
 diff_diff/linalg.py                          | 146 +++++++++++++++++--
 docs/methodology/REGISTRY.md                 |   2 +-
 tests/test_estimators_vcov_type.py           |  34 +++--
 tests/test_linalg_hc2_bm.py                  | 127 ++++++++++++++++
 8 files changed, 389 insertions(+), 55 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0ec74104..278ce96f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- **`MultiPeriodDiD(cluster=..., vcov_type="hc2_bm")` now supported** (`diff_diff/estimators.py:1657`). Pre-PR the combination raised `NotImplementedError` because the cluster-aware CR2 Bell-McCaffrey Satterthwaite DOF for the post-period-average ATT (`avg_att = (1/n_post) Σ_{t ≥ t_treat} β_t`) was not implemented — only the per-coefficient case existed in `_compute_cr2_bm`. New `_compute_cr2_bm_contrast_dof` helper in `diff_diff/linalg.py` generalizes the per-coefficient loop to arbitrary `(k, m)` contrast matrices using the identical Pustejovsky-Tipton 2018 Section 4 algebra; `_compute_cr2_bm` is refactored to call it with `contrasts=eye(k)` so the existing per-coefficient parity to clubSandwich's `coef_test$df_Satt` is preserved (refactor regression at atol=1e-10). `MultiPeriodDiD.fit()` extends its existing avg_att DOF block to branch on `effective_cluster_ids`: one-way `_compute_bm_dof_from_contrasts` when None, cluster-aware `_compute_cr2_bm_contrast_dof` otherwise. Cluster IDs are per-observation length `n` and are NOT subscripted by the rank-deficient column-drop mask. R parity verified at atol=1e-10 against clubSandwich's `Wald_test(constraints=matrix(c, 1), test="HTZ")$df_denom` on the new `mpd_clustered_avg_att_dof` fixture in `benchmarks/data/clubsandwich_cr2_golden.json` (Wald_test's HTZ on a 1-row constraint matrix yields the Satterthwaite t-test DOF). Per-coefficient `period_effects[t].p_value` / `conf_int` and `avg_att` `avg_p_value` / `avg_conf_int` now reflect the correct Satterthwaite DOF rather than the n-k fallback under cluster+hc2_bm. Weighted CR2-BM (`survey_design=` paths) remains a separate gate. New tests: `tests/test_linalg_hc2_bm.py::TestCR2BMContrastDOF` (4 tests: refactor regression, R-parity, shape validation, cluster-count validation); existing `test_multi_period_cluster_plus_hc2_bm_rejected` flipped to behavioral `test_multi_period_cluster_plus_hc2_bm_produces_finite_inference`.
 - **`MultiPeriodDiD(absorb=..., vcov_type in {"hc2", "hc2_bm"})` now supported** (`diff_diff/estimators.py:1476`). Mirrors the DiD-absorb auto-route shipped earlier in this release: when `absorb=` is paired with `vcov_type in {"hc2","hc2_bm"}`, `MultiPeriodDiD.fit()` promotes the absorb columns to `fixed_effects=` internally so the existing full-dummy-design code path computes the algebraically correct vcov on the event-study design (`treated + period_X dummies + treated:period_X interactions + factor(unit)`). Verified at ~1e-10 vs `lm() + sandwich::vcovHC(type="HC2")` and `lm() + clubSandwich::vcovCR(cluster=1:n, type="CR2")` on a 5-cohort × 5-period event-study fixture (new `tests/test_estimators_vcov_type.py::TestMPDAbsorbedFERParity` against `benchmarks/data/clubsandwich_cr2_golden.json` scenario `mpd_absorbed_fe_did`). HC1/CR1 paths on `absorb=` are unchanged (no leverage term). `TwoWayFixedEffects(vcov_type in {"hc2","hc2_bm"})` rejection remains as a follow-up (different fit-path structure — no `fixed_effects=` equivalent inside TWFE). **Behavioral note (full `MultiPeriodDiDResults` surface change under auto-route):** under the auto-route, the entire returned `MultiPeriodDiDResults` reflects the full-dummy fit rather than the within-transformed fit — `result.coefficients`, `result.vcov`, `result.residuals`, `result.fitted_values`, `result.r_squared` all include the FE-dummy entries / un-demeaned values. `result.period_effects[t].effect` / `.se` / `.p_value` / `.conf_int` and `result.avg_att` / `.avg_se` are invariant to this routing (FWL guarantee). MPD requires a time-invariant ever-treated indicator that lies in the span of the intercept and the post-auto-route unit FE dummies (the exact alias depends on the omitted FE reference category under `pd.get_dummies(drop_first=True)`, not just on "the sum of treated-cohort unit dummies"), so `solve_ols` drops one column from that collinear set under R-style rank-deficiency handling. Which specific column is dropped is pivot-order and dummy-coding dependent (in the shipped parity fixture it is a never-treated unit dummy, not the `treated` main effect itself). The per-period interaction coefficients (`treated:period_X`) and `avg_att` are identified and invariant to that choice; parity tests target those rather than the `treated` main effect. **Survey-design scope (replicate weights):** when `survey_design=` uses replicate weights, the auto-route short-circuits the absorb-refit branch at `estimators.py:1693` and routes through the standard `compute_replicate_vcov` path on the fixed full-dummy design — correct because the design does not depend on replicate weights so no per-replicate refit is needed. **Redundant time-FE skip:** when the routed (or directly-supplied) `fixed_effects` list contains the `time` column, MPD silently skips emitting `<time>_<X>` dummies for that entry because the design already absorbs the time dimension via the non-reference period dummies; without the skip, the two blocks would collide on dummy names and the `coefficients` dict would silently collapse duplicates under `var_names`-keyed construction, breaking the coefficients-vs-vcov alignment that downstream consumers rely on. This applies to both the new `absorb=` auto-route and the pre-existing `fixed_effects=[<time_col>]` invocation.
 - **`DifferenceInDifferences(absorb=..., vcov_type in {"hc2", "hc2_bm"})` now supported** (`diff_diff/estimators.py:382`). Previously raised `NotImplementedError` because the HC2 leverage correction and CR2 Bell-McCaffrey DOF depend on the FULL FE hat matrix, while within-transformation (FWL) preserves coefficients and residuals but not the hat. Lift via internal auto-route: when `absorb=` is paired with `vcov_type in {"hc2","hc2_bm"}`, the fit promotes the absorb columns to `fixed_effects=` internally so the existing full-dummy-design code path computes the algebraically correct vcov. Empirically matches `lm() + sandwich::vcovHC(type="HC2")` and `lm() + clubSandwich::vcovCR(cluster=..., type="CR2")` at ~1e-10 (verified via new `tests/test_estimators_vcov_type.py::TestDiDAbsorbedFERParity` against `benchmarks/data/clubsandwich_cr2_golden.json` scenario `absorbed_fe_did`, with the R generator using the singleton-cluster CR2 trick for one-way HC2-BM Satterthwaite DOF). HC1/CR1 paths unchanged. `MultiPeriodDiD(absorb=...)` and `TwoWayFixedEffects` rejections remain as follow-ups (different fit-path structure). **Behavioral note (full `DiDResults` surface change under auto-route):** under the auto-route, the entire returned `DiDResults` reflects the full-dummy fit rather than the within-transformed fit. Specifically, `result.coefficients` and `result.vcov` include the FE-dummy entries (matching the `fixed_effects=` path), `result.residuals` and `result.fitted_values` are on the un-demeaned outcome scale, and `result.r_squared` is computed on the un-demeaned outcome (so it absorbs the FE variance and will typically be higher than the within-R²). `result.att` is invariant to this routing (FWL guarantee). Downstream consumers reading `result.att` are unaffected; consumers reading the broader result surface should expect the full-dummy values. **Survey-design scope:** the auto-route changes the FE handling (and removes the prior absorbed-FE rejection), but `survey_design=` continues to drive its own variance path (Taylor-series linearization or replicate-weight variance, per the existing survey contract) rather than the analytical HC2/HC2-BM sandwich. The auto-route is therefore methodologically meaningful for non-survey fits and for the FE-handling side of survey fits; analytical small-sample inference under `vcov_type in {"hc2","hc2_bm"}` is bypassed when a survey design is supplied.
 - **`SpilloverDiD` Gardner GMM first-stage uncertainty correction across HC1 / Conley / cluster (Wave D).** Closes the documented Wave B/C "SEs biased downward by a few percent" caveat. **Documented synthesis** of Butts (2021) Section 3.1 (the IF construction for spillover-aware DiD) + Gardner (2022) Section 4 (the two-stage GMM sandwich) + Conley (1999) (the spatial kernel). No reference software combines all three — `did2s` (Butts & Gardner) implements the Gardner correction without rings or Conley; `conleyreg` and `acreg` implement Conley without the two-stage correction. Wave D is the synthesis. Applies unconditionally under `vcov_type ∈ {"hc1", "conley", "cluster"}` for both `event_study=False` AND `event_study=True`. **Formula** (Butts 2021 §3.1 + Gardner 2022 §4): `psi_i = gamma_hat' * X_{10,i} * eps_{10,i} - X_{2,i} * eps_{2,i}` where `gamma_hat = (X_10' X_10)^{-1} (X_1' X_2)` is the stage-1-projection-of-stage-2 cross-moment; meat = `Psi' K Psi` with `K` dispatched by `vcov_type` (identity for HC1, block-indicator for cluster, spatial kernel for Conley); vcov = `(X_2' X_2)^{-1} @ meat @ (X_2' X_2)^{-1}`. **Finite-sample multipliers:** `n/(n-p)` for HC1; `G/(G-1) * (n-1)/(n-p)` for cluster CR1; no multiplier for Conley (preserves `conleyreg` / Wave B convention). **Public surface:** `vcov_type="classical"` now raises `NotImplementedError` upfront (the Wave D synthesis has not been derived for the homoskedastic meat structure `sigma_hat^2 * (X_10' X_10)`); REGISTRY's "vcov_type restrictions" block updated accordingly. **Point estimates unchanged** (`tau_total`, `delta_j`, event-study `tau_k` / `delta_jk` are byte-identical to Wave B/C); SE values shift upward by 1-few percent depending on first-stage residual variance. **Implementation:** new module-level helper `_compute_gmm_corrected_meat` in `diff_diff/two_stage.py` (NOT a modification of the existing `_compute_gmm_variance` method — TwoStageDiD's path is unchanged); new module-level helper `_build_butts_fe_design_csr` in `diff_diff/spillover.py`; new module-level helper `_compute_conley_meat` in `diff_diff/conley.py` factored out of `_compute_conley_vcov` so the same kernel-application code path handles both standard sandwich (`X * residuals`) and Wave D IF outer product (`Psi`) cases. **No new public API kwarg** — the correction is unconditional. Wave D variance mode dispatch derives from the public contract: `vcov_type="conley"` → `"conley"`; `cluster=<col>` → `"cluster"` (CR1); otherwise `"hc1"`. **Wave B/C SE goldens re-pinned** at `tests/test_spillover.py::TestSpilloverDiDEventStudyBackwardCompat` (constants renamed `_WAVE_B_GOLDEN_*` → `_WAVE_D_GOLDEN_*`; pre-Wave-D references retained as commented baselines for the directional inflation invariant `_WAVE_B_UNCORRECTED_*`). **Tests:** new test classes `TestSpilloverDiDWaveDGmmCorrectedHc1Hand` (hand-derived `Psi` on a 4-unit × 3-period over-identified panel — matches at `atol=1e-12`), `TestSpilloverDiDWaveDGmmCorrectedEventStudy` (vcov shape on event-study path), `TestSpilloverDiDWaveDGmmCorrectedNanInferenceContract` (rank-deficient column propagation), `TestSpilloverDiDWaveDGmmCorrectedValidatorWiring` (Conley validator fires from the new helper), `TestSpilloverDiDWaveDGmmCorrectedFitIdempotence` (clone + repeat-fit bit-identity per `feedback_fit_does_not_mutate_config`), `TestSpilloverDiDWaveDPublicVarianceContract` (end-to-end public `cluster=<col>` CR1 routing, single-cluster rejection, classical NotImplementedError). Closes the Gardner-GMM follow-up row in `TODO.md`.
diff --git a/benchmarks/R/generate_clubsandwich_golden.R b/benchmarks/R/generate_clubsandwich_golden.R
index cf11f957..74aaf88d 100644
--- a/benchmarks/R/generate_clubsandwich_golden.R
+++ b/benchmarks/R/generate_clubsandwich_golden.R
@@ -172,6 +172,66 @@ output$mpd_absorbed_fe_did <- list(
   target_period = 4L
 )
 
+# --- MPD clustered avg_att DOF scenario (Gate 6 lift PR) ---------------------
+# Pins clubSandwich's compound-contrast Satterthwaite DOF for the post-period-
+# average ATT under cluster-robust CR2. Mirrors MultiPeriodDiD(cluster=unit,
+# vcov_type='hc2_bm', fixed_effects=['unit']) parameterization. Per-coefficient
+# DOFs use coef_test()$df_Satt (the canonical Satterthwaite per-coef API);
+# the compound contrast DOF uses Wald_test(constraints=matrix(c_avg, 1),
+# test='HTZ')$df_denom — on a 1-row constraint matrix HTZ reduces to a
+# Satterthwaite t-test and its df_denom IS the BM Satterthwaite DOF.
+
+d_mpd_cl <- make_mpd_panel(n_total = 15, units_per_cohort = 5, n_periods = 4,
+                           seed = 20260517)
+d_mpd_cl$period_f <- relevel(factor(d_mpd_cl$period), ref = "1")
+for (p in 2:4) {
+  d_mpd_cl[[paste0("treated_period_", p)]] <-
+    d_mpd_cl$treated * (d_mpd_cl$period == p)
+}
+fit_mpd_cl <- lm(y ~ treated + period_f +
+                     treated_period_2 + treated_period_3 + treated_period_4 +
+                     factor(unit),
+                 data = d_mpd_cl)
+vcov_mpd_cr2 <- vcovCR(fit_mpd_cl, cluster = d_mpd_cl$unit, type = "CR2")
+# Per-coefficient DOF via coef_test (canonical Satterthwaite API).
+ct_mpd_cr2 <- coef_test(fit_mpd_cl, vcov = vcov_mpd_cr2)
+# Compound post-period-average contrast: (1/3) * (e_treated_period_2
+# + e_treated_period_3 + e_treated_period_4). Build full-width vector
+# matching coef(fit) order, with zeros on the NA-dropped column.
+all_coef_names <- names(coef(fit_mpd_cl))
+n_coef <- length(all_coef_names)
+c_avg_vec <- setNames(rep(0, n_coef), all_coef_names)
+post_names <- c("treated_period_2", "treated_period_3", "treated_period_4")
+c_avg_vec[post_names] <- 1 / length(post_names)
+# Wald_test ignores NA-dropped coefficients; subset the constraint vector
+# to the non-NA coefficients (clubSandwich's coef_test convention).
+finite_mask <- !is.na(coef(fit_mpd_cl))
+c_avg_kept <- c_avg_vec[finite_mask]
+dof_avg_compound <- Wald_test(
+  fit_mpd_cl,
+  constraints = matrix(c_avg_kept, 1),
+  vcov = vcov_mpd_cr2,
+  test = "HTZ"
+)$df_denom
+output$mpd_clustered_avg_att_dof <- list(
+  unit = d_mpd_cl$unit,
+  period = d_mpd_cl$period,
+  treated = d_mpd_cl$treated,
+  y = d_mpd_cl$y,
+  cluster = d_mpd_cl$unit,
+  coef = as.numeric(coef(fit_mpd_cl)),
+  coef_names = all_coef_names,
+  finite_coef_names = all_coef_names[finite_mask],
+  vcov_cr2 = as.numeric(vcov_mpd_cr2),
+  vcov_cr2_shape = dim(vcov_mpd_cr2),
+  dof_per_coef = as.numeric(ct_mpd_cr2$df_Satt),
+  c_avg = as.numeric(c_avg_kept),
+  dof_avg = unname(dof_avg_compound),
+  post_interaction_names = post_names,
+  reference_period = 1L,
+  n_post_periods = length(post_names)
+)
+
 output$meta <- list(
   source = "clubSandwich",
   clubSandwich_version = as.character(packageVersion("clubSandwich")),
diff --git a/benchmarks/data/clubsandwich_cr2_golden.json b/benchmarks/data/clubsandwich_cr2_golden.json
index 3524bb56..5e406154 100644
--- a/benchmarks/data/clubsandwich_cr2_golden.json
+++ b/benchmarks/data/clubsandwich_cr2_golden.json
@@ -64,11 +64,29 @@
     "reference_period": 1,
     "target_period": 4
   },
+  "mpd_clustered_avg_att_dof": {
+    "unit": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15],
+    "period": [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4],
+    "treated": [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0],
+    "y": [1.475718242198878, 2.464876166312384, 2.035892793748199, 1.879955566349265, 1.693592098244641, 1.57322216497172, 2.134877806109572, 1.587836699503246, 1.360078460155304, 2.408900333261002, 2.684661494195645, 1.579470084731996, 2.117443467913104, 2.54193910725751, 2.89430848785229, 0.9220357582732823, 2.005370923172815, 1.512446789638281, 3.004528500424255, 1.423132996044842, 1.59463900645967, 1.918906489488893, 2.436742089701388, 2.399198105466515, 3.535946723086255, 2.578767123325587, 2.388928731674242, 2.876514265809677, 3.181117266605847, 3.456092529656238, 1.916608742232821, 1.52189364909609, 2.48698665111912, 2.398279253936628, 3.30520605002529, 2.567770225397498, 3.408057609001849, 3.189083203091819, 2.32492044798654, 4.059123729357513, 3.164819482482473, 3.728160022934698, 3.111440547620898, 2.82187895806224, 2.337998035056887, 2.695918867233574, 3.423567186166589, 2.958073707207942, 1.975027301739997, 3.676766732956131, 2.713119374555802, 2.849552712962999, 3.414065876467119, 3.361369664504964, 4.603756942772487, 2.616067663067715, 3.37096210583077, 3.512532311155356, 3.610259114297078, 2.761275269146486],
+    "cluster": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15],
+    "coef": [2.359251305384248, -1.200612010824764, 0.5327194550242089, 0.6692948808413298, 0.8106547643093722, -0.3189197499340118, 0.1870030421977655, 0.4949720392619664, 0.6013565787023286, 0.495779582943745, 0.561877253127895, 0.7721040668330856, 0.3596172903615322, 0.8252782519061873, 0.9043615647062527, 0.6088212670436904, 1.899361529634674, -0.1013396396601202, -0.09553834413504889, 0.04206406769678345, 0.1763800311276931, "NA"],
+    "coef_names": ["(Intercept)", "treated", "period_f2", "period_f3", "period_f4", "treated_period_2", "treated_period_3", "treated_period_4", "factor(unit)2", "factor(unit)3", "factor(unit)4", "factor(unit)5", "factor(unit)6", "factor(unit)7", "factor(unit)8", "factor(unit)9", "factor(unit)10", "factor(unit)11", "factor(unit)12", "factor(unit)13", "factor(unit)14", "factor(unit)15"],
+    "finite_coef_names": ["(Intercept)", "treated", "period_f2", "period_f3", "period_f4", "treated_period_2", "treated_period_3", "treated_period_4", "factor(unit)2", "factor(unit)3", "factor(unit)4", "factor(unit)5", "factor(unit)6", "factor(unit)7", "factor(unit)8", "factor(unit)9", "factor(unit)10", "factor(unit)11", "factor(unit)12", "factor(unit)13", "factor(unit)14"],
+    "vcov_cr2": [0.05086064825955528, -0.05086064825955527, -0.02549586076543726, -0.0921689198357791, -0.0857778124370035, 0.0254958607654373, 0.09216891983577909, 0.0857778124370035, 4.386695528475936e-18, -2.939065356597548e-17, 2.591167677857621e-17, -1.433211621613874e-17, -2.781166956993463e-17, -8.579767187743784e-18, -3.33815465261822e-17, -1.215212036320773e-17, -1.276317097048087e-17, -2.877341252492203e-16, -2.485216692010425e-16, -3.033992195587622e-16, -3.246643760337161e-16, -0.05086064825955527, 0.06934355909557269, 0.02549586076543727, 0.09216891983577907, 0.08577781243700347, -0.04822228462360257, -0.1213975540484044, -0.1077543977102825, 8.259544084399596e-18, 4.849210570704196e-17, -3.851992325453083e-17, 2.783496561506819e-17, 3.968178977941973e-17, 1.635556349950676e-17, 4.710671821445654e-17, 2.23894478167242e-17, 2.390397546167922e-17, 2.805513295267066e-16, 2.405568882933059e-16, 2.957297564913511e-16, 3.164154269740996e-16, -0.02549586076543727, 0.02549586076543729, 0.02739126676589657, 0.02633084307118829, 0.04826133322466428, -0.02739126676589662, -0.02633084307118825, -0.04826133322466423, 3.723022423470681e-18, 2.482083064319797e-17, -9.338377189254e-19, 1.108997467391252e-17, 1.568022335604176e-17, 1.346023029223708e-17, 2.236974208280944e-17, 1.11689627533574e-17, 1.109842425911693e-17, -6.899726911843268e-17, -4.688491658079492e-17, -3.186960995539225e-17, -1.721474236747706e-17, -0.09216891983577911, 0.09216891983577909, 0.02633084307118829, 0.1991759895536381, 0.143168846718287, -0.02633084307118833, -0.1991759895536382, -0.1431688467182871, -1.800189802116327e-17, 3.676150973763307e-17, -6.873468901384877e-17, 1.767745618238533e-17, 4.499952986243048e-17, -1.512741992162097e-18, 4.856783575727061e-17, 1.193196210459719e-17, 1.360900224000064e-17, 7.240616632590146e-16, 5.790134565062263e-16, 7.170832508008937e-16, 7.474178265343089e-16, -0.08577781243700348, 0.08577781243700346, 0.04826133322466427, 0.1431688467182869, 0.1516810698050605, -0.04826133322466433, -0.1431688467182869, -0.1516810698050605, -3.267906516210872e-18, 5.598027388307047e-17, -3.397818038152988e-17, 2.856103400825675e-17, 5.056692506126546e-17, 2.237158045090007e-17, 6.258860826464791e-17, 2.550755659487607e-17, 2.63452573828056e-17, 4.958721068562704e-16, 4.619581368787103e-16, 5.283832373895193e-16, 5.684544199680039e-16, 0.02549586076543735, -0.04822228462360263, -0.02739126676589666, -0.02633084307118838, -0.04826133322466438, 0.08061526905530911, 0.05119111654854618, 0.06108275289055518, -1.876932914038614e-17, -7.212066135619005e-17, 3.304045680507487e-18, -2.872602411242336e-17, -2.672552438581148e-17, -2.495522895730365e-17, -3.708696650034689e-17, -2.581412401256551e-17, -2.62010240919756e-17, 7.523728869290735e-17, 5.614213479272635e-17, 3.943804054512363e-17, 2.622845213492482e-17, 0.0921689198357791, -0.1213975540484045, -0.02633084307118827, -0.1991759895536382, -0.143168846718287, 0.05119111654854617, 0.261940924135656, 0.1724581755094129, 6.524193797473899e-18, -4.702314056972621e-17, 1.085394652559208e-16, -2.625373233463153e-17, -5.583064887718239e-17, 2.598363170737406e-18, -5.814124996738086e-17, -1.622883480739034e-17, -1.951911965161597e-17, -7.108391552920894e-16, -5.651939112713665e-16, -7.032404215904497e-16, -7.331967120209081e-16, 0.08577781243700351, -0.1077543977102826, -0.04826133322466428, -0.1431688467182871, -0.1516810698050605, 0.06108275289055518, 0.1724581755094128, 0.1974766624411603, -2.079304099468661e-17, -7.482462090225148e-17, 4.223618208169396e-17, -5.636010601321775e-17, -7.617098585468455e-17, -4.306538821146097e-17, -9.319865639009792e-17, -4.751483244694094e-17, -4.98957581031253e-17, -4.866034515076169e-16, -4.531757766945564e-16, -5.191166449200513e-16, -5.586934480103878e-16, 4.38669552847595e-18, 8.259544084399581e-18, 3.723022423470661e-18, -1.800189802116328e-17, -3.267906516210876e-18, -1.876932914038609e-17, 6.524193797473937e-18, -2.079304099468658e-17, 4.611102389459619e-32, 4.925795769504484e-32, 3.940832797822716e-32, 4.795786623002699e-32, 4.241482007939216e-32, 4.711126865847382e-32, 4.761923220132377e-32, 4.532881012113933e-32, 4.612196240060957e-32, -8.15672656418635e-32, -5.485184949950333e-32, -7.29802805679375e-32, -7.265920498714174e-32, -2.93906535659755e-17, 4.8492105707042e-17, 2.482083064319799e-17, 3.676150973763307e-17, 5.598027388307052e-17, -7.212066135619004e-17, -4.702314056972619e-17, -7.482462090225151e-17, 4.92579576950449e-32, 1.032059703677637e-31, 3.535305169328944e-32, 6.540067435499116e-32, 6.280712406235152e-32, 6.107568391212733e-32, 7.531816055772046e-32, 6.140193404365326e-32, 6.216929401478822e-32, 6.210054448107855e-32, 7.507209417776685e-32, 8.78969965966184e-32, 1.029807347900446e-31, 2.591167677857622e-17, -3.851992325453079e-17, -9.338377189253873e-19, -6.87346890138487e-17, -3.397818038152988e-17, 3.304045680507462e-18, 1.085394652559208e-16, 4.223618208169393e-17, 3.940832797822717e-32, 3.535305169328942e-32, 9.116862355733825e-32, 3.433693931441738e-32, 2.078871475837578e-32, 4.764942247302648e-32, 2.605021932800794e-32, 3.785286098653788e-32, 3.661104889921864e-32, -2.235514257696816e-31, -1.512758411008191e-31, -2.118971094920678e-31, -2.167724884391275e-31, -1.433211621613873e-17, 2.783496561506817e-17, 1.108997467391251e-17, 1.767745618238532e-17, 2.856103400825675e-17, -2.872602411242332e-17, -2.625373233463149e-17, -5.636010601321773e-17, 4.795786623002699e-32, 6.540067435499108e-32, 3.433693931441737e-32, 5.919872162253473e-32, 5.707057252829587e-32, 5.500956002745244e-32, 6.53332794545155e-32, 5.451456144462013e-32, 5.558636413016834e-32, 7.06323658164929e-32, 8.037596852913411e-32, 8.007264852769692e-32, 8.806373159154387e-32, -2.781166956993461e-17, 3.968178977941973e-17, 1.568022335604175e-17, 4.499952986243047e-17, 5.056692506126545e-17, -2.672552438581144e-17, -5.583064887718234e-17, -7.617098585468452e-17, 4.241482007939219e-32, 6.280712406235145e-32, 2.078871475837576e-32, 5.707057252829588e-32, 6.215033201744175e-32, 5.126985218480241e-32, 6.889972647047276e-32, 5.262111145549584e-32, 5.387822858883626e-32, 1.794030574264836e-31, 1.731146323360398e-31, 1.887414971130578e-31, 2.020409991829524e-31, -8.579767187743776e-18, 1.635556349950675e-17, 1.346023029223708e-17, -1.512741992162083e-18, 2.237158045090007e-17, -2.495522895730362e-17, 2.598363170737428e-18, -4.306538821146096e-17, 4.711126865847382e-32, 6.107568391212731e-32, 4.764942247302648e-32, 5.500956002745246e-32, 5.126985218480241e-32, 5.637478514137718e-32, 6.035326160005507e-32, 5.265874653358727e-32, 5.316083326933704e-32, -8.603774798987926e-33, 2.620315188664936e-32, 8.311143557930366e-33, 1.635911043341937e-32, -3.338154652618217e-17, 4.710671821445653e-17, 2.236974208280943e-17, 4.85678357572706e-17, 6.258860826464789e-17, -3.708696650034683e-17, -5.814124996738081e-17, -9.319865639009783e-17, 4.761923220132378e-32, 7.53181605577204e-32, 2.605021932800792e-32, 6.53332794545155e-32, 6.889972647047274e-32, 6.035326160005503e-32, 8.041028426946574e-32, 6.037394326604204e-32, 6.170716268733255e-32, 1.745775612634519e-31, 1.791329951038338e-31, 1.923394012465904e-31, 2.095303902121773e-31, -1.215212036320771e-17, 2.238944781672419e-17, 1.116896275335739e-17, 1.193196210459718e-17, 2.550755659487605e-17, -2.581412401256547e-17, -1.62288348073903e-17, -4.751483244694092e-17, 4.532881012113933e-32, 6.140193404365322e-32, 3.785286098653789e-32, 5.451456144462016e-32, 5.262111145549583e-32, 5.265874653358725e-32, 6.037394326604206e-32, 5.14728382147454e-32, 5.209124821915356e-32, 4.592202188435321e-32, 6.139615514128817e-32, 5.690082476249406e-32, 6.450703720837188e-32, -1.276317097048086e-17, 2.390397546167922e-17, 1.109842425911692e-17, 1.360900224000063e-17, 2.63452573828056e-17, -2.620102409197557e-17, -1.951911965161594e-17, -4.989575810312531e-17, 4.612196240060959e-32, 6.216929401478814e-32, 3.661104889921864e-32, 5.558636413016833e-32, 5.387822858883625e-32, 5.316083326933705e-32, 6.170716268733257e-32, 5.209124821915355e-32, 5.327047028001424e-32, 5.34918763738996e-32, 6.724306138512336e-32, 6.395003608538436e-32, 7.166652342207051e-32, -2.877341252492204e-16, 2.805513295267067e-16, -6.899726911843262e-17, 7.240616632590145e-16, 4.958721068562705e-16, 7.523728869290728e-17, -7.108391552920895e-16, -4.866034515076168e-16, -8.156726564186349e-32, 6.210054448107864e-32, -2.235514257696816e-31, 7.063236581649293e-32, 1.794030574264836e-31, -8.603774798987821e-33, 1.74577561263452e-31, 4.592202188435326e-32, 5.349187637389968e-32, 7.294476487756018e-30, 7.165035599654263e-30, 7.155163387206988e-30, 7.176755924763191e-30, -2.485216692010428e-16, 2.405568882933062e-16, -4.688491658079486e-17, 5.790134565062261e-16, 4.619581368787103e-16, 5.614213479272619e-17, -5.651939112713665e-16, -4.531757766945565e-16, -5.485184949950333e-32, 7.507209417776689e-32, -1.512758411008191e-31, 8.037596852913412e-32, 1.731146323360398e-31, 2.620315188664941e-32, 1.791329951038337e-31, 6.13961551412882e-32, 6.72430613851234e-32, 7.165035599654267e-30, 7.410450981384384e-30, 7.159925935035454e-30, 7.178920447369508e-30, -3.033992195587624e-16, 2.957297564913511e-16, -3.186960995539227e-17, 7.170832508008936e-16, 5.283832373895193e-16, 3.943804054512361e-17, -7.032404215904494e-16, -5.191166449200513e-16, -7.29802805679375e-32, 8.789699659661838e-32, -2.118971094920679e-31, 8.007264852769691e-32, 1.887414971130578e-31, 8.311143557930378e-33, 1.923394012465903e-31, 5.690082476249404e-32, 6.395003608538435e-32, 7.155163387206991e-30, 7.159925935035451e-30, 7.239296860586662e-30, 7.171085509022199e-30, -3.246643760337163e-16, 3.164154269740997e-16, -1.721474236747707e-17, 7.474178265343089e-16, 5.684544199680037e-16, 2.62284521349249e-17, -7.331967120209081e-16, -5.586934480103879e-16, -7.265920498714173e-32, 1.029807347900446e-31, -2.167724884391277e-31, 8.806373159154382e-32, 2.020409991829523e-31, 1.635911043341936e-32, 2.095303902121772e-31, 6.450703720837188e-32, 7.166652342207051e-32, 7.176755924763191e-30, 7.178920447369511e-30, 7.171085509022198e-30, 7.195798625784396e-30],
+    "vcov_cr2_shape": [21, 21],
+    "dof_per_coef": [3.999999999999998, 8.099999999999973, 3.999999999999999, 3.999999999999995, 4, 8.10000000000003, 8.099999999999975, 8.100000000000019, 7.234852539498682, 12.987409715937888, 7.107153046689751, 9.689867008603835, 8.713557378534674, 9.869789870457211, 11.658110074475884, 8.317990581609115, 8.225205519715464, 1.036361342793117, 1.030992816554161, 1.058519135166968, 1.058900897899314],
+    "c_avg": [0, 0, 0, 0, 0, 0.3333333333333333, 0.3333333333333333, 0.3333333333333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    "dof_avg": 8.100000000000016,
+    "post_interaction_names": ["treated_period_2", "treated_period_3", "treated_period_4"],
+    "reference_period": 1,
+    "n_post_periods": 3
+  },
   "meta": {
     "source": "clubSandwich",
     "clubSandwich_version": "0.7.0",
     "R_version": "R version 4.5.2 (2025-10-31)",
-    "generated_at": "2026-05-17 10:36:19 UTC",
+    "generated_at": "2026-05-18 01:50:55 UTC",
     "note": "CR2 Bell-McCaffrey cluster-robust parity target for diff_diff._compute_cr2_bm"
   }
 }
diff --git a/diff_diff/estimators.py b/diff_diff/estimators.py
index 9d73e968..771a4438 100644
--- a/diff_diff/estimators.py
+++ b/diff_diff/estimators.py
@@ -1647,27 +1647,6 @@ def fit(  # type: ignore[override]
         # Determine if survey vcov should be used
         _use_survey_vcov = resolved_survey is not None and resolved_survey.needs_survey_vcov
 
-        # Reject cluster + vcov_type="hc2_bm": `_compute_cr2_bm` produces CR2
-        # per-coefficient DOF, but the post-period-average contrast needs a
-        # cluster-aware contrast-BM DOF that isn't implemented yet. Pairing
-        # CR2 SEs with one-way BM DOF would be a broken hybrid — reject with
-        # a clear error until the cluster-aware contrast DOF is in place.
-        # Tracked in TODO.md. Users can drop cluster for one-way HC2+BM, or
-        # drop vcov_type for CR1 cluster-robust.
-        if (
-            self.vcov_type == "hc2_bm"
-            and effective_cluster_ids is not None
-            and not _use_survey_vcov
-        ):
-            raise NotImplementedError(
-                "MultiPeriodDiD(cluster=..., vcov_type='hc2_bm') is not yet "
-                "supported: the cluster-aware CR2 Bell-McCaffrey contrast DOF "
-                "for the post-period average has not been implemented. "
-                "Workarounds: use vcov_type='hc2_bm' without cluster (one-way "
-                "HC2 + BM DOF), or use vcov_type='hc1' with cluster (CR1 "
-                "Liang-Zeger cluster-robust)."
-            )
-
         # Remap implicit "classical" + cluster to CR1 (legacy backward compat).
         _fit_vcov_type = self._resolve_effective_vcov_type(effective_cluster_ids)
 
@@ -1870,6 +1849,7 @@ def _refit_mp_absorb(w_r):
         ):
             from diff_diff.linalg import (
                 _compute_bm_dof_from_contrasts,
+                _compute_cr2_bm_contrast_dof,
                 _compute_hat_diagonals,
             )
 
@@ -1880,7 +1860,6 @@ def _refit_mp_absorb(w_r):
                 bread_kept = X_kept.T @ (
                     X_kept * survey_weights[:, np.newaxis] if survey_weights is not None else X_kept
                 )
-                h_diag_kept = _compute_hat_diagonals(X_kept, bread_kept, weights=survey_weights)
                 # Build the contrast matrix: one column per identified coefficient
                 # plus one column for the post-period average contrast (1/n_post
                 # on each post-period interaction column, 0 elsewhere).
@@ -1893,13 +1872,30 @@ def _refit_mp_absorb(w_r):
                         post_contrast_full[interaction_indices[_p]] = 1.0 / _n_post
                 post_contrast_kept = post_contrast_full[_kept]
                 contrasts = np.column_stack([np.eye(n_kept), post_contrast_kept[:, np.newaxis]])
-                _dof_all = _compute_bm_dof_from_contrasts(
-                    X_kept,
-                    bread_kept,
-                    h_diag_kept,
-                    contrasts,
-                    weights=survey_weights,
-                )
+                # Branch on cluster: one-way HC2-BM vs cluster-aware CR2-BM.
+                # Cluster IDs are per-observation length n and are unchanged
+                # by the column-drop applied to X (`_kept` indexes columns
+                # only); pass `effective_cluster_ids` unmodified.
+                if effective_cluster_ids is None:
+                    h_diag_kept = _compute_hat_diagonals(X_kept, bread_kept, weights=survey_weights)
+                    _dof_all = _compute_bm_dof_from_contrasts(
+                        X_kept,
+                        bread_kept,
+                        h_diag_kept,
+                        contrasts,
+                        weights=survey_weights,
+                    )
+                else:
+                    # Cluster-aware CR2 BM Satterthwaite DOF for per-coefficient
+                    # AND post-period-average compound contrast (Gate 6 lift).
+                    # Weighted CR2-BM is a separate gate; survey paths never
+                    # reach this block (outer `not _use_survey_vcov` guard).
+                    _dof_all = _compute_cr2_bm_contrast_dof(
+                        X_kept,
+                        effective_cluster_ids,
+                        bread_kept,
+                        contrasts,
+                    )
                 # Expand per-coefficient DOF back to full width (NaN for dropped).
                 _bm_dof_per_coef = np.full(X.shape[1], np.nan)
                 _bm_dof_per_coef[_kept] = _dof_all[:n_kept]
diff --git a/diff_diff/linalg.py b/diff_diff/linalg.py
index 793e1987..dc4f7db6 100644
--- a/diff_diff/linalg.py
+++ b/diff_diff/linalg.py
@@ -1591,21 +1591,66 @@ def _compute_cr2_bm(
     vcov = np.linalg.solve(bread_matrix, temp.T).T
 
     # --- Per-coefficient Bell-McCaffrey cluster DOF ---
-    # omega_g(c) = A_g @ X_g @ bread_inv @ c  (length n_g)
-    # trace(B) = sum_i (X_i' bread_inv c)^2
-    # trace(B^2) = sum_{g, h} (omega_g' M_{g, h} omega_h)^2
-    dof_vec = np.empty(k)
-    # Precompute X bread_inv (n x k) so contrast-specific q = X_bi[:, j].
-    X_bi = X @ bread_inv
-    # Precompute A_g @ X_g @ bread_inv per cluster (A_g_X_bi shape n_g x k)
-    A_g_Xbi = {g: A_g_matrices[g] @ X[cluster_idx[g]] @ bread_inv for g in unique_clusters}
-    for j in range(k):
-        q = X_bi[:, j]  # length n
+    # Delegate to the contrast-aware helper with `contrasts=I_k` so the
+    # per-coefficient case is `c = e_j` (the j-th basis vector). Bit-identity
+    # vs the prior inline loop holds at machine precision because the same
+    # X_bi / A_g_Xbi precomputes are reused under the same matmul ordering.
+    dof_vec = _cr2_bm_dof_inner(X, M, A_g_matrices, cluster_idx, bread_inv, np.eye(k))
+
+    return vcov, dof_vec
+
+
+def _cr2_bm_dof_inner(
+    X: np.ndarray,
+    M: np.ndarray,
+    A_g_matrices: Dict[Any, np.ndarray],
+    cluster_idx: Dict[Any, np.ndarray],
+    bread_inv: np.ndarray,
+    contrasts: np.ndarray,
+) -> np.ndarray:
+    """Inner DOF loop, parameterized by an arbitrary contrast matrix.
+
+    Computes the CR2 Bell-McCaffrey Satterthwaite DOF for each column of
+    ``contrasts`` (shape ``(k, m)``), using the precomputed residual-maker
+    ``M``, per-cluster adjustment matrices ``A_g_matrices``, cluster index
+    map ``cluster_idx``, and ``bread_inv``. The per-coefficient case is
+    recovered with ``contrasts=np.eye(k)``; compound contrasts (e.g., a
+    post-period-average ATT) are handled by the same algebra without
+    duplication.
+
+    Per-contrast formula (Pustejovsky-Tipton 2018 Section 4 / Appendix A):
+
+      q       = X @ bread_inv @ c                       (length n)
+      omega_g = A_g @ X_g @ bread_inv @ c               (length n_g)
+      trace_B = sum_i q_i**2
+      trace_B2 = sum_{g, h} (omega_g' M_{g, h} omega_h)**2
+      DOF(c)  = trace_B**2 / trace_B2
+
+    Returns
+    -------
+    dof_vec : ndarray of shape (m,)
+        DOF per contrast column. NaN entries indicate degenerate contrasts
+        (trace_B2 ≈ 0 — typically high-collinearity nuisance columns).
+    """
+    m = contrasts.shape[1]
+    unique_clusters = list(cluster_idx.keys())
+    # Precompute once: q-matrix (n, m) and A_g_Xbi (n_g, k) per cluster.
+    # For unit-contrast inputs (contrasts=I_k), this matches the prior
+    # inline implementation exactly: q[:, j] == X_bi[:, j] == X @ bread_inv @ e_j.
+    X_bi = X @ bread_inv  # (n, k)
+    Q = X_bi @ contrasts  # (n, m) — q vectors as columns
+    A_g_Xbi = {
+        g: A_g_matrices[g] @ X[cluster_idx[g]] @ bread_inv for g in unique_clusters
+    }  # each (n_g, k)
+    # Omega per cluster per contrast: (n_g, m) = A_g_Xbi[g] @ contrasts
+    omega_all = {g: A_g_Xbi[g] @ contrasts for g in unique_clusters}
+
+    dof_vec = np.empty(m)
+    for j in range(m):
+        q = Q[:, j]
         trace_B = float(np.sum(q * q))
-        # trace(B^2) = sum_{g, h} (omega_g' M_{g, h} omega_h)^2
         trace_B2 = 0.0
-        # Cache omega_g for this contrast
-        omega_cache = {g: A_g_Xbi[g][:, j] for g in unique_clusters}
+        omega_cache = {g: omega_all[g][:, j] for g in unique_clusters}
         for g in unique_clusters:
             idx_g = cluster_idx[g]
             omega_g = omega_cache[g]
@@ -1617,7 +1662,80 @@ def _compute_cr2_bm(
                 trace_B2 += val * val
         dof_vec[j] = (trace_B * trace_B) / trace_B2 if trace_B2 > 0 else np.nan
 
-    return vcov, dof_vec
+    return dof_vec
+
+
+def _compute_cr2_bm_contrast_dof(
+    X: np.ndarray,
+    cluster_ids: np.ndarray,
+    bread_matrix: np.ndarray,
+    contrasts: np.ndarray,
+) -> np.ndarray:
+    """Per-contrast CR2 Bell-McCaffrey Satterthwaite DOF.
+
+    Generalizes the per-coefficient DOF from :func:`_compute_cr2_bm` to
+    arbitrary linear combinations ``c = sum_j a_j * beta_j``. Used by
+    :class:`MultiPeriodDiD` to compute the Satterthwaite DOF for the
+    post-period-average ATT contrast under cluster-robust CR2 inference.
+
+    Parameters
+    ----------
+    X : ndarray of shape (n, k)
+        Design matrix (post-rank-deficient-column-drop if applicable).
+    cluster_ids : ndarray of shape (n,)
+        Per-observation cluster identifiers. NOT subscripted by any
+        column mask — cluster IDs are unchanged by column drops.
+    bread_matrix : ndarray of shape (k, k)
+        ``X.T @ X`` (or ``X.T @ W @ X`` for survey-weighted — though the
+        weighted CR2 path is deferred to a follow-up PR).
+    contrasts : ndarray of shape (k, m)
+        Each column is a contrast vector ``c`` for the linear combination
+        ``c' beta``. The per-coefficient case is recovered with
+        ``contrasts=np.eye(k)``.
+
+    Returns
+    -------
+    dof_vec : ndarray of shape (m,)
+        Satterthwaite DOF per contrast.
+
+    See Also
+    --------
+    _compute_cr2_bm : per-coefficient DOF (calls this helper internally
+        with ``contrasts=np.eye(k)``).
+    """
+    n, k = X.shape
+    cluster_ids_arr = np.asarray(cluster_ids)
+    unique_clusters = np.unique(cluster_ids_arr)
+    if len(unique_clusters) < 2:
+        raise ValueError(
+            f"Need at least 2 clusters for cluster-robust SEs, got " f"{len(unique_clusters)}"
+        )
+    if contrasts.shape[0] != k:
+        raise ValueError(
+            f"contrasts must have {k} rows (one per coefficient), got " f"{contrasts.shape[0]}"
+        )
+
+    try:
+        bread_inv = np.linalg.solve(bread_matrix, np.eye(k))
+    except np.linalg.LinAlgError as e:
+        if "Singular" in str(e):
+            raise ValueError(
+                "Design matrix is rank-deficient (singular X'X matrix). "
+                "Cannot compute CR2 Bell-McCaffrey variance."
+            ) from e
+        raise
+
+    H = X @ bread_inv @ X.T
+    M = np.eye(n) - H
+    cluster_idx = {g: np.where(cluster_ids_arr == g)[0] for g in unique_clusters}
+    A_g_matrices: Dict[Any, np.ndarray] = {}
+    for g in unique_clusters:
+        idx_g = cluster_idx[g]
+        H_gg = H[np.ix_(idx_g, idx_g)]
+        I_g = np.eye(len(idx_g))
+        A_g_matrices[g] = _cr2_adjustment_matrix(I_g - H_gg)
+
+    return _cr2_bm_dof_inner(X, M, A_g_matrices, cluster_idx, bread_inv, contrasts)
 
 
 def _compute_bm_dof_from_contrasts(
diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md
index c4b49d00..04096f80 100644
--- a/docs/methodology/REGISTRY.md
+++ b/docs/methodology/REGISTRY.md
@@ -2554,7 +2554,7 @@ Shipped in `diff_diff/had_pretests.py` as `stute_joint_pretest()` (residuals-in
         - **`MultiPeriodDiD(absorb=..., vcov_type in {"hc2","hc2_bm"})` — SUPPORTED (auto-route).** Same auto-route pattern as `DifferenceInDifferences`: `MultiPeriodDiD.fit()` internally promotes the absorb columns to `fixed_effects=` for HC2 / HC2-BM callers, so the existing full-dummy code path computes the algebraically correct vcov from the full FE projection on the event-study design (`treated + period_X dummies + treated:period_X interactions + factor(unit)`). Verified at ~1e-10 vs `lm() + sandwich::vcovHC(type="HC2")` and `lm() + clubSandwich::vcovCR(cluster=1:n, type="CR2")` on a 5-cohort × 5-period event-study fixture; the parity target is a per-period interaction `treated:period_X` because MPD requires the `treated` column to be a time-invariant ever-treated indicator, which lies in the span of the intercept and the post-auto-route unit FE dummies (under `pd.get_dummies(..., drop_first=True)` the dropped reference unit is implicit in the intercept, so the exact alias relation depends on the omitted FE category — it is NOT simply "the sum of treated-cohort unit dummies"). `solve_ols` drops one column from the collinear set under R-style rank-deficiency handling; in the shipped parity fixture (4 ever-treated cohorts of 5 units + 1 never-treated cohort of 5 units) it drops a unit dummy from the never-treated cohort (`unit_25`) and the `treated` main effect remains finite, but the specific column that gets NaN'd is pivot-order and dummy-coding dependent. Either way, the slope coefficients (`treated:period_X`) and the post-period-average `avg_att` are identified and invariant to which column was dropped. Same `MultiPeriodDiDResults` surface change as DiD: `vcov`, `residuals`, `fitted_values`, `r_squared`, and `coefficients` reflect the full-dummy fit, with `period_effects[t].effect` / `.se` / `.p_value` / `.conf_int` invariant by FWL. HC1/CR1 paths on `absorb=` are unchanged (no leverage term). Same survey-design scope as DiD: replicate-weight variance routes through the standard `compute_replicate_vcov` path on the fixed full-dummy design rather than the per-replicate refit branch (which targets the demeaning path); since the auto-routed design does not depend on replicate weights, no refit is needed. **Redundant time-FE skip:** when the routed (or directly-supplied) `fixed_effects` list contains the `time` column, MPD silently skips emitting `<time>_<X>` dummies for that entry because the design already absorbs time via the non-reference period dummies. Without the skip, those blocks collide on dummy names and `MultiPeriodDiDResults.coefficients` (built as `{name: coef for name, coef in zip(var_names, coefficients)}`) would silently drop duplicates, breaking the coefficients-vs-vcov alignment that downstream consumers (HonestDiD sub-VCV extraction, BusinessReport, etc.) rely on. The skip applies to BOTH the new `absorb=` auto-route AND the pre-existing `fixed_effects=[<time_col>]` invocation (pre-PR, `fixed_effects=["unit", time]` produced a dict with `len < vcov.shape[0]` and NaN values overwriting the real event-study period coefficients).
         - Workarounds for the still-rejecting paths: use `vcov_type="hc1"` (HC1/CR1 have no leverage term and survive FWL), or switch to `fixed_effects=` dummies so the hat matrix is computed on the full design.
 - [x] Phase 1a: `vcov_type` enum threaded through `DifferenceInDifferences` (`MultiPeriodDiD`, `TwoWayFixedEffects` inherit); `robust=True` <=> `vcov_type="hc1"`, `robust=False` <=> `vcov_type="classical"`. Conflict detection at `__init__`. Results summary prints the variance-family label.
-    - **Note (deviation from the fully-symmetric enum):** `MultiPeriodDiD(cluster=..., vcov_type="hc2_bm")` is intentionally **not supported** and raises `NotImplementedError`. The scalar-coefficient `DifferenceInDifferences` path handles the cluster + CR2 Bell-McCaffrey combination (`_compute_cr2_bm` returns a per-coefficient Satterthwaite DOF that is valid for the single-ATT contrast), but `MultiPeriodDiD` also reports a post-period-average ATT constructed as a *contrast* of the event-study coefficients. The cluster-aware CR2 BM DOF for that contrast (i.e., the Pustejovsky-Tipton 2018 per-cluster adjustment matrices applied to an arbitrary aggregation contrast) is not yet implemented. Pairing CR2 cluster-robust SEs with the one-way Imbens-Kolesar (2016) contrast DOF would be a broken hybrid, so the combination fails fast with a clear workaround message (drop the cluster for one-way HC2+BM, or use `vcov_type="hc1"` with cluster for CR1 Liang-Zeger). Tracked in `TODO.md` under Methodology/Correctness. Applies only to `MultiPeriodDiD`; `DifferenceInDifferences(cluster=..., vcov_type="hc2_bm")` works.
+    - **Note (`MultiPeriodDiD(cluster=..., vcov_type="hc2_bm")` — SUPPORTED via cluster-aware contrast DOF):** the scalar-coefficient `DifferenceInDifferences` path uses `_compute_cr2_bm`'s per-coefficient Satterthwaite DOF directly for the single-ATT contrast, but `MultiPeriodDiD` also reports a post-period-average ATT constructed as a *contrast* of the event-study coefficients (`avg_att = (1/n_post) Σ_{t ≥ t_treat} β_t`). Pre-PR the combination raised `NotImplementedError` because the cluster-aware CR2 Bell-McCaffrey Satterthwaite DOF for an arbitrary linear combination was not implemented — only the per-coefficient case existed. The new `_compute_cr2_bm_contrast_dof` helper in `diff_diff/linalg.py` generalizes the per-coefficient loop to arbitrary `(k, m)` contrast matrices using the identical Pustejovsky-Tipton 2018 Section 4 algebra (`q = X bread_inv c`, `omega_g = A_g X_g bread_inv c`, `DOF = trace(B)² / trace(B²)`), and `_compute_cr2_bm` is refactored to call it with `contrasts=eye(k)` so the per-coefficient case is recovered at machine precision (atol=1e-10, see refactor regression in `tests/test_linalg_hc2_bm.py::TestCR2BMContrastDOF`). `MultiPeriodDiD.fit()` extends its existing avg_att DOF block to branch on cluster presence: one-way `_compute_bm_dof_from_contrasts` for `effective_cluster_ids is None`, cluster-aware `_compute_cr2_bm_contrast_dof` otherwise. R parity verified against clubSandwich's `Wald_test(constraints=matrix(c, 1), test="HTZ")$df_denom` at atol=1e-10 on the `mpd_clustered_avg_att_dof` fixture in `benchmarks/data/clubsandwich_cr2_golden.json` (Wald_test's HTZ on a 1-row constraint matrix yields the Satterthwaite t-test DOF). Cluster IDs are per-observation length `n` and are NOT subscripted by the rank-deficient column-drop mask `_kept` — the helper accepts the full `effective_cluster_ids` array. Weighted CR2-BM (`survey_design=` paths) remains a separate gate.
 - [x] Phase 1a: `clubSandwich::vcovCR(..., type="CR2")` parity harness committed: R script at `benchmarks/R/generate_clubsandwich_golden.R` plus the authoritative R-generated JSON at `benchmarks/data/clubsandwich_cr2_golden.json` (`"source": "clubSandwich"`, with `clubSandwich_version`, `R_version`, and `generated_at` captured in `meta` for forensic traceability). The parity test at `tests/test_linalg_hc2_bm.py::TestCR2BMCluster::test_cr2_parity_with_golden` runs at 1e-6 tolerance and passes at ≤ 7.1e-15 across all three datasets — Python's `_compute_cr2_bm` matches clubSandwich at machine precision.
 - [x] Phase 1b: Calonico-Cattaneo-Farrell (2018) MSE-optimal bandwidth selector. In-house port of `nprobust::lpbwselect(bwselect="mse-dpi")` (nprobust 0.5.0, SHA `36e4e53`) as `diff_diff.mse_optimal_bandwidth` and `BandwidthResult`, backed by the private `diff_diff._nprobust_port` module (`kernel_W`, `lprobust_bw`, `lpbwselect_mse_dpi`). Three-stage DPI with four `lprobust.bw` calls at orders `q+1`, `q+2`, `q`, `p`. Python matches R to `0.0000%` relative error (i.e., bit-parity within float64 precision, ~8-13 digits agreement) on all five stage bandwidths (`c_bw`, `bw_mp2`, `bw_mp3`, `b_mse`, `h_mse`) across three deterministic DGPs (uniform, Beta(2,2), half-normal) via `benchmarks/R/generate_nprobust_golden.R` → `benchmarks/data/nprobust_mse_dpi_golden.json`. **Note:** `weights=` is currently unsupported (raises `NotImplementedError`); nprobust's `lpbwselect` has no weight argument so there is no parity anchor. Weighted-data support deferred to Phase 2 (survey-design adaptation). **Note (public API scope restriction):** the exported wrapper `mse_optimal_bandwidth` hard-codes the HAD Phase 1b configuration (`p=1`, `deriv=0`, `interior=False`, `vce="nn"`, `nnmatch=3`). The underlying port supports a broader surface (`hc0`/`hc1`/`hc2`/`hc3` variance, interior evaluation, higher `p`), but those paths are not parity-tested against `nprobust` and are deferred. Callers needing the broader surface should use `diff_diff._nprobust_port.lpbwselect_mse_dpi` directly and accept that parity has not been verified on non-HAD configurations. **Note (input contract):** the wrapper enforces HAD's support restriction `D_{g,2} >= 0` (front-door `ValueError` on negative doses and empty inputs). `boundary` must equal `0` (Design 1') or `float(d.min())` (Design 1 continuous-near-d_lower) within float tolerance; off-support values raise `ValueError`. When `boundary ~ 0`, the wrapper additionally requires `d.min() <= 0.05 * median(|d|)` as a Design 1' support plausibility heuristic, chosen to pass the paper's thin-boundary-density DGPs (Beta(2,2), d.min/median ~ 3%) while rejecting substantially off-support samples (U(0.5, 1.0), d.min/median ~ 1.0). Detected mass-point designs (`d.min() > 0` with modal fraction at `d.min() > 2%`) raise `NotImplementedError` pointing to the Phase 2 2SLS path per paper Section 3.2.4.
 - [x] Phase 1c: First-order bias estimator `M̂_{ĥ*_G}` and robust variance `V̂_{ĥ*_G}`. Implemented via Calonico-Cattaneo-Titiunik (2014) bias-combined design matrix `Q.q` in the in-house port `diff_diff._nprobust_port.lprobust` (single-eval-point path of `nprobust::lprobust`, npfunctions.R:177-246).
diff --git a/tests/test_estimators_vcov_type.py b/tests/test_estimators_vcov_type.py
index 8ea85b62..2609402b 100644
--- a/tests/test_estimators_vcov_type.py
+++ b/tests/test_estimators_vcov_type.py
@@ -507,13 +507,15 @@ def test_multi_period_fit_honors_classical(self):
         # SEs must differ — vcov_type actually changed the variance family.
         assert r_hc1.avg_se != pytest.approx(r_classical.avg_se, abs=1e-10)
 
-    def test_multi_period_cluster_plus_hc2_bm_rejected(self):
-        """MultiPeriodDiD rejects cluster + hc2_bm until contrast-aware cluster BM lands.
-
-        The CR2 per-coefficient DOF is available, but the post-period-average
-        contrast DOF under cluster-robust Bell-McCaffrey is not yet
-        implemented. Pairing CR2 SEs with one-way BM DOF would be a broken
-        hybrid. Fail fast with a clear workaround.
+    def test_multi_period_cluster_plus_hc2_bm_produces_finite_inference(self):
+        """MultiPeriodDiD(cluster=..., vcov_type='hc2_bm') is now supported.
+
+        The cluster-aware CR2 Bell-McCaffrey contrast DOF for the
+        post-period-average ATT is implemented via the new
+        `_compute_cr2_bm_contrast_dof` helper that generalizes the
+        per-coefficient loop in `_compute_cr2_bm` to arbitrary linear
+        combinations of coefficients. End-to-end smoke test: assert
+        finite avg_att inference (was `NotImplementedError` pre-PR).
         """
         rng = np.random.default_rng(2)
         rows = []
@@ -524,9 +526,21 @@ def test_multi_period_cluster_plus_hc2_bm_rejected(self):
                 rows.append({"unit": i, "time": t, "treated": treated, "y": y})
         data = pd.DataFrame(rows)
 
-        est = MultiPeriodDiD(vcov_type="hc2_bm", cluster="unit")
-        with pytest.raises(NotImplementedError, match="cluster"):
-            est.fit(data, outcome="y", treatment="treated", time="time")
+        res = MultiPeriodDiD(vcov_type="hc2_bm", cluster="unit").fit(
+            data, outcome="y", treatment="treated", time="time"
+        )
+        # Headline contract: finite avg_att inference under cluster+hc2_bm.
+        assert np.isfinite(res.avg_att)
+        assert np.isfinite(res.avg_se)
+        assert np.isfinite(res.avg_t_stat)
+        assert np.isfinite(res.avg_p_value)
+        assert np.isfinite(res.avg_conf_int[0])
+        assert np.isfinite(res.avg_conf_int[1])
+        # Per-period inference should also be finite for the post period.
+        post_pe = res.period_effects[2]
+        assert np.isfinite(post_pe.effect)
+        assert np.isfinite(post_pe.se)
+        assert np.isfinite(post_pe.p_value)
 
     def test_multi_period_fit_honors_hc2_bm(self):
         """MultiPeriodDiD.fit with vcov_type='hc2_bm' uses Bell-McCaffrey DOF.
diff --git a/tests/test_linalg_hc2_bm.py b/tests/test_linalg_hc2_bm.py
index 4d82c545..d4044d8e 100644
--- a/tests/test_linalg_hc2_bm.py
+++ b/tests/test_linalg_hc2_bm.py
@@ -18,6 +18,7 @@
 from diff_diff.linalg import (
     _compute_bm_dof_oneway,
     _compute_cr2_bm,
+    _compute_cr2_bm_contrast_dof,
     _compute_hat_diagonals,
     _cr2_adjustment_matrix,
     compute_robust_vcov,
@@ -616,3 +617,129 @@ def test_hc2_pweight_matches_manual(self, small_ols_dataset):
 
         got = compute_robust_vcov(X, resid, vcov_type="hc2", weights=w, weight_type="pweight")
         np.testing.assert_allclose(got, expected, atol=1e-10)
+
+
+# =============================================================================
+# Cluster-aware CR2 BM contrast-DOF helper (Gate 6 lift)
+# =============================================================================
+
+
+class TestCR2BMContrastDOF:
+    """Tests for `_compute_cr2_bm_contrast_dof`.
+
+    The helper generalizes the per-coefficient Satterthwaite DOF in
+    `_compute_cr2_bm` to arbitrary linear combinations of coefficients
+    (used by `MultiPeriodDiD` to compute the cluster-aware DOF for the
+    post-period-average ATT contrast).
+    """
+
+    def _load_golden_scenario(self):
+        """Load `mpd_clustered_avg_att_dof` scenario from R generator."""
+        import json
+        from pathlib import Path
+
+        golden_path = (
+            Path(__file__).parent.parent / "benchmarks" / "data" / "clubsandwich_cr2_golden.json"
+        )
+        if not golden_path.exists():
+            pytest.skip(
+                "Golden JSON not present; run "
+                "`Rscript benchmarks/R/generate_clubsandwich_golden.R` first."
+            )
+        with open(golden_path) as f:
+            golden = json.load(f)
+        if "mpd_clustered_avg_att_dof" not in golden:
+            pytest.skip(
+                "Golden JSON does not include `mpd_clustered_avg_att_dof` "
+                "scenario; regenerate via the R script."
+            )
+        return golden["mpd_clustered_avg_att_dof"]
+
+    def _build_mpd_design(self, d):
+        """Construct the MPD-style design matrix that mirrors R's lm()
+        formula `treated + period_f + treated_period_X (non-ref) + factor(unit)`."""
+        unit = np.array(d["unit"])
+        period = np.array(d["period"])
+        treated = np.array(d["treated"], dtype=float)
+        n = len(period)
+        n_periods = int(period.max())
+        non_ref = list(range(2, n_periods + 1))
+        const = np.ones(n)
+        period_dummies = np.column_stack([(period == p).astype(float) for p in non_ref])
+        interaction_dummies = np.column_stack(
+            [(treated * (period == p)).astype(float) for p in non_ref]
+        )
+        n_units = int(unit.max())
+        unit_dummies = np.column_stack([(unit == u).astype(float) for u in range(2, n_units + 1)])
+        X_full = np.column_stack(
+            [const, treated, period_dummies, interaction_dummies, unit_dummies]
+        )
+        # Drop the last unit dummy to match R's rank-deficient drop on this
+        # parameterization (never-treated cohort's last unit is collinear with
+        # the intercept + treated + remaining unit dummies).
+        return X_full[:, :-1]
+
+    def test_unit_contrasts_match_compute_cr2_bm(self):
+        """Refactor anchor: calling the helper with `contrasts=eye(k)`
+        produces the same per-coefficient DOFs as `_compute_cr2_bm`.
+
+        Matmul ordering differs (helper applies eye separately, library
+        slices precomputed columns), so use atol=1e-10 not bit-identity.
+        """
+        d = self._load_golden_scenario()
+        X = self._build_mpd_design(d)
+        k = X.shape[1]
+        y = np.array(d["y"], dtype=float)
+        cluster = np.array(d["cluster"])
+        bread = X.T @ X
+        coef, *_ = np.linalg.lstsq(X, y, rcond=None)
+        residuals = y - X @ coef
+
+        _, dof_lib = _compute_cr2_bm(X, residuals, cluster, bread)
+        dof_helper = _compute_cr2_bm_contrast_dof(X, cluster, bread, np.eye(k))
+
+        finite_both = np.isfinite(dof_lib) & np.isfinite(dof_helper)
+        assert finite_both.any(), "expected at least one finite DOF"
+        np.testing.assert_allclose(dof_helper[finite_both], dof_lib[finite_both], atol=1e-10)
+
+    def test_compound_contrast_matches_clubsandwich(self):
+        """R-parity anchor: compound post-period-average contrast DOF
+        matches clubSandwich's `Wald_test(test="HTZ")$df_denom` at 1e-10.
+        """
+        d = self._load_golden_scenario()
+        X = self._build_mpd_design(d)
+        k = X.shape[1]
+        cluster = np.array(d["cluster"])
+        bread = X.T @ X
+
+        # The R golden stores c_avg as a (k_finite,) vector aligned with
+        # finite_coef_names. Our X already has the rank-deficient column
+        # dropped, so c_avg aligns directly.
+        c_avg = np.array(d["c_avg"])
+        assert c_avg.shape == (k,), f"c_avg shape {c_avg.shape} does not match X.shape[1] {k}"
+
+        dof_avg_py = float(_compute_cr2_bm_contrast_dof(X, cluster, bread, c_avg[:, np.newaxis])[0])
+        dof_avg_r = float(d["dof_avg"])
+        np.testing.assert_allclose(dof_avg_py, dof_avg_r, atol=1e-10)
+
+    def test_invalid_contrast_shape_raises(self):
+        """Helper validates that the contrast matrix's row count matches `k`."""
+        rng = np.random.default_rng(20260517)
+        n, k = 30, 3
+        X = rng.standard_normal((n, k))
+        cluster = np.repeat(np.arange(5), 6)
+        bread = X.T @ X
+        bad_contrasts = np.zeros((k + 1, 1))
+        with pytest.raises(ValueError, match=r"rows"):
+            _compute_cr2_bm_contrast_dof(X, cluster, bread, bad_contrasts)
+
+    def test_too_few_clusters_raises(self):
+        """Helper requires at least 2 clusters (matching `_compute_cr2_bm`)."""
+        rng = np.random.default_rng(20260517)
+        n, k = 30, 3
+        X = rng.standard_normal((n, k))
+        # Everyone in cluster 1 -> only 1 unique cluster.
+        cluster = np.ones(n, dtype=int)
+        bread = X.T @ X
+        with pytest.raises(ValueError, match=r"[Nn]eed at least 2 clusters"):
+            _compute_cr2_bm_contrast_dof(X, cluster, bread, np.eye(k))

From 96ce1c54c223f24af2a7b832e1e4ca348b8a16d0 Mon Sep 17 00:00:00 2001
From: igerber <isaac.gerber@gmail.com>
Date: Sun, 17 May 2026 22:03:43 -0400
Subject: [PATCH 2/4] R1 local Codex polish: ndim validation + docstrings +
 R-parity DOF pin
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Local Codex review on commit 79e09624 returned ✅ with 3 P3s (all
documentation/coverage, no actionable P0/P1). Per the test-coverage P3
upgrade rule (feedback_test_coverage_gap_treat_as_actionable.md),
addressing all three:

P3 #1 (Code Quality): `_compute_cr2_bm_contrast_dof` was missing the
`ndim` validation that the parallel one-way `_compute_bm_dof_from_contrasts`
helper has, so a stray `(k,)` 1-D vector would die with a low-level
indexing error instead of a contract error. Added the same shape-tuple
check pattern (`if contrasts.ndim != 2 or contrasts.shape[0] != k`).

P3 #2 (Docs): two stale doc surfaces post-feature-lift —
  - `estimators.py:68-71` base estimator docstring still said MPD did
    NOT support cluster + hc2_bm. Rewrote to describe the new
    cluster-aware contrast-DOF support and flag survey CR2-BM as the
    remaining gate.
  - `tests/test_linalg_hc2_bm.py` module banner still said clustered
    CR2 BM was "deferred to a follow-up". Updated to describe both the
    per-coefficient and the new compound-contrast DOF surfaces, and
    narrow the deferral to the weighted CR2-BM case only.

P3 #3 (Tests): the new MPD test only asserted finite output, so a
regression that silently fell back to the shared n-k DOF would still
pass. Added `test_multi_period_cluster_hc2_bm_avg_att_uses_clubsandwich_dof`
which fits MPD on the new R `mpd_clustered_avg_att_dof` fixture and
recovers the implied Satterthwaite DOF by inverting
`avg_p_value = 2 * (1 - t.cdf(|avg_t_stat|, df))` via scipy.brentq. The
recovered DOF must match the R `Wald_test(test="HTZ")$df_denom` at
atol=1e-6. Also pins that the implied DOF is much smaller than the
n-k fallback (~39 here) — catches a regression to the shared df path.

All 254 tests in tests/test_linalg_hc2_bm.py + test_estimators_vcov_type.py
+ test_estimators.py pass; lint clean.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 diff_diff/estimators.py            |  7 ++-
 diff_diff/linalg.py                |  6 +--
 tests/test_estimators_vcov_type.py | 75 ++++++++++++++++++++++++++++++
 tests/test_linalg_hc2_bm.py        | 13 ++++--
 4 files changed, 92 insertions(+), 9 deletions(-)

diff --git a/diff_diff/estimators.py b/diff_diff/estimators.py
index 771a4438..a10d4e48 100644
--- a/diff_diff/estimators.py
+++ b/diff_diff/estimators.py
@@ -67,8 +67,11 @@ class DifferenceInDifferences:
           ``cluster=``; use ``"hc2_bm"`` for clustered Bell-McCaffrey.
         - ``"hc2_bm"``: one-way HC2 + Imbens-Kolesar (2016) Satterthwaite DOF;
           with ``cluster=``, Pustejovsky-Tipton (2018) CR2 cluster-robust.
-          (Note: ``MultiPeriodDiD`` does NOT yet support ``cluster=`` with
-          ``"hc2_bm"`` — see ``MultiPeriodDiD`` docstring and REGISTRY.md.)
+          ``MultiPeriodDiD(cluster=..., vcov_type="hc2_bm")`` is supported and
+          uses a cluster-aware Bell-McCaffrey contrast DOF for the
+          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.
         - ``"conley"``: Conley 1999 spatial-HAC sandwich. Pass
           ``conley_coords=(lat_col, lon_col)``, ``conley_cutoff_km=<float>``,
           and ``conley_lag_cutoff=<int>`` on the constructor; pass
diff --git a/diff_diff/linalg.py b/diff_diff/linalg.py
index dc4f7db6..d9d9f652 100644
--- a/diff_diff/linalg.py
+++ b/diff_diff/linalg.py
@@ -1710,10 +1710,8 @@ def _compute_cr2_bm_contrast_dof(
         raise ValueError(
             f"Need at least 2 clusters for cluster-robust SEs, got " f"{len(unique_clusters)}"
         )
-    if contrasts.shape[0] != k:
-        raise ValueError(
-            f"contrasts must have {k} rows (one per coefficient), got " f"{contrasts.shape[0]}"
-        )
+    if contrasts.ndim != 2 or contrasts.shape[0] != k:
+        raise ValueError(f"contrasts must have shape (k={k}, m); got {contrasts.shape}")
 
     try:
         bread_inv = np.linalg.solve(bread_matrix, np.eye(k))
diff --git a/tests/test_estimators_vcov_type.py b/tests/test_estimators_vcov_type.py
index 2609402b..cf1db28e 100644
--- a/tests/test_estimators_vcov_type.py
+++ b/tests/test_estimators_vcov_type.py
@@ -542,6 +542,81 @@ def test_multi_period_cluster_plus_hc2_bm_produces_finite_inference(self):
         assert np.isfinite(post_pe.se)
         assert np.isfinite(post_pe.p_value)
 
+    def test_multi_period_cluster_hc2_bm_avg_att_uses_clubsandwich_dof(self):
+        """MPD(cluster=..., hc2_bm) `avg_att` inference uses the new
+        cluster-aware contrast Satterthwaite DOF, not the shared n-k fallback.
+
+        Pins the implied DOF from `avg_p_value` against the R `dof_avg` target
+        on the `mpd_clustered_avg_att_dof` fixture. The R-side compound contrast
+        DOF (from `Wald_test(test="HTZ")$df_denom` on a 1-row constraint matrix
+        — equivalent to the Satterthwaite t-test DOF) is the parity target.
+        Recovers the DOF by inverting `avg_p_value = 2 * (1 - t.cdf(|t|, df))`.
+        """
+        import json
+        from pathlib import Path
+
+        from scipy import stats
+
+        golden_path = (
+            Path(__file__).parent.parent / "benchmarks" / "data" / "clubsandwich_cr2_golden.json"
+        )
+        if not golden_path.exists():
+            pytest.skip(
+                "Golden JSON not present; run "
+                "`Rscript benchmarks/R/generate_clubsandwich_golden.R` first."
+            )
+        with open(golden_path) as f:
+            golden = json.load(f)
+        if "mpd_clustered_avg_att_dof" not in golden:
+            pytest.skip("Golden JSON missing `mpd_clustered_avg_att_dof` scenario.")
+        d = golden["mpd_clustered_avg_att_dof"]
+        dof_avg_r = float(d["dof_avg"])
+
+        data = pd.DataFrame(
+            {
+                "unit": d["unit"],
+                "period": d["period"],
+                "treated": d["treated"],
+                "y": d["y"],
+            }
+        )
+        # MPD parameterization with `fixed_effects=["unit"]` matches the R
+        # generator's `factor(unit)` term (the cluster column is the unit).
+        res = MultiPeriodDiD(vcov_type="hc2_bm", cluster="unit").fit(
+            data,
+            outcome="y",
+            treatment="treated",
+            time="period",
+            fixed_effects=["unit"],
+            reference_period=int(d["reference_period"]),
+            unit="unit",
+        )
+        assert np.isfinite(res.avg_att) and np.isfinite(res.avg_se)
+        # Recover the implied DOF from the reported p_value:
+        # avg_p_value = 2 * (1 - t.cdf(|t|, df))  ->  df = root of
+        # `t.sf(|t|, df) * 2 - p` (Satterthwaite-bounded scalar bisection
+        # via scipy's brentq on a sane interval).
+        t_stat = abs(res.avg_t_stat)
+        p_target = res.avg_p_value
+        # Sanity: BM Satterthwaite DOF is bounded in [1, n-k]. With 60 obs and
+        # ~21 coefficients (per the fixture), DOF is in [1, 39].
+        from scipy.optimize import brentq
+
+        def _residual(df):
+            return 2.0 * stats.t.sf(t_stat, df) - p_target
+
+        implied_dof = brentq(_residual, 1.0, 100.0, xtol=1e-8)
+        # The implied DOF should match the R golden at 1e-6 (small tolerance
+        # accounts for the t.cdf evaluation roundoff, not the DOF computation).
+        np.testing.assert_allclose(implied_dof, dof_avg_r, atol=1e-6)
+        # Pin that the new path is in use, not the n-k fallback: dof_avg_r
+        # is well below n-k for this fixture (60 obs - 21 coefs = 39 > 8.1).
+        assert implied_dof < 30, (
+            f"Implied DOF {implied_dof:.2f} is suspiciously large; expected "
+            f"~{dof_avg_r:.2f} (Satterthwaite-corrected) and the n-k fallback "
+            "would be ~39, so the contrast-DOF helper may not be wired."
+        )
+
     def test_multi_period_fit_honors_hc2_bm(self):
         """MultiPeriodDiD.fit with vcov_type='hc2_bm' uses Bell-McCaffrey DOF.
 
diff --git a/tests/test_linalg_hc2_bm.py b/tests/test_linalg_hc2_bm.py
index d4044d8e..b9619c61 100644
--- a/tests/test_linalg_hc2_bm.py
+++ b/tests/test_linalg_hc2_bm.py
@@ -6,8 +6,15 @@
   ``robust=False`` on ``DifferenceInDifferences``).
 - ``vcov_type="hc2"``: leverage-corrected HC2 one-way.
 - ``vcov_type="hc2_bm"``: HC2 plus Imbens-Kolesar (2016) Satterthwaite DOF.
-
-Cluster-robust CR2 Bell-McCaffrey is deferred to a follow-up Phase 1a commit.
+- ``vcov_type="hc2_bm"`` + ``cluster=``: Pustejovsky-Tipton (2018) CR2 cluster-
+  robust with per-coefficient and (via ``_compute_cr2_bm_contrast_dof``) compound
+  contrast Bell-McCaffrey Satterthwaite DOF. The contrast-DOF helper is the
+  backend for ``MultiPeriodDiD``'s post-period-average ATT inference under the
+  cluster+hc2_bm combination.
+
+Weighted CR2 Bell-McCaffrey (``hc2_bm`` + ``cluster=`` + ``weights=``) remains
+deferred to a follow-up; the corresponding ``NotImplementedError`` gate in
+``_validate_vcov_args`` is exercised by ``TestInvalidInputs``.
 """
 
 from __future__ import annotations
@@ -730,7 +737,7 @@ def test_invalid_contrast_shape_raises(self):
         cluster = np.repeat(np.arange(5), 6)
         bread = X.T @ X
         bad_contrasts = np.zeros((k + 1, 1))
-        with pytest.raises(ValueError, match=r"rows"):
+        with pytest.raises(ValueError, match=r"shape \(k="):
             _compute_cr2_bm_contrast_dof(X, cluster, bread, bad_contrasts)
 
     def test_too_few_clusters_raises(self):

From 41a323e501eddb818c53b4d183785eeb66a5ea34 Mon Sep 17 00:00:00 2001
From: igerber <isaac.gerber@gmail.com>
Date: Sun, 17 May 2026 22:38:56 -0400
Subject: [PATCH 3/4] R2 local Codex polish: pin the SAME compound contrast as
 the R fixture
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Local Codex R2 returned ✅ with 1 substantive P3: the new MPD parity
test was fitting MPD without an explicit `post_periods=` argument, so
MPD's default "last half of periods" rule selected `[3, 4]` on the
4-period fixture while the R generator defines the parity contrast
over `[2, 3, 4]` (per `post_interaction_names` in the JSON). The
Satterthwaite DOFs happened to coincide here (~8.1) which masked the
estimand mismatch — the test would have stayed green if MPD silently
fit the wrong contrast.

Fix: derive `post_periods` from the golden JSON's
`post_interaction_names` field and pass it explicitly to MPD.fit().
The test now asserts that MPD computes `avg_att` over the exact same
contrast vector R uses for the Wald_test DOF target.

This makes the test a genuine estimand-level parity pin rather than
just a DOF-magnitude smell check.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 tests/test_estimators_vcov_type.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/tests/test_estimators_vcov_type.py b/tests/test_estimators_vcov_type.py
index cf1db28e..691b051b 100644
--- a/tests/test_estimators_vcov_type.py
+++ b/tests/test_estimators_vcov_type.py
@@ -582,6 +582,15 @@ def test_multi_period_cluster_hc2_bm_avg_att_uses_clubsandwich_dof(self):
         )
         # MPD parameterization with `fixed_effects=["unit"]` matches the R
         # generator's `factor(unit)` term (the cluster column is the unit).
+        # Derive `post_periods` from the R `post_interaction_names` so the
+        # contrast we compare to is the SAME `c_avg = (1/n_post) Σ e_p` that
+        # the R generator builds. Without this, MPD defaults to the
+        # last-half-of-periods rule and computes avg_att over [3, 4] on this
+        # 4-period panel, but the R fixture's `c_avg` is over [2, 3, 4] —
+        # the DOFs happen to coincide here but the avg_att estimands differ.
+        post_periods = [
+            int(name.rsplit("_", 1)[1]) for name in d["post_interaction_names"]
+        ]
         res = MultiPeriodDiD(vcov_type="hc2_bm", cluster="unit").fit(
             data,
             outcome="y",
@@ -589,6 +598,7 @@ def test_multi_period_cluster_hc2_bm_avg_att_uses_clubsandwich_dof(self):
             time="period",
             fixed_effects=["unit"],
             reference_period=int(d["reference_period"]),
+            post_periods=post_periods,
             unit="unit",
         )
         assert np.isfinite(res.avg_att) and np.isfinite(res.avg_se)

From 59971295ca0b6c794c0d746b48932980e70fafbd Mon Sep 17 00:00:00 2001
From: igerber <isaac.gerber@gmail.com>
Date: Mon, 18 May 2026 06:05:10 -0400
Subject: [PATCH 4/4] R1 CI Codex polish: doc-consistency + perf-debt TODO row
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

CI Codex review on PR #465 (commit 41a323e5) returned ✅ with two
findings:

P2 (Documentation): Two stale doc surfaces post-feature-lift still said
MPD does NOT support cluster + hc2_bm. The HC2/BM scope-limitation Note
at REGISTRY.md:2557 was already updated in 79e09624, but:
- REGISTRY.md:167-176 (main MPD section under HeterogeneousAdoptionDiD
  requirements-checklist) still had the old "not supported" Note.
- estimators.py MPD class docstring's `cluster` and `vcov_type` blocks
  still said `cluster + hc2_bm` raises NotImplementedError.

Both rewritten to describe the new supported path with a pointer to
_compute_cr2_bm_contrast_dof and the clubSandwich Wald_test(HTZ)
parity target. Weighted CR2-BM remains the only documented gate.

P3 (Performance): _compute_cr2_bm_contrast_dof recomputes H, M, and
per-cluster A_g matrices that solve_ols → _compute_cr2_bm already
built for the vcov path. O(n²k) redundant work per clustered hc2_bm
MPD fit; acceptable for typical cluster-robust DiD panel sizes
(n ≤ few thousand). Tracked as a new Performance row in TODO.md;
acknowledged with a Note in REGISTRY.md per the codex deferral rules
(`**Note:**` label + TODO entry downgrades to P3-deferred).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 TODO.md                      |  1 +
 diff_diff/estimators.py      | 22 +++++++++++-----------
 docs/methodology/REGISTRY.md | 25 +++++++++++++++++--------
 3 files changed, 29 insertions(+), 19 deletions(-)

diff --git a/TODO.md b/TODO.md
index ae7550e1..9aa28973 100644
--- a/TODO.md
+++ b/TODO.md
@@ -148,6 +148,7 @@ Deferred items from PR reviews that were not addressed before merge.
 | Rust faer SVD ndarray-to-faer conversion overhead (minimal vs SVD cost) | `rust/src/linalg.rs:67` | #115 | Low |
 | Unrelated label events (e.g., adding `bug` label) re-trigger CI workflows when `ready-for-ci` is already present; filter `labeled`/`unlabeled` events to only `ready-for-ci` transitions | `.github/workflows/rust-test.yml`, `notebooks.yml`, `docs-tests.yml` | #269 | Low |
 | `bread_inv` as a performance kwarg on `compute_robust_vcov` to avoid re-inverting `(X'WX)` when the caller already has it. Deferred from Phase 1a for scope. HC2 and HC2+BM both need the bread inverse, so a shared hint would save one `np.linalg.solve` per sandwich. | `linalg.py::compute_robust_vcov` | Phase 1a | Low |
+| MPD cluster+hc2_bm path computes CR2 precomputes twice — once via `solve_ols` → `_compute_cr2_bm` for vcov + per-coefficient DOF, then again via `_compute_cr2_bm_contrast_dof` from `MultiPeriodDiD.fit()` for the post-period-average contrast DOF. Both rebuild `H = X bread_inv X'`, the residual-maker `M`, and the per-cluster `A_g = (I - H_gg)^{-1/2}` matrices. O(n²k) redundant work; acceptable for typical cluster-robust DiD panel sizes (n ≤ a few thousand). Fix would plumb the contrast DOF through the existing CR2 vcov path (intrusive API change) or share the precomputes via a cached helper. | `linalg.py::_compute_cr2_bm_contrast_dof`, `estimators.py::MultiPeriodDiD.fit` | follow-up | Low |
 | Rust-backend HC2 implementation. Current Rust path only supports HC1; HC2 and CR2 Bell-McCaffrey fall through to the NumPy backend. For large-n fits this is noticeable. | `rust/src/linalg.rs` | Phase 1a | Low |
 | CR2 Bell-McCaffrey DOF uses a naive `O(n² k)` per-coefficient loop over cluster pairs. Pustejovsky-Tipton (2018) Appendix B has a scores-based formulation that avoids the full `n × n` `M` matrix. Switch when a user hits a large-`n` cluster-robust design. | `linalg.py::_compute_cr2_bm` | Phase 1a | Low |
 
diff --git a/diff_diff/estimators.py b/diff_diff/estimators.py
index a10d4e48..e764d3a8 100644
--- a/diff_diff/estimators.py
+++ b/diff_diff/estimators.py
@@ -1108,16 +1108,13 @@ class MultiPeriodDiD(DifferenceInDifferences):
         contradictory (e.g. ``robust=False, vcov_type="hc2"`` raises).
     cluster : str, optional
         Column name for cluster-robust standard errors. With ``vcov_type="hc1"``
-        dispatches to CR1 (Liang-Zeger).
-
-        **Not supported with** ``vcov_type="hc2_bm"``: the cluster-aware CR2
-        Bell-McCaffrey contrast DOF for the post-period-average ATT is not
-        yet implemented, and pairing CR2 SEs with one-way Imbens-Kolesar DOF
-        would be a broken hybrid, so the combination raises
-        ``NotImplementedError`` with a pointer to workarounds. Tracked in
-        ``TODO.md``; also documented as a Note in
-        ``docs/methodology/REGISTRY.md`` under the HeterogeneousAdoptionDiD
-        requirements-checklist block.
+        dispatches to CR1 (Liang-Zeger). With ``vcov_type="hc2_bm"`` dispatches
+        to CR2 cluster-robust SEs with Bell-McCaffrey Satterthwaite DOF on both
+        per-period coefficients and the post-period-average ATT contrast (the
+        latter via the new ``_compute_cr2_bm_contrast_dof`` helper in
+        ``linalg.py``; matches clubSandwich's
+        ``Wald_test(test="HTZ")$df_denom`` at atol=1e-10). Weighted CR2-BM
+        (``survey_design=``) is a separate, still-gated path.
     vcov_type : {"classical", "hc1", "hc2", "hc2_bm", "conley"}, optional
         Variance-covariance family. Defaults to the ``robust`` alias.
 
@@ -1128,7 +1125,10 @@ class MultiPeriodDiD(DifferenceInDifferences):
           ``cluster=``; use ``"hc2_bm"`` without cluster for Bell-McCaffrey.
         - ``"hc2_bm"``: one-way HC2 + Imbens-Kolesar (2016) Satterthwaite DOF
           per coefficient plus a contrast-aware DOF for the post-period-average
-          ATT. **Unsupported with** ``cluster=`` — see ``cluster`` above.
+          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.
         - ``"conley"``: Conley 1999 spatial-HAC sandwich via the panel
           block-decomposed form (matches R ``conleyreg`` with
           ``lag_cutoff > 0``). Pass ``conley_coords=(lat_col, lon_col)``,
diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md
index 04096f80..195601ee 100644
--- a/docs/methodology/REGISTRY.md
+++ b/docs/methodology/REGISTRY.md
@@ -166,14 +166,23 @@ where V is the VCV sub-matrix for post-treatment δ_e coefficients.
 - Alternative: Cluster-robust at unit level via `cluster` parameter (recommended for panel data)
 - `vcov_type="hc2_bm"` (one-way) computes HC2 + Imbens-Kolesar (2016) Satterthwaite DOF
   per coefficient and a contrast-aware DOF for the post-period-average ATT.
-- **Note:** `cluster` + `vcov_type="hc2_bm"` is **not supported** and raises
-  `NotImplementedError`. The cluster-aware CR2 Bell-McCaffrey contrast DOF for the
-  post-period-average ATT (Pustejovsky-Tipton 2018 per-cluster adjustment matrices
-  applied to an arbitrary aggregation contrast) is not yet implemented. Pairing CR2
-  cluster-robust SEs with the one-way Imbens-Kolesar contrast DOF would be a broken
-  hybrid, so the combination fails fast. Workarounds: drop `cluster` for one-way
-  HC2+BM, or keep `cluster` with the default `vcov_type="hc1"` for CR1 (Liang-Zeger).
-  Tracked in `TODO.md` under Methodology/Correctness.
+- `cluster` + `vcov_type="hc2_bm"` is now supported (PR for Gate 6 lift). Both per-period
+  effects and the post-period-average ATT use a cluster-aware Bell-McCaffrey
+  Satterthwaite DOF: the per-coefficient case continues via `_compute_cr2_bm`, and
+  the compound `avg_att` contrast DOF goes through the new
+  `_compute_cr2_bm_contrast_dof` helper in `diff_diff/linalg.py` (Pustejovsky-Tipton
+  2018 §4 algebra generalized to arbitrary `(k, m)` contrast matrices). R parity
+  verified at atol=1e-10 vs clubSandwich's
+  `Wald_test(constraints=matrix(c, 1), test="HTZ")$df_denom`. Weighted CR2-BM
+  (`survey_design=` paths) is still gated separately; see the rows in `TODO.md`
+  under Methodology/Correctness.
+- **Note:** the cluster-aware contrast-DOF path currently recomputes the CR2 hat
+  matrix and per-cluster adjustment matrices that `solve_ols` already built for the
+  vcov dispatch — clustered `hc2_bm` MPD fits pay the O(n²) CR2 setup twice in
+  exchange for a clean call-site contract. Acceptable for typical cluster-robust
+  DiD panel sizes (n ≤ few thousand); tracked in `TODO.md` under Performance for
+  a follow-up that plumbs the contrast DOF through the existing CR2 vcov path or
+  shares precomputes.
 - Optional: Wild cluster bootstrap (complex for multi-coefficient testing;
   requires joint bootstrap distribution)
 - Degrees of freedom adjusted for absorbed fixed effects