Risk-Adjusted Linearizations

Theory

Nonlinear Model

Most dynamic economic models can be formulated as the system of nonlinear equations

\[\begin{aligned} z_{t + 1} & = \mu(z_t, y_t) + \Lambda(z_t)(y_{t + 1} - \mathbb{E}_t y_{t + 1}) + \Sigma(z_t) \varepsilon_{t + 1},\\ 0 & = \log\mathbb{E}_t[\exp(\xi(z_t, y_t) + \Gamma_5 z_{t + 1} + \Gamma_6 y_{t + 1})]. \end{aligned}\]

The vectors $z_t\in \mathbb{R}^{n_z}$ and $y_t \in \mathbb{R}^{n_y}$ are the state and jump variables, respectively. The first vector equation comprise the transition equations of the state variables. The second vector equation comprise the model's expectational equations, which are typically the first-order conditions for the jump variables from agents' optimization problem. The exogenous shocks $\varepsilon_t \in\mathbb{R}^{n_\varepsilon}$ form a martingale difference sequence and can therefore be non-Gaussian. Given some differentiable mapping $\alpha:\mathbb{R}^{n_z}\rightarrow\mathbb{R}^{n_\varepsilon}$, the random variable $X_t = \alpha(z_t)^T \varepsilon_{t + 1}$ has the differentiable, conditional (on $z_t$) cumulant generating function

\[\begin{aligned} ccgf[\alpha(z_t) \mid z_t] = \log\mathbb{E}_t[\exp(\alpha(z_t)^T \varepsilon_{t + 1})]. \end{aligned}\]

The functions

\[\begin{aligned} \xi:\mathbb{R}^{2n_y + 2n_z}\rightarrow \mathbb{R}^{n_y},& \quad \mu:\mathbb{R}^{n_y + n_z}\rightarrow \mathbb{R}^{n_z},\\ \Lambda:\mathbb{R}^{n_z} \rightarrow \mathbb{R}^{n_z \times n_y}, & \quad \Sigma:\mathbb{R}^{n_z}\ \rightarrow \mathbb{R}^{n_z\times n_\varepsilon}, \end{aligned}\]

are differentiable. The first two functions characterize the effects of time $t$ variables on the expectational and state transition equations. The function $\Lambda$ characterizes heteroskedastic endogenous risk that depends on innovations in jump variables while the function $\Sigma$ characterizes exogenous risk. Note that $\Sigma$ is not the variance-covariance matrix of $\varepsilon_t$. The functions $\Lambda$ and $\Sigma$ can also depend on jump variables. Denote the jump-dependent versions as $\tilde{\Lambda}:\mathbb{R}^{n_z\times n_y} \rightarrow \mathbb{R}^{n_z \times n_y}$ and $\tilde{\Sigma}:\mathbb{R}^{n_z \times n_y}\ \rightarrow \mathbb{R}^{n_z\times n_\varepsilon}$. If there exists a mapping $y_t = y(z_t)$, then we define $\Lambda(z_t) = \tilde{\Lambda}(z_t, y(z_t))$ and $\Sigma(z_t) = \tilde{\Sigma}(z_t, y(z_t))$.

The expectational equations can be simplified as

\[\begin{aligned} 0 & = \log\mathbb{E}_t[\exp(\xi(z_t, y_t) + \Gamma_5 z_{t + 1} + \Gamma_6 y_{t + 1})]\\ & = \log[\exp(\xi(z_t, y_t))\mathbb{E}_t[\exp(\Gamma_5 z_{t + 1} + \Gamma_6 y_{t + 1})]]\\ & = \xi(z_t, y_t) + \Gamma_5\mathbb{E}_t z_{t + 1} + \Gamma_6 \mathbb{E}_t y_{t + 1} + \log\mathbb{E}_t[\exp(\Gamma_5 z_{t + 1} + \Gamma_6 y_{t + 1})] - (\Gamma_5\mathbb{E}_t z_{t + 1} + \Gamma_6 \mathbb{E}_t y_{t + 1})\\ & = \xi(z_t, y_t) + \Gamma_5\mathbb{E}_t z_{t + 1} + \Gamma_6 \mathbb{E}_t y_{t + 1} + \mathcal{V}(\Gamma_5 z_{t + 1} + \Gamma_6 y_{t + 1}), \end{aligned}\]

where the last term is defined to be

\[\begin{aligned} \mathcal{V}(x_{t + 1}) = \log\mathbb{E}_t[\exp(x_{t + 1})] - \mathbb{E}_t x_{t + 1}. \end{aligned}\]

As Lopez et al. (2018) describe it, this quantity "is a relative entropy measure, i.e. a nonnegative measure of dispersion that generalizes variance."

Risk-Adjusted Linearizations by Affine Approximation

Many economic models are typically solved by perturbation around the deterministic steady state. To break certainty equivalence so that asset pricing is meaningful, these perturbations need to be at least third order. However, even third-order perturbations can poorly approximate the true global solution. A key problem is that the economy may not spend much time near the deterministic steady state, so a perturbation around this point will be inaccurate.

Instead of perturbing the model's nonlinear equations around the deterministic steady state, we could perturb around the stochastic or "risky" steady state. This point is better for a perturbation because the economy will spend a large amount of time near the stochastic steady state. Lopez et al. (2018) show that an affine approximation of the model's nonlinear equation is equivalent to a linearization around the stochastic steady state. Further, they confirm that in practice this "risk-adjusted" linearization approximates global solutions of canonical economic models very well and outperforms perturbations around the deterministic steady state.

The affine approximation of an dynamic economic model is

\[\begin{aligned} \mathbb{E}[z_{t + 1}] & = \mu(z, y) + \Gamma_1(z_t - z) + \Gamma_2(y_t - y),\\ 0 & = \xi(z, y) + \Gamma_3(z_t - z) + \Gamma_4(y_t - y) + \Gamma_5 \mathbb{E}_t z_{t + 1} + \Gamma_6 \mathbb{E}_t y_{t + 1} + \mathcal{V}(z) + J\mathcal{V}(z)(z_t - z), \end{aligned}\]

where $\Gamma_1, \Gamma_2$ are the Jacobians of $\mu$ with respect to $z_t$ and $y_t$, respectively; $\Gamma_3, \Gamma_4$ are the Jacobians of $\xi$ with respect to $z_t$ and $y_t$, respectively; $\Gamma_5, \Gamma_6$ are constant matrices; $\mathcal{V}(z)$ is the model's entropy; and $J\mathcal{V}(z)$ is the Jacobian of the entropy;

and the state variables $z_t$ and jump variables $y_t$ follow

\[\begin{aligned} z_{t + 1} & = z + \Gamma_1(z_t - z) + \Gamma_2(y_t - y) + (I_{n_z} - \Lambda(z_t) \Psi)^{-1}\Sigma(z_t)\varepsilon_{t + 1},\\ y_t & = y + \Psi(z_t - z). \end{aligned}\]

The unknowns $(z, y, \Psi)$ solve the system of equations

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

Under an affine approximation, the entropy term is a nonnegative function $\mathcal{V}:\mathbb{R}^{n_z} \rightarrow \mathbb{R}_+^{n_y}$ defined such that

\[\begin{aligned} \mathcal{V}(z_t) \equiv \mathcal{V}_t(\exp((\Gamma_5 + \Gamma_6 \Psi)z_{t + 1})) = \vec{ccgf}[(\Gamma_5 + \Gamma_6 \Psi)(I_{n_z} - \Lambda(z_t) \Psi)^{-1} \Sigma(z_t) \mid z_t] \end{aligned}\]

where the notation $\vec{ccgf}$ means that each component $ccgf_i[\cdot \mid \cdot]$ is a conditional cumulant-generating function. Explicitly, define

\[\begin{aligned} A(z_t) = (\Gamma_5 + \Gamma_6 \Psi)(I_{n_z} - \Lambda(z_t) \Psi)^{-1} \Sigma(z_t) = [A_1(z_t), \dots, A_{n_y}(z_t)]^T. \end{aligned}\]

Each $A_i(z_t)$ is a mapping from $z_t$ to the $i$th row vector in $A(z_t)$. Then $ccgf_i[\cdot \mid \cdot]$ is

\[\begin{aligned} ccgf_i[A_i(z_t)\mid z_t] = \log\mathbb{E}_t[\exp(A_i(z_t) \varepsilon_{t + 1})]. \end{aligned}\]

Every $ccgf_i[\cdot \mid \cdot]$ corresponds to an expectational equation and thus acts as a risk correction to each one. In the common case where the individual components of $\varepsilon_{t + 1}$ are independent, $ccgf_i$ simplifies to

\[\begin{aligned} ccgf_i[A_i(z_t)\mid z_t] = \sum_{j = 1}^{n_\varepsilon}\log\mathbb{E}_t[\exp(A_{ij}(z_t) \varepsilon_{j, t + 1})], \end{aligned}\]

i.e. it is the sum of the cumulant-generating functions for each shock $\varepsilon_{j, t + 1}$.

To see why $\mathcal{V}_t(\exp((\Gamma_5 + \Gamma_6 \Psi)z_{t + 1}))$ can be expressed as the conditional cumulant-generating function of $(\Gamma_5 + \Gamma_6 \Psi)(I_{n_z} - \Lambda(z_t) \Psi)^{-1} \Sigma(z_t)$, observe that

\[\begin{aligned} \mathcal{V}_t(\exp((\Gamma_5 + \Gamma_6 \Psi)z_{t + 1})) & = \log\mathbb{E}_t[\exp((\Gamma_5 + \Gamma_6 \Psi)z_{t + 1})] - \mathbb{E}_t[(\Gamma_5 + \Gamma_6 \Psi)z_{t + 1}]\\ & = \log\left(\frac{\mathbb{E}_t[\exp((\Gamma_5 + \Gamma_6 \Psi)z_{t + 1})]}{\exp(\mathbb{E}_t[(\Gamma_5 + \Gamma_6 \Psi)z_{t + 1}])}\right). \end{aligned}\]

Since $\mathbb{E}_t[(\Gamma_5 + \Gamma_6 \Psi)z_{t + 1}]$ is a conditional expectation, it is measurable with respect to the time-$t$ information set. Therefore, we can move the denominator of the fraction within the logarithm inside the numerator's conditional expectation.

\[\begin{aligned} \mathcal{V}_t(\exp((\Gamma_5 + \Gamma_6 \Psi)z_{t + 1})) & = \log\mathbb{E}_t\left[\exp\left((\Gamma_5 + \Gamma_6 \Psi)z_{t + 1} - \mathbb{E}_t[(\Gamma_5 + \Gamma_6 \Psi)z_{t + 1}]\right)\right]\\ & = \log\mathbb{E}_t\left[\exp\left((\Gamma_5 + \Gamma_6 \Psi)(z_{t + 1} - \mathbb{E}_t[z_{t + 1}])\right)\right]. \end{aligned}\]

Using the postulated law of motion for states,

\[\begin{aligned} z_{t + 1} - \mathbb{E}_t[z_{t + 1}] & = \mu(z_t, y_t) + \Lambda(z_t)(y_{t + 1} - \mathbb{E}_t[y_{t + 1}]) + \Sigma(z_t) \varepsilon_{t + 1} - \mu(z_t, y_t)\\ & = \Lambda(z_t) \Psi (z_{t + 1} - \mathbb{E}_t[z_{t + 1}]) + \Sigma(z_t) \varepsilon_{t + 1}\\ (I - \Lambda(z_t) \Psi) (z_{t + 1} - \mathbb{E}_t[z_{t + 1}]) & = \Sigma(z_t) \varepsilon_{t + 1}. \end{aligned}\]

Therefore, the entropy term $\mathcal{V}_t(\cdot)$ becomes

\[\begin{aligned} \mathcal{V}_t(\exp((\Gamma_5 + \Gamma_6 \Psi)z_{t + 1})) & = \log\mathbb{E}_t\left[\exp\left((\Gamma_5 + \Gamma_6 \Psi)(I - \Lambda(z_t) \Psi)^{-1} \Sigma(z_t) \varepsilon_{t + 1}\right)\right]. \end{aligned}\]

The RHS is a vector of logarithms of expected values of linear combinations of the shocks $\varepsilon_{t, + 1}$ with coefficients given by the rows of $(\Gamma_5 + \Gamma_6 \Psi)(I - \Lambda(z_t) \Psi)^{-1} \Sigma(z_t)$. Thus, each element of $\mathcal{V}_t(\cdot)$ is the conditional cumulant-generating function of the random vector $\varepsilon_{t, + 1}$ evaluated at one of the rows of $(\Gamma_5 + \Gamma_6 \Psi)(I - \Lambda(z_t) \Psi)^{-1} \Sigma(z_t)$, as claimed.

Refer to Lopez et al. (2018) "Risk-Adjusted Linearizations of Dynamic Equilibrium Models" for more details about the theory justifying this approximation approach. See Deriving the conditional cumulant generating function for some guidance on calculating the ccgf, which many users may not have seen before.

Implementation as `RiskAdjustedLinearization`

We implement risk-adjusted linearizations of nonlinear dynamic economic models through the wrapper type RiskAdjustedLinearization. The user only needs to define the functions and matrices characterizing the equilibrium of the nonlinear model. Once these functions are defined, the user can create a RiskAdjustedLinearization object, which will automatically create the Jacobian functions needed to compute the affine approximation.

To ensure efficiency in speed and memory, this package takes advantage of a number of features that are easily accessible through Julia.

The Jacobians are calculated using forward-mode automatic differentiation rather than symbolic differentiation.
The Jacobian functions are constructed to be in-place with pre-allocated caches.
Functions provided by the user will be converted into in-place functions with pre-allocated caches.
Calls to nonlinear equation solver nlsolve are accelerated by exploiting sparsity in Jacobians with SparseDiffTools.jl.
Calculation of Jacobians of $\mu$, $\xi$, and $\mathcal{V}$ with automatic differentiation is accelerated by exploiting sparsity with SparseDiffTools.jl.

See the Example for how to use the type. To compare this package's speed with the original MATLAB code, run the wac_disaster.jl or rbc_cc.jl scripts. These scripts assess how long it takes to calculate a risk-adjusted linearization using the two numerical algorithms implemented by this package and by the original authors. The relaxation algorithm is generally around 50x-100x faster while the homotopy algorithm 3x-4x times faster.

The docstring for the constructor is given below. Some of the keywords allow the user to exploit sparsity in the objects comprising a risk-adjusted linearization. For more details on sparse methods, see Spares Arrays and Jacobians.

RiskAdjustedLinearizations.RiskAdjustedLinearization — Type

RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, z, y, Ψ, Nε)
RiskAdjustedLinearization(nonlinear_system, linearized_system, z, y, Ψ, Nz, Ny, Nε)

Creates a first-order perturbation around the stochastic steady state of a discrete-time dynamic economic model.

The first method is the main constructor most users will want, while the second method is the default constructor.

Inputs for First Method

μ::Function: expected state transition function
ξ::Function: nonlinear terms of the expectational equations
ccgf::Function: conditional cumulant generating function of the exogenous shocks
Λ::Function or Λ::AbstractMatrix: function or matrix mapping endogenous risk into state transition equations
Σ::Function or Σ::AbstractMatrix: function or matrix mapping exogenous risk into state transition equations
Γ₅::AbstractMatrix{<: Number}: coefficient matrix on one-period ahead expectation of state variables
Γ₆::AbstractMatrix{<: Number}: coefficient matrix on one-period ahead expectation of jump variables
z::AbstractVector{<: Number}: state variables in stochastic steady state
y::AbstractVector{<: Number}: jump variables in stochastic steady state
Ψ::AbstractMatrix{<: Number}: matrix linking deviations in states to deviations in jumps, i.e. $y_t - y = \Psi(z_t - z)$.
Nε::Int: number of exogenous shocks

Keywords for First Method

sss_vector_cache_init::Function = dims -> Vector{T}(undef, dims): initializer for the cache of steady state vectors.
Λ_cache_init::Function = dims -> Matrix{T}(undef, dims): initializer for the cache of Λ
Σ_cache_init::Function = dims -> Matrix{T}(undef, dims): initializer for the cache of Λ
jacobian_cache_init::Function = dims -> Matrix{T}(undef, dims): initializer for the cache of the Jacobians of μ, ξ, and 𝒱.
jump_dependent_shock_matrices::Bool = false: if true, Λ and Σ are treated as Λ(z, y) and Σ(z, y) to allow dependence on jumps.
sparse_jacobian::Vector{Symbol} = Symbol[]: pass the symbols :μ, :ξ, and/or :𝒱 to declare that the Jacobians of these functions are sparse and should be differentiated using sparse methods from SparseDiffTools.jl
sparsity::AbstractDict = Dict{Symbol, Mtarix}(): a dictionary for declaring the sparsity patterns of the Jacobians of μ, ξ, and 𝒱. The relevant keys are :μz, :μy, :ξz, :ξy, and :J𝒱.
colorvec::AbstractDict = Dict{Symbol, Vector{Int}}(): a dictionary for declaring the the matrix coloring vector. The relevant keys are :μz, :μy, :ξz, :ξy, and :J𝒱.
sparsity_detection::Bool = false: if true, use SparseDiffTools to determine the sparsity pattern. When false (default), the sparsity pattern is estimated by differentiating the Jacobian once with ForwardDiff and assuming any zeros in the calculated Jacobian are supposed to be zeros.

Inputs for Second Method

nonlinear_system::RALNonlinearSystem
linearized_system::RALLinearizedSystem
z::AbstractVector{<: Number}: state variables in stochastic steady state
y::AbstractVector{<: Number}: jump variables in stochastic steady state
Ψ::AbstractMatrix{<: Number}: matrix linking deviations in states to deviations in jumps, i.e. $y_t - y = \Psi(z_t - z)$.
Nz::Int: number of state variables
Ny::Int: number of jump variables
Nε::Int: number of exogenous shocks

source

Helper Types

To organize the functions comprisng a risk-adjusted linearization, we create two helper types, RALNonlinearSystem and RALLinearizedSystem. The first type holds the $\mu$, $\Lambda$, $\Sigma$, $\xi$, and $\mathcal{V}$ functions while the second type holds the $\mu_z$, $\mu_y$, $\xi_z$, $\xi_y$, $J\mathcal{V}$, $\Gamma_5$, and $\Gamma_6$ quantities. The RALNonlinearSystem type holds potentially nonlinear functions, and in particular $\mu$, $\xi$, and $\mathcal{V}$, which need to be linearized (e.g. by automatic differentiation). The RALLinearizedSystem holds both matrices that are only relevant once the model is linearized, such as $\Gamma_1$ (calculated by $\mu_z$), as well as $\Gamma_5$ and $\Gamma_6$ since these latter two quantities are always constant matrices.