diff --git a/CONTEXT.md b/CONTEXT.md index f62ffad..0d97ac8 100644 --- a/CONTEXT.md +++ b/CONTEXT.md @@ -106,6 +106,25 @@ _Avoid_: bare "shrinkage" (ambiguous with cross-variable shrinkage). The conjugate VAR's closed-form log marginal likelihood of the observed response block, `log p(y_{p+1:T} | y_{1:p}, hyperparameters, model)`, attached to `FittedVAR.evidence` by `ConjugateVAR.fit` (`None` on the NUTS path, which has no closed form). Because the value includes the volatility-rescaling Jacobian it is a density over the *observed* data, so a break model and a homoscedastic model on the same observations are directly comparable. It is conditional on the presample and on the hyperparameters it was evaluated at, so with `NIWPrior(select=True)` a ratio of two evidences is an empirical-Bayes Bayes factor. `compare_evidence(**fits)` checks comparability (same variable set, effective sample, window and response digest) and returns an `EvidenceComparison` of Bayes factors and posterior model probabilities. _Avoid_: "log ML" in the API surface (spell out marginal likelihood); "model probability" for a raw Bayes factor — the probability requires prior model weights. +**Predictive pool (`PredictivePool`)**: +A set of fitted models plus one weight each, together defining a single combined predictive distribution. Weights are estimated by `pool_forecasts(fits, holdout)`, which forecasts every model from their shared estimation end and scores those forecasts on the held-out window. The pool is *static*: one weight per model, fixed across horizons and across time. `PredictivePool.combine()` applies the frozen weights to new forecasts (typically from full-sample refits), keeping estimation and application separate. +_Avoid_: "ensemble" (borrowed from ML, and suggests bagging/boosting rather than density combination); bare "model averaging" (ambiguous between BMA over a common likelihood and predictive pooling, which are different objects). + +**Held-out window**: +The stretch of data, supplied as its own `VARData`, on which candidate models' predictive densities are scored. It must genuinely postdate every model's estimation sample — checked against `FittedVAR.data.index[-1]`, which also fixes the shared **forecast origin**. Its length is the number of scored horizons `H`. +_Avoid_: "test set" (implies a train/test split of an exchangeable sample; time ordering is load-bearing here), "validation sample" (used inconsistently for in-sample tuning). + +**Forecast origin**: +The last timestamp of the shared estimation sample — the point every pooled forecast is made from. All models in a pool must share it, so the scores compare forecasts made from one information set. + +**Log-score weights**: +Weights proportional to `exp(total log score)`, i.e. a softmax over each model's summed held-out log predictive density (`method="log_score"`). Equivalent to pseudo-Bayesian model averaging on the held-out score, and they inherit its behaviour: the weight vector collapses onto the single best model as the held-out window lengthens, because total scores diverge linearly. +_Avoid_: "BMA weights" unqualified — genuine BMA uses marginal likelihoods over the estimation sample, not held-out predictive scores. + +**Stacking weights**: +Weights maximising the log score of the *pooled* predictive, `max_w Σ_h log Σ_i w_i p_i(y_h)` on the simplex (`method="stacking"`, the default). Unlike log-score weights they keep complementary models — two models that are each wrong in opposite directions can pool to a weight vector no single model matches — because the objective scores the mixture, not the members. Convex on the simplex, so the SLSQP optimum is global; matches `az.compare(method="stacking")`. +_Avoid_: "optimal weights" — optimal for the held-out log score only, and only for the scored window. + **Deterministic volatility break (`ConjugateVolatility`)**: A volatility process whose per-period scale `s_t` follows a deterministic, hyperparameter-driven path with a known break date — not a stochastic process. Used only by `ConjugateVAR`: the scale enters as data rescaling `ỹ_t = y_t / s_t` with a Jacobian in the marginal likelihood, and its hyperparameters are estimated jointly with λ. `PandemicBreak` (three outbreak scales + geometric decay from March 2020) is the concrete case reproducing Lenza & Primiceri (2020). _Avoid_: "stochastic volatility" — the break is deterministic given its hyperparameters. @@ -138,6 +157,7 @@ _Avoid_: "number of cointegrating vectors" in API surface (fine in prose); "coin - A **stationarity pretest** consumes `VARData` (endogenous block only), a DataFrame, or a Series, and produces a result object — never a modified dataset and never a specification. It sits *beside* the pipeline, not in it: nothing downstream of `VAR.fit()` reads its output. - **Integration order** feeds **cointegration rank**: the Johansen test is only meaningful for series that are individually integrated, and it is conditioned on a lag order (`k_ar_diff = p - 1`) that `select_lag_order` supplies. - A **ConjugateVAR** carries an **NIW prior** and optionally a **deterministic volatility break**; a **VAR** carries a **MinnesotaPrior** and a **PyMC volatility process** (`PyMCVolatilityProcess`, the `build_pymc_latent` extension of the `VolatilityProcess` query surface). Each estimator's fields accept only its compatible components, enforced by types + validators rather than a builder. +- Several **FittedVAR**s sharing a **forecast origin**, scored on a **held-out window**, produce a **predictive pool**; the pool's **stacking** or **log-score weights** then combine any matching set of forecasts. Estimation never refits, and the pool is indifferent to which estimator produced each `FittedVAR`. ## Example dialogue @@ -161,6 +181,8 @@ _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 have three candidate VARs. Which should I forecast with?" +> **Library:** Possibly all three. `pool_forecasts({"var2": f2, "var4": f4, "conj": fc}, holdout=held_out_data)` scores each model's density forecast on the held-out window and returns a `PredictivePool`; `pool.summary()` shows the weights and log scores. Refit on the full sample, then `pool.combine({...})` applies those weights to genuine forecasts. ## Conventions diff --git a/docs/adr/0006-static-predictive-pooling-on-draw-based-gaussian-densities.md b/docs/adr/0006-static-predictive-pooling-on-draw-based-gaussian-densities.md new file mode 100644 index 0000000..0f5cad4 --- /dev/null +++ b/docs/adr/0006-static-predictive-pooling-on-draw-based-gaussian-densities.md @@ -0,0 +1,58 @@ +# Static predictive pooling scores draw-based Gaussian densities on a held-out window + +Combining several fitted VARs into one predictive distribution needs a score. Impulso scores each +candidate model's predictive density on an explicit **held-out window** and turns those scores into +**static weights** — one weight per model, fixed across horizons and across time. `pool_forecasts` +owns estimation (it calls `.forecast()` itself, so the forecasts it scores are provably +density-mode, from a single origin, on the models' own training samples); `PredictivePool.combine()` +owns application (frozen weights applied to *new* forecasts from full-sample refits). Weights come +from either **stacking** (default; the log-score of the pooled predictive, maximised over the +simplex) or **log-score/pseudo-BMA** weights (softmax of per-model total log scores). + +## Considered options + +- **Exact Rao-Blackwellised mixture score** — for each held-out point, average the per-draw Gaussian + predictive densities implied by each posterior draw rather than moment-matching the draws to one + Gaussian. Strictly better (it keeps the mixture's heavy tails) and deferred, not rejected: the + `density=` keyword is the seam, and `density="mixture"` is where it lands. The moment-matched + Gaussian was chosen first because it needs only the forecast draws — no re-entry into the + volatility process — and therefore works identically for every estimator (`VAR`, `ConjugateVAR`) + and every volatility process. +- **Per-variable (diagonal) densities only** — rejected as the default: it throws away the + cross-variable correlation that is the entire point of a VAR. Retained as the `density="diagonal"` + escape hatch for near-singular covariances and small draw counts. +- **Joint-path scoring** — score the whole `H`-step path as one `H·n`-dimensional Gaussian instead + of summing per-horizon scores. Rejected: it needs `S > H·n` draws to be non-singular at realistic + horizons, and the per-horizon sum is the standard density-forecast-evaluation object. +- **Rolling-origin scoring** — re-forecast from each held-out date. Rejected for v1: it multiplies + cost by `H` and requires refitting to stay honest about information sets. One fixed origin keeps + the contract provable from `FittedVAR.data.index[-1]`. +- **Dynamic (time-varying) weights** — Del Negro, Hasegawa & Schorfheide (2016). Out of scope; the + static pool is the object this ADR fixes, and dynamic weights slot in behind the existing + `method=` keyword when they arrive. +- **Pooling over `ForecastResult` objects supplied by the user** — rejected as the estimation entry + point: a `ForecastResult` carries no origin or estimation metadata, so "these were produced from + the same information set" is unenforceable. `pool_forecasts` takes `FittedVAR`s and forecasts them + itself. `combine()` *does* take forecasts, because by then the weights are already fixed. + +## Consequences + +- The **score matrix** — `(H, M)` log predictive densities, one row per held-out date, one column + per model — is the contract between scoring and weighting. Both solvers consume it and nothing + else, which is why they are unit-testable against closed-form optima and why a future + `density="mixture"` is a drop-in. +- Scores are *comparable across models*, not exact predictive densities: the Gaussian + moment-matching approximates a genuinely heavier-tailed posterior-predictive mixture. The + approximation degrades with few draws, stochastic volatility, and heavy tails. +- The pool **never refits**. It cannot, in general — the held-out window postdates every model's + estimation sample by construction, and refitting on it would destroy the held-out property. Users + who want full-sample models for genuine forecasting refit themselves and pass those forecasts to + `combine()`. +- Exponentiation is done after a **per-row maximum shift** (stacking) or a global maximum shift + (log-score weights). The row shift changes the stacking objective by a `w`-independent constant, + so the optimum is unchanged, and it keeps the pool finite at log scores where an unshifted + implementation underflows to zero and reports a degenerate weight vector. +- ArviZ parity is a *property Impulso tests*, not a dependency it takes: with degenerate + (draw-constant) log-likelihoods, `pool_forecasts` reproduces `az.compare(method="stacking")` and + `method="pseudo-BMA"`. The primary regression check is a direct SciPy re-solve of the score + matrix, so ArviZ API drift cannot silently weaken the suite. diff --git a/docs/how-to/index.md b/docs/how-to/index.md index 7ff315b..cb02f1d 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 +pooling ``` diff --git a/docs/how-to/pooling.md b/docs/how-to/pooling.md new file mode 100644 index 0000000..ae5d75f --- /dev/null +++ b/docs/how-to/pooling.md @@ -0,0 +1,172 @@ +# Pooling Several Models + +When two lag orders, two priors, or two estimators all look defensible, you do +not have to pick one. Score them on data they were not fitted to, and let the +scores decide how much of each to keep. + +## Split the sample + +The held-out window must be the periods immediately following the estimation +sample. Split once, keep the split fixed, and fit every candidate on the same +training data: + +```python +import pandas as pd +from impulso import VAR, ConjugateVAR, VARData + +df = pd.read_csv("macro_data.csv", index_col="date", parse_dates=True) +train = VARData.from_df(df.iloc[:-12], endog=["gdp", "inflation", "rate"]) +holdout = VARData.from_df(df.iloc[-12:], endog=["gdp", "inflation", "rate"]) +``` + +## Fit the candidates + +Any mix of estimators works — the pool only needs `FittedVAR` objects that share +their variables and their sample end: + +```python +fits = { + "var2": VAR(lags=2, prior="minnesota").fit(train), + "var4": VAR(lags=4, prior="minnesota").fit(train), + "conjugate": ConjugateVAR(lags=4).fit(train), +} +``` + +## Estimate the weights + +```python +from impulso import pool_forecasts + +pool = pool_forecasts(fits, holdout, method="stacking", seed=0) +pool.summary() +``` + +`summary()` returns one row per model, heaviest weight first, with the total and +per-period held-out log score — the shape of the output, with your own data's +numbers in it: + +| | weight | log_score | mean_log_score | rank | +|---|---|---|---|---| +| var4 | 0.62 | ... | ... | 1 | +| var2 | 0.38 | ... | ... | 2 | +| conjugate | 0.00 | ... | ... | 3 | + +`pool.log_scores` is the full matrix behind those totals — one row per held-out +date, one column per model — and `pool.to_dataframe()` adds the pooled +predictive's own score alongside them. `pool.plot()` draws the weights as a +ranked bar chart. + +## Stacking versus log-score weights + +The two methods answer different questions. When one candidate simply forecasts +better than the others, they agree and both hand it the weight. Take a tightly +shrunk VAR(1) against a loosely shrunk VAR(4), fitted on the same synthetic +series and scored over sixteen held-out quarters: + +```python +import numpy as np +import pandas as pd +from impulso import VARData, pool_forecasts +from impulso.conjugate import ConjugateVAR +from impulso.priors import NIWPrior + +rng = np.random.default_rng(0) +T, H = 240, 16 +A = np.array([[0.55, 0.15], [-0.20, 0.45]]) +y = np.zeros((T + H, 2)) +for t in range(1, T + H): + calm = (t // 4) % 2 == 0 + y[t] = A @ y[t - 1] + rng.standard_normal(2) * ( + np.array([0.3, 1.6]) if calm else np.array([1.6, 0.3]) + ) +index = pd.date_range("1980-01-01", periods=T + H, freq="QS") +frame = pd.DataFrame(y, columns=["output", "prices"], index=index) + +train = VARData.from_df(frame.iloc[:T], endog=["output", "prices"]) +held_out = VARData.from_df(frame.iloc[T:], endog=["output", "prices"]) +candidates = { + "var1_tight": ConjugateVAR(lags=1, prior=NIWPrior(tightness=0.05), draws=500, seed=0).fit(train), + "var4_loose": ConjugateVAR(lags=4, prior=NIWPrior(tightness=1.0), draws=500, seed=1).fit(train), +} + +stacked = pool_forecasts(candidates, held_out, method="stacking", seed=0) +log_scored = pool_forecasts(candidates, held_out, method="log_score", seed=0) +``` + +The two models score -56.7 and -51.2 over the window, and both rules reach the +same verdict: + +| | `var1_tight` | `var4_loose` | pooled log score | +|---|---|---|---| +| `method="stacking"` | 0.000 | 1.000 | -51.2 | +| `method="log_score"` | 0.004 | 0.996 | -51.2 | + +This is the common case: one model dominates, the pool finds it, and the pooled +score matches the winner's. Pooling has told you something useful — that the +second candidate adds nothing — even though it did not improve the forecast. + +The rules come apart when the candidates are genuinely **complementary**: each +one better over a different part of the held-out window, neither better +throughout. Stacking scores the *mixture*, so it can hold both models at +non-trivial weight and reach a score above anything either achieves alone — a +mixture covers regions that no single member covers. Log-score weights compare +the models one at a time and take a softmax of their totals, so they concentrate +on whichever model has the better total even when a blend would score higher. +Whether the resulting pool beats the best single model depends on how close the +totals are; the point is that log-score weights are not trying to maximise the +pooled score, and stacking is. + +Log-score weights also concentrate harder the longer the holdout, because total +scores diverge roughly linearly in the number of held-out periods. That is a +property of the rule, not a bug: it is the behaviour you want if you believe one +candidate is correct and want the evidence to say which. Reach for stacking when +you want the best combined forecast instead. + +## Forecast with the weights + +`pool.holdout_predictive` covers the *held-out window* — useful for plotting the +combination against what actually happened, useless as a forecast, since those +dates have already occurred. For a real forecast, refit on the full sample and +apply the frozen weights: + +```python +full = VARData.from_df(df, endog=["gdp", "inflation", "rate"]) +refits = { + "var2": VAR(lags=2, prior="minnesota").fit(full), + "var4": VAR(lags=4, prior="minnesota").fit(full), + "conjugate": ConjugateVAR(lags=4).fit(full), +} +forecasts = {label: fit.forecast(steps=8) for label, fit in refits.items()} + +combined = pool.combine(forecasts, seed=1) +combined.median() +combined.hdi(0.89) +``` + +`combine` needs the same model labels and variables as the pool, and all +forecasts must run to the same horizon — but that horizon need not be the one +the weights were scored over. It also insists on density-mode forecasts +(`include_shock_uncertainty=True`, the default), because a mean forecast has no +predictive density to pool. + +## What the scores are and are not + +- Each horizon's predictive density is a **Gaussian matched to the forecast + draws**, joint across variables. The true posterior predictive is a + heavier-tailed mixture over draws, so the scores are comparable across models + rather than exact. The approximation degrades with few draws, stochastic + volatility, and fat tails. Pass `density="diagonal"` to score each variable on + its own marginal when the joint covariance is near-singular. +- Scores are **summed over horizons 1 to H from one fixed origin**. This is not + a joint-path density and not a rolling-origin evaluation. +- Weights are **static** — one per model, fixed across horizons and time. A + model that only wins at long horizons cannot be given a horizon-specific + weight. +- A short holdout makes the weights noisy. Twelve quarters is thin; two + observations tell you almost nothing. +- Impulso checks that the holdout postdates every model's estimation sample. It + cannot check that you did not look at the holdout while choosing the + candidates — that part is on you. +- Nothing records how a series was transformed, so pooling models fitted on + differently transformed data will produce meaningless scores. Keep the + transformations identical across candidates. diff --git a/docs/index.md b/docs/index.md index 000a111..3caa8f4 100644 --- a/docs/index.md +++ b/docs/index.md @@ -41,6 +41,7 @@ irf.plot() - **PyMC backend** — full Bayesian estimation with NUTS sampling - **Probabilistic forecasts** — posterior median, HDI credible intervals, tidy DataFrames - **Structural identification** — Cholesky and sign restriction schemes +- **Predictive pooling** — weight competing models by held-out log score, via stacking or pseudo-BMA - **Built-in plotting** — IRF, FEVD, forecast, and historical decomposition plots ## Installation diff --git a/docs/reference/index.md b/docs/reference/index.md index 100e069..383e311 100644 --- a/docs/reference/index.md +++ b/docs/reference/index.md @@ -18,6 +18,7 @@ fitted identified identification scenario +pooling results evidence primitives diff --git a/docs/reference/plotting.md b/docs/reference/plotting.md index ebe6148..6ed8208 100644 --- a/docs/reference/plotting.md +++ b/docs/reference/plotting.md @@ -15,6 +15,7 @@ plot_fevd plot_historical_decomposition plot_counterfactual + plot_pool_weights plot_volatility plot_sv_forecast ``` diff --git a/docs/reference/pooling.md b/docs/reference/pooling.md new file mode 100644 index 0000000..1cdf502 --- /dev/null +++ b/docs/reference/pooling.md @@ -0,0 +1,24 @@ +# Predictive Pooling + +Combine several fitted models into one predictive distribution. `pool_forecasts` +forecasts every candidate from their shared forecast origin, scores those +densities against a held-out window, and turns the resulting log-score matrix +into weights — by stacking (the log score of the *pooled* predictive, maximised +over the simplex) or by a softmax of each model's total score. + +The returned `PredictivePool` carries the weights, the full score matrix, and a +pooled predictive sample over the held-out window. Its weights are frozen once +estimated: `PredictivePool.combine` applies them to new forecasts, which is how +you get a genuine out-of-sample forecast from full-sample refits without +rescoring anything. + +```{eval-rst} +.. currentmodule:: impulso.pooling + +.. autosummary:: + :toctree: generated/ + :nosignatures: + + pool_forecasts + PredictivePool +``` diff --git a/src/impulso/__init__.py b/src/impulso/__init__.py index 8fe3f84..a3d0707 100644 --- a/src/impulso/__init__.py +++ b/src/impulso/__init__.py @@ -19,6 +19,7 @@ from impulso.identification import Cholesky, LongRunRestriction, ProxySVAR, SignRestriction from impulso.identified import IdentifiedVAR from impulso.observation import Gaussian, StudentT + from impulso.pooling import PredictivePool, pool_forecasts from impulso.priors import MinnesotaPrior, NIWPrior from impulso.protocols import ErrorDistribution, VolatilityProcess from impulso.results import ( @@ -75,6 +76,7 @@ "NIWPrior", "NUTSSampler", "PandemicBreak", + "PredictivePool", "ProxySVAR", "SVData", "SVDefaultPrior", @@ -97,6 +99,7 @@ "johansen_test", "kpss_test", "lag_matrices", + "pool_forecasts", "select_lag_order", ] @@ -145,6 +148,8 @@ "VolatilityProcess": "impulso.protocols", "compute_ma_phi": "impulso._ma", "lag_matrices": "impulso._linalg", + "PredictivePool": "impulso.pooling", + "pool_forecasts": "impulso.pooling", } """Map of lazily-exported name to the module that defines it. diff --git a/src/impulso/plotting/__init__.py b/src/impulso/plotting/__init__.py index 94a8474..c4d390c 100644 --- a/src/impulso/plotting/__init__.py +++ b/src/impulso/plotting/__init__.py @@ -7,6 +7,7 @@ from impulso.plotting._forecast import plot_forecast from impulso.plotting._historical_decomposition import plot_historical_decomposition from impulso.plotting._irf import plot_irf +from impulso.plotting._pooling import plot_pool_weights 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 @@ -19,6 +20,7 @@ "plot_forecast", "plot_historical_decomposition", "plot_irf", + "plot_pool_weights", "plot_structural_scenario", "plot_sv_forecast", "plot_volatility", diff --git a/src/impulso/plotting/_pooling.py b/src/impulso/plotting/_pooling.py new file mode 100644 index 0000000..8cd207e --- /dev/null +++ b/src/impulso/plotting/_pooling.py @@ -0,0 +1,55 @@ +"""Predictive-pool plotting.""" + +from typing import TYPE_CHECKING + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.figure import Figure + +if TYPE_CHECKING: + from impulso.pooling import PredictivePool + + +def plot_pool_weights( + pool: "PredictivePool", + figsize: tuple[float, float] | None = None, +) -> Figure: + """Plot pool weights as a ranked horizontal bar chart. + + Bars run heaviest weight first and are annotated with each model's total + held-out log score, so a model carrying little weight can be read against + how badly it actually scored. The title reports the weighting method and + the pooled score alongside the best single model's, which is the number + that says whether pooling bought anything. + + Args: + pool: PredictivePool from `pool_forecasts`. + figsize: Figure size. Defaults to `(8, 1 + 0.5 * n_models)`. + + Returns: + Matplotlib Figure. + """ + summary = pool.summary() + labels = list(summary.index) + weights = summary["weight"].to_numpy(dtype=float) + scores = summary["log_score"].to_numpy(dtype=float) + + fig, ax = plt.subplots(figsize=figsize or (8.0, 1.0 + 0.5 * len(labels))) + positions = np.arange(len(labels))[::-1] + ax.barh(positions, weights, color="tab:blue", alpha=0.8) + ax.set_yticks(positions) + ax.set_yticklabels(labels) + ax.set_xlabel("weight") + ax.set_xlim(0.0, max(1.0, float(weights.max()) * 1.05)) + + offset = 0.01 * ax.get_xlim()[1] + for position, weight, score in zip(positions, weights, scores, strict=True): + ax.text(weight + offset, position, f"log score {score:,.1f}", va="center", fontsize=8) + + pooled = pool.pooled_log_score() + ax.set_title( + f"Pool weights ({pool.method}, {pool.density} density) — " + f"pooled log score {pooled:,.1f} vs best single {scores.max():,.1f}" + ) + fig.tight_layout() + return fig diff --git a/src/impulso/pooling.py b/src/impulso/pooling.py new file mode 100644 index 0000000..b8cac0e --- /dev/null +++ b/src/impulso/pooling.py @@ -0,0 +1,755 @@ +"""Static predictive-density pooling. + +Score several fitted models' predictive densities on a held-out window, turn +those scores into weights, and combine the models into a single predictive +distribution. `pool_forecasts` estimates the weights; `PredictivePool.combine` +applies them to new forecasts. See ADR-0006 for the design and its limits. +""" + +from __future__ import annotations + +import warnings +from typing import TYPE_CHECKING, Literal + +import arviz as az +import numpy as np +import pandas as pd +from matplotlib.figure import Figure +from pydantic import Field, model_validator + +from impulso._base import ImpulsoBaseModel +from impulso.data import VARData +from impulso.results import ForecastResult + +if TYPE_CHECKING: + from collections.abc import Mapping + + from impulso.fitted import FittedVAR + +WeightMethod = Literal["stacking", "log_score"] +"""How held-out log scores become weights.""" + +DensityKind = Literal["gaussian", "diagonal"] +"""How forecast draws become a predictive density.""" + +_COV_JITTER = 1e-10 +_LOG_2PI = float(np.log(2.0 * np.pi)) +_NO_OPTIMISER = "log-score weights are closed-form; no optimiser was run." +_POOLED_COLUMN = "pooled" + + +# -------------------------------------------------------------------------- +# Predictive densities: forecast draws -> (H,) log scores +# -------------------------------------------------------------------------- + + +def _degenerate_density_message(label: str, h: int, reason: str) -> str: + return ( + f"Model {label!r} has a degenerate joint predictive density at holdout step {h + 1}: " + f"{reason}. Pass density='diagonal' to score each variable on its own marginal, " + "or draw more posterior samples." + ) + + +def _joint_log_scores(draws: np.ndarray, y: np.ndarray, label: str) -> np.ndarray: + """Moment-matched joint Gaussian log density, one score per horizon.""" + from scipy.linalg import solve_triangular + + _, n_steps, n_vars = draws.shape + scores = np.empty(n_steps) + eye = np.eye(n_vars) + for h in range(n_steps): + block = draws[:, h, :] + mean = block.mean(axis=0) + cov = np.cov(block, rowvar=False, ddof=1) + mean_var = float(np.mean(np.diag(cov))) + if mean_var <= 0.0: + raise ValueError(_degenerate_density_message(label, h, "every forecast draw is identical")) + # The sample covariance is positive semi-definite by construction, so + # proving the smallest eigenvalue clears the jitter both rejects + # singular densities (where the jitter, not the data, would be doing + # the work) and bounds the jittered condition number at ~n/_COV_JITTER, + # which Cholesky handles comfortably in float64. + if float(np.linalg.eigvalsh(cov)[0]) <= _COV_JITTER * mean_var: + raise ValueError(_degenerate_density_message(label, h, "the predictive covariance is numerically singular")) + chol = np.linalg.cholesky(cov + _COV_JITTER * mean_var * eye) + z = solve_triangular(chol, y[h] - mean, lower=True) + scores[h] = -0.5 * (n_vars * _LOG_2PI + 2.0 * float(np.sum(np.log(np.diag(chol)))) + float(z @ z)) + return scores + + +def _diagonal_log_scores(draws: np.ndarray, y: np.ndarray, label: str) -> np.ndarray: + """Per-variable normal log density summed over variables, one score per horizon.""" + n_vars = draws.shape[2] + mean = draws.mean(axis=0) + var = draws.var(axis=0, ddof=1) + if not bool((var > 0.0).all()): + h, v = (int(i) for i in np.argwhere(var <= 0.0)[0]) + raise ValueError( + f"Model {label!r} has zero forecast variance for variable index {v} at holdout step {h + 1}, " + "so its marginal predictive density is degenerate. Draw more posterior samples, or forecast " + "in density mode (include_shock_uncertainty=True)." + ) + z = (y - mean) / np.sqrt(var) + return -0.5 * (n_vars * _LOG_2PI + np.sum(np.log(var), axis=1) + np.sum(z**2, axis=1)) + + +def _gaussian_log_scores( + draws: np.ndarray, + y: np.ndarray, + density: DensityKind = "gaussian", + label: str = "", +) -> np.ndarray: + """Log predictive density of each held-out row under a draw-based Gaussian. + + Args: + draws: Flattened forecast draws of shape `(S, H, n)` — posterior + draws by horizon by variable. + y: Held-out realisations of shape `(H, n)`. + density: `"gaussian"` for a joint density across variables (moments + matched to the draws), `"diagonal"` for the sum of per-variable + marginals. + label: Model label, used only in error messages. + + Returns: + Array of `H` log scores. + """ + n_sims, n_steps, n_vars = draws.shape + if y.shape != (n_steps, n_vars): + raise ValueError(f"Held-out array has shape {y.shape} but the draws imply {(n_steps, n_vars)}.") + if not bool(np.isfinite(draws).all()): + raise ValueError( + f"Model {label!r} produced non-finite forecast draws (NaN or Inf); its predictive density " + "cannot be scored. This usually means an explosive posterior draw — check the fit." + ) + if n_sims < 2: + raise ValueError(f"Model {label!r} has {n_sims} posterior draw(s); scoring needs at least two.") + if density == "diagonal": + return _diagonal_log_scores(draws, y, label) + if n_sims <= n_vars: + raise ValueError( + f"Model {label!r} has {n_sims} posterior draws for {n_vars} variables; with S <= n the joint " + "predictive covariance is rank-deficient. Sample more draws, or pass density='diagonal'." + ) + return _joint_log_scores(draws, y, label) + + +# -------------------------------------------------------------------------- +# Weight solvers: (H, M) log-score matrix -> (M,) weights +# -------------------------------------------------------------------------- + + +def _check_score_matrix(log_scores: np.ndarray, index: pd.Index | None = None) -> None: + """Reject score matrices no weight vector can be fitted to.""" + if log_scores.ndim != 2: + raise ValueError(f"The log-score matrix must be 2-D (horizons by models), got {log_scores.ndim}-D.") + if log_scores.shape[1] < 2: + raise ValueError(f"Pooling requires at least two fitted models, got {log_scores.shape[1]}.") + if bool(np.isnan(log_scores).any()) or bool(np.isposinf(log_scores).any()): + raise ValueError("The log-score matrix contains NaN or +inf entries; log densities must be finite or -inf.") + dead = ~np.isfinite(log_scores).any(axis=1) + if bool(dead.any()): + pos = int(np.flatnonzero(dead)[0]) + where = f"{index[pos]}" if index is not None else f"holdout step {pos + 1}" + raise ValueError( + f"No model assigns positive predictive density at {where}, so the pooled score is -inf for every " + "weight vector. Check that point for an outlier or level shift, or widen the predictive densities " + "(more posterior draws, or density='diagonal')." + ) + + +def _log_score_weights(log_scores: np.ndarray, index: pd.Index | None = None) -> np.ndarray: + """Softmax of each model's total held-out log score (pseudo-BMA weights). + + Shifting by the maximum total before exponentiating keeps the weights + finite at log scores where an unshifted implementation overflows. + """ + _check_score_matrix(log_scores, index) + totals = log_scores.sum(axis=0) + if not bool(np.isfinite(totals).any()): + raise ValueError( + "Every model scores -inf somewhere in the holdout, so every total log score is -inf and " + "log-score weights are undefined. Use method='stacking', which scores the pooled density " + "row by row rather than model by model." + ) + weights = np.exp(totals - totals.max()) + return weights / weights.sum() + + +def _stacking_weights(log_scores: np.ndarray, index: pd.Index | None = None) -> tuple[np.ndarray, bool, str]: + """Weights maximising the log score of the pooled predictive. + + The objective is convex on the simplex, so the SLSQP solution is the + global optimum. Densities are exponentiated after a per-row maximum + shift, which changes the objective by a `w`-independent constant and + therefore leaves the optimum untouched while keeping the pool finite at + log scores that would otherwise underflow to zero. + + Returns: + Tuple of `(weights, converged, optimiser_message)`. + """ + from scipy.optimize import Bounds, LinearConstraint, minimize + + _check_score_matrix(log_scores, index) + n_models = log_scores.shape[1] + densities = np.exp(log_scores - log_scores.max(axis=1, keepdims=True)) + + def objective(w: np.ndarray) -> float: + return -float(np.sum(np.log(densities @ w))) + + def gradient(w: np.ndarray) -> np.ndarray: + return -np.sum(densities / (densities @ w)[:, None], axis=0) + + solver_kwargs = { + "fun": objective, + "jac": gradient, + "method": "SLSQP", + "bounds": Bounds(0.0, 1.0), + "constraints": LinearConstraint(np.ones((1, n_models)), 1.0, 1.0), + "options": {"ftol": 1e-12, "maxiter": 1000}, + } + with np.errstate(divide="ignore", invalid="ignore"): + result = minimize(x0=np.full(n_models, 1.0 / n_models), **solver_kwargs) + if not result.success: + # Restart from a data-driven point on the simplex; unlike the + # log-score weights this is finite even when every model scores + # -inf somewhere. + column_mass = densities.sum(axis=0) + result = minimize(x0=column_mass / column_mass.sum(), **solver_kwargs) + if not result.success: + raise RuntimeError( + f"Stacking weights failed to converge: {result.message}. The objective is convex on the " + "simplex, so this points at a pathological score matrix — inspect the log scores, or fall " + "back to method='log_score'." + ) + weights = np.clip(np.asarray(result.x, dtype=float), 0.0, None) + return weights / weights.sum(), bool(result.success), str(result.message) + + +def _pooled_row_scores(log_scores: np.ndarray, weights: np.ndarray) -> np.ndarray: + """Log score of the pooled predictive at each held-out point.""" + from scipy.special import logsumexp + + with np.errstate(divide="ignore", invalid="ignore"): + return np.asarray(logsumexp(log_scores + np.log(weights), axis=1), dtype=float) + + +# -------------------------------------------------------------------------- +# Input validation +# -------------------------------------------------------------------------- + + +def _resolve_rng(seed: int | np.random.Generator | None) -> np.random.Generator: + if isinstance(seed, np.random.Generator): + return seed + try: + return np.random.default_rng(seed) + except (TypeError, ValueError) as exc: + raise ValueError(f"seed must be an int, None, or a numpy.random.Generator, got {type(seed).__name__}.") from exc + + +def _spawn(rng: np.random.Generator, count: int) -> list[np.random.Generator]: + try: + return list(rng.spawn(count)) + except (AttributeError, TypeError) as exc: + raise ValueError( + "seed must be an int, None, or a numpy.random.Generator that supports spawning " + "(one built by numpy.random.default_rng does)." + ) from exc + + +def _index_freq(index: pd.DatetimeIndex) -> pd.offsets.BaseOffset | None: + """Best-effort frequency for a DatetimeIndex, or None.""" + if index.freq is not None: + return index.freq + if len(index) >= 3: + inferred = pd.infer_freq(index) + if inferred is not None: + return pd.tseries.frequencies.to_offset(inferred) + return None + + +def _check_reserved_labels(labels: list[str]) -> None: + """Reject model labels that would collide with a generated column.""" + if _POOLED_COLUMN in labels: + raise ValueError( + f"{_POOLED_COLUMN!r} is a reserved model label: PredictivePool.to_dataframe() adds a " + f"{_POOLED_COLUMN!r} column holding the combined predictive's score, which would silently " + "overwrite that model's own scores. Rename the model." + ) + + +def _check_variables(label: str, fit_names: list[str], holdout_names: list[str]) -> None: + if fit_names == holdout_names: + return + message = ( + f"Model {label!r} was fitted on variables {fit_names} but the holdout carries " + f"{holdout_names}; pooled models must share the holdout's variables" + ) + if set(fit_names) == set(holdout_names): + raise ValueError(f"{message} in the same order — reorder the holdout columns to {fit_names}.") + raise ValueError(f"{message}.") + + +def _check_exog(label: str, fit: FittedVAR, holdout: VARData) -> None: + if not fit.has_exog: + return + if holdout.exog is None: + raise ValueError( + f"Model {label!r} was fitted with exogenous regressors, so the holdout must carry their " + "future values too; rebuild it with VARData(..., exog=..., exog_names=...)." + ) + fit_names = list(fit.data.exog_names or []) + holdout_names = list(holdout.exog_names or []) + if fit_names != holdout_names: + raise ValueError( + f"Model {label!r} was fitted with exogenous regressors named {fit_names} but the holdout " + f"carries {holdout_names}; they must match in name and order." + ) + + +def _check_alignment(fits: Mapping[str, FittedVAR], holdout: VARData, origin: pd.Timestamp) -> None: + """Require one shared frequency across the models that the holdout continues. + + Checked against *every* model's estimation index, not just the first: two + models can share a sample end while running at different frequencies, in + which case a forecast step means a different span of time for each and + their held-out scores are not comparable. + """ + freqs = {label: _index_freq(fit.data.index) for label, fit in fits.items()} + known = {label: freq for label, freq in freqs.items() if freq is not None} + if len({freq.freqstr for freq in known.values()}) > 1: + stamps = ", ".join(f"{label}={freq.freqstr}" for label, freq in known.items()) + raise ValueError( + f"Pooled models were estimated at different frequencies: {stamps}. They share a sample " + "end, but a forecast step spans a different amount of time for each, so their held-out " + "scores are not comparable. Resample the candidates onto a common frequency before pooling." + ) + + freq = next(iter(known.values()), None) + if freq is None: + freq = _index_freq(holdout.index) + if freq is None: + warnings.warn( + "Could not infer a frequency for the estimation samples or the holdout, so the held-out " + "dates are assumed to line up positionally with forecast steps 1..H. Pass data with a " + "regular DatetimeIndex if you want that checked.", + UserWarning, + stacklevel=4, + ) + return + if len(known) < len(fits): + unknown = sorted(set(freqs) - set(known)) + warnings.warn( + f"Could not infer a frequency for {unknown}, so their forecast steps are assumed to line " + f"up positionally with the {freq.freqstr} held-out dates.", + UserWarning, + stacklevel=4, + ) + expected = pd.date_range(origin, periods=len(holdout.index) + 1, freq=freq)[1:] + mismatch = np.flatnonzero(expected.to_numpy() != holdout.index.to_numpy()) + if mismatch.size: + first = int(mismatch[0]) + raise ValueError( + f"The holdout does not continue the estimation sample at step {first + 1}: expected " + f"{expected[first].date()} at frequency {freq.freqstr}, got {holdout.index[first].date()}. " + "Pooled scores are indexed by forecast step, so the holdout must be the H periods " + "immediately after the forecast origin." + ) + + +def _validate_pool_inputs(fits: Mapping[str, FittedVAR], holdout: VARData) -> pd.Timestamp: + """Check a pool's inputs and return the shared forecast origin.""" + if len(fits) < 2: + raise ValueError(f"Pooling requires at least two fitted models, got {len(fits)}.") + if len(holdout.index) < 1: + raise ValueError("Pooling needs at least one held-out observation to score; the holdout is empty.") + _check_reserved_labels(list(fits)) + + ends = {} + for label, fit in fits.items(): + _check_variables(label, list(fit.var_names), list(holdout.endog_names)) + _check_exog(label, fit, holdout) + ends[label] = fit.data.index[-1] + if len(set(ends.values())) > 1: + stamps = ", ".join(f"{label}={end.date()}" for label, end in ends.items()) + raise ValueError( + f"Pooled models were estimated to different sample ends: {stamps}; pooling compares " + "forecasts made from a single origin, so every model must end on the same date." + ) + + origin = next(iter(ends.values())) + if holdout.index[0] <= origin: + raise ValueError( + f"The holdout starts at {holdout.index[0].date()} but the models are estimated through " + f"{origin.date()}; the holdout must postdate the estimation sample or the scores are in-sample." + ) + _check_alignment(fits, holdout, origin) + return origin + + +# -------------------------------------------------------------------------- +# Mixture sampling +# -------------------------------------------------------------------------- + + +def _mixture_draws( + stacks: list[np.ndarray], + weights: np.ndarray, + n_draws: int | None, + rng: np.random.Generator, +) -> tuple[np.ndarray, np.ndarray]: + """Draw a pooled sample by picking a model per draw, then a draw within it. + + Returns: + Tuple of `(pooled draws of shape (N, H, n), membership of shape (N,))`. + """ + sizes = np.array([stack.shape[0] for stack in stacks]) + total = int(sizes.min()) if n_draws is None else int(n_draws) + membership = rng.choice(len(stacks), size=total, p=weights) + position = np.floor(rng.random(total) * sizes[membership]).astype(int) + pooled = np.empty((total, *stacks[0].shape[1:])) + for i, stack in enumerate(stacks): + picked = membership == i + if picked.any(): + pooled[picked] = stack[position[picked]] + return pooled, membership + + +def _pooled_forecast_result( + pooled: np.ndarray, + membership: np.ndarray, + labels: list[str], + var_names: list[str], + time_index: pd.DatetimeIndex | None = None, +) -> ForecastResult: + """Wrap pooled draws as a `(chain=1, draw=N, step, variable)` ForecastResult.""" + import xarray as xr + + steps = pooled.shape[1] + coords: dict[str, object] = { + "variable": var_names, + "model": ("draw", np.asarray(labels, dtype=object)[membership]), + } + if time_index is not None: + coords["time"] = ("step", time_index.to_numpy()) + forecast = xr.DataArray( + pooled[np.newaxis], + dims=["chain", "draw", "step", "variable"], + coords=coords, + name="forecast", + ) + idata = az.InferenceData(posterior_predictive=xr.Dataset({"forecast": forecast})) + return ForecastResult(idata=idata, steps=steps, var_names=list(var_names), mode="density") + + +def _flatten(forecast: ForecastResult) -> np.ndarray: + """Collapse a ForecastResult's `(chain, draw, step, variable)` draws to `(S, H, n)`.""" + values = forecast.idata.posterior_predictive["forecast"].values + return values.reshape(-1, *values.shape[2:]) + + +# -------------------------------------------------------------------------- +# Public API +# -------------------------------------------------------------------------- + + +class PredictivePool(ImpulsoBaseModel): + """A set of fitted models with one weight each, defining a pooled predictive. + + Produced by `pool_forecasts`, which scores every model's density forecast + on a held-out window. The weights are *static* — one per model, fixed + across horizons and across time — and they are frozen once estimated: + `combine` applies them to new forecasts without rescoring. + + Attributes: + weights: Weight per model, indexed by label, summing to 1. + log_scores: Log predictive density per held-out date (index) and + model (columns) — the contract between scoring and weighting. + method: Which weight rule produced `weights`. + density: How forecast draws were turned into a predictive density. + var_names: Endogenous variables, in the order every model shares. + steps: Number of held-out periods scored, `H`. + origin: Shared forecast origin — the last date of the estimation + sample every model was fitted on. + holdout_predictive: The pooled predictive over the *held-out window*, + as a `ForecastResult` with a `(chain=1, draw=N, step, variable)` + layout. For genuine forecasts, refit on the full sample and use + `combine`. + membership: Which model produced each pooled draw, as an index into + `labels`. Read-only. + converged: Whether the weight solver converged. Always True for + `method="log_score"`, which is closed-form. + optimiser_message: The solver's own status message. + """ + + weights: pd.Series + log_scores: pd.DataFrame = Field(repr=False) + method: WeightMethod + density: DensityKind + var_names: list[str] + steps: int + origin: pd.Timestamp + holdout_predictive: ForecastResult = Field(repr=False) + membership: np.ndarray = Field(repr=False) + converged: bool = True + optimiser_message: str = "" + + @model_validator(mode="after") + def _validate(self) -> PredictivePool: + weights = self.weights.rename_axis("model") + log_scores = self.log_scores.rename_axis(columns="model") + if list(weights.index) != list(log_scores.columns): + raise ValueError( + f"Weight labels {list(weights.index)} do not match the log-score columns {list(log_scores.columns)}." + ) + _check_reserved_labels(list(weights.index)) + values = weights.to_numpy(dtype=float) + if bool((values < -1e-12).any()): + raise ValueError(f"Pool weights must be non-negative, got {values.tolist()}.") + if not np.isclose(values.sum(), 1.0, atol=1e-8): + raise ValueError(f"Pool weights must sum to 1, got {values.sum()!r}.") + + membership = np.asarray(self.membership) + if membership.ndim != 1: + raise ValueError(f"membership must be 1-D, got {membership.ndim}-D.") + if membership.size and (membership.min() < 0 or membership.max() >= values.size): + raise ValueError( + f"membership indexes {values.size} models but ranges over [{membership.min()}, {membership.max()}]." + ) + membership = membership.astype(int, copy=True) + membership.flags.writeable = False + + object.__setattr__(self, "weights", weights) + object.__setattr__(self, "log_scores", log_scores) + object.__setattr__(self, "membership", membership) + return self + + @property + def labels(self) -> list[str]: + """Model labels, in the order the pool was built.""" + return list(self.weights.index) + + def pooled_log_score(self) -> float: + """Total held-out log score of the pooled predictive.""" + return float(_pooled_row_scores(self.log_scores.to_numpy(), self.weights.to_numpy()).sum()) + + def summary(self) -> pd.DataFrame: + """Per-model weights and held-out scores, heaviest weight first. + + Returns: + DataFrame indexed by label with `weight`, `log_score` (total over + the holdout), `mean_log_score` (per held-out period), and `rank`. + """ + totals = self.log_scores.sum(axis=0) + frame = pd.DataFrame({ + "weight": self.weights, + "log_score": totals, + "mean_log_score": totals / self.steps, + }).sort_values("weight", ascending=False) + frame["rank"] = np.arange(1, len(frame) + 1) + return frame + + def to_dataframe(self) -> pd.DataFrame: + """Per-date log scores for every model plus the pooled predictive. + + The combined predictive occupies a `"pooled"` column, which is why + `"pooled"` is rejected as a model label. + """ + frame = self.log_scores.copy() + frame[_POOLED_COLUMN] = _pooled_row_scores(self.log_scores.to_numpy(), self.weights.to_numpy()) + return frame + + def realised_weights(self) -> pd.Series: + """Empirical share of each model among the pooled draws. + + The Monte Carlo counterpart of `weights`: it converges to them as the + pooled sample grows, and the gap is sampling noise, not a second + estimate. + """ + counts = np.bincount(self.membership, minlength=len(self.labels)) + return pd.Series(counts / counts.sum(), index=self.weights.index, name="realised_weight") + + def combine( + self, + forecasts: Mapping[str, ForecastResult], + n_draws: int | None = None, + seed: int | np.random.Generator | None = None, + ) -> ForecastResult: + """Apply the frozen weights to a new set of forecasts. + + The usual workflow: estimate weights on a held-out window, refit every + model on the full sample, forecast past the end of the data, and + combine. The new forecasts need not have the pool's horizon — only the + same models, the same variables, and the same horizon as each other. + + Args: + forecasts: One density-mode `ForecastResult` per pooled model, + keyed by the pool's labels. + n_draws: Size of the pooled sample. Defaults to the smallest + member's draw count. + seed: RNG seed (int) or Generator for the mixture draws. + + Returns: + ForecastResult holding the pooled draws, with a `model` coordinate + recording which member produced each draw. + + Raises: + ValueError: If the labels, variables, horizons, or forecast mode + do not line up with the pool. + """ + labels = self.labels + if set(forecasts) != set(labels): + raise ValueError(f"combine() needs one forecast per pooled model {labels}, got {sorted(forecasts)}.") + steps = {forecasts[label].steps for label in labels} + if len(steps) > 1: + raise ValueError( + f"All forecasts must run to the same number of steps, got {sorted(steps)}; the pooled " + "draws are a mixture over a common horizon." + ) + for label in labels: + forecast = forecasts[label] + if list(forecast.var_names) != self.var_names: + raise ValueError( + f"Forecast {label!r} forecasts variables {list(forecast.var_names)} but the pool " + f"was estimated on {self.var_names}." + ) + if forecast.mode != "density": + raise ValueError( + f"Forecast {label!r} is a mean forecast; pooling combines predictive densities, so " + "call forecast(include_shock_uncertainty=True) for every member." + ) + if n_draws is not None and n_draws < 1: + raise ValueError(f"n_draws must be at least 1, got {n_draws}.") + + rng = _resolve_rng(seed) + stacks = [_flatten(forecasts[label]) for label in labels] + pooled, membership = _mixture_draws(stacks, self.weights.to_numpy(dtype=float), n_draws, rng) + return _pooled_forecast_result(pooled, membership, labels, self.var_names) + + def plot(self) -> Figure: + """Plot the pool weights as a ranked bar chart.""" + from impulso.plotting import plot_pool_weights + + return plot_pool_weights(self) + + +def pool_forecasts( + fits: Mapping[str, FittedVAR], + holdout: VARData, + *, + method: WeightMethod = "stacking", + density: DensityKind = "gaussian", + n_draws: int | None = None, + seed: int | np.random.Generator | None = None, +) -> PredictivePool: + """Weight several fitted models by their held-out predictive performance. + + Every model is forecast `H = len(holdout.index)` steps from the shared + forecast origin (the last date of the estimation sample, which all models + must share), each forecast's predictive density is scored at the held-out + realisations, and the resulting `(H, M)` log-score matrix is turned into + weights. The pool forecasts the models itself rather than accepting + `ForecastResult`s, because only a `FittedVAR` carries the estimation + metadata that makes "same origin, genuinely held out" checkable. + + Args: + fits: Fitted models keyed by label. At least two, all estimated on + the same variables in the same order and to the same sample end. + holdout: The held-out window — the `H` periods immediately after the + forecast origin. Must carry future exogenous values if any model + was fitted with exogenous regressors. + method: `"stacking"` (default) maximises the log score of the *pooled* + predictive and keeps complementary models; `"log_score"` takes a + softmax of each model's total score and collapses onto the single + best model as `H` grows. + density: `"gaussian"` (default) matches a joint Gaussian to each + horizon's forecast draws; `"diagonal"` scores each variable on its + own marginal, which is the escape hatch when the joint covariance + is near-singular or the draw count is small. + n_draws: Size of the pooled predictive sample. Defaults to the + smallest member's draw count. + seed: RNG seed (int) or Generator. Child generators are spawned in + `fits` insertion order, so results are reproducible but depend on + that order. + + Returns: + PredictivePool holding the weights, the score matrix, and a pooled + predictive sample over the held-out window. + + Raises: + ValueError: On fewer than two models, mismatched variables, differing + estimation ends, a holdout that does not immediately follow the + origin, missing exogenous values, a degenerate predictive density, + or an unknown `method`/`density`. + + Note: + The score is an approximation, and comparable across models rather + than exact. Each horizon's predictive density is a Gaussian matched to + the forecast draws, while the true posterior predictive is a + heavier-tailed mixture over draws; the gap widens with few draws, + stochastic volatility, and fat tails. Scores are summed over horizons + `1..H` from one fixed origin — not a joint-path density, and not a + rolling-origin evaluation. Weights are static, so a short holdout + makes them noisy, and a model that only wins late in the window + cannot be given a horizon-specific weight. Finally, `holdout` is + checked to postdate the estimation sample, but nothing can check that + it was not peeked at while the candidates were being chosen. + """ + if method not in ("stacking", "log_score"): + raise ValueError(f"method must be 'stacking' or 'log_score', got {method!r}.") + if density not in ("gaussian", "diagonal"): + raise ValueError(f"density must be 'gaussian' or 'diagonal', got {density!r}.") + if n_draws is not None and n_draws < 1: + raise ValueError(f"n_draws must be at least 1, got {n_draws}.") + + origin = _validate_pool_inputs(fits, holdout) + labels = list(fits) + steps = len(holdout.index) + realised = np.asarray(holdout.endog, dtype=float) + + rng = _resolve_rng(seed) + children = _spawn(rng, len(labels) + 1) + + stacks: list[np.ndarray] = [] + columns: list[np.ndarray] = [] + for child, label in zip(children[:-1], labels, strict=True): + fit = fits[label] + exog_future = np.asarray(holdout.exog, dtype=float) if fit.has_exog else None + forecast = fit.forecast( + steps=steps, + include_shock_uncertainty=True, + seed=child, + exog_future=exog_future, + ) + draws = _flatten(forecast) + stacks.append(draws) + columns.append(_gaussian_log_scores(draws, realised, density=density, label=label)) + + log_scores = pd.DataFrame( + np.column_stack(columns), + index=pd.DatetimeIndex(holdout.index, name="time"), + columns=pd.Index(labels, name="model"), + ) + matrix = log_scores.to_numpy() + if method == "stacking": + weight_values, converged, message = _stacking_weights(matrix, index=log_scores.index) + else: + weight_values, converged, message = _log_score_weights(matrix, index=log_scores.index), True, _NO_OPTIMISER + + pooled, membership = _mixture_draws(stacks, weight_values, n_draws, children[-1]) + return PredictivePool( + weights=pd.Series(weight_values, index=pd.Index(labels, name="model"), name="weight"), + log_scores=log_scores, + method=method, + density=density, + var_names=list(holdout.endog_names), + steps=steps, + origin=pd.Timestamp(origin), + holdout_predictive=_pooled_forecast_result( + pooled, membership, labels, list(holdout.endog_names), pd.DatetimeIndex(holdout.index) + ), + membership=membership, + converged=converged, + optimiser_message=message, + ) diff --git a/tests/test_plotting.py b/tests/test_plotting.py index 559afb5..21a8959 100644 --- a/tests/test_plotting.py +++ b/tests/test_plotting.py @@ -172,6 +172,42 @@ def test_plot_volatility_returns_figure(synthetic_sv_idata): assert isinstance(fig, Figure) +def test_plot_pool_weights_returns_figure(): + """Smoke test: the pool bar chart renders from a hand-built pool.""" + import pandas as pd + from matplotlib.figure import Figure + + from impulso.plotting import plot_pool_weights + from impulso.pooling import PredictivePool + + labels = ["a", "b"] + index = pd.DatetimeIndex(pd.date_range("2020-01-01", periods=3, freq="QS"), name="time") + log_scores = pd.DataFrame([[-1.0, -2.0], [-1.5, -2.5], [-0.5, -3.0]], index=index, columns=labels) + draws = np.zeros((4, 3, 2)) + da = xr.DataArray( + draws[np.newaxis], + dims=["chain", "draw", "step", "variable"], + coords={"variable": ["y1", "y2"], "model": ("draw", np.array(["a", "a", "b", "b"], dtype=object))}, + name="forecast", + ) + pool = PredictivePool( + weights=pd.Series([0.7, 0.3], index=labels), + log_scores=log_scores, + method="stacking", + density="gaussian", + var_names=["y1", "y2"], + steps=3, + origin=pd.Timestamp("2019-10-01"), + holdout_predictive=ForecastResult( + idata=az.InferenceData(posterior_predictive=xr.Dataset({"forecast": da})), + steps=3, + var_names=["y1", "y2"], + ), + membership=np.array([0, 0, 1, 1]), + ) + assert isinstance(plot_pool_weights(pool), Figure) + + def _make_sv_forecast_result(steps=12, index=None): from impulso.results import SVForecastResult diff --git a/tests/test_pooling.py b/tests/test_pooling.py new file mode 100644 index 0000000..228dae3 --- /dev/null +++ b/tests/test_pooling.py @@ -0,0 +1,1054 @@ +"""Tests for static predictive-density pooling (issue #151). + +All fast: the fixtures hand-build `FittedVAR` posteriors rather than +sampling, following the `test_density_forecast.py` precedent. +""" + +import warnings + +import arviz as az +import matplotlib +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from pydantic import ValidationError + +matplotlib.use("Agg") + +from impulso.data import VARData +from impulso.fitted import FittedVAR +from impulso.pooling import ( + PredictivePool, + _gaussian_log_scores, + _index_freq, + _log_score_weights, + _pooled_row_scores, + _spawn, + _stacking_weights, + pool_forecasts, +) +from impulso.results import ForecastResult +from impulso.volatility import Constant + +A1_DEFAULT = np.array([[0.5, 0.1], [-0.2, 0.3]]) + +# log([[4, 1], [4, 1], [1, 4]]): the analytic reference matrix. +# stacking -> FOC 2*3/(3w + 1) = 3/(4 - 3w) -> w* = 7/9 +# log score -> products 16 vs 4 -> [0.8, 0.2] +REFERENCE_SCORES = np.log(np.array([[4.0, 1.0], [4.0, 1.0], [1.0, 4.0]])) + + +# -------------------------------------------------------------------------- +# Fixtures +# -------------------------------------------------------------------------- + + +def _fitted( + sd, + *, + n_lags: int = 1, + A1: np.ndarray | None = None, + intercept: float = 0.0, + n_chains: int = 2, + n_draws: int = 100, + n_obs: int = 60, + start: str = "2000-01-01", + freq: str = "QS", + var_names: tuple[str, ...] = ("y1", "y2"), + exog: bool = False, +) -> FittedVAR: + """Hand-build a FittedVAR with a point-mass posterior on (B, intercept, L). + + `sd` is the diagonal of the shock Cholesky factor, so the predictive + spread of each variable is controlled exactly. + """ + n_vars = len(var_names) + coefs = np.asarray(A1 if A1 is not None else A1_DEFAULT, dtype=float) + B = np.broadcast_to(coefs, (n_chains, n_draws, n_vars, n_vars * n_lags)).copy() + mu = np.full((n_chains, n_draws, n_vars), float(intercept)) + L = np.zeros((n_chains, n_draws, n_vars, n_vars)) + L[:, :, range(n_vars), range(n_vars)] = np.asarray(sd, dtype=float) + variables = { + "B": (("chain", "draw", "var", "coeff"), B), + "intercept": (("chain", "draw", "var"), mu), + "L": (("chain", "draw", "var1", "var2"), L), + } + if exog: + variables["B_exog"] = (("chain", "draw", "var", "exog"), np.zeros((n_chains, n_draws, n_vars, 1))) + idata = az.InferenceData(posterior=xr.Dataset(variables)) + + index = pd.date_range(start, periods=n_obs, freq=freq) + y = np.zeros((n_obs, n_vars)) + y[0] = 1.0 + for t in range(1, n_obs): + y[t] = intercept + coefs @ y[t - 1] + exog_arr = np.arange(n_obs, dtype=float).reshape(-1, 1) if exog else None + data = VARData( + endog=y, + endog_names=list(var_names), + exog=exog_arr, + exog_names=["x"] if exog else None, + index=index, + ) + return FittedVAR( + idata=idata, + n_lags=n_lags, + data=data, + var_names=list(var_names), + volatility=Constant(), + ) + + +def _mean_path(fit: FittedVAR, steps: int) -> np.ndarray: + """Closed-form conditional mean path A1^h y_T (intercept 0 fixtures).""" + A1 = fit.idata.posterior["B"].values[0, 0] + y = fit.data.endog[-1].copy() + out = np.empty((steps, y.size)) + for h in range(steps): + y = A1 @ y + out[h] = y + return out + + +def _holdout(fit: FittedVAR, values: np.ndarray, *, exog: bool = False) -> VARData: + """Held-out VARData immediately following a fit's estimation sample.""" + steps = values.shape[0] + index = pd.date_range(fit.data.index[-1], periods=steps + 1, freq=fit.data.index.freq)[1:] + return VARData( + endog=np.asarray(values, dtype=float), + endog_names=list(fit.var_names), + exog=np.arange(steps, dtype=float).reshape(-1, 1) if exog else None, + exog_names=["x"] if exog else None, + index=index, + ) + + +@pytest.fixture +def mirrored_pool(): + """Two models with identical means and mirrored shock scales. + + The headline stacking case: neither model dominates, the holdout + alternates between the regions each model covers, and the pooled score + beats both members by a wide margin. + """ + tight_then_wide = _fitted(sd=[0.3, 1.5]) + wide_then_tight = _fitted(sd=[1.5, 0.3]) + steps = 12 + mean = _mean_path(tight_then_wide, steps) + deviation = np.tile(np.array([[0.0, 2.5], [2.5, 0.0]]), (steps // 2, 1)) + holdout = _holdout(tight_then_wide, mean + deviation) + return {"tight_y1": tight_then_wide, "tight_y2": wide_then_tight}, holdout + + +@pytest.fixture +def dominance_pool(): + """Three models where one is unambiguously best on the holdout.""" + good = _fitted(sd=[0.5, 0.5]) + wide = _fitted(sd=[5.0, 5.0]) + biased = _fitted(sd=[0.5, 0.5], intercept=3.0) + steps = 8 + holdout = _holdout(good, _mean_path(good, steps)) + return {"good": good, "wide": wide, "biased": biased}, holdout + + +# -------------------------------------------------------------------------- +# A. Weight solvers against closed-form optima +# -------------------------------------------------------------------------- + + +class TestWeightSolvers: + def test_stacking_interior_optimum(self): + weights, converged, message = _stacking_weights(REFERENCE_SCORES) + assert converged + assert message + np.testing.assert_allclose(weights, [7.0 / 9.0, 2.0 / 9.0], atol=1e-6) + + def test_log_score_closed_form(self): + np.testing.assert_allclose(_log_score_weights(REFERENCE_SCORES), [0.8, 0.2], atol=1e-12) + + def test_dominance_collapses_to_one_model(self): + scores = np.array([[0.0, -50.0], [0.0, -50.0], [0.0, -50.0]]) + np.testing.assert_allclose(_stacking_weights(scores)[0], [1.0, 0.0], atol=1e-8) + np.testing.assert_allclose(_log_score_weights(scores), [1.0, 0.0], atol=1e-8) + + def test_symmetric_rows_split_evenly(self): + scores = np.log(np.array([[4.0, 1.0], [1.0, 4.0]])) + np.testing.assert_allclose(_stacking_weights(scores)[0], [0.5, 0.5], atol=1e-6) + np.testing.assert_allclose(_log_score_weights(scores), [0.5, 0.5], atol=1e-12) + + @pytest.mark.parametrize("method", ["stacking", "log_score"]) + @pytest.mark.parametrize("seed", [0, 1, 2, 3]) + def test_weights_live_on_the_simplex(self, method, seed): + rng = np.random.default_rng(seed) + scores = rng.standard_normal((7, 4)) * 5.0 + weights = _stacking_weights(scores)[0] if method == "stacking" else _log_score_weights(scores) + assert weights.shape == (4,) + assert np.all(weights >= 0.0) + assert np.isclose(weights.sum(), 1.0, atol=1e-10) + + def test_stacking_keeps_a_model_that_dies_at_one_point(self): + """A -inf at one point kills log-score weights but not stacking.""" + scores = np.array([[0.0, -np.inf], [0.0, 0.0], [-3.0, 0.0]]) + weights, _, _ = _stacking_weights(scores) + assert np.all(np.isfinite(weights)) + # Analytic optimum: w1 = 1 / (2 * (1 - exp(-3))). + np.testing.assert_allclose(weights[0], 1.0 / (2.0 * (1.0 - np.exp(-3.0))), atol=1e-6) + assert weights[1] > 0.4 + np.testing.assert_allclose(_log_score_weights(scores), [1.0, 0.0], atol=1e-12) + + def test_dead_point_raises_naming_the_point(self): + scores = np.array([[0.0, 0.0], [-np.inf, -np.inf], [0.0, 0.0]]) + index = pd.DatetimeIndex(["2020-01-01", "2020-04-01", "2020-07-01"]) + with pytest.raises(ValueError, match="2020-04-01"): + _stacking_weights(scores, index=index) + with pytest.raises(ValueError, match="2020-04-01"): + _log_score_weights(scores, index=index) + + def test_dead_point_without_index_names_the_step(self): + scores = np.array([[0.0, 0.0], [-np.inf, -np.inf]]) + with pytest.raises(ValueError, match="holdout step 2"): + _log_score_weights(scores) + + def test_every_total_infinite_rejects_log_score_weights(self): + """Each row survives, but each column dies somewhere.""" + scores = np.array([[0.0, -np.inf], [-np.inf, 0.0]]) + with pytest.raises(ValueError, match="method='stacking'"): + _log_score_weights(scores) + assert np.all(np.isfinite(_stacking_weights(scores)[0])) + + def test_single_model_matrix_rejected(self): + with pytest.raises(ValueError, match="at least two fitted models"): + _log_score_weights(np.zeros((3, 1))) + + def test_non_2d_matrix_rejected(self): + with pytest.raises(ValueError, match="must be 2-D"): + _log_score_weights(np.zeros(3)) + + def test_nan_scores_rejected(self): + with pytest.raises(ValueError, match="NaN or \\+inf"): + _stacking_weights(np.array([[0.0, np.nan], [0.0, 0.0]])) + + def test_optimiser_retries_then_gives_up(self, monkeypatch): + import scipy.optimize + + calls = [] + real = scipy.optimize.minimize + + def always_fails(*args, **kwargs): + calls.append(kwargs.get("x0")) + result = real(*args, **kwargs) + result.success = False + result.message = "synthetic failure" + return result + + monkeypatch.setattr(scipy.optimize, "minimize", always_fails) + with pytest.raises(RuntimeError, match="synthetic failure"): + _stacking_weights(REFERENCE_SCORES) + assert len(calls) == 2 + # The retry starts somewhere other than the uniform point. + assert not np.allclose(calls[0], calls[1]) + + def test_optimiser_retry_can_succeed(self, monkeypatch): + import scipy.optimize + + real = scipy.optimize.minimize + state = {"n": 0} + + def fails_once(*args, **kwargs): + state["n"] += 1 + result = real(*args, **kwargs) + if state["n"] == 1: + result.success = False + return result + + monkeypatch.setattr(scipy.optimize, "minimize", fails_once) + weights, converged, _ = _stacking_weights(REFERENCE_SCORES) + assert converged + assert state["n"] == 2 + np.testing.assert_allclose(weights, [7.0 / 9.0, 2.0 / 9.0], atol=1e-6) + + +# -------------------------------------------------------------------------- +# B. Overflow safety +# -------------------------------------------------------------------------- + + +class TestOverflow: + @pytest.mark.parametrize("shift", [-5e4, 5e4]) + def test_log_score_weights_survive_extreme_shifts(self, shift): + scores = REFERENCE_SCORES + shift + with warnings.catch_warnings(), np.errstate(over="raise", invalid="raise"): + warnings.simplefilter("error") + weights = _log_score_weights(scores) + assert np.all(np.isfinite(weights)) + np.testing.assert_allclose(weights, [0.8, 0.2], atol=1e-12) + + @pytest.mark.parametrize("shift", [-1e5, 1e4]) + def test_stacking_survives_extreme_shifts(self, shift): + """An unshifted exponentiation (as ArviZ does) fails here.""" + scores = REFERENCE_SCORES + shift + with warnings.catch_warnings(): + warnings.simplefilter("error") + weights, converged, _ = _stacking_weights(scores) + assert converged + assert np.all(np.isfinite(weights)) + np.testing.assert_allclose(weights, [7.0 / 9.0, 2.0 / 9.0], atol=1e-6) + + def test_naive_exponentiation_would_have_failed(self): + """Documents why the shift is load-bearing rather than cosmetic.""" + with np.errstate(over="ignore", under="ignore"): + assert not np.isfinite(np.exp(REFERENCE_SCORES + 1e4)).all() + assert np.exp(REFERENCE_SCORES - 1e5).max() == 0.0 + + +# -------------------------------------------------------------------------- +# C. Score matrix +# -------------------------------------------------------------------------- + + +class TestScoreMatrix: + def test_joint_gaussian_matches_scipy(self): + from scipy.stats import multivariate_normal + + rng = np.random.default_rng(7) + draws = rng.standard_normal((500, 4, 3)) @ np.array([[1.0, 0.0, 0.0], [0.6, 0.9, 0.0], [0.2, 0.3, 1.1]]).T + y = rng.standard_normal((4, 3)) + scores = _gaussian_log_scores(draws, y, density="gaussian", label="m") + expected = [ + multivariate_normal.logpdf(y[h], draws[:, h, :].mean(0), np.cov(draws[:, h, :], rowvar=False, ddof=1)) + for h in range(4) + ] + np.testing.assert_allclose(scores, expected, atol=1e-6) + + def test_diagonal_matches_scipy(self): + from scipy.stats import norm + + rng = np.random.default_rng(11) + draws = rng.standard_normal((300, 3, 2)) * np.array([1.0, 2.5]) + y = rng.standard_normal((3, 2)) + scores = _gaussian_log_scores(draws, y, density="diagonal", label="m") + expected = [ + float(np.sum(norm.logpdf(y[h], draws[:, h, :].mean(0), draws[:, h, :].std(0, ddof=1)))) for h in range(3) + ] + np.testing.assert_allclose(scores, expected, atol=1e-10) + + def test_diagonal_ignores_correlation(self): + """The escape hatch really is the product of marginals.""" + rng = np.random.default_rng(3) + base = rng.standard_normal((400, 2, 2)) + correlated = base @ np.array([[1.0, 0.0], [0.95, 0.31]]).T + y = np.zeros((2, 2)) + joint = _gaussian_log_scores(correlated, y, density="gaussian", label="m") + diagonal = _gaussian_log_scores(correlated, y, density="diagonal", label="m") + assert not np.allclose(joint, diagonal) + + def test_singular_covariance_points_at_diagonal(self): + rng = np.random.default_rng(5) + draws = rng.standard_normal((200, 2, 2)) + draws[:, 1, 1] = 3.0 # constant column at horizon 2 + with pytest.raises(ValueError, match="density='diagonal'"): + _gaussian_log_scores(draws, np.zeros((2, 2)), density="gaussian", label="alpha") + + def test_collinear_columns_rejected(self): + rng = np.random.default_rng(5) + draws = rng.standard_normal((200, 1, 2)) + draws[:, 0, 1] = draws[:, 0, 0] + with pytest.raises(ValueError, match="numerically singular"): + _gaussian_log_scores(draws, np.zeros((1, 2)), density="gaussian", label="alpha") + + def test_all_draws_identical_rejected(self): + draws = np.ones((50, 1, 2)) + with pytest.raises(ValueError, match="every forecast draw is identical"): + _gaussian_log_scores(draws, np.zeros((1, 2)), density="gaussian", label="alpha") + + def test_too_few_draws_for_joint_density(self): + rng = np.random.default_rng(2) + with pytest.raises(ValueError, match="S <= n"): + _gaussian_log_scores(rng.standard_normal((3, 2, 3)), np.zeros((2, 3)), label="alpha") + + def test_single_draw_rejected(self): + with pytest.raises(ValueError, match="at least two"): + _gaussian_log_scores(np.zeros((1, 2, 2)), np.zeros((2, 2)), density="diagonal", label="alpha") + + def test_zero_marginal_variance_rejected(self): + rng = np.random.default_rng(9) + draws = rng.standard_normal((40, 2, 2)) + draws[:, 0, 1] = 2.0 + with pytest.raises(ValueError, match="zero forecast variance"): + _gaussian_log_scores(draws, np.zeros((2, 2)), density="diagonal", label="alpha") + + def test_non_finite_draws_rejected(self): + draws = np.ones((10, 2, 2)) + draws[0, 0, 0] = np.nan + with pytest.raises(ValueError, match="non-finite forecast draws"): + _gaussian_log_scores(draws, np.zeros((2, 2)), label="alpha") + + def test_holdout_shape_mismatch_rejected(self): + with pytest.raises(ValueError, match="but the draws imply"): + _gaussian_log_scores(np.zeros((10, 2, 2)), np.zeros((3, 2)), label="alpha") + + def test_pooled_row_scores_tolerate_zero_weights(self): + scores = np.array([[0.0, -1.0], [0.0, -1.0]]) + with warnings.catch_warnings(): + warnings.simplefilter("error") + pooled = _pooled_row_scores(scores, np.array([1.0, 0.0])) + np.testing.assert_allclose(pooled, [0.0, 0.0]) + + +# -------------------------------------------------------------------------- +# D. pool_forecasts end to end +# -------------------------------------------------------------------------- + + +class TestPoolForecasts: + @pytest.mark.parametrize("method", ["stacking", "log_score"]) + def test_dominant_model_takes_all_the_weight(self, dominance_pool, method): + fits, holdout = dominance_pool + pool = pool_forecasts(fits, holdout, method=method, seed=0) + assert pool.weights["good"] > 0.99 + assert pool.log_scores.sum().idxmax() == "good" + + def test_stacking_splits_mirrored_models_evenly(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, method="stacking", seed=0) + np.testing.assert_allclose(pool.weights.to_numpy(), [0.5, 0.5], atol=1e-3) + + def test_pooling_beats_every_single_model(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, method="stacking", seed=0) + assert pool.pooled_log_score() > pool.log_scores.sum().max() + 50.0 + + def test_the_two_methods_disagree(self, mirrored_pool): + """Log-score weights collapse where stacking keeps both models.""" + fits, holdout = mirrored_pool + stacked = pool_forecasts(fits, holdout, method="stacking", seed=0) + log_scored = pool_forecasts(fits, holdout, method="log_score", seed=0) + assert stacked.weights.min() > 0.4 + assert log_scored.weights.min() < 0.05 + assert stacked.pooled_log_score() > log_scored.pooled_log_score() + + def test_stacking_weights_resolve_scipy_independently(self, mirrored_pool): + """Version-proof check: re-solve the published score matrix directly.""" + from scipy.optimize import Bounds, LinearConstraint, minimize + + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, method="stacking", seed=0) + matrix = pool.log_scores.to_numpy() + densities = np.exp(matrix - matrix.max(axis=1, keepdims=True)) + result = minimize( + lambda w: -np.sum(np.log(densities @ w)), + x0=np.array([0.3, 0.7]), + method="SLSQP", + bounds=Bounds(0.0, 1.0), + constraints=LinearConstraint(np.ones((1, 2)), 1.0, 1.0), + options={"ftol": 1e-12, "maxiter": 1000}, + ) + np.testing.assert_allclose(pool.weights.to_numpy(), result.x, atol=1e-5) + + def test_asymmetric_holdout_shifts_the_weights(self, mirrored_pool): + """Three quarters of the holdout favours one model; check against scipy.""" + fits, _ = mirrored_pool + base = next(iter(fits.values())) + steps = 12 + mean = _mean_path(base, steps) + pattern = np.tile(np.array([[0.0, 2.5], [0.0, 2.5], [0.0, 2.5], [2.5, 0.0]]), (steps // 4, 1)) + pool = pool_forecasts(fits, _holdout(base, mean + pattern), method="stacking", seed=0) + assert pool.weights.iloc[0] > pool.weights.iloc[1] + assert pool.weights.min() > 0.0 + + def test_log_scores_frame_is_labelled(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + pd.testing.assert_index_equal(pool.log_scores.index, pd.DatetimeIndex(holdout.index, name="time")) + assert list(pool.log_scores.columns) == list(fits) + assert pool.log_scores.columns.name == "model" + assert pool.labels == list(fits) + + def test_summary_layout(self, dominance_pool): + fits, holdout = dominance_pool + pool = pool_forecasts(fits, holdout, seed=0) + summary = pool.summary() + assert list(summary.columns) == ["weight", "log_score", "mean_log_score", "rank"] + assert summary.index[0] == "good" + assert summary["rank"].tolist() == [1, 2, 3] + assert summary["weight"].is_monotonic_decreasing + np.testing.assert_allclose( + summary["mean_log_score"].to_numpy(), + summary["log_score"].to_numpy() / pool.steps, + ) + + def test_to_dataframe_carries_the_pooled_column(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + frame = pool.to_dataframe() + assert list(frame.columns) == [*fits, "pooled"] + np.testing.assert_allclose(frame["pooled"].sum(), pool.pooled_log_score()) + # A mixture density is bounded by its members: below the best member's + # own score at each point, above that member's score shrunk by its weight. + members = frame[list(fits)] + assert (frame["pooled"] <= members.max(axis=1) + 1e-9).all() + assert (frame["pooled"] >= (members + np.log(pool.weights)).max(axis=1) - 1e-9).all() + # Summed over the window, pooling wins — that is the point. + assert frame["pooled"].sum() > members.sum().max() + + def test_origin_is_the_shared_estimation_end(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + assert pool.origin == next(iter(fits.values())).data.index[-1] + assert pool.steps == len(holdout.index) + assert pool.var_names == holdout.endog_names + assert pool.method == "stacking" + assert pool.density == "gaussian" + assert pool.converged + + def test_diagonal_density_runs(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, density="diagonal", seed=0) + assert pool.density == "diagonal" + assert np.isclose(pool.weights.sum(), 1.0) + + def test_heterogeneous_chain_counts_pool(self, mirrored_pool): + """A conjugate-style (1, 200) posterior pools with a NUTS-style (2, 100).""" + _, holdout = mirrored_pool + fits = { + "conjugate": _fitted(sd=[0.3, 1.5], n_chains=1, n_draws=200), + "nuts": _fitted(sd=[1.5, 0.3], n_chains=2, n_draws=100), + } + pool = pool_forecasts(fits, holdout, seed=0) + assert pool.membership.shape == (200,) + assert np.isclose(pool.weights.sum(), 1.0) + + def test_exogenous_models_consume_the_holdout_exog(self): + fit_a = _fitted(sd=[0.4, 1.2], exog=True) + fit_b = _fitted(sd=[1.2, 0.4], exog=True) + holdout = _holdout(fit_a, _mean_path(fit_a, 6), exog=True) + pool = pool_forecasts({"a": fit_a, "b": fit_b}, holdout, seed=0) + assert np.isclose(pool.weights.sum(), 1.0) + + def test_pool_never_refits(self, mirrored_pool): + """The fits handed in are untouched — the pool holds no reference to them.""" + fits, holdout = mirrored_pool + before = {k: v.data.endog.copy() for k, v in fits.items()} + pool_forecasts(fits, holdout, seed=0) + for label, fit in fits.items(): + np.testing.assert_array_equal(fit.data.endog, before[label]) + + +# -------------------------------------------------------------------------- +# E. ArviZ parity +# -------------------------------------------------------------------------- + + +def _arviz_shim(log_scores: pd.DataFrame) -> dict[str, az.InferenceData]: + """Wrap a score matrix as InferenceData whose loo reproduces the scores. + + Each column becomes a `log_likelihood` broadcast constant across draws, + so PSIS weights are uniform, `p_loo` is zero, and `elpd_loo` is exactly + the summed column. This makes `az.compare` a pure function of the score + matrix — the point of the parity check. Verified against arviz 0.23. + """ + rng = np.random.default_rng(0) + n_chains, n_draws = 2, 200 + out = {} + for label in log_scores.columns: + values = np.broadcast_to(log_scores[label].to_numpy(), (n_chains, n_draws, len(log_scores))).copy() + out[label] = az.InferenceData( + posterior=xr.Dataset({"mu": (("chain", "draw"), rng.standard_normal((n_chains, n_draws)))}), + log_likelihood=xr.Dataset({"y": (("chain", "draw", "obs"), values)}), + ) + return out + + +class TestArviZParity: + def test_stacking_matches_az_compare(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + frame = pd.DataFrame(REFERENCE_SCORES, columns=["a", "b"]) + comparison = az.compare(_arviz_shim(frame), ic="loo", method="stacking", scale="log") + np.testing.assert_allclose(comparison.loc["a", "elpd_loo"], REFERENCE_SCORES[:, 0].sum(), atol=1e-8) + weights = _stacking_weights(REFERENCE_SCORES)[0] + # ArviZ optimises a softmax reparameterisation with its own default + # tolerance, so it lands ~3e-5 from the analytic 7/9; the parity + # claim is "same solution", not "same solver settings". + np.testing.assert_allclose( + [comparison.loc["a", "weight"], comparison.loc["b", "weight"]], weights, atol=1e-4 + ) + np.testing.assert_allclose(weights[0], 7.0 / 9.0, atol=1e-6) + + def test_pseudo_bma_matches_az_compare(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + frame = pd.DataFrame(REFERENCE_SCORES, columns=["a", "b"]) + comparison = az.compare(_arviz_shim(frame), ic="loo", method="BB-pseudo-BMA", scale="log") + assert comparison.loc["a", "weight"] > comparison.loc["b", "weight"] + np.testing.assert_allclose(_log_score_weights(REFERENCE_SCORES), [0.8, 0.2], atol=1e-12) + + +# -------------------------------------------------------------------------- +# F. Validation +# -------------------------------------------------------------------------- + + +class TestValidation: + def test_single_model_rejected(self, mirrored_pool): + fits, holdout = mirrored_pool + with pytest.raises(ValueError, match="at least two fitted models"): + pool_forecasts({"only": next(iter(fits.values()))}, holdout) + + def test_variable_name_mismatch(self, mirrored_pool): + fits, holdout = mirrored_pool + other = _fitted(sd=[0.5, 0.5], var_names=("y1", "z9")) + with pytest.raises(ValueError, match="must share the holdout's variables"): + pool_forecasts({**fits, "odd": other}, holdout) + + def test_variable_order_mismatch_suggests_reordering(self, mirrored_pool): + fits, holdout = mirrored_pool + swapped = VARData( + endog=np.asarray(holdout.endog)[:, ::-1], + endog_names=["y2", "y1"], + index=holdout.index, + ) + with pytest.raises(ValueError, match="reorder the holdout columns"): + pool_forecasts(fits, swapped) + + def test_different_estimation_ends(self, mirrored_pool): + fits, holdout = mirrored_pool + shorter = _fitted(sd=[0.5, 0.5], n_obs=59) + with pytest.raises(ValueError, match="different sample ends"): + pool_forecasts({**fits, "short": shorter}, holdout) + + def test_holdout_must_postdate_the_origin(self, mirrored_pool): + fits, holdout = mirrored_pool + overlapping = VARData( + endog=np.asarray(holdout.endog), + endog_names=holdout.endog_names, + index=pd.date_range("2010-01-01", periods=len(holdout.index), freq="QS"), + ) + with pytest.raises(ValueError, match="must postdate the estimation sample"): + pool_forecasts(fits, overlapping) + + def test_gap_after_the_origin_is_rejected(self, mirrored_pool): + fits, holdout = mirrored_pool + origin = next(iter(fits.values())).data.index[-1] + gapped = VARData( + endog=np.asarray(holdout.endog), + endog_names=holdout.endog_names, + index=pd.date_range(origin, periods=len(holdout.index) + 2, freq="QS")[2:], + ) + with pytest.raises(ValueError, match="does not continue the estimation sample"): + pool_forecasts(fits, gapped) + + def test_unknown_frequency_warns_and_proceeds(self): + index = pd.DatetimeIndex(["2000-01-01", "2000-01-05", "2000-02-11", "2000-04-02"]) + fits = { + "a": _fitted(sd=[0.3, 1.5], n_obs=4), + "b": _fitted(sd=[1.5, 0.3], n_obs=4), + } + fits = { + label: FittedVAR( + idata=fit.idata, + n_lags=fit.n_lags, + data=VARData(endog=np.asarray(fit.data.endog), endog_names=fit.var_names, index=index), + var_names=fit.var_names, + volatility=fit.volatility, + ) + for label, fit in fits.items() + } + holdout = VARData( + endog=np.zeros((3, 2)), + endog_names=["y1", "y2"], + index=pd.DatetimeIndex(["2000-05-09", "2000-06-30", "2000-09-01"]), + ) + with pytest.warns(UserWarning, match="estimation samples or the holdout"): + pool = pool_forecasts(fits, holdout, seed=0) + assert np.isclose(pool.weights.sum(), 1.0) + + def test_missing_holdout_exog(self): + fit_a = _fitted(sd=[0.4, 1.2], exog=True) + fit_b = _fitted(sd=[1.2, 0.4], exog=True) + holdout = _holdout(fit_a, _mean_path(fit_a, 6), exog=False) + with pytest.raises(ValueError, match="exogenous regressors"): + pool_forecasts({"a": fit_a, "b": fit_b}, holdout) + + def test_exog_name_mismatch(self): + fit_a = _fitted(sd=[0.4, 1.2], exog=True) + fit_b = _fitted(sd=[1.2, 0.4], exog=True) + steps = 6 + index = pd.date_range(fit_a.data.index[-1], periods=steps + 1, freq="QS")[1:] + holdout = VARData( + endog=_mean_path(fit_a, steps), + endog_names=["y1", "y2"], + exog=np.arange(steps, dtype=float).reshape(-1, 1), + exog_names=["other"], + index=index, + ) + with pytest.raises(ValueError, match="exogenous regressors named"): + pool_forecasts({"a": fit_a, "b": fit_b}, holdout) + + def test_unknown_method(self, mirrored_pool): + fits, holdout = mirrored_pool + with pytest.raises(ValueError, match="method must be"): + pool_forecasts(fits, holdout, method="magic") + + def test_unknown_density(self, mirrored_pool): + fits, holdout = mirrored_pool + with pytest.raises(ValueError, match="density must be"): + pool_forecasts(fits, holdout, density="student") + + def test_n_draws_must_be_positive(self, mirrored_pool): + fits, holdout = mirrored_pool + with pytest.raises(ValueError, match="n_draws must be at least 1"): + pool_forecasts(fits, holdout, n_draws=0) + + def test_bad_seed_is_reported_as_a_seed_problem(self, mirrored_pool): + fits, holdout = mirrored_pool + with pytest.raises(ValueError, match="seed must be"): + pool_forecasts(fits, holdout, seed="tomorrow") + + @pytest.mark.parametrize("order", [("quarterly", "monthly"), ("monthly", "quarterly")]) + def test_mixed_frequency_fits_rejected_in_either_order(self, order): + """Two models can share a sample end at different frequencies. + + The check must not depend on which one the mapping happens to yield + first: previously only the first fit's index was inspected, so a + quarterly-first mapping pooled a monthly model silently. + """ + candidates = { + "quarterly": _fitted(sd=[0.5, 0.5], n_obs=60, start="2000-01-01", freq="QS"), + "monthly": _fitted(sd=[0.6, 0.6], n_obs=60, start="2009-11-01", freq="MS"), + } + assert candidates["quarterly"].data.index[-1] == candidates["monthly"].data.index[-1] + fits = {label: candidates[label] for label in order} + holdout = _holdout(candidates["quarterly"], np.zeros((6, 2))) + with pytest.raises(ValueError, match="different frequencies"): + pool_forecasts(fits, holdout, seed=0) + + def test_one_unknown_frequency_among_regular_fits_warns(self): + """A model whose index has no inferable frequency is flagged, not ignored.""" + regular = _fitted(sd=[0.5, 0.5], n_obs=60, start="2000-01-01", freq="QS") + irregular_index = pd.DatetimeIndex([ + *pd.date_range("2000-01-01", periods=58, freq="QS"), + "2014-08-13", + "2014-10-01", + ]) + irregular = FittedVAR( + idata=regular.idata, + n_lags=regular.n_lags, + data=VARData( + endog=np.asarray(regular.data.endog), + endog_names=regular.var_names, + index=irregular_index, + ), + var_names=regular.var_names, + volatility=regular.volatility, + ) + holdout = _holdout(regular, _mean_path(regular, 6)) + with pytest.warns(UserWarning, match=r"Could not infer a frequency for \['odd'\]"): + pool = pool_forecasts({"regular": regular, "odd": irregular}, holdout, seed=0) + assert np.isclose(pool.weights.sum(), 1.0) + + def test_non_spawnable_generator_is_reported_as_a_seed_problem(self): + with pytest.raises(ValueError, match="supports spawning"): + _spawn(np.random.RandomState(0), 3) + + def test_frequency_is_inferred_when_the_index_does_not_carry_one(self): + regular = pd.DatetimeIndex(["2000-01-01", "2000-04-01", "2000-07-01", "2000-10-01"]) + assert regular.freq is None + assert _index_freq(regular) is not None + assert _index_freq(pd.DatetimeIndex(["2000-01-01", "2000-01-05"])) is None + + def test_reserved_pooled_label_rejected(self, mirrored_pool): + """'pooled' would overwrite that model's column in to_dataframe().""" + fits, holdout = mirrored_pool + renamed = dict(zip(["pooled", "other"], fits.values(), strict=True)) + with pytest.raises(ValueError, match="reserved model label"): + pool_forecasts(renamed, holdout, seed=0) + + def test_reserved_pooled_label_rejected_on_direct_construction(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + broken = dict(pool) + broken["weights"] = pd.Series(pool.weights.to_numpy(), index=["pooled", "other"]) + broken["log_scores"] = pool.log_scores.set_axis(["pooled", "other"], axis=1) + with pytest.raises(ValidationError, match="reserved model label"): + PredictivePool(**broken) + + def test_empty_holdout_rejected(self, mirrored_pool): + fits, _ = mirrored_pool + origin = next(iter(fits.values())).data.index[-1] + with pytest.raises(ValueError, match="at least one held-out"): + pool_forecasts( + fits, + VARData( + endog=np.zeros((0, 2)), + endog_names=["y1", "y2"], + index=pd.date_range(origin, periods=1, freq="QS")[1:], + ), + ) + + def test_literal_fields_are_enforced_on_direct_construction(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + with pytest.raises(ValidationError): + PredictivePool(**{**dict(pool), "method": "magic"}) + + +# -------------------------------------------------------------------------- +# G. Combined sample +# -------------------------------------------------------------------------- + + +class TestCombinedSample: + def test_pooled_predictive_layout(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + forecast = pool.holdout_predictive + assert isinstance(forecast, ForecastResult) + values = forecast.idata.posterior_predictive["forecast"] + # n_draws defaults to the smallest member's flattened draw count: 2 x 100. + assert values.shape == (1, 200, 12, 2) + assert list(values.coords["variable"].values) == ["y1", "y2"] + pd.testing.assert_index_equal(pd.DatetimeIndex(values.coords["time"].values), pd.DatetimeIndex(holdout.index)) + assert forecast.median().shape == (12, 2) + assert forecast.hdi(0.89).lower.shape == (12, 2) + assert forecast.to_dataframe().shape == (12, 2) + assert forecast.mode == "density" + + def test_realised_weights_track_the_target(self): + fits = {"tight_y1": _fitted(sd=[0.3, 1.5]), "tight_y2": _fitted(sd=[1.5, 0.3])} + base = fits["tight_y1"] + mean = _mean_path(base, 12) + deviation = np.tile(np.array([[0.0, 2.5], [2.5, 0.0]]), (6, 1)) + pool = pool_forecasts(fits, _holdout(base, mean + deviation), n_draws=4000, seed=3) + realised = pool.realised_weights() + assert list(realised.index) == pool.labels + target = pool.weights.to_numpy() + tolerance = 4.0 * np.sqrt(target * (1.0 - target) / 4000) + assert np.all(np.abs(realised.to_numpy() - target) < np.maximum(tolerance, 1e-12)) + + def test_degenerate_weights_draw_from_one_model(self, dominance_pool): + fits, holdout = dominance_pool + pool = pool_forecasts(fits, holdout, method="log_score", seed=0) + assert pool.weights["good"] > 0.999 + winner = pool.labels.index(pool.weights.idxmax()) + assert set(np.unique(pool.membership)) == {winner} + model_coord = pool.holdout_predictive.idata.posterior_predictive["forecast"].coords["model"].values + assert set(model_coord) == {"good"} + + def test_model_coord_agrees_with_membership(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + coord = pool.holdout_predictive.idata.posterior_predictive["forecast"].coords["model"].values + np.testing.assert_array_equal(coord, np.array(pool.labels)[pool.membership]) + + def test_pooled_draws_are_member_draws(self, dominance_pool): + """Every pooled row is verbatim a draw from the model that produced it. + + Also pins the documented RNG contract: child generators are spawned in + `fits` insertion order, one per model, and the mixture takes the last. + """ + fits, holdout = dominance_pool + pool = pool_forecasts(fits, holdout, method="log_score", seed=0) + assert set(np.unique(pool.membership)) == {pool.labels.index("good")} + + children = np.random.default_rng(0).spawn(len(fits) + 1) + winner = fits["good"].forecast(steps=pool.steps, seed=children[pool.labels.index("good")]) + flat = winner.idata.posterior_predictive["forecast"].values.reshape(-1, pool.steps * 2) + pooled = pool.holdout_predictive.idata.posterior_predictive["forecast"].values[0] + known = {row.tobytes() for row in np.ascontiguousarray(flat)} + assert all(row.tobytes() in known for row in np.ascontiguousarray(pooled.reshape(-1, pool.steps * 2))) + + def test_n_draws_override(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, n_draws=37, seed=0) + assert pool.membership.shape == (37,) + assert pool.holdout_predictive.idata.posterior_predictive["forecast"].shape[1] == 37 + + +# -------------------------------------------------------------------------- +# H. Determinism +# -------------------------------------------------------------------------- + + +class TestDeterminism: + def test_same_seed_is_bit_identical(self, mirrored_pool): + fits, holdout = mirrored_pool + a = pool_forecasts(fits, holdout, seed=11) + b = pool_forecasts(fits, holdout, seed=11) + pd.testing.assert_series_equal(a.weights, b.weights) + pd.testing.assert_frame_equal(a.log_scores, b.log_scores) + np.testing.assert_array_equal(a.membership, b.membership) + np.testing.assert_array_equal( + a.holdout_predictive.idata.posterior_predictive["forecast"].values, + b.holdout_predictive.idata.posterior_predictive["forecast"].values, + ) + + def test_generator_seed_is_reproducible(self, mirrored_pool): + fits, holdout = mirrored_pool + a = pool_forecasts(fits, holdout, seed=np.random.default_rng(5)) + b = pool_forecasts(fits, holdout, seed=np.random.default_rng(5)) + pd.testing.assert_frame_equal(a.log_scores, b.log_scores) + np.testing.assert_array_equal(a.membership, b.membership) + + def test_different_seed_changes_the_pooled_draws(self, mirrored_pool): + fits, holdout = mirrored_pool + a = pool_forecasts(fits, holdout, seed=1) + b = pool_forecasts(fits, holdout, seed=2) + assert not np.array_equal( + a.holdout_predictive.idata.posterior_predictive["forecast"].values, + b.holdout_predictive.idata.posterior_predictive["forecast"].values, + ) + + +# -------------------------------------------------------------------------- +# I. combine() +# -------------------------------------------------------------------------- + + +class TestCombine: + def test_applies_weights_to_new_forecasts(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + refits = {label: fit.forecast(steps=20, seed=7) for label, fit in fits.items()} + combined = pool.combine(refits, seed=1) + assert combined.steps == 20 + assert combined.var_names == pool.var_names + values = combined.idata.posterior_predictive["forecast"] + assert values.shape == (1, 200, 20, 2) + assert "time" not in values.coords + assert set(values.coords["model"].values) <= set(pool.labels) + + def test_combine_respects_n_draws(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + refits = {label: fit.forecast(steps=4, seed=7) for label, fit in fits.items()} + assert pool.combine(refits, n_draws=13, seed=1).idata.posterior_predictive["forecast"].shape[1] == 13 + + def test_combine_is_deterministic(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + refits = {label: fit.forecast(steps=4, seed=7) for label, fit in fits.items()} + a = pool.combine(refits, seed=99).idata.posterior_predictive["forecast"].values + b = pool.combine(refits, seed=99).idata.posterior_predictive["forecast"].values + np.testing.assert_array_equal(a, b) + + def test_label_mismatch_rejected(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + refits = {label: fit.forecast(steps=4, seed=7) for label, fit in fits.items()} + refits["extra"] = next(iter(fits.values())).forecast(steps=4, seed=7) + with pytest.raises(ValueError, match="one forecast per pooled model"): + pool.combine(refits) + + def test_variable_mismatch_rejected(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + refits = {label: fit.forecast(steps=4, seed=7) for label, fit in fits.items()} + first = pool.labels[0] + refits[first] = refits[first].model_copy(update={"var_names": ["y2", "y1"]}) + with pytest.raises(ValueError, match="forecasts variables"): + pool.combine(refits) + + def test_step_mismatch_rejected(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + refits = { + label: fit.forecast(steps=steps, seed=7) for steps, (label, fit) in zip([4, 6], fits.items(), strict=True) + } + with pytest.raises(ValueError, match="same number of steps"): + pool.combine(refits) + + def test_mean_mode_rejected(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + refits = {label: fit.forecast(steps=4, include_shock_uncertainty=False) for label, fit in fits.items()} + with pytest.raises(ValueError, match="include_shock_uncertainty=True"): + pool.combine(refits) + + def test_combine_rejects_non_positive_n_draws(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + refits = {label: fit.forecast(steps=4, seed=7) for label, fit in fits.items()} + with pytest.raises(ValueError, match="n_draws must be at least 1"): + pool.combine(refits, n_draws=0) + + def test_combine_rejects_bad_seed(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + refits = {label: fit.forecast(steps=4, seed=7) for label, fit in fits.items()} + with pytest.raises(ValueError, match="seed must be"): + pool.combine(refits, seed="tomorrow") + + +# -------------------------------------------------------------------------- +# J. Frozen contract and plotting +# -------------------------------------------------------------------------- + + +class TestFrozenContract: + def test_attributes_cannot_be_reassigned(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + with pytest.raises(ValidationError): + pool.method = "log_score" + + def test_membership_is_read_only(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + with pytest.raises(ValueError, match="read-only"): + pool.membership[0] = 0 + + def test_weights_must_sum_to_one(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + broken = dict(pool) + broken["weights"] = pool.weights * 2.0 + with pytest.raises(ValidationError, match="sum to 1"): + PredictivePool(**broken) + + def test_weights_must_be_non_negative(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + broken = dict(pool) + broken["weights"] = pd.Series([1.5, -0.5], index=pool.labels) + with pytest.raises(ValidationError, match="non-negative"): + PredictivePool(**broken) + + def test_weight_labels_must_match_the_score_columns(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + broken = dict(pool) + broken["weights"] = pd.Series(pool.weights.to_numpy(), index=["p", "q"]) + with pytest.raises(ValidationError, match="log-score columns"): + PredictivePool(**broken) + + def test_membership_must_index_a_model(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + broken = dict(pool) + broken["membership"] = np.array([0, 1, 7]) + with pytest.raises(ValidationError, match="membership"): + PredictivePool(**broken) + + def test_membership_must_be_one_dimensional(self, mirrored_pool): + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + broken = dict(pool) + broken["membership"] = np.zeros((3, 2), dtype=int) + with pytest.raises(ValidationError, match="must be 1-D"): + PredictivePool(**broken) + + +class TestPlot: + def test_plot_returns_a_figure(self, mirrored_pool): + from matplotlib.figure import Figure + + fits, holdout = mirrored_pool + pool = pool_forecasts(fits, holdout, seed=0) + fig = pool.plot() + assert isinstance(fig, Figure) + assert len(fig.axes) == 1 + labels = [t.get_text() for t in fig.axes[0].get_yticklabels()] + assert set(labels) == set(pool.labels) + assert "stacking" in fig.axes[0].get_title() diff --git a/tests/test_public_api.py b/tests/test_public_api.py index 7a5af57..394fdb8 100644 --- a/tests/test_public_api.py +++ b/tests/test_public_api.py @@ -197,3 +197,34 @@ def test_long_run_restriction_in_all(self): import impulso assert "LongRunRestriction" in impulso.__all__ + + +class TestPoolingPublicAPI: + def test_predictive_pool_importable_from_impulso(self): + from impulso import PredictivePool + from impulso.pooling import PredictivePool as DirectPredictivePool + + assert PredictivePool is DirectPredictivePool + + def test_pool_forecasts_importable_from_impulso(self): + from impulso import pool_forecasts + from impulso.pooling import pool_forecasts as direct_pool_forecasts + + assert pool_forecasts is direct_pool_forecasts + + def test_pooling_names_in_all(self): + import impulso + + assert "PredictivePool" in impulso.__all__ + assert "pool_forecasts" in impulso.__all__ + + def test_every_exported_name_resolves(self): + import impulso + + for name in impulso.__all__: + assert getattr(impulso, name) is not None + + def test_plot_pool_weights_exported(self): + import impulso.plotting + + assert "plot_pool_weights" in impulso.plotting.__all__