8. Storytelling

The storytelling layer takes an estimated, data-fed model and answers the questions an empirical paper actually has to answer:

Which structural shocks were active in each period?
What was the most likely regime path – and how confident are we?
What would the path have looked like without the monetary shock? Without the fiscal regime change?
How much credit do we give each shock for each datapoint?

Three layers do the work, in order: the smoother (filter), the classical historical decomposition (historical_decomposition), and the modern resimulation surface (+resim_bridge).

This chapter is the modern entry-point reference. The legacy chapter Resimulation and counterfactuals remains canonical for the algorithm details and the per-method API and is scheduled to fold into this tree pending review.

8.1. `filter`

Everything downstream starts from the smoother. filter(m) runs the Kalman filter (or, for switching models, the regime-conditional Kalman filter; for nonlinear models, the unscented or particle filter). It returns the smoothed states, smoothed shocks, smoothed regime probabilities, and the filter likelihood:

[filt, LogLik, Incr, retcode, m] = filter(m_est);
assert(retcode == 0, decipher(retcode));

The five outputs:

Output	Meaning
`filt`	Struct of `ts` carrying the filtered, updated, and smoothed quantities (see below).
`LogLik`	Sample log-likelihood at the parameter vector currently bound to `m_est`.
`Incr`	Per-period likelihood increments, useful for diagnostics.
`retcode`	Return code; `0` on success. Use `decipher` on failure.
`m`	The model with the filter output attached. Required by `historical_decomposition` and `+resim_bridge`; capture this output.

The fields of filt follow a predictable naming pattern: each quantity (variables, shocks, measurement errors, regime probabilities, state probabilities) comes in three flavours – filtered (one-step prediction), updated (conditional on the contemporaneous observation), and smoothed (conditional on the whole sample) – plus an Expected_* companion that averages over regime uncertainty. The most-used entries:

Field	Content
`smoothed_variables`	Smoothed endogenous variables (struct of `ts`).
`smoothed_shocks`	Smoothed structural shocks (struct of `ts`).
`smoothed_regime_probabilities`	One `ts` per composite regime; `regime_1` through `regime_h`.
`smoothed_state_probabilities`	One `ts` per chain-state pair.
`filtered_variables`, `updated_variables`, `filtered_shocks`, `updated_shocks`	The one-step and contemporaneous counterparts.
`smoothed_measurement_errors`	When the model declares them.

For a model with a mon chain (active/passive monetary) and a fis chain (active/passive fiscal) the smoothed marginal probability of an active monetary regime is the sum of the composite-regime probabilities consistent with that marginal:

p_active_mon = filt.smoothed_regime_probabilities.regime_1 ...
             + filt.smoothed_regime_probabilities.regime_3;

(The composite-regime numbering follows m.markov_chains.regimes.)

8.2. Historical decomposition 

historical_decomposition splits each datapoint into the contributions of each structural shock plus an initial conditions term:

hd = historical_decomposition(m);

The output is a struct of ts keyed by endogenous-variable name. Each ts carries one column per contributing series; the varnames property lists them. The contributions sum to the realized path: for any variable y and any date t:

y_t  ==  sum_k hd.y(t, k)   over all columns k

The columns include every declared exogenous shock plus an init column that absorbs the effect of initial conditions and the steady state.

For linear constant-parameter models the decomposition is classical and exact. For switching or nonlinear models the decomposition is approximated by a Monte Carlo integration over the regime distribution; the integration size is controlled by the nsim argument to historical_decomposition (default 100 for switching, automatically 1 for linear constant- parameter).

Grouping shocks is supported via the fourth argument:

groups = struct( ...
    demand   = {'eps_g'}, ...
    monetary = {'eps_m'}, ...
    fiscal   = {'eps_f'});

hd = historical_decomposition(m, [], [], groups);

Each group becomes a single column in the output ts; unlisted shocks are pooled under others.

8.3. `+resim_bridge` – the resimulation surface 

The classical historical decomposition tells you who’s responsible for each datapoint. The +resim_bridge package answers the more specific questions: which regime path was the economy really on, what would have happened without a specific shock, what is the exactly-additive credit assignment across shocks.

Five reconstruction methods, three counterfactual modes, and three decomposition utilities sit on top of the same adapter object.

8.3.1. Adapter 

resim_bridge.adapt packages a solved-and-filtered model into the canonical tuple used by every downstream method:

adapter = resim_bridge.adapt(m, filt);

You rarely call this directly – run, counterfactual, shock_decomp, shapley, and var_decomp all call it internally. Direct access is useful when writing a custom reconstruction or counterfactual.

The adapter exposes Tfunc(x, eps, r) (the one-step state propagator), Pfunc(x) (the state-dependent regime transition matrix, constant for exogenous Markov chains), and the smoothed regime-conditional means and probabilities (smoothed.x, smoothed.eps, smoothed.pi).

8.3.2. Reconstruction 

resim_bridge.run(m, filt, method, opts) reconstructs the latent regime path consistent with the smoother output. Four methods:

`method`	What it returns
`'viterbi'`	The most-likely regime path (single deterministic sequence) and the state path obtained by forward-simulating along it.
`'ffbs'`	`M` draws from the posterior regime-path distribution via forward-filter backward-sample. The principled way to get joint uncertainty over the regime history.
`'independent'`	`M` regime paths drawn independently at each date from the smoothed marginals. Ignores the transition law; correct only marginally. Useful for fast exploratory analysis.
`'rao_blackwell'`	Deterministic `pi`-weighted mixture across regimes. Exact for linear models; biased for nonlinear ones with the bias growing in regime uncertainty.

The most common pairing is Viterbi for headline plots and FFBS for uncertainty bands:

viterbi = resim_bridge.run(m, filt, 'viterbi');
ffbs    = resim_bridge.run(m, filt, 'ffbs', struct(M = 1000));

8.3.3. Counterfactuals 

resim_bridge.counterfactual(m, filt, mode, spec) reconstructs a deterministic baseline (Viterbi by default) and re-simulates under one perturbation:

`mode`	`spec` fields	Effect
`'shock'`	`shock_ids`, optional `values`	Zeroes (or replaces) the listed shocks. Regime path held fixed.
`'regime'`	`r_path`	Replaces the reconstructed regime path with a user-supplied sequence. Shocks held fixed.
`'param'`	`parameters` struct	Re-sets parameters, re-solves, re-filters on the same data, and re-runs the reconstruction.

The return struct carries baseline, counter, and delta = counter.x - baseline.x. A what-if on monetary shocks looks like:

cf = resim_bridge.counterfactual(m, filt, 'shock', ...
                                 struct(shock_ids = {'eps_m'}));
plot(cf.baseline.pi, cf.counter.pi);
title('Inflation: actual vs counterfactual without monetary shocks');

8.3.4. Shock decompositions 

Two flavours, both along the reconstructed regime path.

resim_bridge.shock_decomp(m, filt) is the leave-one-out decomposition:

delta_k(t) = x_t^full - x_t^(-k)

It is exact for linear models; for nonlinear models the interaction residual is not allocated.

resim_bridge.shapley(m, filt) is the Shapley decomposition:

phi_k(t) = (1/ne!) sum_sigma [ v(S U {k}, t) - v(S, t) ]

The Shapley decomposition is exactly additive even for nonlinear models:

sum_k phi_k(t) == x_t^full - x_t^base

so that the residual is at machine precision (smoke runs report residual = 1.11e-16).

Exact mode enumerates all \(2^{n_e}\) subsets and runs for ne <= opts.exact_threshold (default 10). For ne above the threshold the sampled estimator of Castro, Gomez, and Tejada (2009) is used.

8.3.5. Ergodic variance decomposition 

resim_bridge.var_decomp(m, kind) works at the model’s ergodic distribution; it does not consume a filter output and is therefore a property of the solved model alone. Three games:

`kind`	Decomposition
`'shocks'`	Shapley allocation of total variance across the `ne` shocks. The regime process runs freely under `Pfunc`.
`'joint'`	Shapley allocation across the shocks and a `regime` player. Removing the regime player fixes the chain at its ergodic distribution.
`'anova'`	Between-vs-within decomposition via the law of total variance: `Var(x) = W + B` with `W = sum_j pi_j * Var(x \| r_t = j)` and `B = sum_j pi_j * (mu_j - mu)(mu_j - mu)'`. Always exact; no Shapley enumeration.

8.4. A worked story 

The empirical section of a paper, in five calls, on the estimated fiscal-monetary switching NK from Estimation:

% 1. filter
[filt, ~, ~, rc, m_est] = filter(m_est);
assert(rc == 0, decipher(rc));

% 2. who's responsible for each datapoint?
hd = historical_decomposition(m_est);

% 3. most-likely regime path
viterbi = resim_bridge.run(m_est, filt, 'viterbi');
plot(viterbi.r_path); ylim([0.5, m_est.markov_chains.regimes_number + 0.5]);
title('Most-likely composite-regime path');

% 4. what if there had been no monetary shocks?
cf = resim_bridge.counterfactual(m_est, filt, 'shock', ...
                                 struct(shock_ids = {'eps_m'}));
plot(cf.baseline.pi, cf.counter.pi);
legend('actual', 'no monetary shocks');

% 5. Shapley credit assignment, exactly additive
sh = resim_bridge.shapley(m_est, filt);

These five calls answer “what drove the data”, “what regime were we in”, “what would have happened without X”, and “how much credit does each shock get” – the standard menu of an empirical paper’s results section.

8.5. Where to look next 

Resimulation and counterfactuals (legacy DSGE chapter) – the long-form algorithm description plus per-method API references for resim_bridge.adapt, run, counterfactual, shock_decomp, shapley, var_decomp. Canonical; fold pending review.
Master filtering (legacy DSGE chapter) – the full surface of filter options, including the choice of nonlinear filter (regime-switching unscented, particle), the initialisation options, and the OBC variants.
Forecasting and simulation – simulate and forecast, used for forward extensions of a Viterbi baseline.