diff --git a/CONTEXT.md b/CONTEXT.md index f62ffad1..cd7cea24 100644 --- a/CONTEXT.md +++ b/CONTEXT.md @@ -22,6 +22,12 @@ The sum of the moving-average coefficients of a stable reduced-form VAR, `C(1) = **Cumulative MA impact matrix (Θ(1))**: The structural counterpart, `Θ(1) = C(1) P`: the total effect of each structural shock on each variable's level. Where `Cholesky` and `SignRestriction` constrain the impact matrix `Θ(0) = P`, `LongRunRestriction` constrains `Θ(1)` to be lower-triangular in a stated variable ordering, with positive diagonal (shock `j` raises variable `j`'s long-run level). _Avoid_: "Blanchard-Quah identification" as an API name — the scheme is named for what it restricts, not for its first users (the prose citation is fine). "Permanent shock" for anything but the first column: only the shock restricted nowhere is unambiguously permanent. +A rule for recovering structural shocks from reduced-form covariance, implemented as adapters of the `IdentificationScheme` Protocol (`Cholesky`, `SignRestriction`, `ZeroSignRestriction`, `ProxySVAR`). The scheme is a pure function: it consumes a Cholesky factor `L` and produces a structural shock matrix `B = identify(L)`. It does not own time iteration. +_Avoid_: "identification strategy" (used colloquially; the Protocol is named "scheme"). + +**Zero-and-sign restrictions (ARW construction)**: +Identification combining exact zeros on the impact matrix with sign restrictions on impulse responses, via the recursive orthogonalisation of Arias, Rubio-Ramírez & Waggoner (2018). With `P = L Q`, a zero is a *linear* condition on one column of `Q` (`Z_j L q_j = 0`), so columns are built one at a time — each drawn uniformly from the unit sphere of the null space of the stacked zero conditions and the already-drawn columns. The zeros hold by construction; only the signs are accept/reject. Admissibility is the Rubio-Ramírez, Waggoner & Zha (2010) counting condition `z_j <= n - j` on shocks sorted by zero count, checked before sampling. Equality throughout is exact identification and reproduces the Cholesky factor; anything looser is set identification. Draws are *unweighted* — no volume-element correction for the ARW uniform-conditional prior — which is documented rather than silent. Failed draws are `NaN`, never a fallback to `L`, because a fallback would violate the zeros the scheme exists to impose. +_Avoid_: "penalty-function zeros" — that is a different (Uhlig-style loss-minimising) construction that only approximates the zeros; these are exact. **Volatility process**: The seam that owns how the structural-shock covariance Σ_t is constructed (in PyMC), evolved over time, and queried. Concrete adapters of the `VolatilityProcess` Protocol: `Constant` (homoscedastic Σ; the default) and `StochasticVolatility` (time-varying). The volatility process owns its downstream computation — forecast covariance paths, time-`t` Cholesky query, per-variable volatility paths — so the pipeline never branches on adapter type. diff --git a/docs/explanation/identification.md b/docs/explanation/identification.md index c74c54fe..69bfa2e1 100644 --- a/docs/explanation/identification.md +++ b/docs/explanation/identification.md @@ -72,3 +72,52 @@ $C(1)$ is the long-run multiplier on the levels of the variables *as modelled*. Triangularity is $n(n-1)/2$ restrictions, exactly the number needed for point identification. With two variables that is the one restriction you wanted. With three it asserts three zeros at once, which is a much stronger joint claim than it looks. Arbitrary (non-recursive) long-run zero patterns are not supported. Two things can go wrong, and Impulso reports them separately. If $I - \sum_j A_j$ is close to singular, $C(1)$ is numerically undefined and those draws are blanked (or, with `on_undefined="raise"`, refused). If a draw is explosive — companion spectral radius above one — the arithmetic is fine but $\sum_h \Phi_h$ diverges, so "long-run effect" has no meaning for it; those draws are always reported and never blanked, because posteriors near a unit root routinely contain them. +## Zero and sign restrictions together + +Cholesky and sign restrictions sit at opposite ends of one spectrum. Cholesky imposes $n(n-1)/2$ zeros, which is exactly enough to pin a single answer. Sign restrictions impose no zeros at all and leave a set of answers. Most credible identification schemes sit in between: a few zeros you can defend on institutional or physical grounds, plus signs on the responses you have a firm prior about. + +`ZeroSignRestriction` covers that middle ground, following {cite:t}`ariasRubioRamirezWaggoner2018`: + +```python +from impulso.identification import ZeroSignRestriction + +scheme = ZeroSignRestriction( + shock_names=["supply", "demand", "policy"], + zero_restrictions={"gdp": ["policy"]}, + sign_restrictions={ + "gdp": {"supply": "+", "demand": "+"}, + "inflation": {"supply": "-", "demand": "+", "policy": "-"}, + }, +) +``` + +This says output does not respond to a policy shock within the period, on top of the usual signs. + +### How the construction works + +Write the structural impact matrix as $P = LQ$, where $L$ is the Cholesky factor of the residual covariance and $Q$ is orthogonal. The response of variable $i$ to shock $j$ on impact is $e_i' L q_j$, so a zero restriction is a *linear* condition on the $j$-th column of $Q$: + +$$Z_j L q_j = 0$$ (eq-zero-condition) + +where $Z_j$ selects the rows carrying zeros on shock $j$. Because {eq}`eq-zero-condition` is linear, it does not have to be searched for. Columns are built one at a time, in decreasing order of how many zeros they carry. At step $k$ the column must satisfy its own zero conditions *and* be orthogonal to the $k-1$ columns already drawn: + +$$R_k = \begin{pmatrix} Z_k L \\ q_1' \\ \vdots \\ q_{k-1}' \end{pmatrix}, \qquad q_k \sim \text{Uniform}\!\left(\mathcal{S} \cap \mathcal{N}(R_k)\right)$$ (eq-arw-recursion) + +with $\mathcal{N}(R_k)$ the null space of $R_k$ and $\mathcal{S}$ the unit sphere. Drawing a standard Gaussian in an orthonormal basis of that null space and normalising gives the uniform draw. The zeros then hold to the precision of the singular value decomposition (SVD) used to find the basis, and orthogonality holds by construction, so the accept/reject step only ever has to test the *signs*. + +### The rank condition + +$R_k$ has $z_k + (k-1)$ rows, where $z_k$ is the number of zeros on the $k$-th shock, so its null space is non-trivial only when + +$$z_j \le n - j, \qquad j = 1, \dots, n$$ (eq-rwz-rank) + +with shocks indexed in decreasing order of $z_j$. This is the counting condition of {cite:t}`rubioRamirezWaggonerZha2010`. It depends on the restriction pattern alone, so Impulso checks it once before any sampling and raises a `ValueError` naming the offending shock. Equality throughout — $z_j = n - j$ for every shock — makes every null space one-dimensional, the answer unique up to column signs, and reproduces the Cholesky factor. That is the exactly-identified end of the spectrum; anything looser leaves a set. + +:::{admonition} Draws are unweighted +:class: warning +Accepted draws are kept as the recursion produced them, with no importance weight. {cite:t}`ariasRubioRamirezWaggoner2018` derive such a weight — a volume-element correction for the manifold the zero restrictions carve out — which their uniform-conditional prior over the identified set requires. Without it, a set-identified posterior from `ZeroSignRestriction` is not that prior exactly. + +Two cases are unaffected. With no zero restrictions the recursion reduces to Gram-Schmidt on Gaussian columns, which is exactly the Haar measure. With zeros that exactly identify the system, the answer is a point up to column signs, and no weight can move a point. In between, treat the spread across draws as reflecting the restrictions and the parameter posterior, not a calibrated prior over rotations. + +Two further differences from the paper are worth stating. Impulso retries rotations *within* each posterior draw, as `SignRestriction` does, rather than rejecting the joint $(\theta, Q)$ pair; and it draws from the orthogonal group $O(n)$ rather than the special orthogonal group $SO(n)$, which is immaterial whenever at least one column's sign is left unpinned. +::: diff --git a/docs/how-to/index.md b/docs/how-to/index.md index 7ff315b0..aaa74ee2 100644 --- a/docs/how-to/index.md +++ b/docs/how-to/index.md @@ -13,4 +13,5 @@ climate-pitfalls sign-restrictions long-run-restrictions heavy-tailed-errors +zero-sign-restrictions ``` diff --git a/docs/how-to/zero-sign-restrictions.md b/docs/how-to/zero-sign-restrictions.md new file mode 100644 index 00000000..e7c3b15e --- /dev/null +++ b/docs/how-to/zero-sign-restrictions.md @@ -0,0 +1,74 @@ +# Combining Zero and Sign Restrictions + +`ZeroSignRestriction` imposes exact zeros on the impact matrix alongside sign restrictions on impulse responses. Use it when part of your identification rests on a timing or exclusion argument you can defend outright, and the rest on the direction a response should take. + +## A worked example + +Take a three-variable system: a temperature anomaly, a measure of economic activity, and emissions. Economic activity moves emissions immediately, and emissions eventually move temperature, but the thermal inertia of the ocean means an activity shock cannot register in this year's temperature anomaly. That is a zero, not a sign. + +```python +from impulso.identification import ZeroSignRestriction + +scheme = ZeroSignRestriction( + shock_names=["climate", "activity", "emissions"], + zero_restrictions={ + # An activity shock has no contemporaneous effect on temperature: + # the physical lag between emissions and warming is far longer + # than the sampling frequency. + "temperature": ["activity"], + }, + sign_restrictions={ + "temperature": {"climate": "+"}, + "activity": {"climate": "-", "activity": "+"}, + "emissions": {"activity": "+", "emissions": "+"}, + }, + random_seed=42, +) + +identified = fitted.set_identification_strategy(scheme) +irf = identified.impulse_response(horizon=20) +``` + +## Specifying restrictions + +- `shock_names` fixes the column order of the returned structural matrix. Name fewer shocks than you have variables and the rest are labelled `unidentified_1`, `unidentified_2`, ... — those columns carry no restrictions and are rotation-arbitrary, so `fevd()` masks their shares. +- `zero_restrictions` maps a variable to the shocks that do not move it on impact. Zeros bind at horizon 0 only; long-run zeros are not supported. +- `sign_restrictions` uses the same format as `SignRestriction`: variable → shock → `"+"` or `"-"`. +- `restriction_horizon=H` imposes the *signs* at horizons `0..H`. The zeros stay at impact. + +A cell cannot be restricted to zero and to a sign at once — that is a contradiction at horizon 0, and construction fails with a `ValueError`. + +## How many zeros are admissible + +Sort the shocks by how many zeros they carry, most first. The shock in position `j` may carry at most `n - j` zeros. Break that and identification is impossible for *any* orthogonal matrix, so `identify()` raises before sampling starts rather than burning through rotations. At the limit — `n - 1`, `n - 2`, ..., `0` — the zeros exactly identify the system and reproduce the Cholesky factor. + +## When draws fail + +A draw fails when no candidate satisfies the sign restrictions within `n_rotations` attempts. + +:::{admonition} Failed draws become NaN, not Cholesky +:class: warning +`ZeroSignRestriction` fills failed draws with `NaN` and warns once with the count and fraction. It deliberately does **not** fall back to the unrotated factor the way `SignRestriction` does: that fallback would silently break the zero restrictions, which are the whole point of the scheme. + +`NaN` draws propagate into impulse-response and variance-decomposition summaries and are rejected by the scenario methods, so treat a non-trivial failure fraction as a result to act on, not a warning to suppress. Raise `n_rotations`, relax the signs, or set `on_failure="raise"` to stop at the first failure. +::: + +## Tuning + +- Each unpinned column sign is effectively a coin flip, so a scheme with `k` sign-restricted shocks accepts roughly one candidate in `2**k` even when the restrictions are otherwise easy. Budget `n_rotations` accordingly. +- Read the acceptance rate off the shock matrix: + + ```python + attrs = identified.shock_matrix().attrs + attrs["zero_sign_acceptance_rate"] # fraction of draws identified + attrs["zero_sign_mean_attempts"] # candidates drawn per draw + attrs["zero_sign_max_zero_violation"] # largest |zero cell|, expect ~1e-14 + ``` + +- Set `random_seed` for reproducibility. +- Naming a shock without giving it any sign restriction leaves its column sign unidentified, so posterior summaries average over both directions. Impulso warns when this happens. + +:::{admonition} Draws are unweighted +:class: note +No importance weight corrects for the volume element of the zero-restricted manifold, so set-identified results are not the uniform-conditional prior of Arias, Rubio-Ramírez and Waggoner (2018). See [the identification explanation](../explanation/identification.md) for what this does and does not affect. +::: diff --git a/docs/reference/identification.md b/docs/reference/identification.md index 65cd7635..acde8504 100644 --- a/docs/reference/identification.md +++ b/docs/reference/identification.md @@ -11,4 +11,5 @@ SignRestriction LongRunRestriction ProxySVAR + ZeroSignRestriction ``` diff --git a/docs/references.bib b/docs/references.bib index 2601e430..9233bfc9 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -156,3 +156,23 @@ @article{blanchardQuah1989 number = {4}, pages = {655--673}, } + +@article{ariasRubioRamirezWaggoner2018, + author = {Arias, Jonas E. and Rubio-Ram{\'i}rez, Juan F. and Waggoner, Daniel F.}, + title = {Inference Based on Structural Vector Autoregressions Identified with Sign and Zero Restrictions: Theory and Applications}, + journal = {Econometrica}, + year = {2018}, + volume = {86}, + number = {2}, + pages = {685--720}, +} + +@article{rubioRamirezWaggonerZha2010, + author = {Rubio-Ram{\'i}rez, Juan F. and Waggoner, Daniel F. and Zha, Tao}, + title = {Structural Vector Autoregressions: Theory of Identification and Algorithms for Inference}, + journal = {The Review of Economic Studies}, + year = {2010}, + volume = {77}, + number = {2}, + pages = {665--696}, +} diff --git a/src/impulso/__init__.py b/src/impulso/__init__.py index 8fe3f846..b1999bb1 100644 --- a/src/impulso/__init__.py +++ b/src/impulso/__init__.py @@ -16,7 +16,7 @@ from impulso.conjugate_volatility import ConjugateVolatility, PandemicBreak from impulso.evidence import EvidenceComparison, ModelEvidence, compare_evidence from impulso.fitted import FittedVAR - from impulso.identification import Cholesky, LongRunRestriction, ProxySVAR, SignRestriction + from impulso.identification import Cholesky, LongRunRestriction, ProxySVAR, SignRestriction, ZeroSignRestriction from impulso.identified import IdentifiedVAR from impulso.observation import Gaussian, StudentT from impulso.priors import MinnesotaPrior, NIWPrior @@ -89,6 +89,7 @@ "VariablePath", "VolatilityProcess", "VolatilityResult", + "ZeroSignRestriction", "adf_test", "compare_evidence", "compute_ma_phi", @@ -108,6 +109,7 @@ "LongRunRestriction": "impulso.identification", "ProxySVAR": "impulso.identification", "SignRestriction": "impulso.identification", + "ZeroSignRestriction": "impulso.identification", "MinnesotaPrior": "impulso.priors", "NIWPrior": "impulso.priors", "ConjugateVAR": "impulso.conjugate", diff --git a/src/impulso/identification.py b/src/impulso/identification.py index a7488f51..1efa3f9c 100644 --- a/src/impulso/identification.py +++ b/src/impulso/identification.py @@ -1,6 +1,7 @@ """Identification schemes for structural VAR analysis.""" import weakref +from dataclasses import dataclass from typing import TYPE_CHECKING, Any, ClassVar, Final, Literal import numpy as np @@ -1037,3 +1038,610 @@ def _first_stage_f(z_c: np.ndarray, u_policy_c: np.ndarray) -> np.ndarray: def shock_coords(self, n_vars: int) -> list[str]: """Identified shock first, then rotation-arbitrary padding.""" return SignRestriction._build_shock_coords([self.shock_name], n_vars) + + +# -------------------------------------------------------------------------- +# Zero-and-sign restrictions (Arias, Rubio-Ramirez & Waggoner, 2018) +# -------------------------------------------------------------------------- + +#: Prefix reserved by the pipeline for rotation-arbitrary shock columns. +_RESERVED_SHOCK_PREFIX = "unidentified_" + + +def _signs_ok(values: np.ndarray, rows: np.ndarray, cols: np.ndarray, signs: np.ndarray) -> bool: + """Check sign restrictions against a response matrix. + + Restrictions are supplied as flat index arrays so the check is one + fancy-index plus one comparison, rather than a Python loop over the + nested restriction dict. + + Args: + values: Response matrix, shape `(n_vars, n_vars)`; rows are + variables, columns are shocks (in data / user order). + rows: Variable indices of the restricted cells. + cols: Shock indices of the restricted cells. + signs: `+1.0` for a `"+"` restriction, `-1.0` for `"-"`. + + Returns: + True if every restricted cell has the required sign. Exact zeros + pass, matching `SignRestriction._check_restrictions`. + """ + if rows.size == 0: + return True + return bool(np.all(values[rows, cols] * signs >= 0.0)) + + +def _signs_ok_at_horizons( + P: np.ndarray, + Phi: np.ndarray, + horizon: int, + rows: np.ndarray, + cols: np.ndarray, + signs: np.ndarray, +) -> bool: + """Check sign restrictions at horizons `1..horizon`. + + Impact (`h = 0`) is *not* checked here: the recursive construction + screens impact signs column by column while it builds `P`, so by the + time this runs the impact restrictions already hold. + + Args: + P: Candidate structural impact matrix, shape `(n_vars, n_vars)`. + Phi: MA coefficient matrices for this draw, shape + `(horizon + 1, n_vars, n_vars)`. + horizon: Highest horizon to check. + rows: Variable indices of the restricted cells. + cols: Shock indices of the restricted cells. + signs: `+1.0` for a `"+"` restriction, `-1.0` for `"-"`. + + Returns: + True if all restrictions hold at every horizon `1..horizon`. + """ + return all(_signs_ok(Phi[h] @ P, rows, cols, signs) for h in range(1, horizon + 1)) + + +@dataclass(frozen=True) +class _ZeroSignLayout: + """Draw-independent bookkeeping compiled once per `identify()` call. + + Two column orderings coexist and mixing them is the easiest way to get + this wrong, so they are named explicitly here: *user order* is the order + of `shock_coords` (what callers see), *construction order* is the + zero-count-descending order the recursion requires. + + Attributes: + labels: Effective shock labels, user order. + order: `order[k]` is the user column built at construction step `k`. + perm_back: Inverse of `order` — permutes construction-order columns + back to user order. + zero_rows: Zero-restricted variable rows per construction step. + sign_rows: Variable indices of sign-restricted cells. + sign_cols: Shock indices (user order) of sign-restricted cells. + sign_vals: `+1.0` / `-1.0` targets matching `sign_rows`/`sign_cols`. + col_rows: Impact-sign-restricted variable rows per construction step. + col_signs: Matching sign targets per construction step. + zero_idx_rows: Variable indices of zero-restricted cells (user order). + zero_idx_cols: Shock indices of zero-restricted cells (user order). + """ + + labels: list[str] + order: list[int] + perm_back: np.ndarray + zero_rows: list[np.ndarray] + sign_rows: np.ndarray + sign_cols: np.ndarray + sign_vals: np.ndarray + col_rows: list[np.ndarray] + col_signs: list[np.ndarray] + zero_idx_rows: np.ndarray + zero_idx_cols: np.ndarray + + +class ZeroSignRestriction(ImpulsoModel): + """Combined zero-and-sign restriction identification. + + Implements the recursive orthogonalisation of Arias, + Rubio-Ramirez & Waggoner (2018). Writing the structural impact matrix + as `P = L Q` with `Q` orthogonal, a zero restriction "variable `i` + does not respond to shock `j` on impact" is the linear condition + `e_i' L q_j = 0` on the `j`-th column of `Q`. Columns are built one at + a time, each drawn uniformly from the unit sphere of the null space of + + R_k = [ Z_k L ; q_1' ; ... ; q_{k-1}' ] + + where `Z_k` selects the rows carrying zero restrictions on shock `k`. + The null-space draw imposes the zeros *exactly* (to SVD precision) and + orthogonality to the earlier columns by construction, so no rejection + step is needed for the zeros — only the sign restrictions are checked + by accept/reject. + + Shocks are ordered internally by their number of zero restrictions, + descending (ties keep the order given in `shock_names`; unnamed padding + columns go last), because the construction requires it. Rows of the + returned matrix are always in data order and columns are permuted back + to `shock_names` order, so the internal ordering is not observable. + + Attributes: + shock_names: Structural shock labels, in the order the columns of + the returned matrix should appear. May be shorter than the + number of variables — remaining columns are labelled + `unidentified_1`, ... and carry no restrictions. + zero_restrictions: Dict mapping variable -> list of shocks that + have zero impact on that variable. Keyed by variable for + consistency with `sign_restrictions`. + sign_restrictions: Dict mapping variable -> {shock: "+" or "-"}, + the same format `SignRestriction` uses. + restriction_horizon: Sign restrictions are imposed at horizons + `0..restriction_horizon`. Zero restrictions are always impact + only (`h = 0`); long-run zeros are not supported. + n_rotations: Maximum candidate draws per posterior draw. + random_seed: Seed for reproducibility. + on_failure: What to do for a posterior draw where no candidate + satisfies the sign restrictions within `n_rotations` attempts. + `"nan"` (default) fills that draw with NaN and warns once at + the end; `"raise"` raises immediately. + + Note: + Candidates are drawn *unweighted*: each accepted draw keeps the + `Q` that the recursion produced, with no importance weight + correcting for the volume element of the zero-restricted manifold. + Arias, Rubio-Ramirez & Waggoner (2018) derive such a weight for + their uniform-conditional prior over the identified set. The + unweighted draws therefore do not represent that prior exactly + when the restrictions leave a set (rather than a point) identified. + Two regimes are unaffected: with no zero restrictions the draws are + exactly Haar, and when the zeros exactly identify the system + (`z_j = n - j` for every shock) the answer is a point up to column + signs. See the explanation page for the full caveat. + """ + + shock_names: list[str] + zero_restrictions: dict[str, list[str]] = Field(default_factory=dict) + sign_restrictions: dict[str, dict[str, str]] = Field(default_factory=dict) + restriction_horizon: int = Field(default=0, ge=0) + n_rotations: int = Field(default=1000, ge=1) + random_seed: int | None = None + on_failure: Literal["nan", "raise"] = "nan" + + # Draws a fresh Q per identify() call, exactly as SignRestriction does — + # forecast-side scenario machinery reads this flag and refuses + # time-varying volatility for such schemes. + _samples_rotations: ClassVar[bool] = True + + # Single-call scratchpad, mirroring ProxySVAR._last_diagnostics: identify() + # writes, IdentifiedVAR.shock_matrix reads it back immediately and attaches + # the entries to the shock-matrix attrs. Not reentrant. + # + # Deliberately *not* `_last_acceptance_rate`: the pipeline surfaces that + # private attribute under the name `sign_restriction_acceptance_rate`, + # which would mislabel this scheme's diagnostics. + _last_diagnostics: dict[str, float] = PrivateAttr(default_factory=dict) + + @model_validator(mode="after") + def _validate_restrictions(self) -> "ZeroSignRestriction": + """Validate everything that does not depend on the number of variables. + + Variable names and the rank condition need `n_vars`, which is only + known at `identify()` time; those are checked there. + + Returns: + The validated instance. + + Raises: + ValueError: On empty/duplicate shock names, reserved shock + labels, unknown shocks, bad sign tokens, a cell restricted + to zero *and* to a sign, or two empty restriction dicts. + """ + import warnings + + self._validate_shock_names() + self._validate_restriction_dicts() + + signed_shocks = {s for signed in self.sign_restrictions.values() for s in signed} + unsigned = [s for s in self.shock_names if s not in signed_shocks] + if unsigned: + warnings.warn( + f"Named shock(s) {unsigned} carry no sign restriction at any horizon, so their " + "column sign is not identified: q and -q are both admissible and posterior " + "summaries will mix the two directions. Add a sign restriction to pin the " + "direction.", + UserWarning, + stacklevel=2, + ) + return self + + def _validate_shock_names(self) -> None: + """Check `shock_names` is a non-empty list of unreserved, unique labels. + + Raises: + ValueError: On an empty list, duplicates, or a reserved prefix. + """ + if not self.shock_names: + raise ValueError("shock_names must name at least one structural shock.") + + duplicates = sorted({s for s in self.shock_names if self.shock_names.count(s) > 1}) + if duplicates: + raise ValueError(f"Duplicate shock names in shock_names: {duplicates}.") + + reserved = [s for s in self.shock_names if s.startswith(_RESERVED_SHOCK_PREFIX)] + if reserved: + raise ValueError( + f"Shock names {reserved} use the reserved prefix {_RESERVED_SHOCK_PREFIX!r}. " + "The pipeline assigns that prefix to rotation-arbitrary columns under partial " + "identification; pick a different label." + ) + + def _validate_restriction_dicts(self) -> None: + """Check the two restriction dicts are non-empty, well-formed, and consistent. + + Raises: + ValueError: If both dicts are empty, a shock is unknown, a sign + token is not `"+"`/`"-"`, or a cell carries a zero and a + sign at once. + """ + if not self.zero_restrictions and not self.sign_restrictions: + raise ValueError( + "ZeroSignRestriction needs at least one of zero_restrictions or sign_restrictions. " + "With neither, every orthogonal Q is admissible and nothing is identified." + ) + + known = set(self.shock_names) + for variable, shocks in self.zero_restrictions.items(): + unknown = sorted(set(shocks) - known) + if unknown: + raise ValueError( + f"zero_restrictions[{variable!r}] references unknown shock(s) {unknown}. " + f"Known shocks: {self.shock_names}." + ) + clash = sorted(set(shocks) & set(self.sign_restrictions.get(variable, {}))) + if clash: + raise ValueError( + f"Restriction conflict on variable {variable!r}: shock(s) {clash} are " + "restricted to zero on impact and simultaneously given a sign. A signed " + "response contradicts a zero response at h = 0." + ) + + for variable, signed in self.sign_restrictions.items(): + unknown = sorted(set(signed) - known) + if unknown: + raise ValueError( + f"sign_restrictions[{variable!r}] references unknown shock(s) {unknown}. " + f"Known shocks: {self.shock_names}." + ) + bad = sorted({s for s, direction in signed.items() if direction not in ("+", "-")}) + if bad: + raise ValueError( + f"sign_restrictions[{variable!r}] has non-sign token(s) for shock(s) {bad}. Use '+' or '-'." + ) + + def identify( + self, + L: np.ndarray, + var_names: list[str], + posterior: "xr.Dataset | None" = None, + data: "VARData | None" = None, + n_lags: int | None = None, + ) -> np.ndarray: + """Apply zero-and-sign-restriction identification. + + Args: + L: Lower-triangular Cholesky factor, shape (chains, draws, n_vars, n_vars). + var_names: Variable names in the data's natural order. + posterior: Required when `self.restriction_horizon > 0`, which + needs the VAR coefficients `B` for the MA recursion. + Ignored for impact-only restrictions. + data: Unused. Accepted for Protocol uniformity. + n_lags: Unused — the lag order is read off `B`. Accepted for + Protocol uniformity. + + Returns: + Structural shock matrix, shape (chains, draws, n_vars, n_vars), + with columns in `shock_names` order (padding last). Draws where + no candidate satisfied the sign restrictions are NaN. Acceptance + diagnostics land on `IdentifiedVAR.shock_matrix()` attrs under + the `zero_sign_` prefix. + + Raises: + ValueError: If a restriction names an unknown variable, more + shocks are named than there are variables, the zero pattern + violates the rank condition, `restriction_horizon > 0` + without a posterior, or `on_failure="raise"` and a draw + found no admissible candidate. + """ + del data, n_lags # unused; lag order comes from B + import warnings + + n_chains, n_draws, n_vars, _ = L.shape + layout = self._compile_layout(var_names, n_vars) + B_all = self._require_coefficients(posterior) + horizon = self.restriction_horizon + + rng = np.random.default_rng(self.random_seed) + eye = np.eye(n_vars) + P = np.full((n_chains, n_draws, n_vars, n_vars), np.nan) + n_total = n_chains * n_draws + n_accepted = 0 + total_attempts = 0 + max_zero_violation = 0.0 + + n_lags_b = B_all.shape[-1] // n_vars if B_all is not None else 0 + for c in range(n_chains): + for d in range(n_draws): + # Hoisted out of the candidate loop: the MA coefficients depend + # on the draw only, not on the rotation being tried. + Phi = compute_ma_phi(lag_matrices(B_all[c, d], n_lags_b), horizon) if B_all is not None else None + accepted, attempts = self._identify_draw(L[c, d], Phi, layout, eye, rng, n_vars) + total_attempts += attempts + + if accepted is None: + if self.on_failure == "raise": + raise ValueError( + f"No admissible rotation for draw (chain={c}, draw={d}) within " + f"n_rotations={self.n_rotations}. Increase n_rotations, relax the sign " + "restrictions, or set on_failure='nan' to keep the draw as NaN." + ) + continue + + P[c, d] = accepted + n_accepted += 1 + if layout.zero_idx_rows.size: + violation = float(np.max(np.abs(accepted[layout.zero_idx_rows, layout.zero_idx_cols]))) + max_zero_violation = max(max_zero_violation, violation) + + failed = n_total - n_accepted + self._last_diagnostics = { + "zero_sign_acceptance_rate": n_accepted / n_total, + "zero_sign_failed_draws": float(failed), + "zero_sign_failed_fraction": failed / n_total, + "zero_sign_mean_attempts": total_attempts / n_total, + "zero_sign_max_zero_violation": max_zero_violation, + } + if failed: + warnings.warn( + f"Sign restrictions not satisfied for {failed}/{n_total} draws " + f"({failed / n_total:.1%}); those draws are NaN. Increase n_rotations or relax the " + "restrictions. NaN draws propagate into IRF/FEVD summaries and are rejected by the " + "scenario methods.", + UserWarning, + stacklevel=2, + ) + return P + + def _require_coefficients(self, posterior: "xr.Dataset | None") -> np.ndarray | None: + """Fetch the `B` draws when horizon restrictions need them. + + Args: + posterior: Posterior Dataset, or None. + + Returns: + The `B` draws, or None for impact-only restrictions. + + Raises: + ValueError: If `restriction_horizon > 0` and `B` is unavailable. + """ + if self.restriction_horizon == 0: + return None + if posterior is None or "B" not in posterior: + raise ValueError( + "restriction_horizon > 0 requires the full posterior with 'B' " + "(VAR coefficients). Pass posterior=fitted.idata.posterior to identify()." + ) + return posterior["B"].values + + def _compile_layout(self, var_names: list[str], n_vars: int) -> _ZeroSignLayout: + """Resolve names to indices and fix the construction order, once per call. + + Everything here depends on the restriction pattern and the variable + names only — not on the posterior draw — so it is hoisted out of the + draw loop. The rank condition is checked here too, before any + sampling happens. + + Args: + var_names: Variable names in data order. + n_vars: Number of endogenous variables. + + Returns: + The compiled layout. + + Raises: + ValueError: If more shocks are named than there are variables, a + restriction names an unknown variable, or the zero pattern + violates the rank condition. + """ + if len(self.shock_names) > n_vars: + raise ValueError( + f"ZeroSignRestriction names {len(self.shock_names)} shocks " + f"({self.shock_names}) but the VAR has only {n_vars} variables. " + "A VAR admits at most n_vars structural shocks." + ) + + referenced = set(self.zero_restrictions) | set(self.sign_restrictions) + unknown_vars = sorted(referenced - set(var_names)) + if unknown_vars: + raise ValueError(f"Restrictions reference unknown variable(s) {unknown_vars}. Variables: {var_names}.") + + labels = self.shock_coords(n_vars) + label_index = {label: j for j, label in enumerate(labels)} + var_index = {name: i for i, name in enumerate(var_names)} + + # Zero-restricted variable rows per shock column, in user order. + zero_rows_user: list[list[int]] = [[] for _ in range(n_vars)] + for variable, shocks in self.zero_restrictions.items(): + i = var_index[variable] + for shock in shocks: + col = label_index[shock] + if i not in zero_rows_user[col]: + zero_rows_user[col].append(i) + zero_rows_user = [sorted(rows) for rows in zero_rows_user] + + # Construction order: most-restricted shock first. Python's sort is + # stable, so ties keep user order and the unrestricted padding + # columns (z = 0) stay last. + order = sorted(range(n_vars), key=lambda u: -len(zero_rows_user[u])) + self._check_rank_condition(order, zero_rows_user, labels, n_vars) + + # Flat sign-restriction index arrays, columns in user order. + rows_list, cols_list, signs_list = [], [], [] + for variable, signed in self.sign_restrictions.items(): + i = var_index[variable] + for shock, direction in signed.items(): + rows_list.append(i) + cols_list.append(label_index[shock]) + signs_list.append(1.0 if direction == "+" else -1.0) + sign_rows = np.asarray(rows_list, dtype=int) + sign_cols = np.asarray(cols_list, dtype=int) + sign_vals = np.asarray(signs_list, dtype=float) + + # Flat zero-cell indices (user order) for the violation diagnostic. + zero_cells = [(i, u) for u in range(n_vars) for i in zero_rows_user[u]] + + return _ZeroSignLayout( + labels=labels, + order=order, + perm_back=np.argsort(order), + zero_rows=[np.asarray(zero_rows_user[u], dtype=int) for u in order], + sign_rows=sign_rows, + sign_cols=sign_cols, + sign_vals=sign_vals, + col_rows=[sign_rows[sign_cols == u] for u in order], + col_signs=[sign_vals[sign_cols == u] for u in order], + zero_idx_rows=np.asarray([i for i, _ in zero_cells], dtype=int), + zero_idx_cols=np.asarray([u for _, u in zero_cells], dtype=int), + ) + + def _identify_draw( + self, + chol: np.ndarray, + Phi: np.ndarray | None, + layout: _ZeroSignLayout, + eye: np.ndarray, + rng: np.random.Generator, + n_vars: int, + ) -> tuple[np.ndarray | None, int]: + """Rejection-sample one posterior draw's structural impact matrix. + + Args: + chol: This draw's Cholesky factor, shape `(n_vars, n_vars)`. + Phi: MA coefficients for this draw, or None when + `restriction_horizon == 0`. + layout: Compiled restriction bookkeeping. + eye: Cached `n_vars` identity. + rng: Random generator. + n_vars: Number of endogenous variables. + + Returns: + Tuple `(P, attempts)`: the accepted impact matrix with columns in + user order, or None if the budget ran out, and the number of + candidates drawn. + """ + for attempt in range(1, self.n_rotations + 1): + candidate = self._draw_candidate(chol, layout, eye, rng, n_vars) + if candidate is None: + continue + P_cand = chol @ candidate[:, layout.perm_back] + if Phi is not None and not _signs_ok_at_horizons( + P_cand, Phi, self.restriction_horizon, layout.sign_rows, layout.sign_cols, layout.sign_vals + ): + continue + return P_cand, attempt + return None, self.n_rotations + + def _check_rank_condition( + self, + order: list[int], + zero_rows_user: list[list[int]], + labels: list[str], + n_vars: int, + ) -> None: + """Verify the Rubio-Ramirez, Waggoner & Zha (2010) rank condition. + + With shocks sorted by zero count descending, the null space at + position `j` (1-based) has dimension `n - z_j - (j - 1)`, so a + non-degenerate column requires `z_j <= n - j`. The check is + deterministic — it depends on the restriction pattern alone — so it + runs once, before any sampling. + + Args: + order: Construction order (user column indices, most-restricted first). + zero_rows_user: Zero-restricted variable rows per user column. + labels: Effective shock labels, user order. + n_vars: Number of endogenous variables. + + Raises: + ValueError: If any sorted position violates the bound. + """ + for k, u in enumerate(order): + z_k = len(zero_rows_user[u]) + bound = n_vars - k - 1 + if z_k > bound: + raise ValueError( + f"Zero restrictions violate the rank condition of Rubio-Ramirez, Waggoner & " + f"Zha (2010): shock {labels[u]!r} carries {z_k} zero restriction(s), but at " + f"position j = {k + 1} of the zero-count ordering a shock may carry at most " + f"n - j = {bound}. With more, the null space for that column is empty and no " + "orthogonal matrix satisfies the restrictions." + ) + + def _draw_candidate( + self, + chol: np.ndarray, + layout: _ZeroSignLayout, + eye: np.ndarray, + rng: np.random.Generator, + n_vars: int, + ) -> np.ndarray | None: + """Draw one candidate orthogonal matrix in construction order. + + Args: + chol: This draw's Cholesky factor, shape `(n_vars, n_vars)`. + layout: Compiled restriction bookkeeping. + eye: Cached `n_vars` identity. + rng: Random generator. + n_vars: Number of endogenous variables. + + Returns: + An orthogonal matrix whose columns are in construction order + and satisfy every zero restriction plus every *impact* sign + restriction, or None if the impact screen failed. + """ + Q = np.empty((n_vars, n_vars)) + for k in range(n_vars): + rows_k = layout.zero_rows[k] + # Rows of R: the z_k zero conditions Z_k L q = 0, plus the k - 1 + # orthogonality conditions against the columns already drawn. + m = rows_k.size + k + if m == 0: + N = eye + else: + R = np.vstack((chol[rows_k, :], Q[:, :k].T)) + # Right singular vectors beyond the m-th span the null space + # of R. Valid even when R is rank-deficient: those vectors + # still have zero singular value, so they lie in the (then + # larger) null space. + _, _, Vh = np.linalg.svd(R, full_matrices=True) + N = Vh[m:].T + + x = rng.standard_normal(n_vars - m) + norm = float(np.linalg.norm(x)) + while norm < 1e-12: + x = rng.standard_normal(n_vars - m) + norm = float(np.linalg.norm(x)) + q = N @ (x / norm) + + if layout.col_rows[k].size: + impact = chol @ q + if not bool(np.all(impact[layout.col_rows[k]] * layout.col_signs[k] >= 0.0)): + # Abandon the WHOLE candidate and restart from column 1. + # Redrawing only this column would be a distribution bug: + # q_k's law is conditional on q_1..q_{k-1}, so retrying + # column k alone conditions the retained prefix on + # "produced a failure here", which is not the marginal of + # the accepted joint draw. Abandoning everything keeps the + # procedure plain rejection sampling over the whole Q. + return None + Q[:, k] = q + return Q + + def shock_coords(self, n_vars: int) -> list[str]: + """Named shocks in user order, then rotation-arbitrary padding.""" + return SignRestriction._build_shock_coords(list(self.shock_names), n_vars) diff --git a/tests/test_zero_sign_restriction.py b/tests/test_zero_sign_restriction.py new file mode 100644 index 00000000..229c55ca --- /dev/null +++ b/tests/test_zero_sign_restriction.py @@ -0,0 +1,523 @@ +"""Tests for ZeroSignRestriction (Arias, Rubio-Ramirez & Waggoner, 2018).""" + +import warnings + +import arviz as az +import numpy as np +import pytest +import xarray as xr +from pydantic import ValidationError + +from impulso.identification import SignRestriction, ZeroSignRestriction +from impulso.protocols import IdentificationScheme + +# -------------------------------------------------------------------------- +# Fixtures +# -------------------------------------------------------------------------- + + +@pytest.fixture +def synthetic_idata_3v(): + """Synthetic InferenceData mimicking a fitted 3-var VAR(1). + + Same recipe as the shared `synthetic_idata_2v` fixture, one variable + wider: 2 chains, 30 draws, 3 variables, 1 lag. + """ + rng = np.random.default_rng(7) + n_chains, n_draws, n_vars, n_lags = 2, 30, 3, 1 + + B = rng.standard_normal((n_chains, n_draws, n_vars, n_vars * n_lags)) * 0.2 + intercept = rng.standard_normal((n_chains, n_draws, n_vars)) * 0.01 + + sigma = np.zeros((n_chains, n_draws, n_vars, n_vars)) + L = np.zeros_like(sigma) + for c in range(n_chains): + for d in range(n_draws): + A = rng.standard_normal((n_vars, n_vars)) * 0.5 + sigma[c, d] = A @ A.T + np.eye(n_vars) + L[c, d] = np.linalg.cholesky(sigma[c, d]) + + posterior = xr.Dataset({ + "B": xr.DataArray(B, dims=["chain", "draw", "var", "coeff"]), + "intercept": xr.DataArray(intercept, dims=["chain", "draw", "var"]), + "Sigma": xr.DataArray( + sigma, + dims=["chain", "draw", "var1", "var2"], + coords={"var1": ["y1", "y2", "y3"], "var2": ["y1", "y2", "y3"]}, + ), + "L": xr.DataArray(L, dims=["chain", "draw", "var1", "var2"]), + }) + return az.InferenceData(posterior=posterior) + + +def _quiet(**kwargs) -> ZeroSignRestriction: + """Build a scheme, suppressing the unsigned-shock UserWarning. + + Most numerics tests deliberately leave column signs unpinned; the + warning is asserted on directly in `TestValidation`. + """ + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + return ZeroSignRestriction(**kwargs) + + +# -------------------------------------------------------------------------- +# Validation (n-independent, at construction) +# -------------------------------------------------------------------------- + + +class TestValidation: + def test_frozen(self): + scheme = _quiet(shock_names=["s1", "s2"], zero_restrictions={"y1": ["s2"]}) + with pytest.raises(ValidationError): + scheme.n_rotations = 5 + + def test_satisfies_protocol(self): + scheme = _quiet(shock_names=["s1", "s2"], zero_restrictions={"y1": ["s2"]}) + assert isinstance(scheme, IdentificationScheme) + + def test_samples_rotations_flag(self): + scheme = _quiet(shock_names=["s1", "s2"], zero_restrictions={"y1": ["s2"]}) + assert scheme._samples_rotations is True + + def test_empty_shock_names(self): + with pytest.raises(ValidationError, match="at least one structural shock"): + ZeroSignRestriction(shock_names=[], zero_restrictions={"y1": ["s2"]}) + + def test_duplicate_shock_names(self): + with pytest.raises(ValidationError, match="Duplicate shock names"): + ZeroSignRestriction(shock_names=["s1", "s1"], zero_restrictions={"y1": ["s1"]}) + + def test_reserved_prefix(self): + with pytest.raises(ValidationError, match="reserved prefix"): + ZeroSignRestriction(shock_names=["s1", "unidentified_1"], zero_restrictions={"y1": ["s1"]}) + + def test_unknown_shock_in_zero_restrictions(self): + with pytest.raises(ValidationError, match="unknown shock"): + ZeroSignRestriction(shock_names=["s1"], zero_restrictions={"y1": ["nope"]}) + + def test_unknown_shock_in_sign_restrictions(self): + with pytest.raises(ValidationError, match="unknown shock"): + ZeroSignRestriction(shock_names=["s1"], sign_restrictions={"y1": {"nope": "+"}}) + + def test_conflicting_zero_and_sign_cell(self): + with pytest.raises(ValidationError, match="conflict"): + ZeroSignRestriction( + shock_names=["s1", "s2"], + zero_restrictions={"y1": ["s2"]}, + sign_restrictions={"y1": {"s2": "+"}}, + ) + + def test_bad_sign_token(self): + with pytest.raises(ValidationError, match="non-sign token"): + ZeroSignRestriction(shock_names=["s1"], sign_restrictions={"y1": {"s1": "positive"}}) + + def test_both_dicts_empty(self): + with pytest.raises(ValidationError, match="at least one of zero_restrictions"): + ZeroSignRestriction(shock_names=["s1"]) + + def test_named_shock_without_sign_warns(self): + with pytest.warns(UserWarning, match="no sign restriction"): + ZeroSignRestriction(shock_names=["s1", "s2"], zero_restrictions={"y1": ["s2"]}) + + def test_fully_signed_shocks_do_not_warn(self): + with warnings.catch_warnings(): + warnings.simplefilter("error") + ZeroSignRestriction( + shock_names=["s1", "s2"], + zero_restrictions={"y1": ["s2"]}, + sign_restrictions={"y2": {"s1": "+", "s2": "+"}}, + ) + + +# -------------------------------------------------------------------------- +# identify() entry checks (need n_vars) +# -------------------------------------------------------------------------- + + +def _chol(n_vars: int, n_draws: int = 5, seed: int = 0) -> np.ndarray: + rng = np.random.default_rng(seed) + sigma = np.zeros((1, n_draws, n_vars, n_vars)) + for d in range(n_draws): + A = rng.standard_normal((n_vars, n_vars)) + sigma[0, d] = A @ A.T + np.eye(n_vars) + return np.linalg.cholesky(sigma) + + +class TestIdentifyEntryChecks: + def test_rank_condition_n2_two_zeros(self): + """n = 2, one shock with 2 zeros: z_1 = 2 > n - 1 = 1.""" + scheme = _quiet(shock_names=["s1", "s2"], zero_restrictions={"y1": ["s2"], "y2": ["s2"]}) + with pytest.raises(ValueError, match="at most"): + scheme.identify(_chol(2), ["y1", "y2"]) + + def test_rank_condition_n3_two_shocks_two_zeros(self): + """n = 3 with z = (2, 2): the second sorted position needs z <= 1.""" + scheme = _quiet( + shock_names=["s1", "s2", "s3"], + zero_restrictions={"y1": ["s2", "s3"], "y2": ["s2", "s3"]}, + ) + with pytest.raises(ValueError, match=r"n - j"): + scheme.identify(_chol(3), ["y1", "y2", "y3"]) + + def test_unknown_variable(self): + scheme = _quiet(shock_names=["s1", "s2"], zero_restrictions={"nope": ["s2"]}) + with pytest.raises(ValueError, match="unknown variable"): + scheme.identify(_chol(2), ["y1", "y2"]) + + def test_too_many_shocks(self): + scheme = _quiet(shock_names=["s1", "s2", "s3"], zero_restrictions={"y1": ["s2"]}) + with pytest.raises(ValueError, match="at most n_vars structural shocks"): + scheme.identify(_chol(2), ["y1", "y2"]) + + def test_horizon_without_posterior(self): + scheme = _quiet( + shock_names=["s1", "s2"], + zero_restrictions={"y1": ["s2"]}, + sign_restrictions={"y2": {"s1": "+", "s2": "+"}}, + restriction_horizon=2, + ) + with pytest.raises(ValueError, match="requires the full posterior"): + scheme.identify(_chol(2), ["y1", "y2"]) + + +# -------------------------------------------------------------------------- +# Numerics — the Cholesky-degeneracy anchor is the go/no-go gate +# -------------------------------------------------------------------------- + + +class TestCholeskyDegeneracyAnchor: + """Full triangular zeros collapse every null space to one dimension. + + With z_j = n - j for every shock the construction has no freedom left: + each column is determined up to its sign, P is lower triangular, and + P P' = Sigma. The only lower-triangular square root of Sigma is the + Cholesky factor up to column signs, so |P| must reproduce + |cholesky(Sigma)| exactly and acceptance must be 1.0. + """ + + def test_matches_cholesky_in_absolute_value(self, synthetic_idata_3v): + L = synthetic_idata_3v.posterior["L"].values + sigma = synthetic_idata_3v.posterior["Sigma"].values + scheme = _quiet( + shock_names=["s1", "s2", "s3"], + zero_restrictions={"y1": ["s2", "s3"], "y2": ["s3"]}, + n_rotations=1, + random_seed=0, + ) + P = scheme.identify(L, ["y1", "y2", "y3"]) + + assert scheme._last_diagnostics["zero_sign_acceptance_rate"] == 1.0 + expected = np.abs(np.linalg.cholesky(sigma)) + max_dev = float(np.max(np.abs(np.abs(P) - expected))) + assert max_dev < 1e-8, f"max deviation from |cholesky(Sigma)| = {max_dev:g}" + + def test_diagonal_signs_pin_the_cholesky_factor_exactly(self, synthetic_idata_3v): + L = synthetic_idata_3v.posterior["L"].values + sigma = synthetic_idata_3v.posterior["Sigma"].values + scheme = ZeroSignRestriction( + shock_names=["s1", "s2", "s3"], + zero_restrictions={"y1": ["s2", "s3"], "y2": ["s3"]}, + sign_restrictions={"y1": {"s1": "+"}, "y2": {"s2": "+"}, "y3": {"s3": "+"}}, + # Each column's sign is an independent coin flip in the degenerate + # (1-dimensional null space) case, so a candidate satisfies all + # three diagonal pins with probability 1/8. The budget has to + # cover that: P(no acceptance in 300 tries) is around 1e-17. + n_rotations=300, + random_seed=0, + ) + P = scheme.identify(L, ["y1", "y2", "y3"]) + + assert scheme._last_diagnostics["zero_sign_acceptance_rate"] == 1.0 + np.testing.assert_allclose(P, np.linalg.cholesky(sigma), atol=1e-8) + + +class TestClosedFormAndInvariants: + def test_n2_analytic_column(self, synthetic_idata_2v): + """One zero on a 2-var VAR identifies that column up to sign.""" + L = synthetic_idata_2v.posterior["L"].values + sigma = synthetic_idata_2v.posterior["Sigma"].values + scheme = _quiet( + shock_names=["s1", "s2"], + zero_restrictions={"y1": ["s2"]}, + n_rotations=1, + random_seed=3, + ) + P = scheme.identify(L, ["y1", "y2"]) + + s11 = sigma[..., 0, 0] + s12 = sigma[..., 0, 1] + s22 = sigma[..., 1, 1] + expected_lower = np.sqrt(s22 - s12**2 / s11) + np.testing.assert_allclose(np.abs(P[..., 0, 1]), 0.0, atol=1e-8) + np.testing.assert_allclose(np.abs(P[..., 1, 1]), expected_lower, atol=1e-8) + + def test_zero_cells_are_exact(self, synthetic_idata_3v): + """A generic single zero restriction is satisfied to SVD precision.""" + L = synthetic_idata_3v.posterior["L"].values + scheme = _quiet( + shock_names=["s1", "s2"], + zero_restrictions={"y2": ["s1"]}, + n_rotations=50, + random_seed=11, + ) + P = scheme.identify(L, ["y1", "y2", "y3"]) + assert np.nanmax(np.abs(P[..., 1, 0])) < 1e-10 + assert scheme._last_diagnostics["zero_sign_max_zero_violation"] < 1e-10 + + def test_reproduces_sigma(self, synthetic_idata_3v): + """P P' = Sigma on every accepted draw — Q is exactly orthogonal.""" + L = synthetic_idata_3v.posterior["L"].values + sigma = synthetic_idata_3v.posterior["Sigma"].values + scheme = _quiet( + shock_names=["s1", "s2"], + zero_restrictions={"y2": ["s1"]}, + sign_restrictions={"y1": {"s1": "+"}, "y3": {"s2": "-"}}, + n_rotations=200, + random_seed=5, + ) + P = scheme.identify(L, ["y1", "y2", "y3"]) + finite = np.isfinite(P).all(axis=(-2, -1)) + assert finite.any() + np.testing.assert_allclose(P[finite] @ np.swapaxes(P[finite], -1, -2), sigma[finite], atol=1e-10) + + def test_columns_land_in_user_order(self, synthetic_idata_3v): + """Internal sorting is invisible: zeros appear under the named shock.""" + L = synthetic_idata_3v.posterior["L"].values + # s3 is the most-restricted shock, so construction reorders — but the + # zero cells must still be at (y1, s3), (y2, s3) and (y1, s2). + scheme = _quiet( + shock_names=["s1", "s2", "s3"], + zero_restrictions={"y1": ["s2", "s3"], "y2": ["s3"]}, + n_rotations=1, + random_seed=0, + ) + P = scheme.identify(L, ["y1", "y2", "y3"]) + assert np.max(np.abs(P[..., 0, 2])) < 1e-10 + assert np.max(np.abs(P[..., 1, 2])) < 1e-10 + assert np.max(np.abs(P[..., 0, 1])) < 1e-10 + assert np.min(np.abs(P[..., 0, 0])) > 1e-8 + + +# -------------------------------------------------------------------------- +# Horizon sign restrictions +# -------------------------------------------------------------------------- + + +class TestHorizonSigns: + def test_signs_hold_at_every_horizon(self, synthetic_idata_3v): + from impulso._linalg import lag_matrices + from impulso._ma import compute_ma_phi + + horizon = 3 + L = synthetic_idata_3v.posterior["L"].values + B = synthetic_idata_3v.posterior["B"].values + scheme = _quiet( + shock_names=["s1", "s2"], + zero_restrictions={"y2": ["s1"]}, + sign_restrictions={"y1": {"s1": "+"}, "y3": {"s2": "-"}}, + restriction_horizon=horizon, + n_rotations=500, + random_seed=17, + ) + # Random synthetic B makes the horizon signs hard, so a fair share of + # draws end as NaN. That is the documented policy, not a failure — + # the loop below only inspects the accepted draws. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + P = scheme.identify(L, ["y1", "y2", "y3"], posterior=synthetic_idata_3v.posterior) + + n_checked = 0 + for c in range(P.shape[0]): + for d in range(P.shape[1]): + if not np.isfinite(P[c, d]).all(): + continue + n_checked += 1 + Phi = compute_ma_phi(lag_matrices(B[c, d], 1), horizon) + for h in range(horizon + 1): + irf = Phi[h] @ P[c, d] + assert irf[0, 0] >= -1e-12 + assert irf[2, 1] <= 1e-12 + assert n_checked > 0 + + def test_zero_restriction_is_impact_only(self, synthetic_idata_3v): + """Zeros bind at h = 0; later horizons are free.""" + L = synthetic_idata_3v.posterior["L"].values + scheme = _quiet( + shock_names=["s1", "s2"], + zero_restrictions={"y2": ["s1"]}, + sign_restrictions={"y1": {"s1": "+"}}, + restriction_horizon=2, + n_rotations=500, + random_seed=23, + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + P = scheme.identify(L, ["y1", "y2", "y3"], posterior=synthetic_idata_3v.posterior) + assert np.nanmax(np.abs(P[..., 1, 0])) < 1e-10 + + +# -------------------------------------------------------------------------- +# Distributional equivalence with SignRestriction when there are no zeros +# -------------------------------------------------------------------------- + + +class TestSignOnlyEquivalence: + def test_matches_sign_restriction_without_zeros(self, synthetic_idata_2v): + """No zeros => the recursion is Gram-Schmidt on Gaussians, i.e. Haar. + + `SignRestriction` draws from SO(n) via `special_ortho_group`; this + scheme draws from O(n). With one named shock the second column is + unidentified, so the two accepted-P laws coincide (a reflection only + flips the unpinned column). Tolerances are coarse on purpose: this + is a two-sample Monte Carlo comparison on 2 x 50 draws, so the + sampling error on an acceptance rate near 0.5 is order 0.05 and the + means of the restricted cells are noisier still. + """ + L = synthetic_idata_2v.posterior["L"].values + zs = ZeroSignRestriction( + shock_names=["s"], + sign_restrictions={"y1": {"s": "+"}, "y2": {"s": "-"}}, + n_rotations=200, + random_seed=101, + ) + sr = SignRestriction( + restrictions={"y1": {"s": "+"}, "y2": {"s": "-"}}, + n_rotations=200, + random_seed=101, + ) + P_zs = zs.identify(L, ["y1", "y2"]) + P_sr = sr.identify(L, ["y1", "y2"]) + + rate_zs = zs._last_diagnostics["zero_sign_acceptance_rate"] + rate_sr = sr._last_acceptance_rate + assert abs(rate_zs - rate_sr) < 0.15 + + for row in (0, 1): + m_zs = float(np.nanmean(P_zs[..., row, 0])) + m_sr = float(np.nanmean(P_sr[..., row, 0])) + assert abs(m_zs - m_sr) < 0.15, f"row {row}: {m_zs:g} vs {m_sr:g}" + + +# -------------------------------------------------------------------------- +# Failure policy and diagnostics +# -------------------------------------------------------------------------- + + +class TestFailurePolicy: + def _impossible(self, **kwargs) -> ZeroSignRestriction: + # y1 and y2 must both rise on impact from s1 while s1 is also + # restricted to zero on y3 — feasible in principle, but with + # n_rotations=1 the single candidate almost never complies. + return ZeroSignRestriction( + shock_names=["s1", "s2"], + zero_restrictions={"y3": ["s1"]}, + sign_restrictions={"y1": {"s1": "+", "s2": "+"}, "y2": {"s1": "+", "s2": "-"}}, + n_rotations=1, + random_seed=1, + **kwargs, + ) + + def test_nan_and_summary_warning(self, synthetic_idata_3v): + L = synthetic_idata_3v.posterior["L"].values + scheme = self._impossible() + with pytest.warns(UserWarning, match="NaN"): + P = scheme.identify(L, ["y1", "y2", "y3"]) + assert np.isnan(P).any() + assert scheme._last_diagnostics["zero_sign_failed_draws"] > 0 + assert 0.0 <= scheme._last_diagnostics["zero_sign_acceptance_rate"] < 1.0 + + def test_raise_policy(self, synthetic_idata_3v): + L = synthetic_idata_3v.posterior["L"].values + scheme = self._impossible(on_failure="raise") + with pytest.raises(ValueError, match="No admissible rotation"): + scheme.identify(L, ["y1", "y2", "y3"]) + + def test_no_fallback_to_cholesky(self, synthetic_idata_3v): + """Failed draws are NaN, never silently the unrotated factor.""" + L = synthetic_idata_3v.posterior["L"].values + scheme = self._impossible() + with pytest.warns(UserWarning): + P = scheme.identify(L, ["y1", "y2", "y3"]) + failed = ~np.isfinite(P).all(axis=(-2, -1)) + assert failed.any() + assert np.isnan(P[failed]).all() + + def test_diagnostics_keys(self, synthetic_idata_3v): + L = synthetic_idata_3v.posterior["L"].values + scheme = _quiet( + shock_names=["s1", "s2"], + zero_restrictions={"y2": ["s1"]}, + n_rotations=10, + random_seed=2, + ) + scheme.identify(L, ["y1", "y2", "y3"]) + assert set(scheme._last_diagnostics) == { + "zero_sign_acceptance_rate", + "zero_sign_failed_draws", + "zero_sign_failed_fraction", + "zero_sign_mean_attempts", + "zero_sign_max_zero_violation", + } + assert scheme._last_diagnostics["zero_sign_mean_attempts"] >= 1.0 + + def test_does_not_set_sign_restriction_acceptance_rate(self, synthetic_idata_3v): + """The scheme must not attach the misnamed SignRestriction attr.""" + L = synthetic_idata_3v.posterior["L"].values + scheme = _quiet(shock_names=["s1", "s2"], zero_restrictions={"y2": ["s1"]}, n_rotations=10) + scheme.identify(L, ["y1", "y2", "y3"]) + assert getattr(scheme, "_last_acceptance_rate", None) is None + + +# -------------------------------------------------------------------------- +# Pipeline integration +# -------------------------------------------------------------------------- + + +class TestPipeline: + def test_importable_from_impulso(self): + import impulso + + assert impulso.ZeroSignRestriction is ZeroSignRestriction + assert "ZeroSignRestriction" in impulso.__all__ + + def test_shock_coords_padding(self): + scheme = _quiet(shock_names=["s1"], zero_restrictions={"y2": ["s1"]}) + assert scheme.shock_coords(3) == ["s1", "unidentified_1", "unidentified_2"] + + def test_diagnostics_reach_shock_matrix_attrs(self, synthetic_idata_3v): + import pandas as pd + + from impulso.data import VARData + from impulso.identified import IdentifiedVAR + from impulso.volatility import Constant + + rng = np.random.default_rng(0) + index = pd.date_range("2000-01-01", periods=60, freq="QS") + data = VARData(endog=rng.standard_normal((60, 3)), endog_names=["y1", "y2", "y3"], index=index) + scheme = _quiet( + shock_names=["s1"], + zero_restrictions={"y2": ["s1"]}, + sign_restrictions={"y1": {"s1": "+"}}, + n_rotations=100, + random_seed=4, + ) + identified = IdentifiedVAR( + idata=synthetic_idata_3v, + n_lags=1, + data=data, + var_names=["y1", "y2", "y3"], + volatility=Constant(), + scheme=scheme, + ) + sm = identified.shock_matrix() + assert "zero_sign_acceptance_rate" in sm.attrs + assert "sign_restriction_acceptance_rate" not in sm.attrs + assert list(sm.coords["shock"].values) == ["s1", "unidentified_1", "unidentified_2"] + + with pytest.warns(UserWarning, match="rotation-arbitrary"): + fevd = identified.fevd(horizon=4) + da = fevd.idata.posterior_predictive["fevd"] + assert np.isnan(da.sel(shock="unidentified_1").values).all() + assert np.isnan(da.sel(shock="unidentified_2").values).all() + assert not np.isnan(da.sel(shock="s1").values).any()