Solving ======== ``solve`` is the gateway between a parsed model and everything else. You pass it a model object and a handful of options; you get back a solved model and a return code. This chapter is the reference for the option surface, the return-code discipline, and the diagnostic protocol when a solve fails. .. contents:: :local: :depth: 2 The signature -------------- :: [m, retcode] = solve(m, name = value, name = value, ...); The first input is a model object. Subsequent arguments are name-value pairs validated by an ``arguments`` block on ``solve``. The first output is the model object with the solution attached; downstream methods (``irf``, ``simulate``, ``estimate``, ``filter``, ...) consume it. The second output is the **retcode** -- always capture it and always check it. The retcode pattern -------------------- The retcode is the single source of truth about whether the solve succeeded:: [m, retcode] = solve(m); assert(retcode == 0, decipher(retcode)); * ``retcode == 0`` -- the model is solved. * ``retcode > 0`` -- the solve failed in a specific, named way. * ``decipher(retcode)`` returns the human-readable explanation. Never read the model object's solution fields without first checking ``retcode``. A non-zero retcode can leave the solution empty or partially populated; downstream methods will then return empty results or raise unrelated errors. A complete list of return codes appears at the end of this chapter. Picking a solver ----------------- The ``solver`` option selects the underlying algorithm. For constant-parameter models RISE picks a sensible default. For switching models it defaults to ``mfi`` (functional iteration). You may force a specific solver by prefixing the name with ``+``:: m = solve(m, solver = '+mn'); % force Newton, no auto-pick The available solvers: .. list-table:: :header-rows: 1 :widths: 20 25 55 * - Name - Family - Notes * - ``mfi`` - Functional iteration (switching default) - Cheapest per iteration. Can miss solutions where the operator is not a contraction; see the diagnostic protocol below. * - ``mn`` - Newton (switching), no Kronecker - Robust to non-contraction. Use when ``mfi`` returns retcode 21. * - ``mnk`` - Newton (switching), with Kronecker - Same robustness as ``mn`` with a different inner solve; sometimes faster on larger systems. * - ``mfi_full``, ``mn_full``, ``mnk_full`` - Full (non-sparse) variants - Use when sparsity exploitation in the standard variants causes numerical issues. * - ``fwz`` - Farmer-Waggoner-Zha Newton - Their algorithm; sometimes finds solutions the others miss, sometimes the reverse. * - ``dsge_groebner`` - Groebner-basis all-solutions - Finds every solution. Can be very slow (minutes to never) on non-trivial models. Do not include in a default diagnostic sweep. * - ``dsge_udc`` - Undetermined coefficients - For switching models. * - ``dsge_schur`` - Schur decomposition - Only with the Maih-Waggoner perturbation strategy. * - ``rise_1`` - Constant-parameter default - Generalised Schur / QZ decomposition. * - ``klein`` - Klein's method - Constant-parameter; classical reference algorithm. * - ``cr``, ``cyclic_reduction`` - Cyclic reduction - Constant-parameter. * - ``AIM`` - AIM solver - Constant-parameter. A *user-defined* solver -- a function handle with a specific signature -- is also accepted; see the ``solve`` doc-comment for the required inputs and outputs. Solver options --------------- ``solver`` may be a plain string or a two-element cell ``{name, opts}`` where ``opts`` is a struct that RISE passes through to the inner iterator. The fields use MATLAB-standard option names:: opts = struct( ... MaxIter = 1e6, ... TolFun = 1e-10, ... Display = 'iter', ... ObjectiveLimit = 1e+12); [m, retcode] = solve(m, solver = {'mn', opts}); .. list-table:: :header-rows: 1 :widths: 25 20 55 * - Field - Default - Role * - ``MaxIter`` - ``1000`` - Maximum number of iterations. * - ``TolFun`` - ``sqrt(eps)`` - Convergence tolerance on the residual. * - ``Display`` - ``'off'`` - Set to ``'iter'`` to print convergence per iteration. * - ``ObjectiveLimit`` - ``1e+12`` - Treat the iteration as divergent if the residual exceeds this. .. note:: The legacy toolbox accepted the same options under the names ``fix_point_maxiter``, ``fix_point_TolFun``, ``fix_point_verbose``, ``fix_point_explosion_limit``. The modern convention is the MATLAB-standard names above; the legacy names have been retired as user-facing options and passing them now errors as an unknown option. Order of approximation ----------------------- ``solve_order`` controls the order of the perturbation expansion:: m = solve(m, solve_order = 3); Supported orders are ``1`` through ``5``. Orders above ``1`` require that the model was *parsed* with at least the matching ``max_deriv_order`` -- this is a parser-time choice:: m = dsge_model('mymodel.rs', max_deriv_order = 3); m = solve(m, solve_order = 3); If you parsed at a lower order than you want to solve, re-parse the model -- you cannot upgrade after the fact. At order :math:`\geq 2` a simulated path can diverge. Pruning is an explicit opt-in for keeping such paths bounded, off by default and enabled through the ``simul_pruned`` option; see :doc:`../../WorkingWithAModel/Forecasting and simulation`. Optimal policy --------------- Two orthogonal options control optimal-policy solutions: .. list-table:: :header-rows: 1 :widths: 35 40 25 * - Option - Values - Default * - ``solve_policy_type`` - ``'ramsey'``, ``'discretion'`` (case-insensitive) - ``'ramsey'`` * - ``solve_policy_equilibrium`` - ``'OLE'``, ``'MPE'`` (case-insensitive) - ``'OLE'`` ``solve_policy_type`` controls commitment versus discretion: * ``'ramsey'`` -- full commitment; the planner re-optimizes only at ``t = 0``. * ``'discretion'`` -- time-consistent / Markov-perfect policy; the planner re-optimizes every period. ``solve_policy_equilibrium`` controls the equilibrium concept in multi-player games: * ``'OLE'`` (open-loop equilibrium) -- the classical Lagrangian formulation in which each player takes the opponent's *path* as given. * ``'MPE'`` (Markov-perfect equilibrium) -- each player takes the opponent's *policy function* as given. All four combinations are supported. Both options are passed to ``solve`` -- they are not declared in the model file. Loose commitment and stochastic replanning are also available through a small extension of this surface; see *Optimal policy* in the Modeling chapter for the model-side syntax. Occasionally-binding constraints --------------------------------- The piecewise-linear (OccBin) solve of an occasionally-binding constraint is provided by the ``occbin.paradigm`` solve paradigm, activated through the ``solve_paradigm`` option:: C = { struct('var','R','dir',-1,'bound',1,'chain','ocb') }; % R >= 1 m = set(m,'solve_paradigm',{@occbin.paradigm,'ref',1,'constraints',C}); m = solve(m); The model-side construction (the binding chain and the convex-combination equation) is covered in :doc:`Occasionally-binding constraints`, and the paradigm interface -- the constraint-cell grammar and its options -- in :doc:`/WorkingWithAModel/Extending RISE through paradigms`. Perturbation strategy ---------------------- ``solve_perturbation_type`` selects how the perturbation is taken around the deterministic steady state of a switching model. The strategies differ in which parameters and steady states are perturbed. .. list-table:: :header-rows: 1 :widths: 25 75 * - Value - Strategy * - ``'maih'`` or ``'m'`` (default) - Maih (2014): regime-specific steady states. * - ``'maih_waggoner'`` or ``'mw'`` - Maih and Waggoner (2018): regime-specific steady states plus perturbed transition probabilities on *all* chains. * - ``{'mw', mkv_list}`` - As above but only the chains in ``mkv_list`` have their probabilities perturbed. * - ``'frwz'`` - Foerster, Rubio-Ramirez, Waggoner and Zha (2016): unique steady state at the ergodic mean of the parameters, perturbing all switching parameters. * - ``{'frwz', part_list}`` - As above but only the parameters in ``part_list`` are perturbed. ``mkv_list`` and ``part_list`` may also be ``{'*none'}``, ``{'*all'}``, ``{'*all_but'}``, or ``{'*all_but', 'v1', 'v2', ...}``. The default ``'maih'`` is the right choice unless you have a specific reason to use one of the others -- see the FRWZ paper for when their strategy is preferred. .. note:: Under the default ``'maih'`` perturbation, RISE auto-declares several auxiliary parameters (``iota_frwz_*``, ``ss_frwz_*``, ``perturbator``, ...) that are not used in the model. Their absence from your calibration triggers a benign ``"1(incomplete!)"`` warning at parse time. This warning is *not* an error; do not chase it. Chase the retcode. Other options worth knowing ---------------------------- .. list-table:: :header-rows: 1 :widths: 30 70 * - Option - Role * - ``solve_linear`` - ``true`` to use more efficient steady-state algorithms when the model is truly linear. * - ``solve_bgp`` - ``true`` to solve the model as non-stationary along a balanced growth path. * - ``solve_derivatives_only`` - ``true`` to compute the symbolic derivatives but skip the solve. Useful for diagnostics. * - ``solve_kill_g`` - ``true`` to kill off the growth element in the first-order solution. * - ``solve_automatic_differentiator`` - Function handle for the automatic differentiator engine (default: ``@adolm.diff``). * - ``solve_user_defined_shocks`` - Struct describing non-Gaussian shock distributions. See the ``solve`` doc-comment. * - ``solve_qz_criterium`` - Eigenvalue cutoff for constant-parameter QZ-based solvers (default: ``sqrt(eps)``). The diagnostic protocol when solve fails ------------------------------------------ When ``solve`` returns a non-zero retcode, work down this list in order; each step is cheap and rules out a class of causes: 1. **Read the retcode.** ``decipher(retcode)`` gives the human-readable reason. If the message is specific (e.g. "Complex solution", "Cannot take a log-expansion of a variable whose steady state is close to 0"), fix the model, not the solver settings. 2. **If retcode is 21 ("Maximum Number of iterations reached or multiple solutions"), raise the iteration cap first.** It is the cheapest test and rules out simple truncation:: opts = struct(MaxIter = 1e8); [m, retcode] = solve(m, solver = {'mfi', opts}); If ``retcode == 0``, you were just being truncated. 3. **If retcode is still 21, switch to a Newton solver.** ``mfi``'s functional iteration is not a contraction on every model; Newton solvers are robust to that and find fixed points ``mfi`` cannot reach:: [m, retcode] = solve(m, solver = '+mn'); % or [m, retcode] = solve(m, solver = '+mnk'); 4. **If Newton returns retcode 25 ("System unstable"), the model is genuinely unstable** for this parameterization. Investigate the calibration -- vary the parameters and recheck the determinacy boundary. The solver is not the issue. 5. **If retcode is 28 ("transition matrix must be diagonal"), an OBC chain has non-zero off-diagonal transition probabilities.** Set them to zero -- OBC switching is endogenous, not stochastic. 6. **If retcode is in the evaluation range (1-11), the failure is upstream of the solver.** Read the message; it points at the steady-state file, the Jacobian, the transition matrix construction, or the parameter restrictions. Fix that before touching ``solve``. Return-code reference ---------------------- The canonical table lives in ``utils.error.registry``; this is the human-readable summary, in three groups. Evaluation (pre-solve) ~~~~~~~~~~~~~~~~~~~~~~~ .. list-table:: :header-rows: 1 :widths: 12 88 * - Code - Meaning * - ``1`` - Steady state could not solve * - ``11`` - Complex or nan steady-state residuals * - ``2`` - NaNs in Jacobian * - ``3`` - Problem in transition matrix * - ``4`` - Parameter restrictions violated * - ``5`` - Definitions are NaN, Inf, or imaginary * - ``6`` - NaNs or Infs in planner objective * - ``7`` - Bounds or restriction violation Solving ~~~~~~~~ .. list-table:: :header-rows: 1 :widths: 12 88 * - Code - Meaning * - ``21`` - Maximum number of iterations reached **or** multiple solutions * - ``210`` - Discretionary outer loop: max iterations reached * - ``22`` - NaNs in solution or no solution * - ``223`` - Complex solution * - ``23`` - Explosion limit reached * - ``230`` - Discretionary outer loop: explosion limit reached * - ``24`` - Explosive solution * - ``25`` - System unstable * - ``251`` - DSGE-VAR solution: system unstable * - ``26`` - VAR approximation to the DSGE failed * - ``261`` - DSGE prior is not proper * - ``262`` - DSGE-VAR solution: covariance matrix not positive definite * - ``263`` - SVAR solution: matrix not invertible * - ``27`` - Cannot take a log-expansion of a variable whose steady state is close to 0 * - ``28`` - For this solver the transition matrix must be diagonal Filtering and likelihood (downstream of solve, listed for completeness) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. list-table:: :header-rows: 1 :widths: 12 88 * - Code - Meaning * - ``301`` - Max iterations reached in Lyapunov solution * - ``302`` - NaNs in Lyapunov solution * - ``303`` - Explosion limit reached in Lyapunov solution * - ``304`` - Failure of covariance in sigma-points determination * - ``305`` - Covariance of forecast errors not positive definite * - ``306`` - Unlikely parameter vector * - ``307`` - NaN/Inf/complex in log prior * - ``308`` - Inconsistent ergodic probabilities * - ``309`` - Error in user-defined endogenous priors * - ``310`` - NaNs or Inf in initial conditions for the state vector * - ``311`` - NaNs or Inf in initial covariance matrix * - ``312`` - NaN/Inf/complex likelihood * - ``313`` - Forecast covariance positive definite but determinant = 0 The retcode-aware pattern -------------------------- The pattern is the same for every solve in a script:: [m, retcode] = solve(m, solver = '+mn'); if retcode ~= 0 warning('solve failed: %s -- trying mfi with raised MaxIter', ... decipher(retcode)); opts = struct(MaxIter = 1e8); [m, retcode] = solve(m, solver = {'mfi', opts}); end assert(retcode == 0, decipher(retcode)); print_solution(m); Always capture the retcode; always decipher on failure; do not read the solution fields until ``retcode == 0``. See *Teaser Example* for the canonical worked solve on the standard documentation model.