From f831339c8dfa615262d6a197bc6fd0e61f9ea22e Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 09:25:48 +0200 Subject: [PATCH 1/6] feat(granger): posterior causal-strength query on FittedVAR (#154) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `FittedVAR.granger_causality(cause, effect)` reports the posterior of the Euclidean norm of the tested lag coefficients of `cause` in the `effect` equation, plus the per-lag posteriors behind it. A magnitude, not a test statistic: nothing divides through by the posterior covariance, so a small effect stays distinguishable from an imprecise one. An optional `rope` adds `p_rope = P(||b|| < rope | data)` — practical negligibility at a threshold the analyst names, deliberately with no default. It is not the probability of no causality: `b = 0` has probability zero under continuous coefficient priors, so that quantity needs a spike-and-slab prior Impulso does not fit. The `GrangerCausalityResult` docstring carries the full statement. The engine lives in a new private `_granger.py` (the `conditional_forecast` delegation precedent), and its extraction is pinned by tests against a hand-built posterior with every entry distinct, cross-checked against `lag_matrices` so the lag-major layout cannot drift between the two consumers. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Egjd7ToFeb9TQqFnfRQZxV --- src/impulso/__init__.py | 3 + src/impulso/_granger.py | 183 ++++++++++++++++++++ src/impulso/fitted.py | 74 ++++++++- src/impulso/results.py | 170 +++++++++++++++++++ tests/test_granger.py | 351 +++++++++++++++++++++++++++++++++++++++ tests/test_public_api.py | 18 ++ 6 files changed, 798 insertions(+), 1 deletion(-) create mode 100644 src/impulso/_granger.py create mode 100644 tests/test_granger.py diff --git a/src/impulso/__init__.py b/src/impulso/__init__.py index b1999bb..a0b71f1 100644 --- a/src/impulso/__init__.py +++ b/src/impulso/__init__.py @@ -28,6 +28,7 @@ DynamicMultiplierResult, FEVDResult, ForecastResult, + GrangerCausalityResult, HDIResult, HistoricalDecompositionResult, IntegrationOrderResult, @@ -63,6 +64,7 @@ "FittedVAR", "ForecastResult", "Gaussian", + "GrangerCausalityResult", "HDIResult", "HistoricalDecompositionResult", "IRFResult", @@ -145,6 +147,7 @@ "StudentT": "impulso.observation", "ErrorDistribution": "impulso.protocols", "VolatilityProcess": "impulso.protocols", + "GrangerCausalityResult": "impulso.results", "compute_ma_phi": "impulso._ma", "lag_matrices": "impulso._linalg", } diff --git a/src/impulso/_granger.py b/src/impulso/_granger.py new file mode 100644 index 0000000..fce9320 --- /dev/null +++ b/src/impulso/_granger.py @@ -0,0 +1,183 @@ +"""Granger-causality strength read off a fitted reduced-form posterior. + +`granger_causality(fitted, ...)` backs `FittedVAR.granger_causality`, +reading the tested lag coefficients straight out of an already-fitted +posterior. It works under any volatility process, because the coefficient +matrix `B` is time-invariant under all of them. + +It reports no probability of *no* causality — see the +`GrangerCausalityResult` docstring for why that quantity does not exist +under continuous coefficient priors. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING, Literal + +import numpy as np + +if TYPE_CHECKING: + from impulso.fitted import FittedVAR + from impulso.results import GrangerCausalityResult, IntegrationOrderResult + + +def _coefficient_indices(n_vars: int, cause_index: int, n_lags_tested: int) -> list[int]: + """Columns of `B` holding one cause's coefficients, lag 1 first. + + `VAR.fit` and `ConjugateVAR.fit` both stack the regressors lag-major — + `X = [y_{t-1}, ..., y_{t-p}]`, each block holding all `n` variables — + and the posterior's `coeff` coordinate labels that layout `L1.v1, + L1.v2, ..., L2.v1, ...`. So column `(k - 1) * n + j` of `B` multiplies + `y_{j, t-k}`, and the tested lags of variable `j` are strided `n` apart + from column `j`. + + Args: + n_vars: Number of endogenous variables `n`. + cause_index: Position `j` of the cause in `var_names`. + n_lags_tested: Number of lags to test, counting from lag 1. + + Returns: + Column indices into the trailing axis of `B`, lag 1 first. + """ + return [lag * n_vars + cause_index for lag in range(n_lags_tested)] + + +def _posterior_coefficients(fitted: FittedVAR) -> np.ndarray: + """Read `B` as `(chain, draw, var, coeff)`. + + Hand-built posteriors may order their dimensions arbitrarily; realign + by name when the canonical labels are present, otherwise trust the + positional convention (the same contract as `dynamic_multiplier`). + """ + B_da = fitted.idata.posterior["B"] + if set(B_da.dims) == {"chain", "draw", "var", "coeff"}: + B_da = B_da.transpose("chain", "draw", "var", "coeff") + return np.asarray(B_da.values, dtype=float) + + +def _validate_pair(cause: str, effect: str, var_names: list[str]) -> tuple[int, int]: + """Resolve the cause/effect names to positions, or explain why not. + + Raises: + ValueError: If either name is unknown, or if they are the same + variable. + """ + unknown = [name for name in (cause, effect) if name not in var_names] + if unknown: + raise ValueError(f"unknown variable(s) {unknown}; this model's variables are {var_names}") + if cause == effect: + raise ValueError( + f"cause and effect must be different variables, got {cause!r} for both; " + "Granger causality compares one variable's past against another's own past." + ) + return var_names.index(cause), var_names.index(effect) + + +def _build_result( + fitted: FittedVAR, + cause: str, + effect: str, + *, + test_lags: int | None, + rope: float | None, + standardize: bool, + augmentation_source: Literal["none", "integration_order", "user"] | None = None, + integration_order_result: IntegrationOrderResult | None = None, +) -> GrangerCausalityResult: + """Extract the tested coefficients and package them with their metadata. + + Sole construction site for `GrangerCausalityResult`, so both entry + points agree on what `n_lags_tested`, `augmentation`, and `scale` mean. + + Args: + fitted: The fitted reduced-form posterior to read `B` from. + cause: Variable whose lags are tested. + effect: Variable whose equation they are tested in. + test_lags: Lags to test, counting from lag 1; `None` tests all + fitted lags. + rope: Region of practical equivalence, or `None`. + standardize: Rescale the draws by `sd(cause) / sd(effect)`. + augmentation_source: Provenance of the untested lags. `None` lets + it be inferred: `"none"` when nothing was held back, `"user"` + when the caller shortened `test_lags` by hand. + integration_order_result: Diagnostics to attach, when consulted. + + Returns: + GrangerCausalityResult in the requested reporting units. + + Raises: + ValueError: On unknown or identical variable names, a `test_lags` + outside `[1, n_lags]`, or a non-positive `rope`. + """ + from impulso.results import GrangerCausalityResult + + var_names = list(fitted.var_names) + cause_index, effect_index = _validate_pair(cause, effect, var_names) + + n_lags_fitted = fitted.n_lags + n_lags_tested = n_lags_fitted if test_lags is None else int(test_lags) + if not 1 <= n_lags_tested <= n_lags_fitted: + raise ValueError(f"test_lags must lie in [1, {n_lags_fitted}] (the fitted lag order), got {test_lags}") + if rope is not None and rope <= 0: + raise ValueError( + f"rope must be positive, got {rope}; it is a magnitude in the reporting units of the coefficients." + ) + + B = _posterior_coefficients(fitted) + columns = _coefficient_indices(len(var_names), cause_index, n_lags_tested) + coef_draws = B[..., effect_index, :][..., columns] # (chain, draw, n_lags_tested) + + scale = 1.0 + if standardize: + sd = np.asarray(fitted.data.endog, dtype=float).std(axis=0, ddof=1) + scale = float(sd[cause_index] / sd[effect_index]) + + augmentation = n_lags_fitted - n_lags_tested + if augmentation_source is None: + augmentation_source = "none" if augmentation == 0 else "user" + + return GrangerCausalityResult( + cause=cause, + effect=effect, + n_lags_tested=n_lags_tested, + n_lags_fitted=n_lags_fitted, + augmentation=augmentation, + augmentation_source=augmentation_source, + standardize=standardize, + scale=scale, + rope=rope, + coef_draws=coef_draws * scale, + integration_order_result=integration_order_result, + ) + + +def granger_causality( + fitted: FittedVAR, + cause: str, + effect: str, + *, + rope: float | None = None, + standardize: bool = True, + test_lags: int | None = None, +) -> GrangerCausalityResult: + """Engine behind `FittedVAR.granger_causality`. + + Args: + fitted: Fitted reduced-form posterior. + cause: Variable whose lags are tested. + effect: Variable whose equation they are tested in. + rope: Region of practical equivalence for `p_rope`. + standardize: Report in `sd(effect)` per `sd(cause)` units. + test_lags: Lags to test; `None` tests every fitted lag. + + Returns: + GrangerCausalityResult for the ordered pair. + """ + return _build_result( + fitted, + cause, + effect, + test_lags=test_lags, + rope=rope, + standardize=standardize, + ) diff --git a/src/impulso/fitted.py b/src/impulso/fitted.py index 7a1c79e..f9d56ff 100644 --- a/src/impulso/fitted.py +++ b/src/impulso/fitted.py @@ -16,7 +16,12 @@ if TYPE_CHECKING: from impulso.identified import IdentifiedVAR - from impulso.results import ConditionalForecastResult, DynamicMultiplierResult, ForecastResult + from impulso.results import ( + ConditionalForecastResult, + DynamicMultiplierResult, + ForecastResult, + GrangerCausalityResult, + ) from impulso.scenario import VariablePath @@ -606,6 +611,73 @@ def dynamic_multiplier(self, horizon: int = 20, cumulative: bool = False) -> "Dy cumulative=cumulative, ) + def granger_causality( + self, + cause: str, + effect: str, + *, + rope: float | None = None, + standardize: bool = True, + test_lags: int | None = None, + ) -> "GrangerCausalityResult": + """Posterior strength of one variable's lags in another's equation. + + Reports the posterior of `‖b‖`, the Euclidean norm of the tested + lag coefficients of `cause` in the `effect` equation, together with + the per-lag posteriors behind it. A magnitude, not a test + statistic: nothing here divides by the posterior covariance, so a + small effect and an imprecise one stay distinguishable. + + Supply a `rope` — a region of practical equivalence, in the + reporting units — to also get `p_rope = P(‖b‖ < rope | data)`. + That is a statement about practical negligibility, **not** the + probability of no causality: under continuous coefficient priors + `b = 0` has probability zero regardless of the data. See + `GrangerCausalityResult` for the full statement. + + Granger causality is conditional predictive precedence within this + set of variables, not intervention. An omitted common driver is + enough to manufacture it. + + No identification scheme is involved — this reads the reduced-form + coefficients directly — and it works under any volatility process, + because `B` is time-invariant under all of them. + + Args: + cause: Variable whose lags are tested. + effect: Variable whose equation they are tested in. Must differ + from `cause`. + rope: Region of practical equivalence for `p_rope`. Must be + positive when given; there is deliberately no default. + standardize: If True (default), scale the draws by + `sd(cause) / sd(effect)` so magnitudes read as standard + deviations of the effect per standard deviation of the + cause. The factor is recorded on the result as `scale`. + test_lags: Number of lags to test, counting from lag 1. + Defaults to every fitted lag. Pass `p` after fitting `p + d` + lags to run the Toda-Yamamoto test by hand; the untested + lags are recorded as `augmentation`, never dropped from the + fit. + + Returns: + GrangerCausalityResult for the ordered `cause -> effect` pair. + + Raises: + ValueError: If either name is unknown, if `cause == effect`, if + `test_lags` is outside `[1, n_lags]`, or if `rope` is not + positive. + """ + from impulso._granger import granger_causality + + return granger_causality( + self, + cause, + effect, + rope=rope, + standardize=standardize, + test_lags=test_lags, + ) + def set_identification_strategy(self, scheme: IdentificationScheme) -> "IdentifiedVAR": """Apply a structural identification scheme. diff --git a/src/impulso/results.py b/src/impulso/results.py index 26a880e..bd4778e 100644 --- a/src/impulso/results.py +++ b/src/impulso/results.py @@ -773,6 +773,176 @@ def d_max(self) -> int: return max(self.order.values(), default=0) +class GrangerCausalityResult(ImpulsoBaseModel): + """Posterior Granger-causal strength for one ordered cause-effect pair. + + The headline quantity is the Euclidean norm of the tested lag + coefficients of `cause` in the `effect` equation — `‖b‖ = sqrt(sum_k + b_k^2)`, where `b_k` multiplies `cause_{t-k}` — evaluated draw by draw, + so the result is a posterior for a *magnitude*. That separates "the + effect is small" from "the effect is imprecisely estimated", which a + Wald-style quadratic form (which divides through by the posterior + covariance) deliberately conflates. Per-lag posteriors are reported + alongside the norm, so a single dominant lag stays visible instead of + being buried in it. + + Granger causality is conditional predictive precedence, not + intervention: the statement is that the past of `cause` improves the + prediction of `effect` beyond `effect`'s own past, *within this system + of variables*. Omitted drivers, temporal aggregation, and simultaneous + feedback each break the step from that statement to a mechanism. + + **What `p_rope` is, and what it is not.** `p_rope` is the posterior + probability that the strength norm falls inside the region of practical + equivalence (ROPE) the analyst supplied: `P(‖b‖ < rope | data)`. It is + NOT `P(no causality)`, and it is not a Bayes factor. Under Impulso's + continuous coefficient priors the event `b = 0` has probability zero + both before and after seeing the data, so no dataset can raise it — a + genuine posterior probability of exact non-causality needs a prior that + puts point mass on the null (spike-and-slab / edge inclusion), which + Impulso does not fit. What `p_rope` does say is that the tested + coefficients are jointly *practically* negligible at the magnitude you + declared negligible. Choosing `rope` is the analyst's job and there is + no default, because there is no data-free notion of "small enough"; + that the choice is explicit and recorded is the honesty of the + statement. With `rope=None` the result reports the distribution only + and `p_rope` is `None`. + + **Reporting units.** With `standardize=True` (the default) the draws are + multiplied by `sd(cause) / sd(effect)`, both sample standard deviations + of the estimation data, so a `rope` is read in standard deviations of + the effect per standard deviation of the cause. The factor is recorded + in `scale`. Under lag augmentation the model is fitted in levels, and + the sample standard deviations of integrated series are inflated by + their trends, so standardised magnitudes are most meaningful compared + within one fitted model rather than across models. + + **Toda-Yamamoto metadata.** Under the lag-augmented procedure the model + is fitted with `n_lags_fitted = n_lags_tested + augmentation` lags and + only the first `n_lags_tested` are tested; the augmented lags are never + reported and `n_lags_tested` is never silently changed to match the + fit. `augmentation_source` records where the augmentation came from: + `"none"` when there is none, `"user"` when it was passed explicitly, + and `"integration_order"` when the integration-order diagnostics were + consulted — including the case where they returned `d_max = 0`, so the + record shows that they were consulted. Those diagnostics are attached + as `integration_order_result` whenever they were consulted. + `IntegrationOrderResult.d_max` is a floor rather than a finding + whenever its `inconclusive` list is non-empty (a variable still + integrated at `max_order` is recorded at `max_order`), which is why + `toda_yamamoto` refuses to run in that case rather than + under-augmenting silently. + + Attributes: + cause: Name of the variable whose lags are tested. + effect: Name of the variable whose equation they are tested in. + n_lags_tested: Number of lags of `cause` entering the strength + norm, counting from lag 1. + n_lags_fitted: Lag order of the model that was fitted. + augmentation: `n_lags_fitted - n_lags_tested`, the lags fitted but + deliberately not tested. + augmentation_source: `"none"`, `"user"`, or `"integration_order"`. + standardize: Whether the draws are in standardised units. + scale: The multiplier applied to the raw coefficient draws — the + standardisation factor, or `1.0` when `standardize` is `False`. + rope: The region of practical equivalence, in the reporting units, + or `None` when none was supplied. + coef_draws: Per-lag coefficient draws in the reporting units, shape + `(chains, draws, n_lags_tested)`, lag 1 first. + integration_order_result: The integration-order diagnostics that + fixed the augmentation, when they were consulted. + """ + + cause: str + effect: str + n_lags_tested: int + n_lags_fitted: int + augmentation: int + augmentation_source: Literal["none", "integration_order", "user"] + standardize: bool + scale: float + rope: float | None = Field(default=None, gt=0) + coef_draws: np.ndarray = Field(repr=False) + integration_order_result: IntegrationOrderResult | None = Field(default=None, repr=False) + + @property + def norm_draws(self) -> np.ndarray: + """Per-draw strength norm `‖b‖`, shape `(chains, draws)`.""" + return np.linalg.norm(self.coef_draws, axis=-1) + + @property + def lag_labels(self) -> list[str]: + """Row labels for the tested lags, `["L1", ..., "Lp"]`.""" + return [f"L{lag}" for lag in range(1, self.n_lags_tested + 1)] + + @property + def p_rope(self) -> float | None: + """Posterior probability that `‖b‖` falls inside the ROPE. + + `None` when no `rope` was supplied. Read the class docstring before + reporting it: this is not the probability of no causality. + """ + if self.rope is None: + return None + return float((self.norm_draws < self.rope).mean()) + + def _stacked(self) -> np.ndarray: + """Per-lag draws with the norm appended, shape `(C, D, p + 1)`. + + Stacking lets one `az.hdi` call cover both, and keeps the array + three-dimensional — `az.hdi` reads a bare 2-D array as + `(draw, shape)` rather than `(chain, draw)` and warns about it. + """ + return np.concatenate([self.coef_draws, self.norm_draws[..., np.newaxis]], axis=-1) + + def median(self) -> float: + """Posterior median of the strength norm. + + Returns: + Median of `‖b‖` across all draws, in the reporting units. + """ + return float(np.median(self.norm_draws)) + + def hdi(self, prob: float = 0.89) -> tuple[float, float]: + """Highest-density interval for the strength norm. + + Args: + prob: Probability mass for the interval. Default 0.89. + + Returns: + Tuple of `(lower, upper)` bounds, in the reporting units. + """ + lower, upper = np.asarray(az.hdi(self._stacked(), hdi_prob=prob))[-1] + return float(lower), float(upper) + + def summary(self, prob: float = 0.89) -> pd.DataFrame: + """Per-lag and overall posterior summary. + + Args: + prob: Probability mass for the HDI columns. Default 0.89. + + Returns: + DataFrame indexed by `["L1", ..., "Lp", "norm"]` with columns + `median`, `hdi_lower`, `hdi_upper`. When a `rope` was supplied a + `p_rope` column is added, filled only on the `norm` row — the + ROPE is a statement about the joint magnitude, not about any one + lag. + """ + stacked = self._stacked() + bounds = np.asarray(az.hdi(stacked, hdi_prob=prob)) + frame = pd.DataFrame( + { + "median": np.median(stacked, axis=(0, 1)), + "hdi_lower": bounds[:, 0], + "hdi_upper": bounds[:, 1], + }, + index=pd.Index([*self.lag_labels, "norm"], name="term"), + ) + if self.rope is not None: + frame["p_rope"] = [np.nan] * self.n_lags_tested + [self.p_rope] + return frame + + class VolatilityResult(VARResultBase): """Result from univariate SV fit — posterior of conditional SD. diff --git a/tests/test_granger.py b/tests/test_granger.py new file mode 100644 index 0000000..c6fb069 --- /dev/null +++ b/tests/test_granger.py @@ -0,0 +1,351 @@ +"""Tests for Bayesian Granger causality. + +The load-bearing group is `TestIndexing`. Everything else in the feature — +directionality, calibration — is downstream of +reading the right columns out of the stacked coefficient matrix `B`, and a +lag-major/variable-major slip there is silent: it still returns plausible +numbers, for the wrong pair. Those tests therefore run against a hand-built +posterior whose every entry is distinct, so any slip changes the answer. +""" + +import arviz as az +import numpy as np +import pandas as pd +import pytest +import xarray as xr + +from impulso import VARData +from impulso._granger import _coefficient_indices +from impulso._linalg import lag_matrices +from impulso.conjugate import ConjugateVAR +from impulso.fitted import FittedVAR +from impulso.priors import NIWPrior +from impulso.volatility import Constant + +# --------------- helpers --------------- + + +def _var_data(endog: np.ndarray, names: list[str]) -> VARData: + """Wrap a raw array in VARData with a monthly index.""" + index = pd.date_range("2000-01-31", periods=endog.shape[0], freq="ME") + return VARData(endog=endog, endog_names=names, index=index) + + +def _hand_built_fitted(B_matrix: np.ndarray, *, seed: int = 3, T: int = 60) -> FittedVAR: + """FittedVAR whose posterior is one constant, hand-chosen `B`. + + No MCMC, no estimation: every draw carries exactly `B_matrix`, so the + extracted coefficients can be asserted against literal values. The + endogenous data is random (only its per-column standard deviations + matter, for the standardisation tests). + """ + n_vars = B_matrix.shape[0] + n_lags = B_matrix.shape[1] // n_vars + B = np.broadcast_to(B_matrix, (2, 10, *B_matrix.shape)).copy() + posterior = xr.Dataset({ + "B": xr.DataArray(B, dims=["chain", "draw", "var", "coeff"]), + "intercept": xr.DataArray(np.zeros((2, 10, n_vars)), dims=["chain", "draw", "var"]), + }) + rng = np.random.default_rng(seed) + # Distinct per-column scales so a swapped standardisation ratio shows up. + endog = rng.standard_normal((T, n_vars)) * np.arange(1, n_vars + 1) + data = _var_data(endog, [f"y{i + 1}" for i in range(n_vars)]) + return FittedVAR( + idata=az.InferenceData(posterior=posterior), + n_lags=n_lags, + data=data, + var_names=list(data.endog_names), + volatility=Constant(), + ) + + +# Lag-major layout: columns are L1.y1, L1.y2, L2.y1, L2.y2. Every entry is +# distinct, so picking the wrong row, lag, or variable cannot go unnoticed. +B_2V_2L = np.array([ + [0.11, 0.12, 0.13, 0.14], + [0.21, 0.22, 0.23, 0.24], +]) + +# 3 variables, 2 lags: L1.y1, L1.y2, L1.y3, L2.y1, L2.y2, L2.y3. +B_3V_2L = np.array([ + [0.11, 0.12, 0.13, 0.14, 0.15, 0.16], + [0.21, 0.22, 0.23, 0.24, 0.25, 0.26], + [0.31, 0.32, 0.33, 0.34, 0.35, 0.36], +]) + + +@pytest.fixture +def hand_built(): + """2-variable, 2-lag hand-built posterior.""" + return _hand_built_fitted(B_2V_2L) + + +# --------------- data-generating processes --------------- + + +def _unidirectional_data(seed: int = 0, T: int = 300) -> VARData: + """y2 Granger-causes y1; y1 never feeds back into y2. + + y1_t = 0.4 y1_{t-1} + 0.4 y2_{t-1} + 0.1 e1_t + y2_t = 0.5 y2_{t-1} + 0.1 e2_t + """ + rng = np.random.default_rng(seed) + y = np.zeros((T, 2)) + for t in range(1, T): + y[t, 0] = 0.4 * y[t - 1, 0] + 0.4 * y[t - 1, 1] + 0.1 * rng.standard_normal() + y[t, 1] = 0.5 * y[t - 1, 1] + 0.1 * rng.standard_normal() + return _var_data(y, ["y1", "y2"]) + + +def _null_data(seed: int, T: int = 200) -> VARData: + """Two independent AR(1)s with rho = 0.5 — no causality either way.""" + rng = np.random.default_rng(seed) + y = np.zeros((T, 2)) + for t in range(1, T): + y[t] = 0.5 * y[t - 1] + 0.1 * rng.standard_normal(2) + return _var_data(y, ["y1", "y2"]) + + +# --------------- 1. indexing (load-bearing) --------------- + + +class TestIndexing: + def test_helper_strides_by_n_vars_from_the_cause_column(self): + # Lag-major: lag k's block starts at column k * n_vars. + assert _coefficient_indices(2, 1, 2) == [1, 3] + assert _coefficient_indices(3, 2, 2) == [2, 5] + assert _coefficient_indices(3, 0, 3) == [0, 3, 6] + + def test_helper_matches_the_coeff_coordinate_labels(self): + # The posterior's `coeff` coord is built lag-major in spec.py as + # [f"L{lag}.{name}" for lag in 1..p for name in var_names]; the + # indices must select exactly the cause's labels, in lag order. + names = ["y1", "y2", "y3"] + coeff = [f"L{lag}.{name}" for lag in (1, 2) for name in names] + indices = _coefficient_indices(len(names), names.index("y2"), 2) + assert [coeff[i] for i in indices] == ["L1.y2", "L2.y2"] + + def test_extracts_the_cause_columns_of_the_effect_equation(self, hand_built): + result = hand_built.granger_causality("y2", "y1", standardize=False) + assert result.coef_draws.shape == (2, 10, 2) + # Row y1 (effect), columns L1.y2 and L2.y2. + np.testing.assert_allclose(result.coef_draws, np.broadcast_to([0.12, 0.14], (2, 10, 2))) + + def test_reverse_direction_reads_the_other_equation(self, hand_built): + result = hand_built.granger_causality("y1", "y2", standardize=False) + # Row y2 (effect), columns L1.y1 and L2.y1. + np.testing.assert_allclose(result.coef_draws, np.broadcast_to([0.21, 0.23], (2, 10, 2))) + + def test_test_lags_keeps_the_leading_lags_only(self, hand_built): + result = hand_built.granger_causality("y2", "y1", standardize=False, test_lags=1) + np.testing.assert_allclose(result.coef_draws, np.broadcast_to([0.12], (2, 10, 1))) + assert result.n_lags_tested == 1 + assert result.n_lags_fitted == 2 + assert result.augmentation == 1 + assert result.augmentation_source == "user" + + def test_agrees_with_the_lag_matrices_split(self, hand_built): + # Drift-proofing: `lag_matrices` is the other consumer of the same + # layout, so the two must never disagree about which entry is + # A_k[effect, cause]. + B = hand_built.idata.posterior["B"].values + expected = np.stack([A[..., 0, 1] for A in lag_matrices(B, 2)], axis=-1) + result = hand_built.granger_causality("y2", "y1", standardize=False) + np.testing.assert_allclose(result.coef_draws, expected) + + def test_three_variable_system_picks_the_right_pair(self): + fitted = _hand_built_fitted(B_3V_2L) + result = fitted.granger_causality("y3", "y2", standardize=False) + # Row y2, columns L1.y3 (index 2) and L2.y3 (index 5). + np.testing.assert_allclose(result.coef_draws, np.broadcast_to([0.23, 0.26], (2, 10, 2))) + + expected = np.stack([A[..., 1, 2] for A in lag_matrices(fitted.idata.posterior["B"].values, 2)], axis=-1) + np.testing.assert_allclose(result.coef_draws, expected) + + +# --------------- 2. validation --------------- + + +class TestValidation: + def test_unknown_cause_lists_the_model_variables(self, hand_built): + with pytest.raises(ValueError, match=r"unknown variable\(s\) \['gdp'\]"): + hand_built.granger_causality("gdp", "y1") + + def test_unknown_effect_is_reported_too(self, hand_built): + with pytest.raises(ValueError, match="y1', 'y2'"): + hand_built.granger_causality("y1", "inflation") + + def test_cause_equal_to_effect_is_refused(self, hand_built): + with pytest.raises(ValueError, match="must be different variables"): + hand_built.granger_causality("y1", "y1") + + @pytest.mark.parametrize("test_lags", [0, -1, 3]) + def test_test_lags_outside_the_fitted_order_is_refused(self, hand_built, test_lags): + with pytest.raises(ValueError, match=r"test_lags must lie in \[1, 2\]"): + hand_built.granger_causality("y2", "y1", test_lags=test_lags) + + @pytest.mark.parametrize("rope", [0.0, -0.1]) + def test_non_positive_rope_is_refused(self, hand_built, rope): + with pytest.raises(ValueError, match="rope must be positive"): + hand_built.granger_causality("y2", "y1", rope=rope) + + +# --------------- 3. standardisation --------------- + + +class TestStandardisation: + def test_scale_is_sd_cause_over_sd_effect(self, hand_built): + endog = np.asarray(hand_built.data.endog) + expected = endog[:, 1].std(ddof=1) / endog[:, 0].std(ddof=1) + result = hand_built.granger_causality("y2", "y1") + assert result.scale == pytest.approx(expected) + np.testing.assert_allclose(result.coef_draws, np.broadcast_to([0.12, 0.14], (2, 10, 2)) * expected) + + def test_reverse_direction_inverts_the_ratio(self, hand_built): + forward = hand_built.granger_causality("y2", "y1").scale + reverse = hand_built.granger_causality("y1", "y2").scale + assert forward * reverse == pytest.approx(1.0) + + def test_disabling_standardisation_leaves_the_raw_draws(self, hand_built): + result = hand_built.granger_causality("y2", "y1", standardize=False) + assert result.standardize is False + assert result.scale == 1.0 + np.testing.assert_allclose(result.coef_draws, np.broadcast_to([0.12, 0.14], (2, 10, 2))) + + +# --------------- 4. result surface --------------- + + +class TestResultSurface: + def test_summary_index_is_the_lags_then_the_norm(self, hand_built): + summary = hand_built.granger_causality("y2", "y1", standardize=False).summary() + assert list(summary.index) == ["L1", "L2", "norm"] + assert list(summary.columns) == ["median", "hdi_lower", "hdi_upper"] + assert summary.loc["L1", "median"] == pytest.approx(0.12) + assert summary.loc["norm", "median"] == pytest.approx(np.hypot(0.12, 0.14)) + + def test_hdi_brackets_the_median(self): + fitted = ConjugateVAR(lags=1, prior=NIWPrior(), draws=200, seed=0).fit(_unidirectional_data()) + result = fitted.granger_causality("y2", "y1") + lower, upper = result.hdi() + assert lower <= result.median() <= upper + summary = result.summary() + assert (summary["hdi_lower"] <= summary["median"]).all() + assert (summary["median"] <= summary["hdi_upper"]).all() + + def test_norm_draws_are_the_per_draw_euclidean_norm(self, hand_built): + result = hand_built.granger_causality("y2", "y1", standardize=False) + assert result.norm_draws.shape == (2, 10) + np.testing.assert_allclose(result.norm_draws, np.hypot(0.12, 0.14)) + + def test_p_rope_is_none_without_a_rope(self, hand_built): + result = hand_built.granger_causality("y2", "y1", standardize=False) + assert result.rope is None + assert result.p_rope is None + assert "p_rope" not in result.summary().columns + + def test_p_rope_lands_on_the_norm_row_only(self, hand_built): + result = hand_built.granger_causality("y2", "y1", standardize=False, rope=0.5) + assert result.p_rope == pytest.approx(1.0) # ||b|| = 0.185 < 0.5 in every draw + summary = result.summary() + assert summary.loc["norm", "p_rope"] == pytest.approx(1.0) + assert np.isnan(summary.loc[["L1", "L2"], "p_rope"]).all() + + def test_p_rope_is_a_probability(self): + fitted = ConjugateVAR(lags=1, prior=NIWPrior(), draws=200, seed=0).fit(_unidirectional_data()) + p_rope = fitted.granger_causality("y2", "y1", rope=0.2).p_rope + assert p_rope is not None + assert 0.0 <= p_rope <= 1.0 + + def test_result_is_frozen(self, hand_built): + result = hand_built.granger_causality("y2", "y1") + with pytest.raises(ValueError, match="frozen"): + result.cause = "y1" + + def test_metadata_defaults_to_no_augmentation(self, hand_built): + result = hand_built.granger_causality("y2", "y1") + assert result.n_lags_tested == result.n_lags_fitted == 2 + assert result.augmentation == 0 + assert result.augmentation_source == "none" + assert result.integration_order_result is None + + +# --------------- 5. directionality --------------- + + +class TestDirectionality: + """One decisive unidirectional system, both directions queried. + + `ConjugateVAR` draws in closed form, so 400 draws on 300 observations + is a fraction of a second and needs no MCMC diagnostics. + """ + + @pytest.fixture(scope="class") + def fitted(self): + return ConjugateVAR(lags=1, prior=NIWPrior(), draws=400, seed=0).fit(_unidirectional_data()) + + def test_true_edge_is_large_and_outside_the_rope(self, fitted): + result = fitted.granger_causality("y2", "y1", rope=0.1) + assert result.median() > 0.15 + assert result.p_rope < 0.05 + + def test_absent_edge_is_small_and_inside_the_rope(self, fitted): + result = fitted.granger_causality("y1", "y2", rope=0.1) + assert result.p_rope > 0.5 + + def test_true_edge_dominates_the_absent_one(self, fitted): + assert fitted.granger_causality("y2", "y1").median() > fitted.granger_causality("y1", "y2").median() + + +# --------------- 6. null calibration --------------- + + +def test_null_system_is_mostly_inside_the_rope(): + """Sanity check, not an exact calibration claim. + + Twenty independent null systems, both directions each. Under the + conjugate Minnesota prior the coefficients are shrunk toward a random + walk, so cross-variable draws are pulled toward zero and `p_rope` at a + rope of 0.1 should sit high nearly everywhere; the thresholds are loose + on purpose, because the prior — not a sampling distribution — is what + sets the exact level. + """ + p_ropes = [] + for seed in range(20): + fitted = ConjugateVAR(lags=1, prior=NIWPrior(), draws=200, seed=seed).fit(_null_data(seed)) + for cause, effect in (("y2", "y1"), ("y1", "y2")): + p_ropes.append(fitted.granger_causality(cause, effect, rope=0.1).p_rope) + + assert len(p_ropes) == 40 + assert float(np.median(p_ropes)) > 0.5 + assert min(p_ropes) > 0.1 + + +# --------------- 8. NUTS smoke --------------- + + +@pytest.mark.slow +def test_granger_causality_on_a_nuts_fit(var_data_2v): + """The reduced-form read must work on a real PyMC posterior too.""" + from impulso import VAR + from impulso.samplers import NUTSSampler + + fitted = VAR(lags=1).fit(var_data_2v, NUTSSampler(draws=100, tune=100, chains=2, cores=1)) + + # Pin the coefficient layout against the real posterior's labels. + coeff = [str(c) for c in fitted.idata.posterior["B"].coords["coeff"].values] + assert coeff == ["L1.y1", "L1.y2"] + + for cause, effect in (("y2", "y1"), ("y1", "y2")): + result = fitted.granger_causality(cause, effect, rope=0.1) + summary = result.summary()[["median", "hdi_lower", "hdi_upper"]] + assert np.isfinite(summary.to_numpy(dtype=float)).all() + lower, upper = result.hdi() + assert lower <= result.median() <= upper + p_rope = result.p_rope + assert p_rope is not None + assert 0.0 <= p_rope <= 1.0 + + index = list(var_data_2v.endog_names).index(cause) + expected = fitted.idata.posterior["B"].values[..., list(var_data_2v.endog_names).index(effect), index] + np.testing.assert_allclose(result.coef_draws[..., 0] / result.scale, expected) diff --git a/tests/test_public_api.py b/tests/test_public_api.py index e348cc1..34e12fe 100644 --- a/tests/test_public_api.py +++ b/tests/test_public_api.py @@ -121,6 +121,24 @@ def test_enable_runtime_checks_drives_pipeline(self): assert "VIOLATION_CAUGHT" in result.stdout, result.stdout +class TestGrangerPublicAPI: + def test_granger_causality_result_importable(self): + from impulso import GrangerCausalityResult + from impulso.results import GrangerCausalityResult as direct + + assert GrangerCausalityResult is direct + + def test_granger_names_in_all(self): + import impulso + + assert "GrangerCausalityResult" in impulso.__all__ + + def test_granger_causality_is_a_fitted_var_method(self): + from impulso.fitted import FittedVAR + + assert callable(FittedVAR.granger_causality) + + class TestVolatilityPublicAPI: def test_constant_importable_from_impulso(self): from impulso import Constant From 1ea5cf64a96e6d5187f89869c57068cf2ebba702 Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 09:26:17 +0200 Subject: [PATCH 2/6] feat(granger): toda_yamamoto() lag-augmented mode (#154) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `toda_yamamoto(data, cause, effect)` runs the Toda-Yamamoto (1995) procedure for possibly-integrated systems: fit the VAR in levels with `p + d` lags, test only the first `p`. The augmented lags are never tested and the reported test lag order is never silently changed to match the fit — the result carries `n_lags_tested` and `n_lags_fitted` separately, with `augmentation` and `augmentation_source` recording where the extra lags came from. `d` comes from `integration_order` unless the caller pins it. Honouring the consumer contract frozen by #140/#197: when the diagnostics leave anything in `inconclusive`, `d_max` is a floor rather than a finding, so this refuses to run — naming the variables, pointing at `.summary()` and at the `d=` override — rather than under-augmenting silently. An explicit `d` skips the diagnostics entirely, so the route also works without statsmodels installed. The fit uses the closed-form conjugate estimator, since augmentation inflates the lag order; exogenous regressors it cannot consume are refused with a message naming the manual `VAR(lags=p + d)` route. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Egjd7ToFeb9TQqFnfRQZxV --- src/impulso/__init__.py | 3 + src/impulso/_granger.py | 209 +++++++++++++++++++++++++++++++++++++-- tests/test_granger.py | 171 +++++++++++++++++++++++++++++++- tests/test_public_api.py | 7 ++ 4 files changed, 379 insertions(+), 11 deletions(-) diff --git a/src/impulso/__init__.py b/src/impulso/__init__.py index a0b71f1..8be0253 100644 --- a/src/impulso/__init__.py +++ b/src/impulso/__init__.py @@ -10,6 +10,7 @@ if TYPE_CHECKING: from types import ModuleType + from impulso._granger import toda_yamamoto from impulso._linalg import lag_matrices from impulso._ma import compute_ma_phi from impulso.conjugate import ConjugateVAR @@ -101,6 +102,7 @@ "kpss_test", "lag_matrices", "select_lag_order", + "toda_yamamoto", ] @@ -148,6 +150,7 @@ "ErrorDistribution": "impulso.protocols", "VolatilityProcess": "impulso.protocols", "GrangerCausalityResult": "impulso.results", + "toda_yamamoto": "impulso._granger", "compute_ma_phi": "impulso._ma", "lag_matrices": "impulso._linalg", } diff --git a/src/impulso/_granger.py b/src/impulso/_granger.py index fce9320..2e4dbb9 100644 --- a/src/impulso/_granger.py +++ b/src/impulso/_granger.py @@ -1,11 +1,19 @@ -"""Granger-causality strength read off a fitted reduced-form posterior. - -`granger_causality(fitted, ...)` backs `FittedVAR.granger_causality`, -reading the tested lag coefficients straight out of an already-fitted -posterior. It works under any volatility process, because the coefficient -matrix `B` is time-invariant under all of them. - -It reports no probability of *no* causality — see the +"""Granger-causality strength and the Toda-Yamamoto lag-augmented mode. + +Two entry points share one result builder, so the metadata a +`GrangerCausalityResult` carries cannot drift between them: + +* `granger_causality(fitted, ...)` — backs `FittedVAR.granger_causality`, + reading the tested lag coefficients straight out of an already-fitted + posterior. Works under any volatility process, because the coefficient + matrix `B` is time-invariant under all of them. +* `toda_yamamoto(data, ...)` — the lag-augmented procedure of Toda and + Yamamoto (1995) for possibly-integrated systems: fit `p + d` lags, test + only the first `p`. It resolves `d` from the integration-order + diagnostics unless the caller pins it, and fits with the closed-form + conjugate estimator so the extra lags cost seconds rather than minutes. + +Neither reports a probability of *no* causality — see the `GrangerCausalityResult` docstring for why that quantity does not exist under continuous coefficient priors. """ @@ -16,10 +24,15 @@ import numpy as np +from impulso.data import VARData + if TYPE_CHECKING: from impulso.fitted import FittedVAR + from impulso.priors import NIWPrior from impulso.results import GrangerCausalityResult, IntegrationOrderResult +_CRITERIA = ("aic", "bic", "hq") + def _coefficient_indices(n_vars: int, cause_index: int, n_lags_tested: int) -> list[int]: """Columns of `B` holding one cause's coefficients, lag 1 first. @@ -181,3 +194,183 @@ def granger_causality( rope=rope, standardize=standardize, ) + + +def _resolve_lag_order(data: VARData, lags: int | str, max_lags: int) -> int: + """Resolve `lags` to a positive integer `p`, selecting if asked.""" + if isinstance(lags, str): + if lags not in _CRITERIA: + raise ValueError(f"lags must be an int or one of {_CRITERIA}, got {lags!r}") + from impulso._lag_selection import select_lag_order + + return int(getattr(select_lag_order(data, max_lags=max_lags), lags)) + p = int(lags) + if p < 1: + raise ValueError(f"lags must be >= 1, got {p}") + return p + + +def _resolve_augmentation( + data: VARData, + d: int | None, + integration_order_result: IntegrationOrderResult | None, + *, + max_order: int, + alpha: float, + regression: Literal["c", "ct"], +) -> tuple[int, Literal["integration_order", "user"], IntegrationOrderResult | None]: + """Fix the augmentation `d`, either from the caller or from diagnostics. + + An explicit `d` skips the diagnostics entirely — deliberately, so the + procedure runs without `statsmodels` installed and so a decision the + analyst has already made is not silently re-litigated. + + Raises: + ValueError: If `d` is negative, or if the diagnostics left any + variable in `inconclusive`, where `d_max` is a floor rather + than a finding and would under-augment the test. + """ + if d is not None: + if d < 0: + raise ValueError(f"d must be non-negative, got {d}") + return int(d), "user", None + + consulted = integration_order_result + if consulted is None: + from impulso._stationarity import integration_order + + consulted = integration_order(data, max_order=max_order, alpha=alpha, regression=regression) + if consulted.inconclusive: + raise ValueError( + f"the integration order of {consulted.inconclusive} is unsettled: each of these is either " + f"still non-stationary at max_order={consulted.max_order} (so its recorded order is a floor, " + "not a finding) or had ADF and KPSS disagree where the search stopped. d_max would then " + "under-augment, and under-augmented Toda-Yamamoto inference is invalid. Inspect the full " + "table with integration_order(...).summary(), then pass the augmentation explicitly as " + "d= once you have decided." + ) + return consulted.d_max, "integration_order", consulted + + +def toda_yamamoto( + data: VARData, + cause: str, + effect: str, + *, + lags: int | Literal["aic", "bic", "hq"] = "aic", + max_lags: int = 12, + d: int | None = None, + integration_order_result: IntegrationOrderResult | None = None, + max_order: int = 2, + alpha: float = 0.05, + regression: Literal["c", "ct"] = "c", + rope: float | None = None, + standardize: bool = True, + prior: NIWPrior | None = None, + draws: int = 1000, + seed: int | None = None, +) -> GrangerCausalityResult: + """Granger causality with Toda-Yamamoto lag augmentation. + + Toda and Yamamoto (1995) make Granger-causality inference valid without + first deciding the integration and cointegration structure: fit the VAR + in levels with `p + d` lags, where `p` is the lag order you would have + chosen and `d` the highest integration order in the system, then test + only the first `p` lags. The extra `d` lags are never tested — they + exist to restore the standard asymptotics — and this function never + silently changes the reported test lag order to match the fitted one: + the result carries `n_lags_tested` and `n_lags_fitted` separately. + + `d` comes from `integration_order` unless it is passed explicitly. When + the diagnostics leave any variable in `inconclusive`, `d_max` is a + floor rather than a finding, so this function refuses to run rather + than under-augment; read the full table and pass `d=` yourself. + + The fit uses the closed-form conjugate estimator (`ConjugateVAR` with + an `NIWPrior`) because augmentation inflates the lag order and the + conjugate path draws in closed form. For the NUTS estimator, a + stochastic-volatility process, or exogenous regressors, run the + procedure by hand — it is three calls: + + ```python + d = integration_order(data).d_max + fitted = VAR(lags=p + d).fit(data) + fitted.granger_causality(cause, effect, test_lags=p) + ``` + + Args: + data: Endogenous data, in levels. Exogenous regressors are not + supported here (the conjugate estimator does not consume them). + cause: Variable whose lags are tested. + effect: Variable whose equation they are tested in. + lags: Test lag order `p`, or an information criterion to select it + with (`"aic"`, `"bic"`, `"hq"`). + max_lags: Upper bound when `lags` is a criterion. + d: Augmentation to use. Passing it skips the diagnostics entirely, + and records `augmentation_source="user"`. + integration_order_result: Diagnostics to reuse instead of running + `integration_order` again. Ignored when `d` is given. + max_order: `max_order` for `integration_order`, when it is run. + alpha: Significance level for `integration_order`, when it is run. + regression: Deterministic terms for `integration_order`'s level + test, when it is run. + rope: Region of practical equivalence for `p_rope`. + standardize: Report in `sd(effect)` per `sd(cause)` units. Note + that the standard deviations of integrated series carry their + trends, so standardised magnitudes compare best within one fit. + prior: Conjugate prior for the fit. Defaults to `NIWPrior()`. + draws: Posterior draws to retain. + seed: Seed for the conjugate sampler. + + Returns: + GrangerCausalityResult with `n_lags_tested = p`, + `augmentation = d`, and the consulted diagnostics attached when + they were run. + + Raises: + ValueError: If `data` carries exogenous regressors, if the names + are unknown or identical, if `lags` or `d` is invalid, or if + the integration-order diagnostics are inconclusive. + """ + from impulso.conjugate import ConjugateVAR + from impulso.priors import NIWPrior + + if data.exog is not None: + raise ValueError( + "toda_yamamoto fits with the conjugate estimator, which estimates endogenous dynamics " + f"only, and this VARData carries exogenous regressors {list(data.exog_names or [])}. Run " + "the procedure manually instead: d = integration_order(data).d_max, then " + "fitted = VAR(lags=p + d).fit(data), then " + "fitted.granger_causality(cause, effect, test_lags=p)." + ) + # Validate the pair before any fitting or diagnostics, so a typo costs + # nothing. + _validate_pair(cause, effect, list(data.endog_names)) + + p = _resolve_lag_order(data, lags, max_lags) + augmentation, source, consulted = _resolve_augmentation( + data, + d, + integration_order_result, + max_order=max_order, + alpha=alpha, + regression=regression, + ) + + fitted = ConjugateVAR( + lags=p + augmentation, + prior=prior if prior is not None else NIWPrior(), + draws=draws, + seed=seed, + ).fit(data) + + return _build_result( + fitted, + cause, + effect, + test_lags=p, + rope=rope, + standardize=standardize, + augmentation_source=source, + integration_order_result=consulted, + ) diff --git a/tests/test_granger.py b/tests/test_granger.py index c6fb069..a56685e 100644 --- a/tests/test_granger.py +++ b/tests/test_granger.py @@ -1,7 +1,7 @@ -"""Tests for Bayesian Granger causality. +"""Tests for Bayesian Granger causality and the Toda-Yamamoto mode. The load-bearing group is `TestIndexing`. Everything else in the feature — -directionality, calibration — is downstream of +directionality, calibration, the augmentation contract — is downstream of reading the right columns out of the stacked coefficient matrix `B`, and a lag-major/variable-major slip there is silent: it still returns plausible numbers, for the wrong pair. Those tests therefore run against a hand-built @@ -14,12 +14,13 @@ import pytest import xarray as xr -from impulso import VARData +from impulso import GrangerCausalityResult, VARData, toda_yamamoto from impulso._granger import _coefficient_indices from impulso._linalg import lag_matrices from impulso.conjugate import ConjugateVAR from impulso.fitted import FittedVAR from impulso.priors import NIWPrior +from impulso.results import IntegrationOrderResult from impulso.volatility import Constant # --------------- helpers --------------- @@ -106,6 +107,21 @@ def _null_data(seed: int, T: int = 200) -> VARData: return _var_data(y, ["y1", "y2"]) +def _i1_unidirectional(seed: int = 7, T: int = 400) -> VARData: + """Both series I(1); x's increments drive y's, never the reverse. + + x is a driftless random walk; y accumulates 0.5 times x's increment + plus its own small noise. Seed chosen so that `integration_order` + settles both series at d = 1 with nothing left inconclusive. + """ + rng = np.random.default_rng(seed) + x = np.cumsum(rng.standard_normal(T)) + y = np.zeros(T) + for t in range(2, T): + y[t] = y[t - 1] + 0.5 * (x[t - 1] - x[t - 2]) + 0.05 * rng.standard_normal() + return _var_data(np.column_stack([x, y]), ["x", "y"]) + + # --------------- 1. indexing (load-bearing) --------------- @@ -321,6 +337,155 @@ def test_null_system_is_mostly_inside_the_rope(): assert min(p_ropes) > 0.1 +# --------------- 7. Toda-Yamamoto contract --------------- + + +def _integration_order_result(order: dict[str, int], inconclusive: list[str]) -> IntegrationOrderResult: + """Hand-built diagnostics, so the refusal contract needs no statsmodels.""" + return IntegrationOrderResult( + order=order, + alpha=0.05, + max_order=2, + regression="c", + inconclusive=inconclusive, + table=pd.DataFrame( + {"joint_status": ["inconclusive"] * len(order)}, + index=pd.MultiIndex.from_tuples([(name, 0) for name in order], names=["variable", "d"]), + ), + ) + + +class TestTodaYamamoto: + def test_inconclusive_diagnostics_are_refused_by_name(self): + diagnostics = _integration_order_result({"y1": 1, "y2": 2}, ["y2"]) + with pytest.raises(ValueError, match="y2") as excinfo: + toda_yamamoto( + _null_data(0), + "y2", + "y1", + lags=1, + integration_order_result=diagnostics, + ) + message = str(excinfo.value) + assert "under-augment" in message + assert "summary()" in message + assert "d=" in message + + def test_explicit_d_splits_tested_from_fitted_lags(self): + result = toda_yamamoto(_null_data(0), "y2", "y1", lags=1, d=1) + assert result.n_lags_tested == 1 + assert result.n_lags_fitted == 2 + assert result.augmentation == 1 + assert result.augmentation_source == "user" + assert result.integration_order_result is None + # The augmented lag is fitted but never reported. + assert list(result.summary().index) == ["L1", "norm"] + assert result.coef_draws.shape[-1] == 1 + + def test_injected_clean_diagnostics_are_consulted_and_attached(self): + diagnostics = _integration_order_result({"y1": 1, "y2": 1}, []) + result = toda_yamamoto( + _null_data(0), + "y2", + "y1", + lags=1, + integration_order_result=diagnostics, + ) + assert result.augmentation == 1 + assert result.augmentation_source == "integration_order" + assert result.integration_order_result is diagnostics + + def test_d_max_of_zero_still_records_that_diagnostics_ran(self): + diagnostics = _integration_order_result({"y1": 0, "y2": 0}, []) + result = toda_yamamoto( + _null_data(0), + "y2", + "y1", + lags=2, + integration_order_result=diagnostics, + ) + assert result.augmentation == 0 + assert result.augmentation_source == "integration_order" + assert result.n_lags_fitted == result.n_lags_tested == 2 + + def test_explicit_d_needs_no_statsmodels(self, monkeypatch): + """The `d=` route must not even try to import the optional extra.""" + import impulso._stationarity as stationarity + + def _absent(module, *, extra): + raise ImportError(f"{module} is not installed") + + monkeypatch.setattr(stationarity, "require", _absent) + result = toda_yamamoto(_null_data(0), "y2", "y1", lags=1, d=1) + assert result.augmentation_source == "user" + assert np.isfinite(result.median()) + + def test_exogenous_data_points_at_the_manual_route(self): + data = _null_data(0) + with_exog = VARData( + endog=np.asarray(data.endog), + endog_names=list(data.endog_names), + exog=np.ones((len(data.index), 1)), + exog_names=["trend"], + index=data.index, + ) + with pytest.raises(ValueError, match="granger_causality") as excinfo: + toda_yamamoto(with_exog, "y2", "y1", lags=1, d=0) + assert "VAR(lags=p + d)" in str(excinfo.value) + + def test_unknown_variable_is_refused_before_fitting(self): + with pytest.raises(ValueError, match=r"unknown variable\(s\)"): + toda_yamamoto(_null_data(0), "gdp", "y1", lags=1, d=0) + + @pytest.mark.parametrize(("lags", "match"), [("nope", "lags must be an int"), (0, "lags must be >= 1")]) + def test_invalid_lag_specification_is_refused(self, lags, match): + with pytest.raises(ValueError, match=match): + toda_yamamoto(_null_data(0), "y2", "y1", lags=lags, d=0) + + def test_negative_d_is_refused(self): + with pytest.raises(ValueError, match="d must be non-negative"): + toda_yamamoto(_null_data(0), "y2", "y1", lags=1, d=-1) + + def test_criterion_string_selects_the_test_lag_order(self): + result = toda_yamamoto(_null_data(0), "y2", "y1", lags="bic", max_lags=4, d=1) + assert result.n_lags_tested >= 1 + assert result.n_lags_fitted == result.n_lags_tested + 1 + + def test_returns_a_granger_causality_result(self): + assert isinstance(toda_yamamoto(_null_data(0), "y2", "y1", lags=1, d=0), GrangerCausalityResult) + + +class TestTodaYamamotoEndToEnd: + """Full path on an I(1) system, diagnostics included.""" + + def test_augmented_fit_recovers_the_direction(self): + pytest.importorskip("statsmodels") + data = _i1_unidirectional() + + forward = toda_yamamoto(data, "x", "y", lags=2, rope=0.1, draws=400, seed=0) + assert forward.augmentation_source == "integration_order" + assert forward.augmentation >= 1 + assert forward.n_lags_fitted == forward.n_lags_tested + forward.augmentation + assert forward.integration_order_result is not None + assert forward.integration_order_result.d_max == forward.augmentation + + # Reuse the diagnostics rather than re-running them for the reverse + # direction; the injected path is asserted separately above. + reverse = toda_yamamoto( + data, + "y", + "x", + lags=2, + rope=0.1, + draws=400, + seed=0, + integration_order_result=forward.integration_order_result, + ) + assert forward.p_rope is not None + assert reverse.p_rope is not None + assert forward.p_rope < reverse.p_rope + + # --------------- 8. NUTS smoke --------------- diff --git a/tests/test_public_api.py b/tests/test_public_api.py index 34e12fe..bc79dfd 100644 --- a/tests/test_public_api.py +++ b/tests/test_public_api.py @@ -122,6 +122,12 @@ def test_enable_runtime_checks_drives_pipeline(self): class TestGrangerPublicAPI: + def test_toda_yamamoto_importable(self): + from impulso import toda_yamamoto + from impulso._granger import toda_yamamoto as direct + + assert toda_yamamoto is direct + def test_granger_causality_result_importable(self): from impulso import GrangerCausalityResult from impulso.results import GrangerCausalityResult as direct @@ -131,6 +137,7 @@ def test_granger_causality_result_importable(self): def test_granger_names_in_all(self): import impulso + assert "toda_yamamoto" in impulso.__all__ assert "GrangerCausalityResult" in impulso.__all__ def test_granger_causality_is_a_fitted_var_method(self): From c26a04001a3f34c3a1baf0ab9d23664e66de0618 Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 09:32:03 +0200 Subject: [PATCH 3/6] docs(granger): how-to, reference, CONTEXT terms, and ADR-0010 (#154) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New how-to covering the fitted-model query, how to read `summary()`, why there is no probability of no causality, the Toda-Yamamoto happy path plus its refusal and the `d=` override, the manual NUTS route, and a worked carbon-dioxide/temperature example with an explicit statement of what it does and does not license (predictive precedence not intervention; omitted forcings; bidirectional physical coupling; annual aggregation). Cross-links with the climate-pitfalls page both ways. New reference page for `toda_yamamoto`; `GrangerCausalityResult` added to the results page. CONTEXT gains three terms — Granger causality, ROPE, Toda-Yamamoto augmentation — plus the two relationships that place the query on `FittedVAR` and wire the augmentation to the integration-order contract. ADR-0010 records the decision: the norm of the tested coefficients as the headline with an analyst-supplied ROPE, against the rejected alternatives (a Wald quadratic headline, a fixed default epsilon, spike-and-slab, and Savage-Dickey). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Egjd7ToFeb9TQqFnfRQZxV --- CONTEXT.md | 14 ++ ...anger-strength-is-a-magnitude-statement.md | 19 ++ docs/how-to/climate-pitfalls.md | 5 +- docs/how-to/granger-causality.md | 238 ++++++++++++++++++ docs/how-to/index.md | 1 + docs/reference/causality.md | 31 +++ docs/reference/index.md | 1 + docs/reference/results.md | 1 + src/impulso/_granger.py | 4 +- 9 files changed, 311 insertions(+), 3 deletions(-) create mode 100644 docs/adr/0010-granger-strength-is-a-magnitude-statement.md create mode 100644 docs/how-to/granger-causality.md create mode 100644 docs/reference/causality.md diff --git a/CONTEXT.md b/CONTEXT.md index 30e35d7..f5291de 100644 --- a/CONTEXT.md +++ b/CONTEXT.md @@ -132,6 +132,18 @@ _Avoid_: "order of differencing" for `d_max` — `d_max` is the maximum across t The number of independent long-run relationships among integrated series, from the Johansen procedure (`johansen_test`). Both sequential tests are reported — `rank_trace` and `rank_max_eigen` — and `rank` is `rank_trace` by documented convention. Decisions rest on **critical values, not p-values** (MacKinnon-Haug-Michelis 1996 tables, as vendored by statsmodels), which is why `alpha` is restricted to 0.10 / 0.05 / 0.01. Rank ≥ 1 means differencing every series discards the long-run relationship. A vector error-correction model (VECM) is **out of scope**; the recommended response is a VAR in levels (the Sims–Stock–Watson stance; the Minnesota prior already shrinks toward random walks). _Avoid_: "number of cointegrating vectors" in API surface (fine in prose); "cointegration test" without saying which statistic, since trace and max-eigen can disagree. +**Granger causality**: +Conditional predictive precedence: the past of one variable improves the prediction of another beyond that other's own past, *within the fitted system of variables*. Reported by `FittedVAR.granger_causality(cause, effect)` as the posterior of the strength norm `‖b‖` over the tested lags of the cause in the effect's equation, with the per-lag posteriors alongside. Ordered and directional — the two orderings are separate queries with unrelated answers. Reduced-form: no identification scheme is involved, and `B` is time-invariant under every volatility process, so it needs no `at`. +_Avoid_: "X causes Y" for a Granger result, and "causal effect" — the finding is about information sets, not interventions; an omitted common driver manufactures it. Reserve "effect" for identified structural objects (IRFs, counterfactuals). + +**ROPE (region of practical equivalence)**: +The magnitude below which the analyst declares a relationship practically negligible, supplied as `rope=` and echoed on the result. `p_rope = P(‖b‖ < rope | data)` is the only probability statement the Granger surface makes. There is deliberately no default: the threshold is the analyst's judgement, and it travels with the number that came from it. Without a `rope` the result reports the distribution and `p_rope` is `None`. +_Avoid_: "probability of no causality" / "probability the coefficient is zero" for `p_rope` — under continuous coefficient priors `P(b = 0) = 0` before and after the data, and an edge-inclusion probability needs a spike-and-slab prior Impulso does not fit (see ADR-0010). Also avoid bare "threshold" — the ROPE is on the magnitude, not on a p-value. + +**Toda-Yamamoto augmentation**: +Fitting a VAR in levels with `p + d` lags and testing only the first `p`, which restores standard Granger inference on possibly-integrated series without differencing (Toda & Yamamoto 1995). `n_lags_tested` and `n_lags_fitted` are separate fields on `GrangerCausalityResult` — the test lag order is never silently changed to match the fit, and the augmented lags never appear in `summary()`. `toda_yamamoto` consumes the frozen `IntegrationOrderResult` contract (`order` / `d_max` / `inconclusive`): it **refuses to run** when `inconclusive` is non-empty, because `d_max` is then a floor and an under-augmented test is invalid rather than imprecise. `augmentation_source` records the provenance — `"user"` for an explicit `d=` (which skips the diagnostics entirely), `"integration_order"` when they were consulted, including when they returned `d_max = 0`. +_Avoid_: "extra lags" without saying they are untested; "corrected for non-stationarity" — nothing is corrected, the asymptotics are restored. + ## Relationships - A **VAR** carries one **prior**, one **volatility process**, and one **observation error distribution**. @@ -148,6 +160,8 @@ _Avoid_: "number of cointegrating vectors" in API surface (fine in prose); "coin - A **VAR** is estimated by NUTS; a **ConjugateVAR** is estimated analytically with a Metropolis step on hyperparameters. Both produce a **FittedVAR**. - 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 **FittedVAR** answers **Granger causality** queries on its own — the coefficients are reduced-form and time-invariant, so neither an identification scheme nor an `at` is involved. The query never refits: `test_lags` selects which of the already-fitted lags are tested. +- **Toda-Yamamoto augmentation** consumes an **integration order**: `toda_yamamoto` reads `d_max` (after checking `inconclusive`), fits `p + d` lags with a **ConjugateVAR**, and produces the same `GrangerCausalityResult` the `FittedVAR` query does — which is why the manual route (`VAR(lags=p + d).fit(...)` then `granger_causality(..., test_lags=p)`) is equivalent and is the documented escape hatch for exogenous regressors, NUTS, and stochastic volatility. - 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. ## Example dialogue diff --git a/docs/adr/0010-granger-strength-is-a-magnitude-statement.md b/docs/adr/0010-granger-strength-is-a-magnitude-statement.md new file mode 100644 index 0000000..67b8209 --- /dev/null +++ b/docs/adr/0010-granger-strength-is-a-magnitude-statement.md @@ -0,0 +1,19 @@ +# Granger-causal strength is a magnitude statement, not a test + +Impulso reports Granger causality as the posterior of `‖b‖`, the Euclidean norm over the tested lags of the cause's coefficients in the effect's equation, evaluated draw by draw — median, highest-density interval (HDI), and the per-lag posteriors alongside it. The only probability statement attached is `p_rope = P(‖b‖ < rope | data)`, against a region of practical equivalence (ROPE) the analyst supplies (Kruschke's formulation). With no `rope` the result reports the distribution and nothing else. There is no default ROPE, and no probability of *no* causality anywhere in the API. + +## Considered options + +- **A Wald-style quadratic form as the headline** — `b'V⁻¹b` over the posterior, the Bayesian echo of the classical test — rejected: dividing through by the posterior covariance conflates "the effect is small" with "the effect is precisely estimated", which is exactly the distinction a posterior is for. A large quadratic can mean a large effect or a tight posterior around a small one, and the reader cannot tell which. It remains a reasonable *supplementary* statistic and may arrive later; it is not the headline. +- **A fixed default epsilon for the ROPE** (0.05, say, in standardised units) — rejected: any default is arbitrary, and a default is precisely what makes a threshold invisible. That the analyst had to name a magnitude, and that the magnitude travels with the number in `GrangerCausalityResult.rope`, is the honesty of the statement. +- **Spike-and-slab / edge-inclusion priors, to report `P(b = 0 | data)` properly** — rejected for this slice, not on principle. It is the only construction that makes the quantity users actually ask for well defined, but it is a different prior and a different estimator, not a post-processing step. Reporting an edge-inclusion probability without fitting one would be a fabrication, and the issue explicitly rules it out. +- **Savage-Dickey density ratio at `b = 0`** — rejected: it gives a Bayes factor for the point null under a continuous prior, but its value depends on the prior density at zero, which under the Minnesota prior is a shrinkage choice rather than a considered statement about the null. It would export the tightness parameter into what reads as evidence. +- **Per-lag summaries only, no aggregate** — rejected: a `p`-lag block needs a single number to be comparable across pairs and models, and users would compute one anyway (usually badly, from the per-lag medians rather than draw by draw). + +## Consequences + +- `p_rope` is `None` unless a `rope` is given. Downstream code must handle that; there is no substitute value. +- The class docstring carries the full statement of what `p_rope` is not, and the how-to page has a section devoted to it. This is documentation load the feature cannot shed: the wrong reading of the number is the natural one. +- The norm is scale-dependent, so `standardize=True` is the default and the applied factor is recorded on the result. Under Toda-Yamamoto augmentation the fit is in levels and the sample standard deviations carry the series' trends, so standardised magnitudes compare within a fit rather than across fits — documented on the result and in the how-to. +- Toda-Yamamoto metadata is kept separable rather than collapsed: `n_lags_tested` and `n_lags_fitted` are distinct fields, and `augmentation_source` records whether the augmentation came from the analyst or from the integration-order diagnostics (including when those returned `d_max = 0`). A procedure that silently tested the lags it happened to fit would be a different, invalid test. +- `toda_yamamoto` refuses to run on inconclusive diagnostics instead of warning. `d_max` is a floor whenever `IntegrationOrderResult.inconclusive` is non-empty, and an under-augmented test is invalid rather than imprecise, so there is nothing useful to return. diff --git a/docs/how-to/climate-pitfalls.md b/docs/how-to/climate-pitfalls.md index 89217f4..70a543c 100644 --- a/docs/how-to/climate-pitfalls.md +++ b/docs/how-to/climate-pitfalls.md @@ -160,6 +160,9 @@ them the way you would record the lag order: without its `alpha` and its `regression` is not reproducible. - **`d_max`**, the highest integration order in the system. It is the augmentation term a Toda-Yamamoto style procedure needs, and it is far - easier to carry forward now than to reconstruct later. + easier to carry forward now than to reconstruct later. See [Granger + causality and Toda-Yamamoto](granger-causality.md) for the procedure that + consumes it — including why it refuses to run when `inconclusive` is + non-empty. - The **cointegration rank** and the lag order it was conditioned on. The rank is not invariant to `k_ar_diff`. diff --git a/docs/how-to/granger-causality.md b/docs/how-to/granger-causality.md new file mode 100644 index 0000000..c58ef55 --- /dev/null +++ b/docs/how-to/granger-causality.md @@ -0,0 +1,238 @@ +# Granger Causality and Toda-Yamamoto + +Granger causality asks whether one variable's past improves the prediction +of another beyond that other variable's own past. Impulso answers it as a +posterior over a *magnitude* rather than as a hypothesis test, and it will +not hand you a probability that there is no causality — the last section +explains why not. + +## Query a fitted model + +The query lives on the reduced-form object. No identification scheme is +involved: it reads the reduced-form coefficients directly. + +```python +from impulso import VAR, VARData + +data = VARData.from_df(df, endog=["co2", "temperature"]) +fitted = VAR(lags=2).fit(data) + +result = fitted.granger_causality("co2", "temperature") +result.summary() +``` + +The pair is ordered. `granger_causality("co2", "temperature")` tests the +lags of carbon dioxide in the temperature equation; swapping the arguments +asks the other question, and the two answers are unrelated. + +## Read the summary + +``` + median hdi_lower hdi_upper +term +L1 0.31 0.18 0.44 +L2 -0.06 -0.19 0.07 +norm 0.32 0.20 0.45 +``` + +One row per tested lag, then the headline. `norm` is the Euclidean norm of +the tested coefficients, `‖b‖ = sqrt(sum_k b_k²)`, computed draw by draw — +so its posterior is a posterior for the joint strength of the whole lag +block, not a summary of the per-lag medians. `hdi_lower` and `hdi_upper` +bound the highest-density interval (HDI), 89% by default; pass +`summary(prob=0.95)` for a different mass. + +Keeping the per-lag rows matters: a strong first lag with an offsetting +second lag is a different finding from two moderate ones, and the norm alone +cannot tell them apart. + +By default the draws are standardised — multiplied by +`sd(cause) / sd(effect)`, both sample standard deviations of the estimation +data — so a magnitude reads as standard deviations of the effect per +standard deviation of the cause. The factor is on the result as `scale`. +Pass `standardize=False` for raw coefficient units. + +## Put a number on "practically zero" + +Supply a region of practical equivalence (ROPE) — the magnitude below which +you would call the relationship negligible — and the result also reports +`p_rope`: + +```python +result = fitted.granger_causality("co2", "temperature", rope=0.05) +result.p_rope # P(||b|| < 0.05 | data) +``` + +There is deliberately no default. The ROPE is where your judgement about +what counts as a small effect enters, and it belongs in the write-up next to +the number it produced. In standardised units it is read as "a shift of one +standard deviation in the cause moves the effect by less than `rope` +standard deviations". + +## Why there is no probability of no causality + +The obvious thing to want is `P(no causality | data)`. Impulso does not +report it, because under the priors it fits that quantity is zero by +construction and would be zero whatever the data said. + +Every coefficient in Impulso has a continuous prior — Normal under +`MinnesotaPrior`, Normal-Inverse-Wishart under `NIWPrior`. A continuous +distribution assigns probability zero to any single point, so +`P(b = 0) = 0` before seeing the data. Conditioning cannot raise a +probability from zero. A model that can answer the question needs a prior +that puts a lump of mass on the null itself — a spike-and-slab, or an +edge-inclusion prior over which coefficients are present at all — which is a +different model, not a different summary of this one. + +So the honest reformulation is the ROPE one: not "is it exactly zero?" but +"is it smaller than I would care about?". That is what `p_rope` answers, and +it is only meaningful because you chose the threshold. Report `p_rope` +together with the `rope` that produced it; alone it is uninterpretable. + +## Toda-Yamamoto for integrated systems + +Standard Granger inference assumes the VAR's asymptotics are the stationary +ones. On integrated series they are not, and the usual fix — difference +everything first — changes the question to one about growth rates and +discards any long-run relationship. + +Toda and Yamamoto (1995) offer a way around it: fit the VAR in levels with +`p + d` lags, where `p` is the lag order you would have chosen and `d` the +highest integration order in the system, then test only the first `p`. The +extra `d` lags are not part of the hypothesis. They exist to restore the +standard asymptotics. + +```python +from impulso import toda_yamamoto + +result = toda_yamamoto(data, "co2", "temperature", lags=2, rope=0.05) + +result.n_lags_tested # 2 — what the answer is about +result.n_lags_fitted # 3 — what was estimated +result.augmentation # 1 +result.augmentation_source # "integration_order" +result.integration_order_result.summary() +``` + +The test lag order is never silently changed to match the fitted one. Both +numbers are on the result, and only the tested lags appear in `summary()`. + +`d` comes from `integration_order` unless you pass it. That call needs the +optional `diagnostics` extra (`pip install "impulso[diagnostics]"`). + +### When it refuses + +`integration_order` lists a variable in `inconclusive` when it is still +non-stationary at `max_order`, or when the Augmented Dickey-Fuller (ADF) and +Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests disagreed where the search +stopped. In that case `d_max` is a floor rather than a finding: the true +augmentation may be higher, and an under-augmented Toda-Yamamoto test is +invalid. + +Rather than guess, `toda_yamamoto` raises, naming the variables. Read the +table, decide, and pass the augmentation yourself: + +```python +from impulso import integration_order + +integration_order(data).summary() # look at every level, per variable +toda_yamamoto(data, "co2", "temperature", lags=2, d=2) +``` + +Passing `d=` skips the diagnostics entirely — the decision is recorded as +`augmentation_source="user"` — so it also works without `statsmodels` +installed. `d=0` is legitimate: it is the plain Granger test, and it is what +the diagnostics themselves return for a stationary system. + +If you already ran the diagnostics, hand them over instead of paying for +them twice: + +```python +diagnostics = integration_order(data) +toda_yamamoto(data, "co2", "temperature", lags=2, integration_order_result=diagnostics) +``` + +### The manual route + +`toda_yamamoto` fits with the conjugate estimator, which draws in closed +form — augmentation inflates the lag order, and this keeps that cheap. It +therefore does not accept exogenous regressors, the NUTS estimator, or a +stochastic-volatility process. For any of those, run the same three steps by +hand: + +```python +from impulso import VAR, integration_order + +d = integration_order(data).d_max # check .inconclusive first +fitted = VAR(lags=2 + d).fit(data) +fitted.granger_causality("co2", "temperature", test_lags=2) +``` + +`test_lags=2` is the whole trick: 2 + `d` lags are estimated, two are +tested. The result records the untested lags as `augmentation`. + +## A worked example, and what it does not license + +Take deseasonalised Mauna Loa carbon dioxide and a global mean surface +temperature series, both annual, both in levels: + +```python +from impulso import VARData, integration_order, toda_yamamoto + +data = VARData.from_df(climate_df, endog=["co2", "temperature"]) +integration_order(data).summary() # both I(1), typically + +forward = toda_yamamoto(data, "co2", "temperature", lags=2, rope=0.05) +reverse = toda_yamamoto(data, "temperature", "co2", lags=2, rope=0.05) + +forward.median(), forward.p_rope +reverse.median(), reverse.p_rope +``` + +Suppose the forward direction comes back with a large norm and a `p_rope` +near zero, and the reverse with a small one. Here is precisely what may be +said: *past carbon dioxide improves the prediction of temperature beyond +temperature's own past, in this bivariate system, at these lags.* Nothing +more. In particular: + +**Granger causality is predictive precedence, not intervention.** It ranks +information sets, not policies. It cannot tell you what temperature would do +under a counterfactual emissions path — that is what +`counterfactual` and `structural_scenario` are for, and they need an +identification scheme. + +**The bivariate system omits the other forcings.** Solar variability, +volcanic and anthropogenic aerosols, and El Niño-Southern Oscillation all +drive temperature and are correlated with the industrial era. A driver +omitted from the system can manufacture apparent causality between the two +variables that remain, or mask it. + +**The physical coupling runs both ways.** Carbon dioxide forces temperature +radiatively; temperature drives carbon dioxide back through ocean solubility +and the response of respiration and the terrestrial carbon sink. A test that +finds one direction "stronger" has found something about the sampling +frequency and the lag structure, not about which mechanism is real. + +**Aggregation distorts lead-lag structure.** Annual means, and the smoothing +inherent in ice-core and other proxy records, compress the timescales the +test measures. A mechanism that operates within a year is invisible to +annual data, and the smoothing can shift apparent leads by whole periods. + +The general warnings about unit-root testing on climate series apply here +too, because Toda-Yamamoto consumes an integration order: see +[Stationarity pitfalls in climate data](climate-pitfalls.md). + +## What to record + +- The **ordered pair** and the **direction**, both ways round if you ran + both. "X and Y are Granger-related" is not a result. +- The **`rope`** alongside every `p_rope`. The number means nothing without + its threshold. +- **`n_lags_tested` and `n_lags_fitted`**, not just the lag order. Under + augmentation these differ, and the difference is the point. +- **`augmentation_source`**, and the integration-order table when the + diagnostics were consulted. +- Whether the magnitudes are **standardised**, and the `scale` if not + obvious. Under lag augmentation the model is in levels, so the standard + deviations carry the series' trends and magnitudes compare best within one + fit. diff --git a/docs/how-to/index.md b/docs/how-to/index.md index 35797ce..c023393 100644 --- a/docs/how-to/index.md +++ b/docs/how-to/index.md @@ -10,6 +10,7 @@ custom-priors lag-selection stationarity-testing climate-pitfalls +granger-causality predictive-checks sign-restrictions long-run-restrictions diff --git a/docs/reference/causality.md b/docs/reference/causality.md new file mode 100644 index 0000000..8035b3f --- /dev/null +++ b/docs/reference/causality.md @@ -0,0 +1,31 @@ +# Granger Causality + +Directional predictive-strength queries on a fitted VAR. The per-pair query +lives on `FittedVAR.granger_causality`; `toda_yamamoto` is the lag-augmented +entry point for possibly-integrated systems. + +`toda_yamamoto` resolves its augmentation from the integration-order +diagnostics unless you pass `d=` yourself, so that path needs the optional +`diagnostics` extra: + +``` +pip install "impulso[diagnostics]" +``` + +```{eval-rst} +.. currentmodule:: impulso + +.. autosummary:: + :toctree: generated/ + :nosignatures: + + toda_yamamoto +``` + +Both entry points return a +{class}`~impulso.results.GrangerCausalityResult`, documented on the +[Results](results.md) page — read its docstring for exactly what `p_rope` +claims and what it does not. + +See [Granger causality and Toda-Yamamoto](../how-to/granger-causality.md) +for the recipes. diff --git a/docs/reference/index.md b/docs/reference/index.md index 100e069..65141cd 100644 --- a/docs/reference/index.md +++ b/docs/reference/index.md @@ -9,6 +9,7 @@ source docstrings by Sphinx autodoc. data spec stationarity +causality priors conjugate volatility diff --git a/docs/reference/results.md b/docs/reference/results.md index f0946bb..a1e7878 100644 --- a/docs/reference/results.md +++ b/docs/reference/results.md @@ -21,6 +21,7 @@ StationarityTestResult CointegrationTestResult IntegrationOrderResult + GrangerCausalityResult SVForecastResult VolatilityResult ``` diff --git a/src/impulso/_granger.py b/src/impulso/_granger.py index 2e4dbb9..a424738 100644 --- a/src/impulso/_granger.py +++ b/src/impulso/_granger.py @@ -244,8 +244,8 @@ def _resolve_augmentation( raise ValueError( f"the integration order of {consulted.inconclusive} is unsettled: each of these is either " f"still non-stationary at max_order={consulted.max_order} (so its recorded order is a floor, " - "not a finding) or had ADF and KPSS disagree where the search stopped. d_max would then " - "under-augment, and under-augmented Toda-Yamamoto inference is invalid. Inspect the full " + "not a finding) or had the two unit-root pretests disagree where the search stopped. d_max " + "would then under-augment, and under-augmented Toda-Yamamoto inference is invalid. Inspect the full " "table with integration_order(...).summary(), then pass the augmentation explicitly as " "d= once you have decided." ) From 4450735954da7c60c256e79c73caf22c9254deef Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 09:33:17 +0200 Subject: [PATCH 4/6] test(granger): pin the posterior dim-realignment fallbacks (#154) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A transposed posterior must be realigned by its canonical dim names, and an unlabelled one must fall back to the positional (chain, draw, var, coeff) convention — the same contract `dynamic_multiplier` relies on. Neither path was exercised. Takes `_granger.py` to full branch coverage. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Egjd7ToFeb9TQqFnfRQZxV --- tests/test_granger.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/tests/test_granger.py b/tests/test_granger.py index a56685e..1cdcb5d 100644 --- a/tests/test_granger.py +++ b/tests/test_granger.py @@ -169,6 +169,23 @@ def test_agrees_with_the_lag_matrices_split(self, hand_built): result = hand_built.granger_causality("y2", "y1", standardize=False) np.testing.assert_allclose(result.coef_draws, expected) + def test_transposed_posterior_is_realigned_by_name(self, hand_built): + # Hand-built posteriors may order their dims arbitrarily; the + # canonical labels are enough to put them back. + transposed = hand_built.idata.posterior["B"].transpose("var", "coeff", "chain", "draw") + hand_built.idata.posterior["B"] = transposed + result = hand_built.granger_causality("y2", "y1", standardize=False) + np.testing.assert_allclose(result.coef_draws, np.broadcast_to([0.12, 0.14], (2, 10, 2))) + + def test_unlabelled_posterior_falls_back_to_the_positional_convention(self): + # No canonical dim names at all — trust (chain, draw, var, coeff), + # the same contract as `dynamic_multiplier`. + fitted = _hand_built_fitted(B_2V_2L) + values = fitted.idata.posterior["B"].values + fitted.idata.posterior["B"] = xr.DataArray(values, dims=["a", "b", "c", "d"]) + result = fitted.granger_causality("y2", "y1", standardize=False) + np.testing.assert_allclose(result.coef_draws, np.broadcast_to([0.12, 0.14], (2, 10, 2))) + def test_three_variable_system_picks_the_right_pair(self): fitted = _hand_built_fitted(B_3V_2L) result = fitted.granger_causality("y3", "y2", standardize=False) From 26ce6ed70635825608f663552d61b1c3f64edb5a Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 09:44:59 +0200 Subject: [PATCH 5/6] docs(granger): state p_rope as the probability of the event, not the event The honesty paragraph's one imprecise sentence read the probability as an assertion of practical negligibility; a p_rope of 0.02 says no such thing. Refs #154 Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Egjd7ToFeb9TQqFnfRQZxV --- src/impulso/results.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/impulso/results.py b/src/impulso/results.py index bd4778e..1221536 100644 --- a/src/impulso/results.py +++ b/src/impulso/results.py @@ -800,9 +800,9 @@ class GrangerCausalityResult(ImpulsoBaseModel): both before and after seeing the data, so no dataset can raise it — a genuine posterior probability of exact non-causality needs a prior that puts point mass on the null (spike-and-slab / edge inclusion), which - Impulso does not fit. What `p_rope` does say is that the tested - coefficients are jointly *practically* negligible at the magnitude you - declared negligible. Choosing `rope` is the analyst's job and there is + Impulso does not fit. What `p_rope` quantifies is the posterior + probability that the tested coefficients are jointly *practically* + negligible at the magnitude you declared. Choosing `rope` is the analyst's job and there is no default, because there is no data-free notion of "small enough"; that the choice is explicit and recorded is the honesty of the statement. With `rope=None` the result reports the distribution only From f69edaf51f98f993e62f8f9b1b68d9bb8754a731 Mon Sep 17 00:00:00 2001 From: Thomas Pinder Date: Wed, 29 Jul 2026 15:51:55 +0200 Subject: [PATCH 6/6] test(granger): use a non-constant exog fixture (#237 rejects constant columns) Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Egjd7ToFeb9TQqFnfRQZxV --- tests/test_granger.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_granger.py b/tests/test_granger.py index 1cdcb5d..0c34728 100644 --- a/tests/test_granger.py +++ b/tests/test_granger.py @@ -442,7 +442,7 @@ def test_exogenous_data_points_at_the_manual_route(self): with_exog = VARData( endog=np.asarray(data.endog), endog_names=list(data.endog_names), - exog=np.ones((len(data.index), 1)), + exog=np.arange(len(data.index), dtype=float).reshape(-1, 1), exog_names=["trend"], index=data.index, )