feat(pooling): static predictive-density pooling

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
 

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.

docs/how-to/index.md

@@ -13,4 +13,5 @@ climate-pitfalls
 sign-restrictions
 long-run-restrictions
 heavy-tailed-errors
+pooling
 ```

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.

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

docs/reference/index.md

@@ -18,6 +18,7 @@ fitted
 identified
 identification
 scenario
+pooling
 results
 evidence
 primitives

docs/reference/plotting.md

@@ -15,6 +15,7 @@
    plot_fevd
    plot_historical_decomposition
    plot_counterfactual
+   plot_pool_weights
    plot_volatility
    plot_sv_forecast
 ```

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
+```