Linear perturbation around steady state


The main functions of this section are in the folder 5_LinearizationFunctions.

The model is linearized with respect to aggregate variables. For this, we write the equilibrium conditions in the form of $F(X,X')=0$, where $X$ and $X'$ are (expected) deviations from steady state in two successive periods. Applying the total differential yields $A*X' = - B*X$, where $A$,$B$ are the first derivatives of $F$ with respect to $X'$,$X$. In the standard setting, we use the generalized Schur decomposition [Klein] to transform this equation into a linearized observation equation $d = gx*k$ and a linearized state transition equation $k' = hx*k$, where $k$ is a vector of the state variables and $d$ is a vector of the control variables ($X = \begin{bmatrix} k \\ d \end{bmatrix}$).

In our code, $F$ is implemented as HANKEstim.Fsys(), while differentiating and solving for $gx$ and $hx$ is done in HANKEstim.SGU(), and linearize_full_model() returns the results as a struct LinearResults:


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


struct LinearResults, containing

  • A::Array{Float64,2},B::Array{Float64,2}: first derivatives of Fsys() with respect to arguments X [B] and XPrime [A]
  • State2Control::Array{Float64,2}: observation equation
  • LOMstate::Array{Float64,2}: state transition equation

Overview of SGU()

SGU(sr, m_par, A, B; estim)

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


  • sr: steady-state structure (variable values, indexes, numerical parameters, ...)
  • A,B: derivative of Fsys() with respect to arguments X [B] and XPrime [A]
  • m_par: model parameters


  • gx,hx: observation equations [gx] and state transition equations [hx]
  • alarm_sgu,nk: alarm_sgu=true when solving algorithm fails, nk number of predetermined variables
  • A,B: first derivatives of Fsys() with respect to arguments X [B] and XPrime [A]

The function executes the following steps:

  • generate devices to retrieve distribution and marginal value functions from compressed states/controls (Γ and DC,IDC)

  • calculate the first derivative of HANKEstim.Fsys() with respect to X and XPrime. We use automatic differentiation (implemented in Julia by the package ForwardDiff). Partial derivatives are calculated using the ForwardDiff.jacobian() function. We exploit that some partial derivatives have known values (contemporaneous marginal value functions and the future marginal distributions) and set them directly instead of calculating them [BL].

  • compute linear observation and state transition equations using the HANKEstim.SolveDiffEq() function

Overview of `SolveDiffEq()'

SolveDiffEq(A, B, n_par; estim)

Calculate the solution to the linearized difference equations defined as P'BP xt = P'AP x{t+1}, where P is the (ntotal x r) semi-unitary model reduction matrix n_par.PRightAll of potentially reduced rank r.


  • A,B: matrices with first derivatives
  • n_par::NumericalParameters: n_par.sol_algo determines the solution algorithm, options are:
    • litx: Linear time iteration (implementation follows Reiter)
    • schur: Klein's algorithm (preferable if number of controls is small)


  • gx,hx: observation equations [gx] and state transition equations [hx]
  • alarm_sgu,nk: alarm_sgu=true when solving algorithm fails, nk number of predetermined variables
  • compute linear observation and state transition equations. The solution algorithm is set in n_par.sol_algo, with the options :schur (mentioned above) and :litx [lit]. The results are matrices that map contemporaneous states to controls [gx], or contemporaneous states to future states [hx]

Overview of Fsys()

Fsys(X, XPrime, Xss, m_par, n_par, indexes, Γ, compressionIndexes, DC, IDC, DCD, IDCD)

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

Split computation into Aggregate Part, handled by Fsys_agg(), and Heterogeneous Agent Part.


  • X,XPrime: deviations from steady state in periods t [X] and t+1 [XPrime]
  • Xss: states and controls in steady state
  • Γ, DC, IDC, DCD,IDCD: transformation matrices to retrieve marginal distributions [Γ], marginal value functions [DC,IDC], and the (linear) interpolant of the copula [DCD,IDCD] from deviations
  • indexes,compressionIndexes: access Xss by variable names (DCT coefficients of compressed $V_m$ and $V_k$ in case of compressionIndexes)


julia> # Solve for steady state, construct Γ,DC,IDC as in SGU()
julia> Fsys(zeros(ntotal),zeros(ntotal),XSS,m_par,n_par,indexes,Γ,compressionIndexes,DC,IDC)
*ntotal*-element Array{Float64,1}:

The function HANKEstim.Fsys() proceeds in the following way:

  1. set up vector F, that contains the errors to all equilibrium conditions. There are as many conditions as deviations from steady state (length of X,XPrime), and conditions are indexed with respective model variable in IndexStruct indexes
  2. generate locally all aggregate variables (for both periods) using @generate_equations
  3. 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 also the perturbed function is a copula)
  4. write all equilibrium condition-errors with respect to aggregate variables to F, using HANKEstim.Fsys_agg()
  5. compute optimal policies with HANKEstim.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 assumed present marginal value functions (in terms of their compressed deviation from steady state) as equilibrium condition-errors (backward iteration of the value function)
  6. compute future marginal distributions and the copula (on the copula grid) from previous distribution and optimal asset policies. Interpolate when necessary. Set difference to assumed future marginal distributions and copula values on the copula nodes as equilibrium condition-errors (forward iteration of the distribution)
  7. compute distribution summary statistics with HANKEstim.distrSummaries() and write equilibrium conditions with their respective (control) variables
  8. return F

Note that the copula is treated as the sum of two interpolants. An interpolant based on the steady-state distribution using the full steady-state marginals as a grid and a "deviations"-function that is defined on the copula grid generated in prepare_linearization(). The actual interpolation is carried out with HANKEstim.myinterpolate3(). Default setting is trilinear interpolation, the code also allows for 3d-Akima interpolation.

Called functions / macros


Write out the expansions around steady state for all variables in aggr_names, i.e. generate code that reads aggregate states/controls from steady state deviations.

Equations take the form of (with variable r as example):

  • r = exp.(Xss[indexes.rSS] .+ X[indexes.r])
  • rPrime = exp.(Xss[indexes.rSS] .+ XPrime[indexes.r])


(module) global aggr_names

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 SGU_estim()), variable-vectors X,XPrime, and Xss only contain the aggregate variables of the model.

myinterpolate3(xgrd1, xgrd2, xgrd3, ygrd, xeval1, xeval2, xeval3)

Trilineary project ygrd on (xgrd1,xgrd2,xgrd3) and use it to interpolate value at (xeval1,xeval2,xeval3).


julia> xgrd = [1.0,6.0];
julia> f((x,y,z)) = x+y+z;
julia> ygrd = f.(collect(Iterators.product(xgrid,xgrid,xgrid));
julia> xeval = [3.0,5.0];
julia> mylinearinterpolate3(xgrd,xgrd,xgrd,ygrd,xeval,xeval,xeval)
2x2x2 Array{Float64,3}:
[:,:,1] =
 9.0 11.0
11.0 13.0
[:,:,2] =
11.0 13.0
13.0 15.0