Introduction to Julia on Computational Economics

Xing Xu, University of Minnesota, 2024 Summer

Overview

The purpose of the class is to get you familiarize with some basic dynamic programming (DP) problems in economics and train you to be comfortable with coding them in Julia.

Why DP?

Dynamic programming is the key to solving modern economic models. Most macro models now are dynamic general equilibrium models. Knowing how to solve recursive models is in additional useful for searching, matching and even game theory (with sequential games).

### Why Julia?

Julia is the most modern tool for computations. It preserves the efficiency and clarity of coding while providing astounding speed. For comparison, c++ and Fortran are as fast but coding in them requires a lot more work as basically everything has to be built from scratch. Another extreme is Stata, which is super easy to code in (say, reg y x), but is absurdly slow and extremely limited in its applications (can’t even perform DP).

Julia reaches a great balance and is perfectly suited for computational economics. You will see its beauty in action once you get comfortable with it.

Background knowledge

I will list the backgrounds used that are considered as known. Additional knowledge will be specified when used.

Math Preliminaries (undergraduate level knowledge is sufficient):

Linear algebra
Real analysis
Probability theory
Basic understanding of Markov Chains
Knowledge of basic optimization theory will be a plus (Sundaram <A First Course in Optimization Theory>)

Economics: * Basic choice theory (Chapter 2 MWG) (Preferences, utility representation)

Basic consumer theory (Chapter 3 MWG) (Utility maximization)
Basic producer theory (Chapter 5 MWG)
An idea of Competitive Equilibrium and Welfare Theorems (Chapter 10 MWG)

Julia (required reading): * Setup of Julia environment, Quantecon 1

Introduction to Julia Quantecon 2
Julia Essentials Quantecon 3
Arrays, Tuples, Ranges, and Other Fundamental Types Quantecon 4
Basic understanding is good: Introduction to Types and Generic Programming Quantecon 5

Other good references:

Notes: * Dirk Krueger’s Macroeconomic Theory notes

Tim Kehoe’s notes on Blackwell sufficient conditions

Books: * Adda & Cooper <Dynamic economics> (strongly recommend)

Stokey, Lucas with Prescott <Recursive methods in economic dynamics>
Ljungqvist & Sargent <Recursive Macroeconomic Theory>
Heer & Maußner <Dynamic General Equilibrium Modeling>
Stachurski <Economic Dynamics: Theory and Computation (Second Edition)>

Lecture 1: Julia Fundamentals

This acts as a review for Lecture 1.1 - 1.5 for Julia on Quantecon. I will only focus on stuffs we will use for the following lectures. This lecture accompanies a succinct assignment to get you familiarized with writing in the language and playing with it on the Jupyter Notebook.

1.1 Packages and Projects

Working with packages is of vital importance in Julia. We use “LinearAlgebra” for matrix operations, “Statistics” for different distributions, “Plots” for, obviously, plots, etc. One thing you always need to keep in mind is how to download some packages.

After downloading the essential packages you need, you have developed something called an “environment”. For compatibility, version control and reproducibility reasons, one often wants to create individual environments for different files. “Projects” manages each individual local environments. It will create a project repository that usually contains a “Project.toml”, a “Manifest.toml” and your code files.

1.1.1 Packages and Projects using the Pkg API

Now what are these exactly? Put this notebook into an individual folder and run the following lines.

import Pkg

Pkg.activate(".")  # Activate a local environment 
Pkg.add("LinearAlgebra") # add the package "LinearAlgebra"

  Activating project at `~/Documents/Teaching/Julia course 2024/Julia/Lectures`
   Resolving package versions...
  No Changes to `~/Documents/Teaching/Julia course 2024/Julia/Lectures/Project.toml`
  No Changes to `~/Documents/Teaching/Julia course 2024/Julia/Lectures/Manifest.toml`

You should find two extra files in the folder. A “Project.toml” and a “Manifest.toml”. Double click on each of them in the folder or just click on them in the editor in Vscode. Take a look. What are in them?

The “Project.toml” file contains the name of the package, its unique UUID, its version, the authors and potential dependencies.

The “Manifest.toml” file shows the information for the resulting dependencies of the added packages.

You can check the status of your package using the following code.

Pkg.status()

Status `~/Documents/Teaching/Julia course 2024/Julia/Lectures/Project.toml`
  [37e2e46d] LinearAlgebra

Sometimes you want to run other people’s code. Say you have stored some people’s project into a folder on your laptop (you just need the “Project.toml” file). To activate the exact environment, the following code will automatically download all packages and their dependencies of the same version and make it the current environment.

Pkg.instantiate()

1.1.2 Packages and Projects using REPL

All the above operations can also be done using REPL (Read-Evaluate-Print-Loop).

To activate REPL in Vscode, press CTRL + Shift + P in Windows or CMD + Shift + P in OS, then type “start REPL” in the command palette and click on it.

Now the Terminal should pop up and you should see “julia>”. This indicates we have successfully started the REPL.

Next type ] to enter the package mode. It should write something like “(Julia course 2023 Summer) pkg>”. To exit, simply press delete.

With the package mode, we can perform the same functionality as before and more.

activate. activate the project at the current directory.

st shows the status.

add LinearAlgebra add the LinearAlgebra package to the local environment.

instantiate automatically downloads all the missing dependencies from the particular “Project.toml”.

1.2 Introductory example for using packages

We will see a lot more packages and their usages along the way. For the moment, I will show some basic functionalities of the LinearAlgebra package as an example. For more information, check here.

First, you need to call the Package with “using”.

using LinearAlgebra

Define a matrix:

A = [1 2 3; 4 5 6; 7 8 9]

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

A couple things you can do. The codes explain themselves.

tr(A)

det(A)

-9.516197353929915e-16

inv(A)

3×3 Matrix{Float64}:
  3.15252e15  -6.30504e15   3.15252e15
 -6.30504e15   1.26101e16  -6.30504e15
  3.15252e15  -6.30504e15   3.15252e15

norm(A)

16.881943016134134

eigvals(A)

3-element Vector{Float64}:
 -1.1168439698070434
 -8.582743335036247e-16
 16.11684396980703

eigvecs(A)

3×3 Matrix{Float64}:
 -0.78583     0.408248  -0.231971
 -0.0867513  -0.816497  -0.525322
  0.612328    0.408248  -0.818673

I represents an identity matrix of any size. For example,

A * I

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

The dot product can also be calculated.

dot(A, I) # Same as tr(A)

A ⋅ I
#⋅ is typed by "\cdot" then press "tab"

Asolves Ax = B. Note that this is strictly preferred to inv(A)*B. Never do the latter.

B = [1, 2, 3]
A\B
# or B/A

3-element Vector{Float64}:
 -0.23333333333333334
  0.46666666666666673
  0.09999999999999994

Element-wise operations are done through . (more discussions on this next lecture).

A .+ 3

3×3 Matrix{Int64}:
  4   5   6
  7   8   9
 10  11  12

You can also perform factorizations on the matrices.

lu(A)

LU{Float64, Matrix{Float64}}
L factor:
3×3 Matrix{Float64}:
 1.0       0.0  0.0
 0.142857  1.0  0.0
 0.571429  0.5  1.0
U factor:
3×3 Matrix{Float64}:
 7.0  8.0        9.0
 0.0  0.857143   1.71429
 0.0  0.0       -1.58603e-16

qr(A)

LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}}
Q factor:
3×3 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}:
 -0.123091   0.904534   0.408248
 -0.492366   0.301511  -0.816497
 -0.86164   -0.301511   0.408248
R factor:
3×3 Matrix{Float64}:
 -8.12404  -9.60114   -11.0782
  0.0       0.904534    1.80907
  0.0       0.0        -8.88178e-16

# Cholesky decomposition doesn't work here, which creates an error
cholesky(A)

PosDefException: PosDefException: matrix is not Hermitian; Cholesky factorization failed.

Not surprisingly, we can also do vector operations.

v1 = [1, 1, 3]
v2 = [2, 2, 3]

3-element Vector{Int64}:
 2
 2
 3

v1 + v2 
# same as v1 .+ v2

3-element Vector{Int64}:
 3
 3
 6

Dot products for vectors are same as for matrices, except for vectors we can simply do:

v1' * v2 #' for transpose

dot(v1, v2)

v1 ⋅ v2

α = 0.3
α * 2

0.6

1.2 Functions

In Julia functions can be built in line or specifically calling function.

f(x, y) =  x^y
f(2,2)

#Or
function f(x, y)
    return x^y
end
f(2,2)

If you want to input an array into f, it gives an error as it is not clear what you mean.

v1 = [1, 2, 3]
f(v1, 2)

MethodError: MethodError: no method matching ^(::Vector{Int64}, ::Int64)

Closest candidates are:
  ^(!Matched::Union{AbstractChar, AbstractString}, ::Integer)
   @ Base strings/basic.jl:733
  ^(!Matched::Diagonal, ::Integer)
   @ LinearAlgebra ~/Applications/Julia-1.9.3.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/diagonal.jl:208
  ^(!Matched::Diagonal, ::Real)
   @ LinearAlgebra ~/Applications/Julia-1.9.3.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/diagonal.jl:207
  ...

You can modify this by defining the function specifically for vectors.

function g(x,y)
    z = similar(x) # similar() generates an size/matrix of the same size and type as the input
    for i in eachindex(x)
        z[i] = x[i]^y
    end
    return z
end
@show g(v1, 2)

g(v1, 2) = [1, 4, 9]

3-element Vector{Int64}:
 1
 4
 9

Or in a more concise fashion,

g(x,y) = [z^y for z in x] # or g(x,y) = [x[i]^y for i in eachindex(x)]
g(v1, 2)

3-element Vector{Int64}:
 1
 4
 9

Instead, you can also apply “.” to broadcast the function onto every element of v1. This comes very handy in a lot of cases.

f(x, y) =  x^y
f.(v1, 2)
# should write f.(v1, Ref(2)) for rigidity. You can think of Ref() as fixing the input here.

3-element Vector{Int64}:
 1
 4
 9

1.3 Types (a short summary)

Types are incredibly important for coding in Julia. It matters for the correctness, consistency and efficiency of the codes. However, my past experience told me that it is almost impossible to learn this from pure memorization. You have to actually work with them and think along the way. Here I will just let you see some commonly used types and give a simple example about declaration of types.

For more info on types, check the documentation.

The default type for (simple) integers is Int64

@show typeof(1)

typeof(1) = Int64

Int64

The default type for (simple) floats is Float64

@show typeof(1.0) #putting ".0" after 1 make the Julia compiler infer it is a float

typeof(1.0) = Float64

Float64

The default type for strings is String

@show typeof("string")

typeof("string") = String

String

Type of Boolean (true of false) variables.

@show typeof(true)
@show typeof(false)

typeof(true) = Bool
typeof(false) = Bool

Bool

# Note that sometimes we use the feature that true equals to 1 and false equals to 0
@show true == 1
@show false == 0

true == 1 = true
false == 0 = true

true

What are the type of arrays?

@show typeof([1, 2, 3])

typeof([1, 2, 3]) = Vector{Int64}

Vector{Int64} (alias for Array{Int64, 1})

@show typeof([1.0, 2.0, 3.0])

typeof([1.0, 2.0, 3.0]) = Vector{Float64}

Vector{Float64} (alias for Array{Float64, 1})

How about the type of [1, 2.0, 3.0], where the first entry is an integer?

@show typeof([1, 2.0, 3.0])

typeof([1, 2.0, 3.0]) = Vector{Float64}

Vector{Float64} (alias for Array{Float64, 1})

Note that it returns a vector of floats. This is subtle but keep in mind that sometimes these type of inconsistencies can create problems. This also indicates that the type within an array is unified, to find the type of the elements in an array. We commonly use the following eltype().

@show eltype([1.0, 2.0, 3.0])
# differentiate with typeof([1.0, 2.0, 3.0])

eltype([1.0, 2.0, 3.0]) = Float64

Float64

Tuples are constructed using small brackets. We can also check the type of a tuple, which is different from the behavior of an array (vector).

@show typeof((1, 2.0, 3.0, "string"))

typeof((1, 2.0, 3.0, "string")) = Tuple{Int64, Float64, Float64, String}

Tuple{Int64, Float64, Float64, String}

Now let’s look at the function in 1.2 again, which is

function g(x,y)
    z = similar(x)
    for i in eachindex(x)
        z[i] = x[i]^y
    end
    return z
end
@show g(v1, 2)

g(v1, 2) = [1, 4, 9]

3-element Vector{Int64}:
 1
 4
 9

In practice, sometimes we want to specify the types of inputs, which can be down using ::.

function g(x :: Vector{Int64},y :: Int64)
    z = similar(x)
    for i in eachindex(x)
        z[i] = x[i]^y
    end
    return z
end
@show g(v1, 2)

g(v1, 2) = [1, 4, 9]

3-element Vector{Int64}:
 1
 4
 9

g([2], 2)

1-element Vector{Int64}:
 4

g([2.0], 2)

1-element Vector{Float64}:
 4.0

g(2.0, 2) # this doesn't work

MethodError: MethodError: no method matching similar(::Float64)

Closest candidates are:
  similar(!Matched::Union{Adjoint{T, var"#s972"}, Transpose{T, var"#s972"}} where {T, var"#s972"<:(AbstractVector)})
   @ LinearAlgebra ~/Applications/Julia-1.9.3.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/adjtrans.jl:329
  similar(!Matched::Union{Adjoint{T, var"#s972"}, Transpose{T, var"#s972"}} where {T, var"#s972"<:(AbstractVector)}, !Matched::Type{T}) where T
   @ LinearAlgebra ~/Applications/Julia-1.9.3.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/adjtrans.jl:330
  similar(!Matched::Union{Adjoint{T, S}, Transpose{T, S}} where {T, S})
   @ LinearAlgebra ~/Applications/Julia-1.9.3.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/adjtrans.jl:333
  ...

However this sometimes creates confusions. We can input a matrix and still get out results:

m1 = [2 2; 2 2]
g(m1, 2)

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

This is why Quantecon stresses the importance for generic programming. See Quantecon 5 and Quantecon 6 for details. In short, you should avoid specifying particular types unless necessary and let your codes be compatible with different types.

In addition, keep in mind that you might specify some types without even knowing that you are doing so. One example is when you want to create a vector of 0 of type Int with length 5. It is tempted to use the built-in function zeros(5). Beware that this implicitly imposes that 0s are of type float.

zeros(5)

5-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0

@show eltype(zeros(5))

eltype(zeros(5)) = Float64

Float64

Suppose you want the zeros with type int, one way is to use Int(). Again you need to broadcast it to each element.

Int.(zeros(5))

5-element Vector{Int64}:
 0
 0
 0
 0
 0

Or you can simply declare within zeros().

zeros(Int,5)

5-element Vector{Int64}:
 0
 0
 0
 0
 0

1.4 Help?

How do you know you can do something as zeros(Int,5)? For most built-in functions, you can find it using Help from the REPL.

If you haven’t quitted Vscode, you are probably still in REPL from 1.1.

If you have, press CTRL + Shift + P in Windows or CMD + Shift + P in OS, then type “start REPL” in the command palette and click on it. You should again see “julia>”, which indicates you are at REPL.

Now simply typing ? enters into the help mode. You should see “help?>”. Let’s try with the function zeros(). Type in the zeros or zeros() and press Return on your keyboard. You should see it shows the first input of function is “[T=Float64,]”, so unless otherwise specified, zeros() will generate an array with type “Float64”.

Help? is usually the first place to go if you encounter some built-in functions that you are not sure with. As we go, you will be introduced to a lot more of them.

1.5 Macros (not Macroeconomics)

Macros are a powerful tool in Julia. Macros are a type of metaprogramming (whatever this means). Just keep in mind that Macros start with an @ and take in general expressions instead of values as inputs. You can also check with help? for documentation of each Macro.

I know it sounds a little absurd but knowing Macros well is a ticket to becoming a Julia master. For the moment, knowing a couple examples is sufficient.

We have used @show to show the output of a function. However, note that the inputs it can take is much more general. Say, it can take in the function itself.

@show f

f = f

f (generic function with 1 method)

@show f(1,2)

f(1, 2) = 1

@time begin ... end works like tic ... tok in Matlab. It gives the duration of running the code between.

@time begin
    1 + 1
end

  0.000001 seconds

Another example is @assert, which can be used for generating some statements if some conditions are not met.

v1 = [1, 2, 3]
v2 = [4, 5, 6]
@assert v1 == v2 "v1 is not equal to v2"

AssertionError: AssertionError: v1 is not equal to v2

# no statement is generated if the condition is met
@assert v1 == v1 "v1 is not equal to v1"

A couple other useful Macros:

@. Broadcasts . onto all arguments in one line.

@inbounds Eliminates array bounds checking within expressions (efficiency).

@inline Performs inline Maths (efficiency).

@view Creates a data structure that acts like an array (it is a subtype of AbstractArray ), but the underlying data is actually part of another array. Often used to create a temporary array for efficiently.

For more, check here.

j(x) = x^2
v1 = [1,2,3]
j.(v1)

3-element Vector{Int64}:
 1
 4
 9

@. (j(v1) + 2)*2

3-element Vector{Int64}:
  6
 12
 22

1.6 Miscellaneous facts about Julia (under construction)

Julia has a particular behavior when assigning an array to another. What is your guess for the result of the following code?

x = [3, 3]
y = x
y[1] = 6
@show x

x = [6, 3]

2-element Vector{Int64}:
 6
 3

One might guess that the output will be “[3, 3]” as x is not changed. It is exactly what will happen in Python or Matlab. However, running it gives a seemingly surprising result.

x = [3, 3]
y = x
y[1] = 6
@show x

x = [6, 3]

2-element Vector{Int64}:
 6
 3

# also
x[2] = 6
@show y

y = [6, 6]

2-element Vector{Int64}:
 6
 6

The reason is that on the background of Python or Matlab. When you run y = x, they create a copy for x in the storage and make it y. However in Julia, y = x simply assigns y to the same array that x is assigned to. That’s why changing either of them changes the other.

There are two alternative ways if you want the same behavior as in Python or Matlab.

x = [3, 3]
y = copy(x)
y[1] = 6
@show x

x = [3, 3]

2-element Vector{Int64}:
 3
 3

x = [3, 3]
y .= x #element-wise assignment dodges the problem
y[1] = 6
@show x

x = [3, 3]

2-element Vector{Int64}:
 3
 3

Either one of them essentially makes y and x to be assigned to “different” arrays in your storage. If you know what a “pointer” is in c++, they now simply point to different places.