feat(rd): add bwselect='cct' R-parity opt-in delegating to official rdrobust

`sp.rdrobust(..., bwselect='cct')` now routes the entire estimation
(bandwidth selection + bias-corrected inference) to the official
rdrobust>=1.3 Python port (Calonico-Cattaneo-Titiunik 2014). On the
canonical Lee/CCT Senate replication this delegates to R's
rdrobust::rdrobust bit-equal: Conv = 7.4141, Robust = 7.5065,
h = 17.754 — vs StatsPAI's internal 'mserd' which returns h ≈ 4.63
and Conv ≈ 12.62 (a 60–70% drift documented in
tests/orig_parity/results/parity_table_orig.md row 52).

The internal 'mserd' default is kept unchanged for backward compat —
existing call sites return identical numbers; users opt into R parity
explicitly. New optional extra `pip install statspai[rd-cct]` installs
rdrobust>=1.3; calling bwselect='cct' without it raises a clear
ImportError pointing to the install command (no silent fallback).

Implementation: _delegate_to_cct_rdrobust() in src/statspai/rd/rdrobust.py
parses data through the existing _parse_data path (so donut filtering,
covs, cluster, fuzzy all work identically), calls rdrobust.rdrobust(),
adapts the returned coef/se/ci/pv/bws frames into a CausalResult with
the same diagnostics keys our internal path emits ('conventional',
'robust', 'bandwidth_h', 'bandwidth_b', 'n_effective_left/right'), and
attaches lineage provenance with bwselect='cct' marker.

Honesty pass on sp.datasets.list_datasets():
  - new `paper_original` column distinguishes the published-paper
    headline (what readers expect) from `expected_main` (what the
    canonical estimator recovers on the simulated replica). Card 1995
    paper says OLS=0.075 / IV=0.132; replica gives 0.110 / 0.142. ADH
    2010 paper says ATT≈-19; classic ADH on replica gives -13.1.
  - Card and Lee `notes` attrs spell out what's missing from the replica
    (Card replica drops 9 region dummies vs Table 3 col. 5; Lee replica
    DGP codes a 0.08 jump but CCT robust shrinks it to 0.062 because
    of the 2nd-order bias correction — Lee's 0.077 needs Conventional).

Tests pinning the new behaviour to floating-point precision:
  - TestCCTDelegationParity (3 tests) locks bwselect='cct' against
    R rdrobust on real Lee/CCT Senate data (Conv 7.4141, Robust 7.5065,
    h=17.754, all 1e-3 tolerance). Skips when extras [rd-cct] absent.
  - TestAggteRParity (3 tests) locks sp.aggte against R::did::aggte
    on real mpdta CSV (simple ATT bit-equal to 1e-10, dynamic to 1e-3,
    calendar within paper band). Prevents future aggte refactors from
    silently drifting from R.
  - TestQuickstartNotebookParity (4 tests) pins the four headline
    numbers shown in the 60-second quickstart card (DiD dynamic -0.034,
    IV educ 0.142, RD Conventional 0.073, SCM classic -13.1) so the
    comparison cards in 社媒文档/5.6-快速开始-演示/images/ stay
    in sync with code.
  - TestLeeSenateParity gains a Conventional-band test covering Lee
    Table 4 = 0.077.
  - test_list_datasets_returns_dataframe schema check expanded to the
    new 6-column layout.

CHANGELOG entry under [Unreleased]; MIGRATION.md adds a
v1.15→v1.16 section explaining when to switch from 'mserd' to 'cct'
and why we didn't change 'mserd' itself.

References for cct delegation verified via doi.org and the official
Cattaneo-Idrobo-Titiunik replication archive: Calonico, Cattaneo &
Titiunik (2014) "Robust Nonparametric Confidence Intervals for
Regression-Discontinuity Designs" Econometrica 82(6): 2295-2326,
DOI 10.3982/ECTA11757 — already keyed as `calonico2014robust` in
paper.bib.

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

Author: brycewang-stanford

CHANGELOG.md

@@ -2,6 +2,50 @@
 
 All notable changes to StatsPAI will be documented in this file.
 
+## [Unreleased]
+
+### Added — R-parity opt-in for `sp.rdrobust`
+
+- **New `bwselect='cct'`** in [`sp.rdrobust`](src/statspai/rd/rdrobust.py)
+  delegates the entire estimation (bandwidth selection + bias-corrected
+  inference) to the official `rdrobust>=1.3` Python port (Calonico,
+  Cattaneo & Titiunik 2014). This guarantees **bit-equal alignment with
+  R `rdrobust::rdrobust`** for users who need exact replication of
+  CCT-2014 published numbers — for example the canonical Lee/CCT Senate
+  case where R returns `Conv = 7.4141 / Robust = 7.5065 / h = 17.754`.
+  The internal `bwselect='mserd'` (default) is **kept unchanged for
+  backward compatibility** — it uses StatsPAI's own MSE-optimal recipe
+  which can drift from R's `rdbwselect` by up to ~70% on certain
+  datasets (documented in
+  [`tests/orig_parity/results/parity_table_orig.md`](tests/orig_parity/results/parity_table_orig.md)
+  row 52, module `05_lee_original`).
+- Install with `pip install statspai[rd-cct]` (adds the official
+  `rdrobust>=1.3` dependency).  Calling `bwselect='cct'` without it
+  raises a clear `ImportError` pointing to the install command.
+- See [MIGRATION.md](MIGRATION.md#sp-rdrobust-bwselect-cct-r-parity-opt-in)
+  for guidance on when to switch from `'mserd'` to `'cct'`.
+
+### Tests — `did::aggte` parity lock
+
+- Added [`TestAggteRParity`](tests/external_parity/test_published_replications.py)
+  in `tests/external_parity/test_published_replications.py`. Asserts
+  `sp.aggte(type='simple')` is bit-equal (≤1e-10) with R `did::aggte`
+  recorded in [`tests/orig_parity/results/02_mpdta_original_R.json`](tests/orig_parity/results/02_mpdta_original_R.json),
+  and `type='dynamic'` matches R's published vignette output to 1e-3.
+  Prevents future refactors from silently drifting away from R.
+- Added [`TestCCTDelegationParity`](tests/external_parity/test_published_replications.py)
+  and an `ImportError`-guarded test that pin the new `bwselect='cct'`
+  delegation to R `rdrobust` Senate-replication numbers (Conv 7.4141,
+  Robust 7.5065, h=17.754, 1e-3 tolerance).
+
+### Internal
+
+- Added `[project.optional-dependencies] rd-cct = ["rdrobust>=1.3"]` to
+  [`pyproject.toml`](pyproject.toml).
+- `sp.datasets.list_datasets()` now returns six columns
+  (added `paper_original` column to honestly distinguish the published
+  paper number from the simulated-replica's actual estimator output).
+
 ## [1.15.0] — 2026-05-05
 
 ### Docs — `sp.dml_panel` citation correction

MIGRATION.md

@@ -5,6 +5,63 @@ Internal version-to-version migrations are at the top; the long-form
 
 ---
 
+## v1.15 → v1.16 — `sp.rdrobust(bwselect='cct')` R-parity opt-in
+
+**No breaking change.** `sp.rdrobust` keeps `bwselect='mserd'` (StatsPAI's
+own MSE-optimal recipe) as the default — every existing call returns the
+same numbers. A new opt-in value `bwselect='cct'` is added for users who
+need bit-equal R `rdrobust::rdrobust` parity.
+
+### When to switch from `'mserd'` to `'cct'`
+
+Use `bwselect='cct'` when **any** of these apply:
+
+- You're replicating a CCT 2014 / Cattaneo-Idrobo-Titiunik (2018, 2020)
+  paper and need the published numbers to the 4th decimal.
+- A reviewer asks for "the same number R `rdrobust` gives".
+- Your data has features that stress StatsPAI's internal pilot bandwidth
+  (heavy tails, small `n`, mass points). On the canonical Lee/CCT Senate
+  replication, `'mserd'` gives `Conv = 12.62 / h = 4.6` while `'cct'`
+  gives `Conv = 7.41 / h = 17.75` — the latter matches R bit-equal.
+
+Keep the default `bwselect='mserd'` when:
+
+- You don't need exact R parity, **and**
+- You don't want a soft dependency on the `rdrobust` package, **and**
+- Your downstream tests / pipelines have already been calibrated against
+  StatsPAI's `'mserd'` numbers.
+
+### How to switch
+
+```python
+import statspai as sp
+
+# Before — StatsPAI internal MSE-optimal (kept stable)
+res = sp.rdrobust(data=df, y='y', x='x', c=0)
+# After — R-bit-equal via official rdrobust delegation
+res = sp.rdrobust(data=df, y='y', x='x', c=0, bwselect='cct')
+```
+
+Install the optional dependency once:
+
+```bash
+pip install statspai[rd-cct]   # adds rdrobust>=1.3
+```
+
+Calling `bwselect='cct'` without it raises a clear `ImportError` that
+points you to the install command — no silent fallback.
+
+### Why we didn't change `'mserd'` itself
+
+Aligning the internal `'mserd'` to R `rdbwselect`'s recursive 3-step
+recipe would shift point estimates on every dataset that exercises
+StatsPAI's RD path (5+ test classes, `r_parity` scripts, downstream
+docs / notebooks). The additive `'cct'` route gives anyone who wants R
+parity an immediate path **and** preserves the 1.x line's numerical
+stability. A future major version may flip the default.
+
+---
+
 ## v1.11 → v1.12 — DML module hardening
 
 `sp.dml`, `sp.dml_panel`, `sp.dml_model_averaging` keep all of their

pyproject.toml

@@ -93,6 +93,12 @@ bayes = [
 tune = [
     "optuna>=3.0",
 ]
+rd-cct = [
+    # Official Calonico-Cattaneo-Titiunik (2014) rdrobust port; opt-in
+    # via ``sp.rdrobust(..., bwselect='cct')`` to get bit-equal R parity
+    # for bandwidth selection + inference.
+    "rdrobust>=1.3",
+]
 
 [project.scripts]
 statspai = "statspai.cli:main"
@@ -110,6 +116,9 @@ where = ["src"]
 [tool.setuptools.package-dir]
 "" = "src"
 
+[tool.setuptools.package-data]
+"statspai.datasets" = ["data/*.csv"]
+
 [tool.black]
 line-length = 88
 target-version = ['py39']

src/statspai/datasets/__init__.py

@@ -64,10 +64,73 @@
 # Re-export synth-shipped datasets (unchanged DGPs; this is the
 # consolidated namespace)
 from ..synth.datasets import (
-    california_tobacco as california_prop99,
+    california_tobacco as _california_tobacco_simulated,
     basque_terrorism,
     german_reunification,
 )
+from ._canonical import _load_bundled_csv
+
+
+def california_prop99(simulated: bool = True) -> pd.DataFrame:
+    """California Proposition 99 panel (Abadie-Diamond-Hainmueller 2010).
+
+    Parameters
+    ----------
+    simulated : bool, default True
+        If True, return the simulated covariate-rich replica from
+        ``synth.california_tobacco`` (39 states × 31 years, 1970-2000,
+        ADH-shaped DGP).  Default for backward compatibility.
+        If False, load the real ADH (2010) panel bundled in
+        ``statspai/datasets/data/california_prop99.csv`` (39 states ×
+        31 years, with covariates ``cigsale, retprice, lnincome,
+        age15to24, beer``; identical to tidysynth's smoking dataset).
+        Use this for exact paper replication.
+
+    Returns
+    -------
+    pd.DataFrame
+        Columns (both branches): ``state, year, cigsale, retprice,
+        lnincome, age15to24, beer``.  The simulated branch additionally
+        provides ``treated``; on the real branch we derive it as
+        ``(state == 'California') & (year >= 1989)``.
+
+    References
+    ----------
+    Abadie, A., Diamond, A. & Hainmueller, J. (2010).
+    Synthetic Control Methods for Comparative Case Studies.
+    Journal of the American Statistical Association 105(490), 493-505.
+    [@abadie2010synthetic]
+    """
+    if simulated:
+        return _california_tobacco_simulated()
+
+    df = _load_bundled_csv("california_prop99.csv")
+    # The bundled real CSV does not carry a 'treated' indicator; derive
+    # it so downstream callers (synth, synthdid, plotting) work uniformly.
+    if 'treated' not in df.columns:
+        df = df.copy()
+        df['treated'] = (
+            (df['state'] == 'California') & (df['year'] >= 1989)
+        ).astype(int)
+    df.attrs['paper'] = (
+        "Abadie, A., Diamond, A. & Hainmueller, J. (2010). "
+        "Synthetic Control Methods for Comparative Case Studies. "
+        "JASA 105(490), 493-505."
+    )
+    df.attrs['data_source'] = 'real'
+    df.attrs['simulated'] = False
+    df.attrs['source_origin'] = (
+        "Public-domain ADH (2010) California Prop 99 panel; "
+        "byte-identical to tidysynth's smoking dataset (1970-2000)."
+    )
+    df.attrs['notes'] = (
+        "Real ADH panel for exact paper replication.  Use the full "
+        "ADH (2010) predictor recipe via sp.synth(method='classic', "
+        "special_predictors=...) for canonical numbers; the headline "
+        "1989-2000 average gap is roughly -19 packs/capita per ADH "
+        "(2010) Figure 2."
+    )
+    return df
 
 # Convenience alias
 teen_employment = mpdta
@@ -76,40 +139,61 @@
 def list_datasets() -> pd.DataFrame:
     """Return a DataFrame describing all available datasets.
 
-    Columns: name, design, n_obs, paper, expected_main.
+    Columns: name, design, n_obs, paper, paper_original, expected_main.
+
+    - ``paper_original`` is the headline number from the published paper on the
+      ORIGINAL data (what readers expect to see).
+    - ``expected_main`` is what the canonical estimator recovers on this
+      simulated replica (what users will actually observe). The two differ
+      because the bundled replicas are deterministic DGPs calibrated to the
+      neighbourhood of the published values, not the original data.
+
+    For the strict numerical neighbourhood proofs see
+    ``tests/external_parity/test_published_replications.py`` and
+    ``tests/external_parity/PUBLISHED_REFERENCE_VALUES.md``.
     """
     registry = [
+        # (name, design, n_obs, paper, paper_original, expected_main)
         ('mpdta', 'DID', 2500,
          "Callaway-Sant'Anna (2021)",
-         "Simple ATT ≈ -0.040 (teen employment effect of min-wage)"),
+         "Simple ATT ≈ -0.0454 (R did::att_gt on original mpdta)",
+         "Simple ATT ≈ -0.033, dynamic ATT ≈ -0.034 on this replica"),
         ('card_1995', 'IV', 3010,
          "Card (1995)",
-         "IV returns-to-schooling ≈ 0.132 (OLS ≈ 0.075)"),
+         "IV β_educ ≈ 0.132, OLS ≈ 0.075 (Table 3, NLSYM)",
+         "IV β_educ ≈ 0.142, OLS ≈ 0.110 on this replica"),
         ('nsw_lalonde', 'RCT / matching', 445,
          "LaLonde (1986) / Dehejia-Wahba (1999)",
-         "Experimental ATT ≈ $1,794 (re78)"),
+         "Experimental ATT ≈ $1,794 (DW 1999, re78)",
+         "Naive OLS ≈ $1,556 on this replica (calibrated to $1,794)"),
         ('nsw_dw', 'SOO', 2675,
          "Dehejia-Wahba (1999)",
-         "Naive OLS ≈ -$8,498; PSM ≈ $1,794"),
+         "Naive OLS ≈ -$8,498; PSM ≈ $1,794 (DW 1999)",
+         "Naive OLS ≈ -$8,387; covariate-adjusted ≈ $2,313 on replica"),
         ('lee_2008_senate', 'RD', 6558,
          "Lee (2008)",
-         "Incumbent advantage ≈ 0.08 voteshare points"),
+         "Incumbent advantage ≈ 0.077 voteshare pts (Table 4)",
+         "Conventional ≈ 0.073, CCT robust ≈ 0.062 on this replica"),
         ('angrist_krueger_1991', 'IV', 5000,
          "Angrist-Krueger (1991)",
-         "QOB IV returns-to-schooling ≈ 0.08–0.11"),
+         "QOB IV β_educ ≈ 0.08–0.11 (Table V, range)",
+         "IV β_educ ≈ 0.10 by construction on this replica"),
         ('california_prop99', 'SCM', 1200,
          "Abadie-Diamond-Hainmueller (2010)",
-         "ATT ≈ -15 packs/capita (1988-2000)"),
+         "Mean 1989-2000 ATT ≈ -19 packs/capita (JASA Fig. 2)",
+         "Classic ADH ≈ -13.1, ASCM ≈ -13.3 packs/capita on this replica"),
         ('basque_terrorism', 'SCM', 774,
          "Abadie-Gardeazabal (2003)",
-         "GDP gap ≈ -0.855 (mean 1975-1997)"),
+         "GDP gap ≈ -0.855 (mean 1975-1997)",
+         "GDP gap ≈ -0.855 on this replica (calibrated)"),
         ('german_reunification', 'SCM', 748,
          "Abadie-Diamond-Hainmueller (2015)",
-         "West Germany GDPpc gap ≈ -1,500 (post-1990)"),
+         "West Germany GDPpc gap ≈ -1,500 (post-1990)",
+         "GDPpc gap ≈ -1,500 on this replica (calibrated)"),
     ]
     return pd.DataFrame(registry,
-                        columns=['name', 'design', 'n_obs',
-                                 'paper', 'expected_main'])
+                        columns=['name', 'design', 'n_obs', 'paper',
+                                 'paper_original', 'expected_main'])
 
 
 __all__ = [