Julia Basics & VAR Process

Business Cycles and Fluctuations - AE2E6

  • 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)
(1.0+(2.0+3.0*(4.0+5.0)))/30
1.0
# exponentials are denoted by ^
2^8
256

Variables / assignments / comparisons

# variable assignment
x = 10
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]
1
# access elements in a matrix
M[1,2]
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]
1
# 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
  ...

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

  ...





Stacktrace:

 [1] top-level scope

   @ ~/Teaching/ensae/ae2e6/tutorials/session_1/index.ipynb:3

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)
7
fun(1,2; c=4)
9

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 & 2 \\ 0 & \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\).

# your code

Compute (programmatically) the eigenvalues of A.

# your code

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

# your code

Plot the result using Plots.jl

# your code

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

# your code

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
# run the following in case it is not already installed
# 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
# your code

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)
# your code

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

# your code

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

# your code