Perturbation Analysis

Computational Economics (ECO309)

Pablo Winant

2025-05-28

Deriving First Order Conditions

Lagrangian

We maximize the lagrangian to get:

$\begin{array}{rclr} \forall t \geq 0 & \frac{\partial L}{\partial i_{t}} & = & 0 \\ \frac{\partial L}{\partial c_{t}} & = & 0 \\ \frac{\partial L}{\partial k_{t + 1}} & = & 0 \end{array}$

It is important to note that we don’t differentiate with respect to a predetermined state
- check that you don’t differentiate w.r.t. $k_{0}$
It looks like the first order condition added four new variables $μ_{t}$ , $ν_{t}$ , $λ_{t}$ , $q_{t}$

Luckily these variables are associated to slackness conditions.

$μ_{t} \geq 0$	$i_{t} \geq 0$
$ν_{t} \geq 0$	$\exp (z_{t}) k_{t}^{α} - i_{t} \geq 0$
$q_{t} \geq 0$	$(1 - δ) k_{t} + i_{t} - k_{t + 1} \geq 0$
$λ_{t} \geq 0$	$\exp (z_{t}) k_{t}^{α} - i_{t} - c_{t} = 0$

The Karush-Kuhn-Tucker states, that for each slackness condition, at any time, either
- the lagrangian is 0 and it disappears from the F.O.C.s
- or it is not 0 and the associated constraint adds another condition

Eliminating constraints

$μ_{t} \geq 0$	$i_{t} \geq 0$
$ν_{t} \geq 0$	$\exp (z_{t}) k_{t}^{α} - i_{t} \geq 0$
$q_{t}$	$(1 - δ) k_{t} + i_{t} - k_{t + 1} = 0$
$λ_{t}$	$\exp (z_{t}) k_{t}^{α} - i_{t} - c_{t} = 0$

In general slackness conditions can be occasionally binding
For perturbation analysis, we need constraints to be always (or never binding)

Let’s review them:

$ν_{t}$ : it is equivalent to $c_{t} \geq 0$
- we necessarily have $c_{t} > 0$ since $U^{'} (0) = \infty$
- hence $ν_{t} = 0$
$μ_{t}$ : it states $k_{t + 1} \geq 0$
- if $k_{t + 1} = 0$ , then $c_{t + 1}$ . We can conclude $k_{t + 1} > 0$
- hence $μ_{t} = 0$
for multipliers associated to an equality constraint, we always keep the system
- multiplier can have any sign
inequality formulation is sometimes found too:
- $c_{t} \leq \exp (z_{t}) k_{t}^{α} - i_{t}$ ( you can destroy production insead of eating or investing it)
- $k_{t + 1} \leq (1 - δ) k_{t} + i_{t}$ (you can destroy capital instead of investing it)

Perturbation Analysis

Write all variables in deviation to the steady-state: $z_{t} = \overset{―}{z} + Δ z_{t}$ $k_{t} = \overset{―}{k} + Δ k_{t}$ $i_{t} = \overset{―}{i} + Δ i_{t}$ $c_{t} = \overset{―}{c} + Δ c_{t}$
- Remark: some smart economists use log-deviations (i.e. $x_{t} = \overset{―}{x} {\hat{x}}_{t}$ to make computations easier)

Replace in the system $β [\frac{{(\overset{―}{c} + Δ c_{t + 1})}^{- γ}}{{(\overset{―}{c} + Δ c_{t})}^{- γ}} ((1 - δ + α (\overset{―}{k} + Δ k_{t + 1})^{α - 1}))] = 1$ $\overset{―}{c} + Δ c_{t} = (\overset{―}{k} + Δ k_{t})^{α} - \overset{―}{i} - Δ i_{t}$ $\overset{―}{k} + Δ k_{t} = (1 - δ) (\overset{―}{k} + Δ k_{t - 1}) + \overset{―}{i} + Δ i_{t - 1}$ $\overset{―}{z} + Δ z_{t} = \overset{―}{z} + Δ ρ z_{t - 1}$
Differentiate…
(if we want to limit the number of equations, we can replace $c_{t}$ by its value)

Finite Differences

Choose small $ϵ > 0$ , typically $\sqrt{machine eps}$
- check eps()
Forward Difference scheme:
- $f^{'} (x) \approx \frac{f (x + ϵ) - f (x)}{ϵ}$
- precision: $O (ϵ)$
- bonus: if $f (x + ϵ)$ can compute $f (x) - f (x - ϵ)$ instead (Backward)
Central Difference scheme:
- Finite $f^{'} (x) \approx \frac{f (x + ϵ) - f (x - ϵ)}{2 ϵ}$
- average of forward and backward
- precision: $O (ϵ^{2})$

Two packages FiniteDiff and FiniteDifferences.

Example:

using FiniteDiff

g(x) = [x[1]^2 x[2]^3]
FiniteDiff.finite_difference_jacobian(g, [0.1, 0.2])

Symbolic Differentiation

manipulate the tree of algebraic expressions
- implements various simplification rules
requires mathematical expression
can produce mathematical insights
sometimes inaccurate:
- cf: ${(\frac{1 + u (x)}{1 + v (x)})}^{100}$

Two main libraries:

SymEngine.jl
- fast symbolic calculation
- mature C++ engine
Symbolics.jl:
- pure julia
- finite difference
- symbolic calculation

Example using Symbolics:

using Symbolics
@variables a x b y
eq = a*sin(x) + b*cos(y)
Symbolics.derivative(eq, x)

Automatic differentiation

There are many flavours of automatic differenation (check JuliaDiff.org)

Forward Accumulation Mode

isomorphic to dual number calculation
compute values and derivatives at the same time
efficient for $f : R^{n} \to R^{m}$ , with $n << m$
- (keeps lots of empty gradients when $n >> m$ )

Reverse Accumulation Mode

Reverse Accumulation / Back Propagation
- efficient for $f : R^{n} \to R^{m}$ , with $n >> m$
- requires data storage (to keep intermediate values)
- graph / example
Very good for machine learning:
- $\nabla_{θ} F (x; θ)$ where $F$ is a scalar objective

Solution

What is the solution of our problem?

At date $t$ controls must be chosen as a function of (predetermined) states
Mathematically speaking, the solution is a function $φ$ such that: $\forall t, x_{t} = φ (s_{t})$
Since the model is linear we look for un unknown matrix $X \in R^{n_{x}} \times R^{n_{s}}$ such that: $Δ x_{t} = X Δ s_{t}$

In the neoclassical model

The states are $k_{t}$ and $z_{t}$
The controls $i_{t}$ and $c_{t}$ must be a function of the states
- there is a decision rule $i ()$ , $c ()$ such that $i_{t} = i (z_{t}, k_{t})$ $c_{t} = c (z_{t}, k_{t})$
In the linearized model: $Δ i_{t} = i_{z} Δ z_{t} + i_{k} Δ k_{t}$ $Δ c_{t} = c_{z} Δ z_{t} + c_{k} Δ k_{t}$

Optimality condition:

Replacing in the system:

$Δ x_{t} = X Δ s_{t}$ $Δ s_{t + 1} = E Δ s_{t} + F X Δ s_{t}$ $Δ x_{t + 1} = X Δ s_{t + 1}$ $A Δ s_{t} + B Δ x_{t} + C Δ s_{t + 1} + D Δ x_{t + 1} = 0$

If we make the full substitution:

$((A + B X) + (D X + C) (E + F X)) Δ s_{t} = 0$

This must be true for all $s_{t}$ . We get the special Ricatti equation:

$(A + B X) + (D X + C) (E + F X) = 0$

This is a quadratic, matrix ( $X$ is 2 by 2 ) equation:

requires special solution method
there are multiple solutions: which should we choose?
today: linear time iteration selects only one solution
- alternative: eigenvalues analysis

Linear Time Iteration (2)

We get the equation: $\begin{array}{rcl} F (X, \tilde{X}) & = & (A + B X) + (C + D \tilde{X}) (E + F X) \\ = & 0 \end{array}$
Consider the linear time iteration algorithm
When the model is well-specified it is guaranteed to converge to the right solution.
- cf linear time iteration by Pontus Rendahl (link)
There are simple criteria to check that the solution is right, and that the model is well specified
$T$ is the time iteration operator… for linear models
- it does forward iteration ( $X_{t}$ as a function of $X_{t + 1}$ )

Algorithm:

choose stopping criteria: $ϵ_{0}$ and $η_{0}$
choose random $X_{0}$
given $X_{n}$ :
- compute $X_{n + 1}$ such that $F (X_{n + 1}, X_{n}) = 0$ $(B + (C + D X_{n}) F) X_{n + 1} + A + (C + D X_{n}) E = 0$ $X_{n + 1} = - (B + (C + D X_{n}) F)^{- 1} (A + (C + D X_{n}) E)$ $X_{n + 1} = T (X_{n})$
- compute:
  - $η_{n} = | X_{n + 1} - X_{n} |$
  - $ϵ_{n} = F (X_{n + 1}, X_{n + 1})$
- if $η_{n} < η_{0}$ and $ϵ_{n} < ϵ_{0}$
  - stop and return $X_{n + 1}$
  - otherwise iterate with $X_{n + 1}$

Simulating the model

Suppose we have found the solution $Δ x_{t} = X Δ s_{t}$
Recall the transition equation: $Δ s_{t + 1} = F Δ s_{t} + G Δ x_{t}$
We can now compute the model evolution following initial deviation in the state: $Δ s_{t} = \underset{P}{\underset{⏟}{(F + G X)}} Δ s_{t - 1}$
- $P$ is the simulation operator
- it is a backward operator (TODO: example of a reaction to a shock)
The system is stable if the biggest eigenvalue of $P$ is smaller than one…
… or if its spectral radius is smaller than 1: $ρ (P) < 1$
This condition is called backward stability
- it rules out explosive solutions
- if $ρ (P) > 1$ one can always find $s_{0}$ such that the model simulation diverges

Spectral radius

How do you compute the spectral radius of matrix P?
- naive approach: compute all eigenvalues, check the value of the biggest one…
```
using LinearAlgebra
M = rand(2,2)
maximum(abs, eigvals(M))
```
- becomes limited when size of matrix growth
- another approach: power iteration method
Power iteration method:
- works for matrices and linear operators
- take a linear operator $L$ over a Banach Space $B$ (vector space with a norm)
- use the fact that for most $u_{0} \in B$ , $\frac{| L^{n + 1} u_{0} |}{| L^{n} u_{0} |} \to ρ (L)$

Algorithm:

choose tolerance criterium: $η > 0$
choose random initial $x_{0}$ and define $u_{0} = \frac{x_{0}}{| x_{0} |}$
- by construction: $| u_{0} | = 1$
given $u_{n}$ , compute
- $x_{n + 1} = L . u_{n}$
- $u_{n + 1} = \frac{x_{n + 1}}{| x_{n + 1} |}$
- compute $η_{n + 1} = | u_{n + 1} - u_{n} |$
- if $η_{n + 1} < η$ :
  - stop and return $| x_{n + 1} |$
  - else iterate with $u_{n + 1}$

Differentials

Consider a Banach Space $B$ .
Consider an operator (i.e. a function): $T : B \to B$ .
Consider $\overset{―}{x} \in B$ .
$T$ is differentiable at $\overset{―}{x}$ if there exists a bounded linear operator $L \in L (B)$ such that: $T (x) = \overset{―}{x} + L . (x - \overset{―}{x}) + o (| x - \overset{―}{x} |)$
- when it exists we denote this operator by $T^{'} (\overset{―}{x})$
Remarks:
- Bounded operator means: $sup_{| x | = 1} | L . x | < + \infty$
- This definition of a derivative is usually referred to as Fréchet-derivative

The derivative of the time-iteration operator

Implicit relation can be differentiated: $F_{X}^{'} (T (X), X) T^{'} (X) + F_{Y}^{'} (T (X), X) = 0$
$F_{X}^{'} (T (X), X)$ being a regular matrix, it is (conceptually) easy to invert: $T^{'} (X) = - (F_{X}^{'} (T (X), X))^{- 1} F_{Y}^{'} (T (X), X)$
Finally, we get the explicit formula for the linear operator $T^{'}$ computed at the steady state: $T^{'} (\overset{―}{X}) . u = ((B + (C + D \overset{―}{X}) F)^{- 1}) D u (E + F \overset{―}{X})$
We can compute the spectral radius of $T^{'} (\overset{―}{X})$ using the power iteration method

Recap

We compute the derivatives of the model
Time iteration algorithm, starting from an initial guess $X_{0}$ and we repeat until convergence: $X_{n + 1} = (B + (C + D X_{n}) F)^{- 1} (A + (C + D X_{n}) E)$
We compute the spectral radius of two operators to ensure the model is well defined and that the solution is the right one.
backward stability: derivative of simulation operator $ρ (F + H \overset{―}{X})$
forward stability: derivative of time iteration operator $ρ (u \mapsto ((B + (C + D \overset{―}{X}) F)^{- 1}) D u (E + F \overset{―}{X}))$

Perturbation Analysis Computational Economics (ECO309) Pablo Winant 2025-05-28

Perturbation Analysis
Today
Deriving First Order Conditions
Neoclassical Growth Model
Deriving First Order Conditions
Lagrangian
Lagrangian
Eliminating constraints
First order model
Steady-State
Perturbation Analysis
Result:
Computing Derivatives
Julia Console
Julia package ecosystem
How do you install pckages?
Main approaches
Manual Differentiation
Finite Differences
Finite Differences: Higher order
Symbolic Differentiation
Automatic Differentiation
Automatic source code rewrite
Dual numbers: operator overloading
Dual numbers
Automatic differentiation
Libraries for AutoDiff
Solution using Linear Time iteration
Our problem
Solution
Optimality condition:
Linear Time Iteration
Linear Time Iteration (2)
Simulating the model
Spectral radius
Stability of the backward operator
Differentials
Examples of linear operators
Back to the time iteration operator
The derivative of the time-iteration operator
Recap