Close SDID placebo R-parity gap: warm-start + R-anchored fixture + test seam

Closes the last queued SDID R-parity follow-up (per
``project_sdid_pr349_followups.md``, now removed from the memory
index — the work is shipped). Symmetric with the existing
``test_jackknife_se_matches_r`` anchor in TestJackknifeSERParity.

Methodology fix — placebo warm-start:
``synthdid:::placebo_se`` (R/vcov.R) seeds Frank-Wolfe per draw with
``weights.boot$omega = sum_normalize(weights$omega[ind[1:N0_placebo]])``
(fit-time ω subsetted + renormalized) and the fit-time
``weights$lambda``, then re-estimates with ``update.omega=TRUE,
update.lambda=TRUE``. Python's ``_placebo_variance_se`` previously
used uniform cold-start, producing finite-iter convergence-pattern
drift on a handful of draws relative to R's reference SE on the
same panel.

Fix: add ``init_omega`` and ``init_lambda`` kwargs to
``_placebo_variance_se``. The dispatcher now passes ``init_omega=
unit_weights, init_lambda=time_weights`` (fit-time outputs); the
loop seeds ``compute_sdid_unit_weights(init_weights=
_sum_normalize(init_omega[pseudo_control_idx]))`` and
``compute_time_weights(init_weights=init_lambda)`` per draw,
mirroring R's warm-start pattern. At the global FW optimum the
two starts are equivalent (strictly convex objective) — this is
a finite-iter parity fix, not a methodology change.

R-parity fixture + test seam:
* ``benchmarks/R/generate_sdid_placebo_parity_fixture.R`` — R 4.5.2
  + synthdid 0.0.9. Reuses the same Y matrix as
  ``TestJackknifeSERParity`` (same R_ATT = 4.980848860060929) so
  jackknife and placebo R-parity tests share an anchor panel.
  Records the 200 per-rep permutations R consumed and the SE from
  R's manual ``placebo_se`` loop (which matches
  ``vcov(method="placebo")`` to machine precision when seeded);
  permutations are 0-indexed for direct numpy consumption.
* ``tests/data/sdid_placebo_indices_r.json`` — committed fixture.
* ``_placebo_variance_se`` gains a private ``_placebo_indices``
  kwarg (underscore-prefixed, test-only). When supplied, each row
  replaces the per-draw ``rng.permutation(n_control)`` so a
  Python fit can consume R's exact permutation sequence and
  produce a bit-identical SE.
* ``test_placebo_se_matches_r`` (in ``TestJackknifeSERParity``)
  intercepts the dispatcher's call to ``_placebo_variance_se`` to
  capture the normalized fit-time inputs, then re-invokes the
  method with R's permutations through the seam. Asserts
  ``|py_se - r_se| < 1e-8`` — Rust FW vs R FW differ at sub-ULP
  on the same warm-start; tight enough to catch real divergences
  without masking BLAS reduction-order tolerance.

Baseline rebase:
``TestScaleEquivariance::test_baseline_parity_small_scale[placebo]``
captured pre-warm-start SE = 0.29385822261006445. New value is
0.293840360160448 (sub-percent shift). The test's bit-identity
contract is preserved per backend; baseline updated with a comment
documenting the warm-start change and pointer to the new R-parity
test that pins the post-fix value to R's reference. p-value (placebo
uses empirical formula, not analytical) is unchanged at
0.004975124378109453.

Verification: ``pytest tests/test_methodology_sdid.py
tests/test_survey_phase5.py -q`` → 230 passed (1 new R-parity
test; existing TestScaleEquivariance baseline rebased; all other
SDID + survey tests unchanged).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: igerber

CHANGELOG.md

@@ -7,6 +7,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [3.3.0] - 2026-04-25
 
+### Fixed
+- **`SyntheticDiD(variance_method="placebo")` SE now uses R-default warm-start** matching `synthdid:::placebo_se`. R's placebo loop seeds Frank-Wolfe per draw with `weights.boot$omega = sum_normalize(weights$omega[ind[1:N0_placebo]])` (fit-time ω subsetted + renormalized) and the fit-time `weights$lambda` — Python previously used uniform cold-start, producing finite-iter convergence-pattern drift on a handful of draws relative to R's reference SE. New `_placebo_variance_se` kwargs `init_omega` / `init_lambda` thread fit-time weights through the existing two-pass FW dispatcher; on the global FW optimum the values are init-independent (strictly convex objective), so the change is a finite-iter parity fix, not a methodology change. Existing placebo SE values shift by sub-percent on most panels; the bit-identity baseline pin in `TestScaleEquivariance::test_baseline_parity_small_scale[placebo]` was rebased from `0.29385822261006445` to `0.293840360160448`. New R-parity test `tests/test_methodology_sdid.py::TestJackknifeSERParity::test_placebo_se_matches_r` asserts SE matches R's `vcov(method="placebo")` to within `< 1e-8` using R's exact permutation sequence (recorded by `benchmarks/R/generate_sdid_placebo_parity_fixture.R` into `tests/data/sdid_placebo_indices_r.json`). The `_placebo_indices` kwarg on `_placebo_variance_se` is the test seam; not part of the public API.
+
 ### Added
 - **`qug_test` and `did_had_pretest_workflow` survey-aware NotImplementedError gates (Phase 4.5 C0 decision gate).** `qug_test(d, *, survey=None, weights=None)` and `did_had_pretest_workflow(..., *, survey=None, weights=None)` now accept the two kwargs as keyword-only with default `None`. Passing either non-`None` raises `NotImplementedError` with an educational message naming the methodology rationale and pointing users to joint Stute (Phase 4.5 C, planned) as the survey-compatible alternative. Mutex guard on `survey=` + `weights=` mirrors `HeterogeneousAdoptionDiD.fit()` at `had.py:2890`. **QUG-under-survey is permanently deferred** — the test statistic uses extreme order statistics `D_{(1)}, D_{(2)}` which are NOT smooth functionals of the empirical CDF, so standard survey machinery (Binder-TSL linearization, Rao-Wu rescaled bootstrap, Krieger-Pfeffermann (1997) EDF tests) does not yield a calibrated test; under cluster sampling the `Exp(1)/Exp(1)` limit law's independence assumption breaks; and the EVT-under-unequal-probability-sampling literature (Quintos et al. 2001, Beirlant et al.) addresses tail-index estimation, not boundary tests. The workflow's gate is **temporary** — Phase 4.5 C will close it for the linearity-family pretests with mechanism varying by test: Rao-Wu rescaled bootstrap for `stute_test` and the joint variants (`stute_joint_pretest`, `joint_pretrends_test`, `joint_homogeneity_test`); weighted OLS residuals + weighted variance estimator for `yatchew_hr_test` (Yatchew 1997 is a closed-form variance-ratio test, not bootstrap-based). Sister pretests (`stute_test`, `yatchew_hr_test`, `stute_joint_pretest`, `joint_pretrends_test`, `joint_homogeneity_test`) keep their closed signatures in this release — Phase 4.5 C will add kwargs and implementation together to avoid API churn. Unweighted `qug_test(d)` and `did_had_pretest_workflow(...)` calls are bit-exact pre-PR (kwargs are keyword-only after `*`; positional path unchanged). New tests at `tests/test_had_pretests.py::TestQUGTest` (5 rejection / mutex / message / regression tests) and the new `TestHADPretestWorkflowSurveyGuards` class (6 tests covering both kwarg paths, mutex, methodology pointer, both aggregate paths, and unweighted regression). See `docs/methodology/REGISTRY.md` § "QUG Null Test" — Note (Phase 4.5 C0) for the full methodology rationale plus a sketch of the (out-of-scope) theoretical bridge that combines endpoint-estimation EVT (Hall 1982, Aarssen-de Haan 1994, Hall-Wang 1999, Beirlant-de Wet-Goegebeur 2006), survey-aware functional CLTs (Boistard-Lopuhaä-Ruiz-Gazen 2017, Bertail-Chautru-Clémençon 2017), and tail-empirical-process theory (Drees 2003) — publishable methodology research, not engineering work.
 - **`HeterogeneousAdoptionDiD` mass-point `survey=` / `weights=` + event-study `aggregate="event_study"` survey composition + multiplier-bootstrap sup-t simultaneous confidence band (Phase 4.5 B).** Closes the two Phase 4.5 A `NotImplementedError` gates: `design="mass_point" + weights/survey` and `aggregate="event_study" + weights/survey`. Weighted 2SLS sandwich in `_fit_mass_point_2sls` follows the Wooldridge 2010 Ch. 12 pweight convention (`w²` in the HC1 meat, `w·u` in the CR1 cluster score, weighted bread `Z'WX`); HC1 and CR1 ("stata" `se_type`) bit-parity with `estimatr::iv_robust(..., weights=, clusters=)` at `atol=1e-10` (new cross-language golden at `benchmarks/data/estimatr_iv_robust_golden.json`, generated by `benchmarks/R/generate_estimatr_iv_robust_golden.R`; `estimatr` added to `benchmarks/R/requirements.R`). `_fit_mass_point_2sls` gains `weights=` + `return_influence=` kwargs and now always returns a 3-tuple `(beta, se, psi)` — `psi` is the per-unit IF on the β̂-scale scaled so `compute_survey_if_variance(psi, trivial_resolved) ≈ V_HC1[1,1]` at `atol=1e-10` (PR #359 IF scale convention applied uniformly; no `sum(psi²)` claims). Event-study per-horizon variance: `survey=` path composes Binder-TSL via `compute_survey_if_variance`; `weights=` shortcut uses the analytical weighted-robust SE (continuous: CCT-2014 `bc_fit.se_robust / |den|`; mass-point: weighted 2SLS pweight sandwich from `_fit_mass_point_2sls` — HC1 / classical / CR1). `survey_metadata` / `variance_formula` / `effective_dose_mean` populated in both regimes (previously hardcoded `None` at `had.py:3366`). New multiplier-bootstrap sup-t: `_sup_t_multiplier_bootstrap` reuses `diff_diff.bootstrap_utils.generate_survey_multiplier_weights_batch` for PSU-level draws with stratum centering + sqrt(n_h/(n_h-1)) small-sample correction + FPC scaling + lonely-PSU handling. On the `weights=` shortcut, sup-t calibration is routed through a synthetic trivial `ResolvedSurveyDesign` so the centered + small-sample-corrected branch fires uniformly — targets the analytical HC1 variance family (`compute_survey_if_variance(IF, trivial) ≈ V_HC1` per the PR #359 IF scale invariant) rather than the raw `sum(ψ²) = ((n-1)/n) · V_HC1` that unit-level Rademacher multipliers would produce on the HC1-scaled IF. Perturbations: `delta = weights @ IF` with NO `(1/n)` prefactor (matching `staggered_bootstrap.py:373` idiom), normalized by per-horizon analytical SE, `(1-alpha)`-quantile of the sup-t distribution. At H=1 the quantile reduces to `Φ⁻¹(1 − alpha/2) ≈ 1.96` up to MC noise (regression-locked by `TestSupTReducesToNormalAtH1`). `HeterogeneousAdoptionDiD.__init__` gains `n_bootstrap: int = 999` and `seed: Optional[int] = None` (CS-parity singular seed); `fit()` gains `cband: bool = True` (only consulted on weighted event-study). `HeterogeneousAdoptionDiDEventStudyResults` extended with `variance_formula`, `effective_dose_mean`, `cband_low`, `cband_high`, `cband_crit_value`, `cband_method`, `cband_n_bootstrap` (all `None` on unweighted fits); surfaced in `to_dict`, `to_dataframe`, `summary`, `__repr__`. Unweighted event-study with `cband=False` preserves pre-Phase 4.5 B numerical output bit-exactly (stability invariant, locked by regression tests). Zero-weight subpopulation convention carries over from PR #359 (filter for design decisions; preserve full ResolvedSurveyDesign for variance). Non-pweight SurveyDesigns (`aweight`, `fweight`, replicate designs) raise `NotImplementedError` on both new paths (reciprocal-guard discipline). Pretest surfaces (`qug_test`, `stute_test`, `yatchew_hr_test`, joint variants, `did_had_pretest_workflow`) remain unweighted in this release — Phase 4.5 C / C0. See `docs/methodology/REGISTRY.md` §HeterogeneousAdoptionDiD "Weighted 2SLS (Phase 4.5 B)", "Event-study survey composition", and "Sup-t multiplier bootstrap" for derivations and invariants.

benchmarks/R/generate_sdid_placebo_parity_fixture.R

@@ -0,0 +1,142 @@
+#!/usr/bin/env Rscript
+# Generate a fixture pinning R's `synthdid::vcov(method="placebo")` SE plus
+# the per-replication permutations R consumed, so the Python R-parity test
+# can feed those exact permutations through `_placebo_variance_se` and
+# assert SE match at machine precision.
+#
+# Usage:
+#   Rscript benchmarks/R/generate_sdid_placebo_parity_fixture.R
+#
+# Output:
+#   tests/data/sdid_placebo_indices_r.json
+#
+# Symmetric with the existing jackknife R-parity test
+# (TestJackknifeSERParity in tests/test_methodology_sdid.py:1410). Reuses
+# the same Y matrix and (N0, N1, T0, T1) shape so the placebo + jackknife
+# parity tests share an anchor panel.
+#
+# R version: 4.5.2; synthdid version: 0.0.9.
+
+library(synthdid)
+library(jsonlite)
+
+# Reconstruct R's panel exactly as TestJackknifeSERParity does (set.seed(42),
+# 23 units × 8 periods, treated = i > N0 with effect 5 in t > T0).
+set.seed(42)
+N0 <- 20
+N1 <- 3
+T0 <- 5
+T1 <- 3
+N <- N0 + N1
+T <- T0 + T1
+Y <- matrix(0, nrow = N, ncol = T)
+for (i in 1:N) {
+  unit_fe <- rnorm(1, sd = 2)
+  for (t in 1:T) {
+    Y[i, t] <- 10 + unit_fe + (t - 1) * 0.3 + rnorm(1, sd = 0.5)
+    if (i > N0 && t > T0) Y[i, t] <- Y[i, t] + 5.0
+  }
+}
+
+# Fit-time ATT (sanity check — must match TestJackknifeSERParity.R_ATT).
+tau_hat <- synthdid_estimate(Y, N0, T0)
+r_att <- as.numeric(tau_hat)
+
+# Reproduce R's placebo_se loop exactly so we can record permutations and
+# the per-rep tau alongside the resulting SE. Mirrors `synthdid:::placebo_se`
+# (R/vcov.R), including the warm-start weights pass-through:
+#
+#   theta = function(ind) {
+#     N0 = length(ind) - N1
+#     weights.boot = weights
+#     weights.boot$omega = sum_normalize(weights$omega[ind[1:N0]])
+#     do.call(synthdid_estimate, c(list(Y = setup$Y[ind, ],
+#         N0 = N0, T0 = setup$T0, X = setup$X[ind, , ],
+#         weights = weights.boot), opts))
+#   }
+#
+# The warm-start `weights.boot$omega` differs from a fresh uniform init
+# at finite FW iterations and is what `vcov(method="placebo")` actually
+# consumes — so reproducing it here is required for bit-identical SE.
+opts_used <- attr(tau_hat, "opts")
+fit_weights <- attr(tau_hat, "weights")
+fit_setup <- attr(tau_hat, "setup")
+replications <- 200
+
+# Use a fresh seed for the placebo loop so the recorded permutations are
+# independent of the fit-time RNG state. Python consumes the recorded
+# permutations directly (no RNG-state matching needed).
+set.seed(42)
+perms <- vector("list", replications)
+taus <- numeric(replications)
+
+for (r in 1:replications) {
+  ind <- sample(1:N0, N0)
+  perms[[r]] <- ind
+  N0_placebo <- N0 - N1
+  weights_boot <- fit_weights
+  weights_boot$omega <- synthdid:::sum_normalize(fit_weights$omega[ind[1:N0_placebo]])
+  # IMPORTANT: R's `placebo_se` uses ONLY the N0 controls (subdivided into
+  # N0-N1 pseudo-controls + N1 pseudo-treated). Real treated rows are NOT
+  # included in the placebo Y matrix — that's what makes the placebo a
+  # null-distribution test. ``Y = setup$Y[ind, ]`` is N0 rows; appending
+  # the real treated rows (i.e., ``setup$Y[c(ind, (N0+1):N), ]``) would
+  # change the test entirely (and produces SE ~0.132 instead of R's 0.226
+  # — a 2× drift on this fixture).
+  est_placebo <- do.call(
+    synthdid_estimate,
+    c(list(
+      Y = fit_setup$Y[ind, ],
+      N0 = N0_placebo,
+      T0 = T0,
+      X = fit_setup$X[ind, , ],
+      weights = weights_boot
+    ), opts_used)
+  )
+  taus[r] <- as.numeric(est_placebo)
+}
+
+r_placebo_se <- sqrt((replications - 1) / replications) * sd(taus)
+
+# Sanity check against R's vcov() entry point. With the warm-start pattern
+# applied explicitly above, the manual loop and `vcov()` should produce
+# the same SE up to MC noise on the seed sequence. Match isn't required
+# for the parity test (we use `r_placebo_se` from our recorded
+# permutations); both values are kept for transparency.
+set.seed(42)
+r_placebo_se_via_vcov <- sqrt(vcov(tau_hat, method = "placebo", replications = replications)[1, 1])
+
+cat(sprintf("R ATT:                       %.15f\n", r_att))
+cat(sprintf("R placebo SE (manual loop):  %.15f\n", r_placebo_se))
+cat(sprintf("R placebo SE (via vcov):     %.15f\n", r_placebo_se_via_vcov))
+cat(sprintf("Replications:                %d\n", replications))
+
+# Convert permutations to 0-indexed for Python (R uses 1-indexed).
+perms_0indexed <- lapply(perms, function(p) as.integer(p - 1L))
+
+payload <- list(
+  metadata = list(
+    R_version = paste(R.version$major, R.version$minor, sep = "."),
+    synthdid_version = as.character(packageVersion("synthdid")),
+    seed = 42L,
+    replications = as.integer(replications),
+    note = paste(
+      "Permutations are 0-indexed for direct numpy consumption.",
+      "R ATT, R placebo SE (manual loop), and per-rep taus are pinned",
+      "for downstream Python parity assertion."
+    )
+  ),
+  N0 = as.integer(N0),
+  N1 = as.integer(N1),
+  T0 = as.integer(T0),
+  T1 = as.integer(T1),
+  R_ATT = r_att,
+  R_PLACEBO_SE = r_placebo_se,
+  R_PLACEBO_SE_VIA_VCOV = r_placebo_se_via_vcov,
+  R_PLACEBO_TAUS = as.numeric(taus),
+  R_PERMUTATIONS = perms_0indexed
+)
+
+out_path <- "tests/data/sdid_placebo_indices_r.json"
+write_json(payload, out_path, auto_unbox = TRUE, digits = 17)
+cat(sprintf("\nWrote %s\n", out_path))

diff_diff/synthetic_did.py

@@ -1153,6 +1153,8 @@ def fit(  # type: ignore[override]
                     min_decrease=min_decrease,
                     replications=self.n_bootstrap,
                     w_control=w_control,
+                    init_omega=unit_weights,
+                    init_lambda=time_weights,
                 )
             se = se_n * Y_scale
             placebo_effects = np.asarray(placebo_effects_n) * Y_scale
@@ -1695,6 +1697,9 @@ def _placebo_variance_se(
         min_decrease: float = 1e-5,
         replications: int = 200,
         w_control=None,
+        init_omega: Optional[np.ndarray] = None,
+        init_lambda: Optional[np.ndarray] = None,
+        _placebo_indices: Optional[np.ndarray] = None,
     ) -> Tuple[float, np.ndarray]:
         """
         Compute placebo-based variance matching R's synthdid methodology.
@@ -1747,6 +1752,21 @@ def _placebo_variance_se(
         rng = np.random.default_rng(self.seed)
         n_pre, n_control = Y_pre_control.shape
 
+        # R-parity test seam (PR follow-up to #349). When
+        # ``_placebo_indices`` is provided, each row replaces the
+        # per-draw ``rng.permutation(n_control)`` so a Python fit can
+        # consume R's exact permutation sequence and produce a
+        # bit-identical SE. Shape: ``(replications, n_control)``,
+        # 0-indexed. Underscore-prefixed → test-only / private API.
+        if _placebo_indices is not None:
+            _placebo_indices = np.asarray(_placebo_indices, dtype=np.int64)
+            if _placebo_indices.ndim != 2 or _placebo_indices.shape[1] != n_control:
+                raise ValueError(
+                    f"_placebo_indices shape {_placebo_indices.shape} does not "
+                    f"match expected (replications, {n_control})"
+                )
+            replications = _placebo_indices.shape[0]
+
         # Ensure we have enough controls for the split
         n_pseudo_control = n_control - n_treated
         if n_pseudo_control < 1:
@@ -1771,10 +1791,17 @@ def _placebo_variance_se(
 
         placebo_estimates = []
 
-        for _ in range(replications):
+        for _rep in range(replications):
             try:
-                # Random permutation of control indices (Algorithm 4, step 1)
-                perm = rng.permutation(n_control)
+                # Random permutation of control indices (Algorithm 4, step 1).
+                # Test seam: when ``_placebo_indices`` is supplied, consume
+                # the externally-provided permutation instead — used by
+                # ``test_placebo_se_matches_r`` to feed R's exact
+                # permutation sequence through Python for bit-identical SE.
+                if _placebo_indices is not None:
+                    perm = _placebo_indices[_rep]
+                else:
+                    perm = rng.permutation(n_control)
 
                 # Split into pseudo-controls and pseudo-treated (step 2)
                 pseudo_control_idx = perm[:n_pseudo_control]
@@ -1805,15 +1832,28 @@ def _placebo_variance_se(
                         Y_post_control[:, pseudo_treated_idx], axis=1
                     )
 
-                # Re-estimate weights on permuted data (matching R's behavior)
-                # R passes update.omega=TRUE, update.lambda=TRUE via opts,
-                # re-estimating weights from uniform initialization (fresh start).
-                # Unit weights: re-estimate on pseudo-control/pseudo-treated data
+                # Re-estimate weights on permuted data (matching R's behavior).
+                # R's `placebo_se` (synthdid:::placebo_se in R/vcov.R) passes
+                # `weights.boot$omega = sum_normalize(weights$omega[ind[1:N0]])`
+                # as a warm-start to `synthdid_estimate` with `update.omega=TRUE`,
+                # then re-estimates ω. Python mirrors that: when fit-time ω
+                # (`init_omega`) is supplied, seed the FW first pass with the
+                # subsetted-renormalized fit-time omega; otherwise use cold-
+                # start (uniform). At the global FW optimum the two are
+                # equivalent, but the warm-start matches R's exact iterates
+                # for bit-identical SE under the R-parity test.
+                if init_omega is not None:
+                    pseudo_omega_init = _sum_normalize(
+                        init_omega[pseudo_control_idx]
+                    )
+                else:
+                    pseudo_omega_init = None
                 pseudo_omega = compute_sdid_unit_weights(
                     Y_pre_pseudo_control,
                     Y_pre_pseudo_treated_mean,
                     zeta_omega=zeta_omega,
                     min_decrease=min_decrease,
+                    init_weights=pseudo_omega_init,
                 )
 
                 # Compose pseudo_omega with control survey weights
@@ -1824,12 +1864,17 @@ def _placebo_variance_se(
                 else:
                     pseudo_omega_eff = pseudo_omega
 
-                # Time weights: re-estimate on pseudo-control data
+                # Time weights: re-estimate on pseudo-control data.
+                # R's `placebo_se` reuses the same `weights.boot` (with
+                # the original ``weights$lambda`` untouched) as warm-
+                # start to ``synthdid_estimate``. Mirror by passing
+                # fit-time λ as ``init_weights`` when supplied.
                 pseudo_lambda = compute_time_weights(
                     Y_pre_pseudo_control,
                     Y_post_pseudo_control,
                     zeta_lambda=zeta_lambda,
                     min_decrease=min_decrease,
+                    init_weights=init_lambda,
                 )
 
                 # Compute placebo SDID estimate (step 4)