Intro to DSGE Model & Julia Basics

Business Cycles and Fluctuations - AE2E6

Objectives

Introduce Julia environment
- run cells
- basic types: numbers, strings, vectors, matrices
- functions
- simple plots
Simulate AR1 models
- impulse response functions
- stochastic simulations
  - conditional / unconditional moments
- develop intuition about eigenvalues / ergodic distributions

Julia: Quick Introduction

Why Julia?

Developped at MIT on top of opensource technologies

linux / git / llvm

Syntax inspired by Matlab but:

more consistent
lots of features from high level languages

Everything is JIT-compiled

no interpreted vs compiled treadeoff -> very fast
most of the base library is written in Julia

Opensource/free + vibrant community

Numbers

# numbers (usual operations usual)
(1.0+(2.0+3.0*(4.0+5.0)))/30

1.0

# exponentials are denoted by ^
2^8

Variables / assignments / comparisons

# variable assignment
x = 10

# Variable names can have Unicode characters
# To get ϵ in the REPL, type \epsilon<TAB>
a = 20
σ = 34
ϵ = 1e-4

0.0001

# comparison 
2 == 3

false

3<=3

true

Strings

# Strings can also contain Unicode characters
fancy_str = "α is a string"

"α is a string"

# double quotes define a character, not a string
'c'

'c': ASCII/Unicode U+0063 (category Ll: Letter, lowercase)

# string interpolation (1/2)
he_who_must_not_be_named = "Voldemort"
"Welcome $(he_who_must_not_be_named)!"

"Welcome Voldemort!"

# string interpolation (2/2)
n = 1999999
println("Iteration ", n, " is still running...")

Iteration 1999999 is still running...

Arrays

# vectors
v = [1,2,3]

3-element Vector{Int64}:
 1
 2
 3

# matrices 
M = [1 2 3 ; 4 5 6 ; 7 8 9]

3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

# matrix multiplication
M*v

3-element Vector{Int64}:
 14
 32
 50

# mutating vectors
x = ["One"]
push!(x, "Two")
push!(x, "Three")
push!(x, "Four")
push!(x, "Five")

# Note how the type of the vector is different from above
# Vectors in Julia hold homoegenous types

# Also note the exclation mark in `push!`: it is a reminder
# of the fact that `push!` mutates its first argument

5-element Vector{String}:
 "One"
 "Two"
 "Three"
 "Four"
 "Five"

# access elements in a vector
v[1]

# access elements in a matrix
M[1,2]

Warning

Indexing in Julia is 1-based. First element of a collection is denoted by 1.

# slice matrices
println(M)

# keep first line
println("First line")
println(M[1,:])

# keep second column
println("Second column")
println(M[:,2])

# extract a submatrix
println(M[1:2,1:2])

[1 2 3; 4 5 6; 7 8 9]
First line
[1, 2, 3]
Second column
[2, 5, 8]
[1 2; 4 5]

# concatenate vectors (horizontally)
vcat( [1,2], [3,4])

4-element Vector{Int64}:
 1
 2
 3
 4

# concatenate vectors
hcat( [1,2], [3,4])

2×2 Matrix{Int64}:
 1  3
 2  4

# transpose matrix
hcat( [1,2], [3,4])'

2×2 adjoint(::Matrix{Int64}) with eltype Int64:
 1  2
 3  4

Tuples

# like in python
# tuples can hold heterogenous data
t = ("This", "is", 1, "tuple")

("This", "is", 1, "tuple")

# access elements in a tuple (still 1-based)
t[3]

# tuples are `immutable`
# The following should raise an exception
push!(t, "not a vector")

MethodError: MethodError: no method matching push!(::Tuple{String, String, Int64, String}, ::String)

Closest candidates are:
  push!(::Any, ::Any, !Matched::Any)
   @ Base abstractarray.jl:3389
  push!(::Any, ::Any, !Matched::Any, !Matched::Any...)
   @ Base abstractarray.jl:3390
  push!(!Matched::Set, ::Any)
   @ Base set.jl:103
  ...

Loops

# loop over any iterable (like in python)
for i in 1:5
    println("Iteration ", i)
end
# note how 1 and 5 are both included.

Iteration 1
Iteration 2
Iteration 3
Iteration 4
Iteration 5

for i ∈ ["Paris", "New-York", "Bogota"]
    println("City: ", i)
end

City: Paris
City: New-York
City: Bogota

# Vector comprehensions (like in python)
# enumerate all squares of even numbers (% computes modulo)
[i^2 for i=1:10 if i%2==1]

5-element Vector{Int64}:
  1
  9
 25
 49
 81

Functions

# matlab like syntax
# with positional and keyword arguments separated by `;`

function fun(a,b; c=3)
    z = a+b*c
    return z
end

fun (generic function with 1 method)

fun(1,2)

fun(1,2; c=4)

Manipulating AR1 Models

We consider here a simple autoregressive model:

\[y_t = A y_{t-1} + B e_t\]

where \(y_t=(y^1_t, y^2_t)\) a vector of variables and \(e_t=(e^1_t, e^2_t)\) a normal i.i.d. multivariate process defined by covariance matrix \(\Sigma \in R^p \times R^p\).

We start by choosing:

\[A = \begin{bmatrix}\rho & 0 \\ 2 & \lambda\end{bmatrix}\]

\[B = \begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}\]

\[\Sigma = \begin{bmatrix}0.05 & 0.005 \\ 0.005 & 0.01\end{bmatrix}\]

Define julia variables for matrices \(A, B, \Sigma\).

Compute (programmatically) the eigenvalues of A.

Simulate the response to a one deviation shock to \(e^1_t\) and \(e^2_t\) over 10 periods.

Plot the result using Plots.jl

What do you get if one eigenvalue is equal to 1? Greater than 1?

Import the Distributions package. Use MvNormal to compute draws from a multivariate distribution with covariance matrix \(\Sigma\)

# we need the distributions package to compute mvnormal
import Pkg; Pkg.add("Distributions")

   Resolving package versions...
  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`

using Distributions

# MvNormal
# check the doc

Perform a stochastic simulation over 20 periods starting with a one standard-deviation in \(e_2\).

function simulate(A,B,Σ,e0=zeros(2); T=20)
    
    # initialize distribution object
    dis = MvNormal(Σ)
    
    # create vector to hold the simulation
    sim = [B*e0]

    for t=1:T

        e1 = rand(dis)

        y0 = sim[end] # last value
        y1 = A*y0 + B*e1

        push!(sim, y1)


    end

    return hcat(sim...)
end

simulate (generic function with 2 methods)

Perform K=500 stochastic simulations over 20 periods starting with a one standard-deviation in \(e_2\). Plot the result.

Ergodic Distribution.

It can be shown that the ergodic distribution of a VAR1 is a multivariate normal law, with covariance matrix \(\Omega\).

This matrix is a solution to the equation \(\Omega = A \Omega A' + B \Sigma B'\).

A simple algorithm to find it consist in applying the recurrence \(\Omega_n = A \Omega_{n-1} A' + B \Sigma B'\) until convergence, starting with \(\Omega_0 =0\).

Implement this algorithm.

function ergodic_steady_state(A,B,Σ; N=1000, tol_η=1e-8)
    Ω0 = Σ*0
    for n = 1:N
        Ω = A*Ω0*A'+B*Σ*B'
        η = maximum( abs.( Ω - Ω0) )
        if η<tol_η
            return Ω
        end
        Ω0 = Ω
    end
    error("No convergence")
end

ergodic_steady_state (generic function with 1 method)

Compare the result with the empirical ergodic distribution obtained from the simulations