Linear perturbation around steady state

Note

The routines described here live mainly in the submodule BASEforHANK.PerturbationSolution (src/SubModules/PerturbationSolution) and are called from top-level helpers such as linearize_full_model.

The model is linearized with respect to aggregate variables around the steady state. We collect all variables in period $t$ in a single vector $X_t$ and write all equilibrium conditions in implicit form

\[F(X_t, X_{t+1}) = 0,\]

where $X_t$ and $X_{t+1}$ contain (expected) deviations from steady state in successive periods. Taking the total differential at the steady state $X^*$ gives

\[A \, \Delta X_{t+1} + B \, \Delta X_t = 0,\]

with

\[A = \left.\frac{\partial F}{\partial X_{t+1}}\right\rvert_{X^*}, \qquad B = \left.\frac{\partial F}{\partial X_t}\right\rvert_{X^*}.\]

We adopt the SGU convention

\[X_t = \begin{bmatrix} d_t \\ k_t \end{bmatrix},\]

whereLinearSolution

$k_t$ is the vector of state / predetermined variables (including shocks),
$d_t$ is the vector of control / jump variables.

The linear solution computes matrices

$gx$: controls as a linear function of states,
$hx$: next-period states as a linear function of current states,

so that we obtain

\[d_t = gx \, k_t, \qquad k_{t+1} = hx \, k_t,\]

with shapes $gx \in \mathbb{R}^{n_d \times n_k}$, $hx \in \mathbb{R}^{n_k \times n_k}$ for $n_k$ states and $n_d$ controls.

In the code, the equilibrium errors $F$ are implemented as BASEforHANK.PerturbationSolution.Fsys. Differentiating and solving for $gx$ and $hx$ is done in BASEforHANK.PerturbationSolution.LinearSolution, which is called by linearize_full_model and returns the results in a LinearResults struct.

In the standard setting, we use the generalized Schur decomposition ^[Klein] to transform the system into the linearized observation equation $d_t = gx \, k_t$ and state transition equation $k_{t+1} = hx \, k_t$.

High-level entry point: `linearize_full_model`

BASEforHANK.linearize_full_model — Function

linearize_full_model(sr, m_par; ss_only = false)

Linearize the full model (i.e. including idiosyncratic states and controls) around the steady state, and solves using LinearSolution().

Arguments

sr::SteadyResults: Output of call_prepare_linearization()
m_par::ModelParameters
ss_only::Bool: If true, only the steady state is checked

Returns

LinearResults, containing

State2Control::Array{Float64,2}: Matrix of observation equation
LOMstate::Array{Float64,2}: Matrix of state transition equation
A::Array{Float64,2}: Jacobian of PerturbationSolution.Fsys() with respect to XPrime.
B::Array{Float64,2}: Jacobian of PerturbationSolution.Fsys() with respect to X.
SolutionError::Bool: indicator whether solution failed
nk::Int: Number of shocks

source

BASEforHANK.PerturbationSolution.LinearSolution — Function

LinearSolution(sr, m_par, A, B; allow_approx_sol=true, ss_only=false)

Calculate the linearized solution to the non-linear difference equations defined by function Fsys(), using Schmitt-Grohé & Uribe (JEDC 2004) style linearization (apply the implicit function theorem to obtain linear observation and state transition equations).

The Jacobian is calculated using the package ForwardDiff

This function computes the first-order derivatives of the equilibrium conditions with respect to states and controls, constructs the linear system, and solves it. It handles the idiosyncratic part of the Jacobian efficiently by exploiting known derivatives (see set_known_derivatives_distr!) and using automatic differentiation for the rest.

Arguments

sr: steady-state structure (variable values, indexes, numerical parameters, ...).
m_par: model parameters.
A,B: derivative of Fsys() with respect to arguments X [B] and XPrime [A]. Can be initialized as zero matrices.
allow_approx_sol: if true (default), the function will attempt to solve the linearized model even if the system is indeterminate (shifting the critical eigenvalues).
ss_only: if true, the function will only check if the steady state is a solution (i.e., if Fsys evaluates to zero at steady state) and return the residuals.

Returns

gx: observation equations matrix (mapping states to controls).
hx: state transition equations matrix (mapping states to future states).
alarm_LinearSolution: true if the solution algorithm fails or is indeterminate.
nk: number of predetermined variables (states).
A,B: first derivatives of Fsys() with respect to arguments X [B] and XPrime [A].

source

linearize_full_model is the user-facing entry point for the perturbation solution. It

allocates empty Jacobian matrices A, B of size sr.n_par.ntotal × sr.n_par.ntotal,
calls BASEforHANK.LinearSolution to
- check steady-state feasibility (option ss_only = true), or
- compute Jacobians and solve the linear model,
wraps the results in a LinearResults struct.

The LinearResults object contains, among other things:

gx: observation matrix mapping states to controls (nd × nk).
hx: transition matrix mapping current states to next states (nk × nk).
A, B: Jacobians of Fsys at steady state, with A = ∂F/∂X', B = ∂F/∂X.
nk: number of predetermined (state) variables; nd is inferred from the length of X.
alarm_LinearSolution: boolean indicator if the solution encountered borderline or approximate handling (e.g. indeterminacy).

Typical usage:

sr = call_prepare_linearization(m_par)
lr = linearize_full_model(sr, m_par)

The resulting lr is then used for impulse responses, variance decompositions, higher-order derivatives, and estimation.

Core linearization: `LinearSolution`

The function BASEforHANK.LinearSolution performs the actual first-order perturbation around the steady state. It executes the following main steps:

Prepare transformations

Construct the transformation elements for compression and decompression:

Γ (shuffle matrices for marginals),
DC / IDC (DCT forward / inverse for value functions),
DCD / IDCD (DCT forward / inverse for copula deviations).

These define the mapping between the compressed representation used in the state vector X and the full-grid distributions and value functions.

Compute Jacobians of Fsys

The equilibrium error function F is implemented as BASEforHANK.PerturbationSolution.Fsys. LinearSolution computes its Jacobians

\[A = \partial F/\partial X', \qquad B = \partial F/\partial X\]

using ForwardDiff.jacobian, while inserting known derivative blocks directly ^[BL]:

contemporaneous marginal value functions → identity block in B over indexes.valueFunction,
future marginal distributions / copula components → pre-filled via set_known_derivatives_distr! in A.

This reduces the dimension of the automatic differentiation problem and speeds up linearization substantially.

BASEforHANK.PerturbationSolution.set_known_derivatives_distr! — Function

set_known_derivatives_distr!(A, dist_indexes, transform_elements, indexes,
                             n_par, transition_type, transform_type)

Pre-fill known derivative blocks of the Jacobian matrix A (with respect to future states) for the distributional part of the state vector. This avoids redundant automatic differentiation and significantly speeds up linearization.

The function exploits the structure of the distribution transition equations: the Jacobian with respect to future distributions follows directly from the state transition rules (Γ matrices) and the distribution representation (copula, CDF, or representative agent).

Multiple methods handle different combinations of:

indexes: Type indicating distribution representation (CopulaTwoAssetsIndexes, CopulaOneAssetIndexes, CDFIndexes, RepAgentIndexes).
transition_type: LinearTransition or NonLinearTransition.
transform_type: LinearTransformation or ParetoTransformation.

Arguments

A::AbstractMatrix: Jacobian matrix (filled in-place).
dist_indexes::Vector: Indices in the state vector corresponding to distributional variables.
transform_elements::TransformationElements: Pre-computed transition operators (Γ matrices).
indexes: Index struct specifying the distribution representation.
n_par::NumericalParameters: Numerical parameters (grid sizes, etc.).
transition_type: Either LinearTransition or NonLinearTransition.
transform_type: Either LinearTransformation or ParetoTransformation.

Behavior

Sets A[dist_indexes, dist_indexes] to -I (negative identity) as the base, then modifies specific blocks corresponding to each asset/income dimension using the transition matrices Γ. For CDF-based states with linear transformation, handles both unit and shuffled derivatives via cumsum logic. For Pareto transformation, uses only the base -I block.

source

Solve the linear system

Given A and B, the function calls BASEforHANK.PerturbationSolution.SolveDiffEq to compute the linear observation and state transition equations. The solver supports the options :schur and :litx; both yield consistent gx, hx mappings.

Return solution

LinearSolution returns

gx, hx, alarm_LinearSolution, nk, A, B

and is typically only called indirectly through linearize_full_model.

Implementation details:

Variable ordering follows an SGU-style convention internally when constructing Jacobians for FO/SO/TO; helper functions convert to the ordering used by Levintal where required for nested AD.
Automatic differentiation uses chunking (ForwardDiff.Chunk) tuned for performance; steady-state feasibility is checked before differentiation.
In verbose mode, maximum residuals per block are printed to aid diagnostics.

Solving the linear system: `SolveDiffEq`

BASEforHANK.PerturbationSolution.SolveDiffEq — Function

SolveDiffEq(A, B, n_par; allow_approx_sol=true)

Calculate the solution to the linearized difference equations defined as

P'*B*P x_t = P'*A*P x_{t+1},

where P is the (ntotal x r) semi-unitary model reduction matrix n_par.PRightAll of potentially reduced rank r.

This function solves the system for the observation equations d_t = gx * k_t and state transition equations k_{t+1} = hx * k_t.

Arguments

A,B: Matrices with first derivatives of the equilibrium conditions with respect to XPrime (A) and X (B).
n_par: Numerical parameters structure (NumericalParameters). n_par.sol_algo determines the solution algorithm. Options are:
- :lit: Linear time iteration (implementation follows Reiter).
- :litx: Optimized Linear time iteration (Reiter + QR + Howard-style improvement).
- :schur: Klein's algorithm (based on generalized Schur decomposition, default). Preferable if the number of controls is small.
allow_approx_sol: Keyword, Bool (default true). If true, approximate solutions are allowed if the system is not unique and determinate (e.g., critical eigenvalues are shifted).

Returns

gx: Observation equations matrix (mapping states to controls).
hx: State transition equations matrix (mapping states to future states).
alarm_LinearSolution: Bool. true when the solving algorithm fails or the solution is indeterminate.
nk: Number of predetermined variables (states) found by the solver.

source

The function BASEforHANK.PerturbationSolution.SolveDiffEq solves the generalized linear difference equation

\[A \, \Delta X_{t+1} + B \, \Delta X_t = 0\]

to obtain gx and hx.

Optional model reduction

If model reduction is active, the system is pre- and post-multiplied by the reduction matrix $\mathcal{P}$ that is computed in model_reduction and stored in n_par.PRightAll. This transforms the full system into a reduced one while preserving its dynamics in the directions that matter for aggregate quantities.

Solution algorithms

SolveDiffEq can use two alternative algorithms:

:schur (default): generalized Schur decomposition ^[Klein] to partition stable / unstable roots and select the stable solution.
:litx: an algorithm based on the Implicit Function Theorem that iterates on the fixed-point conditions for the derivatives Dg, Dh ^[lit].

Both return matrices that map contemporaneous states to controls (gx) and to future states (hx).

Numerical robustness

If allow_approx_sol = true, critical eigenvalues can be slightly shifted to obtain a nearby determinate solution. In such cases alarm_LinearSolution is set to true so that the caller can monitor or discard such draws in estimation.

Equilibrium errors: `Fsys`

BASEforHANK.PerturbationSolution.Fsys — Function

Fsys(X, XPrime, XSS, m_par, n_par, indexes, compressionIndexes, transform_elements; only_F = true)

Equilibrium error function: returns deviations from equilibrium around steady state.

It splits the computation into an Aggregate Part, handled by Fsys_agg(), and a Heterogeneous Agent Part (backward iteration of value functions and forward iteration of distributions).

Arguments

X,XPrime: Deviations from steady state in periods t [X] and t+1 [XPrime]
XSS: States and controls in steady state
m_par: Model parameters
n_par: Numerical parameters
indexes: Index struct to access variables in X and XPrime
compressionIndexes: Indices of DCT coefficients selected for perturbation
transform_elements: A struct containing transformation matrices (Γ, DC, IDC, DCD, IDCD) used to map between compressed states/controls and full grids (distributions and value functions).
only_F::Bool: If true (default), returns only the error vector F. If false, returns additional objects useful for debugging or other computations.

Returns

If only_F = true:
- F: Vector of equilibrium errors.
If only_F = false:
- F: Vector of equilibrium errors.
- pf: Policy functions for the current iteration.
- vf_new: Updated value functions.
- tax_rev: Tax revenue.

source

The function BASEforHANK.PerturbationSolution.Fsys constructs and evaluates the vector of equilibrium errors

\[F(X_t, X_{t+1})\]

for given deviations from steady state X, XPrime.

It proceeds in the following steps:

Initialize residual vector

Set up a vector F that contains the residuals for all equilibrium conditions. There are as many conditions as entries in X / XPrime. Conditions are indexed by the IndexStruct indexes.

Generate aggregate variables

Reconstruct all aggregate variables (for both periods) from the entries in X and XPrime using BASEforHANK.Parsing.@generate_equations.

Rebuild distributions and value functions

Construct the full-grid marginal distributions, marginal value functions, and the copula from the steady-state values and the (compressed) deviations. For the copula, the selection of DCT coefficients that can be perturbed ensures that the perturbed function remains a copula.

Aggregate block

Write all equilibrium errors that depend only on aggregate variables into F using BASEforHANK.PerturbationSolution.Fsys_agg.

Backward iteration of the value function

Compute optimal policies with BASEforHANK.SteadyState.EGM_policyupdate, given future marginal value functions, prices, and individual incomes. Infer present marginal value functions from them (envelope theorem) and set the difference to the assumed present marginal value functions (in terms of their compressed deviation from steady state) as equilibrium errors.

Forward iteration of the distribution

Compute future marginal distributions and the copula (on the copula grid) from the previous-period distribution and optimal asset policies. Interpolate when necessary. Set the difference to the assumed future marginal distributions and copula values on the copula nodes as equilibrium errors.

Distribution summary statistics

Compute distribution summary statistics with BASEforHANK.SteadyState.distrSummaries and write the corresponding equilibrium conditions with their respective control variables.

Return

Return the residual vector F (and, internally, additional objects such as distributions and policy functions when requested).

Additional remarks:

The copula is treated as the sum of two interpolants: one based on thesteady-state distribution using the full steady-state marginals as agrid, and a “deviations” function that is defined on the copula gridgenerated in prepare_linearization. The interpolation is carried outwith BASEforHANK.Tools.myinterpolate3. By default, trilinearinterpolation is used; the code also allows for 3D Akima interpolation.
Copula deviations are stored and manipulated in compressed DCT coefficient space and reconstructed in CDF space on the copula grid for comparisons (see pack_distributions.jl).
Marginal distributions use Γ to respect the unit-sum constraint.

Shapes and storage conventions:

Copula deviations: stored in PDF space as DCT coefficients and reconstructed in CDF space for evaluation. Grids are nb_copula × nk_copula × nh_copula (two-asset) or nb_copula × nh_copula (one-asset).
Marginals: Γ maps N − 1 reduced coordinates to N full probabilities per margin, enforcing unit sum. Income h is discrete and handled as a PDF; it can be converted to a CDF under nonlinear transitions.
Value functions: compressed as deviations in log-inverse-marginal-utility space; uncompressed via IDC and exponentiated back to marginal utility with mutil.

Model reduction

BASEforHANK.PerturbationSolution.model_reduction — Function

model_reduction(sr, lr, m_par)

Produce Model Reduction based on Variance Covariance Matrix of States and Controls.

This function performs the second stage of model reduction. It computes the covariance of the states and controls using the first-stage solution (linearized model lr) and then uses Principal Component Analysis (PCA) via compute_reduction to find a lower-dimensional subspace that captures the most significant dynamics.

The reduction is based on the method described in Bayer, Born, and Luetticke (2024, AER). If n_par.further_compress is true, it updates the projection matrices in n_par to map the full state space to this reduced factor space.

Arguments

sr::SteadyResults: Steady state results containing initial indices and parameters.
lr::LinearResults: Linear solution from the first stage (DCT reduction only).
m_par::ModelParameters: Model parameters.

Returns

SteadyResults: A new SteadyResults struct containing:
- Updated indexes_r (indices for the reduced model).
- Updated n_par (numerical parameters with reduced dimensions and projection matrices).
- Other fields copied from the input sr.

source

BASEforHANK.PerturbationSolution.compute_reduction — Function

compute_reduction(sr, lr, m_par, shock_names)

Compute the second-stage model reduction using Principal Component Analysis (PCA) on the covariance matrices of the linearized model's states and controls.

This function identifies the linear combinations of the first-stage variables (DCT coefficients of marginal value functions and distributions) that contribute most to the model's volatility over the business cycle. It constructs projection matrices (PRightAll, PRightStates) to project the full system onto this reduced subspace, significantly reducing the number of state and control variables for estimation while maintaining approximation quality.

Arguments

sr: Steady state results containing indices and parameters.
lr: Initial linear solution (first-stage reduction) used to compute long-run covariances.
m_par: Model parameters (used for shock variances).
shock_names: Vector of symbols denoting the shocks to include in the covariance calculation.

Returns

indexes_r: Structure containing indices for the reduced model variables.
n_par: Updated numerical parameters object containing:
- Reduced dimensions (nstates_r, ncontrols_r, ntotal_r).
- Projection matrices (PRightAll, PRightStates) mapping full variables to reduced factors.

source

The function model_reduction derives an approximate factor representation from a first solution of the heterogeneous-agent model ^[BBL]. It stores the matrices that map factors to the full set of state and control variables.

Conceptually, it computes the long-run variance–covariance matrix of the state vector and constructs a projection $\mathcal{P}$ such that the reduced variables

\[X_r = \mathcal{P} X\]

capture most of the variation that matters for aggregate quantities. This enables fast re-solves with aggregate-only updates.

At a high level:

Use the linear solution (gx, hx) together with the shock processes in m_par to compute long-run covariances of states and controls.
Perform an eigenvalue decomposition and select eigenvectors associated with the dominant eigenvalues, based on thresholds in n_par.
Construct projection matrices PRightAll, PRightStates that map the original compressed coefficients to factors.
Store these matrices and reduced dimensions (nstates_r, ncontrols_r, ntotal_r) in n_par and return reduced index structures.

If n_par.further_compress == false, the model reduction step is skipped and the identity matrix is used as PRightAll / PRightStates.

Reduced linear system and `update_model`

BASEforHANK.PerturbationSolution.update_model — Function

update_model(sr, lr, m_par)

Efficiently update the linearized solution of the model when only aggregate parameters have changed.

This function re-calculates the Jacobian of the aggregate equilibrium conditions (Fsys_agg) and combines it with the previously computed Jacobian of the heterogeneous agent block (which remains unchanged). It then solves the updated linear system using LinearSolution_reduced_system. This provides a significant speedup during estimation loops where only aggregate parameters are varied.

Arguments

sr::SteadyResults: Steady state results containing the model structure and parameters.
lr::LinearResults: The previous linear solution results, containing the Jacobians A and B from the full model linearization.
m_par::ModelParameters: The new set of model parameters (with updated aggregate parameters).

Returns

LinearResults: A new LinearResults struct containing:
- State2Control: Updated observation matrix.
- LOMstate: Updated state transition matrix.
- A, B: Updated Jacobian matrices of the system.
- SolutionError: Boolean flag indicating if the solution failed.
- nk: Number of stable eigenvalues.

source

BASEforHANK.PerturbationSolution.Fsys_agg — Function

Fsys_agg(X, XPrime, XSS, distrSS, m_par, n_par, indexes)

Return deviations from aggregate equilibrium conditions.

indexes can be both IndexStruct or IndexStructAggr; in the latter case (which is how function is called by LinearSolution_reduced_system()), variable-vectors X,XPrime, and XSS only contain the aggregate variables of the model.

source

BASEforHANK.PerturbationSolution.LinearSolution_reduced_system — Function

LinearSolution_reduced_system(sr, m_par, A, B; allow_approx_sol=false)

Calculate the linearized solution to the non-linear difference equations, but only differentiating with respect to the aggregate part of the model Fsys_agg().

The partial derivatives of the Jacobian belonging to the heterogeneous agent part of the model are taken from the full-model derivatives provided as arguments, A and B (computed by LinearSolution()). This allows for faster updates during estimation when only aggregate parameters change.

Arguments

sr: Steady-state structure (SteadyResults) containing variable values, indexes, numerical parameters, etc.
m_par: Model parameters (ModelParameters).
A,B: Derivative matrices of the full Fsys() with respect to XPrime (A) and X (B). These are used as the base, and only the aggregate parts are updated.
allow_approx_sol: Keyword, Bool (default false). If true, the function will attempt to solve the linearized model even if the system is indeterminate.

Returns

gx: Observation equations matrix (mapping states to controls).
hx: State transition equations matrix (mapping states to future states).
alarm_LinearSolution: Bool. true if the solution algorithm fails.
nk: Number of predetermined variables (states).
A,B: The updated first derivatives of the system (aggregate parts re-computed, heterogeneous parts preserved).

source

During estimation, many parameter updates affect only aggregate parameters (e.g. shock processes, policy rules, some preference parameters). The idiosyncratic block (household problem, transition matrices) often remains unchanged. To exploit this, BASEforHANK provides an efficient re-linearization routine.

`LinearSolution_reduced_system`

LinearSolution_reduced_system(sr, m_par, A, B; allow_approx_sol = false)

Treats the idiosyncratic part of the Jacobians as fixed, reused from a previous full linearization.
Rebuilds only the aggregate part of Fsys and its derivatives by differentiating BASEforHANK.PerturbationSolution.Fsys_agg.
Updates the corresponding rows / columns of A and B.
Calls BASEforHANK.PerturbationSolution.SolveDiffEq on the updated system to obtain new gx, hx, and nk.

It returns updated gx, hx, A, B, and the alarm flag.

Reduced residual function

Conceptually, the reduced system can be written as

\[F_r(X_r, X'_r) = \mathcal{P}' F(\mathcal{P} X_r, \mathcal{P} X'_r),\]

where $\mathcal{P}$ is the projection from model reduction. In this representation:

blocks with known-zero derivatives (e.g. contemporaneous marginal value functions, future marginal distributions) are held constant and not differentiated,
reduced index sets sr.indexes_r track which variables are treated as having nontrivial derivatives,
ForwardDiff is applied only to those blocks when computing the reduced Jacobians A_r, B_r.

The reduced Jacobians are then mapped back into the full A, B by updating only the relevant rows and columns.

`update_model`

update_model is a high-level convenience wrapper used in estimation loops:

Start from steady-state results sr, an existing linear solution lr, and updated parameters m_par.
Call LinearSolution_reduced_system with the Jacobians lr.A, lr.B.
Return a new LinearResults object with updated State2Control, LOMstate, A, B, SolutionError, and nk.

As long as the parameter change only affects aggregate equations, this is much faster than running linearize_full_model again. If parameters change the idiosyncratic structure, a new full call to BASEforHANK.PerturbationSolution.LinearSolution (and optionally model_reduction) is required.

Higher-order solution (SO / TO)

BASEforHANK.PerturbationSolution.compute_derivatives — Function

compute_derivatives(sr, m_par, A, B; estim=false, only_SO=true, H_reduc=nothing, TO_reduc=nothing, fix_known_zeros=false)

Compute the derivatives of the equilibrium conditions F necessary for higher-order perturbation solutions.

This function performs the following steps:

Calculates the first-order Jacobian of the (reduced) system and solves for the linear policy functions gx and hx.
Computes the Hessian (second-order derivatives) of the system, potentially using a reduced-space approach for efficiency.
Optionally computes the third-order derivatives (tensor).

It accounts for model reduction by differentiating the system with respect to the reduced factors and then mapping back to the full space.

Arguments

sr: Steady-state results structure (SteadyResults).
m_par: Model parameters (ModelParameters).
A, B: Initial matrices for the Jacobians (often placeholders).
estim: Boolean (default false). Currently unused in the function body.
only_SO: Boolean (default true). If true, only derivatives up to the second order (Hessian) are computed. If false, third-order derivatives are also computed.
H_reduc: (Optional) Precomputed reduced Hessian matrix. If provided, skips calculation.
TO_reduc: (Optional) Precomputed reduced third-order derivative tensor.
fix_known_zeros: Boolean (default false). If true, utilizes knowledge of zero derivatives for certain distribution variables to speed up computation.
buildpath: String (default ""). Path to save intermediate results like Hessians and third-order tensors.

Returns

gx: Observation matrix (states to controls).
hx: State transition matrix.
alarm_LinearSolution: Boolean indicating if the linear solution failed.
nk: Number of stable eigenvalues.
A, B: First-order Jacobian matrices.
H: Second-order derivative matrix (Hessian).
TO: Third-order derivative matrix (only returned if only_SO=false).

source

BASEforHANK.PerturbationSolution.SolveSylvester — Function

SolveSylvester(gx, hx, B, A, H, n_par, indexes, covariances, shock_names; Xred=nothing, TestOfGoodness=false)

Solve a generalized Sylvester equation to compute second-order approximations of policy functions and risk adjustments.

This function computes the second-order derivatives of the policy functions (Gxx, Hxx) and the risk adjustment terms (Gσσ, Hσσ) for the observation and state transition equations. The notation and method follow Levintal (2017). The system to solve is of the form A*X + B*X*C + D = 0.

Arguments

gx: First-order derivatives of the observation equation (control variables).
hx: First-order derivatives of the state transition equation (state variables).
B, A: Jacobians of the equilibrium conditions (B w.r.t X, A w.r.t XPrime).
H: Hessian of the equilibrium conditions (second-order derivatives).
n_par: Numerical parameters structure.
indexes: Index structure for variables.
covariances: Covariance matrix of the shocks.
shock_names: Vector of shock names.
Xred: (Optional) Initial guess for the solution matrix X. If nothing, it is initialized automatically.
TestOfGoodness: (Optional) If true, computes and prints residuals to test the accuracy of the solution.

Returns

Gxx: Second-order derivatives of the observation equation (controls).
Hxx: Second-order derivatives of the state transition equation (states).
Gσσ: Second-order derivatives with respect to the volatility parameter σ (risk adjustment for controls).
Hσσ: Second-order derivatives with respect to the volatility parameter σ (risk adjustment for states).
Xred: The full solution matrix X (useful for warm-starting subsequent calls).

source

The higher-order solution builds on the linearization to capture local curvature (second order, SO) and optionally third-order effects (TO).

Warning

Solving for the third-order terms of the Taylor expansion is not yet implemented in the toolbox.

Higher-order solutions are needed to compute welfare, aggregate risk-effects, and non-linear impulse responses.

Overview:

FO: compute A, B and solve for gx, hx as described above.
SO: compute the reduced Hessian H_red via nested automatic differentiation on the reduced system, then reconstruct the full Hessian H using known-zero patterns.
TO (optional): compute the reduced third-order derivatives TO_red similarly, then fill the full third-order tensor TO.

Key steps and conventions:

Reduced system Automatic Differentiation: define the reduced residual function F_r(X_r, X'_r) = P' F(P X_r, P X'_r) and differentiate only with respect to variables that are not in constant-derivative blocks (tracked by reduced indexes). This concentrates Automatic Differentiation effort on parts of the model that move with aggregates.
Variable ordering: SO / TO follow an SGU-style ordering X = [controls′; controls; states′; states] for Jacobian / Hessian construction; helper mappings convert to the internal ordering used to solve the higher-order equations, following Levintal (2017).

Caching and reconstruction:

H_red.jld2 and TO_red.jld2 are cached to avoid recomputation during estimation loops. When present, they are loaded and used directly.
fill_hessian(ix_all, ix_const, length_X0, H_red) maps the reduced Hessian back to the full H using index sets of variables with constant derivatives.
fill_thirdDeriv(ix_all, ix_const, length_X0, TO_red) performs the analogous mapping for the third-order tensor TO.

Outputs:

FO: gx, hx, A, B, nk, and alarm_LinearSolution.
SO: additionally H (typically stored sparsely); dimension roughly (length_X0 + 1) × 4 (length_X0 + 1)^2.
TO (optional): additionally TO (typically stored sparsely); dimension roughly (length_X0 + 1) × (2 (length_X0 + 1))^3.

References

KleinPaul Klein (2000), Using the generalized Schur form to solve a multivariate linear rational expectations model, Journal of Economic Dynamics and Control.
BLContemporaneous marginal value functions are irrelevant for optimal individual decisions, so their effect on other model variables is zero. Due to a rich enough set of prices, the future distribution directly only affects the Fokker–Planck equation. For details, see Solving heterogeneous agent models in discrete time with many idiosyncratic states by perturbation methods, Quantitative Economics, Vol. 11(4), November 2020, pp. 1253–1288.
BBL Shocks, Frictions, and Inequality in US Business Cycles, American Economic Review, 2024.
litInvoking the Implicit Function Theorem, there exist functions $g$ and $h$ such that $F\left(\begin{pmatrix} k \\ g(k) \end{pmatrix},\begin{pmatrix} h(k) \\ g(h(k)) \end{pmatrix}\right)=0$. Totally differentiating by $k$ yields $B \begin{pmatrix}\mathbb{I}\\ Dg \end{pmatrix} + A \begin{pmatrix}\mathbb{I}\\ Dg \end{pmatrix} Dh = 0$. The :lit-type algorithms solve this equation for $Dg$ and $Dh$ iteratively.