Statspai

bring it to 215 (27.8%).
- **The index now refuses to build such an entry.** A promotion claiming an
R or Stata side must list a test that loads a reference fixture or calls
R, or record a `provenance` for embedded constants (the live-captured
Stata values behind `did_2x2` and `ddd` are the one legitimate case of
that today). `sides` is also now required: it defaulted to `["py", "R"]`,
so an entry could claim an R comparison by omitting the field. Guarded by
`tests/test_parity_index.py`.

⚠️ Correctness fixes — network

- **`sp.dyadic_regression` double-counted reciprocal dyads in its robust
variance.** The Aronow–Samii–Assenova estimator lets any two dyads that
share a member covary, each pair entering once. The implementation summed
scores by node and took `Σₙ SₙSₙ′ − g′g`, which weights a pair by the
*number* of members it shares. That coincides with the definition when
every unordered pair appears once — so undirected results are unchanged —
but on directed data `(i, j)` and `(j, i)` share both members and their
cross terms entered twice: 1.8% on the standard errors of a 20-node
design. The fix subtracts the within-pair score sums instead of the
per-row ones. **Directed dyadic regressions should be re-run.** Rows with
`i == j` now raise rather than entering with an inconsistent weight. The
repository's own brute-force test already encoded the correct definition;
it only ever ran on undirected data.

Added — network

- Cross-language evidence for 19 network functions against R `igraph`,
`sna`, `ergm` and `dyadRobust` (`tests/reference_parity/test_network_parity.py`):
centrality scores, PageRank, Bonacich power, Katz, Wasserman–Faust
closeness on a disconnected graph, transitivity, assortativity,
reciprocity, modularity, components, `network_summary`, QAP `netlm` /
`netlogit` coefficients, ERGM MPLE, dyadic-robust OLS, and the bundled
karate and Florentine datasets — all at 1e-9 or better, except ERGM MPLE
standard errors at ≤3.2e-7. HITS and the eigenvector column of
`sp.centrality` agree with igraph up to a documented normalisation (the
ratio is constant across nodes); Louvain is compared as a 200-seed
distribution (T3).
- `sp.eigenvector_centrality`'s docstring claimed its L2 normalisation
matched igraph; igraph max-scales. Corrected. `sp.florentine_families`
now documents that the Pucci isolate is omitted relative to `ergm`'s
16-node `flomarriage`.
- Found while building the reference: the R package `dyadRobust` recodes
node ids inside one `dplyr::mutate()` call, and because `mutate`
evaluates sequentially the alter recode is built from the already-recoded
ego column. A node then gets different codes as ego and as alter unless
the ids are already `1..N` in first-appearance order. The fixture uses such
ids and asserts the precondition.

⚠️ Correctness fixes — Mendelian randomisation

Found by comparing every MR function against the R packages that define
it — `MendelianRandomization`, `TwoSampleMR`, `RadialMR`, `MRPRESSO` and
`mr.raps` — on the 28-variant LDL-C → CHD data shipped with
`MendelianRandomization`. The percentages are on that example. **MR
analyses should be re-run**; `MIGRATION.mdmr-sweep` has the table.

- **`sp.mr_ivw` reported the fixed-effect standard error regardless of
heterogeneity** — 0.276 against the reference's 0.530 here (Cochran's
Q p = 3e-10). Both reference packages default to multiplicative random
effects (SE × `max(1, RSE)`) with more than three variants, and so does
`sp.mr_ivw` now; `model="fixed"` restores the old SE. The point estimate
is unchanged. `sp.mr("ivw")`, `sp.mendelian_randomization` and
`sp.mr_leave_one_out` inherit the change (the last used to recompute a
fixed-effect SE inline, 49% off).
- **`sp.mr_egger` did not orient variants**, so its result depended on
which allele each variant was coded against: slope 14% and intercept —
the pleiotropy test — 38% off, with 16 of 28 variants having negative
exposure associations. Its residual standard error was also not floored
at 1. `sp.mr_pleiotropy_egger` (which carried its own copy of the fit)
and Rücker's Q′ in `sp.mr_heterogeneity(method="egger")` are fixed with
it.
- **`sp.mr_median`**: a step weighted median instead of Bowden et al.'s
interpolated one (0.8%); a penalty using the *lower*-tail χ² probability
— up-weighting exactly the heterogeneous variants it exists to suppress —
with no `min(1, 20 q)` cap and centred on IVW instead of the weighted
median (38%); and a bootstrap that redrew the weights. New
`weighting="simple"`.
- **`sp.mr_mode`**: bandwidth, grid and bootstrap all differed from
Hartwig et al.'s (estimate 7% off). The SE is now the bootstrap MAD, as
in the reference. New `phi=` and `refine=`.
- **`sp.mr_cml` penalised its BIC by the number of variants instead of the
GWAS sample size** the method requires (log 28 against log 17723), and
selected six invalid variants where the reference selects two; its SE
held the nuisance exposure effects fixed (~12% small). New `n=` (without
it the old stand-in is used and a warning says so) and
`model_average=`; `K_max` defaults to `n_snps − 2`.
- **`sp.mr_raps` was a different estimator under the name**: a Tukey loss
minimised jointly over `(β, log τ²)` with a one-dimensional sandwich
treating τ² as known — τ² 3.5× too small, the estimate 1.4% off, the SE
47% too small. Rewritten as a port of `mr.raps` 0.4.3 (`loss=`,
`over_dispersion=`, `pruning=`, `niter=`, `tol=`); **the default loss is
now Huber**, the package default. `tuning_c` passed without `loss` warns
that it now applies to Huber; `beta_init` / `tau2_init` are ignored with
a `DeprecationWarning`.
- **`sp.mr_steiger`'s p-value was one-sided and computed as `1 − Φ(z)`**,
exactly 0.0 beyond z ≈ 8.3 (the reference reports 1.8e-73 here). It is
two-sided by default, as `TwoSampleMR` reports it; `alternative="greater"`
gives the one-sided test.
- **`sp.mr_presso` is now a port of `MRPRESSO::mr_presso`.** Outlier
p-values were not Bonferroni-adjusted (seven outliers flagged where the
reference flags two, moving the corrected estimate 2.4%), SEs were
fixed-effect rather than residual-scaled (half the reference's), the
outlier test ran when the global test did not reject, and the
"distortion test" was a z-test rather than PRESSO's resampling test.
**Monte Carlo p-values return to the reference's `k / B`**, reversing
1.5.0's `(k + 1) / (B + 1)`: with Bonferroni the latter's smallest
attainable outlier p-value is `n / (B + 1)`, above 0.05 for 50 variants
at B = 1000, so the test could never flag anything. A 0 now means
"below the Monte Carlo resolution", and a warning fires when `n_boot`
cannot resolve the adjusted threshold. New fields `outlier_pvalues`,
`distortion_coefficient`.
- Every p-value in `sp.mendelian` now comes from a survival function:
`1 − cdf` is accurate only to ~1e-16 absolute, 1e-7 relative on the
heterogeneity p-value of 3e-10 here.

Added — Mendelian randomisation

- Cross-language evidence for 15 MR functions
(`tests/reference_parity/test_mr_R_parity.py`): `mr`, `mr_ivw`,
`mr_egger`, `mr_median`, `mr_mode`, `mr_cml`, `mr_steiger`,
`mr_radial`, `mr_presso`, `mr_leave_one_out`, `mr_pleiotropy_egger`,
`mr_heterogeneity`, `mr_f_statistic` and `mr_multivariable` at 1e-9 or
better on every deterministic quantity (bootstrap SEs and simulated
p-values are Monte Carlo on both sides and not compared). `mr_raps` is
graded `aligned`: the Gaussian fits match at 1e-8, and for the robust
losses the sandwich reproduces `mr.raps`'s SEs at 1e-10 when evaluated at
its own estimates, while the fitted values agree to 5e-5 (β) and 1.5e-3
(τ²) because `mr.raps` stops `uniroot` at its default tolerance and
`integrate` at 1.2e-4 — StatsPAI's root satisfies the estimating
equation more tightly.
- `sp.mr_radial(bonferroni=, alpha=)`: `bonferroni=False` reproduces
`RadialMR::ivw_radial`'s outlier list; the Bonferroni default is kept.
- The Mendelian-family guide is updated for all of the above.

⚠️ Correctness fixes — decomposition

Found by comparing against `DasGuptR`, `ddecompose` and `cdgd` (the last
two maintained by the methods' authors or their groups).

- **`sp.das_gupta` decomposed a different aggregate on stratified data.**
Its docstring defines the aggregate as `Σᵢ ∏_f f_{f,i}` over strata, but
with more than one row per population the code multiplied the factor
*means*. On Das Gupta's own Table 6.5 (four factors × six age groups)
factor A's share of the 1963→1968 change came out 0% where `DasGuptR`
gives 36.8%, and factor D's +333% where it gives −52.5%. It now sums the
per-stratum decompositions; rows are paired by the new `by=` or by
position, and unequal row counts raise. Single-row inputs are unchanged.
- **`sp.gap_closing(method="ipw")` reweighted in the wrong direction**, and
so did the reweighting term of `method="aipw"`: both used the reciprocal
of the density ratio `f_target / f_source`. With `Y = 2X` and no group
effect, the IPW counterfactual gap came out at twice the observed gap
instead of zero; AIPW stayed right only while its outcome model was
right, i.e. it was not doubly robust. IPW now matches
`ddecompose::dfl_decompose` in both directions (1e-9);
`method="regression"` was already right and matches
`ddecompose::ob_decompose`.
- **`sp.yu_elwert_decompose(method="efficient")`'s components did not add
up to the disparity** (0.566 against 0.572 on a 2,000-row example): it
computed selection as a covariance of doubly robust scores rather than
as the residual, and used unnormalised inverse-probability weights. It
now follows `cdgd::cdgd0_manual` — Hajek-normalised potential outcomes,
selection as the residual — and matches it to 1e-9 given the same
nuisance fits. New `inference="analytic"` gives the package's
efficient-influence-function standard errors. `method="plugin"` (the
default) is unchanged.

- **The Gini RIF averaged to neither Gini.** `sp.rifreg(statistic="gini")`,
`sp.rif_values` and the Gini path of `sp.ffl_decompose` combined a
midpoint ECDF with the `n/(n−1)`-corrected Gini; on `cps_wage` the RIF
mean was 0.10405 against a plug-in Gini of 0.10442 and RIF-regression
coefficients were up to 1% off. The RIF is now `dineq::rif(method =
"gini")`'s, to 1e-12 per observation, whose mean is exactly the plug-in
Gini; the Gini inside the decomposition machinery (`dfl_decompose`,
`ffl_decompose`) is that plug-in Gini. `sp.inequality_index("gini")`
keeps its documented bias-corrected default.
- **`sp.ffl_decompose` mislabelled and mis-signed its error terms.**
`spec_error` held the reweighting error and `reweight_error` the
specification error; with `reference=1` the specification error also had
the wrong sign, so the four components did not add up to the gap (0.1184
against 0.1216 for the variance on `cps_wage`). `gap` is now the
difference of RIF means — the plug-in statistics the components add up
to, exactly; for the variance that removes the `n/(n−1)` factor it
carried. Every term now matches `ddecompose::ob_decompose(reweighting =
TRUE)` in both directions for the variance, Gini and 10th/50th/90th
percentiles.

- **`sp.gelbach`'s standard errors dropped a covariance.** They were a
two-term delta method treating the auxiliary regressions and the long
regression as independent. They now follow Gelbach's own Stata command
`b1x2` (each added variable its own group): the stacked-system sandwich,
plus the part due to the long-regression coefficients and both cross
terms, robust by default and homoskedastic with `robust=False`. The
result gains the joint `vcov` of the contributions and `total_se`
(`1′V1`, as `b1x2` reports it). Point estimates are unchanged; on
`cps_wage` the SEs moved by 0.03%.

- **`sp.bauer_sinning` / `sp.yun_nonlinear` reported no inference and no
detailed unexplained part.** They now port `mvdcmp` (Powers, Yoshioka &
Yun): the Yun-weighted detail of both the explained and the unexplained
component (`detailed_unexplained`, including the constant), the joint
delta-method covariance `vcov`, and `se` for the explained, unexplained
and total gap. Point estimates are unchanged. Probit fits now finish with
observed-information Newton steps: Fisher scoring stopped ~1e-8 short of
the MLE (this also touches `sp.fairlie(model="probit")` at that level).

Added — decomposition

- Cross-language evidence for `das_gupta`, `kitagawa_decompose`,
`gap_closing`, `yu_elwert_decompose`, `rifreg`, `ffl_decompose` and
(beyond the mean) `dfl_decompose`
(`tests/reference_parity/test_decomp_R_parity.py`), on Das Gupta's worked
examples, simulated data and `cps_wage`, read by both sides from the same
bytes.
- Stata-side evidence for `gelbach` (`b1x2`, robust and homoskedastic, to
1e-9 on the full covariance), `subgroup_decompose` (`ineqdeco`, GE(0),
GE(1), GE(2) total / within / between to 1e-12) and `source_decompose`
(`descogini`, every share, Gini and Gini correlation to 1e-12).
- Evidence for `bauer_sinning` against Stata `mvdcmp`: logit to 1e-9 on
every term, detail and covariance; probit aligned to 2e-6, because
`mvdcmp` fits its probits at Stata's default tolerance (the fixture
records that Stata's tightly converged probit equals StatsPAI's to
1e-12).
- `sp.source_decompose(gini="population")`: the plug-in Gini of
`descogini` and Lerman & Yitzhaki's covariance formula. The default keeps
the `n/(n−1)`-corrected Gini of `sp.inequality_index`; shares, Gini
correlations and each source's percentage are identical under both.
- `quantile_convention="rifreg"` for `sp.rifreg`, `sp.rif_values` and
`sp.ffl_decompose`: `rifreg` / `ddecompose` use the indicator `y ≤ q`
where `dineq` uses `y < q`.
- `sp.dfl_decompose(stat_convention="hmisc")`: the weighted variance
(`Σw(y − ȳ)² / (Σw − 1)`) and quantile of `Hmisc`, which
`ddecompose::dfl_decompose` uses for the reweighted counterfactual; with
it the variance and quantile decompositions match `ddecompose` to 1e-12.
The default keeps the reliability-weight variance. Neither is wrong —
under non-uniform weights the two are different estimators — and the
choice moves the composition effect in the fourth significant digit.

⚠️ Correctness fixes — dynamic panels

- **`sp.xtdpdsys` did not reproduce Stata's `xtdpdsys`.** It ran
`xtabond2`'s default moment set — strictly exogenous regressors
instrumenting *both* equations in one combined column, and the `h(3)`
one-step weight — whereas Stata's `xtdpdsys` instruments them in the
differenced equation only and weights with `h(2)` (cross quadrants of `H`
zeroed). On `abdata` (`n` on `L.n w k`) the lagged-dependent coefficient
was 0.686 where `xtdpdsys` reports 0.542, and `w`'s −0.20 where it reports
−0.62. `sp.xtdpdsys` now defaults to Stata's convention and matches it to
1e-12 one-step robust, two-step Windmeijer and classical, instrument
count included. The Stata fixture had carried the `xtdpdsys` covariate
spec all along; no test read it. **Re-run system-GMM results obtained
through `sp.xtdpdsys`.** `sp.xtabond(method="system")` keeps
`xtabond2`'s defaults and is unchanged.

Added — dynamic panels

- `h=` (1, 2, 3; `xtabond2`'s `h()`) and `iv_equation=` (`"both"`, `"diff"`,
`"level"`) on `sp.xtabond` and `sp.xtdpdsys`; `sp.xtdpdsys(...,
iv_equation="both", h=3)` is the previous estimate.
- `xtlsdvc` graded bit-exact against Stata `xtlsdvc` (Bruno): its
coefficient comparison over every initialiser and bias order had been in
the suite for months without a registered promotion. A test now also
reads the `abdata` fixture's `B4` spec — Stata's default `xtabond` with
`_cons` — confirming every slope and SE is unaffected by the constant
StatsPAI's difference GMM omits.
- `sp.from_stata("xtdpdsys ...")` now translates to a runnable
`sp.xtdpdsys(...)` call; it failed loud with a stale "system GMM not
implemented" message.

⚠️ Correctness fixes — production functions

- **`sp.olley_pakes` and `sp.levinsohn_petrin` were not the OP / LP
estimators.** They put every input into a degree-3 polynomial and
estimated the labour coefficient in a GMM second stage that used current
labour as its own instrument. OP and LP read the free-input coefficient off
the stage-1 regression of output on labour and a polynomial in (capital,
proxy), then choose the capital coefficient to minimise the squared
productivity innovations under a cubic Markov polynomial. That is what
Stata and R `prodest` compute. On the parity panel (true labour elasticity
0.60) labour moves from 0.999 (OP) / 1.003 (LP) to 0.633 / 0.615. The
stage-1 coefficient now matches both `prodest` implementations to 1e-12;
the capital coefficient is the exact minimiser, where both references stop
their optimisers early.
- **Lags shifted rows instead of periods.** Every production-function
estimator lagged by position within a firm, so a firm that skipped 2002
paired 2003 with 2001. They now lag by calendar year (`time - 1`, missing
across a gap), as Stata's `L.` and `prodest`'s `lagPanel` do. The
"non-consecutive time periods" warning is gone;
`diagnostics["n_calendar_gaps"]` counts the affected rows. Duplicate
(firm, year) rows or a non-numeric `time` now raise.
- **`sp.ackerberg_caves_frazer` returned an arbitrary local optimum:**
Nelder-Mead on the GMM criterion from five fixed starts, with a linear
Markov process by default. It now solves the just-identified moment
conditions exactly (analytic Jacobian) from the stage-1 coefficients and a
grid of starts. It returns the root nearest the stage-1 coefficients,
lists every root in `diagnostics["acf_roots"]` and warns when there is more
than one. On the parity panel R `prodestACF` stops 2.4e-6 from the same
root; Stata `prodest, acf` stops at points that are not roots.
- **`sp.wooldridge_prod` treated labour as exogenous.** It minimised a
stacked sum of squares in which labour was its own instrument in the
productivity equation, where labour is correlated with the innovation. It
is now a joint GMM:
- equation 1 is instrumented by current labour and the control-function
terms;
- equation 2 by capital, lagged labour and lagged control-function terms;
- the Markov polynomial is free, and SEs are analytic and firm-clustered
(previously NaN without `boot_reps`).

`convention="prodest"` reproduces Stata `prodest, method(wrdg)` to 1e-8,
coefficients and its unadjusted variance. That convention imposes a
unit-slope Markov process, which on the parity panel (rho = 0.7) puts
capital at -0.64 against a true 0.30; the default puts it at 0.34. R
`prodestWRDG` also drops the second equation's intercept and leaves the
constant out of its instruments (capital -0.79).

Changed — production functions

- Default `productivity_degree` is 3 (cubic `g`, as in both `prodest`
implementations) for `olley_pakes`, `levinsohn_petrin`,
`ackerberg_caves_frazer` and `wooldridge_prod`. It was 1, and 2 for
`wooldridge_prod`, whose default `polynomial_degree` is now 3 (was 2).
`sp.prod_fn` passes the degrees through only when they are given.
- `functional_form="translog"` raises `NotImplementedError` for OP and LP,
whose free-input coefficients are linear stage-1 coefficients. ACF keeps
translog.
- `sp.wooldridge_prod` gains `vce=` (`"cluster"`, `"robust"`,
`"unadjusted"`) and `convention=` (`"statspai"`, `"prodest"`).
- The firm bootstrap makes each resampled firm its own panel, so a firm
drawn twice no longer collides with itself in the lag.

Added — production functions

- `tests/reference_parity/test_prodest_parity.py` against Stata 18 `prodest`
(SSC) and R `prodest` 1.0.2 on a simulated unbalanced panel with calendar
gaps (`_fixtures/_generate_prodest_{data.py,R.R,stata.do}`).

⚠️ Correctness fixes — small p-values package-wide

- **335 p-values in 180 modules were computed as `1 − cdf(x)`**, including
the coefficient table of every result built on `EconometricResults`
(`sp.regress` and the estimators that share it). In double precision
`1 − cdf` is accurate only to ~1e-16 *absolute*: a p-value of 1e-10
carries a relative error near 1e-6, and anything below ~1e-16 is returned
as exactly 0.0 (for a normal test, any |z| above ~8.3). They are now the
distribution's survival function, `sf(x)`, which is accurate across the
whole range. Moderate p-values move by at most a few units in the 16th
decimal; tiny ones become correct instead of 0. The rewrite was done on
the syntax tree — `1 − <scipy distribution>.cdf(args)` → `.sf(args)`,
including frozen distributions — so operator precedence is never guessed;
two call sites through lazily imported `scipy.stats` were converted by
hand. Found when the MR sweep's heterogeneity p-value (3e-10) disagreed
with `TwoSampleMR` at 1e-7.

full audit.)*

Added — dispatcher + registry wiring

- `sp.mr(method=...)` routes `mr_lap | lap | sample_overlap`,
`mr_clust | clust | clustered`, `grapple | profile_likelihood`,
`mr_cml | cml | constrained_ml`, `mr_raps | raps |
robust_profile_score` to the new estimators.
- All six new functions (5 MR + `dml_panel`) registered in
`registry.py` with full `ParamSpec` metadata, category, tags, and
reference. `sp.describe_function`, `sp.function_schema`, and
`sp.agent_card` cover them.

Added — tests

- `tests/test_mr_frontier.py` — 41 tests covering correctness,
boundary validation, cross-method consistency (`mr_lap` with
`overlap=0` == IVW; `mr_cml` with `K=0` ≈ IVW; `mr_clust`
two-cluster DGP; `mr_raps` outlier-robustness vs IVW), dispatcher
routing, and registry/schema export.
- `tests/test_dml_panel.py` — 13 tests covering recovery under
homogeneous treatment, FE-OLS agreement in the no-confounding
limit, cluster-SE vs iid SE under AR(1) within-unit correlation,
time-FE option, boundary validation, and registry metadata.

(2025, *JSS* 113(1), DOI 10.18637/jss.v113.i01), and the SDID schema
is canonicalised end-to-end so `sp.synth_report(method='sdid', ...)`
produces a full Markdown / text / LaTeX report rather than a row of
N/As. Seven 2022–2025 SCM citations were added to `paper.bib`,
each verified independently via Crossref / arXiv (refs verified
via `crossref` + `arxiv`). Third, a decomposition-module polish — see
the dedicated section below. Fourth, a **ML+causal polish wave (v1.15)**
covering DML / meta-learners / causal forests / causal discovery /
policy learning / mediation — see the dedicated section below. Fifth,
an **RDD module polish (v1.15)** to the 2018–2026 frontier with three
new estimators (`sp.rd_flex`, `sp.rd_bias_aware_fuzzy`, `sp.rd_discrete`),
three reporting helpers (`sp.rd_dashboard`, `sp.rd_compare`,
`sp.rd_robustness_table`), an `sp.rdrobust` polish pass with the CCT-2018
`rho` parameter and discrete-RV / weak-first-stage warnings, and a
`sp.rdplotdensity` upgrade to the Cattaneo-Jansson-Ma (2020) boundary-
adaptive density estimator — see the dedicated section below.

which the Track A goldens were generated with. With the newer port
installed, `bwselect='cct'` therefore reports a bandwidth about 4e-8
(relative) away from the native one, where the two used to agree to
1e-15. No StatsPAI number changed.

The maintainers of `rdlocrand` released 3.0 on 2026-10-04. Every reference
fixture was regenerated against it. The window sequences StatsPAI already
built (those of release 1.0) are the ones 3.0 returns, to the last digit.
Notes are in `docs/dev/2026-10-04-cit-rd-extensions-review.md`.

⚠️ Correctness

- **`sp.rdrandinf`, `sp.rdwinselect` and `sp.rdsensitivity` with `p > 0`
use HC3 for the large-sample test; they used HC2.** `rdlocrand` 3.0 added
`vce = "HC1" | "HC2" | "HC3"` with HC3 as the default, and the same
argument is new here (`vce='hc1' | 'hc2' | 'hc3'`). All three agree with
3.0 to 7e-13 across polynomial orders, kernels, a non-zero null and
evaluation points away from the cutoff. On the U.S. Senate data, for
`termshouse` with `p = 1` in [-1, 1], the large-sample p-value moves from
0.666 to 0.694.
`vce='hc2'` returns the earlier numbers. Calls with `p = 0`, the default,
are unchanged, and so is the randomization p-value.

Changed

- `sp.from_stata` / `sp.stata` read `vce()` on `rdrandinf` and
`rdwinselect`. With `p() > 0` and no `vce()` the translation notes that
releases up to 2.0 used HC2.
- `tests/reference_parity/test_rdlocrand_v1_parity.py` runs against both
1.0 and 3.0. `sp.rdwinselect(approx=True, p=1)` has reference numbers for
the first time (2.0 stopped with an error there).
- With `p > 0`, 3.0 no longer computes a randomization p-value.
`sp.rdrandinf` keeps its own, which permutes outcomes against (score,
assignment) pairs and holds its level in simulation.

What the examples of Wooldridge's *Introductory Econometrics* found

Every example of the book that the companion scripts of Heiss compute was
written with StatsPAI and compared with statsmodels, linearmodels, Stata 18
or R on the book's datasets. Most of it agreed to 1e-12. The notes are in
`docs/dev/2026-10-05-wooldridge-review.md`; the Stata and R numbers are in
`tests/test_wooldridge_examples.py`.

⚠️ Correctness

- **`sp.synth(method='classic')` reported a confidence interval that could
contradict its own p-value.** `pvalue` is the rank of the treated
post/pre RMSPE ratio among the placebos; `ci` was the estimate plus or
minus a normal quantile times the standard deviation of the placebo
ATTs, a different procedure with no guarantee of agreeing with the test.
`ci` is now the interval the test implies: the constant effects that
the same rank test does not reject (Firpo and Possebom 2018), with exact
ends. It excludes zero exactly when `pvalue < alpha`. On Proposition 99
without covariates (p = 3/39 = 0.077) the 95% interval goes from
[-40.7, 1.6] to [-58.6, 19.6] and the 90% one from [-37.3, -1.8] to
[-33.0, -5.3]. With fewer than `1 / alpha - 1` donors no effect can be
rejected and the interval is `(-inf, inf)`; it used to be finite. The
estimate, `se` and `pvalue` are unchanged, the old interval is in
`model_info['ci_normal']`, and `model_info['ci_method']` says which one
`ci` is. `summary()` now prints an infinite end as `-inf` / `inf`
instead of leaving it blank, for every estimator. See `MIGRATION.md`.
- **DFBETAS from `sp.estat(result, 'leverage')` were on the wrong scale.**
The division by the square root of the diagonal of `(X'X)^-1` was
missing, so each value was the change in the coefficient over the
leave-one-out residual standard deviation, not over the standard error of
the coefficient. The `2/sqrt(n)` rule was then applied to that quantity,
and `dfbetas_flagged_obs` listed the wrong observations. The values now
agree with Stata's `dfbeta` and with statsmodels to 1e-8. Leverage and
Cook's distance were right and are unchanged. See `MIGRATION.md`.
- **`predict(new_data)` after `sp.glm` returned the linear index.** With
no argument the same method returned the fitted means, so the two calls
were on different scales for every link but the identity. Both are now
the mean of the outcome; `what='link'` gives the index. After
`sp.poisson` and `sp.nbreg` the call used to fail on the `_cons` column
and now works, with the `offset=` or `exposure=` of the new rows applied.
See `MIGRATION.md`.
- **`sp.from_stata` translated `heckman` to the two-step estimator.**
Stata's `heckman` is maximum likelihood unless `twostep` is given, and
`twostep` itself was reported as an untranslated option. The translation
now writes `method='ml'` or `method='twostep'`. On the `mroz` data the
coefficient on education is 0.10835 by maximum likelihood and 0.10907 by
the two-step estimator; both now match Stata. `sp.heckman` itself is
unchanged. See `MIGRATION.md`.
- **`sp.estat(result, 'endogenous')` and `'overid'` said "p-value not
available"** in their interpretation, next to the p-value they reported.
They now state the conclusion.

Added

- **Formulas.** `log`, `log2`, `log10`, `log1p`, `exp` and `sqrt` can be
written without the `np.` prefix in every formula (`log(wage) ~ educ`).
A transformed outcome is accepted by `sp.ivreg` and `sp.panel`, a
transformed endogenous regressor by `sp.ivreg`
(`hours ~ educ + (log(wage) ~ exper + I(exper**2))`), and `sp.qreg`
takes transformed, factor and interaction terms. A name the data do not
have is reported by name; it used to end in a `KeyError`.
- **`result.predict(data, what=, alpha=)`** for the results of
`sp.regress`, `sp.ivreg`, `sp.glm`, `sp.poisson` and the other
formula-based fits. The design is rebuilt from the formula, so new data
need the raw variables only. `what='confidence'` returns the standard
error and interval of the conditional mean, `what='prediction'` the
interval for a new observation of a linear model, `what='link'` the
index. The results of `sp.logit` and `sp.probit` take a DataFrame.
- **`sp.estat`**: `'archlm'` (Engle's LM test), `'durbinalt'` (Durbin's
alternative test), `'bgodfrey'` with `version='fstat'`, and `'white'`
with `variables='fitted'` (the form of the test on the fitted values and
their squares). `'leverage'` also returns `rstudent` and `rstandard`.
`'beta'` returns the standardized coefficients of Stata's
`regress, beta`.
- **`sp.tobit(formula=)` and `sp.truncreg(formula=)`**, with transformed,
factor and interaction terms; the formula may also come first.
`sp.survreg` accepts the data first.
- **`sp.glm(scale=)`**: `'x2'`, `'dev'` or a number, the factor on the
model-based covariance. `family='poisson', scale='x2'` is the
quasi-Poisson covariance and matches Stata `glm, scale(x2)`.
- **`sp.test`** takes a list of restrictions and coefficient names that
hold operators or blanks (`I(exper ** 2)`, `np.log(sales)`); so does
`sp.lincom`.
- **`sp.lrtest`** accepts the maximum-likelihood regressions (`sp.logit`,
`sp.probit`, `sp.poisson`, `sp.nbreg`, `sp.tobit`, `sp.glm`). It refuses
fits on different samples, different models, and a reversed order.
- **Stata translations**: `estat archlm`, `estat durbinalt`, `truncreg`,
`glm` (`family()`, `link()`, `scale()`, `vce()`), and `vce(robust)` /
`vce(cluster)` on `heckman`.

Fixed

- `sp.panel(method='fd')` with the time variable among the regressors (a
linear trend) raised `KeyError`.
- `sp.panel(method='fe')` and `'twoway'` raised when the effects absorbed a
regressor. The regressor is now omitted with a warning and listed in
`model_info['omitted']`, as Stata and R do.
- `sp.panel(method='mundlak')` raised on a time-invariant regressor or on
period dummies. Unit means that add no rank are left out and listed in
`model_info['cre_means_omitted']`.

⚠️ Correctness

- **`sp.stata` / `sp.from_stata` translated `ivregress gmm` to a
different estimator.** Stata's two-step GMM uses the heteroskedasticity
-robust weight matrix unless told otherwise; the translation called
`sp.ivreg(method='gmm')`, which uses the unadjusted one, and dropped
`wmatrix()`. On Card's schooling data the schooling coefficient came
out 0.1608 against Stata's 0.1588. The translation now follows
`wmatrix()` and `vce()` as Stata 18 resolves them (a `vce(cluster c)`
alone sets the weight matrix too), and nine option combinations
reproduce `ivregress gmm` to 5e-9 through `sp.stata`
(`tests/reference_parity/test_iv_gmm_stata_parity.py`). Combinations
with no counterpart (`wmatrix(hac ...)`, `igmm`, a cluster weight matrix
with a robust variance, `wmatrix(unadjusted) small`) are reported as
untranslated and `sp.stata` refuses them. Calls to `sp.iv` /
`sp.ivreg` themselves return the same numbers as before.
- **`sp.ipw(estimand='ATT' | 'ATC', normalize=False)` was scaled by the
share of the target group.** The Horvitz-Thompson sums were divided by
`n` for every estimand, so the ATT came back multiplied by `P(T=1)` and
the ATC by `P(T=0)`. They are now divided by the number of treated and
of controls. On the NHANES school-meal data of Ding (2024, chapter 13)
the ATT goes from -1.098 to -1.992, the value of the book's own code.
`normalize=True` (the default) and the ATE are unchanged. See
`MIGRATION.md`.
- **The confidence interval of `sp.fisher_exact` was too short.** It was
read off a 101-point grid six outcome standard deviations wide, each
point tested on its own 500 draws, and with `controls` the shifted
outcome was not residualized. It is now the exact inversion of the
randomization test on the assignments that gave the p-value, each end
found by bisection. On the LaLonde experiment the 95% interval goes
from [999, 2589] to [572, 3015]; effects of 1000 and 2600 have
randomization p-values of 0.21 and 0.20 and belong in it. The interval
is infinite, with a warning, when the design has too few assignments to
reject at `alpha`. p-values are unchanged. See `MIGRATION.md`.
- **`sp.gformula_ice_fn` raised `LinAlgError` on a flat confounder list.**
The documented shorthand for "the same confounders at every time", and
any covariate named at more than one time, entered the history twice and
made the default sandwich standard error singular. Repeated names are
now counted once. Point estimates are unchanged.
- **Lee and Zhang-Rubin bounds warn when ties defeat the trimming.** On a
binary or coarsely discrete outcome `trimming='quantile'` keeps every
observation tied at the quantile, trims almost nothing, and the two
bounds collapse (both -0.045 where the sharp bounds are [-0.176,
-0.019] on the trial of Ding's chapter 26). `sp.lee_bounds`,
`sp.principal_strat`, `sp.survivor_average_causal_effect` and
`sp.attrition_bounds` now say so and point to `trimming='exact'`. The
default and its numbers are unchanged.
- **For `statistic='ks'` and `'rank_sum'` the "confidence interval" of
`sp.fisher_exact` was not an interval for the effect.** It was the 2.5%
and 97.5% percentiles of the null distribution of the statistic, centred
on zero whatever the effect. It is now the same constant-effect
inversion as for the difference in means, in outcome units: the
statistic is recomputed on `Y - tau_0 * D` for a stored set of
assignments and each end is bisected. With 12 treated among 60 units
and a constant effect of 0.75 the 95% interval covers 96.0% (rank sum)
and 96.7% (KS) of the time over 150 draws. The same holds for
`statistic='t'`. p-values are unchanged. See `MIGRATION.md`.
- **Arithmetic inside a formula wrapped around on small integer
columns.** `pd.read_stata` returns a Stata `int` as `int16` and a `byte`
as `int8`, and numpy keeps arithmetic in that type, so `I(x**2)` was a
column of garbage for values above 181 (or 11) and the regression was
fitted on it without a warning. On the restaurant-inspection data of
Huntington-Klein's *The Effect* (a count up to 646) the coefficient of
the quadratic model was -0.0171 where Stata and R give -0.0844. It
affected `sp.regress`, `sp.feols`, `sp.fepois`, `sp.feglm`, `sp.glm`,
`sp.logit`, `sp.ivreg`, `sp.panel`, the covariate formulas of the DiD
family and the spatial models. Narrow integer columns are now widened to
`int64` before a formula is evaluated
(`tests/test_narrow_integer_formula.py`). Frames whose integers are
already `int64` or float are unaffected. See `MIGRATION.md`.
- **`sp.ebalance(moments=2)` with a 0/1 covariate returned uniform
weights.** The square of an indicator is the indicator, the dual was
singular, and the estimate reported was the raw difference in means
behind a `ConvergenceWarning`. A moment that repeats another one is now
left out (`model_info['redundant_moments']`) and the weights balance the
rest exactly.
- **`sp.match(method='cem')` follows the `cem` packages.** Sturges' rule
gives the number of cut points, hence one bin fewer than before, and a
bin is closed on the right. The default now reproduces Stata `cem` on
the NSW data (matched sample and estimate). Estimates from the default
change, and so do estimates with integer covariates under an explicit
`n_bins`. See `MIGRATION.md`.
- **`sp.sensemakr` skipped a benchmark name it did not find** and returned
a table without that row. It is now an error.
- **`sp.stata` / `sp.from_stata`, four silent mistranslations of factor
notation.** `abc` was written as the main effects and the three-way
product, without the two-way terms. A variable with no prefix inside ``
or `` was read as continuous where Stata reads a factor. `ab` without
main effects was written as the product of two columns; it is now
refused, because the formula's coefficients would not be Stata's.
`i.pc.xc.x` lost the square. All four now follow Stata and are
pinned against Stata 18
(`tests/reference_parity/test_the_effect_stata_parity.py`).
- **`sp.stata`: `collapse` stores a mean, median or sd as Stata does**
(float unless the source is a long or a double). Regressions on
collapsed data differed from Stata in the sixth digit.

- **`sp.survreg` ignored `robust=` and `cluster=`.** Both were accepted
and written into `model_info`, and the standard errors were the
observed-information ones whatever was asked. They are now the sandwich
and the cluster sandwich, with Stata's factors, and agree with Stata 18
`streg ..., time vce(robust)` and `vce(cluster)` to 1e-6 for the
Weibull, exponential, log-normal and log-logistic models
(`tests/reference_parity/test_survreg_streg_stata_parity.py`). On the
test data the robust standard error of the Weibull scale is 28% larger
than the one that used to be returned under `robust='robust'`. Fits
without these options change in the sixth digit only: the optimiser now
finishes with Newton steps (coefficients were about 1e-5 from the
optimum) and agrees with Stata to 1e-9. See `MIGRATION.md`.
- **The Vuong statistic reported by `sp.zip_model` and `sp.zinb` was too
large.** `diagnostics['vuong_stat']` compared the zero-inflated fit with
a Poisson (negative binomial) density evaluated at the zero-inflated
model's own count coefficients. That is not the maximum of the
comparison model's likelihood, so the comparison model looked worse than
it is. It is now the regression fitted on its own, as in R
`pscl::vuong` and in Vuong (1989). On the test data the statistic goes
from 13.42 to 8.69 for the Poisson pair and from 9.16 to 3.70 for the
negative binomial pair, and both now agree with R and with Stata 18
(`forcevuong`) to 1e-6. Coefficients,
standard errors and log-likelihoods are unchanged. The docstrings also
say that this statistic is not a valid test for zero inflation (Wilson
2015). See `MIGRATION.md`.

Added

From a pass over Huntington-Klein, *The Effect* (2nd ed.): its Stata code
for chapters 4 and 13 to 21 was run in Stata 18 and replayed through
`sp.stata` (836 printed numbers reproduced, 148 before). Notes are in
`docs/dev/2026-10-05-the-effect-2e-review.md`.

- `sp.ebalance(moments=[2, 2, 1])` takes one order per covariate (Stata
`targets(2 2 1)`), and `dof_adjust=True` matches the sample variance
and skewness as Stata `ebalance` does. With it the weights agree with
`ebalance, tolerance(1e-10)` to 1e-10. The default is unchanged.
- `sp.match(method='cem', n_bins=...)` takes a list or a dict, from
covariate to a number of bins or to cut edges.
- `sp.sensemakr(benchmark={'label': [columns]})` benchmarks several
controls as one group (`gbenchmark()` / `gname()` in Stata). Agrees with
Stata `sensemakr` to 1e-12.
- `sp.rdrobust` names a covariate, outcome or running variable that is not
numeric instead of failing inside numpy.
- `sp.rosenbaum_bounds(estimates=True)` adds, for each Gamma, the range of
the Hodges-Lehmann point estimate and the outer ends of its confidence
interval (the `t-hat` and `CI` columns of Stata `rbounds`). The
significance levels and the confidence bounds agree with Rosenbaum's
`DOS2::senWilcox`, the point-estimate bounds with Stata `rbounds`.
`sp.stata` runs `rbounds diff, gamma()` and `bstat, stat()`.
- `sp.stata` runs `xi`, `ebalance`, `cem`, `sensemakr` and
`table g, statistic(mean x)`; `encode, g()`;
`date()`; `egen cut(), group()`; `reghdfe` with `dib3.t` and squares;
`_b[1.d3.t]`; abbreviated variable names; weights that are missing or
zero (the observation is dropped, as in Stata); `[iw = w]` in `regress`
when the weights sum to the number of observations; `i.` covariates in
`psmatch2` and `teffects`; `program def`, `return scalar` of a local
macro and nested programs. Translated `tabstat, by()` and `table` head
the groups by their values; they used to read "Control" and "Treated"
for any 0/1 grouping variable.

Seven functions from a pass over Croissant (2025), *Microeconometrics with
R*, and its companion package `micsr`. The book's chapters on binary,
count and censored outcomes use estimators StatsPAI did not have. Notes
are in `docs/dev/2026-10-04-croissant-microeconometrics-review.md`.

- **`sp.ivprobit` and `sp.ivtobit`** fit a probit and a tobit with
continuous endogenous regressors. Maximum likelihood is the default and
`method='twostep'` gives Newey's minimum chi-squared estimator, as in
Stata's `ivprobit` and `ivtobit`. Both report the Wald test of
exogeneity and take `vce='robust'` or `cluster=`. Coefficients,
standard errors, the log-likelihood and both Wald statistics agree with
Stata 18 to 1e-6 on 14 blocks, with one and two endogenous regressors
and with one and two censoring limits
(`tests/reference_parity/test_ivprobit_ivtobit_stata_parity.py`).
- **`sp.ivpoisson`** fits an exponential-mean model with endogenous
regressors by GMM, as Stata's `ivpoisson gmm`. `errors='multiplicative'`
is Mullahy's (1997) moment condition, the one that is valid when an
omitted variable enters the exponential. One-step, two-step and iterated
GMM, robust, clustered and unadjusted weight matrices, Hansen's J. It
agrees with Stata 18 to 1e-6 on 12 blocks
(`tests/reference_parity/test_ivpoisson_stata_parity.py`).
- **`sp.cmtest`** runs the conditional moment tests of normality,
heteroskedasticity, skewness and kurtosis on a fit from `sp.probit` or
`sp.tobit`. Those models are inconsistent when the latent error is not
normal or not homoskedastic, and they leave no residual to inspect. The
statistics agree with R `micsr::cmtest` to 1e-9 for the tobit and 1e-6
for the probit, and reject 5% of the time under the null in simulation
(`tests/reference_parity/test_cmtest_micsr_parity.py`).
- **`sp.nlogit`** fits the nested logit in its random-utility form on
long-format choice data, with one dissimilarity parameter per nest or a
common one, alternative-specific constants, and the likelihood-ratio
test of IIA against the conditional logit. Coefficients, standard
errors (oim, robust, cluster) and the test agree with Stata 18 `nlogit`
to 1e-6, and coefficients and likelihoods with R `mlogit`. `mlogit`'s
standard errors differ from Stata's by up to 4% on the same estimates;
StatsPAI's are Stata's (`tests/reference_parity/test_nlogit_parity.py`).
- **`sp.etpoisson`** fits a Poisson regression with an endogenous binary
treatment by maximum likelihood, as Stata's `etpoisson`. The outcome
error is integrated out by Gauss-Hermite quadrature (`intpoints=24`).
It reports the correlation of the two errors, the Wald test of
independent equations and the average treatment effect on the count
scale with a delta-method standard error. Coefficients and standard
errors agree with Stata 18 to 1e-6 (observed 1e-10) and the ATE with
`margins r.treat` to 1e-9
(`tests/reference_parity/test_etpoisson_stata_parity.py`).
`sp.from_stata` translates the command.
- **`sp.tobit(method='scls')`** is Powell's symmetrically censored least
squares. It needs a symmetric error and nothing else, so it stays
consistent under heteroskedasticity and non-normality, which is what to
turn to when `sp.cmtest` rejects. In a simulation with heteroskedastic
errors its slope averages 1.01 for a true value of 1 and its 95%
interval covers 97% of the time, where maximum likelihood averages 1.33.
Coefficients agree with `micsr::tobit1(method = "trimmed")` to 1e-8
(`tests/reference_parity/test_binary_r2_and_scls_parity.py`).
- **`sp.probit` and `sp.logit` report a family of goodness-of-fit
measures** in `model_info['r2']`: McFadden, Cox-Snell, Nagelkerke,
Efron, Tjur, Estrella and McKelvey-Zavoina. They agree with R
`DescTools::PseudoR2` and `micsr::rsq` to 1e-6. `model_info['pseudo_r2']`
is unchanged.
- **`sp.survreg(frailty='gamma')`** adds gamma-distributed unobserved
heterogeneity to the four parametric duration models, as Stata's
`streg, frailty(gamma)`. It reports `log(theta)`, the frailty variance
and the likelihood-ratio test of `theta = 0` with its boundary p-value.
Coefficients, standard errors, `theta` and the test agree with Stata 18
to 1e-6. When the variance goes to its boundary the fit is the model
without frailty and a warning says so. `sp.survreg` results also carry
`data_info['llobs']`, so `sp.vuong` can compare two duration
distributions, and `model_info['ll_log_time']`, the log-likelihood in
Stata's convention.
- **`sp.vuong`** compares two non-nested maximum likelihood fits on the
same observations and returns the Vuong statistic with its AIC and BIC
corrections. It takes fits from `sp.poisson`, `sp.nbreg`,
`sp.zip_model`, `sp.zinb`, `sp.hurdle`, `sp.logit` and `sp.probit`,
which now carry their per-observation log-likelihood in
`data_info['llobs']`. Statistics agree with R `pscl::vuong` to 1e-6 on
five model pairs. `pscl` leaves a negative binomial dispersion out of
its parameter count and StatsPAI does not, so the corrected statistics
differ when exactly one model is negative binomial
(`tests/reference_parity/test_vuong_pscl_parity.py`).
- `sp.stata` reads `stset timevar, failure(eventvar)` and runs
`streg ..., distribution() time [frailty(gamma)] [vce()]` through
`sp.survreg`. `streg` for the Weibull or exponential without `time` is
the proportional-hazards metric and is refused, as are `shared()`,
other frailty distributions, and `stset` with `id()`, `enter()` or
`origin()`.
- `sp.ologit`, `sp.oprobit`, `sp.mlogit`, `sp.tobit` and `sp.nlogit` fits
carry their per-observation log-likelihood, so `sp.vuong` takes them.
- `sp.from_stata` and `sp.stata` translate `ivprobit`, `ivtobit` and
`ivpoisson gmm`. `ivtobit` without `ll()` is written as `ll=None`,
because `sp.ivtobit` censors at zero by default and Stata does not.
`ivpoisson cfunction` is refused.

What the syllabus of Xu and Lan's *Causal-Inference Econometrics* found

徐小君、蓝嘉俊《因果推断计量经济学》(清华大学出版社, 2025) publishes no data
and no code, only a sample chapter and the table of contents. The twelve
chapters were therefore used as a syllabus. One Stata do-file per chapter
runs the commands those methods are taught with, on datasets that ship with
Stata, and each log was replayed through `sp.stata` number by number
(`tests/external_parity/xu_lan_syllabus/`). About 1,200 printed numbers
were reproduced at the start and two blocks were not. 85 commands were
refused, and several of the functions behind them had never been compared
with Stata. Three of the findings were made by other passes the same
night and fixed on main first: `ivregress gmm` translating to 2SLS,
`heckman` translating to the two-step estimator, and the missing
`estat durbinalt` / `estat archlm`. Those fixes are theirs (see the
entries above); this pass keeps a second set of Stata numbers for each.
Notes are in
`docs/dev/2026-10-05-xu-lan-causal-econometrics-review.md`; the numbers
below are pinned on committed synthetic data in
`tests/reference_parity/test_textbook_syllabus_stata_parity.py` (40 tests).

**⚠️ Correctness**

- **`sp.arima` fitted a zero-mean model to a series that has a mean.** The
default (`method='statespace'`) estimated no constant, so an AR(1) on a
series with a non-zero mean pushed the autoregressive coefficient towards
one. On the first difference of the log U.S. wholesale price index the
coefficient was 0.763; with the mean estimated it is 0.619 (Stata: 0.615
by exact likelihood). `sp.arima` now has `trend=`. The default estimates
a constant when the series is not differenced and none when it is, as R's
`stats::arima` and statsmodels' `ARIMA` do; `trend='c'` on a differenced
series estimates the drift Stata's `arima` always includes; `trend='n'`
gives the old fit. `method='innovations_mle'` already estimated the
constant and is unchanged. See `MIGRATION.md`.
- **Two-step standard errors of `sp.heckman` now use the observed
information of the first-step probit.** They used the expected
information. The difference is in the fourth digit (0.60627 to 0.60654
for the inverse Mills ratio on Stata's `womenwk`), and it was the whole of
the gap that Track A module 43 carried as a tolerance. The worst relative
gap in a standard error against R `sampleSelection` and Stata goes from
8.6e-5 to 4e-9, and the module's budget is now the default 1e-6.
- **`sp.stata("tobit y x, ul(2)")` censored from below as well.** Stata
censors only on the sides that are written. The translation passed the
upper limit and left `sp.tobit` at its own default of `ll=0`, so every
outcome at or below zero was treated as censored too. `tobit y x` with
no limit (an uncensored normal model in Stata) became a tobit at zero,
and `ll` without a value or with a variable name fell back to zero
silently. Both limits are now written out, `sp.tobit` accepts
`ll=None`, and the forms that cannot be translated are refused. On a
sample censored from above only, the slope goes from 0.965 to 0.998
(truth 1), equal to Stata 18 to 1e-8. See `MIGRATION.md`.
- **`sp.tab` reported a Yates-corrected statistic as "Pearson chi2" on
2 x 2 tables.** The test came from `scipy.stats.chi2_contingency` at its
default, which applies the continuity correction when the table has one
degree of freedom. The label, Stata's `tabulate, chi2` and R's
`chisq.test(correct = FALSE)` all mean the uncorrected statistic. On the
test table it is 11.84, not the 11.11 that was printed. Larger tables
are unchanged. See `MIGRATION.md`.
- **`sp.garch(p >= 1, q=0)` returned numbers for a model that is not
identified.** With no ARCH term the variance never responds to the data
and every value of the GARCH coefficient gives the same likelihood. The
call is now refused. (`p` counts lagged variances and `q` lagged squared
innovations; ARCH(1) is `p=0, q=1`.)

**Added**

- **`sp.sdtest` and `sp.ztest`.** The chi-squared test that a standard
deviation equals a value, the F test that two are equal, and the z test
on a mean or a difference of means when the standard deviation is known.
These are the two worked examples of the book's first chapter. Both take
data or summary statistics (`sp.sdtest(n=10, sd=1.14, sd0=2)`), as
Stata's `sdtest` / `sdtesti` and `ztest` / `ztesti`.
- **`sp.prtest`, `sp.sktest`, `sp.swilk` and `sp.ci`.** The rest of a
first statistics chapter: the z test on one or two proportions, the
skewness-kurtosis and Shapiro-Wilk tests for normality, and confidence
intervals for means, variances, standard deviations and proportions
(exact, Wald, Wilson, Agresti-Coull, Jeffreys). They agree with Stata 18
`prtest`, `sktest`, `swilk` and `ci` to 1e-8 or better (the Shapiro-Wilk
p-value to 1e-5). `sp.stata` runs those commands and `prtesti`.
- **`svar ..., aeq(A) beq(B)` and `svar ..., lreq(C)`** run in `sp.stata`,
with the matrices defined by `matrix A = (...)` before them.
- **`tabulate a b, chi2 exact lrchi2 V`** runs in `sp.stata`. The Pearson
and likelihood-ratio statistics, Cramer's V and Fisher's exact p-values
are in `.attrs['test']` of the table and in `r()`; `sp.tab(...,
output='dataframe')` carries the same dictionary.
- **`sp.nlcom`.** A nonlinear function of the coefficients with a
delta-method standard error: ratios, the long-run effect `b / (1 - c)` of
a partial-adjustment model, turning points. The expression is parsed into
a syntax tree and differentiated node by node; nothing is passed to
`eval`.
- **`sp.svar`.** Structural VAR on a fitted `sp.var`: short-run
restrictions in the AB model (`A=`, `B=`), long-run restrictions
(`long_run=`), or sign restrictions on the impulse responses (`sign=`).
Coefficients, standard errors, the log-likelihood, the test of
over-identifying restrictions, structural impulse responses and variance
decompositions agree with Stata 18 `svar` to 1e-6 in all three
likelihood cases. Sign restrictions return the set of admissible
rotations, and the documentation says what its bands are not: they carry
no sampling uncertainty. `VARResult.fevd()` gives the Cholesky variance
decomposition.
- **`sp.unitroot(test='pp' | 'kpss')`.** Phillips-Perron `Z(t)` and
`Z(rho)` and the KPSS statistic, equal to Stata's `pperron` and `kpss` to
1e-9. KPSS has stationarity as its null; the result says so in `null`.
- **Stata's factor-variable names in `sp.test` and `sp.lincom`.**
`"2.race = 3.race"` and `"x + 1.unionc.x"` find the coefficients the
formula engine calls `C(race)[T.2]` and `x:C(union)[T.1]`. `"i.race"` on
its own tests every level jointly.
- **`sp.heckman(select=None)`** selects where the outcome is observed
(Stata's `heckman y x, select(z1 z2)`, which `sp.stata` now runs), and
`method='mle'` is accepted for `'ml'`. `sp.arima` and `sp.garch` take
`data=` with a column name.
- **`sp.hdfe_ols(absorb_dof='areg')`** (and `sp.absorb_ols`) charges every
absorbed level in a clustered variance, where the default follows
`reghdfe` and does not charge an effect nested in the cluster variable.
With `drop_singletons=False` the fit is then the dummy-variable
regression, slopes and standard errors equal to
`sp.regress('y ~ x + C(g)')` to 1e-10 under the classical, robust and
clustered variances, with and without weights, without building the
dummies. `sp.stata` sends `areg` there: on the NLS panel (4,134 groups)
the command goes from about three minutes to under a second, and its
`_cons` is now Stata's (the dummy-variable translation reported the
first group's level and said so in a note). The result no longer lists
the group dummies. A robust variance with weights and HC2 / HC3 keep
the dummy-variable translation.
- **`sp.ttest(n=, mean=, sd=)`** takes summary statistics for one sample
or two independent ones (Stata's `ttesti`).
- **Abbreviated variable names in `sp.stata`** are spelled out in
`heckman`, `truncreg`, `etregress` and `prais`, and inside the options
that hold a varlist: `select()`, `treat()`, `absorb()`, `by()`,
`cluster()` and `vce(cluster ...)`. An abbreviation that fits two
variables is refused, as in Stata.
- **`sp.stata` runs** `sdtest`, `sdtesti`, `ztest`, `ztesti`, `ttesti`,
`etregress`, `pperron`, `kpss`, `arima`, `arch`, `testparm` and `nlcom`.

**Fixed**

- **`xtreg y x i.year, fe` through `sp.stata` failed on a panel with fewer
units than regressors.** The fixed-effects statistics went through the
between regression of the random-effects variance components, which does
not exist there.
- **`sp.panel(method='mle')` took minutes on a few thousand panels.** It
called the general mixed-model optimiser. The likelihood of the
random-intercept model is now concentrated in one variance ratio. On the
NLS panel (4,134 women, 19,007 rows) the fit goes from more than 150
seconds to 0.3, with the same estimates to 1e-8.
- **`sp.regress` with thousands of dummy columns spent minutes checking
for collinearity.** The scan took one dot product per pair of columns
and regrew its Cholesky factor one row at a time. It now works on blocks
of the cross-product. With about 4,100 dummy columns the two checks go
from more than 400 seconds to about 40; the omitted columns are the
same.

What Gaillac and L'Hour's *Machine Learning for Econometrics* found

A pass over the companion code of Gaillac and L'Hour (2025, OUP): notebooks
and R scripts for chapters 2 to 11. Notes, what was run and what is left
open are in `docs/dev/2026-10-05-gaillac-lhour-ml4econometrics-review.md`.
Most of the book already had a counterpart that reproduces its references
(`hdm::rlassoIV`, `ivmodel` LIML and Fuller on the Card and BLP data to
1e-8; the double-selection simulation of chapter 5 gives bias 0.017 and
coverage 0.949 against the book's 0.012 and 0.942). Three things did not.

- **Cluster-Lasso for panels (chapter 7): `cluster=` on `sp.rlasso`,
`sp.rlasso_effect`, `sp.rlasso_effects` and `sp.rlasso_iv`.** The
penalty loadings of Belloni, Chernozhukov, Hansen and Kozbur (2016) sum
the scores within cluster before squaring, `gamma` defaults to
`0.1 / log(G)`, and the final standard error is cluster-robust. hdm has
no such option and the book sources an unpublished `rlasso_cluster.R`.
Against Stata's lassopack on the same bytes: the penalty level, the
selected variables and the post-Lasso coefficients of `rlasso,
cluster()` agree to 1e-10, and the point estimates and variances of
`pdslasso, cluster()` (double selection and partialling out) and
`ivlasso, cluster()` (selection among instruments, and among both
instruments and controls) agree to 1e-9
(`tests/reference_parity/test_cluster_lasso_stata.py`). Variances carry
no small-sample factor, as in `pdslasso`. What ignoring the clustering
costs: in a within-transformed panel of 100 units and 8 periods with
AR(1) regressors and errors (coefficient 0.8), double selection without
`cluster=` keeps 6.1 controls on average and its 95% interval covers
78.6% of the time; with `cluster=` it keeps 4.6 and covers 93.3% (1,167
draws). On one of the two test panels lassopack and hdm's iteration stop
at different supports. The test shows that both are fixed points of the
loading iteration under StatsPAI's cluster loadings, so the difference
is the path and not the penalty; this also happens without clustering.
- **Classic synthetic control: interval and p-value from one procedure
(chapter 10).** The book's end-of-chapter question 5 is that a test
which rejects and an interval containing zero cannot come from the same
procedure. In `sp.synth(method='classic')` they did not. See the
correctness entry above; the inversion is
`statspai.synth._core.placebo_inversion_ci`, checked against a
brute-force grid. Under `H0: effect = C` the weights do not move, the
treated gap becomes `gap - C` and the gap of a placebo whose donor pool
holds the treated unit with weight `w` becomes `gap + w C`, so each
comparison is a quadratic inequality in `C`.
- **`sp.datasets.nsw_dw()` said it was the Dehejia-Wahba data.** It is a
simulated replica with a built-in latent effect of $1,794 (the treated
mean of `re78` is 5,508; in the real file it is 6,349). The docstring
now opens with that and the frame carries `attrs['simulated'] = True`.
`sp.datasets.nsw_dw(simulated=False)` returns the real sample: the 185
NSW treated and the 2,490 PSID-1 comparison men, bundled, no network.
Its naive difference in 1978 earnings is -$15,205, Table 3 of Dehejia
and Wahba (1999), and it agrees row by row with the files on Dehejia's
NBER page to the rounding of earnings. The bare call is still the
replica, rows unchanged, because parity fixtures are built on it.

Changed

- `sp.survreg` raises `MethodIncompatibility` when a duration is zero or
negative. It used to add `1e-15` to every duration and fit.

Fixed

- **`sp.zip_model` and `sp.zinb` returned NaN for every standard error
when one coefficient of the inflation equation was not identified.**
Under quasi-complete separation that coefficient runs off, the logistic
function overflowed inside the complex-step Hessian, and one NaN spread
to the whole matrix. On the `trips` data of the book, R
`pscl::zeroinfl` reports standard errors for the other 19 coefficients
and StatsPAI's are now within 0.5% of them. `sp.hurdle` was not
affected in its numbers. All three now warn and name the coefficient
the likelihood is flat in (`tests/test_zeroinflated_separation.py`).
Estimates on data without separation do not change.

What Ding's *A First Course in Causal Inference* found

The 28 R programs of the book write every estimator out by hand. Each
deterministic number was recomputed with the matching `sp.*` call
(`tests/external_parity/test_ding_first_course.py`, 19 tests over 17
chapters; review in `docs/dev/2026-10-05-ding-first-course-review.md`,
guide in `docs/guides/ding_first_course.md`). The four corrections are
listed above. What was added:

- `sp.rosenbaum_bounds(method='t')`: Rosenbaum's sensitivity analysis for
the mean pair difference. Agrees with `sensitivitymw::senmw(method =
"t")` to 2e-14.
- `sp.anderson_rubin_test(ar_vcov='HC0' | 'HC1' | 'HC2' | 'HC3')`: the
Anderson-Rubin test and interval with a heteroskedasticity-robust
variance (the Fieller-Anderson-Rubin interval of the book). The
statistic was the homoskedastic F unless `cluster=` was given, and
`vcov` only ever chose the effective F reported beside it. Agrees with
`sandwich::vcovHC` on the reduced form to 1e-9 for one and two
instruments. The default stays `'classic'`.
- `sp.ri_test(stat='rank_sum' | 'lin' | 'lin_t')` and
`sp.fisher_exact(statistic='t')`. `'lin_t'` is Lin's estimate over its
HC2 standard error, the covariate-adjusted statistic whose
randomization test is also valid for the weak null.
- `sp.randomize(rerand_accept=)`: rerandomization by Morgan and Rubin's
criterion. Assignments are redrawn until the Mahalanobis distance,
scaled by `n1 n0 / n`, is below the chi-square quantile of the
acceptance probability; `result.rerandomization` records it. The
`n_rerand` search is unchanged, and its documentation now says that its
threshold is not on that scale.
- `sp.lm_lin(blocks=)`: Lin's estimator fitted in each block of a
stratified experiment and combined by block shares.
- `sp.fisher_exact(...).summary()` prints the confidence level that was
asked for and no longer calls the null-distribution quantiles of
`statistic='ks'` / `'rank_sum'` / `'t'` a confidence interval.
`FisherResult` has an `alpha` attribute.
- `sp.principal_strat(method='principal_score')` recognises one-sided
noncompliance: the arm in which the stratum variable is constant is not
fitted and its empty stratum has share zero.

What de Chaisemartin and D'Haultfoeuille's DiD textbook found

The four applications of the authors' difference-in-differences textbook
(SSC package `cc_xd_didtextbook`) were run in Stata 18 and again in
StatsPAI. Six results were wrong (⚠️ below), the authors' own diagnostic
was missing, and their main estimator refused the book's central example.
Details, one documented departure from the reference command and what is
still open are in `docs/dev/2026-10-05-dcdh-did-textbook-review.md`.

- ⚠️ **`sp.did_multiplegt` with a non-binary treatment.** Every switcher
from a baseline other than 0 was signed as a switch off, whichever way
it moved, and nothing was divided by the size of the change. On the
newspaper panel of Gentzkow, Shapiro and Sinkinson (2011) the estimate
was -0.00082. It is now the effect per unit of treatment, 0.0057791,
which is what Stata `did_multiplegt_old` returns (1e-8), as are the two
placebos. Binary treatments give the same numbers as before.
- ⚠️ **`sp.sun_abraham` on an unbalanced panel.** The weight of a cohort
at a relative time was its share of units. It is now its share of the
observations at that relative time (of their weight under `weights=`),
as in Sun and Abraham (2021) and Stata `eventstudyinteract`. On the
Wolfers (2006) divorce panel, where 52 state-years have no outcome, the
unweighted effect at relative time 0 was -0.1399 and is -0.1054, Stata's
number to nine digits. Weights that vary over time now give Stata's
estimates and standard errors as well. Balanced panels are unchanged.
- ⚠️ **`sp.sun_abraham(event_window=)`.** Relative times outside the
window were left out of the regression, which puts their treated
observations in the reference category. On the Wolfers panel a window of
(-14, 12) moved the effect at relative time 0 from 0.246 to 0.502. The
window now selects what is reported and the regression stays saturated
(`window_rule='report'`, the new default). `window_rule='bin'` pools the
relative times beyond each end into that end, the specification usually
written for `eventstudyinteract`; `window_rule='reference'` is the old
behaviour.
- ⚠️ **`sp.did_multiplegt_dyn(same_switchers=True)`** kept every switcher
observed at all the horizons, including those with no control left at
the long ones, so the switchers were not the same across horizons. It
now keeps the switchers whose effect is estimable at every horizon, as
Stata does.
- ⚠️ **`sp.did_multiplegt_dyn(aggregation='switchers')`** divides the
summed effects by the treatment changes in place at each horizon, which
is the definition of `Av_tot_eff`. Nothing changes for a binary
treatment that stays switched; a treatment that comes back within the
horizons gave a different number before.
- ⚠️ **`sp.did_multiplegt_dyn(controls=)`.** The first-difference
regression behind the option was fitted without the groups that never
switch, whose outcomes were then left unadjusted, and on a panel with
holes a change between two rows several periods apart was used as a
one-period change. The regression also ignored `weights=` and
`trends_nonparam=`. And the analytic variance treated the fitted slopes
as known; it now has the term for their estimation (it was 0.1 to 0.6%
off). Estimates and standard errors agree with Stata to 5e-7 in six
configurations: plain, weighted and clustered, within
`trends_nonparam` cells, normalized, one switch direction, and on an
unbalanced panel.
- ⚠️ **`regress y x, vce(hc2 clustvar, dfadjust)` through `sp.stata` /
`sp.from_stata`** came back as `sp.regress(robust='hc2')`: the cluster
variable and the degrees-of-freedom adjustment were dropped without a
note. It is now `vce='cr2', cluster=, dfadjust=True`.
- **`sp.twowayfeweights`** (new). The weights a two-way fixed effects
regression puts on the treatment effects of the treated cells (de
Chaisemartin and D'Haultfoeuille 2020), for the fixed effects regression
(`type='feTR'`) and the first-difference regression (`'fdTR'`), with
controls, with other treatments in the regression (their contamination
weights are reported too), with observation weights, and with the
regression of chosen variables on the weights. Any treatment, binary or
not, staggered or not. Every count, sum, coefficient and test of the
nine calls in the book agrees with Stata `twowayfeweights`.
- **`sp.did_multiplegt_dyn` takes a non-binary treatment.** Periods are
ranked, so a four-yearly panel works; cells where a group has been on
both sides of its starting treatment are dropped; switchers are matched
on their period-one level. Effects, placebos, analytic standard errors,
normalized effects and `Av_tot_eff` agree with Stata `did_multiplegt_dyn`
on the synthetic test panel and, once one county the command discards by
accident is set aside, on the newspaper panel.
- **Joint tests without the bootstrap.** With `se_method='analytic'`,
`joint_placebo_test`, the new `joint_effects_test`, `joint_overall_test`
and `effects_equal_test` are Wald tests on the analytic joint covariance
and reproduce Stata's p-values. No replicate is drawn any more, so an
analytic fit no longer spends its time on 500 bootstrap draws it did not
use for the standard errors.
- **`design=`, `by_path=` and `normalized_weights=`** in
`sp.did_multiplegt_dyn`: the treatment paths behind an effect, the
effects path by path, and the weights of a normalized effect on the
current and lagged treatments.
- **`sp.regress(dfadjust=True)`**: Bell-McCaffrey degrees of freedom for
`robust='hc2'` and for `vce='cr2'` with `cluster=`, as Stata 18
`vce(hc2 [clustvar], dfadjust)`. `sp.test` now works after
`vce='cr2'` / `'cr3'`, whose full covariance is stored.
- **Translations**: `twowayfeweights`, `did_multiplegt_dyn`, `did_had`
and `did_multiplegt_old` run through `sp.stata`.
- **A sup-t test in `sp.uniform_bands`**, and `terms=`. The result's
`attrs` now carry `max_t` and `supt_pvalue`, the test that every covered
coefficient is zero. It is the test the textbook recommends when many
correlated pre-trend estimates make the Wald test's covariance nearly
singular. `terms=[...]` covers named coefficients of any regression, for
an event study written by hand with lead and lag dummies. On the book's
two examples the p-values are 0.068 and 0.608, against 0.068 and 0.607
from Stata `sotable`.
- **`sp.fect` on large panels.** The initial two-way fit built one dummy
column per unit, so a panel of 7,248 units and 40 periods needed a 17 GB
matrix and did not finish. Above 50 million design cells the same
least-squares fit is solved by alternating projections (equal to the
dense fit to 1e-11 where both run); that panel now takes 0.4 seconds
for `method='fe'`. Smaller panels take the same path as before.
- **`sp.fect` warns when the iterations stop at `max_iter`**
(`ConvergenceWarning`). `model_info['converged']` already said so, and
nothing else did. With weak factors the ATT can keep moving for
thousands of iterations: on the Moser and Voena (2012) panel the
two-factor estimate is 0.277 at `tol=1e-3`, 0.131 at `1e-4` and 0.051 at
`1e-5`.

A causal-ML textbook's notebooks, run against their references

The Python notebooks of Chernozhukov, Hansen, Kallus, Spindler and
Syrgkanis, *Applied Causal Inference Powered by ML and AI*
(arXiv:2403.02467), were rerun in StatsPAI next to R `hdm`, `sensemakr`,
`dagitty`, `rdrobust` and Python `DoubleML`. Double lasso, the DML models,
IRM, the ATT score for difference-in-differences, PLIV, proxy controls and
`rdrobust` with covariates reproduced their references. What did not is
below. Every fix is pinned on simulated data; the book's data are not
redistributed. Reader guide: `docs/guides/applied_causal_ml.md`.

⚠️ Correctness

- **`sp.dml(model='plr')` and `model='pliv'` with a classifier nuisance
used the predicted label, not the predicted probability.** Passing
`ml_m=RandomForestClassifier()` for a 0/1 treatment residualised it as
`D - 1{p > 0.5}`. On the book's 401(k) example the estimate was 6,044
where `DoubleML` and the book give 8,754. A classifier now contributes
`predict_proba` for a 0/1 outcome, treatment or instrument, and is
refused for any other target. Regressors, `model='irm'` and
`model='iivm'` are unchanged. `sp.dml_model_averaging` had the same
defect and the same fix. See MIGRATION `oct2026-causalml-textbook-fixes`.
- **`sp.dml_sensitivity` on an IRM fit used the PLR scaling.** The bound
for the average effect needs `sigma^2 = E[(Y - g(D, X))^2]` and the
Riesz representer `D/m - (1-D)/(1-m)`. The old bound on the 401(k)
example was [-1,276, 16,870]; it is now [2,940, 12,654], equal to
`DoubleML`'s to 1e-10, for the ATE and the ATTE scores.
- **`sp.dml_sensitivity`'s `rv_qa` now accounts for the sampling error of
the bias bound.** It is the confounding strength at which the interval
for the bound reaches zero. It used to subtract `1.96 * se` of the
unadjusted estimate, which ignores that `S` is estimated. On the 401(k)
PLR fit `rv_qa` moves from 0.0486 to 0.0511.

Added

- **`sp.lm_lin`**: the average effect of a randomized experiment by
regression on treatment, centred covariates and their interactions (Lin
2013). Estimate, standard error, degrees of freedom and interval equal R
`estimatr::lm_lin` 2.0.0 to 1e-9 for HC0 to HC3, classical, CR2 (with
Bell-McCaffrey degrees of freedom) and Stata clustering.
`superpopulation=True` adds the variance from estimating the covariate
means, which the book's experiments chapter adds by hand: with random
covariates and heterogeneous effects the regression interval covered
the population effect 83% of the time in simulation and the corrected
one 94%.
- **Standard errors and a uniform band for the TOC curve of
`sp.cate_eval`.** `toc_curve` gains `se`, `ci_lower` / `ci_upper` and
`band_lower` / `band_upper`; the band covers the whole curve, so a
stretch of it above zero is evidence that the score found heterogeneity.
The influence function includes the term for ranking within the
evaluation sample; without it the standard error is more than 10% too
large at the top of the curve. Checked against a known curve
(`tests/reference_parity/test_cate_eval_toc_bands.py`). `m_hat` is no
longer required when `mu1_hat` and `mu0_hat` are given. AUTOC and Qini
are unchanged (equal to grf's to 1e-10).
- **Standard errors and an interval for the sensitivity bounds.**
`sp.dml_sensitivity(...)` returns `se_low`, `se_high`, `ci_low` and
`ci_high` for PLR and IRM fits, with clusters and with repeated
cross-fitting. `DoubleML` 0.11.3 reports the same two standard errors
with the bounds exchanged; a finite-difference check of the influence
functions (`tests/reference_parity/test_dml_sensitivity_bound_scores.py`)
and the book's own code agree with the assignment used here.
- **A weak-instrument-robust confidence set for PLIV.**
`sp.dml(model='pliv')` stores `model_info['anderson_rubin']`: the exact
solution of the Anderson-Rubin inequality on the cross-fitted score. It
is an interval when the first stage is significant and unbounded
otherwise. The book scans a grid, which truncates the set (its AJR
example prints an upper end of 1.99; the set extends to 3.67).
- **`sp.best_linear_projection` accepts `sp.dml(model='irm')` results.** It
projects the cross-fitted doubly-robust scores on covariates or group
indicators, which is `DoubleML`'s `cate()` / `gate()`. Coefficients match
to 1e-10; the joint covariance is in `.attrs['vcov']`, for causal
forests as well, so the projected effect at any covariate value has a
standard error.
- **Joint inference for `sp.rlasso_effects`.** The return value is still a
`dict`, now with `.vcov()` (equal to hdm's to 1e-14) and
`.conf_int(joint=True)`, a sup-t band as in
`confint(<rlassoEffects>, joint = TRUE)`.
- **`sp.sensemakr`** takes string and categorical controls (they enter as
indicators and are benchmarked as a group), takes `kd=` / `ky=`
multiples, and reports the adjusted estimate, standard error, t and
interval for each benchmark, equal to R `sensemakr` 0.1.6 to 1e-9. It
used to stop with a NumPy casting error on the package's own Darfur
example. The treatment standard error is now computed from a QR factor
(1e-14 from R on Track A module 22, was 5e-8).
- **`sp.dag`** accepts named unobserved nodes (`latent=[...]` or
`"F [latent]"`): they are never offered as adjustment variables and
`sp.identify` projects them out (`DAG.latent_projection()`).
`DAG.implied_independencies()` lists the testable implications of a
graph and `DAG.test_implications(data)` tests them; both equal
`dagitty`'s output.

Fixed

- `sp.proximal_regression` refuses a treatment that is not 0/1. It used to
truncate a continuous treatment to integers and return a number.
- `sp.dml_sensitivity` refuses PLIV and IIVM fits with a reason instead of
asking for residuals that those models never store.

What Clarke's *Applied Microeconometrics* found

The Stata chapters of the book's companion site (randomisation inference
and the bootstrap, matching and weighting, clustered inference and the
two-way fixed effects decompositions, synthetic control and synthetic
difference in differences) were run in Stata 18 and replayed through
`sp.stata`. Method and ledger in
`docs/dev/2026-10-04-clarke-applied-microeconometrics-review.md`.

⚠️ Correctness

- **`sp.twfe_decomposition` now decomposes the TWFE coefficient.** Its
2×2 rows were weighted by unit counts and its timing comparisons were
not restricted to their windows, so the rows did not sum to the
coefficient, and its "dCDH" weights were not de Chaisemartin and
D'Haultfoeuille's. On the book's three-unit panel, where every piece is
known by hand, the headline was 1.839 for a coefficient of 27/11. The
rows are now those of `sp.bacon_decomposition` and sum to the
coefficient exactly, the weights are the dCDH (2020) weights, the
estimate is the TWFE coefficient and the standard error is its
unit-clustered one (equal to `xtreg, fe vce(cluster)` to 1e-10). The
comparison with R `bacondecomp` and `twowayfeweights` on `mpdta`, a
strict `xfail` until now, passes. See MIGRATION.
- **Propensity-score matching: units with the same covariates now tie.**
The fitted score of two units with identical covariates could differ in
the last place, because a matrix product rounds each row on its own.
With every tied control kept (`ties='all'`, what `teffects psmatch`
does) that rounding decided which of several identical controls was
matched, and the estimate changed when the rows were shuffled. The
index is now evaluated once per distinct covariate row. On the Lalonde
data with discrete covariates `teffects psmatch, atet` gives 1209.718
(1261.376); StatsPAI gave 1199.666 (1266.607) and now gives 1209.718
(1261.376). Samples with continuous covariates are unaffected. See
MIGRATION.
- **`sp.stata` / `sp.from_stata`: three translations ran something other
than the Stata command.**
- `teffects ipw` was translated to `sp.ipw` with its default bootstrap
standard error. It now asks for the stacked M-estimation standard
error `teffects ipw` reports (`se_method='sandwich'`).
- `teffects aipw` was translated to `sp.aipw` with its default
five-fold cross-fitting, a different (and random) estimator. It now
runs `cross_fit=False, se_method='sandwich'`, and `atet`, which
`teffects aipw` does not have, is refused.
- `boottest` was translated to a call of `sp.wild_cluster_bootstrap`
with arguments that function does not take, so it raised a
`TypeError`. It now runs `sp.wild_cluster_boot` on the regression
before it, with `weight()`, `reps()`, `seed()`, `level()` and a null
value (`boottest x = c`). With few clusters the Rademacher draws are
enumerated and the p-value equals Stata's (0.09179688 on the book's
regression). A `bootcluster()` other than the error cluster, a joint
hypothesis and an unknown weight are refused.

Added

- **A probit propensity score.** `sp.ipw`, `sp.match` and `sp.psmatch2`
take `ps_model='probit'` (default `'logit'`, as before). The sandwich
standard error of `sp.ipw` and the Abadie-Imbens (2016) standard error
of `sp.match` use the probit's score, density and observed Hessian.
Against Stata 18 on the Lalonde data, to the seven digits `teffects`
prints: `teffects ipw (..., probit)` ATET and ATE with their standard
errors, `teffects psmatch (..., probit), atet` with its standard error,
and `psmatch2` (whose default is a probit) to 1e-7.
`tests/reference_parity/test_teffects_probit_stata.py` holds the
numbers. `teffects ... , probit` and `psmatch2` without `logit` are
translated accordingly; they were refused or run as a logit.
- `sp.wild_cluster_boot(h0=)`: the value of the coefficient under the
null. `confidence_set=True` adds `ci_inverted`, the set of null values
the test does not reject, which is what `boottest` prints; the
translated `boottest` asks for it.
- **The Abadie-Imbens (2016) standard error for the ATE.**
`sp.match(estimand='ATE', se_method='abadie_imbens_2016')` reports what
`teffects psmatch, ate` reports, to the seven digits it prints on five
reference fits. The Stata manual prints the score adjustment with a
plus sign; `e(V)` and the paper have a minus (estimating the score can
only lower the variance of the ATE), and that is what is implemented.
`teffects psmatch, ate` is translated with it; it used to be flagged as
not carrying Stata's standard error.
- **`egen` in `sp.stata`**: `count`, `mean`, `median`, `sd`, `min`, `max`,
`total`, `pctile`, `iqr`, `std`, `group`, `tag` and the `row*`
functions, with `if`, `in`, `by()` and behind `by g:`. Each follows the
missing-value rule of `[D] egen`. Run beside Stata 18 on 30 commands:
27 agree to 1e-15 on every row, the other three depend on a row order
Stata's sort does not fix.
- **`sp.stata` runs more of an ordinary do-file.** `quietly { ... }`,
`capture { ... }` and `noisily { ... }` blocks; `capture cmd`, which
swallows what Stata would raise (a variable that is not there,
`restore` without `preserve`) and nothing else; `by g:` / `bysort g
(t):` in front of `generate` and `replace`, with `_n`, `_N` and
subscripts counted within the group; `count [if] [in]`; `_pctile` and
the percentiles, skewness and kurtosis `summarize, detail` leaves in
`r()`; `scalar(name)`; `_b["x"]`; `duplicates drop [varlist, force]`;
`bsample [, cluster()]` (the data are then marked as random);
`predict` after a regression with `i.` variables; the variables
`psmatch2` leaves behind (`_pscore`, `_weight`, `_support` ...); a
`scalar` defined with no data in memory. Loops, frames and matrices
are still refused.
- `scripts/stata_log_replay.py` reads `import delimited` (Stata's
delimiter detection and lower-cased names) and `input` blocks, and
compares `sdid`, `boottest` and `psmatch2` output. It reads a block as
one command (the body from the log's numbered echo, or from the do-file
beside the log when the block ran quietly) and keeps macros, scalars,
matrices and saved datasets across `use`.
- **`sp.stata` runs loops.** `forvalues`, `foreach`, `while`, `if { } else
{ }` and `continue`, with the macros they use (`local x = exp`, `local
++i`, `` `=exp' ``, `` `r(N)' ``, `tempvar`, `tempfile`) and programs
called with arguments. A block's body is run line by line, so each
command in it is translated and refused as it would be outside. Extended
macro functions, `syntax` and `mata` are still refused.
- **`sp.stata` holds more than one dataset.** `save` keeps a copy in the
session (no file is written); `use`, `append using`, `merge 1:1 | m:1 |
1:m ... using` work on those copies and on the
frames passed in the new `files=` argument. Frames (`frame create / copy
/ change / put / post`, `frlink`, `frget`) are run. `merge`, `append`,
and `frlink` / `frget` were run beside Stata 18 and give its
rows in its order.
- **Matrices in `sp.stata`**: `matrix A = J(r, c, .)`, `e(b)`, `e(V)`, a
cell at a time, `A[i, j]` in expressions, `svmat`, `mkmat`, and matrix
expressions with sums, products, a transpose and `inv()`.
- More of `egen`: `skew`, `kurt`, `mad`, `mdev`, `mode`, `rank`, `seq`,
`anycount`, `anymatch`, `cut(, at())`, `rowfirst`, `rowlast`. Nineteen
commands beside Stata 18, all equal on every row.
- `sp.match(estimand='ATE', se_method='abadie_imbens_2016', caliper=)`:
as in Stata the caliper is a condition every unit has to meet, and the
call stops when one does not.
- `sp.causal_forest` accepts `treat=` and `covariates=` for `d=` and `x=`.
- `display "text" exp` evaluates the expression.

With these the three chapter logs of the book run with one command
declined (`regress, vce(bootstrap)`), down from 107.

⚠️ Correctness: local randomization and multi-cutoff RD

Found by running the replication files of Cattaneo, Idrobo and Titiunik
(2024), *A Practical Introduction to Regression Discontinuity Designs:
Extensions*, through StatsPAI and comparing every number with R
`rdlocrand` 2.0, `rdmulti` 2.0 and `rddensity` 2.6. `sp.rdrobust` (20
calls, including clustered, fuzzy and mass-point cases) and the
`sp.rddensity` test agreed to 1e-11 and are unchanged. See MIGRATION.md
for each item below.

- **`sp.rdrandinf(p=1)` estimated the wrong quantity.** It residualized
the outcome on one polynomial pooled across the cutoff. The score is
collinear with treatment inside the window, so the trend absorbed the
effect: 2.90 (p = 0.107) on the Senate data in a ±2.5 window, where the
difference in intercepts is 13.26. It now fits each side separately and
matches `rdlocrand` to 1e-13 in the statistic and the large-sample
p-value, for every kernel, order and evaluation point tried.
- **`sp.rdrandinf(fuzzy=)` reported a p-value that meant nothing.** It
permuted the instrument and recomputed the Wald ratio, so most draws
divided by a first stage near zero. On the book's fuzzy example it
returned p = 0.863 for an effect whose 2SLS p-value is 0.038. The
default is now the randomization test on the reduced form
(`fuzzy_stat='itt'`); `fuzzy_stat='tsls'` gives the large-sample test.
- **`sp.rdrandinf` / `sp.rdrbounds`: `wl` and `wr` are the window's
endpoints**, as in `rdlocrand`, not offsets from the cutoff. The two
agree only at `c = 0`; the book's placebo-cutoff example (`c=1,
wl=0.2348, wr=1.7652`) raised an error. A window that does not contain
the cutoff is now refused with a message that says why.
- **`sp.rdwinselect` did not run the window-selection procedure.** Its
default sequence started at a fraction of the score's range (5, 15.6,
..., 100 points on the Senate data) instead of at the smallest window
with ten observations a side, had no binomial test, did not say which
covariate was least balanced, and without covariates tested quantile
dummies of the score and returned p = 1. It now reproduces Snippet 2.5
of the book: the ten windows, their counts, the binomial p-values and
the recommended window [-0.7652, 0.7652].
- **`sp.rdmc(cutoff_var=)` reported conventional intervals.** The point
estimates matched `rdmulti::rdmc`; the interval and p-value were the
conventional ones, which ignore smoothing bias. On the book's example
the third cutoff went from [0.017, 0.401], p = 0.033 to the robust
[-0.042, 0.408], p = 0.112. The estimate on the normalized score
(`rdmulti`'s "Pooled" row) is now reported as `result.normalized`.

- **`sp.rdrobust`'s automatic manipulation check is now `sp.rddensity`.**
It was McCrary's binned test, which is not built for a score with mass
points. On the book's academic-probation data (429 GPA values) it warned
of sorting with p = 6e-11 where `sp.rddensity` gives 0.082, the number
the book reports; over 60 placebo cutoffs on that data, where nothing
is manipulated, it rejected 43% of the time against 5%. Under simulated
sorting on a continuous score both tests reject every time. Estimates,
standard errors and intervals are unchanged; what changes is when the
warning fires and the value in `model_info['mccrary']['pvalue']` (the
key keeps its name; `'test'` beside it says `'rddensity'`).
`sp.mccrary_test` itself is unchanged.

Added: local randomization options

- `sp.rdrandinf`: `kernel=` (it was accepted and ignored), `nulltau=`,
`bernoulli=`, `ci=` (a grid of effects for the test-inversion
interval), `d=` / `dscale=` (large-sample power, in
`model_info['power']`), `evall=` / `evalr=`, `fuzzy_stat=`.
- `sp.rdwinselect`: `obsmin=`, `wobs=`, `wasymmetric=`, `approx=`,
`n_perms=`, `kernel=`, `dropmissing=`; columns `variable` and
`binom_pvalue`; `attrs['recommended_window']`.
- `sp.rddensity`: `bino_w=`, `bino_wstep=`, `bino_n=`, `bino_nw=`,
`bino_p=` for the binomial tests.
- `sp.rdmcplot`: the RD plot for a design with several cutoffs, one
binned scatter and polynomial fit per cutoff on a shared axis
(`rdmulti::rdmcplot`). The numbers are `sp.rdplot`'s on each cutoff's
units and are returned in `fig.rdmcplot_data`.
- `sp.bitest`: the exact binomial test of Stata's `bitest` / `bitesti`
and R's `binom.test`, from a 0/1 column or from counts. Agrees with
Stata 18 to 1e-12, including `r(k_opp)`. `sp.stata` runs both commands
and stores their `r()` results.
- `sp.rdrandinf(interfci=)`: Rosenbaum's (2007) confidence interval under
interference between units, in `model_info['interf_ci']`.
- `sp.rdwinselect(statistic='hotelling')`: one joint balance test of all
covariates per window (Hotelling's T-squared), with the F p-value under
`approx=True`; agrees with `rdlocrand` to 1e-9.
- `sp.rdwinselect(wmasspoints=True)`: for a discrete score, window `k`
runs from the `k`-th support point below the cutoff to the `k`-th at or
above it, as in `rdlocrand` 1.0 (same windows, counts and p-values).
Releases 1.1 and 2.0 pair the `k`-th on the right with the `(k-1)`-th
on the left, so their first window is empty below the cutoff; on the
book's academic-probation data the first window here is the one the
book itself analyses (208 and 67 observations).
- `sp.from_stata` / `sp.stata`: `rdms` (the boundary points are read from
the `cvar()` variables by `sp.stata`), `bitest`, `bitesti`, and the
`interfci()` and `wmasspoints` options.
- `sp.rdms` without `x2`: one score with cumulative cutoffs
(`rdmulti::rdms(Y, X, C)`), a sharp RD at each cutoff on the whole
sample or, with `ranges=`, on the units in each cutoff's interval.
Agrees with `rdmulti` to 1e-9. `sp.stata` runs `rdms y x, cvar(c)
range(lo hi)` and `rdmcplot y x, cvar(c) pvar() nbinsvar() ...`,
reading the per-cutoff values from the named variables.
- `sp.rdms`: lists for `cutoff1=` / `cutoff2=` estimate several boundary
points in one call and return one row per point; `xnorm=` adds the
pooled estimate on the distance to the boundary. Both agree with
`rdmulti::rdms` to 1e-9 (`test_rdms_points_parity.py`). A call with
scalar cutoffs returns what it did before.

Where StatsPAI departs from `rdlocrand`, and why

- **Randomization p-values with `p > 0`.** `rdlocrand` re-randomizes
treatment labels with the scores held fixed. The observed statistic
extrapolates each side's fit to the cutoff; a relabelled sample fits
through the interior, with about a quarter of the variance. On data
with no effect, 40 observations and `p = 1`, that test rejects a true
null 37% of the time at the 5% level. `sp.rdrandinf` permutes outcomes
against (score, assignment) pairs and rejects 4.7% (3,000 replications)
(`tests/reference_parity/test_rdlocrand_extensions_parity.py`). For
`p = 0`, with or without a kernel, the two schemes coincide and the
60-seed mean p-values agree with R's.
- **The default windows of `sp.rdwinselect`** are those of `rdlocrand`

regression-adjustment ATE (Stata `teffects ra` 1.897784556341). It does not
extrapolate the Angrist–Rokkanen CATE(x) = E[Y1 − Y0 | X = x]. The
`sp.margins_at` / `sp.contrast` / `sp.pwcompare` p-values and intervals use
N(0, 1) instead of Stata's t(df_r) (G − 1 under clustering), `pwcompare`'s
Sidak adjustment underflows to 0 below p ≈ 1e-17, `sp.margins(method="mem")`
evaluates factor dummies at the base level, an `at=` variable not in the
model is ignored, and the margins family averages over all of `data=`
rather than the estimation sample; strict xfails pin Stata's numbers.
Estimates and SEs of those three functions are already bit-exact with
Stata. Not fixed in this release: `src/statspai/postestimation/` is being reworked on a separate line, which the xfails will flag when it lands.

Reference-implementation findings (round 2)

- **Stata `sbcusum.ado`** hard-codes CUSUM boundary constants whose
crossing probabilities are 0.0100000791, 0.0499991147 and 0.0999993370
(relative errors −7.7e-7, +2.5e-6, +1.2e-6). StatsPAI uses the exact root
of strucchange's closed form, which an independent density-propagation
computation confirms to about 1e-12.
- **Stata `_xturfisher.ado`** codes the Fisher L* scale as `3*(5N+3)` while
its df is 5N+4; Stata's own `[XT] xtunitroot` manual prints 5N+4, as do
plm and the logistic-moment identity (asserted for N ∈ {1, 2, 5, 12, 50}).
- **Stata `estat sbsingle`** uses a different trimmed range from strucchange
when 0.15·n is not an integer (n = 250: first regimes 38..212 against
37..213); the p-values agree on the fixtures.
- **Stata `chks`** drops rows with y ≤ 0 even in the linear model (it runs
`gen ln(depvar)` and then `drop if rowmiss`): 90 of 500 rows on the
fixture. The comparison runs on y + 10.
- **R `metafrontier` 0.3.1** (SFA branch): `EC_group` is Shephard-oriented
while `TC_group` is Farrell-oriented, so `MPI_group` mixes the two.
- **R `sfa::zsfm` 1.2.0**: L-BFGS-B stops with a score of 3e-4 to 7e-4;
the T2 comparison is on its own log-likelihood polished to its optimum.
- **Stata `xsmle` 1.4.5**, two-way SAR / SDM: its reported maximum is
below StatsPAI's by 4.85e-5 (SAR) and 5.15e-4 (SDM) on its own objective;
`type(both, leeyu)` returns the `type(ind, leeyu)` fit (N = 768 = 48·16),
so the time effects are dropped.
- **R `rd2d` 1.0.0**, `kernel_type="rad"`: it fits at radius
sqrt(hx² + hy²) but rescales the variance with hx·hy, inflating every
radial SE by sqrt((hx² + hy²)/(hx·hy)) (√2 when hx = hy). Weighted `lm()`
plus `sandwich::vcovHC(HC0)` at that radius equals StatsPAI's SE to 1e-10.
- **R `RDHonest::RDHonestBME`** adds `maximum.bias` (outcome units) to a z
statistic in its p-value; StatsPAI reports the p-value with bias/SE and
keeps R's as `p_value_rdhonest`.
- **Stata `paramed`**, no covariates, with interaction: the NIE gradient
uses θ3·a1 where θ2 + θ3·a1 belongs (the covariate branch is right). Its
SE 0.055580 is reproduced exactly from that gradient; the correct SE is
0.096576.
- **R `mediation::medsens`** stops its default iteration about 1e-4 short
of the fixed point at |ρ| = 0.9; `eps=sqrt(machine eps)` reproduces it.
- **GRAPPLE** (GitHub 317e837) stops its `optim` / `uniroot` early (score
O(1e-3) at its estimate): fitted β agrees to 8e-4, τ² to 6e-5, SE to
3e-5; its sandwich at its own point matches to 8e-16 (l2).
- **R `scinference`** calls `limSolve::lsei(type=1)` without checking
`IsError`. At θ0 ∈ {8.75, 9} lsei fails and returns infeasible weights
(SSR 310.4 against 86.6 at the simplex optimum), and scinference reports
joint p 0.55 / 0.65 where the value is 0.2 / 0.2. StatsPAI equals
scinference's code with quadprog as the solver on every grid point.
- **Stata `synth.ado`** rounds the unit weights to 3 decimals
(`roundmat`), so no Stata SCM tool can reach T2; not run as a reference.

Defaults aligned with the reference implementation

- **`sp.roc_curve(se_method=)` defaults to `"delong"`** (was `"hanley"`, an
approximation no reference uses by default). DeLong is the default of R
`pROC` and Stata `roctab`; `se_method="hanley"` restores the old SE.
- **`sp.direct_standardize(ci_method=)` defaults to `"gamma"`** (was
`"lognormal"`, which no package computes). The Fay-Feuer gamma interval is
the default of R `epitools::ageadjust.direct`; `ci_method="lognormal"`
restores the old interval.
- **`sp.power_case_control(test=)` defaults to `"chi2"`** (was `"wald"`). The
Pearson chi-squared test is what Stata `power twoproportions` and R
`power.prop.test` compute; `test="wald"` restores the old power (the
docstring example moves from 0.9213 to 0.9171, and the required cases from
138 to 141).
- `sp.svydesign(lonely_psu=)` keeps its default: a single-PSU stratum
contributes zero variance (Stata `singleunit(certainty)`) and now warns
unless the rule is given explicitly; R `survey` fails by default and Stata
reports a missing SE. Pass `lonely_psu="fail"` for R's behaviour.

Provenance

- `sp.cuminc` / `sp.finegray`: the module docstring of
`statspai.survival.competing_risks` now records that Gray's CIF variance,
the covariance of Gray's test and the `crr` sandwich variance were
reimplemented from Gray (1988) and Fine & Gray (1999) as NumPy, with
conventions fixed against R `cmprsk` 2.2-12 (Bob Gray, GPL-2 | GPL-3); no
`cmprsk` code is copied, and `cmprsk` is used only to generate test
reference values.

Added

- **Opt-in IRM ATE OOF audit records.** `sp.dml`, `sp.DoubleML.fit`, and
`sp.DoubleMLIRM.fit` accept `store_oof`, `observation_ids`, and in-memory
`external_predictions`. `sp.OOFPredictions` validates row/fold alignment and
declared training scopes; `sp.OOFBundle` exports all cross-fit repeats with
versioned, hash-checked JSON. `result.get_oof()` and `get_residuals()` return
defensive snapshots; `bundle.score_concentration()` reports squared-score
concentration with explicit unavailable states. The supported audit path is
unweighted, unnormalised binary IRM ATE. Normal fits retain their numerical
behavior, and ordinary result serialization excludes individual OOF arrays.
See `docs/dml_oof_audit.md` for scope and provenance limitations.
- **DML diagnostic semantics correction.** The former orthogonality p-value
tested a score centered by construction and was uninformative. Its statistic
and p-value are now unavailable (`None`) with an explicit explanation;
IRM plots label the displayed series as a centered IRM score residual.

- **One Stata `vce()` grammar on every SE-bearing regression estimator.**
`robust=` / `vce=` accept `True`, `"robust"`, `"vce(robust)"`, `"r"`,
`"oim"`, `"hc0"`–`"hc3"`, and a cluster variable written inline —
`vce="cluster firm"`, `"vce(cluster firm)"`, `"cl firm"` — on `regress`,
`ivreg` / `iv`, `glm`, `logit`, `probit`, `cloglog`, `poisson`, `nbreg`,
`ologit`, `oprobit`, `mlogit`, `clogit`, `zip_model`, `zinb`, `hurdle`,
`truncreg`, `biprobit`, `betareg`, `fracreg`, `etregress`, `liml`,
`ppmlhdfe`, `panel_logit`, `panel_probit` and `subgroup_analysis`. One
parser (`statspai.core._vcov_spec.parse_se_request`) resolves every
spelling; a spelling an estimator does not implement raises
`MethodIncompatibility` instead of falling back to other standard errors.
`vce="robust"` on `sp.regress` (previously an error) is HC1, as in Stata.
- **House-style keyword spellings on the causal entry points.** The panel
identifier is accepted as `id=`, covariates as `covariates=`, and adoption
cohorts as `first_treat=`, alongside each estimator's native spelling:
`callaway_santanna` / `staggered_cs` / `staggered_sa` / `sun_abraham`
(`id`/`unit` -> `i`, `time` -> `t`, `first_treat`/`cohort` -> `g`),
`bjs` / `gardner_did` / `wooldridge_did` / `etwfe` / `stacked_did` /
`did_multiplegt(_dyn)` (`id`/`unit` -> `group`, `covariates` ->
`controls`, `treat` -> `treatment`), `event_study` / `synth` / `sdid` /
`gsynth` / `augsynth` / `fect` (`id` -> `unit`), the RD family
(`running` -> `x`, `cutoff` -> `c`, `covariates` -> `covs`; `x` stays the
running variable), `psm` (`treat` -> `d`, `covariates` -> `X`), and
`controls=` on `did`, `ddd`, `did_2x2`, `drdid`, `match`, `psmatch2`,
`ipw`, `aipw`, `tmle`, `ebalance`, `metalearner`, `dml`. Signatures and
schemas are unchanged; passing both spellings raises `TypeError`, and a
misspelt keyword names the closest parameter (`runing=` -> "did you mean
'running'"). Alias maps are listed by `sp.describe_function(name)["aliases"]`.
- **`sp.stata("...", data=df)`** runs Stata command lines through
`sp.from_stata` and returns the fitted result: several lines (or `;`),
`*` / `//` comments, and `test` / `lincom` / `margins, dydx(...)` applied
to the latest fit. A line whose translation is incomplete — `xtreg, fe`
without the `xtset` panel id — raises instead of running a different model.
- **`sp.list_functions(core=True)`** returns the ~30 everyday estimation and
post-estimation verbs (`sp.registry.CORE_FUNCTIONS`), in curated order.
- `inspect.signature(sp.iv)` / `sp.rd` (and `help`) show the dispatcher's
parameters instead of `(*args, **kwargs)`.
- **`result.test(...)` and `result.lincom(...)`** on `EconometricResults`
and `CausalResult`, mirroring Stata's post-estimation commands.
`sp.test` also accepts `test x1 x2`, chained `x1 = x2 = x3`, grouped
`(x1 = 0) (x2 = 1)`, `2*x1`, `_cons` and `_b[x]`.
- **`sp.glm(..., information="observed" | "expected")`.** Observed
information (Stata's default) is now the variance bread; `"expected"`
reproduces R's `glm`. They coincide for canonical links.
- `nbreg` keeps the MASS::`glm.nb` standard error (conditional on the
dispersion) in `model_info["se_conditional_on_dispersion"]`.
- Stata 18 reference fixtures `tests/reference_parity/_fixtures/
_generate_vce_grammar_stata.do` (29 estimators x oim / robust / cluster)
and `_generate_postestimation_stata.do` (test / lincom / p-values /
margins), pinned by `test_vce_grammar_stata_parity.py` and
`test_postestimation_stata_parity.py` at the 1e-6 default budget (three
documented optimiser-stopping exceptions, each asserted by comparing the
objective at both solutions).

- **`sp.causal_forest(..., fe="twoway", id=..., time=...)`: causal forests
with fixed effects for panel data.** Unit (and period) effects are removed
*within every node and every honest leaf*, not once globally
[kattenberg2023causal]. On a staggered panel where adoption selects on the
unit effect (N = 300 units, T = 6, 4 seeds) the pooled forest's CATT bias
was +0.36 to +2.37 and the FE forest's -0.06 to +0.03 (RMSE 0.32-0.42 with
heterogeneous effects). Unbalanced panels are handled by iterating the
two-way within transformation to 1e-10. Doubly-robust averages are refused
for FE forests (no propensity exists for a within-unit design).
- **`clusters=` for `sp.causal_forest`.** Trees are subsampled, split into
honest halves and cross-fitted by cluster, and every average, calibration
test, BLP and RATE reports cluster-robust standard errors;
`equalize_cluster_weights=` gives clusters equal weight. Repeated
observations of a unit must be clustered.
- **`sp.did_forest`: difference-in-differences causal forests for staggered
adoption.** One honest causal forest per clean Callaway-Sant'Anna
group-time comparison on the long-differenced outcome (the DiD causal
forest of [gavrilova2025difference] applied to every group-time cell),
doubly-robust ATT(g,t), event-study / cohort / overall aggregates with
unit- or cluster-level influence-function standard errors, a joint
pre-trend Wald test, unit-level CATEs with little-bag standard errors and a
per-cell heterogeneity test. Monte Carlo (200 replications, N = 800, T = 7,
two cohorts): see `docs/guides/heterogeneity_panel_forests.md`.
- **`sp.calibrate_cate`**: CATE predictions rescaled by the out-of-bag
best-linear-predictor calibration [chernozhukov2025generic]; warns when no
heterogeneity is detected and, for FE forests, that the slope does not
de-attenuate (measured, see the guide).
- `CausalForest.oob_effect()`, `.effect_variance()`, `.split_frequencies()`,
precomputed `Y_hat=` / `W_hat=` nuisances, and a doubly-robust average for
continuous treatments (`average_treatment_effect(target_sample="all")`, the
Riesz representer `(W - W_hat) / Var(W | X)`). On Card (1995) the return to
a year of schooling is 0.0768 (SE 0.0044).

Changed

- **`sp.metalearner` states what its standard error refers to.**
`result.estimate` / `result.se` are the AIPW average and are numerically
identical across `learner='s'|'t'|'x'|'dr'`; the learner-specific
quantity is `model_info['cate_mean']`, which ships no standard error.
The docstring now says so explicitly and `model_info` carries
`se_refers_to` and `cate_mean_se = None`, so a caller cannot pair the
two by accident. No numerical change.
- **Neural estimators (`sp.tarnet`, `sp.dragonnet`, `sp.cfrnet`) flag their
standard error.** `se` is the bootstrap dispersion of the mean of the
fitted individual effects with the network held fixed
(`se_method = 'unit_bootstrap_cate_plugin'`), which ignores estimation
uncertainty in the network; the ML4CI Monte Carlo finds it at 0.3–0.9 of
the sampling dispersion. `model_info` now carries
`se_valid_for_ate = False` and an `se_note`. No numerical change.

- Printed summaries: coefficient tables label the statistic `z` / `P>|z|`
for likelihood-based fits (the p-values already used the normal);
diagnostics print counts as integers (`N instruments : 1`, not `1.0000`);
the `CausalResult` footer omits `None` fields and backend flags and shows
floats to six significant digits. Stored values are unchanged.
- `sp.etable([r1, r2])` accepts a list; a non-result argument or an empty
call raises `TypeError` (both used to return an empty DataFrame), and a
mixed pyfixest / native call no longer drops the native columns.
- `margins_at`, `contrast` and `pwcompare` refuse models whose prediction
scale is not the linear index (logit, poisson, …) and results without a
coefficient covariance matrix, instead of returning index-scale numbers
or assuming zero covariances.
- `betareg(link="cloglog")` fits the cloglog link (it silently fitted logit)
and `link="loglog"` is available; outcomes on or outside (0, 1) raise
instead of being clipped to [1e-6, 1-1e-6]. `fracreg` rejects outcomes
outside [0, 1]. `truncreg` drops observations outside the truncation
limits with a warning, as Stata does. `xtnbreg(model="re")` raises on
`robust=` / `cluster=` instead of warning and ignoring them.
`panel_logit(method="fe")` raises on `robust=` / `cluster=` (Stata rc 198)
and points to `sp.clogit(..., vce="cluster id")`.

- `sp.causal_forest` defaults follow GRF: `n_estimators=2000` (was 100),
subsampling without replacement (`bootstrap=True` now raises),
`ci_group_size=2` (1 with a warning when `max_samples > 0.5`), and
nuisances from out-of-bag regression forests on `[X, W]` (user sklearn
estimators are cross-fitted in cluster-respecting folds).
- For GRF forests, `calibration_test` is the regression of `Y - Y_hat` on
`(W - W_hat) * mean(tau)` and `(W - W_hat) * (tau - mean(tau))` with HC3
(cluster-robust) standard errors and adds `t_vs_zero` / `p_one_sided`;
`best_linear_projection` regresses AIPW scores on unstandardised
covariates with HC3 standard errors; `average_treatment_effect(
target_sample="overlap")` is the partially linear coefficient. These are
the definitions of `grf::test_calibration`, `grf::best_linear_projection`
and `grf::average_treatment_effect`. In-sample inference functions refuse
X/Y/T other than the training data.

Fixed

- **`sp.causal_forest` accepted a keyword its registry entry did not name
and rejected the one it did.** The registry advertised `n_trees` while
the callable only took `n_estimators`, so a caller following
`sp.describe_function("causal_forest")` got a `TypeError`. The registry
now names `n_estimators` and the function accepts `n_trees` as an alias
(through `accepts_aliases`, so passing both raises `TypeError` and
the alias is listed by `sp.describe_function`). Found by the ML4CI
agent-surface audit; it was the only hard schema-versus-signature
failure among 252 audited callables.

- **`sp.from_stata` post-estimation translations run.** `margins, dydx(x)`
emitted `dydx=`, `test x1 x2` emitted `terms=`, and `contrast x` emitted
`terms=` — keywords `sp.margins`, `sp.test` and `sp.contrast` do not take,
so every copied call raised. They now emit `variables=` (plus
`method="mem"` for `atmeans` and an `at={...}` dict), the restriction text
as `hypothesis=`, and `variable=`; `lincom` is translated (with
`level()` → `alpha`), and `marginsplot` targets `sp.marginsplot`. Predictive
margins (`margins` without `dydx`) and options with no `sp.margins`
counterpart (`eyex`, `over()` ...) are refused instead of being mapped to
average marginal effects. The structural translation sweep exempted
post-estimation tools; a new test dispatches each payload on a fitted model.
- `sp.stata` no longer warns "pipe the previous estimator's result_id" on
`margins` / `test` / `lincom` / `contrast` / `marginsplot` lines: that note
is advice for `sp.from_stata` callers, and `sp.stata` does the piping itself.
- README: a sixth beginner example (clustered logit + `margins` / `test` /
`lincom` / `sp.stata`, checked against Stata 18), the `vce()` grammar and
house-style argument names in the Stata/R table, and refreshed tier counts;
the new outputs are pinned by `tests/test_readme_examples.py`.
- `sp.from_stata("xtreg y x, fe")` without a panel id printed
`sp.feols('y ~ x', ...)` — a runnable pooled regression; the formula now
keeps the `| <panel_id>` fixed-effect placeholder.
- Docstring examples of `EconometricResults` / `CausalResult` methods and of
the registry helpers ran on undefined `sp` / `df` (or showed placeholder
output), so `pytest --doctest-modules` failed on 21 of them; each is now
self-contained with checked output.
- `docs/stats.md` (per-module table, source / test LOC) and the README
at-a-glance LOC figures regenerated from `scripts/registry_stats.py`.

- **`sp.from_stata` IV translations were not runnable.** `ivreg2` /
`ivregress` joined multi-variable lists with spaces
(`y ~ x1 x2 + (d ~ z1 z2)`), which `sp.ivreg` rejects, and refused
exogenous regressors written after the parentheses or several endogenous
regressors; `ivreghdfe` shared the last two gaps. The note about Stata's
`small` option was also inverted: `sp.ivreg` SEs match `ivregress ...,
small`, so the note now appears when `small` is absent.
- **`sp.cross_validate(result)`** now explains that the result does not carry
its data and prints a runnable `sp.cross_validate(df, ...)` call rebuilt from
the result's provenance (the old message pointed at a non-existent `data=`
argument). Data mode accepts `sp.ivreg`'s own formula, `y ~ (d ~ z) + x`.
- **`sp.regtable(...).to_text()`** widens columns to fit their widest cell;
a fixed 14-character column ran long labels into their neighbours.
- `WeakIVConfidenceSet` has a compact `repr` instead of dumping every grid
array.
- `examples/rd_lee.py` and the RD snippet in the review-time documentation
crashed (they used columns `lee_2008_senate()` no longer returns);
`tests/test_example_scripts.py` now runs every example script.
- `sp.datasets.lee_2008_senate()` documentation no longer attributes the
bundled U.S. Senate extract to Lee (2008), which studies House elections;
the data are Cattaneo, Frandsen & Titiunik (2015), distributed with R
`rdrobust`.
- `docs/guides/migration-from-r.md` mapped R post-estimation calls to four
methods that do not exist (`result.robust()`, `result.summary(vcov=)`,
`result.test()`, `result.marginal_effects()`) and to a multiway-cluster
spelling that raises; they now point at `sp.regress(robust=)`, `sp.test`,
`sp.margins` and `vcov={"CRV1": "c1 + c2"}`.

- Citation keys: the best-linear-predictor calibration test was attributed to
`chernozhukov2020generic` (a paper on quantile effect functions) and the
generalized random forest to `athey2019surrogate` (the surrogate index).
They now point to `chernozhukov2025generic` and `athey2019generalized`;
`souto2025forests` lists its second author's family name as recorded by
arXiv and DataCite. New verified entries: `kattenberg2023causal`,
`gavrilova2025difference`, `aytug2026attenuated`, `aytug2026fixed`,
`aytug2026causalfe`, `chernozhukov2025generic`.

⚠️ Correctness

- **`sp.dml(model="plr", score="IV-type")` now computes DoubleML's IV-type
score.** The option fits the third nuisance `g(X) = E[Y − θ̃D | X]`
(cross-fitted on the same partition with a clone of `ml_g`, `θ̃` the
preliminary partialling-out estimate) and solves
`E[(Y − Dθ − g(X))(D − m(X))] = 0`, as `doubleml.DoubleMLPLR` does in both
its ports. Since the option was introduced it had substituted
`l(X) = E[Y|X]` for `g(X)` while keeping the IV-type denominator `Σ v̂·D`.
That variant is unbiased at the truth but is **not Neyman-orthogonal in
the propensity nuisance** — the Gateaux derivative of its moment in the
direction of `m̂` is `θ₀·E[m₀(X)δ(X)]`, so regularisation error in `m̂`
entered the estimate at first order — and it reproduced neither DoubleML
port: on the ML4CI benchmark (n = 2,000, p = 20, five shared folds, OLS
learners) it returned 1.0232 where both ports return 1.0094, a gap of
0.40 standard errors that the unpinned parity test's 0.05 tolerance
could not see. With linear nuisance learners the corrected option
coincides with partialling out to machine precision (the ML4CI companion
paper's Proposition 1); with regularised learners it is a distinct
estimator, and its standard error changes accordingly.
`tests/test_dml_ivtype_reference.py` pins it to DoubleML-for-Python on a
shared fold vector to 1e-10 (linear) and 1e-8 (lasso). Found by the
fold-controlled reconciliation in the ML4CI companion paper. See
MIGRATION.md.

- **`sp.causal_forest(...).ate()` and `.att()` now return the doubly-robust
estimate their own standard error describes.** The float value was the
plug-in average of the fitted CATE predictions while `.se`, `.ci` and
`.pvalue` came from the GRF-style AIPW score, whose point estimate is a
different number: on a misspecified binary design (n = 1,500, 500 trees)
the two differ by 0.049 for the ATE and 0.120 for the ATT against a
standard error of 0.081, so the printed interval was not centred on the
number beside it and a caller who compared the point estimate with a
benchmark was not comparing the quantity the interval covered. The float
is now the doubly-robust estimate — what `grf::average_treatment_effect`
returns — and the plug-in average is kept in
`detail["plug_in_estimate"]`, with the difference in
`detail["plug_in_minus_aipw"]` as a diagnostic (it grows where the
fitted propensity approaches its bounds). A continuous treatment, where
the score cannot be formed and the aggregation falls back to the plug-in
average, is unaffected. The package's own reference-parity test
`test_causal_forest_aipw_recovery.py::test_plugin_is_more_biased_than_aipw`
asserted that the plug-in average is the more biased aggregation while
reading it off `cf.ate()`, which is what made the defect visible. Found
while tracing the causal forest's miss on the LaLonde benchmark in the
ML4CI companion paper. See MIGRATION.md.

- **Rows with a missing cluster variable leave the estimation sample, on
every estimator that takes the Stata `vce()` grammar** (`regress`, `ivreg`,
`liml`, `glm`, `logit`, `probit`, `cloglog`, `poisson`, `nbreg`,
`ppmlhdfe`, `ologit`, `oprobit`, `mlogit`, `clogit`, `zip_model`, `zinb`,
`hurdle`, `truncreg`, `biprobit`, `etregress`, `fracreg`, `betareg`,
`panel_logit`, `panel_probit`), as Stata's `vce(cluster v)` and R
`fixest` do. A `StatsPAIWarning` gives the count and
`model_info["n_missing_cluster_dropped"]` records it. Before, the
estimators disagreed: `regress`, `glm`, `logit` and the count / zero-inflated
models raised, the ordered and multinomial logits dropped the rows, and
`ivreg` / `liml` **kept them in the fit** -- coefficients and clustered
standard errors then differed from Stata's (and from the same call on the
complete rows). One decorator (`core._vcov_spec.markout_clusters`) now does
it for all of them; `tests/test_cluster_markout.py` pins every estimator to
its fit on the complete rows.
- **`sp.margins(fit, data=df)` returned NaN** when `df` had a missing value in
a model variable (for example the rows the fit itself had dropped). It now
raises `MethodIncompatibility` naming the row count; omit `data=` to average
over the fit's own estimation sample (Stata's `e(sample)`). `sp.stata`'s
`margins` line no longer passes the raw data when the fit carries its
design, so it averages over the estimation sample as Stata does.

- **ML robust / cluster standard errors now follow Stata.** `vce="robust"`
on likelihood-based estimators omitted Stata's N/(N-1) factor (SEs
~0.04% small at N = 1,200); `glm`, `ologit`, `oprobit`, `mlogit`,
`clogit`, `zip`, `zinb` and `hurdle` cluster SEs carried the regress
family's (N-1)/(N-K) factor that Stata's ML commands do not apply.
`vce="hc0"` keeps the unscaled R `sandwich` HC0.
- **Options that were accepted and ignored now take effect:** `robust=` /
`cluster=` on `truncreg`, `biprobit`, `betareg`, `panel_logit` /
`panel_probit` (RE), and `robust=` on `subgroup_analysis`. `logit` /
`probit` / `cloglog` answered `vce="cluster g"` with heteroskedasticity-
robust SEs; `poisson` and `nbreg` answered it (and `nbreg` `"hc3"`) with
model-based SEs; `ologit` family answered any unknown string with OIM.
- **`nbreg` standard errors.** The robust / cluster sandwich used the Poisson
score `y - mu` (robust SEs up to ~150% off against Stata on overdispersed
data) and the OIM bread ignored the beta / dispersion cross-information;
`dispersion="constant"` (NB1) coefficients were 8.6% off. Both now come
from the joint (beta, ln dispersion) likelihood, polished by Newton, and
match Stata `nbreg` to ~1e-12 (NB2) / 1e-7 (NB1).
- **`sp.test` / `sp.lincom` used a diagonal covariance.** Any restriction on
two correlated coefficients was wrong (`test x1 = x2` after robust OLS:
F 154.4 vs Stata 167.0); they now use the full covariance and the fit's own
reference distribution (F / t or chi2 / z). A misspelled coefficient used
to be read as zero and now raises.
- **p-values and confidence intervals.** Clustered `regress` / `ivreg` used
t(N-K); they now use t(G-1) as Stata does (a p-value of 7e-32 vs Stata's
9e-15 on 40 clusters). Likelihood-based estimators used t(N-K) and now use
z, as Stata reports.
- **`sp.margins` returned index coefficients.** After `logit` it reported
beta (0.568) where Stata's `margins, dydx(*)` reports the average marginal
effect on Pr(y) (0.126), and with interaction terms its "SE" was
std(dydx)/sqrt(n). It now computes AMEs on the prediction scale (identity,
logit, probit, cloglog, log links) with delta-method SEs from the full
covariance, matching Stata to ~1e-10.
- **LIML robust / cluster meat** uses the instrument-projected regressors
P_Z X, as Stata `ivregress liml` and `linearmodels` do, instead of
(I - kappa M_Z) X (~0.05% SE shift); `sp.liml(robust=...)` now applies the
same N/(N-K) small-sample factor as its classical and cluster SEs and as
`sp.iv(method="liml")`.
- **`clogit` standard errors** came from BFGS's quasi-Newton `hess_inv`
(off by up to 4% against Stata); they now use the analytic information.
- **`glm` with a non-canonical link** uses the observed information (see
Added) and iterates until the coefficients, not only the deviance, have
converged (gaussian/log coefficients were ~1e-4 short of the MLE).
- Newton-polished, complex-step observed information for `truncreg`,
`biprobit`, `betareg`, `fracreg`, `mlogit`, `ologit`, `oprobit`, `zip`,
`zinb`, `hurdle` and panel logit / probit replaces second-difference
Hessians (1e-6 – 3e-5 relative SE error).

- **`sp.ivreg` / `sp.iv` first-stage F now uses the requested variance
estimator.** With `robust=` or `cluster=`, the first-stage F (the
`diagnostics` entry, `model_info['first_stage_f']`, the weak-instrument
warning and `sp.regtable`'s "First-stage F" row) was still the classical
homoskedastic F while every other number used the robust / clustered vcov.
It is now the Wald test of the excluded instruments under that vcov divided
by the number of instruments, on `F(m, N-k)` or `F(m, G-1)` -- what Stata's
`estat firststage` reports after `ivregress ..., vce(robust)` /
`vce(cluster)`. On the Card (1995) extract with HC1 it moves from 16.72 to
17.51. Checked against Stata 18 to ~1e-13 for one and two instruments, HC1
and clustered (`tests/reference_parity/test_iv_first_stage_f_stata_parity.py`).
Non-robust fits are unchanged. The classical F stays available as
`First-stage F, non-robust (<endog>)` and
`model_info['first_stage'][j]['f_statistic_nonrobust']`, since Stock-Yogo
critical values are tabulated for it. `sp.lasso_iv` defaults to HC1, so its
reported first-stage F changes too.

- **`sp.cross_validate` could compare different models or different variance
estimators.** Three defects, each able to produce a wrong verdict:
(1) with an IV formula (`y ~ w1 + w2 | x ~ z`) the first exogenous regressor
was promoted to "treatment" and dropped from the statspai and linearmodels
engines' models, which then agreed with each other against pyfixest;
(2) the statspai engine ignored `cluster=` for OLS / FE / IV and `vcov=` for
IV, and the R engine passed fixest the string `"cluster"`, which clusters on
the first fixed effect (or errors without one);
(3) linearmodels reported large-sample SEs. Every engine now fits the same
model with the same variance estimator and small-sample convention; on
OLS / FE / IV with clustered and HC1 errors the statspai, pyfixest, R and
linearmodels SEs agree to ~1e-10 (linearmodels with absorbed FE differs by
a documented FE-dof convention, ~0.15%).

- **`sp.regtable(..., model_labels=[...])` with a repeated label silently
dropped columns** from `to_dataframe()`, `to_excel()`, `to_word()`,
`save('.csv')`, `to_dict()` / `to_json()` and `sp.collect` exports: each row
was keyed by label, so `["OLS", "OLS", "2SLS"]` exported two columns, the
surviving `OLS` holding the second model. Text / LaTeX / Markdown / HTML were
unaffected. All exports now keep every model; `to_dict()["columns"]`
disambiguates a repeated label as `"OLS (1)"`, `"OLS (2)"` while
`model_labels` keeps the labels as given.

- **Weak-IV confidence sets printed inward-biased endpoints.**
`WeakIVConfidenceSet.summary()` and `.as_intervals()` (Anderson-Rubin, CLR,
K) reported the extreme grid points inside the set, not the root-found
`lower` / `upper`, so the printed interval was too narrow by up to a grid
step (Card 1995 AR set: `[0.0389, 0.2601]` printed, `[0.0384, 0.2612]`
correct). The outer endpoints now come from `lower` / `upper`; interior
boundaries of a disconnected set stay at grid resolution.

- **`sp.causal_forest` now grows generalized random forests.** Since it was
introduced the estimator grew scikit-learn regression trees *on the
treatment residual* -- splits chased variation in the propensity residual,
not treatment-effect heterogeneity -- kept the sklearn leaf mean of that
residual in leaves with fewer than two honest rows, and ran every in-sample
diagnostic (ATE, calibration test, RATE, BLP) on predictions from trees that
had seen the rows. The new default engine implements the GRF algorithm
[athey2019generalized]: gradient-based causal splits with stabilised
child constraints, honest cluster-level subsampling with leaf pruning,
out-of-bag nuisance forests and CATE predictions, and bootstrap-of-little-
bags variances. On five designs with tau(x) = 1 + 2 x0 (n = 2000, 500 trees)
the correlation of the CATE with the truth rose from 0.922 to 0.990, RMSE
fell from 1.02 to 0.31 and the predicted spread went from 59% to 93% of the
true spread, matching R grf 2.6.1 (0.990 / 0.305 / 93%). Against grf on a
shared fixture, RMSE ratios are 0.998 (i.i.d.) and 1.017 (clustered),
median pointwise-variance ratios 1.01 / 0.94 and CI coverage within one
point (tier T3, `tests/reference_parity/test_grf_engine_statistical_parity.py`);
the calibration test, AIPW ATE/ATT/ATC/overlap, BLP and `sandwich::vcovCL`
HC0-HC3 given a forest match grf / sandwich to 8.5e-15. Track A module 13:
ATE SE now within 0.22% of grf (was 7.7%). `split_rule="legacy"`
reproduces the old estimator for one release. See MIGRATION.md.
- **`CausalForest.effect_interval()` returned percentiles of individual tree
predictions**, which are not confidence intervals for the forest's CATE.
It now returns `tau(x) +/- z sqrt(V(x))` with the little-bag variance.
- **`sp.honest_variance` measured no sampling uncertainty**: it re-split one
fixed vector of predictions and divided by sqrt(n_splits), so its SE went
to zero as n_splits grew. It is deprecated and, for GRF forests, returns
the doubly-robust ATE and its influence-function SE.
- **`sp.callaway_santanna(control_group="notyettreated")` admitted
already-treated cohorts as controls** for pre-treatment cells under
`base_period="universal"`, and ignored `anticipation` in the control set:
the cutoff was `G > t` where R `did` uses `G > max(t, base) + anticipation`
(`get_did_cohort_index`). On a three-cohort panel cohort 6's placebo
ATT(6, 1..3) were 0.54-0.58 with covariates (R: -0.03 to -0.02) and
1.21-1.43 without (R: 0.27-0.47). All cells now match R `did` 2.3.0 to
3.6e-12 (ATT) / 8.9e-11 (SE) across panel and repeated cross-sections, both
base periods, `anticipation` 0 and 1, with and without covariates
(`tests/reference_parity/test_cs_notyet_cutoff_parity.py`). The old rule is
what Stata `csdid, asinr` does, and the docs had called it the R
convention: on `mpdta` R `did` gives ATT(2007, 2004) = 0.033813 (the csdid
default number), not asinr's 0.032971. It remains available as
`notyet_cutoff="asinr"`, still pinned to Stata. `notyet_cutoff` is now also
honoured for repeated cross-sections, where it was ignored.
- **`sp.did_bcf` compared treated and never-treated units over different
calendar windows**: never-treated units were split at the median period
whatever the cohort's adoption date, so any non-linear common trend
entered the ATT (common quadratic trend, true ATT 1: cohort CATTs 0.01 and
1.64, now 1.07 and 1.03). Each cohort is now compared with never-treated
units over its own pre and post windows. The covariate-path standard error
(the average posterior spread divided by sqrt(n)) is replaced by the
bootstrap standard deviation of the cohort ATT, combined conservatively
across cohorts; the silent fallback to a difference in means on any
exception is removed. The docstring no longer attributes the estimator to
[souto2025forests], whose model is specified in outcome levels.

⚠️ JSS Track A closure: the last open cross-language rows (StatsPAI vs R / Stata)

Closes the items the JSS manuscript's tolerance ledger still carried as open
or convention-only after phase 3. Every number below is on the committed
same-byte fixtures of `tests/r_parity/` and `tests/stata_parity/`; both
reference sides re-derive with zero drift (R 89/89, Stata 85/85, Python 89/89).

⚠️ Correctness

- **`sp.gardner_did` / `sp.did_2stage`: the default standard error is now the
did2s corrected two-stage clustered variance.** The previous default
(`vce='analytic'`) clustered only the Stage-2 residuals, treating the imputed
counterfactual as known, and understated uncertainty (~26% on `mpdta`; the
docstring warned of ~0.78 coverage at a nominal 95%). The default now builds
the Stage-2 cluster sandwich from the influence function of *both* stages,
exactly as R `did2s::did2s` 1.2.1 and Stata `did2s` do
[gardner2022twostage; butts2022stage]: with `X1` the Stage-1 design, `X10`
that design with treated rows zeroed, `e1` the Stage-1 residual and `X2` /
`e2` the Stage-2 design and residual, `gamma = (X10'X10)^{-1} X1'X2`,
`s_g = sum_{i in g} (x2_i e2_i - gamma' x10_i e1_i)`,
`V = (X2'X2)^{-1} [sum_g s_g s_g'] (X2'X2)^{-1}`, with no small-sample cluster
factor (neither reference applies one). On the mpdta replica the static-ATT
SE moves from 0.0051177 to 0.0069036 and agrees with R `did2s` at rel
2.7e-10 and Stata `did2s` at 1.3e-14; all 26 event-study horizon SEs match
`did2s` at 2.1e-14. Point estimates are unchanged in every mode. The
event-study overall ATT SE is now the delta-method SE through the full
cross-horizon covariance (`model_info['event_study']['vcov']`) instead of
treating horizons as independent. `vce='stage2'` reproduces the pre-1.29
number bit for bit (with a warning); `model_info['se_convention']` records
which variance was used. Track A module `73_did2s` registers `rel_se: 1e-6`
and its SE is a strict three-way parity row rather than a documented gap.
See MIGRATION.md.

Added

- **`sp.aggte(..., share_variance=True)`.** The aggregated Callaway--Sant'Anna
variance carries the sampling variability of the estimated cohort shares
(R `did:::wif`, the authors' convention) by default; `share_variance=False`
holds the shares fixed, which is what Stata `csdid`'s `estat group`
`GAverage` does. Point estimates are identical either way; only the
`'group'` overall (and any multi-cohort cell) SE moves. On the `mpdta`
replica the fixed-share group-overall SE joins `csdid` at rel 6.2e-15, so
Track A module 04 now pins the Stata convention as a joining row
(`group_overall_fixedshare`) instead of a 0.27% note.
`result.model_info["share_variance"]` records the choice.
- **`sp.match` / `sp.psm`: `se_method='abadie_imbens_2016'`**, the Abadie &
Imbens (2016) variance for nearest-neighbour ATT matching on an *estimated*
logit propensity score [abadie2016matching] -- what Stata
`teffects psmatch (y) (t x, logit), atet` reports:
`var = sigma^2_delta - c'V_gamma c + d'V_gamma d` with `sigma^2_delta` the
population-ATT base term (`[sum_treated (y_i - yhat0_i - delta)^2 +
sum_i xi2_i (K_i^2 - K'_i)] / N1^2`), `c`, `d` the score-derivative vectors
of the Stata manual ([CAUSAL] teffects nnmatch, "PSM, ATE, and ATET variance
adjustment") and `V_gamma` the logit observed-information covariance (new
module `src/statspai/matching/_ai2016.py`; the terms are returned in
`model_info['ai2016_components']`). Conventions pinned against Stata 18: the
within-arm conditional-variance set is the unit plus its `ai_matches`
nearest same-arm units, ties included (`ai_matches=J` is Stata `nn(J+1)`);
the extra clustering behind `d(delta)/d(gamma)` is Euclidean on the raw
covariates with ties; neighbour sets are found on the unclipped score. On
`sp.datasets.nsw_dw()` (Track A `11_psm`) the SE is 621.79328 against
Stata's 621.79325 (rel 6.4e-8; with `ai_matches=2`, Stata `nn(3)`, rel
< 1e-6), and plugging Stata's own `e(bps)`/`e(Vps)` into the formula
reproduces `e(V)` to 1e-12, so the residual is Stata's ML stopping rule.
Term by term: base SE 673.42, `c'Vc = 74,905`, `d'Vd = 8,036` -- the
estimated-score term is *negative* for the ATT here, which is why the
corrected SE sits below the psmatch2 `ai(1)` SE (643.35). The default
(`'abadie_imbens'`, Abadie--Imbens 2006 sample-ATT, = `psmatch2, ai(J)`) is
unchanged; the 3.4% between the two Stata commands is now a documented
convention with two named sources (population- vs sample-ATT target, and
the estimated-score term only `teffects` charges for) rather than an open
item. Requires `distance='propensity'`, `method='nearest'`,
`estimand='ATT'`, `bias_correction=False`; anything else raises
`MethodIncompatibility`.
- **`sp.fect(cv=True)`**: fect's cross-validation for the factor number
(`method='ife'`, grid `r_range=(0, 5)`) and the nuclear-norm penalty
(`method='mc'`, fect's `nlambda=10` log grid or `lambda_grid=`), with fect's
rolling / block / treated-units holdouts (`cv_method`, `cv_prop`, `cv_nobs`,
`cv_donut`, `cv_buffer`, `k`), criteria `mspe` / `gmspe` / `moment` / `pc`
and the `cv_rule` (`1se` / `min` / `1pct`) re-pick; the CV table is in
`model_info["cv"]`, the selection in `r_cv` / `lambda_cv`; the final fit is
refitted at the selected value, so `sp.fect(cv=True)` equals
`sp.fect(r=<selected>)` bit for bit. Evidence: the CV scoring equals fect's
on identical folds (rel 5.9e-8 ife, 6.5e-9 mc) and the curve and selection
agree with R across 20 fold seeds (T3;
`tests/reference_parity/test_fect_interflex_cv_parity.py`). Note that
`fect()`'s user-facing default is `cv.method = "rolling"` (the `"all_units"`
in `fect.default`'s signature is overridden), which the port follows.
- **`sp.interflex(estimator='kernel', cv=True)`**: interflex's bandwidth
cross-validation (`kfold`, `bw_grid`, `metric`, `random_state`); the CV table
is in `model_info["cv"]`, the selection in `bw_cv`. Evidence: per-fold
losses equal interflex's `getError.CV` on identical folds (rel 2.8e-13);
curve and selection agree with R across 20 fold seeds (T3). The JSS
manuscript's "CV selectors not ported" limitation is closed.
- **Stata bridges for Track A `18_augsynth` and `19_gsynth`** (audited
Stata/Mata algorithm bridges, the same class as `08_dml`). `18_augsynth.do`
reproduces `augsynth` 0.2.0's Ridge ASCM (simplex-QP SCM weights on
donor-mean-centred pre-period outcomes, the 21-point lambda path with
leave-one-period-out folds and the 1-SE rule recomputed in Mata, ridge
correction): ATT within 1.7e-8 of StatsPAI and 7.9e-6 of R, where the R
residual is `augsynth::synth_qp`'s OSQP stopping tolerance (eps 1e-8; at
1e-13 augsynth agrees with both exact solvers to 1e-11). `allsynth` (SSC
1.32) refuses the fixture (`K + 2` donors) and, where it runs, is a
different estimand. `19_gsynth.do` implements gsynth's control-only factor
convention (two-way demeaned control panel, rank-r SVD, treated-unit
pre-period projection on `[1, F]`, r fixed at the R `r.cv`): ATT within
2.2e-15 of StatsPAI and 7.4e-14 of R. `fect_stata method(ife)` is a
different estimator (IFEct EM over all untreated cells; converged R
`fect(method="ife")` -0.33206 vs gsynth -0.32417) and is emitted as the
unjoined `ife_att_avg_fect_stata` diagnostic; `fect_stata` has no
`method(gsynth)`. `compare.py::STATA_SKIP_REASON` is down to four modules
(13, 77, 79, 80); 85 of 89 Track A modules carry a Stata reference.

Fixed

- **`sp.rdplot` raised `IndexError` when the largest running-variable value
fell one ulp past the last evenly spaced bin edge** (the phase-3
`_rdplot_core` port of `rdrobust::rdplot`; on the Lee-Senate replica
`x_max = 0.8635116005610454` while `seq(c, x_max, range/J)` ended at
`0.8635116005610453`). The R port's `seq` now snaps its last element to
the endpoint, so the maximum observation is inside the last bin as
`findInterval(rightmost.closed = TRUE)` intends; bin means and counts are
unchanged for every other observation. Found by the JSS figure script
(`ex04_lee_rd`), which is why the worked examples are executed rather than
summarised.

Changed (parity harness)

- Track A `27_glmm_aghq` moves from the machine tier to the iterative tier
(`rel_est` 1e-6 -> 5e-6). Not a convention gap: the `logLik` row joins
`lme4` at 1e-13 and `lme4`'s own deviance function reproduces StatsPAI's
optimum to 5e-11, but the 8-point AGQ objective is flat along (intercept,
sqrt(var)) and `lme4`'s optimisers (`bobyqa`, `nloptwrap`, `optimx/nlminb`)
scatter by 1.1e-6 to 1.5e-6 in the intercept at the same deviance; Stata
`melogit` (`mvaghermite`) sits 2e-8 from StatsPAI and all sides agree to
3e-8 at 30 quadrature points. The contract test now also guards `logLik` at
1e-9 so quadrature drift cannot hide behind the optimiser-noise budget.
Registered in `docs/dev/r_parity_tolerances.md` ("Loosening applied
2026-09-22"). Strictness tiers are now 79 machine / 8 iterative / 1
moderate / 1 T4.
- Track A `18_augsynth`'s 2e-5 point budget is re-graded A (mechanism: OSQP
tolerance on the R side), unchanged in value.
- `sp.gardner_did` no longer emits the "0/n bootstrap replicates succeeded"
`RuntimeWarning` in static mode (the per-coefficient series was never
filled there; the reported SE was always the correct overall bootstrap SE),
and warns when the cluster variable has a single level under the default
`vce`.
- `paper.bib` gains `abadie2016matching` (Abadie & Imbens 2016, Econometrica
84(2), 781-807, DOI 10.3982/ECTA11293; refs verified via Crossref and
doi.org).