StatsPAI 1.17.0 — correctness fixes + agent-UX hardening

Added

• Registry-example bind guard (tests/test_registry_examples_bind.py).
  A parametrized test now statically parses every registered example string
  and binds the keyword arguments of its sp.<name>(...) call against the real
  signature, failing if an example references a keyword the function does not
  accept or does not parse. This locks down the agent copy-paste path — an
  agent reads sp.describe_function(name) and runs the example verbatim — so
  the registry/example drift fixed below cannot silently return. 373 examples
  bind green.

• parity optional extra — opt-in DoubleML reference pin for sp.dml.
  pip install -e ".[dev,parity]" now installs doubleml-for-py (the Python
  DoubleML reference of Bach, Chernozhukov, Kurz & Spindler, JMLR 23(53),
  2022), so tests/external_parity/test_dml_python_parity.py runs instead
  of silently skipping. Under identical scikit-learn learners and folds,
  sp.dml(model='plr') reproduces doubleml-for-py to machine precision on
  the seed-42 fixture — |Δ coefficient| = 1.1e-16 and |Δ standard error| = 1.4e-17, i.e. one float64 unit in the last place. doubleml remains not
  a runtime dependency. The measured numbers, software versions, and the
  divergence discussion are recorded in a new Double Machine Learning Parity
  section of docs/joss_validation_dossier.md, with a one-command reproduce
  path added to docs/joss_reviewer_guide.md. Verified by installing the
  extra and running both tests/external_parity/test_dml_python_parity.py and
  tests/reference_parity/test_dml_parity.py (55 DML tests green).

Fixed

• ⚠️ Correctness: sp.drdid(method='trad') returned ~half the true ATT.
  The traditional doubly-robust DiD branch of sp.drdid (Sant'Anna & Zhao
  2020) normalised each of its four cell terms (treated/control × post/pre) by
  the full sample size n instead of by that cell's weight mass. This
  multiplied every term by the cell's sample share (~0.25 each on a balanced
  2×2), biasing the ATT toward zero by ~50%: on a 2×2 with a true ATT of 2.0
  (raw DiD 1.96) the traditional estimator returned ≈1.04. Each term is now
  normalised by its own weight total, so method='trad' reduces exactly to the
  raw 2×2 DiD when no covariates are supplied and recovers the true ATT with
  covariates. The improved (locally efficient) method='imp' — the default
  — already normalised correctly and is unchanged, so no default-path or
  parity/dossier numbers move. sp.drdid now also raises ValueError on an
  unknown method instead of silently treating it as traditional (previously
  e.g. method='ipw' ran the traditional branch). See MIGRATION.md.

• ⚠️ Correctness: sp.multiway_cluster_vcov undercounted intersection
  clusters, biasing multiway-cluster-robust standard errors. The
  Cameron-Gelbach-Miller inclusion-exclusion builds an intersection cluster
  (the unique combinations of the clustering dimensions); its key was formed by
  joining the dimensions into a single string with a "\0" separator, but
  NumPy fixed-width unicode strips the embedded NUL, so e.g. (1, 23) and
  (12, 3) both collapsed to "123". On a 40×50 crossed-cluster DGP this
  merged 1733 true intersection clusters into 1639, biasing the two-way SE by
  ~0.2% (~0.5% at three-way) versus the canonical estimator. The intersection
  key is now built collision-free via np.unique(axis=0) on per-dimension
  integer codes. sp.multiway_cluster_vcov now reproduces sandwich::vcovCL
  and sp.twoway_cluster to machine precision (two-way exact; three-way rel
  ~4e-7), pinned by the new tests/r_parity module 56 and a direct
  twoway-vs-multiway regression test. Propagates to did.harvest and
  panel.feols multiway-clustered SEs. sp.twoway_cluster itself was already
  correct (distinct collision-free key) and is unchanged. See MIGRATION.md.

• .glance() crashed (OverflowError: cannot convert float infinity to integer) on Cox and parametric-survival (survreg) results. Those
  estimators deliberately store df_resid = inf to signal a large-sample
  (normal) reference distribution, but glance() cast the residual degrees of
  freedom with int() unconditionally. The cast now passes non-finite degrees
  of freedom through unchanged; finite results keep an integer df_resid
  (no change). A crash-hunt across ~48 fitted results confirmed the rest of
  the §3 unified-result export surface (summary/to_latex/to_markdown/
  to_word/to_excel/cite/tidy/for_agent/plot) is otherwise clean
  across 20+ estimator families. Covered by tests/test_glance_survival.py.

• sp.event_study results crashed the library's own exporters, plotters and
  pre-trend tools (canonical-column mismatch). sp.event_study emitted its
  coefficient table under the column name estimate, but the rest of the DID
  family — and every downstream consumer — keys on the canonical att column
  (did._core.EVENT_STUDY_COLUMNS). So the canonical event-study estimator
  was incompatible with its own tooling: .tidy() (and the .to_markdown() /
  .to_excel() / .to_word() exporters that delegate to it) raised
  TypeError: unsupported operand type(s) for /: 'NoneType' and 'float';
  .plot() / .event_study_plot() / sp.enhanced_event_study_plot raised
  KeyError: 'att'; and sp.honest_did / sp.breakdown_m raised
  ValueError: missing {'att'}. The event-study table now carries the
  canonical att column (with estimate retained as a backward-compatible
  alias), fixing every consumer at the source. No numerical change.

• sp.pretrends_test / sp.pretrends_summary crashed (LinAlgError: Singular matrix) on every standard sp.event_study result — the same
  reference-period defect already fixed in pretrends_power: the SE = 0
  omitted period made the diagonal VCV singular. It is now dropped before
  inversion, with a clear ValueError on a genuinely collinear pre-period set.

• sp.diagnose_result crashed (TypeError: bad operand type for abs(): 'str') on sp.synth results. The donor-pool check iterated the synthetic
  weights, but sp.synth stores them as a ['unit', 'weight'] DataFrame, so
  iteration yielded column-name strings. The weights are now coerced to their
  numeric values regardless of container (DataFrame / Series / dict / array).

• sp.pretrends_power crashed (LinAlgError: Singular matrix) on every
  standard sp.event_study result. Roth's (2022) pre-trend power calculation
  inverts the pre-period variance–covariance matrix, but the omitted reference
  period (relative time −1) is reported with a standard error of exactly zero,
  so the diagonal VCV was singular and np.linalg.inv raised on the exact
  workflow shown in the function's own docstring. The reference period (and any
  other mechanically-normalised, zero-SE period) is now dropped before
  inversion — it is the baseline, not an estimated coefficient — so the joint
  pre-trend test runs on the estimated pre-periods only. A full-rank
  model_info['vcv_pre'] and a full-length delta are aligned to the retained
  periods, and a still-singular VCV now raises a clear ValueError (collinear
  pre-periods) instead of an opaque NumPy error. No output changes for any call
  that previously succeeded. Covered by tests/test_pretrends_power.py.

• ⚠️ Correctness — sp.structural_break sup-F p-value used the wrong null
  distribution. The Chow/sup-F statistic is a supremum of the F statistic
  over candidate break points, so under H0 it follows the Andrews (1993)
  sup-F law — not F(k, n-2k). The previous code referred the maximised
  statistic to the ordinary F CDF, which ignored the maximisation and
  massively over-rejected: on pure Gaussian white noise at the 5% level the
  test flagged a spurious structural break in 33–37% of series (measured,
  n ∈ {100, 200, 400}). The p-value is now computed from the Andrews (1993)
  limiting null — a q-vector Brownian-bridge functional sampled by a
  deterministic, cached simulation on a grid tied to the sample size —
  restoring nominal size (~0.05) while retaining power (1.00 / 0.88 to
  detect a one-/half-σ mean shift at n=200). The same correct threshold now
  drives the Bai-Perron sequential supF(l+1|l) stopping rule (previously the
  same naive-F over-detection), so method='bai-perron' no longer
  over-segments noise. As a side benefit the Bai-Perron result now populates
  f_stats / p_values (one sup-F statistic and Andrews p-value per detected
  break, chronologically aligned) instead of returning None. Reference
  verified via Crossref / Econometric Society / RePEc: Andrews, D.W.K. (1993),
  Econometrica 61(4), 821-856, doi:10.2307/2951764. See MIGRATION.md.

• 33 registered example strings were statically broken (agent-UX). Six
  failed to parse (stray/unmatched parens, a positional-after-keyword
  ... placeholder, an unclosed call) and 27 passed a keyword the function
  does not accept — a deterministic TypeError/SyntaxError on the exact
  agent copy-paste path. Root causes were two long-standing parameter-name
  drifts: the Mendelian-randomization family (mr_egger/mr_ivw/mr_raps/
  mr_presso) used the short b_exp/b_out/se_exp/se_out names instead
  of the implemented beta_exposure/beta_outcome/se_exposure/se_outcome,
  and the Bayesian family (bayes_rd/bayes_fuzzy_rd/bayes_mte/bayes_did/
  bayes_iv) plus metalearner/tmle/causal_impact/sensemakr/
  spec_curve/qreg/tobit/heckman/cluster_cross_interference/
  causal_dqn/pci_mtp/bartik/ffl_decompose drifted from their
  signatures. qreg/tobit/heckman/spec_curve additionally had wrong
  call shapes (formula passed positionally into data; flat controls
  where a list-of-lists was required) and were rebuilt to runnable form and
  executed to confirm. Fixed acros...

StatsPAI v1.16.1 — synth default restored to classic SCM + simplex weight projection (⚠️ correctness)

StatsPAI 1.16.1 is a patch release with two ⚠️ correctness fixes to the synthetic-control path, a Python 3.10+ agent-schema fix, and a reviewer-doc refresh. Re-run any affected synthetic-control analyses.

📦 PyPI: https://pypi.org/project/StatsPAI/1.16.1/ · pip install -U StatsPAI

⚠️ Correctness — sp.synth() default restored to canonical classic SCM

The bare sp.synth(...) entry point again defaults to method='classic' (Abadie, Diamond & Hainmueller 2010), so a default call returns convex, non-negative, sum-to-one donor weights. The signature default had silently drifted to method='augmented' (Augmented SCM, Ben-Michael, Feller & Rothstein 2021), whose ridge correction is designed to allow negative donor weights by extrapolating outside the donor convex hull — surprising for the canonical synth() entry point and inconsistent with the docstring, the migration-from-R mapping, and the Prop99 examples. Augmented SCM remains fully available via method='augmented' (or 'ascm'); every non-default method is unchanged.

⚠️ Correctness — synthetic-control weights projected back onto the simplex

solve_simplex_weights — the inner W solver shared by sp.synth and the SCM/sdid/augsynth/gsynth family — now projects the SLSQP solution back onto the unit simplex (clipping sub-tolerance negative weights to zero and renormalising to sum 1). SLSQP enforces w_j ≥ 0, Σw = 1 only up to its own tolerance, so the raw solution could carry small negative donor weights (observed down to ≈ −7.5e-4). Donor weights change by the solver's sub-tolerance noise only; the projection moves the native output toward the reference clean-simplex solution, so R Synth / gsynth parity is preserved.

Fixed — agent schema generation preserves full typing shapes

sp.function_schema / the registry schema generator now keep parametrised typing annotations intact across Python 3.9–3.13. registry._stringify_annotation previously collapsed aliases such as Optional[Dict[str, Any]] to the bare origin name (Optional, Dict) on Python 3.10+ because those aliases expose __name__; it now resolves typing.-prefixed and __origin__-bearing annotations first. Machine-readable parameter shapes are now stable and version-independent. No estimator numbers change.

Docs

Reviewer-facing validation docs refreshed (focused regression-test command, compatibility test-path pointer, 2026-06-01 activity/measurement dates). Live docs/stats.md counts re-measured against the 1.16.1 source tree (source 269,010 LOC).

Verification

• Full suite: 5848 passed, 46 skipped, 1 xfailed (exit 0)
• Synth reference parity: 3/3 passed · registry drift check: OK · twine check: PASSED
• Clean-venv install from PyPI verified: 1020 registered functions

Full notes in CHANGELOG.md under [1.16.1].

StatsPAI v1.16.0 — Arellano–Bond GMM + qreg SE correctness fixes; parity harness → 50 R / 43 Stata modules

⚠️ Correctness fix

• sp.qreg Powell sandwich SE was wrong by a factor of √n — every
  pre-fix p-value, z-statistic, and confidence interval emitted by
  sp.qreg was unusable. The closed-form Koenker (2005, eq. 3.7) iid
  kernel sandwich is V = τ(1−τ) / f̂(0)² · (X'X)⁻¹. _qreg_se had an
  extra factor of n in the denominator (/ (n * f0**2)), so the
  reported SE was the correct SE divided by √n — on n = 500 the SE
  was ~20× too small. The fix removes the spurious n. After the fix
  the three-way parity at the median tolerance (tests/r_parity/40_qreg)
  matches quantreg::rq within 1.4–6.8 % and Stata qreg within 2.9 %,
  consistent with the documented kernel-vs-Koenker-Bassett SE method
  gap. Action: any analysis that previously used sp.qreg SE,
  z-statistic, p-value, or CI must be re-run; point estimates are
  unaffected. See MIGRATION.md § sp-qreg-se-fix
  for the per-call impact and rerun recipe.

• sp.xtabond (Arellano-Bond difference GMM) point estimates AND SEs
  were wrong — finding #12. The estimator built a flat, fixed set of
  lagged-level instrument columns (gmm_lags=(2,5)) and then dropped
  every row missing any of them, which on a short panel discards most of
  the sample; it also used W = (Z'Z)⁻¹ as the one-step weight. The
  correct Arellano-Bond estimator uses a block-diagonal GMM
  instrument matrix (every available deeper lag is a period-specific
  moment, missing lags filled with 0, no rows dropped) and the
  one-step weight W = (Σᵢ Zᵢ'H Zᵢ)⁻¹ where H carries the MA(1)
  structure of the differenced errors (2 on the diagonal, −1 on the
  first off-diagonals). On the parity DGP the old code gave
  β_{y₋₁}=0.264 (se 0.224) vs Stata's 0.391 (se 0.046) — a 48 %
  estimate gap and an 80 % SE gap. After the rewrite the one-step robust
  estimates match Stata's xtabond y x, lags(1) vce(robust) to machine
  precision (tests/r_parity/50_xtabond, rel ≈ 1e-15 on both β and
  SE). The default gmm_lags is now (2, None) (all available deeper
  lags, matching Stata's default; pass an explicit max to cap). Two-step
  GMM now applies the Windmeijer (2005) finite-sample SE correction.
  Action: re-run any analysis that used sp.xtabond — both point
  estimates and SEs change. See
  MIGRATION.md § sp-xtabond-fix.

• sp.xtabond(method='system') / sp.panel(method='system') now raise
  NotImplementedError instead of returning an unvalidated (and, after
  the difference-GMM rewrite, badly distorted) estimate. Proper
  Blundell-Bond system GMM requires a stacked level equation and its own
  Stata xtdpdsys parity reference, which is planned for a future
  release. Action: use method='difference' (Arellano-Bond), now
  validated to machine precision.

Added — Parity coverage expansion (2026-05-28 session)

• 15 net-new parity modules (tests/r_parity/{37–51}_*) covering
  sp.ppmlhdfe, sp.drdid, sp.arima, sp.qreg, sp.tobit,
  sp.nbreg, sp.heckman, sp.mlogit, sp.ologit, sp.clogit,
  sp.probit, sp.oprobit, sp.xtabond, sp.newey, and a 3-FE PPML
  variant. The 3-way Track A table
  (tests/r_parity/results/parity_table_3way.md) now covers 50
  R-joined modules versus 36 previously, with a Stata reference for 43
  versus 21 (50_xtabond is a Py-Stata-only migration check omitted
  from the R-joined table). The expansion surfaced the qreg and newey
  SE fixes above and further P1/P2 findings recorded in
  tests/r_parity/PARITY_SESSION_2026-05-28.md.

Fixed

• Cleaned up JOSS review follow-ups: removed two uncited duplicate BibTeX
  entries that caused editorialbot DOI suggestions, aligned the AKM
  shift-share citation key / DOI metadata, and refreshed v1.15.6 wording in
  reviewer-facing docs and README release callouts.
• tools/audit_citations.py now treats transient HTTP/socket/SSL timeouts as
  unresolved citation lookups instead of leaking Python tracebacks.
• tests/r_parity/36_mediation.py referenced model_info["n_boot"],
  but sp.mediation's schema renamed this to n_boot_requested /
  n_boot_successful / n_boot_failed. The parity script crashed
  before producing JSON; pinned it to the new key.

StatsPAI 1.15.6 — Scott Rozelle co-author + JOSS readiness

[1.15.6] — 2026-05-24

Changed — Co-authorship, JOSS submission readiness

• Added Scott Rozelle as co-author across all package metadata:
  pyproject.toml, src/statspai/__init__.py (__author__), CITATION.cff,
  .zenodo.json, mkdocs.yml, the package citation templates in
  src/statspai/_citation.py (BibTeX / APA / plain), and the README BibTeX
  snippets (English and Chinese).
• ⚠️ Downstream-facing rename: unified the package BibTeX key to
  wang2026statspai (CLAUDE.md §10 lastnameYEARkeyword convention).
  Previous keys emitted or documented in earlier versions
  (wang_statspai_2026, wang_rozelle_statspai_2026, statspai2026software,
  bare statspai) are removed in favor of a single canonical key. Downstream
  .tex files that cite the previous key need a one-line rename to
  \cite{wang2026statspai}. The impact surface is small — only users who
  literally copied the previous BibTeX entry into their own .bib are
  affected; users who regenerate via sp.citation("bibtex") get the new key
  automatically.
• sp.citation("bibtex") now emits the unified key and the updated author
  list. sp.citation("apa") and sp.citation("plain") already reflected the
  co-author; both surfaces now also carry the 1.15.6 version string.
• CITATION.cff version / date-released bumped to 1.15.6 / 2026-05-24.

Added — JOSS reviewer-facing documentation

• docs/joss_reviewer_guide.md — install, smoke test, representative offline
  examples, targeted tests, and build check, intended as a short reviewer
  path. All five smoke-test API calls are verified against the current
  registry (ivreg, callaway_santanna + aggte, rdrobust, synth,
  describe_function / function_schema).
• docs/joss_validation_dossier.md — project status, registry counts,
  validation tracks (R-parity / Stata-parity / reference-parity / Monte Carlo
  coverage / snapshot tests / citation audits), parity anchors, research-use
  statement (working-paper use; no published peer-reviewed article yet),
  open-core / commercial-downstream disclosure (StatsPAI Inc. + CoPaper.AI),
  and reproducible-check commands.
• Both pages added to the MkDocs navigation.

Changed — paper.md (JOSS manuscript)

• Repo URL casing corrected to canonical StatsPAI; added Zenodo archive
  reference (@wang2026statspai).
• Research-impact paragraph rephrased to match the actual current state:
  StatsPAI is used in working-paper workflows connected to Stanford REAP; no
  peer-reviewed research article using the package has yet been published.
• AI Usage Disclosure rewritten to spell out exactly what generative AI was
  used for (code generation, refactoring, test scaffolding, documentation
  drafting, manuscript copy-editing), to note that exact model identifiers
  were not retained for all exploratory sessions, and to confirm that
  generative AI will not be used to produce substantive responses to JOSS
  editors or reviewers.
• Acknowledgements split: explicit Author Contributions subsection
  attributing roles to each author, and an open-core / commercial-downstream
  disclosure (StatsPAI Inc. is the legal entity; CoPaper.AI is a commercial
  downstream product that may call the MIT-licensed StatsPAI package; the
  package itself remains permanently open source under MIT).
• paper.bib adds the wang2026statspai software entry pointing at the
  Zenodo concept DOI.

StatsPAI 1.15.5 — Agent-card coverage ratchet

StatsPAI 1.15.5

Released to PyPI and TestPyPI on 2026-05-21.

Added

• Agent-card coverage audit and CI ratchet: scripts/agent_card_coverage.py, docs/agent_cards_spec.md, and tests/test_agent_card_coverage.py.
• Generated baseline cards for the 1,018-function registry via scripts/gen_baseline_cards.py and src/statspai/_baseline_cards.py.
• FunctionSpec.inherits_from support so canonical estimator variants can inherit parent assumptions, preconditions, failure modes, alternatives, and typical_n_min while keeping method-specific metadata.

Changed

• Bumped package, docs, README, changelog, and citation metadata to 1.15.5.
• Refreshed registry statistics: 1,018 registered public functions across 80 submodules.
• Clarified DiD docs for continuous_did(method='cgs') and did_multiplegt_dyn as experimental MVP paths, not paper-parity estimators.

Validation

• pytest tests/test_agent_card_coverage.py tests/test_stability.py tests/test_mcp_protocol.py::TestResources::test_read_per_function_returns_agent_card -q --no-cov: 30 passed.
• python scripts/agent_card_coverage.py --check: passed.
• python scripts/registry_stats.py --check: passed.
• python scripts/schema_quality.py: passed.
• python -m build and twine check: passed.
• TestPyPI and PyPI install smoke tests for StatsPAI==1.15.5: passed.

PyPI: https://pypi.org/project/StatsPAI/1.15.5/
TestPyPI: https://test.pypi.org/project/StatsPAI/1.15.5/

StatsPAI 1.15.0 — Five polish waves (IV / synth / decomposition / ML+causal / RDD)

A bundled minor release on top of v1.13.1 covering five module-level polish waves. (1.14.0 was an internal cut never released to PyPI; PyPI publishes 1.13.1 → 1.15.0 directly.)

Headline

🔧 IV polish

• New sp.iv.iv_diag(data, y, endog, instruments, exog, ...) reporting bundle with Olea–Pflueger (2013) effective F, Lee–McCrary–Moreira–Porter (2022) tF adjusted CI, Anderson–Rubin / optional Moreira CLR / optional Kleibergen K weak-IV-robust sets, Kleibergen–Paap (2006) rk LM/Wald F, Conley–Hansen–Rossi (2012) plausibly-exogenous LTZ sensitivity, Blandhol–Bonney–Mogstad–Torgovitsky (2025) and Słoczyński (2024) TSLS-as-LATE caveat, and an OLS comparator (Young 2022).
• sp.iv.iv_compare(formula, data, methods=...) for k-class / JIVE side-by-side.
• 4 IV diagnostic plots (plot_iv_forest, plot_iv_forest_from_diag, plot_weak_iv_ci_overlay, plot_iv_diagnostics).
• sp.iv(absorb=...) HDFE 2SLS path.
• 13 verified paper.bib entries added; 4 enriched with volume/issue/pages.

🔧 Synthetic-control polish

• Every estimator now ships a publication-grade .to_latex() / .to_excel() / .to_word() pipeline.
• Trajectory + gap plots gain prediction-interval / pre-RMSPE ribbon options (Cattaneo-Feng-Titiunik 2021 JASA; Cattaneo-Feng-Palomba-Titiunik 2025 JSS).
• SDID schema canonicalised — sp.synth_report(method='sdid', ...) produces a full Markdown / text / LaTeX report instead of a row of N/As.
• 7 new 2022–2025 SCM citations Crossref/arXiv-verified.

🔧 Decomposition polish

• New sp.decomposition.yu_elwert distributional-decomposition module (Yu & Elwert 2024).
• Unified sp.decompose(...) dispatcher across the family (RIF / FFL / inequality / Oaxaca / Yu-Elwert).
• Shared influence-function / WLS / statistic-value backbone in decomposition/_common.py.
• Method-citations registry on every result class (.cite()).

🔧 ML+causal polish

• sp.dml_sensitivity + DMLSensitivityResult — Chernozhukov-Cinelli-Newey-Sharma-Syrgkanis (2022) "Long Story Short" framework: robustness value RV_q, RV_{q,α}, scenario bias bounds, benchmark-covariate comparisons, sensemakr-style contour plot().
• sp.dml_diagnostics — DoubleML-style 2×2 publication panel (overlap / score density / residual balance / orthogonality test) (Bach-Kurz-Chernozhukov-Spindler-Klaassen 2024 JSS).
• sp.cate_eval — Yadlowsky-Fleming-Shah-Brunskill-Wager (2025) RATE / AUTOC / Qini with closed-form influence-function SEs, decoupled from the forest backbone.
• ⚠️ Correctness fix — forest.CausalForest.best_linear_projection rewritten to AIPW pseudo-outcome Γ_i with HC1 SEs (Semenova-Chernozhukov 2021); the legacy plug-in OLS path was anti-conservative. Re-fit and report new HC1 numbers.
• ⚠️ Correctness fix — mediation.mediate no longer silently substitutes the point estimate for failed bootstrap replicates (which artificially shrunk SEs). Now retries up to 5x and surfaces boot_failure_rate in model_info. Regenerate SEs from prior versions if heavy bootstrap failure.
• Causal-discovery DAG visualisation (.to_networkx() / .to_dot() / .plot() / .edge_list()) on every result class plus sp.causal_discovery.{to_networkx, to_dot, plot_dag, edge_list, shd}.
• PolicyTreeResult promotion (Athey-Wager 2021): IF-SE on policy value, 95% CI from AIPW scores, Graphviz plot_tree().
• Mediation sensitivity plot upgrade (publication-style ACME(ρ) curve).
• 22 new ML-causal polish tests + DTR / QTE coverage closure.

🔧 RDD polish (the v1.15 cut name)

• 3 new estimators: sp.rd_flex (Noack-Olma-Rothe 2025 cross-fit ML adjustment), sp.rd_bias_aware_fuzzy (Noack-Rothe 2024 Econometrica AR-style weak-IV-robust fuzzy CI), sp.rd_discrete (Kolesár-Rothe 2018 AER honest CIs for discrete RVs).
• 3 reporting helpers: sp.rd_dashboard (4-panel CCT diagnostic), sp.rd_compare (side-by-side methods table), sp.rd_robustness_table (kernel × bandwidth × poly × donut sweep with to_latex() / to_excel()).
• sp.rdrobust polish: new rho parameter (Calonico-Cattaneo-Farrell 2018), discrete-RV warning when running variable has < 30 distinct values, weak-first-stage warning when fuzzy F < 10, first-stage F exposed at result.model_info['first_stage_F'].
• sp.rdplotdensity upgrade: boundary-adaptive local-polynomial CDF-regression density (Cattaneo-Jansson-Ma 2020 JASA).
• Dispatcher: sp.rd(..., method='flex' | 'bias_aware' | 'discrete') with full alias coverage.
• 21 new RD polish tests; all 156 RD tests pass.

Other

• fix(did): BJS imputation (Borusyak-Jaravel-Spiess 2024) DiD support repaired.
• New neural-causal / synth / spatial-DID export modules.
• tests/perf/05_feols_jax_bootstrap_bench.py — accelerator benchmark harness (CPU seq / CPU JAX vmap / GPU T4 / GPU A100).
• paper.bib cleanup: dropped duplicate semenova2023debiased entry; corrected sp.dml_panel attribution from misidentified Semenova-Chernozhukov 2023 to Clarke-Polselli 2025 Econometrics Journal 29(1).
• sp.synth(method='cluster') ClusterSC (Rho-Tang-Bergam-Cummings-Misra 2025, arXiv:2503.21629) author list corrected.

Stats

• 5284 tests passing on Python 3.13 / macOS arm64 (full suite, no skips of new code).
• 1017 functions in the agent registry.

Citations

All new citations DOI/arXiv-verified via Crossref + publisher + arXiv (2026-05-05 / 2026-05-06). See paper.bib and CHANGELOG.md for the full list.

Install

```bash
pip install StatsPAI==1.15.0
```

Full notes: CHANGELOG.md.

StatsPAI 1.13.1 — Stability tiers + external-validity dossier + cold-start surgery

Headline

Stability tiers, external-validity dossier, and cold-start surgery in a
single release. Every FunctionSpec now carries a stability field
plus per-function limitations, surfaced through
sp.describe_function, sp.help, sp.list_functions(stability=...),
the statspai list CLI, and the LLM-facing sp.function_schema
description; sp.recommend / sp.causal / sp.paper default to
dropping experimental / deprecated entries unless
allow_experimental=True is passed — closing a path where an agent
could silently land on a frontier MVP. Eight high-impact estimators
(aipw, aggte, pretrends_test, sensitivity_rr, mccrary_test,
oster_bounds, wild_cluster_bootstrap, rd_honest) are upgraded
from auto-registered stubs to hand-written specs with full assumption
/ failure-mode / alternative metadata. A weak-instrument preflight
gate in sp.preflight(data, "ivreg", formula=...) raises a structured
warning row when the first-stage F falls below the Staiger–Stock
(1997) or Stock–Yogo (2005) thresholds, and sp.recommend(... design='iv') adaptively reorders LIML / AR ahead of 2SLS on weak
first stages. A 36-module R parity harness, 21-module Stata parity
harness, 4-dataset original-paper replay (Card 1995, Callaway–Sant'Anna
mpdta, Abadie Basque, LaLonde NSW + PSID-1), Track-C performance
harness (HDFE / CS-DiD / SCM / DML log-log scaling), B=1000
Monte-Carlo coverage run, and a 900-trial CausalAgentBench prompt
suite all ship under tests/r_parity/, tests/stata_parity/,
tests/orig_parity/, tests/perf/, and tests/agent_bench/ with
paired R/Python/Stata drivers, JSON results, and 3-way Markdown +
LaTeX parity tables suitable for direct paper inclusion. A new
sp.validation_report() / sp.coverage_matrix() /
sp.reproduce_jss_tables() meta-API summarises the live registry,
materialises the parity / coverage / agent-bench artifacts as JSON,
and optionally re-runs the harnesses end-to-end so referees can
verify StatsPAI's external-validity claims without leaving Python.
Cold-start surgery in three steps brings sklearn submodules pulled by
import statspai from 245 to 0 (statspai.forest lazy-loaded —
Step 1B; 18 estimator files import sklearn lazily inside function
bodies — Step 1C; HAL TMLE classes drop sklearn class inheritance —
Step 1D), pinned by a new
test_sklearn_budget_ceiling_on_bare_import_statspai contract. The
workflow / paper orchestration layer replaces silent except: pass
paths with WorkflowDegradedWarning + structured degradations
records on the result object, so optional-stage failures surface in
PaperDraft.to_dict() and the rendered Pipeline notes section
instead of disappearing. sp.principal_strat(instrument=...) ships a
proper Angrist-Imbens-Rubin Wald-LATE estimator (the kwarg was
previously stubbed); sp.hal_tmle(variant='projection') keeps its
NotImplementedError but now points at a written-out RFC
(docs/rfc/hal_tmle_projection.md) instead of raising in silence.
Lazy-loading of optional families via __getattr__ keeps import statspai fast without breaking same-name function/subpackage
collisions (bartik, deepiv, proximal, …) — pinned by a
late-bind / post-import-shadow contract test and a committed
__init__.pyi stub generator so IDE / mypy see lazy-loaded names. A
latent Callaway–Sant'Anna REG inference scaling bug — discovered
because the parity harness flagged it — is fixed in
did/callaway_santanna.py.

Added

• Weak-instrument preflight gate in sp.preflight(data, "ivreg", formula=...). The new first_stage_strength check parses the
  Wilkinson IV formula, runs the first-stage OLS, and emits a
  warning row when the partial F-statistic falls below either the
  Staiger–Stock (1997) rule of thumb (F < 10, "very weak") or the
  Stock–Yogo (2005) 10% maximum-size critical value of 16.38 for
  one endog / one instrument (F < 16.38, "weak"). The warning
  payload includes structured recovery_hints pointing at
  method='liml', inference='ar', and sp.anderson_rubin_ci(...)
  so an LLM agent can branch on the typed envelope without parsing
  prose. Closes the §5.3 robustness DGP follow-up: Track B's
  test_iv_weak_instrument_undercoverage documents that 2SLS+HC1
  under-covers at the 0.88 level on a pi=0.10 first stage; the
  preflight now flags this before the user pays for the 2SLS fit.

• Adaptive IV ranking in sp.recommend(... design='iv'). When
  the live first-stage F is below 10 the recommendation list is
  reordered: LIML moves to the top, an Anderson–Rubin row is
  inserted, and the 2SLS row is annotated with
  very_weak_iv=True plus a rationale that explains the
  HC1-coverage failure mode. When F is in [10, 16.38) the order
  stays 2SLS → LIML but the rationales reference the
  Stock–Yogo threshold. The 2SLS row now also carries
  first_stage_F as a numeric field so downstream tooling can
  consume it without regex.

• Preflight / recommend tests. tests/test_preflight.py adds
  TestIVFirstStageStrength (6 cases covering strong / weak /
  borderline / non-IV / missing-columns / JSON-safe payloads).
  tests/test_smart_workflow.py::TestRecommend adds two cases
  pinning the 2SLS-first ordering on strong instruments and the
  LIML-first / AR-included ordering on weak instruments.

• R parity harness — 36 paired R/Python modules. New
  tests/r_parity/ ships 36 paired scripts (one R, one Python) that
  replay the same DGPs through fixest / did / csdid / gsynth /
  MatchIt / DoubleML / rdrobust / Synth / lme4 / plm /
  frontier / MR-PRESSO / lavaan / mediation / WeightIt /
  cobalt / mlogit / nlme / ordinal / … and StatsPAI's matching
  estimator. Each module emits <id>_R.json and <id>_py.json;
  compare.py produces a 3-way parity table (parity_table.md /
  .tex / parity_table_3way.md / .tex) tightened with a small
  ID-column + longtable layout for direct paper inclusion. Two
  parallel tracks — "orig" (canonical-dataset replays) and "perf"
  (timing under matched DGPs) — share the same compare-tooling.

• Stata parity harness — 21-module StatsPAI ↔ Stata 3-way
  compare. New tests/stata_parity/ ships 21 paired
  .do/.py scripts covering reghdfe, xtreg, csdid, did_imputation,
  synth, synth_runner, ivreg2, xtivreg, rdrobust,
  psmatch2, teffects, xtfrontier, mixed / melogit /
  mepoisson, bayes, gmm, boottest, plus the Stata→Python
  translator round-trip. Drivers write Stata results to JSON via the
  stata-mcp stata_do tool; Python drivers run StatsPAI through
  import statspai as sp; compare_stata.py joins on (module, estimator, statistic) and emits parity_table_stata.md /
  parity_table_stata.tex plus the 3-way StatsPAI ↔ R ↔ Stata table.

• Canonical-dataset original-paper replays. New
  tests/orig_parity/ adds 4 module pairs that replay each paper's
  headline number bit-equal to the published value: Card (1995)
  returns-to-schooling on wooldridge::card; Callaway–Sant'Anna
  (2021) staggered DiD on did::mpdta (the package's vendored
  minimum-wage panel); Abadie–Diamond–Hainmueller Basque on
  Synth::basque; LaLonde (1986) NSW on MatchIt::lalonde plus a
  4b sub-module on causalsens::lalonde.psid (true Dehejia–Wahba
  NSW + PSID-1, 2675 obs) where sp.regress + sp.psm recover the
  published −15,205 to relative tolerance 1.5e-05. Drivers + bundled
  data + JSON results + parity_table_orig.md are all committed.

• Track-C performance harness — log-log timing. New
  tests/perf/ adds matched R/Python timing runs for HDFE
  (fixest::feols vs sp.feols), CS-DiD (did::att_gt vs
  sp.callaway_santanna), SCM (Synth::synth vs sp.synth), and
  DML (DoubleML::DoubleMLPLR vs sp.dml) at log-spaced N. Drivers
  plus compare_perf.py produce perf_table.md / .tex and a
  track_c_loglog.{pdf,png} log-log scaling figure.

• Coverage Monte Carlo at B=1000. New
  tests/coverage_monte_carlo/run_b1000.py measures 95% CI coverage
  for OLS (0.952), 2×2 DiD (0.955), and strong-Z IV (0.962) — all
  inside the 99% Wilson band [0.935, 0.967] around nominal 0.95. The
  full slow pytest tests/coverage_monte_carlo/ -m slow sweep at
  B=1000 also passes 8/8 (753.51 s). Frozen run lives at
  results_b1000/coverage_b1000.json; previous fast-track results
  remain at the default B=200.

• CausalAgentBench scaffolding (mock-mode shipped, API run gated).
  New tests/agent_bench/ ships a 50-prompt × L1/L2/L3 difficulty
  × 6 cells × 3 reps = 900-trial agent bench with a deterministic
  mock-LLM runner (runners/mock_llm.py), a frozen OSF
  pre-registration protocol (prompts/_protocol.md), and a grader
  (runners/grader.py) emitting an H1–H5 directional results table.
  Mock dry-run completes in <1 s and produces
  results/headline.md + results/scores.csv + results/trials.jsonl;
  the production --api flag is one switch away once the OSF
  pre-registration and API budget clear.

• sp.validation_report / sp.coverage_matrix /
  sp.reproduce_jss_tables — JSS-grade validation meta-API. New
  src/statspai/validation.py (863 LOC, 58-test battery in
  tests/test_jss_validation_api.py) exposes three top-level
  functions for the paper-submission audit trail:
  sp.validation_report() summarises the live sp.registry plus
  the materialised parity / coverage / agent-bench artifacts as a
  structured ValidationReport (registry.total_functions,
  evidence.r_parity_modules, evidence.stata_parity_modules,
  evidence.coverage_b1000, evidence.agent_bench_trials, …) with a
  one-paragraph .summary() and full JSON .to_dict();
  sp.coverage_matrix() enumerates every reference-implementation
  parity claim with its expected tolerance, observed gap, and the
  driver script that produced the JSON; sp.reproduce_jss_tables()
  returns a ReproductionResult enumerating the exact `Rscr...

StatsPAI 1.12.1 — sp.citation() + Zenodo DOI

Citation metadata polish — no numerical or API changes to any estimator. Safe drop-in for 1.12.0.

⚠️ DOI correction (post-publish): the original notes for this release accidentally cited DOI 10.5281/zenodo.18636688, which actually belongs to a different project (constitutional-alignment). The correct StatsPAI Zenodo concept DOI is 10.5281/zenodo.19933900 (versioned DOI for v1.12.1: 10.5281/zenodo.19934144). The PyPI wheel for 1.12.1 still embeds the wrong DOI in sp.citation(); a 1.12.2 patch with the correct DOI will follow shortly.

Highlights

• sp.citation() — package-level citation helper:

  import statspai as sp
  sp.citation()            # BibTeX (default)
  sp.citation("apa")       # APA-style human-readable
  sp.citation("plain")     # minimal plain text
  sp.__citation__          # default BibTeX as a str

  Distinct from sp.cite(), which formats individual coefficients inline.
• CITATION.cff at the repo root — GitHub now renders a "Cite this repository" button. Bundled in the sdist via MANIFEST.in.
• Zenodo concept DOI 10.5281/zenodo.19933900 — always resolves to the latest archived release. Surfaces in sp.citation() output (1.12.2+), the README citation block, and a DOI badge next to the JOSS-pending status badge.
• .zenodo.json — future GitHub Releases mint version-specific DOIs with consistent metadata.

A JOSS paper for StatsPAI is currently under review; once accepted, the journal article will become the preferred citation and sp.citation() will be updated to return it.

Citing this release

@software{wang_statspai_2026,
  author       = {Wang, Biaoyue},
  title        = {StatsPAI: The Agent-Native Causal Inference \& Econometrics Toolkit for Python},
  year         = {2026},
  version      = {1.12.1},
  doi          = {10.5281/zenodo.19934144},
  url          = {https://doi.org/10.5281/zenodo.19934144},
  license      = {MIT},
}

For the always-latest-version DOI, use the concept DOI 10.5281/zenodo.19933900 instead of the v1.12.1 versioned DOI shown above.

Install

pip install --upgrade statspai

⚠️ Prefer pip install statspai==1.12.2 once it's published — 1.12.1 ships an incorrect Zenodo DOI in sp.citation() output.

StatsPAI 1.12.0 — DML hardening + TMLE correctness pass

Headline

The whole dml/ module got a careful audit. sp.dml / sp.dml_panel
/ sp.dml_model_averaging all stay backwards-compatible at the
call-site level (existing scripts keep working) but several internal
numerical behaviours change — see the ⚠️ Correctness section and
MIGRATION.md.

⚠️ Correctness

• sp.dml(model='irm') and sp.dml(model='iivm') now use
  StratifiedKFold (stratified by D and Z respectively) — the old
  KFold could produce a fold whose subgroup mask was empty, in which
  case the AIPW score for that fold's test rows was silently filled
  with zeros (biased point estimate, biased SE). Empty subgroups now
  raise IdentificationFailure with a clear remedy. Estimates may
  shift slightly on data sets where the old KFold happened to
  produce extreme folds.
• sp.dml_panel(binary_treatment=True) is now a deprecated no-op. The
  previous classifier path fit a propensity on within-demeaned features
  but raw {0,1} labels — there is no clean interpretation as
  E[D̃ | X̃] for the result. The estimator now always uses a regressor
  on D̃ (PLR-with-FE is agnostic to D's type). A DeprecationWarning
  is emitted, and D ∈ {0,1} is validated when the flag is True.
• sp.dml_model_averaging now drops rows with NaN in y / treat /
  covariates / sample_weight (matching every other DML class);
  previously NaNs propagated into sklearn fits and could produce NaN
  estimates undetected by the existing denom < 1e-12 guard.
• sp.dml_model_averaging: the default weight_rule is now
  "short_stacking" — Ahrens, Hansen, Schaffer & Wiemann (2025, JAE)
  eq. 7 — which solves a constrained least squares stacking problem on
  cross-fitted nuisance predictions and plugs the stacked nuisance
  into a single PLR moment equation. The previous "inverse_risk"
  default (heuristic 1/MSE-weighted average of per-candidate θ̂_k) was
  not in the cited paper and is preserved as a clearly labelled
  baseline. New "single_best" matches the paper's footnote 8
  formulation. Per-nuisance stacking weights are exposed as
  model_info["weights_g"] / weights_m.
• sp.dml(model='pliv') raises RuntimeError when the
  ML-residualised partial correlation |corr(z̃, d̃)| falls below
  1e-3 (was 1e-6, too lenient to catch genuine weak-IV collapse).
  A new model_info["diagnostics"] block reports the partial
  correlation and an approximate first-stage F.

Added

• All four sp.dml(model=…) variants now accept a random_state=
  argument (default 42) controlling fold assignment. Repeated splits
  use random_state + rep so a single seed fully determines the
  result.
• sample_weight= support on sp.dml(model='plr'), sp.dml(model='irm'),
  sp.dml_panel, and sp.dml_model_averaging (any weight rule). The
  weighted estimator uses a Z-estimator sandwich variance throughout.
  sp.dml(model='pliv') and sp.dml(model='iivm') raise
  NotImplementedError if a non-trivial weight is supplied — the
  weighted Wald-ratio variance derivation is non-trivial and lands
  in a follow-up. sample_weight may be passed as a 1-D array, a
  pandas Series, or a column name string.
• New model_info["diagnostics"] block on every variant:
  • PLR: residual scales, partial correlation y_resid·d_resid,
    within-R² of each nuisance.
  • IRM: propensity p01/p99/min/max, n clipped below/above the
    [0.01, 0.99] overlap clip, n times the subgroup g̃₁/g̃₀ fit
    fell back to the subgroup mean.
  • IIVM: instrument-propensity p01/p99/min/max, clipping counts,
    subgroup fallbacks for both g(z, X) and r(z, X), and
    E[ψ_b] (the LATE Wald-ratio denominator — proximity to zero
    indicates a weak first stage).
  • PLIV: first-stage partial correlation, approximate first-stage
    F, residual scales.
  • panel_dml: y/d residual std, within-R², cluster Ω, weighted flag.
• sp.dml_panel(sample_weight=…) does a weighted within transform
  (subtract weighted unit / time means) and reports a weighted
  Liang-Zeger cluster SE.

Changed

• Internal flag rename _BINARY_TREATMENT → _ML_M_TARGET_BINARY and
  _BINARY_INSTRUMENT → _ML_R_TARGET_BINARY on the per-model DML
  classes. The new names describe the nuisance-target shape
  (the IIVM ml_m actually models the instrument propensity, not D).
  These flags are private (underscore-prefixed); no public API change.
• paper.bib: filled in the missing volume / number / pages
  fields on @ahrens2025model (40(3):249–269), verified via the Wiley
  Online Library record and the JAE issue listing.

Internal

• Per-rep diagnostics now flow back to model_info["diagnostics"]
  via a new _aggregate_diagnostics helper on _DoubleMLBase. Each
  subclass populates self._last_rep_diagnostics inside
  _fit_one_rep; the base merges across reps (sum for counts, mean
  for floats, OR for booleans, concat for lists).

⚠️ Correctness — TMLE module audit pass

• sp.tmle.SuperLearner previously ran NNLS and post-hoc-normalised
  weights to sum to 1, which is not the simplex-constrained
  optimum (rescaling an unconstrained NNLS solution gives the simplex
  optimum only when the unconstrained sum already equals 1, a
  measure-zero event). Replaced with a direct SLSQP QP on the
  simplex; ensemble predictions are now genuinely the convex
  combination minimising squared loss. Affects every downstream
  caller — sp.tmle, sp.hal_tmle, and any user code that builds a
  Super Learner directly. Numerical results will shift slightly on
  data sets where the old NNLS solution did not happen to be on the
  simplex.
• sp.tmle.ltmle censoring half-implementation: the regime-following
  indicator now includes & (C_k_obs == 1) so censored units are
  excluded from the targeting equation rather than continuing to
  contribute with 1/p_c-inflated weights. (sp.tmle.ltmle_survival
  was already correct on this; ltmle.py was the regression.)
• sp.tmle.ltmle_survival influence function: previously used
  -H * (T_k - h_star_regime) summed across intervals as the
  influence function for both the RMST contrast and the
  terminal risk difference at K. The proper EIF for :math:E[S^a(t)]
  (Cai & van der Laan 2020) needs the survival-product factor
  :math:S^a(t)/S^a(j) and the IC for the terminal RD at K is the
  EIF of :math:S^a(K) alone (NOT the cumulative-across-K RMST IC).
  Refactored _run_regime to expose the per-subject sequences
  S_seq, h_star_seq, H_seq, T_seq; the call site now
  computes the RMST and terminal-RD EIFs separately via
  _eif_rmst and _eif_survival_at_k. SE estimates change —
  generally smaller for RMST (was conservative), and the terminal-RD
  SE is now correctly tied to its target functional rather than
  picking up RMST's cross-time aggregation.
• sp.hal_tmle(variant='projection') was a no-op in v1.11.x and
  earlier. The projection variant ran an ad-hoc shrinkage on
  model_info["eps"] after the point estimate had already been
  computed; the variant flag did not change the estimate. The path
  now raises :class:NotImplementedError honestly until the proper
  Riesz-projection step (Li-Qiu-Wang-vdL 2025 §3.2) is ported.
• sp.hal_tmle docstring previously claimed the basis was "rich
  enough to approximate any càdlàg function of bounded variation",
  the property of full HAL (Benkeser & van der Laan 2016). The
  implementation only builds main-effects indicator basis
  functions :math:\\mathbb 1\\{x_j \\le a_j\\} — i.e.
  L1-penalised additive piecewise-constant regression, NOT full HAL.
  Docstring is corrected; numerical behaviour unchanged.

Fixed — TMLE convergence + overlap diagnostics

• sp.tmle._fit_epsilon now emits a UserWarning when the Newton
  iteration on the fluctuation parameter fails to converge in
  max_iter steps, instead of silently returning the last value
  (which yields a non-targeted plug-in). The warning includes the
  final score magnitude and ε for diagnosis.
• sp.tmle now reports model_info['propensity_diagnostics'] (min,
  max, p01, p99, n clipped below/above, clip share) and emits a
  UserWarning when ≥ 5 % of propensities hit the
  propensity_bounds clip — same overlap convention as
  sp.metalearner. AIPW scores blow up at e≈0/1, so heavy clipping
  silently changes the estimand from ATE in the population to ATE
  on the trimmed sample.
• sp.tmle.SuperLearner(task='classification') validates that the
  target is binary (was silently dropping non-{0,1} columns of
  predict_proba); switches to StratifiedKFold so every fold has
  both classes; predict() clips to (1e-6, 1-1e-6) for
  classification (was inconsistent with predict_proba which
  already clipped).

Fixed — TMLE / HAL-TMLE citations (§10 verification pass)

• paper.bib now records three previously-uncatalogued HAL-TMLE
  references with full Crossref/arXiv-verified metadata (added
  2026-04-30):
  • @li2025regularized — arXiv:2506.17214, verified via
    arxiv.org. Earlier inline-cited title in hal_tmle.py was
    "Highly Adaptive Lasso Implementations"; the paper's actual
    title is "Highly Adaptive Lasso Implied Working Models" —
    fixed in docstring + model_info['citation'].
  • @vanderlaan2023efficient — IJB 19(1):261–289,
    doi 10.1515/ijb-2019-0092, verified via degruyterbrill.com.
  • @benkeser2016highly — IEEE DSAA 2016, pp. 689–696,
    doi 10.1109/DSAA.2016.93, verified via Crossref API.
• tmle.py:_CITATIONS['tmle'] now includes the vanderlaan2006targeted
  reference that the docstring already cites (was missing — docstring
  promised it via [@vanderlaan2006targeted] but the inline BibTeX
  registered only vanderlaan2007super). Author punctuation /
  capitalisation aligned to paper.bib.
• ltmle_survival.py cai2020step reference reformatted to match
  paper.bib (year 2020 vs the previous docstring's 2019; the IJB
  volume's nominal year is 2020)...

StatsPAI 1.9.0 — Agent-native API surface (12 modules)

The 1.9.0 line ships StatsPAI's first deliberately agent-shaped API
surface — 12 new top-level entry points designed for Claude Code /
Cursor / Copilot CLI workflows where the LLM, not a human, is doing
the calling. No estimator numerical paths changed; all
additions are new functions or strictly additive parameters with
"agent" as the default so existing behaviour is byte-identical.

Added — Agent serialization & error envelope (Phase 1)

• CausalResult.to_dict(detail=...) and
  EconometricResults.to_dict(detail=...) — unified payload
  control with three documented levels:

  • "minimal" (~150 tokens) — bare answer; no diagnostics.
  • "standard" (~250 tokens) — current default; coefficients +
    scalar diagnostics + detail_head rows. Byte-identical to
    legacy to_dict().
  • "agent" (~620 tokens) — adds violations / warnings
    / next_steps / suggested_functions so an LLM can plan
    its next call without another round-trip.

  for_agent() is now a thin alias for to_dict(detail="agent");
  to_agent_summary() is unchanged but its docstring now points
  at to_dict(detail="agent") as the canonical flat form.

• execute_tool MCP error envelope — when an estimator raises
  a structured StatsPAIError subclass, the MCP tools/call
  response now surfaces error_kind (e.g.
  "method_incompatibility") plus the full error_payload
  dict (code / recovery_hint / diagnostics /
  alternative_functions). Legacy error / remediation
  fields preserved.

Added — MCP server polish (Phase 1)

• statspai-mcp console script wired in pyproject.toml so
  pip install statspai exposes it on PATH.
• statspai://function/{name} per-function resources surfacing
  the registry's full agent-card (description, signature,
  assumptions, failure_modes, alternatives, typical_n_min, example).
  Listed via the new resources/templates/list handler.
• statspai://functions machine-readable JSON index for
  one-shot tool discovery.
• Typed JSON-RPC errors mapped to canonical MCP codes:
  -32002 (resource not found), -32602 (invalid params),
  -32000 (server fallback). Replaces the previous blanket
  -32000.
• notifications/* silenced — Claude Desktop / Cursor send
  notifications/initialized after the handshake; the server now
  drops any method whose name starts with notifications/ per
  the MCP spec, instead of replying with -32601 noise on every
  session.
• MCP-level detail parameter on tools/call — agents pick
  detail="minimal" | "standard" | "agent" per call to control
  token cost. Validation rejects invalid values with -32602.

Added — Workflow primitives (Phases 2-4)

• sp.audit(result) — missing-evidence checklist (the
  read-only counterpart to sp.assumption_audit): inspects what
  robustness / sensitivity diagnostics are stored on a fitted
  result and surfaces which method-family checks are still
  missing. Returns {checks: [{name, question, status, severity, importance, suggest_function, ...}], summary, coverage} with
  18 curated checks across DID/RD/IV/synth/matching/OLS.

• sp.detect_design(data, **hints) — heuristic design
  identifier: returns {design, confidence, identified, candidates, n_obs, columns} with design ∈ {"panel", "rd", "cross_section"}. Symmetric (unit, time) pair dedup; RD
  confidence capped at 0.30 without explicit hint to avoid
  noise-data false positives.

• sp.preflight(data, method, **kwargs) — method-specific
  pre-estimation diagnostics distinct from
  sp.check_identification (design-level) and
  sp.assumption_audit (re-runs tests). Cheap shape / column /
  treatment-binarity / sample-size checks per method family;
  returns {verdict: "PASS" | "WARN" | "FAIL", checks, summary, known_method}.

• CausalResult.cite(format=...) and
  sp.bib_for(result) — multi-format citations:
  "bibtex" (default, byte-identical to legacy cite()),
  "apa" (parsed prose), "json" (structured {type, key, authors, year, title, journal, volume, number, pages, publisher, fields}). LaTeX-diacritic normalisation
  ({\\"o} → ö); multi-entry BibTeX strings (e.g.
  twfe_decomposition cites both Goodman-Bacon 2021 AND de
  Chaisemartin & D'Haultfœuille 2020) round-trip both authors —
  zero hallucination per CLAUDE.md §10.

• sp.examples(name) — runnable code snippets for any
  registered function; 10 hand-curated flagship snippets, falls
  back to registry.example for the rest.

• sp.session(seed=42) — deterministic-RNG context manager
  snapshotting Python random and NumPy's legacy global MT19937
  generator; restores prior state on exit even when an exception
  is raised inside the block. Lazy torch / jax interop — never
  auto-imports. Documented escape hatch for
  np.random.default_rng() (which is not covered — pass
  state.seed explicitly).

• result.brief() / sp.brief(result) — one-line
  dashboard string (~95 chars typical, ≤ 140 hard cap) for
  multi-result agent loops.

• MCP prompts/list + prompts/get — three curated
  workflow prompt templates (audit_did_result /
  design_then_estimate / robustness_followup) surfaced as
  one-click buttons in MCP-compliant clients.

Changed

• CausalResult.to_dict / EconometricResults.to_dict now
  accept a keyword-only detail parameter. Default "standard"
  preserves the legacy shape exactly. CausalResult's
  detail_head is also keyword-only now (was positional-or-
  keyword) to close the to_dict("agent") foot-gun.

• CausalResult.cite() now accepts format= keyword; zero-arg
  call still returns BibTeX, byte-identical to
  cite(format="bibtex").

Tests

+422 targeted tests across the agent stack, all passing.
Token-budget assertions pin the size of every detail level so
future changes can't accidentally bloat the LLM tool-result channel.

No numerical changes

Every estimator's coefficient / SE / CI / p-value path is byte-
identical to 1.8.0. The 12 new modules are introspection,
serialization, prompt-rendering, and RNG-management primitives —
they read from existing result state, never recompute it.