15. Failure diagnosis: diagnose(m) and friends
Models fail: a steady state that will not solve, a solution that does not exist or is not stable, a simulation that quietly turns into NaNs, a perfect-foresight path the solver cannot reach. RISE’s failure-diagnosis suite turns each of those from a terse return code into a diagnosis: what check failed, observed versus required, who is implicated, and the next command to run. Two principles run through it:
no silent failures — anything that used to fail quietly now says where and why it died (governable, see the option below);
the right vocabulary per model class — eigenvalue (Blanchard-Kahn) counting is only ever used for constant-parameter models; for regime-switching models stability is judged by mean-square stability (MSS), and the diagnosis speaks that language exclusively.
15.1. One entry point
diagnose(m) % prints the diagnosis for whatever fails first
report = diagnose(m) % same, as a struct
diagnose checks the steady state first, then the solution. What you
get, by failure class:
Steady state fails. The worst equations ranked against the engine’s acceptance floor, labeled by their tag (or by the equation text itself for untagged models), the variables and parameters entering them, and a structural-rank check that names declared-but-never-used variables outright. When most equations are above the floor, the report says the failure is widespread — one global culprit (a bad discount/depreciation /growth number, a bad starting point) is more likely than the top-ranked equations — instead of overattributing.
No solution / unstable solution. For a constant-parameter model: root counting, with the root table (the failed solve retains its eigenvalues), the stable-versus-predetermined comparison, and — for the indeterminate case — a clearly-labeled heuristic (in New-Keynesian models, check the Taylor principle). For a regime-switching model: the MSS spectral radius against the criterion, each regime’s conditional spectral radius and persistence, and the culprit regime:
No mean-square-stable solution: ... The MSS spectral radius is 1.666,
above the criterion 1 ...
regime 1: conditional radius 0.9402, persistence 0.95
regime 2: conditional radius 1.322, persistence 0.95
Regime(s) 2 are conditionally explosive and, with the persistence
above, the process spends long enough there for second moments to
explode on average.
Note what that example teaches: a regime may be conditionally explosive and the model still be MSS — the mix is what matters, and only the MSS radius decides.
Simulation diverges. No more NaN-filled databases without a word: the failure site reports the first non-finite value (period, variable, regime), the growth profile of the run-up, and the pruning hint when simulating above order 1 unpruned.
Perfect foresight does not converge. The stacked residual at the solver’s final iterate is reshaped to (equation x period) and the worst cells are reported, plus the periods carrying 80% of the residual mass — “equation 3 around periods 12–15” localizes what “did not converge” cannot.
15.2. Fix-actions
Diagnosis ends with something to do, and two of the actions are automated:
Steady-state homotopy. Every successful steady-state solve remembers its calibration. When a later calibration fails:
[report, mfixed] = diagnose(m, 'fix', "homotopy");
walks the parameters gradually from the last working calibration to the
failing one, warm-starting each step from the previous solution. On
success, mfixed comes back solved at the target. On a stall, the
breakpoint is the diagnosis: “solvable up to 47% of the way — the
target calibration likely admits no steady state”, with mfixed solved
at the breakpoint mix so you can inspect how the steady state degenerates.
Perfect-foresight shock continuation. With simul_homotopy=true, a
failed stacked solve is retried by walking the shock size up from zero,
each step warm-started from the previous path. A stall short of full size
is reported as “solvable up to X% of the shock size” — the shock, not the
solver, is the problem.
15.3. The volume knob
m = set(m, 'diagnostics_on_failure', 'warn'); % default
% 'error' % throw instead
% 'off' % legacy silence
This gates every failure-site diagnostic (simulation divergence, PF worst cells, continuation summaries). Estimation is untouched by design — the likelihood loop stays silent at any setting.
15.4. The mean-square stability check and its algorithms
A regime-switching first-order solution x_t = T(s_t) x_{t-1} is
mean-square stable when the operator on the regime-indexed second moments,
X_j <- T_j (sum_i Q(i,j) X_i) T_j' with Q(i,j) the probability of
moving from regime i to regime j, has spectral radius below one (Costa,
Fragoso and Marques 2005, Theorem 3.9). RISE offers three ways of
computing that radius, selected by stability_algorithm:
'gmh'(the default) and'cfm'build the matrix of the operator,(Q' kron I) diag(kron(T_i,T_i)), of dimension h n^2 for n states and h regimes, and take its eigenvalues. The cost grows with the sixth power of the number of states: 280 seconds at 80 states with two regimes.'matrix_free'applies the operator itself inside an Arnoldi iteration, started in the cone of positive semidefinite matrices where the spectral radius is an eigenvalue. It costs 2 h n^3 operations per application: two seconds at 500 states with two regimes, where the matrix would take two terabytes. Use it on large models instead of switching the check off (solve_check_stability = false), and under estimation, where the Kronecker algorithms are skipped above a size threshold and this one is not. The self-consistent linearization ({'m','scl'}) checks the first-order block with the same algorithm before it computes the higher orders: at 40 states and two regimes an order-2 solve takes about 2 seconds under'matrix_free'and 15 to 19 under'gmh'.
Until September 2026 the two Kronecker algorithms placed Q(i,j) next to
the block of regime j, which gives an operator similar to that of the
time-reversed chain (whose transition matrix is diag(pi)^-1 Q' diag(pi),
pi the stationary distribution, not Q').
The radius is the same when the chain is reversible, which every two-state
chain and every product of independent two-state chains is; on a
non-reversible chain with three or more regimes it differs, and the
verdict could be wrong in either direction. Both routines now build the
operator of the chain itself, and diagnose reports its radius.
15.5. Is the stable solution the only one? Determinacy
Mean-square stability says that the solution solve returned has
bounded second moments. It does not say that it is the only such
solution: a regime-switching model can have several MSS minimal-state
solutions, and sunspot solutions with bounded second moments. The
determinacy check answers that question for the first-order solution:
[ms,rc,sm] = solve(m);
d = rise.engine.dsge_tools.stability.determinacy(ms, sm)
Write the first-order system of regime i as
sum_j A+_ij E[x(t+1) | s(t+1)=j] + A0_i x(t) + A-_i x(t-1) = 0, the
blocks carrying the transition probabilities p_ij, and let T be
the returned solution. Any other solution differs from it by a process
y(t) = sum_j G_ij E[y(t+1) | s(t+1)=j] with
G_ij = -(A0_i + sum_k A+_ik T_k)^{-1} A+_ij, driven by its
forward-looking block. d.no_bubble_radius is the spectral radius r
of the operator U_i <- sum_j (1/p_ij) G_ij U_j G_ij' in the chain’s own
direction of time. The verdict d.verdict is
'determinate'whenris below one: every other candidate has unbounded second moments, and the returned solution is the unique MSS solution;'indeterminate'whenris above one: an MSS sunspot solution exists, and the check builds it (below);'not MSS'when the returned solution is not itself strictly MSS, so the question does not arise;'unresolved'when a radius is withinMargin(default1e-9) of one, or when a check of the input fails;d.messagesays which. At a radius of exactly one the answer depends on the class of solutions admitted (bounded second moments, or second moments that converge), and floating-point numbers cannot place a radius that close to one on either side. The checks, ind.checks: the solution and the derivatives are finite,A0_i + sum_j A+_ij T_jis invertible (rcond), the solution solves the first-order system of the derivatives (solution_residual), and the Arnoldi iterations converged. A solution with a unit root (a balanced-growth model solved in levels, for instance) is'unresolved'too, but whenris below one the message says what remains true: no bubble has second moments that grow less than exponentially, so the solution is the only one among those whose second moments grow subexponentially.
The probabilities p_ij are those the derivative blocks carry. They
cannot be read off the blocks, because the same equations can be spread
over the regime pairs in more than one way, so the check builds them from
the model: the true chain, except under the Maih-Waggoner perturbation,
which expands at sigma = 0 with the chains it perturbs frozen (their
transition matrices become the identity). d.reference says which law
was used, d.reference_law holds it, and the option 'ReferenceLaw'
supplies another. Under Maih-Waggoner the verdict therefore concerns the
economy in which the perturbed chains never switch; the stability of the
solution simulated under the true chain is a separate question, reported
separately: d.mss_radius is the second-moment radius under the true
chain, d.reference_mss_radius the one under the reference law. The
verdict is uniform over initial regimes. When the reference chain is
reducible (a frozen chain, for instance) d.no_bubble_by_regime gives
the radius of the regimes reachable from each initial regime.
Above one the check does not stop at a failed condition: it builds an MSS
sunspot solution and returns it in d.sunspot. With U the Perron
eigenvector of the operator (positive semidefinite, L(U) = r U),
factored as U_i = W_i W_i', and F_ij = G_ij / p_ij on the forward
rows, the sunspot is y(t) = W_s(t) z(t) with
z(t+1) = K_ij z(t) + e(t+1), e a martingale difference given the
current information and the next regime. The transitions K_ij solve
the full forward equations sum_j p_ij F_ij W_j K_ij = W_i in minimum
norm, and alpha = max_i lambda_max(sum_j p_ij K_ij' K_ij) bounds the
growth of the sunspot’s second moments: 1/r in exact arithmetic, the
smallest rate a sunspot can have. The rank of each W_i is a numerical
decision: a computed U carries noise eigenvalues, and a noise direction
fails the full forward equations. The check tries eigenvalue thresholds of
1e-10, 1e-8 and 1e-6 relative to the largest and keeps the
first support that passes (support_threshold); when the Perron
eigenvector of the chosen Method does not pass, the Arnoldi
iteration at a tolerance of 1e-13, then the Kronecker route, are
tried. A computed K leaves a residual E_i in the forward
equations. Where C_i = [sqrt(p_ij) F_ij W_j]_j has full row rank, a
right inverse corrects it exactly at a distance eta it bounds, and the
test budgets that correction against the margin: (|X_i| + eta)^2 < 1,
with X_i = [sqrt(p_ij) K_ij]_j. On a singular support (C_i rank
deficient, the common case) no correction can put W_i in the range of
C_i, and the residual is a numerical check only. The fields are
verified (the forward equations hold to 1e-8 and the budgeted
contraction is below one: a numerical verification in floating point,
not a validated proof), contraction (alpha),
contraction_bound (alpha with the correction budgeted),
singular_support, residual (the largest relative error in the
forward equations), bound (1/r), radius (the second-moment
radius of the process built, computed on z: at most alpha),
support_threshold, loadings
(B_ij = W_j K_ij W_i^+ = U_j F_ij' U_i^+ / r, in the coordinates of the
forward variables), variables (the forward-looking variables that
carry it; auxiliary variables pinned by a static equation carry none),
dimension (the rank of W_i in each regime: one on a
constant-parameter model, the unstable direction), regimes (where it
lives) and initial_regimes (the regimes from which those can be
reached). On a constant-parameter model the check reduces to the
Blanchard-Kahn condition; on the Farmer, Waggoner and Zha (2011) example
with a unique minimal-state solution but a passive policy in both regimes
it reports 'indeterminate' with a verified sunspot: a unique
minimal-state solution is not determinacy.
The check is Cho’s (2016) no-bubble condition with the orientation of the
chain made explicit (it matters on non-reversible chains with three or
more regimes), and it applies to every perturbation type, including M
with the self-consistent linearization, where it concerns the first-order
system at the rest configuration returned. Large models use the
matrix-free route ('Method','eigs'), as the stability check does.
The default solver can miss the stable solution. Functional iteration from
the backward start (mfi, which is Cho’s forward method) is attracted by
the unique MSS solution of a determinate model, but only locally: from its
start it can converge to a solution that is not MSS, or not converge at
all, and Newton from the default start can do the same. solve then
reports no MSS solution for a model that has exactly one. The search
function tries Newton from random starting points when that happens and
checks what it finds:
[ms,d,info] = rise.engine.dsge_tools.stability.find_stable_solution(m)
[ms,d,info] = rise.engine.dsge_tools.stability.find_stable_solution(m, ...
'Starts',50,'SolveOptions',{'solve_perturbation_type',{'m','scl'}})
d.verdict = 'determinate' means the solution found is the unique MSS
solution, whichever start found it; 'KeepSearching', true tries every
start and reports every distinct MSS solution (two or more prove
indeterminacy). When no start returns an MSS solution the function says
so, which is not a proof that none exists.
15.6. Beyond mean-square stability: higher moments
Mean-square stability bounds the second moments of the first-order
solution. A pruned simulation of order k is driven by k-fold
products of the first-order state, so its second moments need the
2k-th moments of that state, and mean-square stability does not
deliver them: a regime that is explosive but short-lived can leave the
second moments bounded and the fourth moments unbounded (with
x(t) = 1.15 x(t-1) in a regime left with probability 0.4 and
0.3 x(t-1) in the other, left with probability 0.1, the second-moment
radius is 0.80 and the fourth-moment radius 1.05). The check:
[ms,rc,sm] = solve(m);
r = rise.engine.dsge_tools.stability.moment_stability(ms, sm) % fourth moments
r = rise.engine.dsge_tools.stability.moment_stability(ms, sm, 'Order', 6)
reports the spectral radius of the operator of the moments of the
requested order (the mean-square-stability operator of x Kronecker
x) next to the mean-square-stability radius. By Jensen’s inequality
the fourth-moment radius is at least the square of the second-moment
radius. On the published models we ran it is far below one (0.22 to
0.86). The fourth-moment operator acts on n^2-by-n^2 matrices; the
default route applies the state block to each index of the moment tensor
(about 3 seconds at 20 states, a minute at 40, four to five minutes at 60,
with two regimes), so it is a diagnostic for models with a few dozen
states.