feat(diagnostics): E-BFMI, treedepth saturation, stability plot, reserved-prefix doc

CONTEXT.md

@@ -122,7 +122,7 @@ _Avoid_: "order of differencing" for `d_max` — `d_max` is the maximum across t
 The number of independent long-run relationships among integrated series, from the Johansen procedure (`johansen_test`). Both sequential tests are reported — `rank_trace` and `rank_max_eigen` — and `rank` is `rank_trace` by documented convention. Decisions rest on **critical values, not p-values** (MacKinnon-Haug-Michelis 1996 tables, as vendored by statsmodels), which is why `alpha` is restricted to 0.10 / 0.05 / 0.01. Rank ≥ 1 means differencing every series discards the long-run relationship. A vector error-correction model (VECM) is **out of scope**; the recommended response is a VAR in levels (the Sims–Stock–Watson stance; the Minnesota prior already shrinks toward random walks).
 _Avoid_: "number of cointegrating vectors" in API surface (fine in prose); "cointegration test" without saying which statistic, since trace and max-eigen can disagree.
 **Convergence report**:
-The VAR-aware verdict on whether a fitted posterior is usable, produced by `convergence_report()` (or the delegating `FittedVAR.convergence_report()` / `IdentifiedVAR.convergence_report()`). It reports R-hat and both effective sample sizes per *parameter block* with the worst coordinate named, the global divergence count, and the posterior distribution of the spectral radius, and carries machine-readable `DiagnosticMessage` codes for the two VAR-specific failure modes. Its `status` — `"passed"` / `"warnings"` / `"failed"` — reserves `"failed"` for sampler pathology; explosive draws warn but never fail. See `docs/adr/0008-convergence-report-blocks-and-thresholds.md`.
+The VAR-aware verdict on whether a fitted posterior is usable, produced by `convergence_report()` (or the delegating `FittedVAR.convergence_report()` / `IdentifiedVAR.convergence_report()`). It reports R-hat and both effective sample sizes per *parameter block* with the worst coordinate named, the global divergence, energy (E-BFMI) and max-treedepth statistics, and the posterior distribution of the spectral radius, and carries machine-readable `DiagnosticMessage` codes for the two VAR-specific failure modes. Its `status` — `"passed"` / `"warnings"` / `"failed"` — reserves `"failed"` for sampler pathology; explosive draws, low E-BFMI and treedepth saturation warn but never fail. See `docs/adr/0008-convergence-report-blocks-and-thresholds.md`.
 _Avoid_: "diagnostics" as a synonym for this one object — the diagnostics family is wider, and the name `.diagnostics()` is reserved.
 
 **Parameter block**:
@@ -191,4 +191,4 @@ _Avoid_: "unstable" for a whole posterior — stability is a per-draw property,
 - "Σ" now means the *scale* matrix under `StudentT` errors and the covariance under `Gaussian` errors. `sigma()` returns the same object either way; when the number has to be a variance, say so and use `innovation_covariance()`.
 - "Counterfactual" in the wider literature spans shock-path edits (Impulso's meaning), policy-rule replacement (Sims–Zha style; out of scope), and Lucas-robust constructions (McKay–Wolf; out of scope). When comparing with external work, say which one is meant.
 - "Companion" is overloaded: the *companion matrix* is the stacked first-order form of a VAR(p), while the ADPRR "calibrated companion" `q_cal` is the plausibility statistic's partner quantity. They share nothing. Always write "companion matrix" in full; never shorten it to "the companion".
-- `StabilitySummary` (the convergence report's spectral-radius block) is distinct from the `StabilityResult` planned for the ecological-stability work: the former summarises one scalar per draw for a diagnostic verdict, the latter will carry the full complex eigenvalue spectrum and the reactivity/return-rate measures derived from it. Both read `companion_eigenvalues`; neither subsumes the other.
+- `StabilitySummary` (the convergence report's spectral-radius block) is distinct from the `StabilityResult` planned for the ecological-stability work: the former summarises one scalar per draw for a diagnostic verdict, keeping a capped subset of the complex roots only so `plot()` has something to scatter, whereas the latter will carry the full spectrum for every draw and the reactivity/return-rate measures derived from it. Both read `companion_eigenvalues`; neither subsumes the other.

docs/adr/0008-convergence-report-blocks-and-thresholds.md

@@ -23,20 +23,31 @@ The `identification` block exists but is normally empty: the structural shock ma
 | R-hat | 1.01 | 1.05 | Vehtari et al. (2021); classic Gelman–Rubin |
 | Effective sample size | 400 | 100 | 100 per chain at four chains |
 | Divergence rate | any divergence | 1% | Betancourt (2017) |
+| E-BFMI | 0.3 | *never* | Betancourt (2016), arXiv:1604.00695 |
+| Max-treedepth saturation rate | 1% | *never* | — |
 | Explosive draw fraction | 5% | *never* | — |
 
-R-hat and ESS comparisons are strict, so a metric sitting exactly on a threshold passes; the two rate thresholds (divergence rate, explosive fraction) trigger at the boundary. Thresholds live in a frozen `ConvergenceThresholds` model rather than as module constants so a caller can tighten them for a specific study and the report echoes back what it used.
+R-hat, ESS and E-BFMI comparisons are strict, so a metric sitting exactly on a threshold passes; the three rate thresholds (divergence rate, treedepth saturation rate, explosive fraction) trigger at the boundary. Thresholds live in a frozen `ConvergenceThresholds` model rather than as module constants so a caller can tighten them for a specific study and the report echoes back what it used.
+
+The treedepth threshold has no literature source, unlike the others. Stan and PyMC surface any saturation at all, but a handful of hits in a long run costs wall-clock time and nothing else, so reporting them would train users to skim past the message. One transition in a hundred is the point at which the sampler is spending real effort on trajectories it never gets to finish. Both backends record the flag under different names — `reached_max_treedepth` in PyMC, `maxdepth_reached` in nutpie — so this is the one statistic the report resolves through a name map rather than reading directly.
 
 ## Why explosive draws never fail
 
 Posterior mass on parameter draws whose companion matrix has spectral radius at or above 1 is reported prominently, with its consequences (impulse responses that diverge with the horizon, unbounded forecast fans, uninterpretable long-horizon FEVD shares, drifting historical-decomposition baselines) and its remedies. It is still only a warning, and at fractions below `explosive_warn` only informational.
 
 The reason is that explosiveness is a property of the *model*, not of the sampler. Macroeconomic data in levels under a Minnesota prior centred on a random walk puts substantial mass near the unit circle by construction; that is the prior doing its job, and a fraction of draws crossing it is expected rather than pathological. Failing the report there would train users to ignore `"failed"`, which must keep meaning "these draws do not describe the posterior". Convergence and stability are different questions and are reported as such.
 
+## Why E-BFMI and max-treedepth never fail either
+
+Both are statements about *efficiency*, not about wrongness. A saturated tree depth means NUTS stopped a trajectory before it turned back on itself, so the draws are more autocorrelated than they need to be; a low E-BFMI means momentum resampling is exploring the energy distribution slowly, so the tails are undersampled relative to the bulk. Neither says the retained draws come from the wrong distribution — unlike a divergence, which says the sampler could not follow the geometry at all, or an unmixed R-hat, which says the chains are not describing one distribution.
+
+Both are also remediable by changing the sampler or the parameterisation without touching the model, so failing on them would block work that is merely slower than it should be. Both metrics are on the report regardless, so a caller who wants them fatal reads `ebfmi` and `treedepth_saturation_rate` and decides.
+
 ## Rejected alternatives
 
 - **A thin wrapper over `az.summary`.** Rejected: it produces one row per coordinate with no block structure, no stability, and no VAR-specific interpretation — exactly the output users already have and cannot act on.
 - **Living in `results.py` alongside the other result objects.** Rejected: `VARResultBase` contracts for `median`/`hdi`/`to_dataframe`/`plot` over a posterior-predictive DataArray, and a convergence report has no such array. Following `LagOrderResult`'s precedent would have forced a fake `plot` and a fake `median`. A dedicated `diagnostics.py` also gives the diagnostics family (issue #57's umbrella) somewhere to grow.
 - **Per-block divergence attribution.** Rejected: a divergence is a property of a trajectory through the whole parameter space. Splitting the count by block would invent an attribution the sampler never made.
 - **Warning through `warnings.warn`.** Rejected: the report object carries `status`, `messages`, and the per-block table, so the caller decides whether to print, raise, or ignore. A diagnostic that emits warnings cannot be used inside a loop over model specifications.
-- **A `.plot()` method in v1.** Deferred: issue #57 owns diagnostic visuals, and the raw `(chain, draw)` radius array is exposed so a histogram or unit-circle scatter is a few lines away.
+- **A `.plot()` method in v1.** Deferred, then added on `StabilitySummary` alone (issue #179): `report.stability.plot()` gives the spectral-radius histogram beside the companion-root scatter, matching how every other result object in the library is plotted. `ConvergenceReport` itself still has no `plot` — a block table of R-hat and ESS is a table, and issue #57 owns whatever diagnostic visuals go beyond stability.
+- **Retaining every companion eigenvalue on `StabilitySummary`.** Rejected: the array grows as `draws × n_vars × n_lags` complex numbers — roughly 15 MB for 4000 draws of a 240×240 companion matrix — on an object whose every statistic is derived from the radii. The summary keeps a chain-pooled, deterministically strided subset of at most 200 draws, computed from the single eigendecomposition the radii already require, which is more points than the scatter panel can distinguish and a fixed ceiling regardless of posterior size.

docs/reference/diagnostics.md

@@ -2,11 +2,13 @@
 
 Convergence and dynamic-stability diagnostics for a fitted VAR posterior.
 `convergence_report` reports R-hat and effective sample size *per parameter
-block* with the offending coordinate named, counts divergences globally, and
-summarises the posterior distribution of the companion-matrix spectral
-radius. Reach for it through `FittedVAR.convergence_report()` or
+block* with the offending coordinate named; counts divergences, energy
+pathologies (E-BFMI) and max-treedepth saturation globally; and summarises
+the posterior distribution of the companion-matrix spectral radius. Reach
+for it through `FittedVAR.convergence_report()` or
 `IdentifiedVAR.convergence_report()`; the free function is the entry point
-for posteriors built by hand.
+for posteriors built by hand. `report.stability.plot()` draws the radius
+posterior beside the companion roots on the unit circle.
 
 ```{eval-rst}
 .. currentmodule:: impulso.diagnostics

docs/reference/plotting.md

@@ -17,4 +17,5 @@
    plot_counterfactual
    plot_volatility
    plot_sv_forecast
+   plot_stability
 ```