From e6c5b22e024a62b76199b1ec1daf290a2fc17b10 Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 04:01:31 +0200 Subject: [PATCH 1/6] feat(diagnostics): companion-matrix stability primitives (#142) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add `_stability.py` as the single source of truth for the VAR companion form. `B` already stores lag blocks concatenated in lag order along its trailing axis — exactly the companion matrix's top block row — so the coefficients are copied in verbatim with no slicing. `companion_eigenvalues` returns the raw complex roots rather than only their moduli, because downstream reactivity and return-rate measures need the imaginary parts. The batch is processed in chunks: cost scales as O(N * (n * p)^3), so a large posterior over a large system allocates heavily without a bound on peak memory. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Egjd7ToFeb9TQqFnfRQZxV --- src/impulso/_stability.py | 143 ++++++++++++++++++++++++++++++++++++++ tests/test_stability.py | 124 +++++++++++++++++++++++++++++++++ 2 files changed, 267 insertions(+) create mode 100644 src/impulso/_stability.py create mode 100644 tests/test_stability.py diff --git a/src/impulso/_stability.py b/src/impulso/_stability.py new file mode 100644 index 0000000..a05990d --- /dev/null +++ b/src/impulso/_stability.py @@ -0,0 +1,143 @@ +"""Companion-form stability primitives for reduced-form VAR coefficients. + +A VAR(p) is stable when every root of its companion matrix lies strictly +inside the unit circle. Stacking the lag matrices `A_1, ..., A_p` into the +first block row of a `(n*p, n*p)` matrix turns the p-th order system into a +first-order one:: + + F = [[A_1, A_2, ..., A_{p-1}, A_p], + [ I , 0 , ..., 0 , 0 ], + [ 0 , I , ..., 0 , 0 ], + [ ... ... I , 0 ]] + +Impulso stores the lag coefficients as a single matrix `B` of shape +`(n, n*p)` whose trailing axis concatenates the lag blocks in lag order — +exactly the layout of the companion matrix's top block row — so `B` is +copied in verbatim, with no slicing. The intercept and any exogenous block +live in separate posterior variables and play no part in stability. + +This module is the single source of truth for the companion form. It is +deliberately rank-agnostic in the same way as `lag_matrices` and +`sigma_from_cholesky`: only the trailing two axes are interpreted, so a +single draw `(n, n*p)` and a posterior tensor `(chain, draw, n, n*p)` both +work without branching. `companion_eigenvalues` returns the raw complex +roots rather than only their moduli, because downstream consumers (return +rates and reactivity measures borrowed from theoretical ecology) need the +imaginary parts too. +""" + +from __future__ import annotations + +import numpy as np + + +def _companion_dims(B: np.ndarray, n_lags: int) -> tuple[int, int]: + """Validate the coefficient layout and return `(n_vars, n_vars * n_lags)`. + + Raises: + ValueError: If `n_lags` is not positive, if the trailing axis of *B* + is not divisible by `n_lags`, if *B* has fewer than two axes, or + if the resulting block row is not `n_vars` rows tall. + """ + if n_lags < 1: + raise ValueError(f"n_lags must be positive, got {n_lags}") + if B.ndim < 2: + raise ValueError(f"B must have at least 2 dimensions (n, n * n_lags), got shape {B.shape}") + n_coeffs = B.shape[-1] + if n_coeffs % n_lags != 0: + raise ValueError(f"B trailing axis {n_coeffs} is not divisible by n_lags {n_lags}") + n_vars = n_coeffs // n_lags + if B.shape[-2] != n_vars: + raise ValueError( + f"B must have shape (..., n, n * n_lags); trailing axis {n_coeffs} with " + f"n_lags {n_lags} implies n = {n_vars}, but B has {B.shape[-2]} rows" + ) + return n_vars, n_coeffs + + +def companion_matrix(B: np.ndarray, n_lags: int) -> np.ndarray: + """Build the first-order companion matrix of a VAR(p). + + Args: + B: Stacked lag coefficients with trailing shape `(n, n * n_lags)`, + lag blocks concatenated in lag order (lag 1 first). Leading axes + are arbitrary batch dimensions, typically `(chains, draws)`. + n_lags: Number of lag blocks stacked along the trailing axis. + + Returns: + Companion matrix with trailing shape `(n * n_lags, n * n_lags)` and + the same leading batch axes as *B*. The top `n` rows are *B* + verbatim; the sub-diagonal identity blocks shift the lag state. + + Raises: + ValueError: If the coefficient layout is inconsistent with `n_lags`. + """ + B = np.asarray(B) + n_vars, m = _companion_dims(B, n_lags) + F = np.zeros((*B.shape[:-2], m, m), dtype=np.result_type(B.dtype, np.float64)) + F[..., :n_vars, :] = B + if n_lags > 1: + shift = np.arange(m - n_vars) + F[..., n_vars + shift, shift] = 1.0 + return F + + +def companion_eigenvalues(B: np.ndarray, n_lags: int, *, chunk_size: int = 256) -> np.ndarray: + """Eigenvalues of the companion matrix, one set per batch element. + + Cost scales as `O(N * (n * p)**3)` for `N` batch elements, so a large + posterior over a large system is expensive: 4000 draws of an 8x8 + companion matrix take about 0.03 s, but 200 draws of a 240x240 one take + about 2 s and allocate roughly 90 MB. The batch is therefore processed in + chunks, and callers holding many draws should thin them first. + + Args: + B: Stacked lag coefficients with trailing shape `(n, n * n_lags)`. + n_lags: Number of lag blocks stacked along the trailing axis. + chunk_size: Number of companion matrices to materialise at once. + Bounds peak memory; it does not change the result. + + Returns: + Complex array with trailing axis of length `n * n_lags` and the same + leading batch axes as *B*. Eigenvalues are unordered within each set, + as returned by `numpy.linalg.eigvals`. + + Raises: + ValueError: If the coefficient layout is inconsistent with `n_lags`, + or if `chunk_size` is not positive. + """ + B = np.asarray(B) + n_vars, m = _companion_dims(B, n_lags) + if chunk_size < 1: + raise ValueError(f"chunk_size must be positive, got {chunk_size}") + + leading = B.shape[:-2] + flat = B.reshape(-1, n_vars, m) + out = np.empty((flat.shape[0], m), dtype=np.complex128) + for start in range(0, flat.shape[0], chunk_size): + stop = start + chunk_size + out[start:stop] = np.linalg.eigvals(companion_matrix(flat[start:stop], n_lags)) + return out.reshape(*leading, m) + + +def spectral_radius(B: np.ndarray, n_lags: int, *, chunk_size: int = 256) -> np.ndarray: + """Largest companion-matrix eigenvalue modulus, one value per batch element. + + A draw is stable when its spectral radius is strictly below 1 and + explosive at or above it. + + Args: + B: Stacked lag coefficients with trailing shape `(n, n * n_lags)`. + n_lags: Number of lag blocks stacked along the trailing axis. + chunk_size: Forwarded to `companion_eigenvalues`. + + Returns: + Real array with the same shape as *B*'s leading batch axes — a scalar + (0-d) array for a single draw, `(chains, draws)` for a posterior. + + Raises: + ValueError: If the coefficient layout is inconsistent with `n_lags`, + or if `chunk_size` is not positive. + """ + eigenvalues = companion_eigenvalues(B, n_lags, chunk_size=chunk_size) + return np.max(np.abs(eigenvalues), axis=-1) diff --git a/tests/test_stability.py b/tests/test_stability.py new file mode 100644 index 0000000..3ac2c3c --- /dev/null +++ b/tests/test_stability.py @@ -0,0 +1,124 @@ +"""Tests for the companion-form stability primitives.""" + +import numpy as np +import pytest + +from impulso._stability import companion_eigenvalues, companion_matrix, spectral_radius + + +class TestCompanionMatrix: + def test_shape_single_draw(self): + B = np.zeros((3, 6)) + assert companion_matrix(B, n_lags=2).shape == (6, 6) + + def test_shape_batched(self): + B = np.zeros((4, 50, 3, 9)) + assert companion_matrix(B, n_lags=3).shape == (4, 50, 9, 9) + + def test_top_block_is_B_verbatim(self): + rng = np.random.default_rng(0) + B = rng.standard_normal((2, 6)) + F = companion_matrix(B, n_lags=3) + np.testing.assert_array_equal(F[:2, :], B) + + def test_subdiagonal_is_identity(self): + B = np.zeros((2, 6)) + F = companion_matrix(B, n_lags=3) + np.testing.assert_array_equal(F[2:, :4], np.eye(4)) + + def test_trailing_block_column_is_zero_below_top(self): + B = np.ones((2, 6)) + F = companion_matrix(B, n_lags=3) + np.testing.assert_array_equal(F[2:, 4:], np.zeros((4, 2))) + + def test_single_lag_is_B_itself(self): + rng = np.random.default_rng(1) + B = rng.standard_normal((3, 3)) + np.testing.assert_array_equal(companion_matrix(B, n_lags=1), B) + + def test_batch_agnostic(self): + rng = np.random.default_rng(2) + B = rng.standard_normal((2, 5, 3, 6)) + batched = companion_matrix(B, n_lags=2) + for c in range(2): + for d in range(5): + np.testing.assert_array_equal(batched[c, d], companion_matrix(B[c, d], n_lags=2)) + + def test_rejects_zero_lags(self): + with pytest.raises(ValueError, match="n_lags must be positive"): + companion_matrix(np.zeros((2, 4)), n_lags=0) + + def test_rejects_indivisible_trailing_axis(self): + with pytest.raises(ValueError, match="not divisible"): + companion_matrix(np.zeros((2, 5)), n_lags=2) + + def test_rejects_non_matching_row_count(self): + with pytest.raises(ValueError, match="rows"): + companion_matrix(np.zeros((3, 4)), n_lags=2) + + def test_rejects_one_dimensional_input(self): + with pytest.raises(ValueError, match="at least 2 dimensions"): + companion_matrix(np.zeros(4), n_lags=1) + + +class TestCompanionEigenvalues: + def test_known_scalar_ar2_roots(self): + # Scalar AR(2): y_t = 0.5 y_{t-1} + 0.2 y_{t-2}. The companion + # eigenvalues are the roots of lambda^2 - 0.5 lambda - 0.2. + B = np.array([[0.5, 0.2]]) + eigs = companion_eigenvalues(B, n_lags=2) + expected = np.roots([1.0, -0.5, -0.2]) + np.testing.assert_allclose(np.sort(eigs.real), np.sort(expected), atol=1e-12) + np.testing.assert_allclose(np.abs(eigs).max(), 0.7623475, atol=1e-6) + + def test_diagonal_var1_eigenvalues_are_diagonal_entries(self): + B = np.diag([0.9, 0.3]) + eigs = companion_eigenvalues(B, n_lags=1) + np.testing.assert_allclose(np.sort(eigs.real), [0.3, 0.9], atol=1e-12) + + def test_complex_case_returns_complex_dtype_and_real_radius(self): + # A rotation-like block has a conjugate pair of complex eigenvalues. + B = np.array([[0.0, -0.8], [0.8, 0.0]]) + eigs = companion_eigenvalues(B, n_lags=1) + assert np.iscomplexobj(eigs) + assert np.abs(eigs.imag).max() > 0.5 + radius = spectral_radius(B, n_lags=1) + assert not np.iscomplexobj(radius) + np.testing.assert_allclose(radius, 0.8, atol=1e-12) + + def test_chunking_does_not_change_result(self): + rng = np.random.default_rng(3) + B = rng.standard_normal((3, 7, 2, 4)) * 0.3 + small = spectral_radius(B, n_lags=2, chunk_size=1) + large = spectral_radius(B, n_lags=2, chunk_size=10_000) + np.testing.assert_allclose(small, large) + + def test_rejects_non_positive_chunk_size(self): + with pytest.raises(ValueError, match="chunk_size must be positive"): + companion_eigenvalues(np.zeros((2, 2)), n_lags=1, chunk_size=0) + + +class TestSpectralRadius: + def test_matches_naive_loop(self): + rng = np.random.default_rng(4) + B = rng.standard_normal((2, 6, 3, 6)) * 0.2 + vectorised = spectral_radius(B, n_lags=2) + naive = np.zeros((2, 6)) + for c in range(2): + for d in range(6): + F = companion_matrix(B[c, d], n_lags=2) + naive[c, d] = np.abs(np.linalg.eigvals(F)).max() + np.testing.assert_allclose(vectorised, naive) + + def test_batch_shape_preserved(self): + B = np.zeros((5, 11, 2, 4)) + assert spectral_radius(B, n_lags=2).shape == (5, 11) + + def test_single_draw_returns_scalar(self): + radius = spectral_radius(np.diag([0.5, 0.5]), n_lags=1) + assert radius.shape == () + np.testing.assert_allclose(float(radius), 0.5) + + def test_explosive_draw_detected(self): + radius = spectral_radius(np.diag([1.2, 0.5]), n_lags=1) + assert float(radius) > 1.0 From 8506f78ae26e1cf66f09a5de7562cf63471742cf Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 04:01:43 +0200 Subject: [PATCH 2/6] feat(diagnostics): VAR-aware convergence report (#142) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `convergence_report(idata, n_lags=...)` answers what a generic `arviz.summary` cannot: is this VAR posterior usable? It reports R-hat and both effective sample sizes *per parameter block* — coefficient, intercept, exog, covariance, volatility, identification, other — each with the coordinate attaining the worst value, so a mixing problem is attributed to a part of the model rather than to the model as a whole. Divergences are counted once, globally: a divergent transition is a property of a trajectory, not of any single parameter. Two failure modes get named messages with remedies. `rhat_without_ divergences` is the characteristic VAR pathology — near-collinear lag regressors make the posterior ill-conditioned, diagonal mass-matrix adaptation mixes badly across it, and NUTS never has to reject a trajectory to do so, which is why zero divergences is not reassurance here. `explosive_draws` reports posterior mass on non-stationary parameter draws and what it breaks downstream. Explosive draws warn but never fail: mass near a unit root is a legitimate posterior statement about level data under a random-walk prior mean, not evidence the sampler misbehaved. "failed" is reserved for sampler pathology. Block resolution is three-tiered — a static map of Impulso's own variable names, the `v{i}_` prefix on stochastic-volatility latents, then the new optional `posterior_var_names()` capability on `VolatilityProcess` (documented as optional, in the `_samples_rotations` mould, not a protocol requirement). Anything unrecognised lands in `other`; an unknown variable is never an error. Conjugate fits are supported with honest gaps: single chain, so R-hat is None throughout and the report says why; no sampler statistics, so divergences are None at info severity; stability computed in full. The report never calls `warnings.warn` — the returned object is the channel. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Egjd7ToFeb9TQqFnfRQZxV --- src/impulso/__init__.py | 27 + src/impulso/conjugate_volatility.py | 14 + src/impulso/diagnostics.py | 870 ++++++++++++++++++++++++++++ src/impulso/fitted.py | 63 ++ src/impulso/identified.py | 41 ++ src/impulso/protocols.py | 9 + src/impulso/sv/spec.py | 12 + src/impulso/volatility.py | 12 + tests/conftest.py | 124 ++++ tests/test_diagnostics.py | 539 +++++++++++++++++ tests/test_public_api.py | 48 ++ 11 files changed, 1759 insertions(+) create mode 100644 src/impulso/diagnostics.py create mode 100644 tests/test_diagnostics.py diff --git a/src/impulso/__init__.py b/src/impulso/__init__.py index 8fe3f84..8e4619c 100644 --- a/src/impulso/__init__.py +++ b/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", @@ -83,6 +96,7 @@ "ShockPath", "SignRestriction", "StationarityTestResult", + "StabilitySummary", "StochasticVolatility", "StudentT", "VARData", @@ -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, diff --git a/src/impulso/conjugate_volatility.py b/src/impulso/conjugate_volatility.py index f6eb6fe..ffcd206 100644 --- a/src/impulso/conjugate_volatility.py +++ b/src/impulso/conjugate_volatility.py @@ -118,6 +118,20 @@ def _forecast_indices(self, steps: int) -> np.ndarray: """Absolute in-sample-equivalent time indices for forecast steps `0..steps-1`.""" raise NotImplementedError + # --- optional diagnostics capability --- + def posterior_var_names(self) -> tuple[str, ...]: + """Posterior variables this adapter is responsible for. + + The optional `VolatilityProcess` capability read by + `impulso.diagnostics.assign_blocks`. The break's free + hyperparameters *are* its posterior variables — the conjugate + sampler packs one entry per `hyperparameter_priors` key — so the + two stay in sync by construction. They route to the report's + volatility block, since `is_time_varying` is True: they parameterise + the scale path `s_t`, not the base covariance. + """ + return tuple(self.hyperparameter_priors()) + # --- query surface (shared) --- def cholesky_at(self, posterior: xr.Dataset, t: int | None) -> np.ndarray: """Cholesky factor `L_t = s_t * L_base` at time `t` for every draw. diff --git a/src/impulso/diagnostics.py b/src/impulso/diagnostics.py new file mode 100644 index 0000000..d0a4255 --- /dev/null +++ b/src/impulso/diagnostics.py @@ -0,0 +1,870 @@ +"""VAR-aware convergence and stability diagnostics. + +`convergence_report` answers the question a generic `arviz.summary` cannot: +*is this particular VAR posterior usable?* It differs from a plain summary +table in three ways. + +**Parameter blocks.** A VAR posterior mixes quantities with very different +sampling behaviour — lag coefficients, intercepts, the covariance +parameterisation, volatility latents. A single worst-R-hat over all of them +hides which part of the model is failing, so every metric is reported per +block (see `assign_blocks`) with the offending coordinate named. + +**Dynamic stability.** Convergence is necessary but not sufficient: a +perfectly-mixed posterior can put mass on explosive parameter draws, whose +impulse responses diverge with the horizon and whose forecast fans are +unbounded. The report computes the companion-matrix spectral radius of every +draw and reports the explosive fraction alongside the sampling metrics. + +**The VAR failure mode.** Elevated R-hat with *zero* divergences is common in +VARs and unusual elsewhere: near-collinear lag regressors give an +ill-conditioned posterior that a diagonal mass matrix explores badly, and +NUTS never has to reject a trajectory to mix poorly. The report says so +explicitly, with remedies, rather than leaving the user to conclude that no +divergences means no problem. + +The report never calls `warnings.warn`. The returned object *is* the +channel: `status`, `messages`, and the per-block table carry everything, so +callers decide whether to print, raise, or ignore. +""" + +from __future__ import annotations + +import re +import warnings +from typing import TYPE_CHECKING, Literal, Self + +import arviz as az +import numpy as np +import pandas as pd +from pydantic import Field, model_validator + +from impulso._base import ImpulsoBaseModel, ImpulsoModel +from impulso._stability import spectral_radius + +if TYPE_CHECKING: + import xarray as xr + + from impulso.protocols import VolatilityProcess + +# -------------------------------------------------------------------------- +# Block taxonomy +# -------------------------------------------------------------------------- + +# Canonical report order. Blocks with no variables are omitted entirely. +_BLOCK_ORDER: tuple[str, ...] = ( + "coefficient", + "intercept", + "exog", + "covariance", + "volatility", + "identification", + "other", +) + +# Tier 1 of block resolution: posterior variables Impulso itself registers. +_BLOCK_MAP: dict[str, str] = { + "B": "coefficient", + "intercept": "intercept", + "B_exog": "exog", + "sigma_sd": "covariance", + "tril_offdiag": "covariance", + "L": "covariance", + "Sigma": "covariance", + "h": "volatility", + "R_chol": "volatility", + "R_chol_offdiag": "volatility", + "structural_shock_matrix": "identification", + "P": "identification", +} + +# Tier 2: per-variable stochastic-volatility latents are registered with a +# `v{i}_` prefix (see `impulso.sv.spec`), so the whole family maps by pattern. +_SV_PREFIX = re.compile(r"^v\d+_") + +# The only sampler statistic PyMC and nutpie agree on the name of. Every other +# stat (tree depth, acceptance rate, energy) is spelled differently by the two +# backends, so the report reads this one and nothing else. +_DIVERGING_KEY = "diverging" + + +def assign_blocks( + posterior: xr.Dataset, + volatility: VolatilityProcess | None = None, +) -> dict[str, list[str]]: + """Group posterior variables into diagnostic blocks. + + Resolution is three-tiered, first match wins: + + 1. A static map of the variable names Impulso's own estimators register. + 2. The `v{i}_` prefix carried by every per-variable stochastic-volatility + latent, present and future. + 3. The optional `posterior_var_names()` capability on the volatility + process, letting a custom adapter claim the names it registered. + Claimed names join `volatility` if the adapter is time-varying and + `covariance` otherwise. + + Anything left over lands in `other`. Unknown variables are never an + error — a hand-built or third-party posterior still gets a report, and + the block's variable list makes clear what was not recognised. + + Args: + posterior: The posterior Dataset (`idata.posterior`). + volatility: Volatility process used at fit time, consulted for the + optional `posterior_var_names()` capability. Optional. + + Returns: + Mapping from block name to its sorted variable names, in canonical + block order. Blocks with no variables are absent. + """ + claimed: dict[str, str] = {} + hook = getattr(volatility, "posterior_var_names", None) + if hook is not None: + target = "volatility" if getattr(volatility, "is_time_varying", False) else "covariance" + claimed = dict.fromkeys(hook(), target) + + grouped: dict[str, list[str]] = {} + for raw in posterior.data_vars: + name = str(raw) + block = _BLOCK_MAP.get(name) + if block is None and _SV_PREFIX.match(name): + block = "volatility" + if block is None: + block = claimed.get(name, "other") + grouped.setdefault(block, []).append(name) + return {block: sorted(grouped[block]) for block in _BLOCK_ORDER if block in grouped} + + +# -------------------------------------------------------------------------- +# Result objects +# -------------------------------------------------------------------------- + + +class ConvergenceThresholds(ImpulsoModel): + """Cut-offs separating a passing report from warnings and failures. + + Comparisons are strict, so a metric sitting exactly on a threshold + passes: `max_rhat == 1.01` does not warn. + + Attributes: + rhat_warn: R-hat above this warns. Default 1.01, the rank-normalised + split-R-hat cut-off of Vehtari et al. (2021). + rhat_fail: R-hat above this fails. Default 1.05, the classic + Gelman-Rubin rule of thumb. + ess_warn: Effective sample size below this warns. Default 400 — + 100 per chain at the default four chains. + ess_fail: Effective sample size below this fails. Default 100. + divergence_fail_rate: Divergence rate at or above which the report + fails. Default 0.01; any divergence at all warns. + explosive_warn: Fraction of explosive draws at or above which the + explosive-draw message is raised from informational to a + warning. Default 0.05. Explosive draws never fail a report. + """ + + rhat_warn: float = 1.01 + rhat_fail: float = 1.05 + ess_warn: float = 400.0 + ess_fail: float = 100.0 + divergence_fail_rate: float = 0.01 + explosive_warn: float = 0.05 + + +class DiagnosticMessage(ImpulsoModel): + """A single machine-readable finding. + + Attributes: + code: Stable identifier. Prose may be reworded between releases; + the code is the contract programmatic callers match on. + severity: `"info"`, `"warning"`, or `"failure"`. The report's + `status` is the worst severity present. + message: Human-readable explanation, including remedies where + remedies exist. + block: Parameter block the finding concerns, or None if global. + """ + + code: str + severity: Literal["info", "warning", "failure"] + message: str + block: str | None = None + + +class BlockDiagnostics(ImpulsoModel): + """Convergence metrics for one parameter block. + + Every metric is the worst value over the block's coordinates, paired + with a label naming where it occurred (`"B[y2, L1.y1]"`). Metrics are + None when undefined for every coordinate — R-hat from a single chain, + for instance, or a deterministic that is constant across draws. + + Attributes: + block: Block name (see `assign_blocks`). + var_names: Posterior variables in this block. + n_variables: Number of posterior variables in this block. + n_coordinates: Number of scalar coordinates across those variables. + max_rhat: Worst rank-normalised split R-hat, or None. + max_rhat_coord: Coordinate label attaining `max_rhat`, or None. + min_ess_bulk: Smallest bulk effective sample size, or None. + min_ess_bulk_coord: Coordinate label attaining `min_ess_bulk`. + min_ess_tail: Smallest tail effective sample size, or None. + min_ess_tail_coord: Coordinate label attaining `min_ess_tail`. + """ + + block: str + var_names: list[str] + n_variables: int + n_coordinates: int + max_rhat: float | None = None + max_rhat_coord: str | None = None + min_ess_bulk: float | None = None + min_ess_bulk_coord: str | None = None + min_ess_tail: float | None = None + min_ess_tail_coord: str | None = None + + +class StabilitySummary(ImpulsoBaseModel): + """Posterior distribution of the companion-matrix spectral radius. + + A draw is explosive when its spectral radius reaches 1: the implied + system has no stationary solution, its impulse responses grow without + bound in the horizon, and its forecast fan widens indefinitely. Some + explosive mass is normal on level data under a random-walk prior mean, + which is why it never fails a report on its own. + + Attributes: + radius: Read-only spectral radii with shape `(chain, draw)` — after + thinning, if `stability_draws` was used. + p_explosive: Fraction of draws with radius >= 1. + max_radius: Largest radius over all draws. + n_vars: Number of endogenous variables. + n_lags: Lag order. + hdi_prob: Default probability mass for `hdi`, also used by + `to_dataframe`. + thinned_from: Original number of draws per chain when the radii were + computed on a thinned subset, else None. + """ + + radius: np.ndarray = Field(repr=False) + p_explosive: float + max_radius: float + n_vars: int + n_lags: int + hdi_prob: float = 0.89 + thinned_from: int | None = None + + @model_validator(mode="after") + def _make_readonly(self) -> Self: + radius = np.asarray(self.radius).copy() + radius.flags.writeable = False + object.__setattr__(self, "radius", radius) + return self + + def median(self) -> float: + """Posterior median spectral radius.""" + return float(np.median(self.radius)) + + def hdi(self, prob: float | None = None) -> tuple[float, float]: + """Highest-density interval of the spectral radius. + + Args: + prob: Probability mass. Defaults to `hdi_prob` (0.89). + + Returns: + `(lower, upper)` bounds. + """ + # Pooled over chains: the interval is a statement about the posterior, + # not about any one chain, and a flat array sidesteps ArviZ's pending + # reinterpretation of 2-D input as (chain, draw). + pooled = np.asarray(self.radius).reshape(-1) + bounds = az.hdi(pooled, hdi_prob=self.hdi_prob if prob is None else prob) + return float(bounds[0]), float(bounds[1]) + + def to_dataframe(self) -> pd.DataFrame: + """Single-row frame of the stability summary.""" + lower, upper = self.hdi() + return pd.DataFrame( + [ + { + "median_radius": self.median(), + "hdi_lower": lower, + "hdi_upper": upper, + "max_radius": self.max_radius, + "p_explosive": self.p_explosive, + "n_vars": self.n_vars, + "n_lags": self.n_lags, + } + ], + index=pd.Index(["stability"], name="quantity"), + ) + + +class ConvergenceReport(ImpulsoBaseModel): + """VAR-aware convergence and stability diagnostics for one posterior. + + Produced by `convergence_report`, or by the delegating + `FittedVAR.convergence_report` / `IdentifiedVAR.convergence_report`. + + Attributes: + blocks: Per-block metrics in canonical order. + stability: Spectral-radius summary over the posterior draws. + divergences: Number of divergent transitions, or None when the + sampler recorded no statistics. + n_transitions: Total post-warmup transitions, or None. + divergence_rate: `divergences / n_transitions`, or None. + sampler_stats_available: Whether divergence statistics were found. + n_chains: Number of chains in the posterior. + n_draws: Number of post-warmup draws per chain. + thresholds: Thresholds used to derive `status`. + messages: Findings, each with a stable `code`. + status: `"passed"`, `"warnings"`, or `"failed"` — the worst message + severity. `"failed"` is reserved for sampler pathology. + """ + + blocks: list[BlockDiagnostics] + stability: StabilitySummary + divergences: int | None + n_transitions: int | None + divergence_rate: float | None + sampler_stats_available: bool + n_chains: int + n_draws: int + thresholds: ConvergenceThresholds + messages: list[DiagnosticMessage] + status: Literal["passed", "warnings", "failed"] + + @property + def max_rhat(self) -> float | None: + """Worst R-hat across all blocks, or None if undefined everywhere.""" + return _extreme([block.max_rhat for block in self.blocks], "max") + + @property + def min_ess_bulk(self) -> float | None: + """Smallest bulk effective sample size across all blocks.""" + return _extreme([block.min_ess_bulk for block in self.blocks], "min") + + @property + def min_ess_tail(self) -> float | None: + """Smallest tail effective sample size across all blocks.""" + return _extreme([block.min_ess_tail for block in self.blocks], "min") + + def to_dataframe(self) -> pd.DataFrame: + """Per-block metric table, indexed by block in canonical order.""" + rows = [ + { + "n_variables": block.n_variables, + "n_coordinates": block.n_coordinates, + "max_rhat": block.max_rhat, + "max_rhat_coord": block.max_rhat_coord, + "min_ess_bulk": block.min_ess_bulk, + "min_ess_bulk_coord": block.min_ess_bulk_coord, + "min_ess_tail": block.min_ess_tail, + "min_ess_tail_coord": block.min_ess_tail_coord, + } + for block in self.blocks + ] + return pd.DataFrame(rows, index=pd.Index([block.block for block in self.blocks], name="block")) + + def summary(self) -> str: + """Multi-line human-readable rendering of the whole report.""" + divergences = "unavailable" if self.divergences is None else str(self.divergences) + header = [ + f"Convergence report: {self.status.upper()}", + f" {self.n_chains} chains x {self.n_draws} draws | divergences: {divergences}", + ] + table = self.to_dataframe()[["max_rhat", "min_ess_bulk", "min_ess_tail", "max_rhat_coord"]] + lower, upper = self.stability.hdi() + stability = ( + f"Stability: median spectral radius {self.stability.median():.3f} " + f"[{lower:.3f}, {upper:.3f}] | max {self.stability.max_radius:.3f} | " + f"explosive draws {self.stability.p_explosive:.1%}" + ) + lines = [*header, "", table.to_string(), "", stability] + if self.messages: + lines += ["", "Messages:"] + lines += [f" [{msg.severity}] {msg.code}: {msg.message}" for msg in self.messages] + return "\n".join(lines) + + def __repr__(self) -> str: + """One-line status plus headline numbers.""" + rhat = "n/a" if self.max_rhat is None else f"{self.max_rhat:.3f}" + ess = "n/a" if self.min_ess_bulk is None else f"{self.min_ess_bulk:.0f}" + divergences = "n/a" if self.divergences is None else str(self.divergences) + return ( + f"ConvergenceReport(status={self.status!r}, max_rhat={rhat}, " + f"min_ess_bulk={ess}, divergences={divergences}, " + f"p_explosive={self.stability.p_explosive:.3f})" + ) + + +# -------------------------------------------------------------------------- +# Metric computation +# -------------------------------------------------------------------------- + + +def _extreme(values: list[float | None], mode: Literal["max", "min"]) -> float | None: + """Worst of the non-None entries, or None when every entry is None.""" + present = [value for value in values if value is not None] + if not present: + return None + return max(present) if mode == "max" else min(present) + + +def _fallback_coords(var_names: list[str] | None, n_lags: int) -> dict[str, list[str]]: + """Labels for dimensions a hand-built or conjugate posterior leaves bare. + + `VAR.fit` stamps coords on the posterior, so its labels are read + directly. `ConjugateVAR` builds its Dataset from dims alone, and this + supplies the same names so both estimators produce identical coordinate + labels for the same model. + """ + if not var_names: + return {} + names = list(var_names) + coeff = [f"L{lag}.{name}" for lag in range(1, n_lags + 1) for name in names] + return {"var": names, "var1": names, "var2": names, "variable": names, "response": names, "coeff": coeff} + + +def _worst_in_variable( + name: str, + da: xr.DataArray, + mode: Literal["max", "min"], + fallback: dict[str, list[str]], +) -> tuple[float, str] | None: + """Extreme value and its coordinate label within one variable's metrics. + + Returns None when the metric is NaN at every coordinate — a single-chain + R-hat, or a deterministic that never varies across draws. + """ + values = np.asarray(da.values, dtype=float) + if values.ndim == 0: + scalar = float(values) + return None if np.isnan(scalar) else (scalar, name) + flat = values.reshape(-1) + if bool(np.all(np.isnan(flat))): + return None + index = int(np.nanargmax(flat)) if mode == "max" else int(np.nanargmin(flat)) + position = np.unravel_index(index, values.shape) + labels = [] + for axis, (dim, i) in enumerate(zip(da.dims, position, strict=True)): + if dim in da.coords: + labels.append(str(da.coords[dim].values[i])) + elif dim in fallback and len(fallback[dim]) == values.shape[axis]: + labels.append(fallback[dim][i]) + else: + labels.append(str(i)) + return float(flat[index]), f"{name}[{', '.join(labels)}]" + + +def _worst_across( + metrics: xr.Dataset, + names: list[str], + mode: Literal["max", "min"], + fallback: dict[str, list[str]], +) -> tuple[float | None, str | None]: + """Extreme metric value and label across every variable in a block.""" + best: tuple[float, str] | None = None + for name in names: + found = _worst_in_variable(name, metrics[name], mode, fallback) + if found is None: + continue + if best is None or (found[0] > best[0] if mode == "max" else found[0] < best[0]): + best = found + return best if best is not None else (None, None) + + +def _block_metrics( + posterior: xr.Dataset, + block: str, + names: list[str], + fallback: dict[str, list[str]], +) -> BlockDiagnostics: + """Compute R-hat and both effective sample sizes for one block.""" + subset = posterior[names] + # A deterministic that is constant across draws divides by a zero + # variance inside ArviZ and yields NaN, which is the honest answer and is + # handled downstream. Silence the numpy notice so a report never emits a + # warning of its own. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + rhat = az.rhat(subset) + ess_bulk = az.ess(subset, method="bulk") + ess_tail = az.ess(subset, method="tail") + + max_rhat, max_rhat_coord = _worst_across(rhat, names, "max", fallback) + min_bulk, min_bulk_coord = _worst_across(ess_bulk, names, "min", fallback) + min_tail, min_tail_coord = _worst_across(ess_tail, names, "min", fallback) + n_coordinates = sum(int(np.prod(rhat[name].shape, dtype=int)) for name in names) + return BlockDiagnostics( + block=block, + var_names=list(names), + n_variables=len(names), + n_coordinates=n_coordinates, + max_rhat=max_rhat, + max_rhat_coord=max_rhat_coord, + min_ess_bulk=min_bulk, + min_ess_bulk_coord=min_bulk_coord, + min_ess_tail=min_tail, + min_ess_tail_coord=min_tail_coord, + ) + + +def _divergences(idata: az.InferenceData) -> tuple[int | None, int | None, float | None, bool]: + """Global divergence counts read from `sample_stats` only. + + Divergences are a property of a *trajectory*, not of any one parameter, + so they are never attributed to a block. Only the post-warmup group is + read: nutpie also emits `warmup_sample_stats`, whose divergences belong + to adaptation and say nothing about the retained draws. + """ + if "sample_stats" not in idata.groups() or _DIVERGING_KEY not in idata.sample_stats: + return None, None, None, False + diverging = np.asarray(idata.sample_stats[_DIVERGING_KEY].values) + count = int(diverging.sum()) + total = int(diverging.size) + return count, total, (count / total if total else 0.0), True + + +# -------------------------------------------------------------------------- +# Messages and status +# -------------------------------------------------------------------------- + +_NUTPIE_REMEDY = 'NUTSSampler(nuts_sampler="nutpie", nuts_sampler_kwargs={"low_rank_modified_mass_matrix": True})' + + +def _rhat_messages( + report_max_rhat: float | None, + coord: str | None, + divergences: int | None, + stats_available: bool, + thresholds: ConvergenceThresholds, +) -> list[DiagnosticMessage]: + """R-hat findings, including the VAR-specific zero-divergence case.""" + if report_max_rhat is None or report_max_rhat <= thresholds.rhat_warn: + return [] + where = f" (worst: {coord})" if coord else "" + if report_max_rhat > thresholds.rhat_fail: + messages = [ + DiagnosticMessage( + code="rhat_high", + severity="failure", + message=( + f"R-hat reaches {report_max_rhat:.3f}{where}, above the failure " + f"threshold of {thresholds.rhat_fail}. The chains have not mixed; " + "the posterior draws do not describe a single distribution and no " + "downstream quantity should be reported." + ), + ) + ] + else: + messages = [ + DiagnosticMessage( + code="rhat_elevated", + severity="warning", + message=( + f"R-hat reaches {report_max_rhat:.3f}{where}, above the warning " + f"threshold of {thresholds.rhat_warn}. Mixing is imperfect; run " + "longer chains before trusting tail quantities." + ), + ) + ] + if stats_available and divergences == 0: + messages.append( + DiagnosticMessage( + code="rhat_without_divergences", + severity="warning", + message=( + f"R-hat up to {report_max_rhat:.3f} with zero divergences is common in " + "VARs: near-collinear lag regressors make the posterior ill-conditioned " + "and diagonal mass-matrix adaptation mixes poorly across it, without the " + "sampler ever having to reject a trajectory. Divergences are not a " + "reliable alarm here. Remedies, in order: switch to nutpie's low-rank " + f"mass matrix — {_NUTPIE_REMEDY} — then lengthen `tune`, then tighten the " + "Minnesota prior (smaller `tightness`), then reduce the number of " + "variables or lags." + ), + ) + ) + return messages + + +def _ess_messages( + metric: float | None, + coord: str | None, + kind: str, + thresholds: ConvergenceThresholds, +) -> list[DiagnosticMessage]: + """Effective-sample-size findings for one flavour of ESS.""" + if metric is None or metric >= thresholds.ess_warn: + return [] + where = f" (worst: {coord})" if coord else "" + if metric < thresholds.ess_fail: + return [ + DiagnosticMessage( + code=f"ess_{kind}_low", + severity="failure", + message=( + f"{kind.capitalize()} effective sample size falls to {metric:.0f}{where}, " + f"below the failure threshold of {thresholds.ess_fail:.0f}. Posterior " + "summaries carry more Monte Carlo error than signal." + ), + ) + ] + return [ + DiagnosticMessage( + code=f"ess_{kind}_marginal", + severity="warning", + message=( + f"{kind.capitalize()} effective sample size falls to {metric:.0f}{where}, " + f"below the warning threshold of {thresholds.ess_warn:.0f}. Draw more " + "samples before reporting intervals from this block." + ), + ) + ] + + +def _divergence_messages( + divergences: int | None, + n_transitions: int | None, + rate: float | None, + stats_available: bool, + thresholds: ConvergenceThresholds, +) -> list[DiagnosticMessage]: + """Divergence findings, or the informational note when none were recorded.""" + if not stats_available: + return [ + DiagnosticMessage( + code="sampler_stats_missing", + severity="info", + message=( + "No sampler statistics were found, so divergences could not be counted. " + "This is expected for `ConjugateVAR`, which draws coefficients in closed " + "form and has no trajectories to diverge, and for hand-built posteriors. " + "It does not by itself indicate a problem." + ), + ) + ] + if not divergences: + return [] + # `stats_available` guarantees both counts are present; default defensively + # so the message code path never depends on a nullable arithmetic operand. + observed_rate = rate if rate is not None else 0.0 + total = n_transitions if n_transitions is not None else 0 + severity = "failure" if observed_rate >= thresholds.divergence_fail_rate else "warning" + return [ + DiagnosticMessage( + code="divergences_present", + severity=severity, + message=( + f"{divergences} of {total} transitions diverged ({observed_rate:.2%}). " + "Divergent trajectories mean the sampler could not follow the posterior's " + "geometry, so the draws are biased toward the regions it could reach. " + "Raise `target_accept` toward 0.95, lengthen `tune`, or tighten the prior." + ), + ) + ] + + +def _chain_messages(n_chains: int) -> list[DiagnosticMessage]: + """The single-chain note: R-hat is a between-chain statistic.""" + if n_chains >= 2: + return [] + return [ + DiagnosticMessage( + code="single_chain", + severity="warning", + message=( + "Only one chain is present, so R-hat is undefined and is reported as None " + "for every block; effective sample size is still computed. This is the " + "normal shape of a `ConjugateVAR` posterior, where the coefficient and " + "Cholesky draws are exact conditional draws (their effective sample size " + "is nominal by construction) and only the hyperparameters — sampled by " + "random-walk Metropolis — carry meaningful autocorrelation. For a NUTS " + "fit, sample at least two chains before trusting any convergence claim." + ), + ) + ] + + +def _stability_messages( + stability: StabilitySummary, + thresholds: ConvergenceThresholds, +) -> list[DiagnosticMessage]: + """The explosive-draw finding. Never a failure — see ADR-0008.""" + if stability.p_explosive <= 0: + return [] + severity = "warning" if stability.p_explosive >= thresholds.explosive_warn else "info" + return [ + DiagnosticMessage( + code="explosive_draws", + severity=severity, + block="coefficient", + message=( + f"{stability.p_explosive:.1%} of draws are explosive (companion-matrix " + f"spectral radius >= 1; largest {stability.max_radius:.3f}). Their impulse " + "responses diverge with the horizon, their forecast fans are unbounded, " + "their long-horizon FEVD shares are uninterpretable, and their historical " + "decomposition baselines drift. This is a property of the model, not of the " + "sampler, and near-unit-root mass is legitimate on level data under a " + "random-walk prior mean. If it is not intended: difference the data (or " + "test for cointegration), tighten the shrinkage, or restrict reported " + "horizons to where the responses are still meaningful." + ), + ) + ] + + +def _derive_status(messages: list[DiagnosticMessage]) -> Literal["passed", "warnings", "failed"]: + """Worst severity present. Only sampler pathology reaches `"failed"`.""" + severities = {message.severity for message in messages} + if "failure" in severities: + return "failed" + if "warning" in severities: + return "warnings" + return "passed" + + +# -------------------------------------------------------------------------- +# Entry point +# -------------------------------------------------------------------------- + + +def _stability_summary( + posterior: xr.Dataset, + n_lags: int, + hdi_prob: float, + stability_draws: int | None, +) -> StabilitySummary: + """Spectral-radius summary over the posterior's lag coefficients.""" + B_da = posterior["B"] + if set(B_da.dims) == {"chain", "draw", "var", "coeff"}: + B_da = B_da.transpose("chain", "draw", "var", "coeff") + B = np.asarray(B_da.values, dtype=float) + + thinned_from: int | None = None + if stability_draws is not None: + if stability_draws < 1: + raise ValueError(f"stability_draws must be positive, got {stability_draws}") + n_draws = B.shape[1] + if stability_draws < n_draws: + # Deterministic stride, never an RNG: two calls on one posterior + # must return identical numbers. + stride = -(-n_draws // stability_draws) + B = B[:, ::stride] + thinned_from = n_draws + + radius = spectral_radius(B, n_lags) + return StabilitySummary( + radius=radius, + p_explosive=float(np.mean(radius >= 1.0)), + max_radius=float(np.max(radius)), + n_vars=B.shape[-2], + n_lags=n_lags, + hdi_prob=hdi_prob, + thinned_from=thinned_from, + ) + + +def convergence_report( + idata: az.InferenceData, + *, + n_lags: int, + var_names: list[str] | None = None, + volatility: VolatilityProcess | None = None, + thresholds: ConvergenceThresholds | None = None, + hdi_prob: float = 0.89, + stability_draws: int | None = None, +) -> ConvergenceReport: + """Build a VAR-aware convergence and stability report for a posterior. + + Reports R-hat and both effective sample sizes per parameter block, each + with the coordinate that attains the worst value; the global divergence + count; and the posterior distribution of the companion-matrix spectral + radius. See the module docstring for why a VAR needs its own report. + + `"failed"` is reserved for sampler pathology — R-hat above + `rhat_fail`, effective sample size below `ess_fail`, or a divergence + rate at or above `divergence_fail_rate`. Explosive draws warn but never + fail: mass near a unit root is a legitimate posterior statement about + level data, not evidence that the sampler misbehaved. + + Cost is dominated by the eigendecomposition of one `(n * p, n * p)` + companion matrix per draw. Pass `stability_draws` to compute the radii + on a deterministically strided subset when the posterior is large. + + Args: + idata: Posterior to diagnose. Must carry a `posterior` group with + reduced-form lag coefficients `B`; `sample_stats` is read when + present and its absence is reported, not an error. + n_lags: Lag order, needed to build the companion matrix. + var_names: Endogenous variable names, used to label coordinates on + posteriors that carry no coords of their own (`ConjugateVAR` + builds one such). Optional. + volatility: Volatility process used at fit time. Consulted for the + optional `posterior_var_names()` capability when assigning + blocks. Optional. + thresholds: Cut-offs driving `status`. Defaults to + `ConvergenceThresholds()`. + hdi_prob: Default probability mass for the spectral radius interval. + stability_draws: Approximate number of draws per chain to retain for + the spectral-radius computation. Thinning is a deterministic + stride, so repeated calls agree exactly. + + Returns: + A `ConvergenceReport`. Nothing is warned or raised on a bad + posterior — the report is the channel. + + Raises: + ValueError: If `idata` has no `posterior` group, if the posterior + carries no `B` (a univariate `FittedSV` posterior has no lag + coefficients and is not supported), or if `stability_draws` is + not positive. + """ + if "posterior" not in idata.groups(): + raise ValueError("convergence_report requires an InferenceData with a `posterior` group.") + posterior = idata.posterior + if "B" not in posterior: + raise ValueError( + "convergence_report requires reduced-form VAR lag coefficients `B` in the " + "posterior, and this one has none. A univariate stochastic-volatility fit " + "(`FittedSV`) has no lag coefficients and therefore no companion matrix, so " + "it is not supported; diagnose it with ArviZ directly." + ) + + thresholds = thresholds or ConvergenceThresholds() + fallback = _fallback_coords(var_names, n_lags) + blocks = [ + _block_metrics(posterior, block, names, fallback) + for block, names in assign_blocks(posterior, volatility).items() + ] + stability = _stability_summary(posterior, n_lags, hdi_prob, stability_draws) + divergences, n_transitions, rate, stats_available = _divergences(idata) + + max_rhat = _extreme([block.max_rhat for block in blocks], "max") + rhat_coord = next((block.max_rhat_coord for block in blocks if block.max_rhat == max_rhat), None) + min_bulk = _extreme([block.min_ess_bulk for block in blocks], "min") + bulk_coord = next((block.min_ess_bulk_coord for block in blocks if block.min_ess_bulk == min_bulk), None) + min_tail = _extreme([block.min_ess_tail for block in blocks], "min") + tail_coord = next((block.min_ess_tail_coord for block in blocks if block.min_ess_tail == min_tail), None) + + n_chains = int(posterior.sizes["chain"]) + messages = [ + *_rhat_messages(max_rhat, rhat_coord, divergences, stats_available, thresholds), + *_ess_messages(min_bulk, bulk_coord, "bulk", thresholds), + *_ess_messages(min_tail, tail_coord, "tail", thresholds), + *_divergence_messages(divergences, n_transitions, rate, stats_available, thresholds), + *_chain_messages(n_chains), + *_stability_messages(stability, thresholds), + ] + + return ConvergenceReport( + blocks=blocks, + stability=stability, + divergences=divergences, + n_transitions=n_transitions, + divergence_rate=rate, + sampler_stats_available=stats_available, + n_chains=n_chains, + n_draws=int(posterior.sizes["draw"]), + thresholds=thresholds, + messages=messages, + status=_derive_status(messages), + ) diff --git a/src/impulso/fitted.py b/src/impulso/fitted.py index f329a61..c38bd7a 100644 --- a/src/impulso/fitted.py +++ b/src/impulso/fitted.py @@ -15,6 +15,7 @@ from impulso.protocols import ErrorDistribution, IdentificationScheme, VolatilityProcess if TYPE_CHECKING: + from impulso.diagnostics import ConvergenceReport, ConvergenceThresholds from impulso.identified import IdentifiedVAR from impulso.results import ConditionalForecastResult, DynamicMultiplierResult, ForecastResult from impulso.scenario import VariablePath @@ -507,6 +508,68 @@ def dynamic_multiplier(self, horizon: int = 20, cumulative: bool = False) -> "Dy cumulative=cumulative, ) + def convergence_report( + self, + thresholds: "ConvergenceThresholds | None" = None, + hdi_prob: float = 0.89, + stability_draws: int | None = None, + ) -> "ConvergenceReport": + """Diagnose this posterior's convergence and dynamic stability. + + Reports R-hat and both effective sample sizes *per parameter block* + — lag coefficients, intercept, exogenous coefficients, covariance, + volatility latents — each with the coordinate attaining the worst + value, so a mixing problem is attributed to a part of the model + rather than to the model as a whole. Divergences are reported once, + globally: a divergent transition is a property of a trajectory, not + of any single parameter. + + Alongside the sampling metrics the report computes the + companion-matrix spectral radius of every draw, built from the + posterior's `B` on the assumption Impulso's own estimators satisfy + — lag blocks concatenated in lag order along the trailing axis. + Draws whose radius reaches 1 are explosive; the report says how many + there are and what breaks because of them, but explosive draws never + fail a report, because mass near a unit root is a legitimate + posterior statement about level data. + + Two failure modes get named messages with remedies: + `rhat_without_divergences` (the characteristic VAR pathology — an + ill-conditioned posterior from near-collinear lag regressors that + NUTS mixes badly across without ever diverging) and + `explosive_draws`. + + A `ConjugateVAR` fit is supported with honest gaps: its posterior + has a single chain, so R-hat is None throughout and the report says + why, and it carries no sampler statistics, so divergences are None. + Stability is computed in full. + + Args: + thresholds: Cut-offs driving the report's status. Defaults to + `ConvergenceThresholds()`. + hdi_prob: Default probability mass for the spectral-radius + interval. + stability_draws: Approximate number of draws per chain to retain + for the spectral-radius computation, thinned by a + deterministic stride. Use on large posteriors, where one + eigendecomposition per draw dominates the cost. + + Returns: + A `ConvergenceReport`. Nothing is warned or raised for a bad + posterior — the returned object is the channel. + """ + from impulso.diagnostics import convergence_report + + return convergence_report( + self.idata, + n_lags=self.n_lags, + var_names=self.var_names, + volatility=self.volatility, + thresholds=thresholds, + hdi_prob=hdi_prob, + stability_draws=stability_draws, + ) + def set_identification_strategy(self, scheme: IdentificationScheme) -> "IdentifiedVAR": """Apply a structural identification scheme. diff --git a/src/impulso/identified.py b/src/impulso/identified.py index c4d9283..fa80d51 100644 --- a/src/impulso/identified.py +++ b/src/impulso/identified.py @@ -24,6 +24,7 @@ ) if TYPE_CHECKING: + from impulso.diagnostics import ConvergenceReport, ConvergenceThresholds from impulso.scenario import ShockPath, VariablePath # Type alias for the `at=` parameter used by query methods. @@ -764,3 +765,43 @@ def structural_scenario( adjusting=adjusting if adjusting is None else list(adjusting), shocks=list(shocks or []), ) + + def convergence_report( + self, + thresholds: "ConvergenceThresholds | None" = None, + hdi_prob: float = 0.89, + stability_draws: int | None = None, + ) -> "ConvergenceReport": + """Diagnose the underlying reduced-form posterior. + + Identical to `FittedVAR.convergence_report` — an `IdentifiedVAR` + shares its `idata` with the `FittedVAR` it came from, and + identification adds nothing that needs diagnosing. The structural + shock matrix is deliberately excluded: under `Cholesky` it is a + deterministic function of draws already diagnosed in the covariance + block, and under `SignRestriction` it is resampled per call, so its + R-hat would describe the rotation sampler rather than the posterior. + + Args: + thresholds: Cut-offs driving the report's status. Defaults to + `ConvergenceThresholds()`. + hdi_prob: Default probability mass for the spectral-radius + interval. + stability_draws: Approximate number of draws per chain to retain + for the spectral-radius computation, thinned by a + deterministic stride. + + Returns: + A `ConvergenceReport` for the reduced-form posterior. + """ + from impulso.diagnostics import convergence_report + + return convergence_report( + self.idata, + n_lags=self.n_lags, + var_names=self.var_names, + volatility=self.volatility, + thresholds=thresholds, + hdi_prob=hdi_prob, + stability_draws=stability_draws, + ) diff --git a/src/impulso/protocols.py b/src/impulso/protocols.py index e6a99dc..a1b25b2 100644 --- a/src/impulso/protocols.py +++ b/src/impulso/protocols.py @@ -212,6 +212,15 @@ class VolatilityProcess(Protocol): Adapters own their downstream computation: time-`t` query and forward simulation for forecasts. + + Optional capability: `posterior_var_names(self) -> tuple[str, ...]`, + naming the posterior variables the adapter is responsible for. It is + *not* a required protocol method — `impulso.diagnostics.assign_blocks` + reads it via `getattr(volatility, "posterior_var_names", None)` (the + same pattern as `IdentificationScheme._samples_rotations`) to route an + adapter's own variables into the report's covariance or volatility + block. Adapters that omit it lose nothing beyond block attribution: + unrecognised variables land in the report's `other` block. """ name: str diff --git a/src/impulso/sv/spec.py b/src/impulso/sv/spec.py index faa69b7..7390a3b 100644 --- a/src/impulso/sv/spec.py +++ b/src/impulso/sv/spec.py @@ -228,6 +228,18 @@ def _clark_reconstruct(h: np.ndarray, R_chol: np.ndarray) -> np.ndarray: R = R_chol.shape[:-2] + (1,) * extra + R_chol.shape[-2:] return np.exp(h / 2)[..., :, np.newaxis] * R_chol.reshape(R) + def posterior_var_names(self) -> tuple[str, ...]: + """Posterior variables this adapter is responsible for. + + The optional `VolatilityProcess` capability read by + `impulso.diagnostics.assign_blocks`. Only the shared quantities are + named: the per-variable latents (`v0_h`, `v0_sigma_eta`, …) carry + the `v{i}_` prefix that block assignment already recognises by + pattern, so listing them here would duplicate that rule and go + stale whenever a new dynamics adapter adds a parameter. + """ + return ("h", "R_chol", "R_chol_offdiag") + def cholesky_at(self, posterior: "xr.Dataset", t: int | None) -> np.ndarray: """Return L_t = diag(exp(h_t / 2)) @ R_chol for the requested t. diff --git a/src/impulso/volatility.py b/src/impulso/volatility.py index fa941b2..e1afd2c 100644 --- a/src/impulso/volatility.py +++ b/src/impulso/volatility.py @@ -91,6 +91,18 @@ def build_pymc_latent( # from the posterior instead of re-decomposing Σ on every call. return pm.Deterministic("L", L) + def posterior_var_names(self) -> tuple[str, ...]: + """Posterior variables this adapter is responsible for. + + The optional `VolatilityProcess` capability read by + `impulso.diagnostics.assign_blocks`. `Sigma` appears here even + though `spec.py` registers it (the deterministic belongs to the + constant parameterisation), and `tril_offdiag` is absent from the + posterior when `n_vars == 1` — claiming a name the posterior lacks + is harmless, since block assignment iterates over what is there. + """ + return ("sigma_sd", "tril_offdiag", "L", "Sigma") + def cholesky_at(self, posterior: "xr.Dataset", t: int | None) -> np.ndarray: """Return the lower-triangular Cholesky factor of Σ for every draw. diff --git a/tests/conftest.py b/tests/conftest.py index 0b3f523..b78780c 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -191,6 +191,130 @@ def permanent_transitory_2v(): "L": L_draws, "idata": az.InferenceData(posterior=posterior), } +# --------------- Diagnostics posterior factory --------------- + +# PyMC and nutpie spell almost every sampler statistic differently; only +# `diverging` is common. The factory can emit either shape so tests can prove +# the report reads that one key and ignores the rest. +_PYMC_EXTRA_STATS = ("reached_max_treedepth", "tree_depth", "acceptance_rate", "lp", "energy", "step_size") +_NUTPIE_EXTRA_STATS = ("maxdepth_reached", "depth", "mean_tree_accept", "logp", "energy", "step_size", "tuning") + + +def _coefficient_draws(rng, shape, explosive_frac, bad_coord): + """Lag-coefficient draws centred on `0.5 * I`, with optional pathologies.""" + n_chains, n_draws, n_vars, n_coeff = shape + base = np.zeros((n_vars, n_coeff)) + base[:, :n_vars] = 0.5 * np.eye(n_vars) + if explosive_frac > 0: + # Noiseless two-point mixture: every radius is exactly 0.5 or 1.2, so + # the reported statistics are exact rather than approximate. + B = np.broadcast_to(base, shape).copy() + n_explosive = round(explosive_frac * n_draws) + for chain in range(n_chains): + B[chain, rng.permutation(n_draws)[:n_explosive], 0, 0] = 1.2 + else: + B = base + 0.02 * rng.standard_normal(shape) + if bad_coord is not None: + row, col = bad_coord + B[:, :, row, col] += np.arange(n_chains)[:, None] * 2.0 + return B + + +def _cholesky_draws(rng, sigma_sd, n_chains, n_draws, n_vars): + """Lower-triangular Cholesky draws; the upper triangle is a structural zero.""" + L = np.zeros((n_chains, n_draws, n_vars, n_vars)) + diag = np.arange(n_vars) + L[..., diag, diag] = sigma_sd + for i in range(1, n_vars): + for j in range(i): + L[..., i, j] = 0.1 * rng.standard_normal((n_chains, n_draws)) + return L + + +def _sampler_stats(rng, n_chains, n_draws, divergences, nutpie_shaped): + """A `sample_stats` group in either backend's shape, plus nutpie's warmup group.""" + flat = np.zeros(n_chains * n_draws, dtype=bool) + flat[:divergences] = True + stats = {"diverging": (("chain", "draw"), flat.reshape(n_chains, n_draws))} + for name in _NUTPIE_EXTRA_STATS if nutpie_shaped else _PYMC_EXTRA_STATS: + stats[name] = (("chain", "draw"), rng.standard_normal((n_chains, n_draws))) + groups = {"sample_stats": xr.Dataset(stats)} + if nutpie_shaped: + # Warmup divergences belong to adaptation and must not be counted. + groups["warmup_sample_stats"] = xr.Dataset({ + "diverging": (("chain", "draw"), np.ones((n_chains, n_draws), dtype=bool)) + }) + return groups + + +@pytest.fixture +def make_var_posterior(): + """Factory building synthetic VAR posteriors with controlled pathologies. + + Returns a callable accepting: + + * `n_chains`, `n_draws`, `n_vars`, `n_lags` — posterior shape. + * `bad_coord` — `(row, col)` index into `B`; that coordinate gets a + chain-dependent offset, so R-hat blows up while every other + coordinate stays healthy. + * `explosive_frac` — fraction of draws per chain whose lag block is + `diag(1.2, 0.5, ...)` instead of `diag(0.5, ...)`. Switches `B` to a + noiseless two-point mixture so the spectral-radius statistics are + exact, with the explosive draws scattered within each chain so the + indicator is neither autocorrelated nor chain-dependent. + * `divergences` — divergent-transition count, placed at fixed positions. + `None` omits the `sample_stats` group entirely. + * `extra_vars` — mapping of posterior variable name to trailing shape. + * `coords` — attach `var`/`coeff` coords (as `VAR.fit` does) or leave + the dims bare (as `ConjugateVAR` does). + * `nutpie_shaped` — emit nutpie's sampler-stat names plus a + `warmup_sample_stats` group full of divergences that must be ignored. + """ + + def _make( + *, + n_chains=4, + n_draws=200, + n_vars=2, + n_lags=1, + bad_coord=None, + explosive_frac=0.0, + divergences=0, + extra_vars=None, + coords=True, + nutpie_shaped=False, + seed=0, + ): + rng = np.random.default_rng(seed) + n_coeff = n_vars * n_lags + var_names = [f"y{i + 1}" for i in range(n_vars)] + coeff_names = [f"L{lag}.{name}" for lag in range(1, n_lags + 1) for name in var_names] + + B = _coefficient_draws(rng, (n_chains, n_draws, n_vars, n_coeff), explosive_frac, bad_coord) + sigma_sd = 0.5 + 0.05 * np.abs(rng.standard_normal((n_chains, n_draws, n_vars))) + L = _cholesky_draws(rng, sigma_sd, n_chains, n_draws, n_vars) + + data_vars = { + "B": (("chain", "draw", "var", "coeff"), B), + "intercept": (("chain", "draw", "var"), 0.01 * rng.standard_normal((n_chains, n_draws, n_vars))), + "sigma_sd": (("chain", "draw", "var"), sigma_sd), + "L": (("chain", "draw", "var1", "var2"), L), + } + for name, shape in (extra_vars or {}).items(): + dims = tuple(f"{name}_dim_{k}" for k in range(len(shape))) + data_vars[name] = (("chain", "draw", *dims), rng.standard_normal((n_chains, n_draws, *shape))) + + posterior_coords = ( + {"var": var_names, "coeff": coeff_names, "var1": var_names, "var2": var_names} if coords else None + ) + groups = {"posterior": xr.Dataset(data_vars, coords=posterior_coords)} + + if divergences is not None: + groups |= _sampler_stats(rng, n_chains, n_draws, divergences, nutpie_shaped) + + return az.InferenceData(**groups) + + return _make # --------------- SV fixtures --------------- diff --git a/tests/test_diagnostics.py b/tests/test_diagnostics.py new file mode 100644 index 0000000..f3f1065 --- /dev/null +++ b/tests/test_diagnostics.py @@ -0,0 +1,539 @@ +"""Tests for the VAR-aware convergence and stability report.""" + +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from arviz import InferenceData + +from impulso.conjugate import ConjugateVAR +from impulso.conjugate_volatility import PandemicBreak +from impulso.diagnostics import ( + BlockDiagnostics, + ConvergenceReport, + ConvergenceThresholds, + StabilitySummary, + assign_blocks, + convergence_report, +) +from impulso.fitted import FittedVAR +from impulso.identification import Cholesky +from impulso.priors import NIWPrior +from impulso.sv.spec import StochasticVolatility +from impulso.volatility import Constant + + +def codes(report: ConvergenceReport) -> list[str]: + return [message.code for message in report.messages] + + +def blocks_by_name(report: ConvergenceReport) -> dict[str, BlockDiagnostics]: + return {block.block: block for block in report.blocks} + + +class TestHealthyPosterior: + def test_status_passed_with_no_messages(self, make_var_posterior): + report = convergence_report(make_var_posterior(), n_lags=1, var_names=["y1", "y2"]) + assert report.status == "passed" + assert report.messages == [] + + def test_headline_metrics_are_sane(self, make_var_posterior): + report = convergence_report(make_var_posterior(), n_lags=1, var_names=["y1", "y2"]) + assert report.max_rhat is not None + assert report.max_rhat < 1.01 + assert report.min_ess_bulk is not None + assert report.min_ess_bulk > 400 + assert report.n_chains == 4 + assert report.n_draws == 200 + + def test_stability_is_centred_on_half(self, make_var_posterior): + report = convergence_report(make_var_posterior(), n_lags=1) + assert report.stability.p_explosive == 0.0 + assert report.stability.median() == pytest.approx(0.5, abs=0.05) + assert report.stability.n_vars == 2 + assert report.stability.n_lags == 1 + + def test_blocks_present_in_canonical_order(self, make_var_posterior): + report = convergence_report(make_var_posterior(), n_lags=1) + assert [block.block for block in report.blocks] == ["coefficient", "intercept", "covariance"] + + def test_structural_zeros_do_not_break_the_covariance_block(self, make_var_posterior): + # `L` is lower triangular, so L[y1, y2] is constant at zero and its + # R-hat is NaN. The block still reports a finite worst value. + block = blocks_by_name(convergence_report(make_var_posterior(), n_lags=1))["covariance"] + assert block.max_rhat is not None + assert np.isfinite(block.max_rhat) + assert block.n_coordinates == 6 # sigma_sd (2) + L (2x2) + + def test_report_emits_no_python_warnings(self, make_var_posterior): + # The report object is the channel; it never uses `warnings.warn`, + # and it swallows the divide-by-zero notice ArviZ raises on the + # structural zeros in `L`. + import warnings + + idata = make_var_posterior() + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + convergence_report(idata, n_lags=1) + assert caught == [] + + +class TestNonConvergedPosterior: + @pytest.fixture + def report(self, make_var_posterior): + idata = make_var_posterior(bad_coord=(1, 0), divergences=0) + return convergence_report(idata, n_lags=1, var_names=["y1", "y2"]) + + def test_status_failed(self, report): + assert report.status == "failed" + assert report.max_rhat > 1.05 + + def test_rhat_without_divergences_code_present(self, report): + assert "rhat_without_divergences" in codes(report) + + def test_message_quotes_the_low_rank_mass_matrix_remedy(self, report): + message = next(m for m in report.messages if m.code == "rhat_without_divergences") + assert "low_rank_modified_mass_matrix" in message.message + assert 'nuts_sampler="nutpie"' in message.message + + def test_worst_coordinate_named_with_coords(self, report): + assert report.blocks[0].max_rhat_coord == "B[y2, L1.y1]" + + def test_worst_coordinate_named_without_coords(self, make_var_posterior): + idata = make_var_posterior(bad_coord=(1, 0), coords=False) + report = convergence_report(idata, n_lags=1, var_names=["y1", "y2"]) + assert report.blocks[0].max_rhat_coord == "B[y2, L1.y1]" + + def test_worst_coordinate_falls_back_to_positions(self, make_var_posterior): + idata = make_var_posterior(bad_coord=(1, 0), coords=False) + report = convergence_report(idata, n_lags=1, var_names=None) + assert report.blocks[0].max_rhat_coord == "B[1, 0]" + + def test_healthy_blocks_are_not_implicated(self, report): + intercept_rhat = blocks_by_name(report)["intercept"].max_rhat + assert intercept_rhat is not None + assert intercept_rhat < 1.01 + + def test_divergences_present_suppresses_the_zero_divergence_message(self, make_var_posterior): + idata = make_var_posterior(bad_coord=(1, 0), divergences=40) + report = convergence_report(idata, n_lags=1) + assert "rhat_without_divergences" not in codes(report) + assert "divergences_present" in codes(report) + assert report.status == "failed" + + +class TestExplosiveDraws: + @pytest.fixture + def report(self, make_var_posterior): + idata = make_var_posterior(explosive_frac=0.15) + return convergence_report(idata, n_lags=1, var_names=["y1", "y2"]) + + def test_exact_stability_statistics(self, report): + assert report.stability.p_explosive == pytest.approx(0.15) + assert report.stability.max_radius == pytest.approx(1.2) + assert report.stability.median() == pytest.approx(0.5) + + def test_explosive_code_raised_to_warning(self, report): + message = next(m for m in report.messages if m.code == "explosive_draws") + assert message.severity == "warning" + assert message.block == "coefficient" + + def test_explosive_draws_never_fail_the_report(self, report): + # Reserved-status decision: `failed` means sampler pathology only. + # Near-unit-root mass is a legitimate posterior statement. + assert report.status == "warnings" + + def test_below_threshold_explosive_mass_is_informational(self, make_var_posterior): + idata = make_var_posterior(explosive_frac=0.02) + report = convergence_report(idata, n_lags=1) + message = next(m for m in report.messages if m.code == "explosive_draws") + assert message.severity == "info" + assert report.status == "passed" + + def test_no_explosive_message_when_all_draws_stable(self, make_var_posterior): + report = convergence_report(make_var_posterior(), n_lags=1) + assert "explosive_draws" not in codes(report) + + +class TestDivergences: + def test_counts_are_exact(self, make_var_posterior): + report = convergence_report(make_var_posterior(divergences=7), n_lags=1) + assert report.divergences == 7 + assert report.n_transitions == 800 + assert report.divergence_rate == pytest.approx(7 / 800) + assert report.sampler_stats_available is True + + def test_low_rate_warns_but_does_not_fail(self, make_var_posterior): + report = convergence_report(make_var_posterior(divergences=3), n_lags=1) + assert report.status == "warnings" + assert next(m for m in report.messages if m.code == "divergences_present").severity == "warning" + + def test_rate_at_one_percent_fails(self, make_var_posterior): + report = convergence_report(make_var_posterior(divergences=8), n_lags=1) + assert report.divergence_rate == pytest.approx(0.01) + assert report.status == "failed" + + def test_nutpie_shaped_stats_match_pymc_shaped(self, make_var_posterior): + nutpie = convergence_report(make_var_posterior(divergences=5, nutpie_shaped=True), n_lags=1) + pymc = convergence_report(make_var_posterior(divergences=5), n_lags=1) + assert nutpie.divergences == pymc.divergences == 5 + assert nutpie.n_transitions == pymc.n_transitions == 800 + + def test_warmup_divergences_are_ignored(self, make_var_posterior): + idata = make_var_posterior(divergences=0, nutpie_shaped=True) + assert "warmup_sample_stats" in idata.groups() + report = convergence_report(idata, n_lags=1) + assert report.divergences == 0 + assert report.status == "passed" + + +class TestMissingSamplerStats: + @pytest.fixture + def report(self, make_var_posterior): + return convergence_report(make_var_posterior(divergences=None), n_lags=1) + + def test_availability_flag_and_none_counts(self, report): + assert report.sampler_stats_available is False + assert report.divergences is None + assert report.n_transitions is None + assert report.divergence_rate is None + + def test_message_is_informational_only(self, report): + message = next(m for m in report.messages if m.code == "sampler_stats_missing") + assert message.severity == "info" + + def test_status_not_degraded_by_missing_stats_alone(self, report): + assert report.status == "passed" + + def test_zero_divergence_message_requires_stats(self, make_var_posterior): + idata = make_var_posterior(bad_coord=(1, 0), divergences=None) + report = convergence_report(idata, n_lags=1) + assert "rhat_without_divergences" not in codes(report) + + +class TestSingleChain: + @pytest.fixture + def report(self, make_var_posterior): + idata = make_var_posterior(n_chains=1, n_draws=800) + return convergence_report(idata, n_lags=1, var_names=["y1", "y2"]) + + def test_rhat_is_none_everywhere(self, report): + assert report.max_rhat is None + assert all(block.max_rhat is None and block.max_rhat_coord is None for block in report.blocks) + + def test_ess_still_computed(self, report): + assert report.min_ess_bulk is not None + assert report.min_ess_bulk > 0 + + def test_single_chain_message_and_capped_status(self, report): + assert "single_chain" in codes(report) + assert report.status != "failed" + + +class TestBlockAssignment: + @pytest.mark.parametrize( + ("name", "expected"), + [ + ("B", "coefficient"), + ("intercept", "intercept"), + ("B_exog", "exog"), + ("sigma_sd", "covariance"), + ("tril_offdiag", "covariance"), + ("L", "covariance"), + ("Sigma", "covariance"), + ("h", "volatility"), + ("R_chol", "volatility"), + ("R_chol_offdiag", "volatility"), + ("v0_h", "volatility"), + ("v1_sigma_eta", "volatility"), + ("v12_phi", "volatility"), + ("structural_shock_matrix", "identification"), + ("P", "identification"), + ("weird_thing", "other"), + ("lambda_", "other"), + ], + ) + def test_variable_maps_to_expected_block(self, name, expected): + posterior = xr.Dataset({name: (("chain", "draw"), np.zeros((2, 3)))}) + assert assign_blocks(posterior) == {expected: [name]} + + def test_unknown_variable_never_raises(self, make_var_posterior): + idata = make_var_posterior(extra_vars={"weird_thing": (3,)}) + report = convergence_report(idata, n_lags=1) + assert blocks_by_name(report)["other"].var_names == ["weird_thing"] + + def test_absent_blocks_are_omitted(self, make_var_posterior): + report = convergence_report(make_var_posterior(), n_lags=1) + assert "exog" not in blocks_by_name(report) + assert "volatility" not in blocks_by_name(report) + + def test_canonical_order_in_to_dataframe(self, make_var_posterior): + idata = make_var_posterior(extra_vars={"weird_thing": (), "B_exog": (1,), "h": (4, 2)}) + frame = convergence_report(idata, n_lags=1).to_dataframe() + assert list(frame.index) == ["coefficient", "intercept", "exog", "covariance", "volatility", "other"] + + def test_variables_sorted_within_a_block(self, make_var_posterior): + idata = make_var_posterior(extra_vars={"Sigma": (2, 2), "tril_offdiag": (1,)}) + assert blocks_by_name(convergence_report(idata, n_lags=1))["covariance"].var_names == [ + "L", + "Sigma", + "sigma_sd", + "tril_offdiag", + ] + + +class TestAdapterHook: + def test_constant_claims_its_own_variables(self): + assert Constant().posterior_var_names() == ("sigma_sd", "tril_offdiag", "L", "Sigma") + + def test_stochastic_volatility_claims_shared_variables(self): + assert StochasticVolatility().posterior_var_names() == ("h", "R_chol", "R_chol_offdiag") + + def test_pandemic_break_claims_its_hyperparameters(self): + assert PandemicBreak(start=3).posterior_var_names() == ("s_march", "s_april", "s_may", "rho") + + def test_claimed_names_route_to_volatility_for_a_time_varying_adapter(self): + posterior = xr.Dataset({ + name: (("chain", "draw"), np.zeros((2, 3))) for name in ("s_march", "s_april", "s_may", "rho", "lambda_") + }) + blocks = assign_blocks(posterior, volatility=PandemicBreak(start=3)) + assert blocks["volatility"] == ["rho", "s_april", "s_march", "s_may"] + assert blocks["other"] == ["lambda_"] + + def test_claimed_names_route_to_covariance_for_a_constant_adapter(self): + posterior = xr.Dataset({"custom_scale": (("chain", "draw"), np.zeros((2, 3)))}) + + class _CustomConstant(Constant): + def posterior_var_names(self) -> tuple[str, ...]: + return ("custom_scale",) + + assert assign_blocks(posterior, volatility=_CustomConstant()) == {"covariance": ["custom_scale"]} + + def test_adapter_without_the_hook_is_fine(self): + posterior = xr.Dataset({"custom_scale": (("chain", "draw"), np.zeros((2, 3)))}) + + class _Bare: + is_time_varying = False + + assert assign_blocks(posterior, volatility=_Bare()) == {"other": ["custom_scale"]} + + +class TestThresholds: + def test_custom_thresholds_flip_status(self, make_var_posterior): + idata = make_var_posterior() + strict = ConvergenceThresholds(rhat_warn=1.0, rhat_fail=1.001) + assert convergence_report(idata, n_lags=1).status == "passed" + assert convergence_report(idata, n_lags=1, thresholds=strict).status == "failed" + + def test_comparison_is_strict_at_the_threshold(self, make_var_posterior): + idata = make_var_posterior() + observed = convergence_report(idata, n_lags=1).max_rhat + assert observed is not None + on_threshold = ConvergenceThresholds(rhat_warn=observed) + just_below = ConvergenceThresholds(rhat_warn=observed - 1e-12) + assert convergence_report(idata, n_lags=1, thresholds=on_threshold).status == "passed" + assert convergence_report(idata, n_lags=1, thresholds=just_below).status == "warnings" + + def test_thresholds_echoed_on_the_report(self, make_var_posterior): + custom = ConvergenceThresholds(ess_warn=10.0, explosive_warn=0.5) + report = convergence_report(make_var_posterior(), n_lags=1, thresholds=custom) + assert report.thresholds == custom + + def test_defaults_match_the_documented_values(self): + thresholds = ConvergenceThresholds() + assert (thresholds.rhat_warn, thresholds.rhat_fail) == (1.01, 1.05) + assert (thresholds.ess_warn, thresholds.ess_fail) == (400.0, 100.0) + assert thresholds.divergence_fail_rate == 0.01 + assert thresholds.explosive_warn == 0.05 + + +class TestRendering: + def test_to_dataframe_columns_and_index(self, make_var_posterior): + frame = convergence_report(make_var_posterior(), n_lags=1).to_dataframe() + assert frame.index.name == "block" + assert list(frame.columns) == [ + "n_variables", + "n_coordinates", + "max_rhat", + "max_rhat_coord", + "min_ess_bulk", + "min_ess_bulk_coord", + "min_ess_tail", + "min_ess_tail_coord", + ] + + def test_summary_mentions_blocks_divergences_and_headlines(self, make_var_posterior): + report = convergence_report(make_var_posterior(divergences=3), n_lags=1, var_names=["y1", "y2"]) + text = report.summary() + assert "coefficient" in text + assert "covariance" in text + assert "divergences: 3" in text + assert "explosive draws" in text + assert "divergences_present" in text + + def test_summary_reports_unavailable_stats(self, make_var_posterior): + text = convergence_report(make_var_posterior(divergences=None), n_lags=1).summary() + assert "divergences: unavailable" in text + + def test_repr_is_one_line(self, make_var_posterior): + text = repr(convergence_report(make_var_posterior(), n_lags=1)) + assert "\n" not in text + assert text.startswith("ConvergenceReport(status='passed'") + + def test_stability_to_dataframe_is_a_single_row(self, make_var_posterior): + frame = convergence_report(make_var_posterior(), n_lags=1).stability.to_dataframe() + assert isinstance(frame, pd.DataFrame) + assert len(frame) == 1 + assert frame.loc["stability", "p_explosive"] == 0.0 + + def test_stability_hdi_brackets_the_median(self, make_var_posterior): + stability = convergence_report(make_var_posterior(), n_lags=1).stability + lower, upper = stability.hdi() + assert lower <= stability.median() <= upper + wide_lower, wide_upper = stability.hdi(prob=0.99) + assert wide_lower <= lower and wide_upper >= upper + + +class TestStabilityArray: + def test_radius_is_read_only(self, make_var_posterior): + radius = convergence_report(make_var_posterior(), n_lags=1).stability.radius + with pytest.raises(ValueError, match="read-only"): + radius[0, 0] = 99.0 + + def test_radius_shape_matches_the_posterior(self, make_var_posterior): + stability = convergence_report(make_var_posterior(), n_lags=1).stability + assert stability.radius.shape == (4, 200) + assert stability.thinned_from is None + + def test_thinning_is_deterministic_and_recorded(self, make_var_posterior): + idata = make_var_posterior(explosive_frac=0.15) + first = convergence_report(idata, n_lags=1, stability_draws=50).stability + second = convergence_report(idata, n_lags=1, stability_draws=50).stability + assert first.thinned_from == 200 + assert first.radius.shape == (4, 50) + np.testing.assert_array_equal(first.radius, second.radius) + assert first.p_explosive == pytest.approx(0.15, abs=0.1) + + def test_thinning_larger_than_the_posterior_is_a_no_op(self, make_var_posterior): + stability = convergence_report(make_var_posterior(), n_lags=1, stability_draws=10_000).stability + assert stability.thinned_from is None + assert stability.radius.shape == (4, 200) + + def test_rejects_non_positive_thinning(self, make_var_posterior): + with pytest.raises(ValueError, match="stability_draws must be positive"): + convergence_report(make_var_posterior(), n_lags=1, stability_draws=0) + + +class TestDeterminism: + def test_two_calls_agree(self, make_var_posterior): + idata = make_var_posterior(bad_coord=(1, 0), divergences=3) + first = convergence_report(idata, n_lags=1, var_names=["y1", "y2"]) + second = convergence_report(idata, n_lags=1, var_names=["y1", "y2"]) + assert first.to_dataframe().equals(second.to_dataframe()) + assert codes(first) == codes(second) + assert first.status == second.status + np.testing.assert_array_equal(first.stability.radius, second.stability.radius) + + +class TestUnsupportedPosteriors: + def test_missing_posterior_group_raises(self): + idata = InferenceData(prior=xr.Dataset({"B": (("chain", "draw"), np.zeros((2, 3)))})) + with pytest.raises(ValueError, match="`posterior` group"): + convergence_report(idata, n_lags=1) + + def test_missing_B_names_fitted_sv(self): + posterior = xr.Dataset({"h": (("chain", "draw", "time"), np.zeros((2, 3, 4)))}) + with pytest.raises(ValueError, match="FittedSV"): + convergence_report(InferenceData(posterior=posterior), n_lags=1) + + +class TestPipelineIntegration: + @pytest.fixture + def fitted(self, make_var_posterior, var_data_2v): + return FittedVAR( + idata=make_var_posterior(), + n_lags=1, + data=var_data_2v, + var_names=["y1", "y2"], + volatility=Constant(), + ) + + def test_fitted_var_method_delegates(self, fitted): + report = fitted.convergence_report() + assert isinstance(report, ConvergenceReport) + assert report.status == "passed" + assert report.blocks[0].max_rhat_coord.startswith("B[y") + + def test_fitted_var_method_forwards_options(self, fitted): + report = fitted.convergence_report( + thresholds=ConvergenceThresholds(rhat_warn=1.0, rhat_fail=1.0), + hdi_prob=0.5, + stability_draws=50, + ) + assert report.status == "failed" + assert report.stability.hdi_prob == 0.5 + assert report.stability.thinned_from == 200 + + def test_identified_var_matches_fitted_var(self, fitted): + identified = fitted.set_identification_strategy(Cholesky(ordering=["y1", "y2"])) + assert identified.convergence_report().to_dataframe().equals(fitted.convergence_report().to_dataframe()) + + def test_legacy_shock_matrix_lands_in_identification_block(self, make_var_posterior, var_data_2v): + fitted = FittedVAR( + idata=make_var_posterior(extra_vars={"structural_shock_matrix": (2, 2)}), + n_lags=1, + data=var_data_2v, + var_names=["y1", "y2"], + volatility=Constant(), + ) + report = fitted.convergence_report() + assert blocks_by_name(report)["identification"].var_names == ["structural_shock_matrix"] + + def test_conjugate_fit_is_supported_with_honest_gaps(self, var_data_2v): + fitted = ConjugateVAR(lags=1, prior=NIWPrior(), draws=1000, tune=20, seed=0).fit(var_data_2v) + report = fitted.convergence_report() + assert report.n_chains == 1 + assert report.max_rhat is None + assert report.blocks[0].max_rhat_coord is None + assert report.sampler_stats_available is False + assert report.divergences is None + assert {"single_chain", "sampler_stats_missing"} <= set(codes(report)) + # Coefficient and Cholesky draws are exact conditional draws, so their + # effective sample size is nominal by construction. + assert report.min_ess_bulk is not None + assert report.min_ess_bulk > 400 + assert report.stability.n_vars == 2 + assert np.isfinite(report.stability.max_radius) + # Nothing here is sampler pathology: no R-hat, no divergences, so the + # single-chain note is the only thing keeping it off "passed". + assert not any(message.severity == "failure" for message in report.messages) + assert report.status == "warnings" + + def test_conjugate_posterior_labels_coordinates_without_coords(self, var_data_2v): + fitted = ConjugateVAR(lags=1, prior=NIWPrior(), draws=1000, tune=20, seed=0).fit(var_data_2v) + coord = fitted.convergence_report().blocks[0].min_ess_bulk_coord + assert coord is not None + assert coord.startswith("B[y") + assert "L1." in coord + + +@pytest.mark.slow +class TestRealNUTSFit: + def test_report_on_a_real_posterior(self, var_data_2v): + from impulso.samplers import NUTSSampler + from impulso.spec import VAR + + fitted = VAR(lags=1).fit( + var_data_2v, + sampler=NUTSSampler(draws=100, tune=100, chains=2, cores=1, random_seed=42, progressbar=False), + ) + report = fitted.convergence_report() + + assert report.status in {"passed", "warnings", "failed"} + assert isinstance(report.divergences, int) + assert report.sampler_stats_available is True + assert {"coefficient", "intercept", "covariance"} <= set(blocks_by_name(report)) + assert np.isfinite(report.stability.p_explosive) + assert isinstance(report.stability, StabilitySummary) + assert report.blocks[0].max_rhat_coord.startswith("B[y") + assert "Convergence report" in report.summary() diff --git a/tests/test_public_api.py b/tests/test_public_api.py index 7a5af57..848b93a 100644 --- a/tests/test_public_api.py +++ b/tests/test_public_api.py @@ -197,3 +197,51 @@ def test_long_run_restriction_in_all(self): import impulso assert "LongRunRestriction" in impulso.__all__ +class TestDiagnosticsPublicAPI: + def test_convergence_report_importable_from_impulso(self): + from impulso import convergence_report + from impulso.diagnostics import convergence_report as direct + + assert convergence_report is direct + + def test_result_types_importable_from_impulso(self): + import impulso + from impulso import diagnostics + + for name in ( + "BlockDiagnostics", + "ConvergenceReport", + "ConvergenceThresholds", + "DiagnosticMessage", + "StabilitySummary", + ): + assert getattr(impulso, name) is getattr(diagnostics, name) + + def test_stability_primitives_importable_from_impulso(self): + from impulso import companion_matrix, spectral_radius + from impulso._stability import companion_matrix as direct_companion + from impulso._stability import spectral_radius as direct_radius + + assert companion_matrix is direct_companion + assert spectral_radius is direct_radius + + def test_diagnostics_names_in_all(self): + import impulso + + for name in ( + "BlockDiagnostics", + "ConvergenceReport", + "ConvergenceThresholds", + "DiagnosticMessage", + "StabilitySummary", + "companion_matrix", + "convergence_report", + "spectral_radius", + ): + assert name in impulso.__all__ + + def test_all_entries_resolve(self): + import impulso + + for name in impulso.__all__: + assert getattr(impulso, name) is not None From 69c51a716aa02e2c7892f29731355ef8c1c96825 Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 04:05:15 +0200 Subject: [PATCH 3/6] docs(diagnostics): reference page, ADR-0008 and domain terms (#142) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add the `diagnostics` reference page to the API toctree, list `companion_matrix` / `spectral_radius` under primitives, and record the two load-bearing judgements in ADR-0008: why metrics are reported per parameter block (a single worst-R-hat over a VAR posterior tells the user nothing about what to change), and why explosive draws warn but never fail (explosiveness is a property of the model, not the sampler — failing there would train users to ignore "failed"). CONTEXT.md gains the four terms, two relationship bullets, and two flagged ambiguities: "companion" is overloaded against ADPRR's calibrated companion `q_cal`, so always write "companion matrix" in full; and `StabilitySummary` here is not the `StabilityResult` planned for the eigenvalue-spectrum work. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Egjd7ToFeb9TQqFnfRQZxV --- CONTEXT.md | 19 +++++++++ ...onvergence-report-blocks-and-thresholds.md | 42 +++++++++++++++++++ docs/reference/diagnostics.md | 25 +++++++++++ docs/reference/index.md | 1 + docs/reference/primitives.md | 8 ++-- docs/references.bib | 15 +++++++ 6 files changed, 107 insertions(+), 3 deletions(-) create mode 100644 docs/adr/0008-convergence-report-blocks-and-thresholds.md create mode 100644 docs/reference/diagnostics.md diff --git a/CONTEXT.md b/CONTEXT.md index f62ffad..5437040 100644 --- a/CONTEXT.md +++ b/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. diff --git a/docs/adr/0008-convergence-report-blocks-and-thresholds.md b/docs/adr/0008-convergence-report-blocks-and-thresholds.md new file mode 100644 index 0000000..3f8d515 --- /dev/null +++ b/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* | — | + +Comparisons are strict, so a metric sitting exactly on a threshold passes. 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. diff --git a/docs/reference/diagnostics.md b/docs/reference/diagnostics.md new file mode 100644 index 0000000..4d5a64b --- /dev/null +++ b/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 +``` diff --git a/docs/reference/index.md b/docs/reference/index.md index 100e069..4a08930 100644 --- a/docs/reference/index.md +++ b/docs/reference/index.md @@ -20,6 +20,7 @@ identification scenario results evidence +diagnostics primitives protocols plotting diff --git a/docs/reference/primitives.md b/docs/reference/primitives.md index 45a90d5..0824adb 100644 --- a/docs/reference/primitives.md +++ b/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 ``` diff --git a/docs/references.bib b/docs/references.bib index 2601e43..358ba3a 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -155,4 +155,19 @@ @article{blanchardQuah1989 volume = {79}, 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}, } From bc37e51d1b186dd4a3e542a4ebfc952dff864209 Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 04:15:07 +0200 Subject: [PATCH 4/6] docs(diagnostics): correct threshold-boundary prose to match the code R-hat/ESS comparisons are strict, but the divergence-rate and explosive- fraction thresholds trigger at the boundary; the docstring and ADR-0008 claimed blanket strictness, contradicting the tested behaviour. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Egjd7ToFeb9TQqFnfRQZxV --- docs/adr/0008-convergence-report-blocks-and-thresholds.md | 2 +- src/impulso/diagnostics.py | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/adr/0008-convergence-report-blocks-and-thresholds.md b/docs/adr/0008-convergence-report-blocks-and-thresholds.md index 3f8d515..5faa75a 100644 --- a/docs/adr/0008-convergence-report-blocks-and-thresholds.md +++ b/docs/adr/0008-convergence-report-blocks-and-thresholds.md @@ -25,7 +25,7 @@ The `identification` block exists but is normally empty: the structural shock ma | Divergence rate | any divergence | 1% | Betancourt (2017) | | Explosive draw fraction | 5% | *never* | — | -Comparisons are strict, so a metric sitting exactly on a threshold passes. 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 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 diff --git a/src/impulso/diagnostics.py b/src/impulso/diagnostics.py index d0a4255..244b03c 100644 --- a/src/impulso/diagnostics.py +++ b/src/impulso/diagnostics.py @@ -143,8 +143,11 @@ def assign_blocks( class ConvergenceThresholds(ImpulsoModel): """Cut-offs separating a passing report from warnings and failures. - Comparisons are strict, so a metric sitting exactly on a threshold - passes: `max_rhat == 1.01` does not warn. + R-hat and ESS comparisons are strict, so a metric sitting exactly on a + threshold passes: `max_rhat == 1.01` does not warn. The rate thresholds + trigger at the boundary: a divergence rate of exactly + `divergence_fail_rate` fails, and an explosive fraction of exactly + `explosive_warn` escalates to a warning. Attributes: rhat_warn: R-hat above this warns. Default 1.01, the rank-normalised From daf27e5575472fa64df64a7dd9fff533b6b97449 Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 23:06:17 +0200 Subject: [PATCH 5/6] style: normalise blank lines at rebase keep-both junctions Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Egjd7ToFeb9TQqFnfRQZxV --- src/impulso/__init__.py | 2 +- tests/conftest.py | 2 ++ tests/test_public_api.py | 2 ++ 3 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/impulso/__init__.py b/src/impulso/__init__.py index 8e4619c..2fd3868 100644 --- a/src/impulso/__init__.py +++ b/src/impulso/__init__.py @@ -95,8 +95,8 @@ "ScenarioResult", "ShockPath", "SignRestriction", - "StationarityTestResult", "StabilitySummary", + "StationarityTestResult", "StochasticVolatility", "StudentT", "VARData", diff --git a/tests/conftest.py b/tests/conftest.py index b78780c..9941d77 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -191,6 +191,8 @@ def permanent_transitory_2v(): "L": L_draws, "idata": az.InferenceData(posterior=posterior), } + + # --------------- Diagnostics posterior factory --------------- # PyMC and nutpie spell almost every sampler statistic differently; only diff --git a/tests/test_public_api.py b/tests/test_public_api.py index 848b93a..4294b50 100644 --- a/tests/test_public_api.py +++ b/tests/test_public_api.py @@ -197,6 +197,8 @@ def test_long_run_restriction_in_all(self): import impulso assert "LongRunRestriction" in impulso.__all__ + + class TestDiagnosticsPublicAPI: def test_convergence_report_importable_from_impulso(self): from impulso import convergence_report From 388ec52188768746c38881ae8b1abdbe23f3dcbe Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 23:25:00 +0200 Subject: [PATCH 6/6] fix(docs): restore closing brace lost resolving the references.bib rebase conflict The rebase conflict resolution merged the new bibliography entry into the blanchardQuah1989 entry, dropping its closing brace. sphinxcontrib-bibtex then failed to parse references.bib, breaking build-docs and docs-linkcheck. --- docs/references.bib | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/references.bib b/docs/references.bib index 358ba3a..fcba56a 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -155,6 +155,8 @@ @article{blanchardQuah1989 volume = {79}, 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},