19. Extending RISE through paradigms

This chapter describes how RISE separates what a model is from how it is solved, through the notion of a paradigm. A paradigm is a self-contained computational strategy – perturbation, recursive local linearization (RLL), \(k\)-th order Taylor projection (TP), occasionally-binding-constraint (occbin) piecewise linearization, and so on – that plugs into mainland RISE through a small, fixed interface. The aim is to add advanced solution methods without modifying RISE’s core and without paradigms interfering with one another.

This chapter is primarily for developers extending RISE. It first describes what a paradigm is and the interface a paradigm must satisfy, and only then gives the concrete paradigms that ship with RISE as worked examples. Three paradigms are currently provided (each under paradigms/):

occbin – piecewise-linear, OccBin-style occasionally-binding constraints (occbin.paradigm);
recursive local linearization (RLL) (rll.paradigm);
\(k\)-th order Taylor projection (TP) (rtp.paradigm).

Each is described in Examples, after the general interface.

19.1. Goals

The paradigm architecture is designed to deliver:

a clean separation between model structure and computational method;
support for multiple parameterizations of the same model;
extensibility to advanced methods without touching RISE core;
prevention of cross-talk between parameterizations;
strictly read-only access to RISE internals by paradigm code.

A model contains exactly one paradigm at a time, but may carry multiple parameter vectors, each producing its own solution object.

19.2. What a paradigm is

A paradigm \(\mathcal{P}\) is a computational universe that

solves a parameterized model \(m(\theta)\),
defines the internal representation of the solution \(S(\theta)\), and
provides operators acting on that solution (e.g. a one-step forecast).

It is not the model, not a parameter vector, and not a container for multiple solutions. Formally,

\[S(\theta) = \mathcal{P}.\texttt{solve\_engine}\bigl(m(\theta),\, \texttt{regular\_solver}\bigr),\]

where \(S(\theta)\) is partly transparent to RISE (through the standard envelope below) and partly opaque (through the S.solution payload).

19.3. Responsibilities of a paradigm

A paradigm must:

operate on a single parameter vector at a time;
return a single solution object S conforming to the standard envelope;
provide operators derived from that solution;
remain unaware of parameter multiplicity;
treat the model object M as read-only – it may call regular_solver(M, ...) to obtain derivatives, but must never modify M nor return a modified M to RISE;
avoid modifying RISE core structures or global options.

A paradigm may use RISE’s baseline solver internally, use RISE’s constraint enforcement, and keep mutable internal state inside S.solution (see Mutable solution objects).

A paradigm must not store global state across parameterizations, return a modified M, or manage arrays of solutions.

19.4. The standard solution envelope

Every paradigm’s solve_engine returns a struct S with the fields below. This standard envelope lets RISE orchestrate the workflow (convergence checking, re-solve decisions, filter routing) without knowing anything about the paradigm’s internals.

Field	Type	Meaning
`is_solved`	`[1x2 logical]`	`[ss_ok, dynamics_ok]`: first entry, the steady state / baseline solve succeeded; second entry, the dynamics / projection solve succeeded.
`warrant_resolving`	`logical`	`true` if the solve failed or degraded and RISE should attempt a re-solve.
`retcode`	`integer`	`0` = success (RISE’s usual `retcode` convention).
`H`	matrix or `[]`	measurement-error covariance for the Kalman filter; `[]` if the paradigm provides none.
`solution`	handle object	opaque, paradigm-owned payload. RISE stores it but never inspects it. All paradigm-specific state lives here.

The envelope is deliberately minimal. Fields specific to perturbation theory (eigenvalues, number of solutions, Tz, …) are not required, because they are meaningless or misleading for paradigms such as RLL and TP; those paradigms keep such quantities inside S.solution.

A paradigm should start from a failed default and populate fields as each stage succeeds, so partial failures are reported correctly:

function S = make_envelope()
S.is_solved         = [false, false];
S.warrant_resolving = true;
S.retcode           = 1;
S.H                 = [];
S.solution          = [];
end

19.5. The minimum interface contract

A paradigm exposes an engines struct with these fields.

name (required)

A short string identifier, e.g. 'perturbation', 'rll', 'rtp'.

solve_engine (required)

S = solve_engine(M, regular_solver)

M is read-only: the paradigm may call regular_solver(M, ...) internally, but the solved model is used internally only and never returned. S must conform to the standard envelope, and no global or shared state may be modified.

one_step_engine (required)

one_step = one_step_engine(S)

The only input is S – the model M is not passed here, so everything the closure needs must have been captured in S.solution during solve_engine. It returns a closure

[x_t, retcode, extra...] = one_step(x_tm1, eps_t, regime)

with x_tm1 the lagged endogenous vector \(x_{t-1}\), eps_t the shock vector \(\varepsilon_t\), and regime the integer index of the current Markov regime. It returns a nonzero retcode if the step fails; extra... are optional outputs such as the local linear coefficients (k, H, G) consumed by the smoother.

filter_engine (optional)

Future extension for filtering / smoothing.

use_rise_constraint_enforcer (optional)

Logical, default true.

19.6. Mutable solution objects

For RLL and TP the effective policy function is state-dependent: it is updated at each simulation step as the economy evolves, whereas RISE stores the solution S once. The recommended resolution is a MATLAB handle class that can be copied (matlab.mixin.Copyable): one_step_engine(S) gives its closure a copy of S.solution, which the closure then updates in place, step after step. Each simulation, forecast or filter run builds its own closure, so each starts from the solved policy and none sees the state another run left behind; RISE’s stored S.solution stays as solved. Updating S.solution itself in place makes a result depend on what ran before it: two identical simulations of a re-anchored Taylor projection gave two paths.

classdef SolutionHandle < matlab.mixin.Copyable
    properties
        % paradigm-specific fields populated by solve_engine, e.g.
        T        % TP: current polynomial coefficients
        k, H, G  % RLL: current local linear policy
        % ...
    end
end

The workflow is:

solve_engine constructs a fresh SolutionHandle, populates it, and stores it as S.solution;
one_step_engine(S) takes payload = copy(S.solution) and builds a closure capturing that copy by reference;
at each step the closure updates payload in place; S.solution is left as solved;
RISE reads only the envelope fields and never the payload.

rtp.SolutionHandle and rll.SolutionHandle follow this pattern.

Note

Because RISE never reads S.solution directly, this pattern is safe. The only contract is that one_step(x_tm1, eps_t, regime) always returns a valid x_t and retcode, whatever the payload’s internal state.

19.7. Responsibilities of RISE

RISE manages model structure and parameter multiplicity; calls solve_engine(M, regular_solver) for each \(\theta_i\) (passing M read-only); stores each \(S_i\) in model.solution(i); and routes simulation, forecasting and filtering through the closures returned by one_step_engine (and filter_engine). The workflow is

\[\theta_i \;\rightarrow\; m(\theta_i) \;\rightarrow\; S_i = \texttt{solve\_engine}(m(\theta_i),\,\texttt{regular\_solver}) \;\rightarrow\; \text{store } S_i \;\rightarrow\; \text{simulation / filtering}.\]

RISE must not inspect or modify S.solution, assume any format inside it, or reuse a one_step closure across different parameter vectors.

19.8. Isolation, safety and persistence

To prevent cross-talk between parameterizations, each call to solve_engine must construct a fresh SolutionHandle (copying a handle copies only the reference, not the data), no paradigm may keep persistent state outside S.solution, and RISE must not reuse engines across parameter vectors unless the paradigm is stateless.

For saving models to disk, paradigms should provide serialization hooks

blob = serialize_solution(S)
S    = deserialize_solution(blob)

which RISE stores without inspecting.

19.9. Examples

19.9.1. The OccBin paradigm

The piecewise-linear (Guerrieri–Iacoviello, “OccBin”) method for an arbitrary number of occasionally-binding constraints is shipped as a ready-to-use paradigm, occbin.paradigm. The whole method lives in the paradigm package: RISE exchanges only the standard solution envelope and never needs a dedicated occbin option or any in-core occbin state.

Write the model as a Markov-switching model with one reference regime in which no constraint binds, plus, for each constraint, a chain whose second state is the binding branch (e.g. the policy-rate equation becomes R = bind*r_zlb + (1-bind)*R_TAYLOR). Then activate the paradigm at solve time:

m = dsge_model('mymodel');

C = { {'R>=1','ocb'} };   % constraint R >= 1 on the 'ocb' chain

m = set(m,'solve_paradigm',{@occbin.paradigm,'ref',1,'constraints',C});
m = solve(m);
sims = simulate(m,'simul_periods',1000);

RISE stays agnostic about what the paradigm needs: the paradigm itself imposes a single common (unique, imposed) steady state and pins the constraint chains’ transition probabilities to zero – the binding regime is chosen by the constraints, not by chance – so you do not pass sstate_imposed / sstate_unique and do not set those transition probabilities in the calibration.

constraints is a cell with one entry per constraint, each an intuitive {inequality, chain} pair:

element	meaning
`inequality`	`'VAR OP BOUND'` – `VAR` an endogenous variable, `OP` one of `>= > <= <`, `BOUND` a number or a parameter name (e.g. `'R>=r_zlb'`). `>=`/`>` binds when `VAR` falls below the bound; `<=`/`<` binds when it rises above.
`chain`	name of the Markov chain that switches this constraint (char)

Several constraints are just several pairs:

C = { {'maxlev<=0','debt'}, ...   % borrowing limit, 'debt' chain
      {'r>=1',     'zlb'} };       % ZLB on the policy rate, 'zlb' chain

The explicit struct form – struct('var','R','dir',-1,'bound',1,'chain','ocb'), with dir=-1 ⇔ >= and dir=+1 ⇔ <= – is still accepted for each entry, so existing scripts keep working.

ref is the reference regime index (where every chain is in its non-binding state). Further options: maxspell (default 200, the longest binding spell tried before giving up), use_pinv (default true) and anchor_at_bound (default true: the constant of a binding regime is taken where the binding constrained variables sit on their bounds, so they land on them even when the binding branch is nonlinear in them; false is the classical constant, the residual at the steady state; see Occasionally-binding constraints). The paradigm enforces the constraint itself, so use_rise_constraint_enforcer is false.

At solve time the paradigm imposes the (common) steady state, evaluates the structural system under each regime, and solves only the reference regime. At simulation time it takes a reference-regime step each period; if a constraint is violated it grows a binding spell (backward recursion with the reference rule as the terminal condition) until the period after the spell is slack again. For several constraints the active set at each date is the set of violated constraints, and the matching regime is read from the model’s regime table, so the method scales to any number of chains.

There is no filter_engine: OccBin has no associated filtering procedure, so estimation/filtering under this paradigm is not provided.

The model-side construction (the binding chain and the convex-combination equation) is described in Occasionally-binding constraints.

19.9.2. Recursive local linearization (RLL)

solve_engine calls regular_solver read-only, extracts the baseline local policy \((k_0, H_0, G_0)\), stores everything in an rll.SolutionHandle, and returns the standard envelope with S.solution = payload.
one_step_engine(S) returns a closure that works on its own copy of the payload: it builds a zero-shock stance, re-evaluates the structural derivatives at that stance, solves the local linear system for updated \((k_t, H_t, G_t)\), updates its copy in place, and returns \(x_t\) together with \((k_t, H_t, G_t)\) for the smoother. Each simulation therefore starts from the solved policy.
Options, passed as {@rll.paradigm, Name, Value, ...}:
- 'solver' solves each local system: 'rise_qz' (the default) solves it exactly by a generalized Schur decomposition, which the constant-parameter local systems of RLL allow; 'fi' (functional iteration), 'cr' (cyclic reduction) and 'newton' iterate, at most 'MaxIter' (100) times to the tolerance 'TolFun' (1e-8), and far from the steady state functional iteration may need many more iterations than that; 'shalva_qz' is another QZ variant.
- 'use_steady_state': with false (the default) the step is \(x_t = k_t + H_t x_{t-1} + G_t \varepsilon_t\), where the constant \(k_t\) comes from the model evaluated at the current stance; with true it is \(x_t = \bar{x} + H_t (x_{t-1} - \bar{x}) + G_t \varepsilon_t\), the local slopes applied around the global steady state \(\bar{x}\). Far from the steady state the latter is usually less accurate than first order: on the growth model from capital at 0.4 times its steady state, the largest error in consumption against perfect foresight is 3.3% of the steady state with the defaults, 5.9% at first order and 24% with 'use_steady_state', true.
- A re-linearization that fails (an iterative solver that does not converge, or a local system with no unique stable solution) or that the guardrail rejects ('guardrail', 'max_stable_root': a complex, non-finite or explosive local solution) keeps the previous local solution, with one warning per simulation. At an occasionally binding bound such as R = max(0, ...) the local system holds the rate at the bound, has no unique stable solution, and RLL keeps the previous local solution there; the constraint enforcer still holds the rate at the bound.
These were 'fi' and true before the defaults changed; code that passed 'solver', 'rise_qz', 'use_steady_state', false gets the same results as before.

19.9.3. \(k\)-th order Taylor projection (TP)

solve_engine calls regular_solver read-only, uses the perturbation solution as a warm start, runs rtp.solve order by order to obtain \(T^* = \{T_1^*,\dots,T_h^*\}\), and stores \(T^*\), the topology and get_F in an rtp.SolutionHandle. The solved coefficients are converted to divided differences (Taylor coefficients) before they are stored, so that

\[x_t = \sum_{q=0}^{K} T^*_{\text{regime}}\{q{+}1\}\, \bigl(z_t - z_{\text{anchor}}\bigr)^{\otimes q}.\]
one_step_engine(S) returns a closure with two schemes, selected by the reanchor option:
- re-anchored (reanchor=true, the default — Levintal’s scheme): at each step it rebuilds get_F with the anchor at the current state \(x_{t-1}\) (rtp.from_rise(m, x_{t-1}, opts), which evaluates the model derivatives along the trajectory, with \(x_{t+1}=g(\text{forward state})\)), re-solves rtp.solve there, and evaluates the polynomial. This re-centers the approximation on the current state and is what makes TP more accurate than perturbation away from the steady state, at the cost of a projection re-solve per step. Only the state entries of \(x_{t-1}\) define the anchor: the re-solve starts from the current polynomial, whose value at the new state is the guess for \(x_t\), so the other entries of the date-0 vector, which a simulation plan may leave at any value, play no role. Each closure re-anchors its own copy of the solution, so every simulation starts from \(T^*\).
- fixed polynomial (reanchor=false): evaluates the stored steady-state anchored \(T^*\) directly. Fast, but at the steady-state anchor this coincides with perturbation (single-anchor TP and perturbation impose the same conditions there).
In both cases \(z_t = (p_{t-1}, b_{t-1}, \sigma{=}1, \varepsilon_t)\) and \(z_{\text{anchor}} = (p_{t-1}, b_{t-1}, \sigma{=}0, \varepsilon{=}0)\).
The projection equations hold in expectation over next period’s shocks. Next period’s shock enters the next period’s state as \(\sigma\,\varepsilon_{t+1}\), and the residual of regime \(i\) is

\[R_i(z) = \mathbb{E}_{\varepsilon_{t+1}} \sum_j f_{ij}\bigl(v_{ij}(z;\varepsilon_{t+1})\bigr), \qquad \varepsilon_{t+1} \sim N(0, I).\]

At the anchor \(\sigma = 0\), so the \(q\)-th derivative of the residual there is a polynomial of degree at most \(q\) in \(\varepsilon_{t+1}\), and a rule exact for the moments of \(N(0,I)\) up to the order \(K\) gives the expectation exactly. rtp.shock_nodes builds it: the node 0 at order 1, the \(n_e+1\) vertices of a regular simplex at order 2, the \(2n_e\) nodes \(\pm\sqrt{n_e}\,e_i\) at order 3, Stroud’s degree-5 rule (\(2n_e^2+1\) nodes) at orders 4 and 5, and tensor-product Gauss–Hermite beyond. rtp.residual and rtp.jacobian_analytical sum over the nodes. This gives the policy its risk correction at order 2 and above (the \(\sigma^2\) terms): on the RBC model at a shock volatility of 0.05, a shock-free path from the steady state takes capital up by 0.127 under order-2 perturbation and by 0.126 under re-anchored order-2 Taylor projection. The option 'expectation','none' sets \(\varepsilon_{t+1} = 0\) instead: certainty equivalence, which is what the projection did before the option existed. With a news horizon \(k > 0\) the shock drawn next period is the one announced \(k\) periods ahead, the last shock slot of the next period’s state.
The projection residual of regime \(i\) composes the next period’s policy with the state transition, \(T_j(\Psi_i(z))\), so it needs \(T_j\) and its derivatives at \(\Psi_i(z_{\text{anchor}})\), where the state lands one period after the anchor, not at \(T_j\)’s own anchor (the two coincide only at a steady-state anchor of a constant-parameter model). rtp.residual moves \(T_j\) there (rtp.recenter, exact for a polynomial; the displacement comes from rtp.forward_shift), the lead point is evaluated for each pair of regimes, and rtp.jacobian_analytical has the matching terms. Since the model’s derivatives are evaluated at that landing point, which every coefficient of \(T_j\) moves, the analytical Jacobian also differentiates them through it.
The projection equations are solved anew at every state, and they need not have a root there. Far from the steady state, the root near the policy can disappear at a fold (a point where the Jacobian is singular): on the growth model written in levels, the order-2 equations lose it below 0.52 of the steady-state capital. Newton then either fails or converges to another root (capital 80 times its value, negative consumption). A converged re-solve is therefore kept only if today’s policy at the new state, the constant of each regime’s polynomial, is within reanchor_max_jump (default 0.5, relative to \(\max(|x|, 1)\), variable by variable) of the last converged policy evaluated at that state; good roots stay below 0.09 on that model, at orders 1 to 3 and from 0.9 to 0.3 of the steady-state capital. A failed or rejected re-solve warns (rtp:paradigm:stepFallback), and the step reuses the last converged policy (in levels at order 2 from 0.4 of the steady-state capital this is no more accurate than order-2 perturbation, and slow).

Write the model in logs. The same growth model with consumption and capital in logs (@endogenous(log) C K) keeps the root of its projection equations at every distance, and Taylor projection beats perturbation from every starting point. Largest error in consumption against perfect foresight, % of the steady state:

capital at the start / steady state	0.9	0.7	0.5	0.4	0.3
perturbation, order 2, logs	0.0008	0.026	0.160	0.329	0.648
Taylor projection, order 2, logs	0.00002	0.0003	0.0012	0.0017	0.0015
Taylor projection, order 3, logs	0.00001	0.00003	0.0002	0.0007	0.0018
perturbation, order 2, levels	0.006	0.188	1.05	2.05	3.74
Taylor projection, order 2, levels	0.003	0.137	1.74 (2 steps fall back)	2.05 (10)	3.74 (16)
Taylor projection, order 3, levels	0.0002	0.030	0.389	0.921	1.75

In levels, order 3 keeps its root down to 0.3 of the steady-state capital and is about as accurate as order-3 perturbation; in logs, order 2 is 40 to 450 times more accurate than order-2 perturbation.

Solution accuracy can be quantified with the model’s accuracy method; on the RBC model re-anchored TP attains an Euler error 8–13 times smaller than perturbation at \(\pm 30\%\) capital.

19.9.4. Risk-adjusted linear (RAL)

The risk-adjusted linear paradigm, ral.paradigm, turns a higher-order perturbation solution into a linear one that still carries the effect of uncertainty. solve_engine calls regular_solver read-only at an arbitrary internal order \(k\) (the user’s choice), then takes the linearization of that policy function: it keeps every term that fits a linear representation in the states and shocks – the constant (the precautionary \(\sigma^a\) level shifts) and the linear and risk-adjusted slope terms – and discards the genuine curvature (\(\text{state}\times\text{state}\), \(\text{state}\times\text{shock}\), …) and every higher-order block. The result is the affine policy

\[x_t = \bar{x} + c + A\,(x_{t-1}-\bar{x}) + B\,\varepsilon_t,\]

where \(c\) is the precautionary level correction and \(A, B\) are the risk-adjusted slopes. Under regime switching the linearization is taken per regime, so each regime gets its own precautionary level and slopes.

Unlike the other paradigms, RAL dissolves after solving. Because its output is a genuine first-order solution, solve_engine hands the native state_space back to RISE (signalling S.dissolve = true in the envelope) and steps aside: the dispatcher adopts the native solution and clears the paradigm, so from then on the model is an ordinary first-order model. Every downstream path – simulate, irf, and in particular the analytic linear Kalman filter that RISE auto-selects when no filter is named – runs natively, with no bespoke RAL one_step or filter. The precautionary level is carried in the sigma (perturbation) column of the first-order solution, so it enters simulation and the filter constant and correctly cancels out of impulse responses (which are deviations).

Activate it at solve time:

m = dsge_model('mymodel','max_deriv_order',2);
m = set(m,'solve_paradigm',{@ral.paradigm,'order',2});
m = solve(m);                       % m is now a native first-order model
sims = simulate(m,'simul_periods',1000);

The order option (default 2) is the internal perturbation order from which the risk correction is read: order 2 captures the leading precautionary level and the risk-adjusted slopes; higher orders additionally fold in the \(\sigma^a\) (\(a>2\)) corrections. The output is first order regardless. Evaluation is first-order cheap and needs no pruning – the only nonlinear content, the curvature, has been removed – and the residual gap to the full order-\(k\) solution is exactly that discarded curvature.

19.10. Design principles

A model has exactly one paradigm.
A model may have multiple parameterizations; each has one solution S.
M is always read-only inside a paradigm.
S.solution is always opaque to RISE.
The envelope gives RISE just enough to orchestrate workflows.
Mutable simulation-time state lives inside S.solution via handle semantics.
Paradigm-specific solution concepts (eigenvalues, number of solutions, Taylor coefficients) live inside S.solution, not in the envelope.
A dissolving paradigm (e.g. ral.paradigm) is the one exception to the first and fourth points: when its result is itself a standard solution it sets S.dissolve = true and hands back a native state_space; the dispatcher adopts it and clears the paradigm, leaving a paradigm-free model that RISE filters and simulates with its ordinary machinery.