Perturbation Method

This is a technical documentation page, one should read the tutorial to get started.

The toolbox has support for working with perturbation methods as summarized in [Ahn et al., 2017]. For the people who are coming from the DYNARE/RISE background, one can consider this as the pre-parser portion of analyzing the model files. Unfortunately, the differentiations are coded in DYNARE/RISE does not generalize to problems at the scale of heterogeneous agent models (see MATLABAutoDiff ). In this computation, the steady-state values need to be found before the differentiation.

For most applications, one can follow the recipe of

  1. Population var2ss, var2ischoice, and var2isshock

  2. Write the model equation sin setup_eqns (see Writing the setup_eqns Function ).

  3. Solve the model with some linear system solve_perturbation

  4. Compute IRFs with simulate_irf

Class

class MFG.ADHelper

Functions

ADHelper.get_gensys_form()

Convert the equations into the GENSYS form

Parameters:eqns (in class) – eqns
Returns:[G0, G1, psi, pi]
\[G0\cdot y_{t+1} = G1\cdot y_t + psi \cdot \eta_t + pi \cdot \varepsilon_t\]

Hence, the standard variable names for GENSYS.
ADHelper.get_sspace_form()

From derivatives, get the state-space representation

Parameters:eqns (matrix in class) – derivative values from dynamics equations
Returns:[G1, psi, pi]

Note

Intuitively, it “inverts” G0 matrix. For usages, it does not hurt to treat it as a black-box, and most instances, one would not need to call this function directly.

See also

ADHelper.krylov_reduction(varargin)

Reduces the dimensionality of the state space using Krylov subspace method

Parameters:
  • G1 (in class) – G1
  • n_choice (in class) – n_choice
  • n_g (in class) – n_g
  • n_p (in class) – n_p
  • n_arnoldi (in class) – n_arnoldi
Returns:

(in class)

Note

This helper is written so that one would not have to worry about the consistency of the inputs for complex functions. Check out other functions and tutorials.

See also

References

ADHelper.schur_solver()

Solver based on Schur decomposition.

Note

This is based on Sims’s gensys, but uses pre-clean+ Schur decomposition rather than taking the QZ-decomposition directly. QZ-decomposition scales worse with dimensionality of the problem.

Parameters:
  • G1 (in class) – G1
  • psi (in class) – psi
  • pi (in class) – pi
  • n_choie (in class) – n_choice
Returns:

(in class)

Warning

This class is written because users do not have to worry about setting the inputs to this function. Hence, make sure to check other functions (or Tutorial: Using the Perturbation Helper ).

ADHelper.setup_eqns()

User-Defined function of model equations

Returns:[eqns]
  • eqns (matrix): derivative values from vertically stacked equations

Note

Warning

The expected output is the derivative matrix, NOT the dual number. Hence, one would likely need to call getderivs function of the MATLABAutoDiff toolbox. For example

eqns = [eqn1; eqn2; eqn3];
eqns = getderivs(eqns);

See Tutorial: Using the Perturbation Helper for more details.

Todo

  • Consider changing syntax so that user would not have to call getderivs
ADHelper.set_var2advars_dvars()

Initialize automatic differentiation

Parameters:
Returns:

(set in class)

See also

References

ADHelper.set_var2loc()

Given the vector of positions, assign it locations

Parameters:var2ss (in class) – var2ss
Returns:(in class)

Note

This function puts variables in a nice order so that less reorganization is required, but users can just consider this as a blackbox.

ADHelper.simulate_irf(T, N, shock_loc)

Simulate impulse response functions

Parameters:
  • T (double) – time length for simulation
  • N (double) – number of steps
  • shock_loc (int) – which shock to compute impulse response of
Returns:

[output, time]

  • output (matrix): simulated path
  • time (vector): time vector

See also

ADHelper.solve_perturbation()

Solves perturbation solution

Parameters:
Returns:

(in class)

Example

This function contains many things. Refer to the Tutorial: Using the Perturbation Helper .

Properties

ADHelper.var2ss
Type:struct
Description:This struct set the variables.
ADHelper.var2ischoice
Type:struct of true/false
Description:Whether the variable is a choice/non-predetermined variable.
ADHelper.var2isshock
Type:struct of true/false
Description:Whether the variable is a shock/exogenous variable.
ADHelper.var2loc
Type:struct
Description:The variable location of the stacked variable
ADHelper.var2advars
Type:

struct
Description:

Struct containing the variables declared as a “dual” number. See https://github.com/sehyoun/MATLABAutoDiff .
Example:

One can get the dual number of consumption by

helper.var2advars.consumption
ADHelper.shock2loc
Type:struct
Description:A struct containing the index location of shocks. This is used with impact .
ADHelper.eqns
Type:matrix
Description:A matrix of the first order equation
ADHelper.n_vars
Type:integer
Description:Total number of variables
ADHelper.n_choice
Type:integer
Description:Number of choice/not-predetermined variables
ADHelper.n_choice_cached
Type:integer
Description:Cache location of for n_choice .
ADHelper.n_choice_red
Type:integer
Description:Number of choice variables after model reductions
ADHelper.n_g
Type:integer
Description:Number of predetermined variables
ADHelper.n_g_cached
Type:integer
Description:Cache location for n_g .
ADHelper.n_g_red
Type:integer
Description:Number of predetermined variables after model reductions
ADHelper.n_p
Type:

integer
Description:

Number of “redundant” variables. This is clear with an example, consider a dynamics given by

\[\begin{split}\frac{dx}{dt} = 2 y\\ y = 2x\end{split}\]

Though the equation has 2 variables, in fact this only has one “intrinsic” dimensionality. In this case, \(y\) can be solved out to get

\[\frac{dx}{dt} = 4 x\]

The class automatically cleans away these variables without dynamics \(\left(\frac{d}{dt}\right)\)-terms, and n_p refers to the number of these variables.

With small problems where QZ-decomposition can be used, a more general form of the problem can be setup, but in this class, it is assumed that the variables can be solved out directly, i.e., G1(redundant, redundant) is invertible.
ADHelper.n_pred
Type:integer
Description:Number of predetermined variables
ADHelper.n_arnoldi
Type:integer
Description:Number of arnoldi iteration step to take
ADHelper.is_cont_time
Type:true/false
Description:Whether a continuous time problem
ADHelper.state_reductions
Type:complicated
Description:Instruction for state reductions
Todo:Will be documented further
ADHelper.choice_reductions
Type:complicated
Description:Instructions for choice variable reductions.
ADHelper.G1
Type:

matrix
Description:

Unfortunately, G1 stands for different things at different point of computation. It can be

  • G1 matrix after G0 has been inverted
  • G1 dynamics matrix coming from the solution of schur_solver
Todo:

Separate out G1 so that it does not stand for different things at different time.
ADHelper.psi
Type:matrix
Description:\(\Psi\) in gensys notation.
ADHelper.pi
Type:matrix
Description:\(\Pi\) in gensys notation.
ADHelper.G1_cached
Type:matrix
Description:Cache location for G1
ADHelper.psi_cached
Type:matrix
Description:Cache location for psi
ADHelper.pi_cached
Type:matrix
Description:Cache location for pi
ADHelper.impact
Type:matrix
Description:impact matrix of gensys
ADHelper.eu
Type:

existence flags
Description:

existence flag

  • eu(1): existence (1 = solution exists)
  • eu(2): uniqueness (1 = solution is unique)
ADHelper.F
Type:matrix
Description:How solutions are linked
Todo:Will be documented
ADHelper.reduction
Type:matrix
Description:Reduction matrix which collects the instructions from state_reductions and choice_reductions
See also:inv_reduction
ADHelper.inv_reduction
Type:matrix
Description:Unwind the model reduction
See also:reduction

Citation Requests

If you use the method presented here, please consider [1] citing:

[Reiter, 2009] introduced the idea of the perturbation, and [Reiter, 2010] the model reduction. Hence, people coming from economics literature should consider citing them. However, the model reduction method as introduced in [Reiter, 2010] is numerically unstable. [2] Therefore, a stabilization needs to be made, and this is done by using [Freund, 2003]. I, jointly with Greg Kaplan, Benjamin Moll , Thomas Winberry , and Christian Wolf, put everything together in [Ahn et al., 2017]. Of course, none of this is possible if we can not take the derivatives in our life-time. [Rall, 1981] introduced the idea of automating differentiation, which is what allows us to differentiate within a reasonable time (which I found out via [Domke, 2009] . I requested the citation toward the blog since it explains the intuition clearly and concisely.).

Footnotes

[1]The citation requirements can balloon if we were to cite everything from the past, so use your discretion with the citations.
[2]The proof presented in the paper is based on the Cayley-Hamilton theorem, and hence does not apply when there are rounding errors.