SDID methodology review: rewrite to match R synthdid + Rust parallel variance

CLAUDE.md

@@ -68,6 +68,10 @@ cross-platform compilation - no OpenBLAS or Intel MKL installation required.
 
 - **`diff_diff/synthetic_did.py`** - Synthetic DiD estimator:
   - `SyntheticDiD` - Synthetic control combined with DiD (Arkhangelsky et al. 2021)
+  - Frank-Wolfe solver matching R's `synthdid::sc.weight.fw()` for unit and time weights
+  - Auto-computed regularization from data noise level (`zeta_omega`, `zeta_lambda`)
+  - Two-pass sparsification for unit weights (100 iters → sparsify → 1000 iters)
+  - Bootstrap SE uses fixed weights (matches R's `bootstrap_sample`)
 
 - **`diff_diff/staggered.py`** - Staggered adoption DiD main module:
   - `CallawaySantAnna` - Callaway & Sant'Anna (2021) estimator for heterogeneous treatment timing
@@ -163,6 +167,9 @@ cross-platform compilation - no OpenBLAS or Intel MKL installation required.
     - `compute_unit_distance_matrix()` - Parallel pairwise RMSE distance computation (4-8x speedup)
     - `loocv_grid_search()` - Parallel LOOCV across tuning parameters (10-50x speedup)
     - `bootstrap_trop_variance()` - Parallel bootstrap variance estimation (5-15x speedup)
+  - **`rust/src/sdid_variance.rs`** - SDID variance estimation acceleration:
+    - `placebo_variance_sdid()` - Parallel placebo SE computation (~8x speedup)
+    - `bootstrap_variance_sdid()` - Parallel bootstrap SE computation (~6x speedup)
   - Uses pure-Rust `faer` library for linear algebra (no external BLAS/LAPACK dependencies)
   - Cross-platform: builds on Linux, macOS, and Windows without additional setup
   - Provides 4-8x speedup for SyntheticDiD, 5-20x speedup for TROP
@@ -373,6 +380,7 @@ See `docs/benchmarks.rst` for full methodology and validation results.
 
 Tests mirror the source modules:
 - `tests/test_estimators.py` - Tests for DifferenceInDifferences, TWFE, MultiPeriodDiD, SyntheticDiD
+- `tests/test_methodology_sdid.py` - Methodology tests for SDID: Frank-Wolfe solver, regularization, sparsification, edge cases
 - `tests/test_staggered.py` - Tests for CallawaySantAnna
 - `tests/test_sun_abraham.py` - Tests for SunAbraham interaction-weighted estimator
 - `tests/test_imputation.py` - Tests for ImputationDiD (Borusyak et al. 2024) estimator

METHODOLOGY_REVIEW.md

@@ -25,7 +25,7 @@ Each estimator in diff-diff should be periodically reviewed to ensure:
 | TwoWayFixedEffects | `twfe.py` | `fixest::feols()` | **Complete** | 2026-02-08 |
 | CallawaySantAnna | `staggered.py` | `did::att_gt()` | **Complete** | 2026-01-24 |
 | SunAbraham | `sun_abraham.py` | `fixest::sunab()` | Not Started | - |
-| SyntheticDiD | `synthetic_did.py` | `synthdid::synthdid_estimate()` | Not Started | - |
+| SyntheticDiD | `synthetic_did.py` | `synthdid::synthdid_estimate()` | **Complete** | 2026-02-10 |
 | TripleDifference | `triple_diff.py` | (forthcoming) | Not Started | - |
 | TROP | `trop.py` | (forthcoming) | Not Started | - |
 | BaconDecomposition | `bacon.py` | `bacondecomp::bacon()` | Not Started | - |
@@ -314,14 +314,37 @@ variables appear to the left of the `|` separator.
 | Module | `synthetic_did.py` |
 | Primary Reference | Arkhangelsky et al. (2021) |
 | R Reference | `synthdid::synthdid_estimate()` |
-| Status | Not Started |
-| Last Review | - |
+| Status | **Complete** |
+| Last Review | 2026-02-10 |
 
 **Corrections Made:**
-- (None yet)
+1. **Time weights: Frank-Wolfe on collapsed form** (was heuristic inverse-distance).
+   Replaced ad-hoc inverse-distance weighting with the Frank-Wolfe algorithm operating
+   on the collapsed (N_co x T_pre) problem as specified in Algorithm 1 of
+   Arkhangelsky et al. (2021), matching R's `synthdid::fw.step()`.
+2. **Unit weights: Frank-Wolfe with two-pass sparsification** (was projected gradient
+   descent with wrong penalty). Replaced projected gradient descent (which used an
+   incorrect penalty formulation) with Frank-Wolfe optimization followed by two-pass
+   sparsification, matching R's `synthdid::sc.weight.fw()` and `sparsify_function()`.
+3. **Auto-computed regularization from data noise level** (was `lambda_reg=0.0`,
+   `zeta=1.0`). Regularization parameters `zeta_omega` and `zeta_lambda` are now
+   computed automatically from the data noise level (N_tr * sigma^2) as specified in
+   Appendix D of Arkhangelsky et al. (2021), matching R's default behavior.
+4. **Bootstrap SE uses fixed weights matching R's `bootstrap_sample`** (was
+   re-estimating all weights). The bootstrap variance procedure now holds unit and time
+   weights fixed at their point estimates and only re-estimates the treatment effect,
+   matching the approach in R's `synthdid::bootstrap_sample()`.
+5. **Default `variance_method` changed to `"placebo"`** matching R's default. The R
+   package uses placebo variance by default (`synthdid_estimate` returns an object whose
+   `vcov()` uses the placebo method); our default now matches.
+6. **Deprecated `lambda_reg` and `zeta` params; new params are `zeta_omega` and
+   `zeta_lambda`**. The old parameters had unclear semantics and did not correspond to
+   the paper's notation. The new parameters directly match the paper and R package
+   naming conventions. `lambda_reg` and `zeta` are deprecated with warnings and will
+   be removed in a future release.
 
 **Outstanding Concerns:**
-- (None yet)
+- (None)
 
 ---
 

README.md

@@ -1163,11 +1163,12 @@ Use Synthetic DiD instead of standard DiD when:
 
 ```python
 SyntheticDiD(
-    lambda_reg=0.0,     # Regularization toward uniform weights (0 = no reg)
-    zeta=1.0,           # Time weight regularization (higher = more uniform)
-    alpha=0.05,         # Significance level
-    n_bootstrap=200,    # Bootstrap iterations for SE (0 = placebo-based)
-    seed=None           # Random seed for reproducibility
+    zeta_omega=None,        # Unit weight regularization (None = auto-computed from data)
+    zeta_lambda=None,       # Time weight regularization (None = auto-computed from data)
+    alpha=0.05,             # Significance level
+    variance_method="placebo",  # "placebo" (default, matches R) or "bootstrap"
+    n_bootstrap=200,        # Replications for SE estimation
+    seed=None               # Random seed for reproducibility
 )
 ```
 
@@ -1872,11 +1873,12 @@ MultiPeriodDiD(
 
 ```python
 SyntheticDiD(
-    lambda_reg=0.0,     # L2 regularization for unit weights
-    zeta=1.0,           # Regularization for time weights
-    alpha=0.05,         # Significance level for CIs
-    n_bootstrap=200,    # Bootstrap iterations for SE
-    seed=None           # Random seed for reproducibility
+    zeta_omega=None,        # Unit weight regularization (None = auto from data)
+    zeta_lambda=None,       # Time weight regularization (None = auto from data)
+    alpha=0.05,             # Significance level for CIs
+    variance_method="placebo",  # "placebo" (R default) or "bootstrap"
+    n_bootstrap=200,        # Replications for SE estimation
+    seed=None               # Random seed for reproducibility
 )
 ```
 

benchmarks/R/benchmark_synthdid.R

@@ -82,6 +82,15 @@ se_time <- as.numeric(difftime(Sys.time(), se_start, units = "secs"))
 
 total_time <- estimation_time + se_time
 
+# Compute noise level and regularization (to match Python's auto-computed values)
+N0 <- setup$N0
+T0 <- setup$T0
+N1 <- nrow(setup$Y) - N0
+T1 <- ncol(setup$Y) - T0
+noise_level <- sd(apply(setup$Y[1:N0, 1:T0], 1, diff))
+zeta_omega <- ((N1 * T1)^(1/4)) * noise_level
+zeta_lambda <- 1e-6 * noise_level
+
 # Format output
 results <- list(
   estimator = "synthdid::synthdid_estimate",
@@ -90,10 +99,15 @@ results <- list(
   att = as.numeric(tau_hat),
   se = se,
 
-  # Weights
+  # Weights (full precision)
   unit_weights = as.numeric(unit_weights),
   time_weights = as.numeric(time_weights),
 
+  # Regularization parameters
+  noise_level = noise_level,
+  zeta_omega = zeta_omega,
+  zeta_lambda = zeta_lambda,
+
   # Timing
   timing = list(
     estimation_seconds = estimation_time,
@@ -105,10 +119,10 @@ results <- list(
   metadata = list(
     r_version = R.version.string,
     synthdid_version = as.character(packageVersion("synthdid")),
-    n_control = setup$N0,
-    n_treated = nrow(setup$Y) - setup$N0,
-    n_pre_periods = setup$T0,
-    n_post_periods = ncol(setup$Y) - setup$T0
+    n_control = N0,
+    n_treated = N1,
+    n_pre_periods = T0,
+    n_post_periods = T1
   )
 )
 

benchmarks/python/benchmark_synthdid.py

@@ -106,9 +106,13 @@ def main():
         # Point estimate and SE
         "att": float(results.att),
         "se": float(results.se),
-        # Weights
+        # Weights (full precision)
         "unit_weights": unit_weights_df["weight"].tolist(),
         "time_weights": time_weights_df["weight"].tolist(),
+        # Regularization parameters
+        "noise_level": float(results.noise_level) if results.noise_level is not None else None,
+        "zeta_omega": float(results.zeta_omega) if results.zeta_omega is not None else None,
+        "zeta_lambda": float(results.zeta_lambda) if results.zeta_lambda is not None else None,
         # Timing
         "timing": {
             "estimation_seconds": total_time,

diff_diff/_backend.py

@@ -30,6 +30,14 @@
         # TROP estimator acceleration (joint method)
         loocv_grid_search_joint as _rust_loocv_grid_search_joint,
         bootstrap_trop_variance_joint as _rust_bootstrap_trop_variance_joint,
+        # SDID weights (Frank-Wolfe matching R's synthdid)
+        compute_sdid_unit_weights as _rust_sdid_unit_weights,
+        compute_time_weights as _rust_compute_time_weights,
+        compute_noise_level as _rust_compute_noise_level,
+        sc_weight_fw as _rust_sc_weight_fw,
+        # SDID variance estimation (parallel placebo and bootstrap)
+        placebo_variance_sdid as _rust_placebo_variance_sdid,
+        bootstrap_variance_sdid as _rust_bootstrap_variance_sdid,
     )
     _rust_available = True
 except ImportError:
@@ -46,6 +54,14 @@
     # TROP estimator acceleration (joint method)
     _rust_loocv_grid_search_joint = None
     _rust_bootstrap_trop_variance_joint = None
+    # SDID weights (Frank-Wolfe matching R's synthdid)
+    _rust_sdid_unit_weights = None
+    _rust_compute_time_weights = None
+    _rust_compute_noise_level = None
+    _rust_sc_weight_fw = None
+    # SDID variance estimation (parallel placebo and bootstrap)
+    _rust_placebo_variance_sdid = None
+    _rust_bootstrap_variance_sdid = None
 
 # Determine final backend based on environment variable and availability
 if _backend_env == 'python':
@@ -63,6 +79,14 @@
     # TROP estimator acceleration (joint method)
     _rust_loocv_grid_search_joint = None
     _rust_bootstrap_trop_variance_joint = None
+    # SDID weights (Frank-Wolfe matching R's synthdid)
+    _rust_sdid_unit_weights = None
+    _rust_compute_time_weights = None
+    _rust_compute_noise_level = None
+    _rust_sc_weight_fw = None
+    # SDID variance estimation (parallel placebo and bootstrap)
+    _rust_placebo_variance_sdid = None
+    _rust_bootstrap_variance_sdid = None
 elif _backend_env == 'rust':
     # Force Rust mode - fail if not available
     if not _rust_available:
@@ -89,4 +113,12 @@
     # TROP estimator acceleration (joint method)
     '_rust_loocv_grid_search_joint',
     '_rust_bootstrap_trop_variance_joint',
+    # SDID weights (Frank-Wolfe matching R's synthdid)
+    '_rust_sdid_unit_weights',
+    '_rust_compute_time_weights',
+    '_rust_compute_noise_level',
+    '_rust_sc_weight_fw',
+    # SDID variance estimation (parallel placebo and bootstrap)
+    '_rust_placebo_variance_sdid',
+    '_rust_bootstrap_variance_sdid',
 ]

diff_diff/results.py

@@ -605,8 +605,10 @@ class SyntheticDiDResults:
     pre_periods: List[Any]
     post_periods: List[Any]
     alpha: float = 0.05
-    variance_method: str = field(default="bootstrap")
-    lambda_reg: Optional[float] = field(default=None)
+    variance_method: str = field(default="placebo")
+    noise_level: Optional[float] = field(default=None)
+    zeta_omega: Optional[float] = field(default=None)
+    zeta_lambda: Optional[float] = field(default=None)
     pre_treatment_fit: Optional[float] = field(default=None)
     placebo_effects: Optional[np.ndarray] = field(default=None)
     n_bootstrap: Optional[int] = field(default=None)
@@ -650,8 +652,12 @@ def summary(self, alpha: Optional[float] = None) -> str:
             f"{'Post-treatment periods:':<25} {len(self.post_periods):>10}",
         ]
 
-        if self.lambda_reg is not None:
-            lines.append(f"{'Regularization (lambda):':<25} {self.lambda_reg:>10.4f}")
+        if self.zeta_omega is not None:
+            lines.append(f"{'Zeta (unit weights):':<25} {self.zeta_omega:>10.4f}")
+        if self.zeta_lambda is not None:
+            lines.append(f"{'Zeta (time weights):':<25} {self.zeta_lambda:>10.6f}")
+        if self.noise_level is not None:
+            lines.append(f"{'Noise level:':<25} {self.noise_level:>10.4f}")
 
         if self.pre_treatment_fit is not None:
             lines.append(f"{'Pre-treatment fit (RMSE):':<25} {self.pre_treatment_fit:>10.4f}")
@@ -731,7 +737,9 @@ def to_dict(self) -> Dict[str, Any]:
             "n_pre_periods": len(self.pre_periods),
             "n_post_periods": len(self.post_periods),
             "variance_method": self.variance_method,
-            "lambda_reg": self.lambda_reg,
+            "noise_level": self.noise_level,
+            "zeta_omega": self.zeta_omega,
+            "zeta_lambda": self.zeta_lambda,
             "pre_treatment_fit": self.pre_treatment_fit,
         }
         if self.n_bootstrap is not None: