# 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`

:

`HANKEstim.linearize_full_model`

— Function`linearize_full_model()`

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

.

**Returns**

`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()`

`HANKEstim.SGU`

— Function`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`

**Arguments**

`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

**Returns**

`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()'

`HANKEstim.SolveDiffEq`

— Function`SolveDiffEq(A, B, n_par; estim)`

Calculate the solution to the linearized difference equations defined as P'*B*P x*t = 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.**Arguments**

`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)

**Returns**

`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()`

`HANKEstim.Fsys`

— Function`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*.

**Arguments**

`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`

)

**Example**

```
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}:
0.0
0.0
...
0.0
```

The function `HANKEstim.Fsys()`

proceeds in the following way:

- 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`

- generate locally all aggregate variables (for both periods) using
`@generate_equations`

- 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)
- write all equilibrium condition-errors with respect to
*aggregate*variables to`F`

, using`HANKEstim.Fsys_agg()`

- 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*) - 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*) - compute distribution summary statistics with
`HANKEstim.distrSummaries()`

and write equilibrium conditions with their respective (control) variables - 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

`HANKEstim.@generate_equations`

— Macro`@generate_equations()`

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])`

**Requires**

(module) global `aggr_names`

`HANKEstim.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 `SGU_estim()`

), variable-vectors `X`

,`XPrime`

, and `Xss`

only contain the aggregate variables of the model.

`HANKEstim.myinterpolate3`

— Function`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`

).

**Example**

```
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
```

- KleinSee the paper Using the generalized Schur form to solve a multivariate linear rational expectations model by Paul Klein (JEDC 2000)
- BLContemporaneous marginal value functions are irrelevant for optimal decisions, so its effect on other model variables is 0. Due to a rich enough set of prices, the future distribution directly only affects the Fokker-Planck equation. For details, see the paper Solving heterogeneous agent models in discrete time with many idiosyncratic states by perturbation methods
- 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`

-algorithm solves this equation for $Dg$ and $Dh$ iteratively.