Value Function Iteration with Grid Search
In this notebook, I will solve a simple Dynamic Programming Problem: The Neoclassical Growth Model without uncertainty. I will solve the DP problem by discretizing the Bellman equation and performing value function iteration. The maximization problem is carried on using grid search (brute force).
This note aims to be pedagogical, and does not try to achieve maximum performance.
Neoclassical growth model
Consider the Neoclassical growth model: the production function is
The representative household chooses a consumption sequence to maximize the following utility function:
It is easy to show that the steady-state capital stock is given by:
Bellman Equation
The Bellman Equation associated to this problem is:
with the policy function:
Numerical Implementation
First, define the parameters of the model:
xxxxxxxxxx
md"""
## Value Function Iteration with Grid Search
In this notebook, I will solve a simple Dynamic Programming Problem: The Neoclassical Growth Model without uncertainty. I will solve the DP problem by discretizing the Bellman equation and performing value function iteration. The maximization problem is carried on using grid search (brute force).
This note aims to be pedagogical, and does not try to achieve maximum performance.
##### Neoclassical growth model
Consider the Neoclassical growth model: the production function is $k^\alpha$ and capital depreciates at rate $\delta$. The law of motion for the capital stock is:
$k_{t+1} = (1-\delta) k_t + k^\alpha_t - c_t.$
The representative household chooses a consumption sequence to maximize the following utility function:
$\sum_{t=0}^\infty \beta^t \log(c_t).$
It is easy to show that the steady-state capital stock is given by:
$k_{ss} = \left(\dfrac{\alpha}{1/\beta - (1-\delta)} \right)^{\frac{1}{1-\alpha}}.$
##### Bellman Equation
The Bellman Equation associated to this problem is:
$V(k) = \max_{k'} \{\log(k^\alpha + (1-\delta) - k') + \beta V(k')\},$
with the policy function: $k' = g(k)$.
##### Numerical Implementation
First, define the parameters of the model: $\alpha = 0.3$, $\beta = 0.96$, and $\delta = 0.1$. Second, discretize the capital space in $n_k = 200$ points. The points are defined in a grid (a vector) linearly spaced between $2\times k_{ss}/n_k$ and $2\times k_{ss}$.
"""
0.0292082
0.0584164
0.0876247
0.116833
0.146041
0.175249
0.204458
0.233666
0.262874
0.292082
0.32129
0.350499
0.379707
0.408915
0.438123
0.467332
0.49654
0.525748
0.554956
0.584164
5.57877
5.60798
5.63719
5.66639
5.6956
5.72481
5.75402
5.78323
5.81244
5.84164
xxxxxxxxxx
begin
alpha = 0.3
beta = 0.96
delta = 0.1
k_ss = (alpha/((1/beta)-(1-delta)))^(1/(1-alpha)) # SS capital
nK = 200 # grid points
gK = collect(range(2*k_ss/nK, 2*k_ss, length = nK)) # array (see output above!)
end
Value Function Iteration
I will now perform the iteration of the Value Function. Define the Value Function in the grid
The general procedure involve to guess a value function,
where
Afterwise, compute the (maximum) distance between the
If this distance is smaller than some tolerance,
In the code below, I will initialize the Value function (
It will be useful to pre-compute the consumption array (matrix
I will also define the tolerance and the maximum number of iterations. See the consumption matrix below:
xxxxxxxxxx
md"""
##### Value Function Iteration
I will now perform the iteration of the Value Function. Define the Value Function in the grid $n_k$: $V(k_i) = V_i$. That is, $V(k_i)$ is the associated value at capital $k_i$.
The general procedure involve to guess a value function, $V^n(k_i)$, solve the maximization problem and compute the next iteration, $V^{n+1}(k_i)$:
$V^{n+1}(k_i) = \max_{k_j} \{\log(c_{ij}) + \beta V^{n}(k_j) \},$
where $i$ is the subscript associated with the state $k_i$, while $j$ is the grid point associated with control $k'_j$. $c_{ij}$ is the consumption associated with state $i$ and control $j$.
Afterwise, compute the (maximum) distance between the $V^n$ and $V^{n+1}$:
$d = \max_{i = 1,..,n_k} | V^{n+1}_i - V^n_i|.$
If this distance is smaller than some tolerance, $d < \varepsilon$, stop. Otherwise, use $V^{n+1}(k_i)$ as the new guess and repeat the procedure until convergence.
In the code below, I will initialize the Value function ($V^0_i = 0.0$, but it could be different) and the Policy function in the grid $n_k$. Note that the array $g$ will tell us that if we are in state $i$, we will choose the point $j$ in the capital grid, so $ g_i = j$ (so $j$ is the grid point NOT the capital value).
It will be useful to pre-compute the consumption array (matrix $n_k\times n_k$) with elements $c_{ij}$ defined as:
$c_{ij} = k_i^\alpha + (1-\delta)k_i - k_j.$
I will also define the tolerance and the maximum number of iterations. See the consumption matrix below:
"""
200×200 Matrix{Float64}:
0.343538 0.31433 0.285121 0.255913 0.226705 … -5.41048 -5.43969 -5.4689
0.449907 0.420699 0.391491 0.362283 0.333074 -5.30411 -5.33332 -5.36253
0.531366 0.502158 0.47295 0.443742 0.414534 -5.22265 -5.25186 -5.28107
0.601075 0.571866 0.542658 0.51345 0.484242 -5.15294 -5.18215 -5.21136
0.663719 0.634511 0.605303 0.576095 0.546887 -5.0903 -5.11951 -5.14872
0.721574 0.692366 0.663157 0.633949 0.604741 … -5.03245 -5.06165 -5.09086
0.775931 0.746723 0.717515 0.688307 0.659098 -4.97809 -5.0073 -5.0365
⋮ ⋱
6.78208 6.75287 6.72366 6.69445 6.66524 1.02806 0.998847 0.969639
6.81095 6.78174 6.75253 6.72333 6.69412 … 1.05693 1.02772 0.998515
6.83982 6.81061 6.7814 6.75219 6.72298 1.0858 1.05659 1.02738
6.86867 6.83947 6.81026 6.78105 6.75184 1.11465 1.08545 1.05624
6.89752 6.86831 6.83911 6.8099 6.78069 1.1435 1.11429 1.08509
6.92636 6.89715 6.86794 6.83874 6.80953 1.17234 1.14313 1.11392
xxxxxxxxxx
begin
V = fill(0.0, nK) # note this is a float array
gp = fill(0, nK) # note this is an integer array
maxiter = 1000
tol = 10^(-6)
c = @. gK^alpha + (1-delta)*gK - gK' # note the transpose. rows: state; cols: cont
end
Note that some
In the code below, I will implement the value function iteration. The outer loop iterate over
Then, in the most inner loop, I choose the
Note: The loop can be very inefficient in some languages. Nevertheless, it is pretty easy to incorporate extensions where we exploit the monotonicity/concavity of the VF.
xxxxxxxxxx
md"""
Note that some $c_{ij}$ are negative. Be careful to rule out this possibility during the maximization.
In the code below, I will implement the value function iteration. The outer loop iterate over $n$. Then, in the first inner loop, I start to solve for every value function $V_i$. In short, I iterate over the state $i$.
Then, in the most inner loop, I choose the $j$ that maximizes the VF. There are many ways you can do this step. Here, I am just looping over all $j$ values (but skipping the ones that have negative consumption), sometimes people refer to this approach as grid search. The procedure is summarized as following: first, initialize $V_{max}$ to large negative number; second, for every $j$, compute the VF associated with the control $j$, and in case $V_j > V_{max}$, set $V_{max} = V_j$.
_Note_: The loop can be very inefficient in some languages. Nevertheless, it is pretty easy to incorporate extensions where we exploit the monotonicity/concavity of the VF.
"""
xxxxxxxxxx
begin
for iter = 1:maxiter
V_next = fill(0.0, nK) # here we store the Vn+1
# solve value function
for i = 1:nK # this iterates for every state i
# perform maximization problem
V_g = -1.0*10e100 # These are temporary values for the value and policy
g_g = 0
for j = 1:nK
if c[i, j] <= 0; # if consumption is negative we skip this iteration
continue
end
vj = log(c[i, j]) + beta * V[j] # compute the VF
if vj > V_g # perform maximization
V_g = vj
g_g = j
end
end
V_next[i] = V_g # once we finish the maximization we assign to Vn+1
gp[i] = g_g
end
# Check point here: time limit / iteration / sup
d = maximum(abs.(V - V_next))
if d < tol
println("Tol. achieved: $d")
break # break the loop in case we found the VF!
end
if iter == maxiter;
println("Max iterations achieved. VF did not converge")
end
println("Iter: $iter")
println("Tol: $d") # ps these prints will not show in the notebook.
V .= V_next # Update the guess
end
end
Let us now see the array V and g:
xxxxxxxxxx
md"""Let us now see the array V and g:"""
-4.30336
-3.78435
-3.4481
-3.18961
-2.97673
-2.79305
-2.62913
-2.48173
-2.34686
-2.22205
-2.10546
-1.9957
-1.89155
-1.79267
-1.69834
-1.60818
-1.52171
-1.43841
-1.35806
-1.28039
4.1403
4.15891
4.1775
4.19604
4.21448
4.23289
4.25125
4.2695
4.28769
4.30586
xxxxxxxxxx
V
6
9
10
12
13
15
16
17
18
19
20
22
23
24
25
26
27
28
29
30
176
177
178
179
180
180
181
182
183
184
xxxxxxxxxx
gp
It is useful to construct the policy function with capital values instead of grid points.
xxxxxxxxxx
md"""It is useful to construct the policy function with capital values instead of grid points."""
0.175249
0.262874
0.292082
0.350499
0.379707
0.438123
0.467332
0.49654
0.525748
0.554956
0.584164
0.642581
0.671789
0.700997
0.730206
0.759414
0.788622
0.81783
0.847038
0.876247
5.14065
5.16986
5.19906
5.22827
5.25748
5.25748
5.28669
5.3159
5.3451
5.37431
xxxxxxxxxx
g = gK[gp]
xxxxxxxxxx
begin
using Plots # load plot package
plot(gK, V)
xlabel!("k")
ylabel!("V(k)")
end
xxxxxxxxxx
begin
plot(gK, g)
xlabel!("k")
ylabel!("g(k)")
end