Numerical Algorithms

To calculate the risk-adjusted linearization, we need to solve a system of nonlinear equations. These equations are generally solvable using Newton-type methods. The package currently has two available algorithms, relaxation and homotopy continuation

`solve!`

The primary interface for calculating a risk-adjusted linearization once a RiskAdjustedLinearization object is created is the function solve!. The user selects the desired numerical algorithm through algorithm keyword of solve!.

All of the available algorithms need to solve a system of nonlinear equations. We use nlsolve for this purpose, and all keyword arguments for nlsolve can be passed as keyword arguments to solve!, e.g. autodiff and ftol. The user can also exploit sparsity in the Jacobian of the system of nonlinear equations to accelerate nlsolve by using the keywords sparse_jacobian, sparsity, colorvec, jac_cache, and/or sparsity_detection. For details, see Exploiting Sparsity.

RiskAdjustedLinearizations.solve! — Function

solve!(m; algorithm = :relaxation, autodiff = :central, verbose = :high, kwargs...)
solve!(m, z0, y0; kwargs...)
solve!(m, z0, y0, Ψ0; kwargs...)

computes the risk-adjusted linearization of the dynamic economic model described by m and updates m with the solution, e.g. the coefficients $(z, y, \Psi)$.

The three available solve! algorithms are slight variations on each other.

Method 1: uses the z, y, and Ψ fields of m as initial guesses for $(z, y, \Psi)$ and proceeds with the numerical algorithm specified by algorithm
Method 2: uses z0 and y0 as initial guesses for the deterministic steady state, which is then used as the initial guess for $(z, Y, \Psi)$ for the numerical algorithm specified by algorithm.
Method 3: uses z0, y0, and Ψ0 as initial guesses for $(z, Y, \Psi)$ and proceeds with the numerical algorithm specified by algorithm.

Inputs

m::RiskAdjustedLinearization: object holding functions needed to calculate the risk-adjusted linearization
z0::AbstractVector{S1}: initial guess for $z$
y0::AbstractVector{S1}: initial guess for $y$
Ψ0::AbstractVector{S1}: initial guess for $\Psi$
S1 <: Real

Keywords

algorithm::Symbol = :relaxation: which numerical algorithm to use? Can be one of [:relaxation, :homotopy, :deterministic]
autodiff::Symbol = :central: use autodiff or not? This keyword is the same as in nlsolve
use_anderson::Bool = false: use Anderson acceleration if the relaxation algorithm is applied. Defaults to false
step::Float64 = .1: size of step from 0 to 1 if the homotopy algorithm is applied. Defaults to 0.1
sparse_jacobian::Bool = false: if true, exploit sparsity in the Jacobian in calls to nlsolve using SparseDiffTools.jl. If jac_cache and sparsity are nothing, then solve! will attempt to determine the sparsity pattern.
sparsity::Union{AbstractArray, Nothing} = nothing: sparsity pattern for the Jacobian in calls to nlsolve
colorvec = nothing: matrix coloring vector for sparse Jacobian in calls to nlsolve
jac_cache = nothing: pre-allocated Jacobian cache for calls to nlsolve during the numerical algorithms
sparsity_detection::Bool = false: If true, use SparsityDetection.jl to detect sparsity pattern (only relevant if both jac_cache and sparsity are nothing). If false, then the sparsity pattern is determined by using finite differences to calculate a Jacobian and assuming any zeros will always be zero. Currently, SparsityDetection.jl fails to work.

The solution algorithms all use nlsolve to calculate the solution to systems of nonlinear equations. The user can pass in any of the keyword arguments for nlsolve to adjust the settings of the nonlinear solver.

For the keywords relevant to specific methods, see the docstring for the underlying method being called. Note these methods are not exported.

:relaxation -> relaxation!
:homotopy -> homotopy!
:deterministic -> deterministic_steadystate!

source

Relaxation

The first and default numerical algorithm is a relaxation algorithm. The key problem in solving the equations characterizing $(z, y, \Psi)$ is that it is difficult to jointly solve the nonlinear matrix equation for $\Psi$ along with the steady-state equations for $z$ and $y$ due to the presence of the entropy term. The relaxation algorithm splits the solution of these equations into two steps, which allows us to calculate guesses of $\Psi$ using linear algebra. It is in this sense that this iterative algorithm is a relaxation algorithm.

The system of equations characterizing the coefficients $(z, y, \Psi)$ are solved iteratively in two separate steps. Given previous guesses $(z_{n - 1}, y_{n - 1}, \Psi_{n - 1})$, we calculate $(z_n, y_n)$ such that

\[\begin{aligned} 0 & = \mu(z_n, y_n) - z_n,\\ 0 & = \xi(z_n, y_n) + \Gamma_5 z_n + \Gamma_6 y_n + \mathcal{V}(z_{n - 1}),\\ \end{aligned}\]

is satisfied. In other words, we hold the entropy term constant and update $(z_n, y_n)$ in the remaining terms. The coefficients are solved efficiently through nlsolve with $(z_{n - 1}, y_{n - 1})$ as initial guesses.

Then we compute $\Psi_n$ by solving

\[\begin{aligned} 0 & = \Gamma_3 + \Gamma_4 \Psi_n + (\Gamma_5 + \Gamma_6 \Psi_n)(\Gamma_1 + \Gamma_2 \Psi_n) + J\mathcal{V}(z_{n - 1}). \end{aligned}\]

with a Generalized Schur decomposition (also known as QZ decomposition). Notice that we also hold the Jacobian of the entropy constant. Only after we have a new round of $(z_n, y_n, \Psi_n)$ do we update the entropy-related terms.

Convergence is achieved once $(z_n, y_n, \Psi_n)$ are sufficiently close under some norm. By default, we use the $L^\infty$ norm (maximum absolute error).

Homotopy Continuation

When the deterministic steady state exists, it is typically an easy problem to solve numerically. We can therefore use the equations characterizing the deterministic steady state for a homotopy continuation method. Let $q$ be the embedding parameter. Then the homotopy continuation method iteratively solves

\[\begin{aligned} 0 & = \mu(z, y) - z,\\ 0 & = \xi(z, y) + \Gamma_5 z + \Gamma_6 y + q \mathcal{V}(z),\\ 0 & = \Gamma_3 + \Gamma_4 \Psi + (\Gamma_5 + \Gamma_6 \Psi)(\Gamma_1 + \Gamma_2 \Psi) + q J\mathcal{V}(z) \end{aligned}\]

for the coefficients $(z_q, y_q, \Psi_q)$ by increasing $q$ from 0 to 1.

Blanchard-Kahn Conditions

At the end of solve!, we check the stochastic steady state found is locally unique and saddle-path stable by checking what are known as the Blanchard-Kahn conditions. Standard references for computational macroeconomics explain what these conditions are, so we defer to them (e.g. Blanchard-Kahn (1980), Klein (2000), and Sims (2002)). For the stochastic steady state, these conditions are essentially identical to the conditions for the deterministic steady state, but the Jacobian of the expectational equations to $z_t$ also includes the Jacobian of the entropy. In the deterministic steady state, the entropy is zero, hence the Jacobian of the entropy is zero. In the stochastic steady state, the entropy is no longer zero and varies with $z_t$, hence the Jacobian of the expectational equations to $z_t$ depends on entropy.

Exploiting Sparsity

When solving for the deterministic or stochastic steady state, this package solves a system of nonlinear equations by calling nlsolve, whose underlying algorithms typically require calculating the Jacobian of the system of nonlinear equations. For many economic models, this system is sparse because each individual equation usually depends on a small subset of the coefficients $(z, y, \Psi)$. To exploit this sparsity and further accelerate computation time, we can use methods implemented by SparseDiffTools.jl. For an example, please see this script.

We automate the setup process by letting the user pass the keyword sparse_jacobian = true to solve!. If this keyword is true, then there are three ways to exploit sparsity.

If no other keywords are passed, then solve! will attempt to determine the sparsity pattern. By default, the sparsity pattern is determined by using finite differences to calculate a Jacobian and assuming any zeros will always be zero. If the keyword sparsity_detection = true, then solve! will try using SparsityDetection.jl. Currently, the latter approach does not work with RiskAdjustedLinearizations.jl.
The keyword sparsity can be used to specify the sparsity pattern of the Jacobian. If colorvec is not also passed, then the matrix coloring vector is computed based on sparsity.
The keyword jac_cache allows the user to specify the sparsity pattern of the Jacobian and additionally pre-allocate the Jacobian's cache, which potentially achieves speed gains by avoiding extra allocations when the Jacobian function is repeatedly constructed.

If solve! is called once, then the first two approaches are essentially the same. If solve! is repeatedly called (e.g. if the model's parameters are changed), then the second two approaches are strictly faster because computing the sparsity pattern or pre-allocating the Jacobian's cache only needs to be done once, as long as the system of equations does not change.

To simplify using the sparsity, colorvec, and jac_cache keywords, we implement two helper functions, compute_sparsity_pattern and preallocate_jac_cache. The first function calculates sparsity and colorvec while the second ones computes jac_cache. See the docstrings below and this example for more details.

Some additional caveats on these methods:

Creating a cached Jacobian with automatic differentiation via ForwardColorJacCache will not work because the objective function changes in each loop of the homotopy and relaxation algorithms, so the cached Dual matrices will have information on the wrong function after a loop completes. Currently, RiskAdjustedLinearizations.jl has not implemented a way to update the information on the function required by the Dual matrices.
If automatic differentiation does not work with dense Jacobians due to problems with reinterpreting the chunk size, then it will also not work when using sparse Jacobians.

Docstrings

RiskAdjustedLinearizations.relaxation! — Function

relaxation!(ral, xₙ₋₁, Ψₙ₋₁; tol = 1e-10, max_iters = 1000, damping = .5, pnorm = Inf,
            schur_fnct = schur!, autodiff = :central, use_anderson = false, m = 5,
            verbose = :none, kwargs...)

solves for the coefficients $(z, y, \Psi)$ of a risk-adjusted linearization by the following relaxation algorithm:

Initialize guesses for $(z, y, \Psi)$
Do until convergence
a) Solve for $(z, y)$ using the expectational and state transition equations and fixing $\Psi$.
b) Use a QZ decomposition to solve for $\Psi$ while fixing $(z, y)$.

Types:

S1 <: Number
S2 <: Real
S3 <: Real

Inputs

m::RiskAdjustedLinearization: object holding functions needed to calculate the risk-adjusted linearization
xₙ₋₁::AbstractVector{S1}: initial guess for $(z, y)$
Ψₙ₋₁::AbstractVector{S1}: initial guess for $\Psi$

Keywords

tol::S2: convergence tolerance of residual norm for relaxation algorithm
max_iters::Int: maximumm number of iterations
damping::S2: guesses are updated as the weighted average xₙ = damping * proposal + (1 - damping) * xₙ₋₁.
pnorm::S3: norm for residual tolerance
schur_fnct::Function: function for calculating the Schur factorization during QZ decomposition
autodiff::Symbol: specifies whether to use autoamtic differentiation in nlsolve (and is the same keyword as the autodiff keyword for nlsolve)
use_anderson::Bool: set to true to apply Anderson acceleration to the fixed point iteration of the relaxation algorithm
m::Int: m coefficient if using Anderson acceleration
sparse_jacobian::Bool = false: if true, exploit sparsity in the Jacobian in calls to nlsolve using SparseDiffTools.jl. If jac_cache and sparsity are nothing, then relaxation! will attempt to determine the sparsity pattern.
sparsity::Union{AbstractArray, Nothing} = nothing: sparsity pattern for the Jacobian in calls to nlsolve
colorvec = nothing: matrix coloring vector for sparse Jacobian in calls to nlsolve
jac_cache = nothing: pre-allocated Jacobian cache for calls to nlsolve during the numerical algorithms
sparsity_detection::Bool = false: If true, use SparsityDetection.jl to detect sparsity pattern (only relevant if both jac_cache and sparsity are nothing). If false, then the sparsity pattern is determined by using finite differences to calculate a Jacobian and assuming any zeros will always be zero. Currently, SparsityDetection.jl fails to work.
verbose::Symbol: verbosity of information printed out during solution. a) :low -> statement when homotopy continuation succeeds b) :high -> statement when homotopy continuation succeeds and for each successful iteration

source

RiskAdjustedLinearizations.homotopy! — Function

homotopy!(m, xₙ₋₁; step = .1, pnorm = Inf, verbose = :none, kwargs...)

solves the system of equations characterizing a risk-adjusted linearization by a homotopy method with embedding parameter $q$, which steps from 0 to 1, with $q = 1$ obtaining the true solution.

Currently, the only algorithm for choosing $q$ is a simple uniform step search. Given a step size $\Delta$, we solve the homotopy starting from $q = \Delta$ and increase $q$ by $\Delta$ until $q$ reaches 1 or passes 1 (in which case, we force $q = 1$).

Types:

S1 <: Number

Inputs

m::RiskAdjustedLinearization: object holding functions needed to calculate the risk-adjusted linearization
xₙ₋₁::AbstractVector{S1}: initial guess for $(z, y, \Psi)$

Keywords

step::Float64: size of the uniform step from step to 1.
pnorm::Float64: norm under which to evaluate the errors after homotopy succeeds.
sparse_jacobian::Bool = false: if true, exploit sparsity in the Jacobian in calls to nlsolve using SparseDiffTools.jl. If jac_cache and sparsity are nothing, then homotopy! will attempt to determine the sparsity pattern.
sparsity::Union{AbstractArray, Nothing} = nothing: sparsity pattern for the Jacobian in calls to nlsolve
colorvec = nothing: matrix coloring vector for sparse Jacobian in calls to nlsolve
jac_cache = nothing: pre-allocated Jacobian cache for calls to nlsolve during the numerical algorithms
sparsity_detection::Bool = false: If true, use SparsityDetection.jl to detect sparsity pattern (only relevant if both jac_cache and sparsity are nothing). If false, then the sparsity pattern is determined by using finite differences to calculate a Jacobian and assuming any zeros will always be zero. Currently, SparsityDetection.jl fails to work.
verbose::Symbol: verbosity of information printed out during solution. a) :low -> statement when homotopy continuation succeeds b) :high -> statement when homotopy continuation succeeds and for each successful iteration

source

RiskAdjustedLinearizations.blanchard_kahn — Function

blanchard_kahn(m::RiskAdjustedLinearization; deterministic::Bool = false, verbose::Symbol = :high)

checks the Blanchard-Kahn conditions for whether a first-order perturbation is saddle-path stable or not.

If verbose is :low or :high, a print statement will be shown if the Blanchard-Kahn conditions are satisfied.

source

RiskAdjustedLinearizations.compute_sparsity_pattern — Function

compute_sparsity_pattern(f::Function, x::AbstractVector{<: Number}, nrow::Int;
                         sparsity::Union{AbstractArray, Nothing} = nothing,
                         sparsity_detection::Bool = false)

calculates the sparsity pattern of the Jacobian of the functions μ, ξ, and 𝒱.

Inputs

f: is the function to be differentiated, e.g. z -> 𝒱(z, Ψ, (1, 2))
x: the vector at which differentiation occurs
nrow: specifies the number of rows of the Jacobian

Keywords

sparsity: sparsity pattern of the Jacobian
sparsity_detection: if true, use SparsityDetection.jl to determine the sparsity pattern. If false, then the sparsity pattern is determined by using automatic differentiation to calculate a Jacobian and assuming any zeros are supposed to be zero.

source

compute_sparsity_pattern(m::RiskAdjustedLinearization, algorithm::Symbol; q::Float64 = .1,
                         sparsity::Union{AbstractArray, Nothing} = nothing,
                         sparsity_detection::Bool = false)

calculates the sparsity pattern and matrix coloring vector of the Jacobian of the nonlinear system of equations for either the deterministic or stochastic steady state, depending on which algorithm is called.

Keywords

q: step size for homotopy. Should satisfy 0 < q < 1 and is only required to ensure that the sparsity pattern is correctly determined when algorithm = :homotopy and thus the dependence of the entropy 𝒱 on the coefficients (z, y, Ψ) matters.
sparsity: sparsity pattern of the Jacobian of the nonlinear system of equations
sparsity_detection: if true, use SparsityDetection.jl to determine the sparsity pattern. If false, then the sparsity pattern is determined by using finite differences to calculate a Jacobian and assuming any zeros are supposed to be zero.

source

RiskAdjustedLinearizations.preallocate_jac_cache — Function

preallocate_jac_cache(m::RiskAdjustedLinearization, algorithm::Symbol; q::Float64 = .1,
                      sparsity::Union{AbstractArray, Nothing} = nothing,
                      sparsity_detection::Bool = false)

pre-allocates the cache for the Jacobian of the nonlinear system of equations for either the deterministic or stochastic steady state, depending on which algorithm is called.

Keywords

q: step size for homotopy. Should satisfy 0 < q < 1 and is only required to ensure that the sparsity pattern is correctly determined when algorithm = :homotopy and thus the dependence of the entropy 𝒱 on the coefficients (z, y, Ψ) matters.
sparsity: the sparsity pattern of the Jacobian of the nonlinear system of equations
sparsity_detection: if true, use SparsityDetection.jl to determine the sparsity pattern. If false, then the sparsity pattern is determined by using finite differences to calculate a Jacobian and assuming any zeros are supposed to be zero.

source