This post, and other posts with similar tags and headers, are mainly directed at the students who are following the Dynamic Programming course at Dept. of Economics, University of Copenhagen in the Spring 2017. The course is taught using Matlab, but I will provide a few pointers as to how you can use Julia to solve the same problems. So if you are an outsider reading this I welcome you, but you won’t find all the explanations and theory here. Then you’ll have to come visit us at UCPH and enroll in the class!
We start out with very simple examples and Julia code, and gradually ramp up the difficulty level.
Simple model
Consider the simplest consumption-savings model:
\(V^*(M) = \max_{C_1,C_2,\ldots,C_T}\beta^0\sqrt{C_1},\beta^1\sqrt{C2},\ldots,\beta^T\sqrt{C_T}\)
subject to
$$\bar{M} = C_1+C_2+\ldots+C_T\equiv m_1$$
$$m_{t+1} = m_t – C_t$$
$$C_t\in\mathbb{N}$$
$$\beta=0.9,\quad\bar{M}=5$$
Let us try to solve such models.
1) Brute force
This is not a very clever way, but let us start somewhere. We don’t have a better method – yet! Set T=2. Note, there is only one free choice to vary. Given some consumption in period t, the other is given due to the “No cake left” boundary condition.
# Primitives β = 0.90 α = 0.5 u(x, α) = x^α ## State M̄ = 5.0 Cᵒᵖᵗ = u(0.0, α) + β*u(M̄, α) Vᵒᵖᵗ = -Inf for C1 in 1.0:1.0:M̄ V = u(C1, α) + β*u(M̄-C1, α) if Vᵒᵖᵗ <= V Vᵒᵖᵗ = V Cᵒᵖᵗ = C1 end end println("Results") println("* optimimum: ", Vᵒᵖᵗ) println("* optimizer: ", Cᵒᵖᵗ)
This tells us that our agent should consume three units in the first period and two in the last period.
2) Vary β and M
Check that different M’s and β’s give sensible predictions.
We don’t want to rerun that code manually each time – we’re going to need something smarter. Let’s wrap it in a function. In Julia, a typical function (method) definition looks something like the following.
function brute_force(M̄, α, β) Vᵒᵖᵗ = u(0.0, α) + β*u(M̄, α) Cᵒᵖᵗ = 0 for C1 in 1.0:1.0:M̄ V = u(C1, α) + β*u(M̄-C1, α) if Vᵒᵖᵗ <= V Vᵒᵖᵗ = V Cᵒᵖᵗ = C1 end end Vᵒᵖᵗ, Cᵒᵖᵗ end
Notice that whatever is on the last line will be returned. For example, you can verify that it returns the same thing as the loop above, or check what happens if β is 0:
brute_force(M̄, α, 0.0)
The output tells us that the agent should consume everything in the first period. This is consistent with the intuition that the agent doesn’t care about the future at all. What if we have no units at the first period?
brute_force(0, α, β)
We see that the agent should consume nothing. This might seem like a stupid thing to test for, but it is always important to sanity check your code. If it can’t predict things that are obvious, it will probably also give wrong answers to more complex cases. It can be hard to come up with such cases, but the last two examples will always hold. Full discounting -> full consumption today. No ressources -> no consumption today (unless you can borrow).
3) Solve a model using backwards induction
It requires little convincing to realize, that it is wasteful to check all possible combinations – especially when the number of periods increases. Instead, we use the principle of optimality to solve the problem using dynamic programming. To solve the model using dynamic programming, we need to create a few variables to hold the value and policy functions. These will hold the solutions from each (t,m) subproblem. We start with T=3
T = 3 M = 0.0:M̄ Vᵒᵖᵗ = [zeros(M) for t = 1:T] Cᵒᵖᵗ = [zeros(M) for t = 1:T]
We’ve used vectors of vectors here, but could just as well have created a $5\times 3$ matrix in this case. However, in a more general class of models, the number of states might vary over time, and there a vector of vector becomes convenient. With these variables in place, we can solve the last period by realizing that given some resource level in period T, we need to consume everything in order fulfill the constraints of the problem (no cake left on the table!).
Vᵒᵖᵗ[end] = sqrt.(M) Cᵒᵖᵗ[end] = M
Notice the dot after sqrt. This is Julia syntax for “broadcast the square root function over M”. Then, we solve for the period t=T-1, then t=T-2, and then we’re done. We implement a straight forward loop over all possible m’s and feasible consumption given these m’s, by using the fact that the feasible values for consumption make up the set
$$\{0,1,2,…,m_t\}$$.
for t = T-1:-1:1 for imₜ = 1:length(M) Vᵒᵖᵗ[t][imₜ] = -Inf # some value lower than sqrt(0) for (icₜ,cₜ) = enumerate(0:M[imₜ]) v = u(cₜ, α) + β*Vᵒᵖᵗ[t+1][1+imₜ-icₜ] if v > Vᵒᵖᵗ[t][imₜ] Vᵒᵖᵗ[t][imₜ] = v Cᵒᵖᵗ[t][imₜ] = cₜ end end end end println("Results") println("* optimimum: ", Vᵒᵖᵗ) println("* optimizer: ", Cᵒᵖᵗ)
Again, you should try to see if some of all this shouldn’t really be separate functions.
Julia bits and pieces
To close off this post, let me just briefly touch upon some of the things we used above, that might be different from Matlab. To allocate an empty array, write
V = Vector(10)
or a vector with some specific element value (here, a Float64 NaN)
W = fill(NaN, 10)
Indexing into arrays is done with square brackets. For example, let’s assign the first element of W to a variable WW
WW = W[1]
And the same goes for setting the value at some index
W[1] = 0.4; W
Functions are define using either short or long form
function f(x,y) x+y+rand() # returns the sum of the two inputs and a random number end g(x,y) = x+y+rand()
and comments are inserted using # or the #= =# blocks
# This is a comment #= This is a comment =#
For-loops are written as
for a = some_iterator # do something end
or to get both the value and a counter use enumerate
for (i, a) = enumerate(some_iterator) # do stuff with a and index something with i end
Code that should be run conditional on some statement being true is written as
if condition # do something if condition is satisfied end
and that’s all for now, folks!