diff --git a/CONTEXT.md b/CONTEXT.md index f62ffad1..a68096a2 100644 --- a/CONTEXT.md +++ b/CONTEXT.md @@ -74,6 +74,22 @@ _Avoid_: "scenario" for an all-shocks-adjust conditional forecast — scenarios A conditional forecast with structural attribution (Antolín-Díaz, Petrella & Rubio-Ramírez 2021): pinned variable paths must be absorbed by a named `adjusting` set of shocks (non-adjusting shocks keep their unconditional draws), and/or future shock paths are prescribed directly. Computed by `IdentifiedVAR.structural_scenario()`. _Avoid_: "conditional forecast" when an adjusting set is named — the restriction is the point. +**Entropic tilting**: +Imposing a distributional statement on a forecast by *reweighting the draws it already produced* — never re-solving, never re-simulating. The weights are the ones closest to uniform in relative entropy subject to the requested moments (Robertson, Tallman & Whiteman 2005), so the tilt moves the forecast as little as the target allows. Entry points are `ForecastResult.tilt()` and `ConditionalForecastResult.tilt()` (inherited by `ScenarioResult`), returning a `TiltedForecastResult` that holds the parent draws by reference. Hard and soft conditioning *chain* (`conditional_forecast(...).tilt(...)`) rather than mixing in one call: hard pins hold pathwise on every draw and reweighting never moves a draw, so preservation is a theorem, not a code path. +_Avoid_: "importance weights" / "importance sampling" — there is no proposal distribution here; the draws come from the posterior predictive itself. "Reweighting" and "tilting weights" are the terms. + +**Probability target (`ProbabilityTarget`), moment target (`MomentTarget`)**: +The soft counterparts of the condition vocabulary — statements about the forecast *distribution* rather than about every draw. `ProbabilityTarget(variable, horizon, threshold, probability, direction)` requires an event to carry a given probability after tilting (`probability=1.0` is exact conditioning on it); `MomentTarget(variable, horizon, mean)` fixes a tilted mean. `horizon` is 1-based, matching `VariablePath`. A target no draw satisfies is refused with the draw counts — reweighting cannot move mass where there is none. +_Avoid_: "soft condition" as an API term — the objects are targets; "soft conditioning" stays prose. + +**Effective sample size (ESS)**: +The Kish quantity `1 / Σ wᵢ²` reported on every tilted result: how many draws the tilt is effectively using. Equals the draw count under uniform weights, `1` when all mass sits on one draw. `ess_fraction` is its share of the sample, and a tilt below 10% of the sample warns — tilted summaries are honest posterior quantities only to the extent that enough draws carry weight. +_Avoid_: confusing this with the MCMC effective sample size in `arviz.ess` — that measures chain autocorrelation, this measures weight concentration. + +**Reverse stress**: +Scenario analysis run backwards: name the outcome, get the shocks. `IdentifiedVAR.reverse_stress(variable, threshold, steps, ...)` draws an unconditional forecast together with the structural shocks behind it, tilts onto the stress event, and reports the **shock cocktail** — the tilted-weighted mean of the retained shocks, in one-standard-deviation units. Because it averages realised draws rather than solving a projection, the cocktail inherits the model's own shock correlations and needs no norm choice. Its magnitude `q = ‖E_w[ε]‖²` is in the same units as the plausibility statistic. +_Avoid_: describing the cocktail as "the cause" of the outcome — it is a conditional mean, an association under the estimated model; other configurations in the retained set may look nothing like it. + **Condition vocabulary (`ShockPath`, `VariablePath`)**: Frozen spec objects expressing scenario content. `ShockPath(shock, values, start, end)` sets a structural shock's path — values in one-standard-deviation units, scalar broadcast (`0.0` = "switch the shock off") or explicit array; `start`/`end` timestamps window it in-sample (counterfactual), while on the forecast axis values run from step 1 with `NaN` marking free entries. `VariablePath(variable, values)` pins a future endogenous path the same NaN-masked way. A scalar stays scalar until *application* time, where it broadcasts to the full resolved window — deliberately not the same as a length-1 array, which pins exactly one period. Each method accepts only the condition types legal for it — illegal combinations are unrepresentable rather than validated away. _Avoid_: "Conditions object" / "Scenario object" for these primitives — a bundling `Scenario` container may arrive later for connector round-trips; the primitives are paths. @@ -84,7 +100,7 @@ _Avoid_: "driving" / "offsetting" shocks in API surface (fine in prose, where AD **Plausibility statistic (q)**: Per-draw squared Mahalanobis distance of a scenario's binding restriction values from their unconditional law, `q = c̄′(C C′)⁻¹ c̄ = ‖μ*‖²`, plus `‖v_S‖²` for prescribed shock paths — distributed `χ²_r` under the model when all shocks adjust (`r` = number of binding restrictions), reported with `r` and the tail probability `P(χ²_r ≥ q)`. The Leeper–Zha "modest interventions" check in the ADPRR lineage: large `q` (tiny tail probability) means the scenario demands incredible shocks and the model's answer should not be trusted. Stored per draw as a `plausibility` variable on `ConditionalForecastResult` and `ScenarioResult`, alongside the ADPRR-calibrated companion `q_cal ∈ [0.5, 1]` (`plausibility_calibrated`; McCulloch binomial matching with `z = q/2`) — finite only under the unconditional-variance mode, pegged at its ceiling of 1 under hard pins, floored at 0.5 with no conditions. -_Avoid_: calling it a Kullback–Leibler divergence under hard conditions — the conditional law is singular there and that KL is infinite; the KL form applies only to future *soft* conditioning. "Modesty statistic" stays prose-only. +_Avoid_: calling it a Kullback–Leibler divergence under hard conditions — the conditional law is singular there and that KL is infinite. The KL form belongs to soft conditioning, where it is real: a tilted result reports `kl_divergence` alongside its ESS, and `reverse_stress` calibrates *that* divergence into its `q_cal`. Name which one you mean. "Modesty statistic" stays prose-only. **at**: The time-index parameter on time-varying queries (`impulse_response(at=...)`, `fevd(at=...)`). Accepts an integer `t`, the literal `"last"` (most recent), `"all"` (full T-axis returned in the result), or `None` (default; resolves to `"last"` for stochastic volatility, ignored for constant volatility). @@ -134,6 +150,8 @@ _Avoid_: "number of cointegrating vectors" in API surface (fine in prose); "coin - A **FittedVAR** computes **conditional forecasts** on its own — all shocks adjust, so no identification scheme is involved (the dynamic-multiplier placement logic). - An **IdentifiedVAR** computes **historical counterfactuals** and **structural scenarios** through the four-layer scenario engine (back out → constrain → solve → propagate); the propagate layer is shared with `forecast()` and the **historical decomposition**. - The **condition vocabulary** is consumed by all three scenario methods; each method accepts only the condition types legal for it. +- A **forecast result** (plain, conditional, or scenario) is **entropically tilted** onto **targets** by `.tilt()`, producing a `TiltedForecastResult`. Tilting is a layer *on top of* the scenario engine, not a stage inside it: it consumes draws and returns weights. +- An **IdentifiedVAR** runs **reverse stress** by drawing forecast paths together with their structural shocks and tilting onto the stress event; the **shock cocktail** is the weighted mean of the shocks that survive. - 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. @@ -161,6 +179,11 @@ _Avoid_: "number of cointegrating vectors" in API surface (fine in prose); "coin > > **User:** "2020Q2 is wrecking my estimates. Can I stop dummying it out?" > **Library:** `VAR(lags=4, error_dist="student_t").fit(VARData(...))`. The t likelihood downweights the observation automatically; the degrees of freedom come back in the posterior as `nu`. Pass `StudentT(nu=5.0)` to fix them instead — the robust choice on short samples. +> **User:** "I don't have a path — I just want a 40% chance of recession next year." +> **Library:** `fitted.forecast(steps=8).tilt([ProbabilityTarget(variable="gdp", horizon=4, threshold=0.0, probability=0.4)])`. The draws are reweighted, not re-solved; check `result.summary()["ess_fraction"]` before reading the tilted bands. +> +> **User:** "What shocks would push inflation below target by 2027?" +> **Library:** `identified.reverse_stress(variable="inflation", threshold=1.0, steps=12, horizon=8)`. `result.shock_cocktail()` is the average structural configuration among the draws that got there. ## Conventions diff --git a/docs/adr/0009-entropic-tilting-post-hoc-reweighting.md b/docs/adr/0009-entropic-tilting-post-hoc-reweighting.md new file mode 100644 index 00000000..7c962ca6 --- /dev/null +++ b/docs/adr/0009-entropic-tilting-post-hoc-reweighting.md @@ -0,0 +1,25 @@ +# Soft conditioning is a post-hoc reweighting layer, not a second solver + +Distributional statements about a forecast — "give recession odds of 40%", "assume the market's mean path" — are imposed by *entropic tilting*: reweight the draws an existing forecast already produced so the requested moments hold, choosing the weights that minimise relative entropy to the untilted forecast (Robertson, Tallman & Whiteman 2005, *Journal of Money, Credit and Banking*). The draws never move. `ForecastResult.tilt()` and `ConditionalForecastResult.tilt()` return a `TiltedForecastResult` carrying the parent draws by reference plus a weight per draw, and every summary on it is weighted. + +Hard and soft conditioning are **chained, never mixed in one call**: `fitted.conditional_forecast(...).tilt(targets)`. This is not a limitation dodged by API design — it is exact. Hard pins hold pathwise on *every* draw, and reweighting a set of draws that all satisfy a constraint leaves them all satisfying it. Preservation is a theorem, so there is no interaction term to solve for. + +The solver splits on structure. A single `ProbabilityTarget` is a two-mass problem whose solution is `p / N_A` inside the event and `(1 - p) / (N - N_A)` outside, with `p = 1` collapsing to exact conditioning; no optimiser runs. Anything else goes through the convex dual `lambda* = argmin log((1/N) sum_i exp(lambda'(g_i - t)))`, minimised with the analytic gradient on a log-sum-exp-stabilised objective. + +`IdentifiedVAR.reverse_stress()` is the structural application: draw an unconditional forecast together with the structural shocks behind it, condition on the stress event, and report the **shock cocktail** — the tilted-weighted mean of the retained shocks. + +## Considered options + +- **Re-solve the forecast as importance sampling from a tilted proposal** — rejected: it needs a proposal distribution and a re-simulation pass for every target, and its weights carry proposal noise on top of the posterior's. Tilting reweights the *existing* posterior-predictive sample, so the untilted forecast and every tilt of it share the same draws exactly and can be compared line for line. +- **Kalman-smoother soft conditioning inside the scenario engine** — rejected for now, on the same grounds ADR-0005 deferred the smoother backend: it is the more general object (it would handle soft conditions on latent states, not just on forecast functionals) but it is a second solver to build and validate, and it cannot express a target on an arbitrary functional of the path. The tilting layer is distribution-free and composes with whatever the engine produces. +- **Cocktail as a minimum-norm projection onto the constraint set** — rejected: it needs an arbitrary norm choice, and the minimum-norm answer ignores the model's own shock correlations. The tilted conditional mean of *realised* draws inherits those correlations for free and is exactly the plain sample mean over the retained draws when `probability = 1.0`. +- **A `soft=` flag on `conditional_forecast`** — rejected: it invites the mixed hard/soft call that chaining already answers, and it would put a reweighting concept inside a solve. + +## Consequences + +- The relative entropy `KL = sum_i w_i log(N w_i)` is now a real, finite number. ADR-0005 reserved the Kullback–Leibler form for future soft conditioning precisely because the hard-conditioned shock law is singular; this is that future. The plausibility vocabulary therefore carries two divergences that must not be conflated — `q` (Mahalanobis, hard conditions) and `KL` (entropic, soft targets). +- Diagnostics are mandatory, not optional. A tilt that concentrates its weight on a handful of draws produces summaries that *look* like posterior quantities but rest on almost no sample, so the Kish effective sample size `1 / sum_i w_i^2` is reported on every result and a tilt below 10% of the draw count warns. +- Targets can be infeasible in a way solves never are: no reweighting can give positive probability to an event no draw satisfies. Empty-support and out-of-hull failures are raised at moment-construction time with the draw counts, before any optimiser sees them, and joint infeasibility is caught post-solve by comparing achieved against requested. +- `reverse_stress`'s `q_cal` applies the ADPRR binomial calibration to the entropic divergence, `q_cal = (1 + sqrt(1 - exp(-2·KL/d)))/2`. This is an **extension** of that calibration, not a result from the paper: it substitutes the now-finite entropic divergence for the `z = q/2` that hard conditioning makes infinite. It is documented as such on the method. +- Getting the shocks out of a forecast required a new engine entry point, `structural_forecast_draws`, which reproduces `structural_scenario_engine`'s no-ingredient branch exactly — same deterministic path, same forecast factors, same RNG stream — so matched seeds nest across the scenario family. Its cost is one forecast-sized array of shocks retained alongside the paths. +- Tilting is distribution-free: it consumes draws and returns weights, so it composes unchanged with heavier-tailed forecast errors when those arrive. Nothing in this layer assumes Gaussian innovations. diff --git a/docs/how-to/index.md b/docs/how-to/index.md index 7ff315b0..309b855f 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 +probabilistic-conditioning ``` diff --git a/docs/how-to/probabilistic-conditioning.md b/docs/how-to/probabilistic-conditioning.md new file mode 100644 index 00000000..4f4feb33 --- /dev/null +++ b/docs/how-to/probabilistic-conditioning.md @@ -0,0 +1,113 @@ +# Conditioning on Probabilities + +Some scenarios are not paths. "Assume a 40% chance of recession next year", "assume the market's mean rate path", "what shocks would push inflation below target?" — none of these pin a variable to a number on every draw. They are statements about the forecast *distribution*, and Impulso imposes them by **entropic tilting** {cite:p}`robertsonTallmanWhiteman2005,kruegerClarkRavazzolo2017`: reweight the draws you already have so the statement holds, using the weights closest (in relative entropy) to leaving them alone. + +## State a probability target + +Tilt a density forecast onto the probability you want: + +```python +from impulso.scenario import ProbabilityTarget + +forecast = fitted.forecast(steps=8, seed=0) + +tilted = forecast.tilt([ + ProbabilityTarget(variable="gdp", horizon=4, threshold=0.0, probability=0.4) +]) + +tilted.median() # weighted median +tilted.hdi() # weighted HDI +tilted.plot() # tilted fan against the untilted median +``` + +`horizon` is 1-based, the same convention as `VariablePath`. The event is `gdp < 0.0` at step 4; pass `direction="above"` for the other side. Use `MomentTarget(variable, horizon, mean)` to fix a mean instead of a probability. + +## Read the diagnostics before the forecast + +Tilting cannot create draws. If you ask for a probability the sample barely supports, the weights pile onto a handful of draws and every summary silently becomes a summary of those few draws. `summary()` reports what happened: + +```python +tilted.summary() +# {'ess': 612.4, 'ess_fraction': 0.153, 'kl_divergence': 0.44, 'n_draws': 4000, +# 'targets': [{'target': 'P(gdp[h=4] < 0)', 'requested': 0.4, +# 'achieved': 0.4, 'draws_in_event': 380}]} +``` + +- **`ess`** is the Kish effective sample size, `1 / Σ wᵢ²` — the number of draws the tilt is effectively using. `ess_fraction` expresses it as a share of `n_draws`, and a tilt below 10% of the sample warns. +- **`kl_divergence`** is how far the tilt moved the forecast. Zero means the target was already true; large values mean the model disagrees with you. +- **`achieved`** against **`requested`** confirms the solver hit the target. `draws_in_event` is how many of the original draws satisfied it before reweighting. + +A target no draw satisfies is refused outright, naming the counts: + +```python +forecast.tilt([ProbabilityTarget(variable="gdp", horizon=4, threshold=-50.0, probability=0.4)]) +# ValueError: 0 of 4000 draws satisfy P(gdp[h=4] < -50); the target is +# unachievable by reweighting — widen the threshold or increase the number of draws. +``` + +## Chain hard pins with soft targets + +Hard conditioning and probability targets combine by chaining, not by a mixed call: + +```python +from impulso.scenario import VariablePath + +conditional = fitted.conditional_forecast( + steps=8, + conditions=[VariablePath(variable="rate", values=rate_path)], + seed=0, +) +both = conditional.tilt([ + ProbabilityTarget(variable="gdp", horizon=4, threshold=0.0, probability=0.4) +]) +``` + +The order does not matter to the result and the pins are safe: they hold pathwise on *every* draw, and reweighting never moves a draw, so no weighting can break them. `ScenarioResult` inherits `tilt()` too, so a structural scenario can carry a probability target the same way. + +:::{admonition} Tilting needs a density forecast +:class: warning +`tilt()` refuses a mean forecast (`include_shock_uncertainty=False`). Mean-mode draws carry parameter uncertainty only, so probabilities read off them are not predictive probabilities. +::: + +## Run the scenario backwards + +Reverse stress testing asks the inverse question: not "what happens if this shock hits?" but "what shocks would deliver this outcome?" + +```python +result = identified.reverse_stress( + variable="inflation", + threshold=1.0, + steps=12, + horizon=8, + direction="below", + seed=0, +) + +result.shock_cocktail() # step x shock, in one-standard-deviation units +result.summary() +result.plot() +``` + +Impulso draws an unconditional forecast together with the structural shocks behind it, conditions on the event (`probability=1.0` by default — exact conditioning), and averages the shocks of the draws that produced the outcome. That average is the **shock cocktail**: the structural configuration most associated with the stress event, in one-standard-deviation units. + +Reading the output: + +- `baseline_probability` is how likely the event was before conditioning. If it is already 0.4, the "stress" is not stressful. +- The cocktail's entries are shock sizes. A `-2.1` on a demand shock at step 3 means the outcome is associated with a 2.1-standard-deviation negative demand shock there. +- `q = ‖cocktail‖²` is the cocktail's total magnitude in the same units as the scenario plausibility statistic, so `q = 9` reads as "a 3-standard-deviation configuration". +- `q_cal` is a separate reading: it applies the ADPRR binomial calibration to the *tilt's relative entropy* — how far conditioning on the event moved the forecast — mapping it onto `[0.5, 1]`, with 0.5 meaning the event cost nothing to impose. + +Softening the conditioning keeps more of the sample: + +```python +result = identified.reverse_stress( + variable="inflation", threshold=1.0, steps=12, probability=0.6, seed=0 +) +``` + +At `probability=0.6` the event carries 60% of the weight instead of all of it, the effective sample size rises, and the cocktail shrinks toward zero — the outcome is being made likely rather than certain. + +:::{admonition} Cocktails are associations, not causes +:class: note +The cocktail is a conditional mean over draws, so it reports which shock configurations *accompany* the outcome under the estimated model. It is not a unique cause: other configurations in the retained set may look nothing like the average, and a wide `hdi()` around the conditioned path is the signal to look at them. +::: diff --git a/docs/reference/plotting.md b/docs/reference/plotting.md index ebe61488..d6b8623f 100644 --- a/docs/reference/plotting.md +++ b/docs/reference/plotting.md @@ -10,6 +10,8 @@ plot_forecast plot_conditional_forecast plot_structural_scenario + plot_tilted_forecast + plot_reverse_stress plot_dynamic_multiplier plot_irf plot_fevd diff --git a/docs/reference/results.md b/docs/reference/results.md index f0946bbd..9107f59e 100644 --- a/docs/reference/results.md +++ b/docs/reference/results.md @@ -11,6 +11,8 @@ ForecastResult ConditionalForecastResult ScenarioResult + TiltedForecastResult + ReverseStressResult DynamicMultiplierResult IRFResult FEVDResult diff --git a/docs/reference/scenario.md b/docs/reference/scenario.md index 363c7273..96e71f22 100644 --- a/docs/reference/scenario.md +++ b/docs/reference/scenario.md @@ -1,10 +1,13 @@ # Scenario Conditions -The condition vocabulary for scenario analysis: `ShockPath` sets a -structural shock's path (historical counterfactuals today; prescribed -scenario shocks when `structural_scenario` arrives), `VariablePath` pins a -future endogenous path (for the forthcoming conditional-forecast and -structural-scenario methods). +The condition vocabulary for scenario analysis. `ShockPath` sets a +structural shock's path (in-sample counterfactual edits; forecast-side +prescriptions for `structural_scenario`) and `VariablePath` pins a future +endogenous path (`conditional_forecast`, `structural_scenario`) — both +hard conditions, holding on every draw. The *targets* are soft: they state +a fact about the forecast distribution that entropic tilting imposes by +reweighting draws (`ProbabilityTarget` for an event probability, +`MomentTarget` for a mean). ```{eval-rst} .. currentmodule:: impulso.scenario @@ -15,4 +18,6 @@ structural-scenario methods). ShockPath VariablePath + ProbabilityTarget + MomentTarget ``` diff --git a/docs/references.bib b/docs/references.bib index 2601e430..196c16d1 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -156,3 +156,23 @@ @article{blanchardQuah1989 number = {4}, pages = {655--673}, } + +@article{robertsonTallmanWhiteman2005, + author = {Robertson, John C. and Tallman, Ellis W. and Whiteman, Charles H.}, + title = {Forecasting Using Relative Entropy}, + journal = {Journal of Money, Credit and Banking}, + year = {2005}, + volume = {37}, + number = {3}, + pages = {383--401}, +} + +@article{kruegerClarkRavazzolo2017, + author = {Krueger, Fabian and Clark, Todd E. and Ravazzolo, Francesco}, + title = {Using Entropic Tilting to Combine {BVAR} Forecasts with External Nowcasts}, + journal = {Journal of Business \& Economic Statistics}, + year = {2017}, + volume = {35}, + number = {3}, + pages = {470--485}, +} diff --git a/src/impulso/__init__.py b/src/impulso/__init__.py index 8fe3f846..036747bc 100644 --- a/src/impulso/__init__.py +++ b/src/impulso/__init__.py @@ -33,13 +33,15 @@ IntegrationOrderResult, IRFResult, LagOrderResult, + ReverseStressResult, ScenarioResult, StationarityTestResult, SVForecastResult, + TiltedForecastResult, VolatilityResult, ) from impulso.samplers import NUTSSampler - from impulso.scenario import ShockPath, VariablePath + from impulso.scenario import MomentTarget, ProbabilityTarget, ShockPath, VariablePath from impulso.sv.data import SVData from impulso.sv.fitted import FittedSV from impulso.sv.priors import SVDefaultPrior @@ -72,10 +74,13 @@ "LongRunRestriction", "MinnesotaPrior", "ModelEvidence", + "MomentTarget", "NIWPrior", "NUTSSampler", "PandemicBreak", + "ProbabilityTarget", "ProxySVAR", + "ReverseStressResult", "SVData", "SVDefaultPrior", "SVForecastResult", @@ -85,6 +90,7 @@ "StationarityTestResult", "StochasticVolatility", "StudentT", + "TiltedForecastResult", "VARData", "VariablePath", "VolatilityProcess", @@ -121,9 +127,13 @@ "ConditionalForecastResult": "impulso.results", "CounterfactualResult": "impulso.results", "ScenarioResult": "impulso.results", + "TiltedForecastResult": "impulso.results", + "ReverseStressResult": "impulso.results", "DynamicMultiplierResult": "impulso.results", "ShockPath": "impulso.scenario", "VariablePath": "impulso.scenario", + "ProbabilityTarget": "impulso.scenario", + "MomentTarget": "impulso.scenario", "IRFResult": "impulso.results", "FEVDResult": "impulso.results", "HistoricalDecompositionResult": "impulso.results", diff --git a/src/impulso/_scenario.py b/src/impulso/_scenario.py index 130ee3b5..930d9652 100644 --- a/src/impulso/_scenario.py +++ b/src/impulso/_scenario.py @@ -808,3 +808,66 @@ def _calibrate_scenario_q( if path_uncertainty == "unconditional" and n_prescribed == 0: return (1.0 + np.sqrt(1.0 - np.exp(-q / d_total))) / 2.0 return np.ones_like(q) + + +# --- reverse stress (structural forecast draws; ADR-0009) --- + + +def structural_forecast_draws( + identified: IdentifiedVAR, + steps: int, + seed: int | np.random.Generator | None = None, + exog_future: np.ndarray | None = None, +) -> tuple[np.ndarray, np.ndarray]: + """Unconditional forecast draws *and* the structural shocks behind them. + + Reverse stress testing needs the shocks that generated each forecast + draw, and neither `forecast()` nor the scenario engines return them. + This reproduces `structural_scenario_engine`'s no-ingredient branch + (no conditions, no prescriptions, every shock adjusting) exactly — + same deterministic path, same scheme-identified forecast factors, and + the same RNG stream contract: one generator, `_forecast_shock_matrices` + first (which consumes randomness only under time-varying volatility), + then per-step `standard_normal((C, D, n))` draws in `forecast()`'s + order. Under a shared seed the paths match `structural_scenario(steps, + seed=...)` draw for draw. + + Memory note: the returned shocks are forecast-sized `(C, D, steps, n)` + — the same footprint as the paths themselves. + + Args: + identified: The identified VAR. + steps: Forecast horizon. + seed: RNG seed or Generator. + exog_future: Future exogenous values, validated by the caller. + + Returns: + Tuple `(paths, eps)`: forecast paths `(C, D, steps, n)` and the + structural shocks `(C, D, steps, n)` in one-standard-deviation + units, aligned with `identified.shock_names` on the last axis. + """ + from impulso._linalg import lag_matrices + from impulso._propagate import propagate + + posterior = identified.idata.posterior + B_draws = posterior["B"].values + n_lags = identified.n_lags + n_chains, n_draws, n_vars, _ = B_draws.shape + + intercept = posterior["intercept"].values + forcing = np.broadcast_to(intercept[:, :, np.newaxis, :], (n_chains, n_draws, steps, n_vars)).copy() + if exog_future is not None: + forcing += np.einsum("cdij,hj->cdhi", posterior["B_exog"].values, exog_future) + A_lags = lag_matrices(B_draws, n_lags) + b_path = propagate(A_lags, forcing, identified.data.endog[-n_lags:]) + + rng = seed if isinstance(seed, np.random.Generator) else np.random.default_rng(seed) + P_path = _forecast_shock_matrices(identified, steps, rng) + + eps = np.empty((n_chains, n_draws, steps, n_vars)) + for h in range(steps): # per-step draws, forecast()'s stream order + eps[:, :, h, :] = rng.standard_normal((n_chains, n_draws, n_vars)) + + u = np.einsum("cdhij,cdhj->cdhi", P_path, eps) + deviation = propagate(A_lags, u, np.zeros((n_lags, n_vars))) + return b_path + deviation, eps diff --git a/src/impulso/_tilting.py b/src/impulso/_tilting.py new file mode 100644 index 00000000..86d2563f --- /dev/null +++ b/src/impulso/_tilting.py @@ -0,0 +1,500 @@ +"""Entropic tilting — post-hoc reweighting of forecast draws (ADR-0009). + +Solves the minimum-relative-entropy problem of Robertson, Tallman and +Whiteman (2005), + + min_w sum_i w_i log(N w_i) + s.t. sum_i w_i g_k(y_i) = t_k (k = 1..K), sum_i w_i = 1, w >= 0, + +over the draws of a forecast that has already been produced. The draws +never move — only their weights — so anything every draw already +satisfies (a hard pin from `conditional_forecast`, say) survives any +reweighting untouched. + +Two solvers share one entry point. A single `ProbabilityTarget` reduces +to a two-mass problem with a closed form (`p / N_A` inside the event, +`(1 - p) / (N - N_A)` outside), so no optimiser runs. Anything else goes +through the convex dual: `lambda* = argmin log((1/N) sum_i +exp(lambda'(g_i - t)))`, minimised by BFGS with the analytic gradient +`E_lambda[g] - t` and a log-sum-exp-stabilised objective. + +This module is pure numpy/scipy: impulso types are imported lazily or +under `TYPE_CHECKING` only. +""" + +from __future__ import annotations + +import warnings +from typing import TYPE_CHECKING + +import numpy as np + +if TYPE_CHECKING: + from impulso.results import ConditionalForecastResult, ForecastResult, TiltedForecastResult + from impulso.scenario import MomentTarget, ProbabilityTarget + + Target = ProbabilityTarget | MomentTarget + +# Default warning threshold on ESS as a fraction of the draw count. +ESS_WARN_FRACTION = 0.1 +# Post-solve tolerance on |achieved - requested| and the dual's norm guard. +ACHIEVED_TOL = 1e-6 +LAMBDA_GUARD = 1e6 + + +# --- moment construction ------------------------------------------------- + + +def target_label(target: Target) -> str: + """Human-readable label for a target, used as the `target` coordinate.""" + from impulso.scenario import ProbabilityTarget + + if isinstance(target, ProbabilityTarget): + sign = "<" if target.direction == "below" else ">" + return f"P({target.variable}[h={target.horizon}] {sign} {target.threshold:g})" + return f"E[{target.variable}[h={target.horizon}]]" + + +def build_moments( + forecast: np.ndarray, + targets: list[Target], + var_names: list[str], + steps: int, +) -> tuple[np.ndarray, np.ndarray]: + """Build the moment matrix `G` and requested moment vector `t`. + + Column `k` holds the moment function of target `k` evaluated on every + draw: the event indicator for a `ProbabilityTarget`, the level itself + for a `MomentTarget`. Draws are flattened over `(chain, draw)` in C + order, so `G[i]` lines up with `forecast.reshape(N, steps, n)[i]`. + + Feasibility that does not depend on the solver is checked here: an + event no draw satisfies (or, for `p < 1`, one every draw satisfies) + and a requested mean outside the draws' range are unachievable by + reweighting, whatever the optimiser does. + + Args: + forecast: Forecast draws of shape `(C, D, steps, n)`. + targets: The targets to impose (must be non-empty). + var_names: Endogenous variable names. + steps: Forecast horizon. + + Returns: + Tuple `(G, t)` with `G` of shape `(N, K)` and `t` of shape `(K,)`. + + Raises: + TypeError: If a target is neither a `ProbabilityTarget` nor a + `MomentTarget`. + ValueError: On an empty target list, an unknown variable, a + horizon beyond the forecast, or an unachievable target. + """ + from impulso.scenario import MomentTarget, ProbabilityTarget + + if not targets: + raise ValueError("Tilting requires at least one target; an empty target list would leave the weights uniform.") + n_chains, n_draws, _, n_vars = forecast.shape + n_total = n_chains * n_draws + flat = forecast.reshape(n_total, steps, n_vars) + + G = np.empty((n_total, len(targets))) + t = np.empty(len(targets)) + for k, target in enumerate(targets): + if not isinstance(target, ProbabilityTarget | MomentTarget): + raise TypeError(f"tilt() accepts ProbabilityTarget or MomentTarget, got {type(target).__name__}") + i, h = _resolve_target(target, var_names, steps) + y = flat[:, h, i] + if isinstance(target, ProbabilityTarget): + event = (y < target.threshold) if target.direction == "below" else (y > target.threshold) + _check_event_support(target, int(event.sum()), n_total) + G[:, k] = event.astype(float) + t[k] = target.probability + else: + _check_moment_support(target, y, n_total) + G[:, k] = y + t[k] = target.mean + return G, t + + +def _resolve_target(target: Target, var_names: list[str], steps: int) -> tuple[int, int]: + """Resolve a target to `(variable_index, 0-based step index)`.""" + if target.variable not in var_names: + raise ValueError(f"Unknown variable {target.variable!r}; available variables: {var_names}") + if target.horizon > steps: + raise ValueError( + f"Target for {target.variable!r} refers to horizon {target.horizon} " + f"but the forecast has only {steps} steps." + ) + return var_names.index(target.variable), target.horizon - 1 + + +def _check_event_support(target: ProbabilityTarget, n_event: int, n_total: int) -> None: + """Reject probability targets no reweighting of these draws can hit.""" + label = target_label(target) + if n_event == 0: + raise ValueError( + f"0 of {n_total} draws satisfy {label}; the target is unachievable by reweighting " + "— widen the threshold or increase the number of draws." + ) + if n_event == n_total and target.probability < 1.0: + raise ValueError( + f"All {n_total} of {n_total} draws satisfy {label}, so the only achievable " + f"probability is 1.0, not {target.probability:g}; tighten the threshold or " + "increase the number of draws." + ) + + +def _check_moment_support(target: MomentTarget, y: np.ndarray, n_total: int) -> None: + """Reject mean targets outside the convex hull of the draws.""" + lo, hi = float(y.min()), float(y.max()) + if not lo < target.mean < hi: + raise ValueError( + f"{target_label(target)} = {target.mean:g} lies outside the range spanned by the " + f"{n_total} draws ([{lo:g}, {hi:g}]); reweighting cannot move mass where there is " + "none — relax the target or increase the number of draws." + ) + + +# --- solvers ------------------------------------------------------------- + + +def solve_tilt(G: np.ndarray, t: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Solve for the minimum-relative-entropy weights. + + A single binary moment column (one `ProbabilityTarget`) takes the + closed-form two-mass solution; everything else takes the convex dual. + + Args: + G: Moment matrix `(N, K)` from `build_moments`. + t: Requested moments `(K,)`. + + Returns: + Tuple `(weights, achieved)`: weights `(N,)` summing to 1 and the + achieved moments `G' w` `(K,)`. + + Raises: + ValueError: If the dual fails to reproduce the requested moments + (jointly infeasible targets). + """ + if G.shape[1] == 1 and _is_binary(G[:, 0]): + return _closed_form_two_mass(G[:, 0], float(t[0])) + return _dual_solve(G, t) + + +def _is_binary(column: np.ndarray) -> bool: + """Whether a moment column is a 0/1 indicator.""" + return bool(np.all((column == 0.0) | (column == 1.0))) + + +def _closed_form_two_mass(indicator: np.ndarray, probability: float) -> tuple[np.ndarray, np.ndarray]: + """Closed-form tilt for a single event probability. + + Relative entropy is minimised by keeping the weights uniform *within* + the event and within its complement, so the whole problem collapses + to splitting mass `p` over `N_A` draws and `1 - p` over the rest. + At `p = 1` this is exact conditioning on the event. + """ + in_event = indicator > 0.5 + n_total = indicator.size + n_event = int(in_event.sum()) + weights = np.empty(n_total) + if probability >= 1.0: + weights[in_event] = 1.0 / n_event + weights[~in_event] = 0.0 + else: + weights[in_event] = probability / n_event + weights[~in_event] = (1.0 - probability) / (n_total - n_event) + return weights, np.array([float(weights[in_event].sum())]) + + +def _dual_solve(G: np.ndarray, t: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Minimise the convex dual with BFGS, then polish with Newton steps. + + BFGS on the log-sum-exp objective with the analytic gradient gets + close; its default gradient tolerance is looser than the moment + tolerance we then check against, so a handful of Newton steps on the + same objective (whose Hessian is the tilted covariance of `g`) take + the residual to machine precision. Both stages are on the same convex + dual, so the polish cannot move to a different optimum. + """ + from scipy.optimize import minimize + from scipy.special import logsumexp + + n_total, n_targets = G.shape + centred = G - t + + def objective(lam: np.ndarray) -> tuple[float, np.ndarray]: + z = centred @ lam + normaliser = logsumexp(z) + w = np.exp(z - normaliser) + return float(normaliser - np.log(n_total)), centred.T @ w + + result = minimize(objective, np.zeros(n_targets), jac=True, method="BFGS", options={"gtol": 1e-12}) + lam = _newton_polish(centred, np.asarray(result.x, dtype=float)) + weights = _weights_from_dual(centred, lam) + achieved = G.T @ weights + _check_dual_solution(achieved, t, lam) + return weights, achieved + + +def _newton_polish(centred: np.ndarray, lam: np.ndarray, max_steps: int = 50) -> np.ndarray: + """Newton iterations on the dual until the moment residual is machine-small.""" + for _ in range(max_steps): + w = _weights_from_dual(centred, lam) + gradient = centred.T @ w + if np.max(np.abs(gradient)) <= 1e-14: + break + centred_mean = centred - gradient + hessian = (centred_mean * w[:, np.newaxis]).T @ centred_mean + step, *_ = np.linalg.lstsq(hessian, gradient, rcond=None) + if not np.all(np.isfinite(step)): + break + lam = lam - step + return lam + + +def _weights_from_dual(centred: np.ndarray, lam: np.ndarray) -> np.ndarray: + """Normalised tilting weights `w_i ∝ exp(lambda'(g_i - t))`.""" + from scipy.special import logsumexp + + z = centred @ lam + weights = np.exp(z - logsumexp(z)) + return weights / weights.sum() + + +def _check_dual_solution(achieved: np.ndarray, t: np.ndarray, lam: np.ndarray) -> None: + """Reject a dual solution that misses the requested moments.""" + gap = float(np.max(np.abs(achieved - t))) + norm = float(np.linalg.norm(lam)) + if gap > ACHIEVED_TOL or norm > LAMBDA_GUARD: + raise ValueError( + "The requested targets are not jointly achievable by reweighting these draws: the " + f"entropic dual reached requested={np.array2string(t, precision=6)} vs " + f"achieved={np.array2string(achieved, precision=6)} (max gap {gap:.3g}, " + f"|lambda| = {norm:.3g}). Relax a target or increase the number of draws." + ) + + +# --- diagnostics --------------------------------------------------------- + + +def ess(weights: np.ndarray) -> float: + """Kish effective sample size `1 / sum_i w_i^2`. + + Equals `N` for uniform weights and `1` when all mass sits on one + draw, so `ess / N` reads directly as the fraction of the sample the + tilt actually uses. + """ + return float(1.0 / np.sum(weights**2)) + + +def kl_divergence(weights: np.ndarray) -> float: + """Relative entropy `sum_i w_i log(N w_i)` of the tilt from uniform. + + This is the objective the tilt minimises, and (unlike the hard-pin + case of ADR-0005) a genuinely finite divergence — soft conditioning + keeps the tilted law absolutely continuous with respect to the + untilted one. Zero-weight draws contribute nothing. + """ + n_total = weights.size + positive = weights[weights > 0.0] + return float(np.sum(positive * np.log(n_total * positive))) + + +# --- weighted summaries -------------------------------------------------- + + +def _sorted_columns(x: np.ndarray, weights: np.ndarray) -> tuple[np.ndarray, np.ndarray, tuple[int, ...]]: + """Sort each column of `x` and carry the weights along. + + Returns `(sorted values, normalised sorted weights, trailing shape)` + with both arrays shaped `(N, M)` where `M` is the product of the + trailing dimensions of `x`. + """ + x = np.asarray(x, dtype=float) + n_total = weights.size + if x.shape[0] != n_total: + raise ValueError(f"x must have {n_total} draws on its leading axis, got {x.shape[0]}") + trailing = x.shape[1:] + flat = x.reshape(n_total, -1) + order = np.argsort(flat, axis=0, kind="stable") + values = np.take_along_axis(flat, order, axis=0) + w = weights[order] + return values, w / w.sum(axis=0, keepdims=True), trailing + + +def weighted_quantile(x: np.ndarray, weights: np.ndarray, q: float) -> np.ndarray | float: + """Weighted quantile with the mid-cumulative-weight interpolation. + + Each sorted draw `j` is placed at `c_j = W_j - w_j / 2` (its + cumulative weight less half its own mass) and the quantile is read + off by linear interpolation between those knots, clamped at the + extremes. Under uniform weights this reproduces + `np.quantile(..., method="hazen")` exactly. + + Args: + x: Values with draws on the leading axis, shape `(N, ...)`. + weights: Non-negative weights `(N,)`; normalised internally. + q: Quantile level in `[0, 1]`. + + Returns: + Array shaped like `x.shape[1:]`, or a float when `x` is 1-D. + """ + values, w, trailing = _sorted_columns(x, weights) + knots = np.cumsum(w, axis=0) - w / 2.0 + out = np.array([np.interp(q, knots[:, m], values[:, m]) for m in range(values.shape[1])]) + return float(out[0]) if not trailing else out.reshape(trailing) + + +def weighted_hdi(x: np.ndarray, weights: np.ndarray, prob: float = 0.89) -> tuple[np.ndarray, np.ndarray]: + """Weighted highest-density interval by minimum width. + + Scans every sorted draw as a candidate left endpoint, takes the + first right endpoint carrying at least `prob` of the weight, and + keeps the narrowest such interval. Under uniform weights and a + non-integer `prob * N` this reproduces `arviz.hdi`. + + Args: + x: Values with draws on the leading axis, shape `(N, ...)`. + weights: Non-negative weights `(N,)`; normalised internally. + prob: Probability mass the interval must carry. + + Returns: + Tuple `(lower, upper)`, each shaped like `x.shape[1:]` (0-D + arrays when `x` is 1-D). + + Raises: + ValueError: If `prob` is not in `(0, 1]`. + """ + if not 0.0 < prob <= 1.0: + raise ValueError(f"prob must satisfy 0 < prob <= 1, got {prob}") + values, w, trailing = _sorted_columns(x, weights) + cumulative = np.cumsum(w, axis=0) + below = cumulative - w # weight strictly left of each candidate start + n_cols = values.shape[1] + lower = np.empty(n_cols) + upper = np.empty(n_cols) + for m in range(n_cols): + lower[m], upper[m] = _min_width_interval(values[:, m], cumulative[:, m], below[:, m], prob) + return lower.reshape(trailing), upper.reshape(trailing) + + +def _min_width_interval( + values: np.ndarray, + cumulative: np.ndarray, + below: np.ndarray, + prob: float, +) -> tuple[float, float]: + """Narrowest interval of sorted draws carrying at least `prob` mass.""" + # searchsorted on the cumulative weights gives, for each start i, the + # first index j with cumulative[j] >= below[i] + prob (the two-pointer + # sweep, vectorised — cumulative is non-decreasing). + ends = np.searchsorted(cumulative, below + prob - 1e-12, side="left") + feasible = ends < values.size + if not feasible.any(): + return float(values[0]), float(values[-1]) + starts = np.flatnonzero(feasible) + ends = ends[feasible] + widths = values[ends] - values[starts] + best = int(np.argmin(widths)) + return float(values[starts[best]]), float(values[ends[best]]) + + +# --- result assembly ----------------------------------------------------- + + +def tilt_result( + result: ForecastResult | ConditionalForecastResult, + targets: list[Target], + ess_warn_fraction: float = ESS_WARN_FRACTION, +) -> TiltedForecastResult: + """Tilt an existing forecast result onto the requested targets. + + The parent's `"forecast"` DataArray is carried into the new result by + reference — tilting adds weights, it does not copy or move draws. + + Args: + result: A density-mode `ForecastResult` or + `ConditionalForecastResult` (`ScenarioResult` included). + targets: `ProbabilityTarget` / `MomentTarget` list. + ess_warn_fraction: Warn when the effective sample size falls + below this fraction of the draw count. + + Returns: + A frozen `TiltedForecastResult`. + + Raises: + ValueError: If the parent result is a mean forecast, if + `ess_warn_fraction` is outside `[0, 1]`, or if the targets + are unachievable. + """ + import arviz as az + import xarray as xr + + from impulso.results import TiltedForecastResult + from impulso.scenario import ProbabilityTarget + + if result.mode != "density": + raise ValueError( + f"Entropic tilting needs a density forecast, but this result is a {result.mode!r} " + "forecast: mean-mode draws carry parameter uncertainty only, so probabilities read " + "off them are not predictive probabilities. Re-run with " + "include_shock_uncertainty=True." + ) + if not 0.0 <= ess_warn_fraction <= 1.0: + raise ValueError(f"ess_warn_fraction must lie in [0, 1], got {ess_warn_fraction}") + + targets = list(targets) + da = result.idata.posterior_predictive["forecast"] + G, t = build_moments(da.values, targets, result.var_names, result.steps) + weights, achieved = solve_tilt(G, t) + + n_chains, n_draws = da.shape[:2] + labels = [target_label(target) for target in targets] + counts = np.array([ + float(G[:, k].sum()) if isinstance(target, ProbabilityTarget) else np.nan for k, target in enumerate(targets) + ]) + diagnostics = tilt_diagnostics(weights, ess_warn_fraction, stacklevel=4) + + ds = xr.Dataset({ + "forecast": da, + "tilting_weights": xr.DataArray(weights.reshape(n_chains, n_draws), dims=["chain", "draw"]), + "achieved": xr.DataArray(achieved, dims=["target"], coords={"target": labels}), + "requested": xr.DataArray(t, dims=["target"], coords={"target": labels}), + "event_draws": xr.DataArray(counts, dims=["target"], coords={"target": labels}), + }) + ds.attrs.update(diagnostics) + return TiltedForecastResult( + idata=az.InferenceData(posterior_predictive=ds), + steps=result.steps, + var_names=list(result.var_names), + targets=targets, + ) + + +def tilt_diagnostics( + weights: np.ndarray, + ess_warn_fraction: float = ESS_WARN_FRACTION, + stacklevel: int = 3, +) -> dict[str, float]: + """Effective sample size and relative entropy, with the degeneracy warning. + + Args: + weights: Normalised tilting weights `(N,)`. + ess_warn_fraction: Warn below this ESS fraction. + stacklevel: `warnings.warn` stacklevel targeting the public caller. + + Returns: + Dict with `ess`, `ess_fraction`, and `kl_divergence`. + """ + n_total = weights.size + ess_value = ess(weights) + fraction = ess_value / n_total + if fraction < ess_warn_fraction: + warnings.warn( + f"The tilt is concentrated: effective sample size {ess_value:.1f} of {n_total} draws " + f"({fraction:.1%}, below the {ess_warn_fraction:.0%} threshold). Tilted summaries rest " + "on few draws and are noisy — relax the targets or draw a larger forecast sample.", + UserWarning, + stacklevel=stacklevel, + ) + return {"ess": ess_value, "ess_fraction": fraction, "kl_divergence": kl_divergence(weights)} diff --git a/src/impulso/identified.py b/src/impulso/identified.py index c4d92836..e9bddbd7 100644 --- a/src/impulso/identified.py +++ b/src/impulso/identified.py @@ -20,6 +20,7 @@ FEVDResult, HistoricalDecompositionResult, IRFResult, + ReverseStressResult, ScenarioResult, ) @@ -601,6 +602,170 @@ def counterfactual( idata = az.InferenceData(posterior_predictive=xr.Dataset({"counterfactual": cf_da, "actual": actual_da})) return CounterfactualResult(idata=idata, var_names=self.var_names) + def _validate_forecast_exog(self, steps: int, exog_future: np.ndarray | None) -> np.ndarray | None: + """Validate a forecast-side `exog_future` against the posterior. + + Shared by every forecast-side method on this object: the data must + not carry exogenous regressors the estimator never consumed, an + explicit `exog_future` needs a `B_exog` to multiply and the right + shape, and a model with exogenous data cannot forecast without one. + + Args: + steps: Forecast horizon. + exog_future: Future exogenous values, or None. + + Returns: + The coerced float array, or None. + + Raises: + ValueError: On any of the mismatches above. + """ + posterior = self.idata.posterior + if self.data.exog is not None and "B_exog" not in posterior: + raise ValueError( + "This IdentifiedVAR's data carries exogenous regressors the estimator " + "never consumed (no B_exog in the posterior); refit with an estimator " + "that supports them before scenario analysis." + ) + if exog_future is not None: + if "B_exog" not in posterior: + raise ValueError("exog_future provided but the posterior carries no B_exog.") + exog_future = np.asarray(exog_future, dtype=float) + n_exog = posterior["B_exog"].shape[-1] + if exog_future.shape != (steps, n_exog): + raise ValueError(f"exog_future must have shape ({steps}, {n_exog}), got {exog_future.shape}.") + if self.data.exog is not None and exog_future is None: + raise ValueError("exog_future is required when the model includes exogenous variables") + return exog_future + + def reverse_stress( + self, + variable: str, + threshold: float, + steps: int, + horizon: int | None = None, + probability: float = 1.0, + direction: Literal["below", "above"] = "below", + seed: int | np.random.Generator | None = None, + exog_future: np.ndarray | None = None, + ) -> "ReverseStressResult": + """Reverse stress test: which shocks would deliver this outcome? + + Ordinary scenario analysis runs forwards — you name the shocks and + read off the outcome. This runs backwards: you name the outcome + (`variable` crossing `threshold` at `horizon`) and read off the + *shock cocktail* that delivers it. Draw an unconditional density + forecast together with the structural shocks behind it + (`ADR-0009`), reweight the draws so the stress event carries + `probability` (entropic tilting; the default 1.0 is exact + conditioning on the event), and report the tilted-weighted mean of + the retained shocks — the average structural configuration among + the draws that produced the outcome. + + Because the cocktail averages *realised* draws rather than solving + a projection problem, it inherits the model's own shock + correlations and needs no arbitrary norm choice. Its magnitude + `q = ‖E_w[ε]‖²` is in the same one-standard-deviation units as + the scenario plausibility statistic, so a cocktail of total size 9 + is "a 3-sd configuration". + + The result's `q_cal` applies the ADPRR binomial calibration to + the tilt's relative entropy, `q_cal = (1 + sqrt(1 - exp(-2·KL/d))) + / 2` with `d = steps · n_vars`. This is an *extension* of that + calibration to soft conditioning, not a result from the paper: it + substitutes the (now finite) entropic divergence for the + hard-conditioning `z = q/2` that ADR-0005 documents as infinite. + + Note: + Under time-varying volatility the forecast factors are built + per simulated volatility path (conditional-on-path, as in + `structural_scenario`), and `SignRestriction` is not supported + there — the scheme re-samples rotations per call, so no single + structural coordinate system spans the forecast steps. + + Args: + variable: The endogenous variable to stress. + threshold: The threshold it must cross, in its own units. + steps: Number of forecast steps to simulate. + horizon: The 1-based step the event refers to. Defaults to + `steps` (the end of the forecast). + probability: Requested probability of the stress event, + `0 < p <= 1`. The default 1.0 conditions on it outright; + a smaller value softens the conditioning and keeps more + of the sample. + direction: `"below"` for `variable < threshold` (default) or + `"above"`. + seed: RNG seed (int) or Generator. Matched seeds reproduce + `structural_scenario`'s draws exactly. + exog_future: Future exogenous values, shape `(steps, k)`. + Required if the posterior carries `B_exog`. + + Returns: + ReverseStressResult with the conditioned forecast draws, the + structural shocks, the tilting weights, and the cocktail. + + Raises: + ValueError: On an unknown variable, a horizon outside + `1..steps`, a probability outside `(0, 1]`, a stress event + no draw satisfies, `SignRestriction` under time-varying + volatility, or exogenous-data mismatches. + """ + from impulso._scenario import structural_forecast_draws + from impulso._tilting import build_moments, solve_tilt, tilt_diagnostics + from impulso.results import ReverseStressResult + from impulso.scenario import ProbabilityTarget + + horizon = steps if horizon is None else horizon + if horizon > steps: + raise ValueError(f"horizon must lie in 1..steps, got horizon={horizon} with steps={steps}") + target = ProbabilityTarget( + variable=variable, + horizon=horizon, + threshold=threshold, + probability=probability, + direction=direction, + ) + exog_future = self._validate_forecast_exog(steps, exog_future) + + paths, eps = structural_forecast_draws(self, steps, seed=seed, exog_future=exog_future) + G, t = build_moments(paths, [target], self.var_names, steps) + weights, achieved = solve_tilt(G, t) + diagnostics = tilt_diagnostics(weights, stacklevel=3) + + n_chains, n_draws = paths.shape[:2] + n_total = n_chains * n_draws + cocktail = np.einsum("i,ihj->hj", weights, eps.reshape(n_total, *eps.shape[2:])) + q = float(np.sum(cocktail**2)) + d_total = steps * len(self.var_names) + q_cal = float((1.0 + np.sqrt(1.0 - np.exp(-2.0 * diagnostics["kl_divergence"] / d_total))) / 2.0) + + ds = xr.Dataset({ + "forecast": xr.DataArray( + paths, dims=["chain", "draw", "step", "variable"], coords={"variable": self.var_names} + ), + "structural_shocks": xr.DataArray( + eps, dims=["chain", "draw", "step", "shock"], coords={"shock": self.shock_names} + ), + "tilting_weights": xr.DataArray(weights.reshape(n_chains, n_draws), dims=["chain", "draw"]), + "shock_cocktail": xr.DataArray(cocktail, dims=["step", "shock"], coords={"shock": self.shock_names}), + }) + ds.attrs.update(diagnostics) + ds.attrs["baseline_probability"] = float(G[:, 0].mean()) + ds.attrs["achieved_probability"] = float(achieved[0]) + ds.attrs["q"] = q + ds.attrs["q_cal"] = q_cal + return ReverseStressResult( + idata=az.InferenceData(posterior_predictive=ds), + steps=steps, + var_names=self.var_names, + shock_names=self.shock_names, + variable=variable, + threshold=float(threshold), + horizon=horizon, + direction=direction, + probability=float(probability), + ) + def structural_scenario( self, steps: int, @@ -710,22 +875,7 @@ def structural_scenario( ) if path_uncertainty not in ("none", "unconditional"): raise ValueError(f"path_uncertainty must be 'none' or 'unconditional', got {path_uncertainty!r}") - posterior = self.idata.posterior - if self.data.exog is not None and "B_exog" not in posterior: - raise ValueError( - "This IdentifiedVAR's data carries exogenous regressors the estimator " - "never consumed (no B_exog in the posterior); refit with an estimator " - "that supports them before scenario analysis." - ) - if exog_future is not None: - if "B_exog" not in posterior: - raise ValueError("exog_future provided but the posterior carries no B_exog.") - exog_future = np.asarray(exog_future, dtype=float) - n_exog = posterior["B_exog"].shape[-1] - if exog_future.shape != (steps, n_exog): - raise ValueError(f"exog_future must have shape ({steps}, {n_exog}), got {exog_future.shape}.") - if self.data.exog is not None and exog_future is None: - raise ValueError("exog_future is required when the model includes exogenous variables") + exog_future = self._validate_forecast_exog(steps, exog_future) paths, q, q_cond, q_cal, r = structural_scenario_engine( self, diff --git a/src/impulso/plotting/__init__.py b/src/impulso/plotting/__init__.py index 94a84743..5fefe29d 100644 --- a/src/impulso/plotting/__init__.py +++ b/src/impulso/plotting/__init__.py @@ -10,6 +10,7 @@ from impulso.plotting._structural_scenario import plot_structural_scenario from impulso.plotting._sv_forecast import plot_sv_forecast from impulso.plotting._sv_volatility import plot_volatility +from impulso.plotting._tilted_forecast import plot_reverse_stress, plot_tilted_forecast __all__ = [ "plot_conditional_forecast", @@ -19,7 +20,9 @@ "plot_forecast", "plot_historical_decomposition", "plot_irf", + "plot_reverse_stress", "plot_structural_scenario", "plot_sv_forecast", + "plot_tilted_forecast", "plot_volatility", ] diff --git a/src/impulso/plotting/_tilted_forecast.py b/src/impulso/plotting/_tilted_forecast.py new file mode 100644 index 00000000..c5e1432f --- /dev/null +++ b/src/impulso/plotting/_tilted_forecast.py @@ -0,0 +1,146 @@ +"""Plotting for entropically tilted forecasts and reverse stress tests.""" + +from typing import TYPE_CHECKING + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.figure import Figure + +if TYPE_CHECKING: + from impulso.results import ReverseStressResult, TiltedForecastResult + + +def _tilt_subtitle(result: "TiltedForecastResult | ReverseStressResult") -> str: + """Effective-sample-size and relative-entropy summary for a suptitle.""" + attrs = result.idata.posterior_predictive.attrs + n_draws = result.idata.posterior_predictive["tilting_weights"].size + return f"ESS {attrs['ess']:.0f} of {n_draws} draws, KL {attrs['kl_divergence']:.3f}" + + +def plot_tilted_forecast( + result: "TiltedForecastResult", + prob: float = 0.89, + figsize: tuple[float, float] | None = None, +) -> Figure: + """Plot the tilted forecast against the untilted median, one panel per variable. + + The weighted median and its weighted HDI band show the tilted + distribution; the untilted median is overlaid dashed so the effect of + the targets is visible directly. + + Args: + result: TiltedForecastResult. + prob: Probability mass for the HDI band. Default 0.89. + figsize: Figure size. Defaults to `(12, 3 * n_vars)`. + + Returns: + Matplotlib Figure. + """ + med = result.median() + base = result.base_median() + hdi = result.hdi(prob=prob) + steps_axis = np.arange(1, result.steps + 1) + n_vars = len(result.var_names) + + if figsize is None: + figsize = (12, 3 * n_vars) + + fig, axes = plt.subplots(n_vars, 1, figsize=figsize, sharex=True, squeeze=False) + axes = axes[:, 0] + fig.suptitle(f"Tilted Forecast ({_tilt_subtitle(result)})") + + for i, var in enumerate(result.var_names): + axes[i].plot(steps_axis, med[var].values, color="C0", linewidth=1.2, label="tilted median") + axes[i].fill_between( + steps_axis, + hdi.lower[var].values, + hdi.upper[var].values, + color="C0", + alpha=0.25, + linewidth=0, + label=f"{int(prob * 100)}% HDI (tilted)", + ) + axes[i].plot( + steps_axis, + base[var].values, + color="grey", + linewidth=1.0, + linestyle="--", + label="untilted median", + ) + axes[i].set_ylabel(var) + if i == 0: + axes[i].legend(fontsize=8, loc="upper right") + + axes[-1].set_xlabel("Step") + fig.tight_layout() + return fig + + +def plot_reverse_stress( + result: "ReverseStressResult", + prob: float = 0.89, + figsize: tuple[float, float] | None = None, +) -> Figure: + """Plot the stressed variable's conditioned fan and the shock cocktail. + + The top panel shows the stressed variable under the event-conditioned + weights, with the threshold marked and the untilted median dashed for + reference. One bar panel per structural shock then shows the + cocktail — the tilted-weighted mean shock path, in + one-standard-deviation units. + + Args: + result: ReverseStressResult. + prob: Probability mass for the HDI band. Default 0.89. + figsize: Figure size. Defaults to `(12, 3 * (1 + n_shocks))`. + + Returns: + Matplotlib Figure. + """ + cocktail = result.shock_cocktail() + n_shocks = len(result.shock_names) + steps_axis = np.arange(1, result.steps + 1) + + if figsize is None: + figsize = (12, 3 * (1 + n_shocks)) + + fig, axes = plt.subplots(1 + n_shocks, 1, figsize=figsize, sharex=True, squeeze=False) + axes = axes[:, 0] + attrs = result.idata.posterior_predictive.attrs + sign = "<" if result.direction == "below" else ">" + fig.suptitle( + f"Reverse Stress: P({result.variable}[h={result.horizon}] {sign} {result.threshold:g}) " + f"= {attrs['achieved_probability']:.2f} (baseline {attrs['baseline_probability']:.2f}; " + f"{_tilt_subtitle(result)}; q = {attrs['q']:.2f})" + ) + + med = result.median() + base = result.base_median() + hdi = result.hdi(prob=prob) + var = result.variable + axes[0].plot(steps_axis, med[var].values, color="C3", linewidth=1.2, label="conditioned median") + axes[0].fill_between( + steps_axis, + hdi.lower[var].values, + hdi.upper[var].values, + color="C3", + alpha=0.25, + linewidth=0, + label=f"{int(prob * 100)}% HDI (conditioned)", + ) + axes[0].plot(steps_axis, base[var].values, color="grey", linewidth=1.0, linestyle="--", label="untilted median") + axes[0].axhline(result.threshold, color="black", linewidth=0.8, linestyle=":", label="threshold") + axes[0].scatter([result.horizon], [result.threshold], color="black", marker="x", s=30, zorder=3) + axes[0].set_ylabel(var) + axes[0].legend(fontsize=8, loc="upper right") + + for j, shock in enumerate(result.shock_names): + ax = axes[1 + j] + ax.bar(steps_axis, cocktail[shock].values, color="C1", width=0.7) + ax.axhline(0.0, color="black", linewidth=0.8) + ax.set_ylabel(f"{shock}\n(sd units)") + + axes[-1].set_xlabel("Step") + fig.tight_layout() + return fig diff --git a/src/impulso/results.py b/src/impulso/results.py index 26a880e6..4d586e11 100644 --- a/src/impulso/results.py +++ b/src/impulso/results.py @@ -11,7 +11,10 @@ from pydantic import Field, model_validator from impulso._base import ImpulsoBaseModel -from impulso.scenario import ShockPath, VariablePath +from impulso.scenario import MomentTarget, ProbabilityTarget, ShockPath, VariablePath + +# Targets accepted by the tilting entry points. +Target = ProbabilityTarget | MomentTarget def _wide_frame(da: xr.DataArray, row_dim: str, col_dim: str = "shock") -> pd.DataFrame: @@ -160,6 +163,30 @@ def plot(self) -> Figure: return plot_forecast(self) + def tilt(self, targets: list[Target], ess_warn_fraction: float = 0.1) -> "TiltedForecastResult": + """Reweight these draws to satisfy distributional targets (entropic tilting). + + See `TiltedForecastResult` for what comes back and + `ADR-0009` for why this is a post-hoc reweighting layer rather + than a re-solve. + + Args: + targets: `ProbabilityTarget` / `MomentTarget` list. + ess_warn_fraction: Warn when the effective sample size falls + below this fraction of the draw count. Default 0.1. + + Returns: + TiltedForecastResult carrying these draws by reference plus + the tilting weights and diagnostics. + + Raises: + ValueError: If this is a mean forecast, or if the targets are + unachievable by reweighting these draws. + """ + from impulso._tilting import tilt_result + + return tilt_result(self, list(targets), ess_warn_fraction) + class IRFResult(VARResultBase): """Result from impulse response function computation. @@ -475,6 +502,31 @@ def plot(self) -> Figure: return plot_conditional_forecast(self) + def tilt(self, targets: list[Target], ess_warn_fraction: float = 0.1) -> "TiltedForecastResult": + """Reweight these draws to satisfy distributional targets (entropic tilting). + + Chaining hard conditioning with soft targets is the supported way + to mix them: the pins hold pathwise on *every* draw here, and + reweighting never moves a draw, so the pins survive the tilt + exactly — a theorem, not a code path. + + Args: + targets: `ProbabilityTarget` / `MomentTarget` list. + ess_warn_fraction: Warn when the effective sample size falls + below this fraction of the draw count. Default 0.1. + + Returns: + TiltedForecastResult carrying these draws by reference plus + the tilting weights and diagnostics. + + Raises: + ValueError: If this is a mean forecast, or if the targets are + unachievable by reweighting these draws. + """ + from impulso._tilting import tilt_result + + return tilt_result(self, list(targets), ess_warn_fraction) + class ScenarioResult(ConditionalForecastResult): """Result from structural scenario analysis. @@ -504,6 +556,199 @@ def plot(self) -> Figure: return plot_structural_scenario(self) +class _WeightedResultMixin: + """Shared weighted summaries for tilting-derived results. + + Both tilted forecasts and reverse-stress results summarise the same + `"forecast"` draws under a `"tilting_weights"` variable, so the + weighted median / HDI / DataFrame surface is written once here. + """ + + def _weights_flat(self) -> np.ndarray: + """Normalised tilting weights flattened over `(chain, draw)`.""" + return self.idata.posterior_predictive["tilting_weights"].values.ravel() + + def _forecast_flat(self) -> np.ndarray: + """Forecast draws reshaped to `(N, steps, n_vars)`.""" + da = self.idata.posterior_predictive["forecast"] + n_chains, n_draws = da.shape[:2] + return da.values.reshape(n_chains * n_draws, *da.shape[2:]) + + @property + def weights(self) -> np.ndarray: + """Tilting weights, shape `(chain, draw)`, summing to 1.""" + return self.idata.posterior_predictive["tilting_weights"].values + + def median(self) -> pd.DataFrame: + """Weighted posterior median forecast (step-indexed).""" + from impulso._tilting import weighted_quantile + + med = weighted_quantile(self._forecast_flat(), self._weights_flat(), 0.5) + df = pd.DataFrame(np.asarray(med), columns=self.var_names) + df.index.name = "step" + return df + + def base_median(self) -> pd.DataFrame: + """Untilted posterior median, for comparison against `median()`.""" + med = self.idata.posterior_predictive["forecast"].median(dim=("chain", "draw")).values + df = pd.DataFrame(med, columns=self.var_names) + df.index.name = "step" + return df + + def hdi(self, prob: float = 0.89) -> HDIResult: + """Weighted highest-density interval for the forecast. + + Args: + prob: Probability mass for the HDI. Default 0.89. + + Returns: + HDIResult whose `lower` / `upper` DataFrames mirror `median()`. + """ + from impulso._tilting import weighted_hdi + + lower, upper = weighted_hdi(self._forecast_flat(), self._weights_flat(), prob) + return HDIResult( + lower=pd.DataFrame(lower, columns=self.var_names), + upper=pd.DataFrame(upper, columns=self.var_names), + prob=prob, + ) + + def to_dataframe(self) -> pd.DataFrame: + """Weighted posterior median as a DataFrame (passthrough to `median()`).""" + return self.median() + + +class TiltedForecastResult(_WeightedResultMixin, VARResultBase): + """Forecast draws reweighted by entropic tilting (ADR-0009). + + The posterior-predictive Dataset carries `"forecast"` — the parent + result's draws, held by reference, never copied or moved — plus + `"tilting_weights"` `(chain, draw)` and the per-target + `"requested"` / `"achieved"` / `"event_draws"` vectors over a + `target` coordinate. Dataset attrs hold `ess`, `ess_fraction`, and + `kl_divergence`. + + Every summary on this object is weighted: `median()` and `hdi()` + read the tilted distribution, while `base_median()` returns the + untilted median so the two can be compared directly. + + Attributes: + idata: InferenceData with the parent draws and the weights. + steps: Number of forecast steps. + var_names: Names of forecasted variables. + targets: The targets echoed from the call. + """ + + steps: int + var_names: list[str] + targets: list[Target] = Field(default_factory=list, repr=False) + + def summary(self) -> dict[str, object]: + """Diagnostics and per-target achievement. + + Returns: + Dict with `ess`, `ess_fraction`, `kl_divergence`, `n_draws`, + and a `targets` list of per-target dicts holding `target`, + `requested`, `achieved`, and `draws_in_event` (`None` for + moment targets). + """ + pp = self.idata.posterior_predictive + rows = [] + for k, label in enumerate(pp["target"].values.tolist()): + count = float(pp["event_draws"].values[k]) + rows.append({ + "target": label, + "requested": float(pp["requested"].values[k]), + "achieved": float(pp["achieved"].values[k]), + "draws_in_event": None if np.isnan(count) else int(count), + }) + return { + "ess": float(pp.attrs["ess"]), + "ess_fraction": float(pp.attrs["ess_fraction"]), + "kl_divergence": float(pp.attrs["kl_divergence"]), + "n_draws": int(pp["tilting_weights"].size), + "targets": rows, + } + + def plot(self) -> Figure: + """Plot the tilted fan chart against the untilted median.""" + from impulso.plotting import plot_tilted_forecast + + return plot_tilted_forecast(self) + + +class ReverseStressResult(_WeightedResultMixin, VARResultBase): + """Shock cocktail behind a stress event, from reverse stress testing. + + The posterior-predictive Dataset carries `"forecast"` + (chain, draw, step, variable), the structural shocks that generated + those draws (`"structural_shocks"`, chain, draw, step, shock), the + `"tilting_weights"` that condition on the event, and the + `"shock_cocktail"` (step, shock) — the tilted-weighted mean of the + retained structural shocks, in one-standard-deviation units. Dataset + attrs hold `baseline_probability`, `achieved_probability`, `ess`, + `ess_fraction`, `kl_divergence`, `q`, and `q_cal`. + + Attributes: + idata: InferenceData with the draws, weights, and cocktail. + steps: Number of forecast steps. + var_names: Names of forecasted variables. + shock_names: Structural shock coordinate labels. + variable: The stressed variable. + threshold: The stress threshold, in the variable's units. + horizon: The 1-based forecast step the event refers to. + direction: `"below"` or `"above"`. + probability: Requested probability of the stress event. + """ + + steps: int + var_names: list[str] + shock_names: list[str] + variable: str + threshold: float + horizon: int + direction: Literal["below", "above"] = "below" + probability: float = 1.0 + + def shock_cocktail(self) -> pd.DataFrame: + """The shock cocktail as a step-indexed DataFrame. + + Returns: + DataFrame indexed by forecast step (1-based) with one column + per structural shock, in one-standard-deviation units. + """ + da = self.idata.posterior_predictive["shock_cocktail"] + df = pd.DataFrame(da.values, columns=self.shock_names, index=pd.RangeIndex(1, self.steps + 1, name="step")) + return df + + def summary(self) -> dict[str, float]: + """Event probabilities, tilt diagnostics, and cocktail plausibility. + + Returns: + Dict with `baseline_probability`, `requested_probability`, + `achieved_probability`, `ess`, `ess_fraction`, + `kl_divergence`, `n_draws`, `q`, and `q_cal`. + """ + pp = self.idata.posterior_predictive + return { + "baseline_probability": float(pp.attrs["baseline_probability"]), + "requested_probability": float(self.probability), + "achieved_probability": float(pp.attrs["achieved_probability"]), + "ess": float(pp.attrs["ess"]), + "ess_fraction": float(pp.attrs["ess_fraction"]), + "kl_divergence": float(pp.attrs["kl_divergence"]), + "n_draws": int(pp["tilting_weights"].size), + "q": float(pp.attrs["q"]), + "q_cal": float(pp.attrs["q_cal"]), + } + + def plot(self) -> Figure: + """Plot the stressed variable's tilted fan and the shock cocktail.""" + from impulso.plotting import plot_reverse_stress + + return plot_reverse_stress(self) + + class CounterfactualResult(VARResultBase): """Historical counterfactual paths alongside the actual data. diff --git a/src/impulso/scenario.py b/src/impulso/scenario.py index 5187a67a..013e46e1 100644 --- a/src/impulso/scenario.py +++ b/src/impulso/scenario.py @@ -2,18 +2,22 @@ Typed, frozen spec objects expressing scenario content (ADR-0005): `ShockPath` sets a structural shock's path, `VariablePath` pins a future -endogenous path. Each scenario method accepts only the condition types +endogenous path, and the *targets* — `ProbabilityTarget`, `MomentTarget` — +state distributional facts a forecast should satisfy after entropic +tilting (ADR-0009). Each scenario method accepts only the condition types that are legal for it, so illegal combinations are unrepresentable rather than validated away. """ from __future__ import annotations +from typing import Literal + import numpy as np import pandas as pd from pydantic import field_validator, model_validator -from impulso._base import ImpulsoBaseModel +from impulso._base import ImpulsoBaseModel, ImpulsoModel def _coerce_values(value: float | np.ndarray) -> float | np.ndarray: @@ -105,3 +109,94 @@ class VariablePath(ImpulsoBaseModel): @classmethod def _validate_values(cls, value: float | np.ndarray) -> float | np.ndarray: return _coerce_values(value) + + +def _check_finite(value: float, field: str) -> float: + """Reject NaN/inf on a target field.""" + if not np.isfinite(value): + raise ValueError(f"{field} must be finite, got {value}") + return float(value) + + +def _check_horizon(value: int) -> int: + """Reject non-positive horizons (steps are 1-based, as on `VariablePath`).""" + if value < 1: + raise ValueError(f"horizon is 1-based (step 1 is the first forecast step) and must be >= 1, got {value}") + return value + + +class ProbabilityTarget(ImpulsoModel): + """Require an event to carry a given probability after tilting. + + The event is `{y[horizon] < threshold}` (`direction="below"`, the + default) or `{y[horizon] > threshold}` — strict inequalities, so a + draw sitting exactly on the threshold is outside the event. Tilting + reweights the existing forecast draws to hit `probability` while + staying as close as possible (in relative entropy) to the untilted + forecast; it never moves a draw, so nothing the draws already satisfy + can be broken. `probability=1.0` is pure conditioning: draws outside + the event get weight zero. + + Attributes: + variable: Name of the endogenous variable the event refers to. + horizon: Forecast step the event refers to, 1-based (step 1 is + the first forecast step, matching `VariablePath`). + threshold: The event's threshold, in the variable's own units. + probability: Requested probability of the event, `0 < p <= 1`. + direction: `"below"` for `y < threshold` (default) or `"above"` + for `y > threshold`. + """ + + variable: str + horizon: int + threshold: float + probability: float + direction: Literal["below", "above"] = "below" + + @field_validator("horizon") + @classmethod + def _validate_horizon(cls, value: int) -> int: + return _check_horizon(value) + + @field_validator("threshold") + @classmethod + def _validate_threshold(cls, value: float) -> float: + return _check_finite(value, "threshold") + + @field_validator("probability") + @classmethod + def _validate_probability(cls, value: float) -> float: + value = _check_finite(value, "probability") + if not 0.0 < value <= 1.0: + raise ValueError(f"probability must satisfy 0 < p <= 1, got {value}") + return value + + +class MomentTarget(ImpulsoModel): + """Require a variable's tilted forecast mean at one horizon. + + The moment condition is `E_w[y[horizon]] = mean`, imposed on the + existing draws by reweighting. The requested mean must lie strictly + inside the range the draws already span — tilting cannot move mass + where there is none. + + Attributes: + variable: Name of the endogenous variable. + horizon: Forecast step, 1-based (step 1 is the first forecast + step, matching `VariablePath`). + mean: Requested tilted mean, in the variable's own units. + """ + + variable: str + horizon: int + mean: float + + @field_validator("horizon") + @classmethod + def _validate_horizon(cls, value: int) -> int: + return _check_horizon(value) + + @field_validator("mean") + @classmethod + def _validate_mean(cls, value: float) -> float: + return _check_finite(value, "mean") diff --git a/tests/test_reverse_stress.py b/tests/test_reverse_stress.py new file mode 100644 index 00000000..3a051171 --- /dev/null +++ b/tests/test_reverse_stress.py @@ -0,0 +1,294 @@ +"""Invariant tests for IdentifiedVAR.reverse_stress (issue #150). + +Core identities: the base draws nest inside `structural_scenario` under a +matched seed (the RNG stream contract of `structural_forecast_draws`); the +cocktail at `probability=1.0` is exactly the plain mean of the retained +shocks; and with a deterministic posterior the cocktail reproduces the +textbook truncated-normal mean. +""" + +import matplotlib + +matplotlib.use("Agg") + +import arviz as az +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from matplotlib.figure import Figure + +from impulso._scenario import structural_forecast_draws +from impulso.data import VARData +from impulso.fitted import FittedVAR +from impulso.identification import Cholesky +from impulso.volatility import Constant + + +@pytest.fixture +def identified_2v(synthetic_idata_2v, var_data_2v): + """Cholesky-identified 2-var VAR(1) from the synthetic posterior.""" + fitted = FittedVAR( + idata=synthetic_idata_2v, + n_lags=1, + data=var_data_2v, + var_names=["y1", "y2"], + volatility=Constant(), + ) + return fitted.set_identification_strategy(Cholesky(ordering=["y1", "y2"])) + + +def _deterministic_fitted(n_draws=4000): + """A posterior whose draws all carry the same parameters. + + Every draw shares one `B`, `intercept` and `L`, so the only randomness + left in a forecast is the structural shocks themselves — which turns + the reverse-stress cocktail into a textbook truncated-normal mean. + """ + n_vars = 2 + B = np.tile(np.array([[0.5, 0.1], [-0.2, 0.3]]), (1, n_draws, 1, 1)) + intercept = np.tile(np.array([0.1, -0.05]), (1, n_draws, 1)) + L_single = np.linalg.cholesky(np.array([[1.0, 0.3], [0.3, 0.8]])) + L = np.tile(L_single, (1, n_draws, 1, 1)) + posterior = xr.Dataset({ + "B": (("chain", "draw", "var", "coeff"), B), + "intercept": (("chain", "draw", "var"), intercept), + "L": (("chain", "draw", "var1", "var2"), L), + }) + data = VARData( + endog=np.random.default_rng(1).standard_normal((12, n_vars)), + endog_names=["y1", "y2"], + index=pd.date_range("2000-01-01", periods=12, freq="QS"), + ) + return FittedVAR( + idata=az.InferenceData(posterior=posterior), + n_lags=1, + data=data, + var_names=["y1", "y2"], + volatility=Constant(), + ) + + +def _median_threshold(identified, steps, horizon, seed, variable_index=0): + """Empirical median of the matched-seed forecast at `horizon`.""" + paths, _ = structural_forecast_draws(identified, steps, seed=seed) + return float(np.median(paths[:, :, horizon - 1, variable_index])), paths + + +class TestMatchedSeedNesting: + def test_base_draws_match_structural_scenario_under_a_shared_seed(self, identified_2v): + paths, eps = structural_forecast_draws(identified_2v, 6, seed=123) + scenario = identified_2v.structural_scenario(6, seed=123) + assert np.allclose(paths, scenario.idata.posterior_predictive["forecast"].values, atol=1e-8) + assert eps.shape == paths.shape + + def test_reverse_stress_carries_those_same_draws(self, identified_2v): + threshold, paths = _median_threshold(identified_2v, 6, 6, seed=123) + result = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=6, seed=123) + assert np.allclose(result.idata.posterior_predictive["forecast"].values, paths, atol=1e-8) + + def test_shocks_reproduce_the_paths_through_the_scenario_engine(self, identified_2v): + # The scenario engine with every shock prescribed at the drawn + # values must reproduce the same paths, which pins the returned + # eps to the ones that actually generated the forecast. + from impulso.scenario import ShockPath + + paths, eps = structural_forecast_draws(identified_2v, 3, seed=7) + prescribed = identified_2v.structural_scenario( + 3, + shocks=[ShockPath(shock=name, values=eps[0, 0, :, j]) for j, name in enumerate(identified_2v.shock_names)], + seed=7, + ) + # Draw (0, 0) had exactly those shocks prescribed, so its path is unchanged. + assert np.allclose(prescribed.idata.posterior_predictive["forecast"].values[0, 0], paths[0, 0], atol=1e-10) + + +class TestCocktail: + def test_cocktail_is_the_plain_mean_of_the_retained_shocks(self, identified_2v): + threshold, paths = _median_threshold(identified_2v, 5, 5, seed=42) + _, eps = structural_forecast_draws(identified_2v, 5, seed=42) + result = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=5, seed=42) + + n_total = paths.shape[0] * paths.shape[1] + retained = paths[:, :, 4, 0].reshape(n_total) < threshold + oracle = eps.reshape(n_total, 5, 2)[retained].mean(axis=0) + assert np.allclose(result.shock_cocktail().values, oracle, atol=1e-14) + + def test_q_is_the_squared_norm_of_the_cocktail(self, identified_2v): + threshold, _ = _median_threshold(identified_2v, 5, 5, seed=42) + result = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=5, seed=42) + cocktail = result.shock_cocktail().values + assert result.summary()["q"] == pytest.approx(float(np.sum(cocktail**2)), rel=1e-13) + + def test_cocktail_is_labelled_by_shock_names(self, identified_2v): + threshold, _ = _median_threshold(identified_2v, 4, 4, seed=5) + result = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=4, seed=5) + cocktail = result.shock_cocktail() + assert list(cocktail.columns) == identified_2v.shock_names + assert list(cocktail.index) == [1, 2, 3, 4] + assert cocktail.index.name == "step" + + def test_soft_conditioning_shrinks_the_cocktail_toward_zero(self, identified_2v): + threshold, _ = _median_threshold(identified_2v, 4, 4, seed=8) + hard = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=4, seed=8) + soft = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=4, probability=0.6, seed=8) + assert soft.summary()["q"] < hard.summary()["q"] + assert soft.summary()["ess"] > hard.summary()["ess"] + + +class TestProbabilities: + def test_threshold_at_the_empirical_median_gives_a_half_baseline(self, identified_2v): + threshold, _ = _median_threshold(identified_2v, 5, 5, seed=13) + result = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=5, seed=13) + summary = result.summary() + assert summary["baseline_probability"] == pytest.approx(0.5, abs=0.02) + assert summary["achieved_probability"] == pytest.approx(1.0, abs=1e-12) + + def test_a_certain_event_leaves_the_weights_uniform(self, identified_2v): + result = identified_2v.reverse_stress(variable="y1", threshold=-1e6, steps=4, direction="above", seed=14) + summary = result.summary() + assert summary["baseline_probability"] == pytest.approx(1.0, abs=1e-12) + assert summary["kl_divergence"] == pytest.approx(0.0, abs=1e-12) + assert summary["q_cal"] == pytest.approx(0.5, abs=1e-12) + assert summary["ess"] == pytest.approx(float(result.weights.size), rel=1e-12) + # Uniform weights average all the shocks, so the cocktail is Monte + # Carlo noise around zero rather than a stress configuration. + _, eps = structural_forecast_draws(identified_2v, 4, seed=14) + assert np.allclose(result.shock_cocktail().values, eps.mean(axis=(0, 1)), atol=1e-14) + assert summary["q"] < 0.5 + + def test_horizon_defaults_to_the_final_step(self, identified_2v): + threshold, _ = _median_threshold(identified_2v, 5, 5, seed=15) + default = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=5, seed=15) + explicit = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=5, horizon=5, seed=15) + assert default.horizon == 5 + assert np.allclose(default.shock_cocktail().values, explicit.shock_cocktail().values) + + +class TestSignSanity: + """A deterministic posterior turns the cocktail into a textbook oracle.""" + + def test_impact_shock_matches_the_truncated_normal_mean(self): + fitted = _deterministic_fitted() + identified = fitted.set_identification_strategy(Cholesky(ordering=["y1", "y2"])) + # With identical parameters across draws, y1 at step 1 is + # b1 + P[0, 0] * eps_1[y1]; thresholding at b1 selects the draws + # with a negative y1 shock on impact. + b1 = float(fitted.forecast(4, include_shock_uncertainty=False).median()["y1"].iloc[0]) + result = identified.reverse_stress(variable="y1", threshold=b1, steps=4, horizon=1, seed=99) + cocktail = result.shock_cocktail() + + assert cocktail.loc[1, "y1"] == pytest.approx(-np.sqrt(2.0 / np.pi), abs=0.05) + # Under a natural-order Cholesky, y1 on impact loads on its own + # shock only, so nothing else is selected. + assert abs(cocktail.loc[1, "y2"]) < 0.1 + assert np.all(np.abs(cocktail.loc[2:].values) < 0.1) + assert result.summary()["baseline_probability"] == pytest.approx(0.5, abs=0.02) + + def test_direction_above_flips_the_cocktail_sign(self): + fitted = _deterministic_fitted() + identified = fitted.set_identification_strategy(Cholesky(ordering=["y1", "y2"])) + b1 = float(fitted.forecast(4, include_shock_uncertainty=False).median()["y1"].iloc[0]) + result = identified.reverse_stress(variable="y1", threshold=b1, steps=4, horizon=1, direction="above", seed=99) + assert result.shock_cocktail().loc[1, "y1"] == pytest.approx(np.sqrt(2.0 / np.pi), abs=0.05) + + +class TestGuards: + def test_empty_support_raises_with_the_draw_counts(self, identified_2v): + with pytest.raises(ValueError, match=r"0 of 100 draws satisfy"): + identified_2v.reverse_stress(variable="y1", threshold=-1e6, steps=4, seed=1) + + def test_horizon_beyond_steps_raises(self, identified_2v): + with pytest.raises(ValueError, match=r"horizon must lie in 1\.\.steps"): + identified_2v.reverse_stress(variable="y1", threshold=0.0, steps=4, horizon=9, seed=1) + + def test_zero_horizon_raises(self, identified_2v): + with pytest.raises(ValueError, match="1-based"): + identified_2v.reverse_stress(variable="y1", threshold=0.0, steps=4, horizon=0, seed=1) + + def test_unknown_variable_raises(self, identified_2v): + with pytest.raises(ValueError, match="Unknown variable"): + identified_2v.reverse_stress(variable="nope", threshold=0.0, steps=4, seed=1) + + def test_probability_outside_the_unit_interval_raises(self, identified_2v): + with pytest.raises(ValueError, match="0 < p <= 1"): + identified_2v.reverse_stress(variable="y1", threshold=0.0, steps=4, probability=1.5, seed=1) + + def test_exog_future_required_when_the_model_carries_exogenous_data(self, rng, var_data_2v, synthetic_idata_2v): + data = VARData( + endog=var_data_2v.endog, + endog_names=["y1", "y2"], + index=var_data_2v.index, + exog=rng.standard_normal((len(var_data_2v.index), 1)), + exog_names=["z"], + ) + posterior = synthetic_idata_2v.posterior.copy() + posterior["B_exog"] = xr.DataArray( + np.zeros((2, 50, 2, 1)), + dims=["chain", "draw", "var", "exog"], + ) + fitted = FittedVAR( + idata=az.InferenceData(posterior=posterior), + n_lags=1, + data=data, + var_names=["y1", "y2"], + volatility=Constant(), + ) + identified = fitted.set_identification_strategy(Cholesky(ordering=["y1", "y2"])) + with pytest.raises(ValueError, match="exog_future is required"): + identified.reverse_stress(variable="y1", threshold=0.0, steps=4, seed=1) + + def test_ess_warning_propagates(self, identified_2v): + paths, _ = structural_forecast_draws(identified_2v, 4, seed=21) + threshold = float(np.quantile(paths[:, :, 3, 0], 0.03)) + with pytest.warns(UserWarning, match="tilt is concentrated"): + result = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=4, seed=21) + assert result.summary()["ess_fraction"] < 0.1 + + +class TestResultSurface: + def test_summary_keys_are_pinned(self, identified_2v): + threshold, _ = _median_threshold(identified_2v, 4, 4, seed=31) + result = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=4, seed=31) + assert set(result.summary()) == { + "baseline_probability", + "requested_probability", + "achieved_probability", + "ess", + "ess_fraction", + "kl_divergence", + "n_draws", + "q", + "q_cal", + } + + def test_median_and_hdi_are_weighted(self, identified_2v): + threshold, _ = _median_threshold(identified_2v, 4, 4, seed=32) + result = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=4, seed=32) + med = result.median() + assert med.shape == (4, 2) + assert list(med.columns) == ["y1", "y2"] + # Conditioning on y1 being low pushes its conditioned median below + # the untilted one at the conditioning horizon. + assert med["y1"].iloc[3] < result.base_median()["y1"].iloc[3] + hdi = result.hdi(prob=0.5) + assert np.all(hdi.upper.values >= hdi.lower.values) + assert result.to_dataframe().equals(med) + + def test_structural_shocks_are_carried_on_the_result(self, identified_2v): + threshold, _ = _median_threshold(identified_2v, 4, 4, seed=33) + result = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=4, seed=33) + shocks = result.idata.posterior_predictive["structural_shocks"] + assert shocks.dims == ("chain", "draw", "step", "shock") + assert list(shocks.coords["shock"].values) == identified_2v.shock_names + + def test_plot_returns_a_figure(self, identified_2v): + threshold, _ = _median_threshold(identified_2v, 4, 4, seed=34) + result = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=4, seed=34) + assert isinstance(result.plot(), Figure) + + def test_result_is_frozen(self, identified_2v): + threshold, _ = _median_threshold(identified_2v, 4, 4, seed=35) + result = identified_2v.reverse_stress(variable="y1", threshold=threshold, steps=4, seed=35) + with pytest.raises(ValueError, match="frozen"): + result.threshold = 0.0 diff --git a/tests/test_structural_scenario.py b/tests/test_structural_scenario.py index abd90603..4a68ad21 100644 --- a/tests/test_structural_scenario.py +++ b/tests/test_structural_scenario.py @@ -17,7 +17,7 @@ from impulso._linalg import lag_matrices from impulso._propagate import propagate -from impulso._scenario import _resolve_adjusting +from impulso._scenario import _resolve_adjusting, structural_forecast_draws from impulso.fitted import FittedVAR from impulso.identification import Cholesky, SignRestriction from impulso.identified import IdentifiedVAR @@ -293,6 +293,19 @@ def test_matched_seed_nesting_under_generator_consuming_volatility(self, pair_sv atol=1e-8, ) + def test_structural_forecast_draws_nest_under_generator_consuming_volatility(self, pair_sv): + """structural_forecast_draws holds its RNG contract when the volatility path draws too. + + `_forecast_shock_matrices` consumes the generator only under + time-varying volatility, so this is the branch where a mis-ordered + stream would silently desynchronise the two code paths. + """ + _, identified = pair_sv + paths, eps = structural_forecast_draws(identified, 5, seed=23) + scn = identified.structural_scenario(steps=5, seed=23) + np.testing.assert_allclose(paths, scn.idata.posterior_predictive["forecast"].values, atol=1e-12) + assert eps.shape == paths.shape + def _single_draw_identified_exog(): """Single-draw exog posterior, Cholesky-identified.""" diff --git a/tests/test_tilting.py b/tests/test_tilting.py new file mode 100644 index 00000000..251761f2 --- /dev/null +++ b/tests/test_tilting.py @@ -0,0 +1,508 @@ +"""Invariant tests for entropic tilting (issue #150). + +The solver half is checked against analytic oracles: the closed-form +two-mass tilt and its ESS/KL, the Gaussian mean tilt whose dual solution +is known in closed form, and the primal problem solved independently by +SLSQP. The result half checks the guards, the weighted summaries, and the +theorem that hard pins survive any reweighting. +""" + +import matplotlib + +matplotlib.use("Agg") + +import numpy as np +import pytest + +from impulso._tilting import ( + build_moments, + ess, + kl_divergence, + solve_tilt, + weighted_hdi, + weighted_quantile, +) +from impulso.scenario import MomentTarget, ProbabilityTarget + + +def _forecast_from_series(x: np.ndarray, n_chains: int = 1) -> np.ndarray: + """Wrap a 1-D sample as a `(C, D, 1, 1)` forecast array.""" + return x.reshape(n_chains, -1, 1, 1) + + +class TestVocabulary: + def test_horizon_must_be_positive(self): + with pytest.raises(ValueError, match="1-based"): + ProbabilityTarget(variable="y1", horizon=0, threshold=0.0, probability=0.5) + + def test_probability_bounds(self): + with pytest.raises(ValueError, match="0 < p <= 1"): + ProbabilityTarget(variable="y1", horizon=1, threshold=0.0, probability=0.0) + with pytest.raises(ValueError, match="0 < p <= 1"): + ProbabilityTarget(variable="y1", horizon=1, threshold=0.0, probability=1.5) + + def test_non_finite_threshold_rejected(self): + with pytest.raises(ValueError, match="finite"): + ProbabilityTarget(variable="y1", horizon=1, threshold=np.inf, probability=0.5) + + def test_moment_target_requires_finite_mean(self): + with pytest.raises(ValueError, match="finite"): + MomentTarget(variable="y1", horizon=1, mean=np.nan) + + def test_targets_are_frozen(self): + target = MomentTarget(variable="y1", horizon=2, mean=0.5) + with pytest.raises(ValueError, match="frozen"): + target.mean = 1.0 + + +class TestClosedForm: + """Single probability target: the two-mass solution and its oracles.""" + + @staticmethod + def _setup(n_total=100, n_event=20, probability=0.3): + x = np.arange(float(n_total)) + # The first n_event draws sit below the threshold. + forecast = _forecast_from_series(x) + target = ProbabilityTarget(variable="y1", horizon=1, threshold=float(n_event) - 0.5, probability=probability) + G, t = build_moments(forecast, [target], ["y1"], steps=1) + return G, t, n_total, n_event, probability + + def test_weights_match_the_analytic_two_mass_solution(self): + G, t, n_total, n_event, p = self._setup() + w, achieved = solve_tilt(G, t) + expected = np.where(G[:, 0] > 0.5, p / n_event, (1.0 - p) / (n_total - n_event)) + assert np.allclose(w, expected, atol=1e-15, rtol=0.0) + assert achieved[0] == pytest.approx(p, abs=1e-15) + assert w.sum() == pytest.approx(1.0, abs=1e-14) + + def test_ess_matches_the_analytic_value(self): + G, t, n_total, n_event, p = self._setup() + w, _ = solve_tilt(G, t) + oracle = 1.0 / (p**2 / n_event + (1.0 - p) ** 2 / (n_total - n_event)) + assert ess(w) == pytest.approx(oracle, rel=1e-13) + + def test_kl_matches_the_analytic_value(self): + G, t, n_total, n_event, p = self._setup() + w, _ = solve_tilt(G, t) + p_hat = n_event / n_total + oracle = p * np.log(p / p_hat) + (1.0 - p) * np.log((1.0 - p) / (1.0 - p_hat)) + assert kl_divergence(w) == pytest.approx(oracle, rel=1e-13) + + def test_probability_one_is_pure_conditioning(self): + G, t, _, n_event, _ = self._setup(probability=1.0) + w, achieved = solve_tilt(G, t) + in_event = G[:, 0] > 0.5 + assert np.allclose(w[in_event], 1.0 / n_event, atol=1e-15, rtol=0.0) + assert np.all(w[~in_event] == 0.0) + assert achieved[0] == pytest.approx(1.0, abs=1e-15) + assert ess(w) == pytest.approx(float(n_event), rel=1e-13) + + def test_direction_above_flips_the_event(self): + x = np.arange(100.0) + target = ProbabilityTarget(variable="y1", horizon=1, threshold=79.5, probability=0.5, direction="above") + G, _ = build_moments(_forecast_from_series(x), [target], ["y1"], steps=1) + assert G[:, 0].sum() == 20.0 + assert np.all(G[80:, 0] == 1.0) + + def test_uniform_weights_when_the_target_matches_the_sample(self): + G, t, _, _, _ = self._setup(probability=0.2) + w, _ = solve_tilt(G, t) + assert np.allclose(w, 1.0 / 100.0, atol=1e-15, rtol=0.0) + assert kl_divergence(w) == pytest.approx(0.0, abs=1e-15) + + +class TestGaussianMeanTilt: + """Tilting N(0,1) draws to a target mean — the RTW textbook case.""" + + @staticmethod + def _solve(mean=0.5, n_total=50_000): + x = np.random.default_rng(0).standard_normal(n_total) + target = MomentTarget(variable="x", horizon=1, mean=mean) + G, t = build_moments(_forecast_from_series(x), [target], ["x"], steps=1) + w, achieved = solve_tilt(G, t) + return x, w, achieved + + def test_achieved_mean_hits_the_target(self): + _, _, achieved = self._solve() + assert achieved[0] == pytest.approx(0.5, abs=1e-8) + + def test_log_weights_are_affine_in_the_variable(self): + x, w, _ = self._solve() + # log w = lambda x - log Z exactly, so a two-point fit reproduces + # every other draw. + design = np.column_stack([x, np.ones_like(x)]) + coef, *_ = np.linalg.lstsq(design, np.log(w), rcond=None) + residual = np.max(np.abs(design @ coef - np.log(w))) + assert residual < 1e-10 + + def test_dual_solution_matches_an_independent_root_find(self): + from scipy.optimize import brentq + + x, w, _ = self._solve() + design = np.column_stack([x, np.ones_like(x)]) + lam_hat = np.linalg.lstsq(design, np.log(w), rcond=None)[0][0] + + def gap(lam): + z = np.exp(lam * (x - x.max())) + return float((x * z).sum() / z.sum() - 0.5) + + lam_star = brentq(gap, -5.0, 5.0, xtol=1e-14, rtol=1e-15) + assert lam_hat == pytest.approx(lam_star, abs=1e-8) + + def test_dual_solution_is_near_the_population_value(self): + # Tilting a standard normal to mean m has population lambda = m. + x, w, _ = self._solve() + design = np.column_stack([x, np.ones_like(x)]) + lam_hat = np.linalg.lstsq(design, np.log(w), rcond=None)[0][0] + assert lam_hat == pytest.approx(0.5, abs=3.0 / np.sqrt(x.size)) + + +class TestMultiTargetDual: + def test_dual_matches_the_primal_slsqp_solution(self): + from scipy.optimize import minimize + + rng = np.random.default_rng(3) + n_total = 200 + y = rng.standard_normal((1, n_total, 2, 2)) + targets = [ + MomentTarget(variable="y1", horizon=1, mean=0.3), + ProbabilityTarget(variable="y2", horizon=2, threshold=0.0, probability=0.7), + ] + G, t = build_moments(y, targets, ["y1", "y2"], steps=2) + w, achieved = solve_tilt(G, t) + assert np.allclose(achieved, t, atol=1e-6) + + def entropy(v): + v = np.clip(v, 1e-300, None) + return float(np.sum(v * np.log(n_total * v))) + + primal = minimize( + entropy, + np.full(n_total, 1.0 / n_total), + method="SLSQP", + bounds=[(0.0, 1.0)] * n_total, + constraints=[ + {"type": "eq", "fun": lambda v: v.sum() - 1.0}, + {"type": "eq", "fun": lambda v, G=G, t=t: G.T @ v - t}, + ], + options={"maxiter": 500, "ftol": 1e-12}, + ) + assert np.allclose(w, primal.x, atol=1e-5) + + def test_multiple_probability_targets_take_the_dual_path(self): + rng = np.random.default_rng(11) + y = rng.standard_normal((2, 500, 3, 1)) + targets = [ + ProbabilityTarget(variable="y1", horizon=1, threshold=0.0, probability=0.8), + ProbabilityTarget(variable="y1", horizon=3, threshold=0.5, probability=0.6), + ] + G, t = build_moments(y, targets, ["y1"], steps=3) + w, achieved = solve_tilt(G, t) + assert np.allclose(achieved, t, atol=1e-8) + assert w.min() >= 0.0 + + +class TestInfeasibility: + def test_empty_event_names_the_draw_counts(self): + x = np.arange(100.0) + target = ProbabilityTarget(variable="y1", horizon=1, threshold=-5.0, probability=0.5) + with pytest.raises(ValueError, match=r"0 of 100 draws satisfy"): + build_moments(_forecast_from_series(x), [target], ["y1"], steps=1) + + def test_full_event_with_partial_probability_raises(self): + x = np.arange(100.0) + target = ProbabilityTarget(variable="y1", horizon=1, threshold=1e6, probability=0.5) + with pytest.raises(ValueError, match=r"only achievable probability is 1\.0"): + build_moments(_forecast_from_series(x), [target], ["y1"], steps=1) + + def test_moment_target_beyond_the_sample_range_raises(self): + x = np.arange(100.0) + target = MomentTarget(variable="y1", horizon=1, mean=500.0) + with pytest.raises(ValueError, match="outside the range spanned"): + build_moments(_forecast_from_series(x), [target], ["y1"], steps=1) + + def test_jointly_infeasible_targets_report_achieved_values(self): + # Two disjoint events cannot both carry probability 0.7. + x = np.arange(100.0) + forecast = _forecast_from_series(x) + targets = [ + ProbabilityTarget(variable="y1", horizon=1, threshold=49.5, probability=0.7), + ProbabilityTarget(variable="y1", horizon=1, threshold=49.5, probability=0.7, direction="above"), + ] + G, t = build_moments(forecast, targets, ["y1"], steps=1) + with pytest.raises(ValueError, match="not jointly achievable"): + solve_tilt(G, t) + + def test_unknown_variable_raises(self): + x = np.arange(10.0) + target = MomentTarget(variable="nope", horizon=1, mean=1.0) + with pytest.raises(ValueError, match="Unknown variable"): + build_moments(_forecast_from_series(x), [target], ["y1"], steps=1) + + def test_horizon_beyond_the_forecast_raises(self): + x = np.arange(10.0) + target = MomentTarget(variable="y1", horizon=4, mean=1.0) + with pytest.raises(ValueError, match="only 1 steps"): + build_moments(_forecast_from_series(x), [target], ["y1"], steps=1) + + def test_empty_target_list_raises(self): + with pytest.raises(ValueError, match="at least one target"): + build_moments(_forecast_from_series(np.arange(10.0)), [], ["y1"], steps=1) + + def test_wrong_target_type_raises(self): + from impulso.scenario import VariablePath + + with pytest.raises(TypeError, match="ProbabilityTarget or MomentTarget"): + build_moments( + _forecast_from_series(np.arange(10.0)), + [VariablePath(variable="y1", values=0.0)], + ["y1"], + steps=1, + ) + + +class TestWeightedSummaries: + def test_uniform_weights_reproduce_numpy_hazen_quantiles(self): + rng = np.random.default_rng(5) + x = rng.standard_normal(501) + w = np.full(x.size, 1.0 / x.size) + for q in (0.05, 0.25, 0.5, 0.75, 0.95): + assert weighted_quantile(x, w, q) == pytest.approx(float(np.quantile(x, q, method="hazen")), abs=1e-12) + + def test_two_mass_quantile_matches_the_closed_form(self): + x = np.array([0.0, 1.0]) + w = np.array([0.25, 0.75]) + # Knots at c = [0.125, 0.625]; q = 0.5 interpolates between them. + expected = (0.5 - 0.125) / (0.625 - 0.125) + assert weighted_quantile(x, w, 0.5) == pytest.approx(expected, abs=1e-14) + assert weighted_quantile(x, w, 0.05) == pytest.approx(0.0, abs=1e-14) + assert weighted_quantile(x, w, 0.99) == pytest.approx(1.0, abs=1e-14) + + def test_quantile_ignores_weight_normalisation_and_ordering(self): + rng = np.random.default_rng(6) + x = rng.standard_normal(200) + w = rng.random(200) + base = weighted_quantile(x, w / w.sum(), 0.4) + assert weighted_quantile(x, 7.0 * w, 0.4) == pytest.approx(base, abs=1e-12) + perm = rng.permutation(200) + assert weighted_quantile(x[perm], w[perm], 0.4) == pytest.approx(base, abs=1e-12) + + def test_quantile_is_monotone_in_q(self): + rng = np.random.default_rng(7) + x = rng.standard_normal(300) + w = rng.random(300) + levels = np.linspace(0.01, 0.99, 25) + values = [weighted_quantile(x, w, q) for q in levels] + assert np.all(np.diff(values) >= -1e-12) + + def test_quantile_vectorises_over_trailing_dimensions(self): + rng = np.random.default_rng(8) + x = rng.standard_normal((400, 3, 2)) + w = rng.random(400) + out = weighted_quantile(x, w, 0.5) + assert out.shape == (3, 2) + assert out[1, 1] == pytest.approx(weighted_quantile(x[:, 1, 1], w, 0.5), abs=1e-12) + + def test_uniform_weighted_hdi_matches_arviz(self): + import arviz as az + + rng = np.random.default_rng(9) + x = rng.standard_normal(501) # prob * n is not an integer + w = np.full(x.size, 1.0 / x.size) + lower, upper = weighted_hdi(x, w, prob=0.89) + expected = az.hdi(x, hdi_prob=0.89) + assert float(lower) == pytest.approx(float(expected[0]), abs=1e-12) + assert float(upper) == pytest.approx(float(expected[1]), abs=1e-12) + + def test_hdi_concentrates_when_the_weights_concentrate(self): + rng = np.random.default_rng(10) + x = rng.standard_normal(1000) + uniform = np.full(x.size, 1.0 / x.size) + concentrated = np.where(np.abs(x) < 0.25, 1.0, 1e-8) + concentrated = concentrated / concentrated.sum() + wide = weighted_hdi(x, uniform, prob=0.89) + narrow = weighted_hdi(x, concentrated, prob=0.89) + assert float(narrow[1] - narrow[0]) < float(wide[1] - wide[0]) + + def test_hdi_rejects_an_out_of_range_probability(self): + with pytest.raises(ValueError, match="0 < prob <= 1"): + weighted_hdi(np.arange(10.0), np.full(10, 0.1), prob=1.5) + + def test_summaries_reject_mismatched_draw_counts(self): + with pytest.raises(ValueError, match="draws on its leading axis"): + weighted_quantile(np.arange(10.0), np.full(5, 0.2), 0.5) + + +@pytest.fixture +def fitted_2v(synthetic_idata_2v, var_data_2v): + """Reduced-form 2-var VAR(1) from the synthetic posterior.""" + from impulso.fitted import FittedVAR + from impulso.volatility import Constant + + return FittedVAR( + idata=synthetic_idata_2v, + n_lags=1, + data=var_data_2v, + var_names=["y1", "y2"], + volatility=Constant(), + ) + + +def _median_target(forecast, variable="y1", horizon=3, probability=0.8): + """A probability target at the forecast's own empirical median.""" + da = forecast.idata.posterior_predictive["forecast"] + threshold = float(np.median(da.sel(variable=variable).isel(step=horizon - 1).values)) + return ProbabilityTarget(variable=variable, horizon=horizon, threshold=threshold, probability=probability) + + +class TestTiltEntryPoints: + def test_mean_mode_forecast_is_refused(self, fitted_2v): + forecast = fitted_2v.forecast(4, include_shock_uncertainty=False) + target = ProbabilityTarget(variable="y1", horizon=2, threshold=0.0, probability=0.5) + with pytest.raises(ValueError, match="needs a density forecast"): + forecast.tilt([target]) + + def test_mean_mode_conditional_forecast_is_refused(self, fitted_2v): + from impulso.scenario import VariablePath + + result = fitted_2v.conditional_forecast( + 4, + conditions=[VariablePath(variable="y1", values=np.array([0.1]))], + include_shock_uncertainty=False, + ) + with pytest.raises(ValueError, match="needs a density forecast"): + result.tilt([MomentTarget(variable="y2", horizon=2, mean=0.0)]) + + def test_unknown_variable_raises_through_tilt(self, fitted_2v): + forecast = fitted_2v.forecast(4, seed=0) + with pytest.raises(ValueError, match="Unknown variable"): + forecast.tilt([MomentTarget(variable="nope", horizon=1, mean=0.0)]) + + def test_horizon_beyond_the_forecast_raises_through_tilt(self, fitted_2v): + forecast = fitted_2v.forecast(4, seed=0) + with pytest.raises(ValueError, match="only 4 steps"): + forecast.tilt([MomentTarget(variable="y1", horizon=9, mean=0.0)]) + + def test_ess_warn_fraction_must_be_a_fraction(self, fitted_2v): + forecast = fitted_2v.forecast(4, seed=0) + with pytest.raises(ValueError, match=r"ess_warn_fraction must lie in \[0, 1\]"): + forecast.tilt([_median_target(forecast)], ess_warn_fraction=2.0) + + def test_scenario_result_inherits_tilt(self, fitted_2v): + from impulso.identification import Cholesky + + identified = fitted_2v.set_identification_strategy(Cholesky(ordering=["y1", "y2"])) + scenario = identified.structural_scenario(4, seed=2) + tilted = scenario.tilt([_median_target(scenario, horizon=2)]) + assert tilted.summary()["targets"][0]["achieved"] == pytest.approx(0.8, abs=1e-12) + + +class TestHardConditioningSurvivesTilting: + def test_pins_hold_on_every_draw_after_reweighting(self, fitted_2v): + from impulso.scenario import VariablePath + + pinned_path = np.array([0.2, 0.15, np.nan, np.nan]) + conditional = fitted_2v.conditional_forecast( + 4, + conditions=[VariablePath(variable="y1", values=pinned_path)], + seed=11, + ) + before = conditional.idata.posterior_predictive["forecast"].sel(variable="y1").values + assert np.allclose(before[:, :, 0], 0.2, atol=1e-10) + + tilted = conditional.tilt([_median_target(conditional, variable="y2", horizon=4, probability=0.75)]) + after = tilted.idata.posterior_predictive["forecast"].sel(variable="y1").values + # Reweighting never moves a draw, so the pins are untouched. + assert np.allclose(after[:, :, 0], 0.2, atol=1e-10) + assert np.allclose(after[:, :, 1], 0.15, atol=1e-10) + assert np.array_equal(before, after) + # The tilted median at a pinned step is the pin itself. + assert tilted.median()["y1"].iloc[0] == pytest.approx(0.2, abs=1e-10) + + def test_tilted_result_shares_the_parent_draws_by_reference(self, fitted_2v): + forecast = fitted_2v.forecast(4, seed=3) + tilted = forecast.tilt([_median_target(forecast)]) + parent = forecast.idata.posterior_predictive["forecast"].values + child = tilted.idata.posterior_predictive["forecast"].values + assert np.shares_memory(parent, child) + + +class TestDegeneracy: + def test_far_tail_target_warns_and_reports_a_small_ess(self, fitted_2v): + forecast = fitted_2v.forecast(4, seed=4) + da = forecast.idata.posterior_predictive["forecast"].sel(variable="y1").isel(step=3).values + # Keep only the most extreme few per cent of draws. + threshold = float(np.quantile(da, 0.03)) + target = ProbabilityTarget(variable="y1", horizon=4, threshold=threshold, probability=1.0) + with pytest.warns(UserWarning, match="tilt is concentrated"): + tilted = forecast.tilt([target]) + assert tilted.summary()["ess_fraction"] < 0.1 + + def test_no_warning_when_the_tilt_is_mild(self, fitted_2v): + import warnings + + forecast = fitted_2v.forecast(4, seed=4) + with warnings.catch_warnings(): + warnings.simplefilter("error") + forecast.tilt([_median_target(forecast, probability=0.55)]) + + +class TestTiltedResultSurface: + def test_median_and_hdi_mirror_the_forecast_result_shapes(self, fitted_2v): + forecast = fitted_2v.forecast(5, seed=5) + tilted = forecast.tilt([_median_target(forecast, horizon=2)]) + assert tilted.median().shape == forecast.median().shape + assert list(tilted.median().columns) == ["y1", "y2"] + hdi = tilted.hdi(prob=0.5) + assert hdi.lower.shape == (5, 2) + assert hdi.upper.shape == (5, 2) + assert np.all(hdi.upper.values >= hdi.lower.values) + assert tilted.to_dataframe().equals(tilted.median()) + + def test_uniform_tilt_reproduces_the_untilted_median(self, fitted_2v): + forecast = fitted_2v.forecast(4, seed=6) + da = forecast.idata.posterior_predictive["forecast"].sel(variable="y1").isel(step=1).values + n_below = int((da < 0.0).sum()) + # Requesting exactly the empirical probability leaves weights uniform. + target = ProbabilityTarget(variable="y1", horizon=2, threshold=0.0, probability=n_below / da.size) + tilted = forecast.tilt([target]) + assert tilted.summary()["kl_divergence"] == pytest.approx(0.0, abs=1e-12) + assert np.allclose(tilted.median().values, tilted.base_median().values, atol=1e-12) + + def test_summary_surfaces_requested_achieved_and_event_counts(self, fitted_2v): + forecast = fitted_2v.forecast(4, seed=7) + prob_target = _median_target(forecast, horizon=2, probability=0.65) + moment_target = MomentTarget( + variable="y2", + horizon=3, + mean=float(np.mean(forecast.idata.posterior_predictive["forecast"].sel(variable="y2").isel(step=2))), + ) + tilted = forecast.tilt([prob_target, moment_target]) + summary = tilted.summary() + assert set(summary) == {"ess", "ess_fraction", "kl_divergence", "n_draws", "targets"} + assert summary["n_draws"] == 100 + rows = summary["targets"] + assert rows[0]["requested"] == pytest.approx(0.65, abs=1e-12) + assert rows[0]["achieved"] == pytest.approx(0.65, abs=1e-8) + assert rows[0]["draws_in_event"] == 50 + assert rows[1]["draws_in_event"] is None + assert tilted.targets == [prob_target, moment_target] + + def test_weights_are_chain_draw_shaped_and_normalised(self, fitted_2v): + forecast = fitted_2v.forecast(4, seed=8) + tilted = forecast.tilt([_median_target(forecast)]) + assert tilted.weights.shape == (2, 50) + assert tilted.weights.sum() == pytest.approx(1.0, abs=1e-12) + + def test_plot_returns_a_figure(self, fitted_2v): + from matplotlib.figure import Figure + + forecast = fitted_2v.forecast(4, seed=9) + tilted = forecast.tilt([_median_target(forecast)]) + assert isinstance(tilted.plot(), Figure) + + def test_result_is_frozen(self, fitted_2v): + forecast = fitted_2v.forecast(4, seed=10) + tilted = forecast.tilt([_median_target(forecast)]) + with pytest.raises(ValueError, match="frozen"): + tilted.steps = 9