test: Update reference parity tests for all estimators

Author: brycewang-stanford

tests/r_parity/03_hdfe.py

@@ -44,7 +44,7 @@ def main() -> None:
     df = make_panel()
     dump_csv(df, MODULE)
 
-    fit = sp.fast.feols(FORMULA, data=df, vcov="iid")
+    fit = sp.fast.feols(FORMULA, data=df, vcov="iid", ssc="fixest")
 
     # FeolsResult exposes .coef() and .se() as methods returning a
     # pandas Series indexed by the non-absorbed regressor names.
@@ -74,17 +74,11 @@ def main() -> None:
             "n_years": int(df["year"].nunique()),
             "df_resid": int(fit.df_resid),
             "rss": float(fit.rss),
-            "df_convention": (
-                "sp.fast.feols uses df_resid = N - sum(FE_levels). "
-                "fixest::feols uses df_resid = N - sum(FE_levels) - 1, "
-                "subtracting one extra for the cross-FE collinearity. "
-                "The point estimates are identical to ~1e-15; the IID "
-                "SE ratio fixest/sp.fast is exactly "
-                "sqrt(sp_df / fixest_df) = sqrt(9730 / 9729) = "
-                "1.0000514, so the rel SE diff is 5.14e-5. Reviewers "
-                "should treat this as a small-sample-correction "
-                "convention difference, not as a numerical bug; "
-                "fixest's SE = sp.fast's SE * sqrt(9730/9729)."
+            "ssc": "fixest",
+            "parity_note": (
+                "sp.fast.feols uses ssc='fixest' so the residual FE "
+                "rank and IID standard errors match fixest::feols and "
+                "reghdfe on this two-way HDFE fixture."
             ),
         },
     )

tests/r_parity/05_sunab.py

@@ -1,8 +1,9 @@
 """StatsPAI Sun-Abraham event-study parity (Python side) -- Module 05.
 
-Runs sp.sun_abraham on the mpdta replica and emits the simple
-weighted average ATT and the dynamic event-study coefficients at
-relative times -3, -2, 0, 1, 2 (post-period anchor reference). The
+Runs sp.sun_abraham(..., aggregation="fixest_att") on the mpdta
+replica and emits the R/Stata-compatible weighted average ATT plus
+the dynamic event-study coefficients at relative times -3, -2, 0, 1,
+2 (post-period anchor reference). The
 companion R script runs fixest::feols(... | sunab(g, t)) on the
 same CSV.
 
@@ -23,7 +24,10 @@ def main() -> None:
     df = sp.datasets.mpdta()
     dump_csv(df, MODULE)
 
-    fit = sp.sun_abraham(df, y="lemp", g="first_treat", t="year", i="countyreal")
+    fit = sp.sun_abraham(
+        df, y="lemp", g="first_treat", t="year", i="countyreal",
+        aggregation="fixest_att",
+    )
 
     rows: list[ParityRecord] = [
         ParityRecord(
@@ -35,6 +39,14 @@ def main() -> None:
             n=int(len(df)),
         )
     ]
+    rows.append(
+        ParityRecord(
+            module=MODULE, side="py", statistic="event_time_avg_ATT",
+            estimate=float(fit.model_info["att_event_time"]),
+            se=float(fit.model_info["se_event_time"]),
+            n=int(len(df)),
+        )
+    )
 
     # Per-relative-time event-study coefficients.
     es = fit.model_info.get("event_study")
@@ -58,15 +70,14 @@ def main() -> None:
         extra={
             "control_group": "nevertreated",
             "method": fit.method,
-            "aggregation_note": (
-                "sp.sun_abraham reports the unweighted mean of the "
-                "post-treatment per-relative-time ATTs as its summary "
-                "estimate. fixest::feols with agg='att' uses a "
-                "cohort-share-weighted average that recovers the "
-                "Callaway-Sant'Anna simple ATT (-0.03298). On this "
-                "DGP the two summaries differ by ~2.4% (sp: -0.0338; "
-                "fixest: -0.0330). Per-relative-time event-study "
-                "coefficients agree at rel < 1e-11."
+            "aggregation": fit.model_info["summary_aggregation"],
+            "aggregation_parity_note": (
+                "The weighted_avg_ATT row uses "
+                "sp.sun_abraham(..., aggregation='fixest_att'), which "
+                "weights post-treatment cohort-time cells by treated "
+                "cohort size and matches fixest::summary(..., agg='att') "
+                "on the mpdta fixture. The historical equal-weighted "
+                "post-event-time summary is retained as event_time_avg_ATT."
             ),
         },
     )

tests/r_parity/06_rd.R

@@ -2,7 +2,7 @@
 #
 # Reads data/06_rd.csv (the StatsPAI Lee 2008 senate replica) and
 # runs rdrobust::rdrobust with package defaults (kernel = triangular,
-# p = 1, q = p + 1 = 2, bwselect = mserd). Tolerance: rel < 1e-3.
+# p = 1, q = p + 1 = 2, bwselect = mserd). Tolerance: rel < 1e-6.
 
 .args <- commandArgs(trailingOnly = FALSE)
 .file_arg <- grep("^--file=", .args, value = TRUE)

tests/r_parity/06_rd.py

@@ -1,12 +1,12 @@
 """StatsPAI RD CCT bias-corrected parity (Python side) -- Module 06.
 
-Runs sp.rdrobust on the Lee 2008 senate replica with the package
-defaults (kernel = triangular, p = 1, q = 2, bwselect = mserd) and
-emits both the conventional and the robust bias-corrected estimates.
-The companion R script runs rdrobust::rdrobust with identical
-defaults.
+Runs sp.rdrobust(..., bwselect="cct") on the Lee 2008 senate replica
+so the Track A row exercises the same Calonico-Cattaneo-Titiunik
+``rdrobust`` bandwidth selector used by R and Stata.  The legacy
+StatsPAI internal ``mserd`` selector is still recorded as a diagnostic
+row, but the parity headline is the canonical CCT path.
 
-Tolerance: rel < 1e-3 (iterative bandwidth selection).
+Tolerance: rel < 1e-6 against R/Stata CCT defaults.
 """
 from __future__ import annotations
 
@@ -22,8 +22,9 @@ def main() -> None:
     df = sp.datasets.lee_2008_senate()
     dump_csv(df, MODULE)
 
-    # Default mserd bandwidth selector.
-    fit = sp.rdrobust(df, y="voteshare_next", x="margin", c=0.0)
+    # Canonical R/Stata rdrobust bandwidth selector via the official
+    # rdrobust Python port.
+    fit = sp.rdrobust(df, y="voteshare_next", x="margin", c=0.0, bwselect="cct")
 
     rows: list[ParityRecord] = []
     for label in ("conventional", "robust"):
@@ -52,11 +53,31 @@ def main() -> None:
         )
     )
 
-    # Forced-bandwidth replicate at h = b = 0.042287 so the
-    # bandwidth-selector convention difference is isolated from the
-    # local-polynomial estimator math itself. R-side mirrors this
-    # with rdrobust(..., h = 0.042287, b = 0.042287).
-    H_FORCED = 0.042287
+    # Legacy internal-selector diagnostic.  This keeps the old default
+    # visible without mixing it with the R/Stata default-h parity rows.
+    legacy = sp.rdrobust(df, y="voteshare_next", x="margin", c=0.0)
+    rows.append(
+        ParityRecord(
+            module=MODULE, side="py",
+            statistic="legacy_internal_mserd_bandwidth_h",
+            estimate=float(legacy.model_info["bandwidth_h"]), n=int(len(df)),
+        )
+    )
+    rows.append(
+        ParityRecord(
+            module=MODULE, side="py",
+            statistic="legacy_internal_mserd_robust_est",
+            estimate=float(legacy.model_info["robust"]["estimate"]),
+            se=float(legacy.model_info["robust"]["se"]),
+            ci_lo=float(legacy.model_info["robust"]["ci"][0]),
+            ci_hi=float(legacy.model_info["robust"]["ci"][1]),
+            n=int(len(df)),
+        )
+    )
+
+    # Forced-bandwidth replicate at the legacy h = b = 0.042287 so the
+    # local-polynomial estimator math remains separately pinned.
+    H_FORCED = float(legacy.model_info["bandwidth_h"])
     fit_forced = sp.rdrobust(df, y="voteshare_next", x="margin", c=0.0,
                               h=H_FORCED, b=H_FORCED)
     for label in ("conventional", "robust"):
@@ -80,23 +101,14 @@ def main() -> None:
             "p": fit.model_info["polynomial_p"],
             "q": fit.model_info["polynomial_q"],
             "bwselect": fit.model_info["bwselect"],
-            "bandwidth_selector_gap": (
-                "On this Lee-2008 replica, sp.rdrobust's MSE-RD "
-                "bandwidth selector returns h=0.042 while "
-                "rdrobust::rdrobust returns h=0.176 -- a ~4x "
-                "difference that propagates into the headline "
-                "estimate (sp default ~0.062; R default ~0.076). "
-                "When BOTH implementations are forced to the same "
-                "bandwidth h=0.042287 the bias-corrected point "
-                "estimate matches at rel ~ 1e-6 and the conventional "
-                "SE at ~6%. The discrepancy is in the regularisation "
-                "term of the optimal-bandwidth formula and is "
-                "documented in tests/reference_parity/test_rd_parity.py "
-                "as a known gap; reviewers should treat the bandwidth "
-                "selector as a calibrated package-default choice that "
-                "differs across implementations, not as evidence of a "
-                "numerical bug in the underlying local-polynomial "
-                "estimator."
+            "bandwidth_parity_note": (
+                "Track A uses sp.rdrobust(..., bwselect='cct'), which "
+                "delegates to the official rdrobust Python port and "
+                "matches R/Stata rdrobust default mserd bandwidths on "
+                "the Lee-2008 fixture. The legacy StatsPAI internal "
+                "mserd selector is retained as legacy_internal_mserd_* "
+                "diagnostic rows and should not be used as the "
+                "cross-language default-h parity claim."
             ),
         },
     )

tests/r_parity/07_scm.py

@@ -15,13 +15,13 @@
 52_scm_unique; users who need exact R numbers can call the optional
 ``backend='synth'`` bridge.
 """
+
 from __future__ import annotations
 
 import statspai as sp
 
 from _common import ParityRecord, dump_csv, write_results
 
-
 MODULE = "07_scm"
 
 
@@ -53,13 +53,17 @@ def main() -> None:
 
     rows: list[ParityRecord] = [
         ParityRecord(
-            module=MODULE, side="py", statistic="avg_post_gap",
+            module=MODULE,
+            side="py",
+            statistic="avg_post_gap",
             estimate=float(fit.estimate),
             se=float(fit.se),
             n=int(len(df)),
         ),
         ParityRecord(
-            module=MODULE, side="py", statistic="pre_treatment_rmse",
+            module=MODULE,
+            side="py",
+            statistic="pre_treatment_rmse",
             estimate=float(fit.model_info["pre_treatment_rmse"]),
             n=int(len(df)),
         ),
@@ -77,15 +81,18 @@ def main() -> None:
     for unit in donor_pool:
         rows.append(
             ParityRecord(
-                module=MODULE, side="py",
+                module=MODULE,
+                side="py",
                 statistic=f"weight_{unit}",
                 estimate=float(weight_map.get(unit, 0.0)),
                 n=int(len(df)),
             )
         )
 
     write_results(
-        MODULE, "py", rows,
+        MODULE,
+        "py",
+        rows,
         extra={
             "method": "classic",
             "backend": fit.model_info.get("backend", "native"),
@@ -96,15 +103,29 @@ def main() -> None:
             "n_donors": int(fit.model_info["n_donors"]),
             "placebo": False,
             "tier": "T4",
+            "solver_best_start": fit.model_info.get("solver_best_start"),
+            "solver_near_best_start_count": int(
+                fit.model_info.get("solver_near_best_start_count", 0)
+            ),
+            "solver_near_best_weight_class_count": int(
+                fit.model_info.get("solver_near_best_weight_class_count", 0)
+            ),
+            "solver_near_best_weight_l1_max": float(
+                fit.model_info.get("solver_near_best_weight_l1_max", 0.0)
+            ),
+            "weight_solution_nonunique": bool(
+                fit.model_info.get("weight_solution_nonunique", False)
+            ),
             "native_note": (
                 "Headline row uses backend='native'. The Basque donor-weight "
                 "solution is not unique under the ADH special-predictor "
-                "nested-V specification: native tracks Stata synth on the "
-                "same CSV, while R Synth and Stata synth choose measurably "
-                "different local optima. Native correctness is separately "
-                "certified on a uniquely identified DGP in module "
-                "52_scm_unique; backend='synth' is available when exact R "
-                "Synth numbers are required."
+                "nested-V specification: deterministic multi-start diagnostics "
+                "find multiple near-best donor-weight classes, native tracks "
+                "Stata synth on the same CSV, and R Synth and Stata synth "
+                "choose measurably different local optima. Native correctness "
+                "is separately certified on a uniquely identified DGP in "
+                "module 52_scm_unique; backend='synth' is available when exact "
+                "R Synth numbers are required."
             ),
         },
     )

tests/r_parity/08_dml.R

@@ -4,10 +4,9 @@
 # DoubleML::DoubleMLPLR with mlr3::regr.lm nuisance learners. Five
 # folds, single repetition. Tolerance: rel < 1e-3.
 #
-# Note: cross-fitting fold splits are seeded so that R sklearn-style
-# KFold and mlr3 ResamplingCV(folds=5) DO NOT produce bit-equal
-# splits. The gap on theta_DML_PLR therefore reflects fold-noise
-# alone (the score function is identical on both sides).
+# The CSV carries a deterministic fold_id column written by the Python
+# side.  Both implementations consume those exact folds through their
+# explicit sample-splitting APIs.
 
 .args <- commandArgs(trailingOnly = FALSE)
 .file_arg <- grep("^--file=", .args, value = TRUE)
@@ -51,6 +50,11 @@ dml_obj <- DoubleML::DoubleMLPLR$new(
   score   = "partialling out",
   dml_procedure = "dml2"
 )
+fold_levels <- sort(unique(dt$fold_id))
+test_ids <- lapply(fold_levels, function(k) which(dt$fold_id == k))
+train_ids <- lapply(fold_levels, function(k) which(dt$fold_id != k))
+smpls <- list(list(train_ids = train_ids, test_ids = test_ids))
+dml_obj$set_sample_splitting(smpls)
 dml_obj$fit()
 
 theta <- as.numeric(dml_obj$coef["educ"])
@@ -74,5 +78,7 @@ write_results(MODULE, rows,
                            n_folds = 5L,
                            ml_g = "regr.lm",
                            ml_m = "regr.lm",
+                           fold_source = "user",
+                           fold_column = "fold_id",
                            score = "partialling out",
                            dml_procedure = "dml2"))

tests/r_parity/08_dml.py

@@ -24,6 +24,7 @@
 
 def main() -> None:
     df = sp.datasets.card_1995()
+    df["fold_id"] = np.arange(len(df)) % 5
     dump_csv(df, MODULE)
 
     # NOTE: pass numpy seed so sklearn's KFold split is deterministic;
@@ -34,6 +35,8 @@ def main() -> None:
         model="plr",
         model_y=LinearRegression(),
         model_d=LinearRegression(),
+        n_folds=5,
+        fold_indices=df["fold_id"].to_numpy(),
     )
 
     rows: list[ParityRecord] = [
@@ -55,18 +58,12 @@ def main() -> None:
             "ml_g": "LinearRegression",
             "ml_m": "LinearRegression",
             "seed": PARITY_SEED,
-            "fold_split_note": (
-                "DML PLR with linear-regression nuisance learners is "
-                "equivalent to OLS on the partialled-out residuals, so "
-                "both implementations should converge to OLS beta(educ) "
-                "= 0.10999 (Module 01) up to fold-split Monte Carlo "
-                "noise. sklearn.KFold(shuffle=True, seed=42) and mlr3 "
-                "ResamplingCV(folds=5, set.seed(42)) produce different "
-                "partitions, so the residual gap (~2e-4 absolute, "
-                "~2.4e-3 relative) is fold-split RNG drift, not a "
-                "score-function or orthogonalisation bug. Forcing the "
-                "same fold indices on both sides closes the gap to "
-                "machine precision."
+            "fold_source": str(fit.model_info["fold_source"]),
+            "fold_column": "fold_id",
+            "fold_parity_note": (
+                "Python and R read the same fold_id column and use "
+                "explicit sample-splitting APIs before fitting the "
+                "linear-nuisance PLR score."
             ),
         },
     )

tests/r_parity/10_honest_did.R

@@ -4,7 +4,9 @@
 # runs HonestDiD::createSensitivityResults under "FLCI" (the
 # finite-sample length-optimised CI under the smoothness restriction
 # Delta^SD), which is the closest match to sp.honest_did's smoothness
-# restriction. Tolerance: abs < 0.05 on the CI bounds.
+# restriction. Tolerance: abs < 5e-4 on CI bounds; Python-R is
+# machine-equal through the R backend, while the Stata port is kept
+# inside this package-port envelope.
 
 .args <- commandArgs(trailingOnly = FALSE)
 .file_arg <- grep("^--file=", .args, value = TRUE)