6. Solution accuracy

A successful solve (retcode == 0) only tells you that the chosen method’s own residual is small at its anchor. It does not tell you how good the policy is as a global approximation to the true equilibrium. The standard, method-agnostic check is the Euler-equation error: substitute the candidate policy into the model’s optimality conditions on a set of states, integrate the conditional expectation over next-period shocks, and report the residual.

RISE measures this through a single, paradigm-free engine. Whatever solution method produced the policy — perturbation, Taylor projection, occasionally binding constraints, restricted log-linear, … — accuracy is evaluated through the same one-step forecast engine that simulate, irf and forecast use. The accuracy routine never asks how a forecast is produced.

6.1. The accuracy method

ee = accuracy(m)                 % at the deployed solve order
ee = accuracy(m, order)          % report orders 1..order
ee = accuracy(m, order, nxcuts, nsims)

Arguments:

  • order (default solve_order): the largest approximation order to assess. Accuracy is reported for orders 1..max(order). Each order is measured on the solution actually deployed at that order (the model is re-solved per order and assessed through the one-step engine), so the number reflects what you would simulate with. Under a paradigm only the deployed order is assessed.

  • nxcuts (default 10): number of Gauss–Hermite quadrature nodes per shock (made odd so that 0 is a node).

  • nsims (default 1000): number of ergodic draws of the state at which the Euler error is evaluated.

The result is a struct with one field per equation (eqtn_1, eqtn_2, …), each carrying one field per regime (regime_1, …), with

  • mean, min, max — the raw |residual| aggregated across the ergodic draws, per order;

  • log10_mean, log10_max — their base-10 logs, as the projection literature usually reports.

ee = accuracy(m, 5);
ee.eqtn_1.regime_1.log10_max     % log10 max Euler error, orders 1..5

Note

The errors are raw |residuals|, not normalized; take logs or any other transformation as desired. The smaller the standard deviation of the shocks, the higher the apparent accuracy — even if the solution is not accurate.

6.2. Visualizing where a solution is accurate

accuracy_plot evaluates the Euler error of the deployed solution on a regular grid of one or two endogenous variables (the others held at the steady state, shocks set to zero) and plots log10 of the largest residual across equations.

accuracy_plot(m, 'K')                 % line: log10|error| vs capital
accuracy_plot(m, {'K','A'})           % heatmap over (K, A)
[data, ax] = accuracy_plot(m, 'K', 'spread', 0.30, 'ngrid', 25, 'nxcuts', 7);

The same exhibit is discoverable under the plotting package as rise.plot.accuracy(m, ...), a thin alias of the accuracy_plot method.

For a perturbation solution the profile is a deep, narrow well: essentially exact at the steady state (its expansion point) and degrading away from it. This is the standard exhibit for “where in the state space is the solution reliable”.

6.3. Where to evaluate: ergodic cloud versus grid

accuracy and accuracy_plot compute the same Euler residual and integrate the conditional expectation with the same Gauss–Hermite quadrature over next-period shocks. They differ only in where in the state space the residual is evaluated.

  • accuracy draws the state from an ergodic simulation of the model (nsims draws from its stationary distribution). The reported error is weighted by how often the model actually visits each region and respects the co-movement among states (you never evaluate at an economically impossible combination), and it scales painlessly to high-dimensional state spaces. The cost: it under-samples the rare, extreme regions — deep recessions, near a bound, disaster tails — which is exactly where a low-order approximation is worst, so the verdict can flatter the solution. It is also a Monte-Carlo estimate (seed- and length-dependent).

  • accuracy_plot evaluates on a deterministic grid of one or two variables around the steady state, the rest held fixed. It is reproducible and lets you deliberately probe regions the simulation rarely reaches — push a state \(\pm x\%\) and watch the error grow. The cost: the curse of dimensionality (you grid one or two variables and freeze the rest, ignoring co-movement), grid points can be economically inconsistent, and the multiplicative spread anchor*(1 +/- spread) is awkward for variables whose steady state is zero or sign-indefinite.

Use the ergodic measure for “is this solution good enough for what I simulate?”; use the grid measure to see where and how fast the solution degrades away from the steady state. They are complementary.

Note

nxcuts builds a full Gauss–Hermite tensor over all shocks (\(\texttt{nxcuts}^{\,n_{\text{shocks}}}\) nodes). For models with many shocks this is infeasible (e.g. \(7^{11}\approx 2\times10^{9}\)); lower nxcuts, or integrate over only the shocks that matter.

6.4. Pruning and the accuracy verdict

At approximation order \(\geq 2\) the two measures part ways on pruning (see Forecasting and simulation).

  • The grid measure is pruning-free. Each evaluation applies the one-step forecaster at most two steps deep (\(x_{t-1}\to x_t\to x_{t+1}\)) at fixed, externally supplied states; nothing is iterated forward, so nothing compounds and nothing explodes. You get the pointwise accuracy of the order-\(k\) policy with no stand on pruning.

  • The ergodic measure is pruning-conditional. It simulates the model for nsims periods to build the cloud, and an unpruned order-\(\geq 2\) simulation can diverge. So the ergodic verdict is implicitly conditioned on the pruning scheme in force.

To compare high-order solutions (across orders, or perturbation versus Taylor projection) without taking a position on pruning, prefer the grid measure. Use the ergodic measure when you specifically want accuracy over the region the (pruned) model actually visits, and state that the number is conditional on the pruning choice.

6.5. How it works (paradigm-free)

The engine touches exactly two black boxes, neither of which it questions:

  • the one-step forecaster [x, y] = ff(y0, shks, r) (state, shocks, regime \(\to\) next; the forecast is the first output), built by set_model and shared with the simulators;

  • the structural residual f (the model’s equations, probs_times_dynamic, which already weights the future regime by its transition probability).

For each state \(x_{t-1}\) and current regime \(r\),

\[\text{resid}(x_{t-1}) = \sum_k w_k \sum_j f\bigl(r, j,\; x_{t+1}, x_t, x_{t-1}, e_t\bigr), \qquad x = \texttt{ff}(\cdot),\]

where the inner sums are Gauss–Hermite over next-period shocks and over future regimes \(j\). Both ingredients are paradigm-free, so one routine serves every solution method with no branch on solution type. The core lives in rise.engine.dsge_tools.accuracy (one_step_forecaster, euler_residuals, structural_residual, grid_states, shock_quadrature).

6.6. Perturbation versus Taylor projection

Solution accuracy is the natural way to compare solution methods. Two facts about \(k\)-th order Taylor projection (the rtp paradigm) come straight out of it:

  • At the steady-state anchor, single-anchor TP is perturbation. Both impose the same conditions (the residual and its derivatives up to order \(K\) vanish at the anchor), and RISE’s rtp reuses the perturbation machinery, so the Euler-error profiles coincide. There is no accuracy gain from a fixed, steady-state-anchored TP.

  • Re-anchored TP beats perturbation away from the steady state. With reanchor=true (the default, see Extending RISE through paradigms) the one-step re-solves the projection at the current state each step, re-centering the approximation where the model actually is. On the RBC model at order 3 the re-anchored Euler error is 8–13 times smaller than perturbation at \(\pm 30\%\) capital, while perturbation is more accurate only in a narrow band at the steady state. This is the reason to use Taylor projection.