feat(diagnostics): VAR-aware convergence and stability report

CONTEXT.md

@@ -121,6 +121,21 @@ _Avoid_: "order of differencing" for `d_max` — `d_max` is the maximum across t
 **Cointegration rank**:
 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`.
+_Avoid_: "diagnostics" as a synonym for this one object — the diagnostics family is wider, and the name `.diagnostics()` is reserved.
+
+**Parameter block**:
+A group of posterior variables that share a role in the model and tend to share sampling behaviour: `coefficient`, `intercept`, `exog`, `covariance`, `volatility`, `identification`, `other`. The unit of attribution in the convergence report — a mixing problem belongs to a block, not to the model as a whole. Assignment is three-tiered (static name map, then the `v{i}_` stochastic-volatility prefix, then the optional `posterior_var_names()` capability on the volatility process), with unrecognised variables falling to `other` rather than raising.
+_Avoid_: "parameter group" / "variable family" — the report's field is `block`.
+
+**Companion matrix**:
+The `(n·p, n·p)` matrix that rewrites a VAR(p) as a first-order system: the lag coefficients `B` form its top block row verbatim and sub-diagonal identity blocks shift the lag state. Built by `companion_matrix(B, n_lags)`; the intercept and exogenous block play no part in it.
+_Avoid_: bare "companion" — always say "companion matrix" in full (see Flagged ambiguities).
+
+**Spectral radius / explosive draw**:
+The largest companion-matrix eigenvalue modulus of a single posterior draw, computed by `spectral_radius(B, n_lags)`. A draw is *stable* below 1 and *explosive* at or above it: an explosive draw's impulse responses diverge with the horizon, its forecast fan is unbounded, its long-horizon FEVD shares are uninterpretable, and its historical-decomposition baseline drifts. Some explosive mass is legitimate on level data under a random-walk prior mean, which is why the convergence report warns on it and never fails.
+_Avoid_: "unstable" for a whole posterior — stability is a per-draw property, summarised by the explosive *fraction* (`p_explosive`).
 
 ## Relationships
 
@@ -138,6 +153,8 @@ _Avoid_: "number of cointegrating vectors" in API surface (fine in prose); "coin
 - A **stationarity pretest** consumes `VARData` (endogenous block only), a DataFrame, or a Series, and produces a result object — never a modified dataset and never a specification. It sits *beside* the pipeline, not in it: nothing downstream of `VAR.fit()` reads its output.
 - **Integration order** feeds **cointegration rank**: the Johansen test is only meaningful for series that are individually integrated, and it is conditioned on a lag order (`k_ar_diff = p - 1`) that `select_lag_order` supplies.
 - A **ConjugateVAR** carries an **NIW prior** and optionally a **deterministic volatility break**; a **VAR** carries a **MinnesotaPrior** and a **PyMC volatility process** (`PyMCVolatilityProcess`, the `build_pymc_latent` extension of the `VolatilityProcess` query surface). Each estimator's fields accept only its compatible components, enforced by types + validators rather than a builder.
+- A **FittedVAR** produces a **convergence report** on its own; identification adds nothing that needs diagnosing, so `IdentifiedVAR.convergence_report()` returns the same reduced-form answer. The report partitions the posterior into **parameter blocks**, consulting the **volatility process** for the variables it registered.
+- A **convergence report** carries the **spectral radius** of every draw, computed from the **companion matrix** built out of the same `B` that drives the moving-average recursion — so convergence and dynamic stability are answered from one object rather than two.
 
 ## Example dialogue
 
@@ -173,3 +190,5 @@ _Avoid_: "number of cointegrating vectors" in API surface (fine in prose); "coin
 - "Minnesota prior" now denotes two distinct encodings: the independent-Normal `MinnesotaPrior` (NUTS path) and the conjugate `NIWPrior` (`ConjugateVAR`). Name the estimator when it matters.
 - "Σ" 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.

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

@@ -0,0 +1,42 @@
+# The convergence report is block-structured, and explosive draws never fail it
+
+`convergence_report` is a VAR-specific diagnostic object rather than a wrapper over `arviz.summary`. It makes three commitments: every sampling metric is reported per *parameter block* with the offending coordinate named; dynamic stability is reported alongside convergence, computed from the companion matrix of every draw; and the two VAR-specific failure modes get named, machine-readable messages carrying remedies. A `"failed"` status is reserved for sampler pathology — R-hat above 1.05, effective sample size below 100, or a divergence rate at or above 1% — and is never triggered by explosive draws.
+
+## Block taxonomy
+
+Blocks are `coefficient`, `intercept`, `exog`, `covariance`, `volatility`, `identification`, `other`, reported in that order and omitted when empty. A single worst-R-hat over a whole VAR posterior hides which part of the model is failing: the lag coefficients, the covariance parameterisation, and the stochastic-volatility latents mix at very different rates, and a user staring at `max_rhat = 1.4` learns nothing about what to change.
+
+Resolution is three-tiered, first match wins:
+
+1. A static map of the posterior variable names Impulso's own estimators register (`B`, `intercept`, `B_exog`, `sigma_sd`, `tril_offdiag`, `L`, `Sigma`, `h`, `R_chol`, `R_chol_offdiag`, `structural_shock_matrix`, `P`).
+2. The `v{i}_` prefix carried by every per-variable stochastic-volatility latent, which covers the whole family — present adapters and future ones — without enumerating parameter names that change whenever a dynamics adapter gains a field.
+3. The optional `posterior_var_names()` capability on `VolatilityProcess`, letting an adapter claim the variables it registered. It is documented as an optional capability in the mould of `IdentificationScheme._samples_rotations`, read through `getattr`, and is deliberately *not* a protocol requirement — a third-party adapter that omits it still works.
+
+Anything unresolved lands in `other`, never an error. A hand-built or third-party posterior still gets a full report, and the block's variable list makes plain what was not recognised. Refusing to diagnose a posterior because one variable is unfamiliar would be the wrong trade.
+
+The `identification` block exists but is normally empty: the structural shock matrix is memoised lazily on `IdentifiedVAR` and never written back to the posterior. Excluding it is deliberate rather than incidental. Under `Cholesky` it is a deterministic function of draws already diagnosed in the covariance block, so its R-hat adds nothing; under `SignRestriction` a fresh rotation is drawn per call, so its R-hat would describe the rotation sampler rather than the posterior — actively misleading. The block is kept for legacy and hand-built posteriors that do carry the variable.
+
+## Thresholds
+
+| Metric | Warn | Fail | Source |
+| --- | --- | --- | --- |
+| 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) |
+| 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.
+
+## 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.
+
+## 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.

docs/reference/diagnostics.md

@@ -0,0 +1,25 @@
+# Diagnostics
+
+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
+`IdentifiedVAR.convergence_report()`; the free function is the entry point
+for posteriors built by hand.
+
+```{eval-rst}
+.. currentmodule:: impulso.diagnostics
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   convergence_report
+   ConvergenceReport
+   BlockDiagnostics
+   StabilitySummary
+   DiagnosticMessage
+   ConvergenceThresholds
+   assign_blocks
+```

docs/reference/index.md

@@ -20,6 +20,7 @@ identification
 scenario
 results
 evidence
+diagnostics
 primitives
 protocols
 plotting

docs/reference/primitives.md

@@ -1,8 +1,8 @@
 # Primitives
 
-Moving-average building blocks shared by the IRF, FEVD, and
-dynamic-multiplier machinery, published for downstream libraries that
-compose with Impulso posteriors.
+Moving-average and companion-form building blocks shared by the IRF, FEVD,
+dynamic-multiplier, and stability machinery, published for downstream
+libraries that compose with Impulso posteriors.
 
 ```{eval-rst}
 .. currentmodule:: impulso
@@ -13,4 +13,6 @@ compose with Impulso posteriors.
 
    compute_ma_phi
    lag_matrices
+   companion_matrix
+   spectral_radius
 ```

docs/references.bib

@@ -156,3 +156,20 @@ @article{blanchardQuah1989
   number  = {4},
   pages   = {655--673},
 }
+
+@article{vehtari2021,
+  author  = {Vehtari, Aki and Gelman, Andrew and Simpson, Daniel and Carpenter, Bob and B\"urkner, Paul-Christian},
+  title   = {Rank-Normalization, Folding, and Localization: An Improved $\widehat{R}$ for Assessing Convergence of MCMC},
+  journal = {Bayesian Analysis},
+  year    = {2021},
+  volume  = {16},
+  number  = {2},
+  pages   = {667--718},
+}
+
+@misc{betancourt2017,
+  author = {Betancourt, Michael},
+  title  = {A Conceptual Introduction to Hamiltonian Monte Carlo},
+  year   = {2017},
+  note   = {arXiv:1701.02434},
+}

src/impulso/__init__.py

@@ -12,8 +12,17 @@
 
     from impulso._linalg import lag_matrices
     from impulso._ma import compute_ma_phi
+    from impulso._stability import companion_matrix, spectral_radius
     from impulso.conjugate import ConjugateVAR
     from impulso.conjugate_volatility import ConjugateVolatility, PandemicBreak
+    from impulso.diagnostics import (
+        BlockDiagnostics,
+        ConvergenceReport,
+        ConvergenceThresholds,
+        DiagnosticMessage,
+        StabilitySummary,
+        convergence_report,
+    )
     from impulso.evidence import EvidenceComparison, ModelEvidence, compare_evidence
     from impulso.fitted import FittedVAR
     from impulso.identification import Cholesky, LongRunRestriction, ProxySVAR, SignRestriction
@@ -48,13 +57,17 @@
 
 __all__ = [
     "VAR",
+    "BlockDiagnostics",
     "Cholesky",
     "CointegrationTestResult",
     "ConditionalForecastResult",
     "ConjugateVAR",
     "ConjugateVolatility",
     "Constant",
+    "ConvergenceReport",
+    "ConvergenceThresholds",
     "CounterfactualResult",
+    "DiagnosticMessage",
     "DynamicMultiplierResult",
     "ErrorDistribution",
     "EvidenceComparison",
@@ -82,6 +95,7 @@
     "ScenarioResult",
     "ShockPath",
     "SignRestriction",
+    "StabilitySummary",
     "StationarityTestResult",
     "StochasticVolatility",
     "StudentT",
@@ -90,14 +104,17 @@
     "VolatilityProcess",
     "VolatilityResult",
     "adf_test",
+    "companion_matrix",
     "compare_evidence",
     "compute_ma_phi",
+    "convergence_report",
     "enable_runtime_checks",
     "integration_order",
     "johansen_test",
     "kpss_test",
     "lag_matrices",
     "select_lag_order",
+    "spectral_radius",
 ]
 
 
@@ -145,6 +162,14 @@
     "VolatilityProcess": "impulso.protocols",
     "compute_ma_phi": "impulso._ma",
     "lag_matrices": "impulso._linalg",
+    "companion_matrix": "impulso._stability",
+    "spectral_radius": "impulso._stability",
+    "convergence_report": "impulso.diagnostics",
+    "ConvergenceReport": "impulso.diagnostics",
+    "ConvergenceThresholds": "impulso.diagnostics",
+    "BlockDiagnostics": "impulso.diagnostics",
+    "DiagnosticMessage": "impulso.diagnostics",
+    "StabilitySummary": "impulso.diagnostics",
 }
 """Map of lazily-exported name to the module that defines it.
 
@@ -229,6 +254,7 @@ def enable_runtime_checks() -> None:
     from beartype.roar import BeartypeDecorHintPep484585Exception
 
     import impulso.data
+    import impulso.diagnostics
     import impulso.fitted
     import impulso.identified
     import impulso.spec
@@ -242,6 +268,7 @@ def enable_runtime_checks() -> None:
         impulso.spec,
         impulso.fitted,
         impulso.identified,
+        impulso.diagnostics,
         impulso.sv.data,
         impulso.sv.spec,
         impulso.sv.fitted,