Block Jacobian Storage for DSGE

General Usage
,
1

Summary

When working with large system matrices, it is often convenient to express them in block form; moreover, since I consider a particular class of dynamic models, this system often yields a block Toeplitz form which is handy for sparse storage of system Jacobians.

Motivation

Consider the class of discrete time DSGE models with targets Y (these are the typically zeroed out), exogenous shocks Z, and endogenous states X in the following form:

F(X,Z) = Y

In most modeling languages, this system is treated as a single function call in order to calculate derivatives. However, an efficient modification of this paradigm was proposed in (Auclert et al. 2021) using a modular structure that leverages variable substitution via the model structure to reduce the dimensionality.

RBC Model Example

To justify this approach, consider the typical RBC model with unknowns C, K, Z, an exogenous shock to technology ε, and targets euler, goods_mkt, and shock_res.

The following is implemented in my modeling language DynamicMacroeconomics which is still a huge work in progress.

@simple function market_clearing(C, Y, I, r, β, γ)
    euler = (C ^ -γ) - (β * (1 + lead(r)) * (lead(C) ^ -γ))
    goods_mkt = Y - C - I
    return euler, goods_mkt
end

@simple function firms(Z, K, α, δ)
    r = α * Z * lag(K) ^ (α - 1) - δ
    Y = Z * lag(K) ^ α
    return r, Y
end

@simple function households(K, δ)
    I = K - (1 - δ) * lag(K)
    return I
end

@simple function shocks(Z, ρ, ε)
    shock_res = log(Z) - ρ * log(lag(Z)) - ε
    return shock_res
end

Here, we notice the intermediate states r, Y, and I need not be part of the Jacobian as long as their partials are accumulated via chain rules.

This not only reduces the dimension of the endogenous states U \subset X, but also implies faster AD for simple blocks (since it can be evaluated in the recursive form). In these simple cases, Jacobians are block Toeplitz and must be chain ruled according to each named block. Since these objects are recycled and interact, it would be beneficial to create a consistent container.

Existing Approach

To start sequence_jacobian defines a custom dictionary type which overloads matrix multiplication to occur consistently across named dimensions.

Pros Cons
getattr/setattr methods are super convenient Not an array type
Before conversion to matrix, this is technically sparse Does not admit a sparse matrix during conversion
Can store either matrix or AccumulatedDerivative Designed only to work in the sequence space

This is currently the approach that I use, albeit in a very unsophisticated way. Given a Jacobian for a single block M, I define the following constructor:

function SparseJacobian(M::AbstractMatrix, unknowns, targets)
    jac_dict = Dict(
        target => Dict(
            unknown => OffsetArrays.centered(
                M[o, (3 * i - 2) : (3 * i)]
            ) for (i, unknown) in enumerate(unknowns)
        ) for (o, target) in enumerate(targets)
    )
    return SparseJacobian(jac_dict, unknowns, targets)
end

Since I would rather nix this approach, I will spare further details on this post. If you are curious, I have the full implementation (with the modeling language) hosted here.

Related Packages

  • OffsetArrays.jl for time indexing (although I could write a custom sparse array to handle this).
  • ComponentArrays.jl for naming and casting to block matrix, but this may be a thing of the past.
  • ToeplitzMatrices.jl which serves as a starting point for efficient storage of Toeplitz operators.

Questions

To facilitate discussion, I’ll ask the following questions:

  • How would you handle this problem?
  • Why haven’t Toeplitz operators been better leveraged for dynamic modeling?
  • Are there better standards for discrete time forward/backward looking systems? It seems largely ignored besides some work done years ago on DifferentiableStateSpaceModels.jl