Problem Set 2: Problem 4 Coding SolutionsĀ¶

Goal: solve the following Bellman equation using value function iteration:

$V (x) = \max_{s \in [0,x]} \{\log (x-s) + \delta \mathbb{E} V (s + \tilde{y}_{+1})\}$

InĀ [1]:
## Load packages ##
#using Pkg
#Pkg.add("Plots")
using LinearAlgebra
using Plots

Precompiling Plots
ā”Œ Warning: attempting to remove probably stale pidfile
ā”‚   path = "/Users/sarahrobinson/.julia/compiled/v1.10/PlotUtils/YveHG_Kzo3C.ji.pidfile"
ā”” @ FileWatching.Pidfile /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/FileWatching/src/pidfile.jl:244
  āœ“ PlotUtils
  āœ“ PlotThemes
  āœ“ RecipesPipeline
  āœ“ Plots
  āœ“ Plots ā†’ UnitfulExt
  5 dependencies successfully precompiled in 43 seconds. 149 already precompiled.
Precompiling IJuliaExt
  āœ“ Plots ā†’ IJuliaExt
  1 dependency successfully precompiled in 2 seconds. 161 already precompiled.

Part 1Ā¶

InĀ [54]:
## Part 1: Specify parameters values ##
Ī“ = 0.9

# Make vector with \tilde{y} values and probabilities (Py)
Y = [y for y in 0.1:0.1:0.9]
Py = ones(9)./9 # uniform distribution probabilities

# Wealth level vector #
X = [x for x in 0.01:0.01:10];

Part 2Ā¶

Take one iteration of the Bellman equation using the initial guess $V(x_{+1}) = x_{+1} $

InĀ [3]:
## Specify initial guess ##
V0 = X;

Part (a)Ā¶

InĀ [4]:
## Part (a): Solve for continuation value function $V (s + \tilde{y}_{+1})$) for each s value 

#### Subpart (i): Create two matricies for each set of pairs (s, y) \in (X x Y) ####
# Note that this gives us two 1000x9 matrities #
s_vals = repeat(X,1,9)
y_vals = repeat(Y',1000,1);

{Note} a unique feature of how we have defined the grids. For the s_vals matrix, the first row corresponds with s = 0.01, the second row with s = 0.02, the third row with s = 0.03, ...., the 999th row with 9.99, and the 1000th row with 10.00. The presents an easy way to go back and forth between row or column {index} and value of $s$. Specifically, for any value of $s$, this value of $s$ lands on the grid is in row $s*100$. Do you see this? This is an important feature that we will use below for both values of $s$ and $x$.

InĀ [5]:
s_vals
Out[5]:
1000Ɨ9 Matrix{Float64}:
  0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01
  0.02   0.02   0.02   0.02   0.02   0.02   0.02   0.02   0.02
  0.03   0.03   0.03   0.03   0.03   0.03   0.03   0.03   0.03
  0.04   0.04   0.04   0.04   0.04   0.04   0.04   0.04   0.04
  0.05   0.05   0.05   0.05   0.05   0.05   0.05   0.05   0.05
  0.06   0.06   0.06   0.06   0.06   0.06   0.06   0.06   0.06
  0.07   0.07   0.07   0.07   0.07   0.07   0.07   0.07   0.07
  0.08   0.08   0.08   0.08   0.08   0.08   0.08   0.08   0.08
  0.09   0.09   0.09   0.09   0.09   0.09   0.09   0.09   0.09
  0.1    0.1    0.1    0.1    0.1    0.1    0.1    0.1    0.1
  0.11   0.11   0.11   0.11   0.11   0.11   0.11   0.11   0.11
  0.12   0.12   0.12   0.12   0.12   0.12   0.12   0.12   0.12
  0.13   0.13   0.13   0.13   0.13   0.13   0.13   0.13   0.13
  ā‹®                                  ā‹®                   
  9.89   9.89   9.89   9.89   9.89   9.89   9.89   9.89   9.89
  9.9    9.9    9.9    9.9    9.9    9.9    9.9    9.9    9.9
  9.91   9.91   9.91   9.91   9.91   9.91   9.91   9.91   9.91
  9.92   9.92   9.92   9.92   9.92   9.92   9.92   9.92   9.92
  9.93   9.93   9.93   9.93   9.93   9.93   9.93   9.93   9.93
  9.94   9.94   9.94   9.94   9.94   9.94   9.94   9.94   9.94
  9.95   9.95   9.95   9.95   9.95   9.95   9.95   9.95   9.95
  9.96   9.96   9.96   9.96   9.96   9.96   9.96   9.96   9.96
  9.97   9.97   9.97   9.97   9.97   9.97   9.97   9.97   9.97
  9.98   9.98   9.98   9.98   9.98   9.98   9.98   9.98   9.98
  9.99   9.99   9.99   9.99   9.99   9.99   9.99   9.99   9.99
 10.0   10.0   10.0   10.0   10.0   10.0   10.0   10.0   10.0
InĀ [6]:
#### Subpart (ii): For each pair (s, y), get x_{+1} ####

Note from our wealth equation we have that $x_{+1} = s + \tilde{y}_{+1}$. This results in a 1000x9 matrix that gives us every possible value of $x_{+1}$ given a value of $s$ and $\tilde{y}_{+1}$. The matrix corresponds to the following ($s, \tilde{y}_{+1}$) pairs:

$$ \begin{bmatrix} (0.01, 0.1) & (0.01, 0.2) & \dots & (0.01, 0.9) \\ (0.02, 0.1) & (0.02, 0.2) & \dots & (0.02, 0.9) \\ \vdots & \vdots & \ddots & \vdots \\ (10, 0.1) & (10, 0.2) & \dots & (10, 0.9) \end{bmatrix} $$

So each row corresponds to a given savings value $s$, and each column corresponds to a given value of $\tilde{y}_{+1}$.

InĀ [7]:
# Calculate x_{+1} = s + \tilde{y}. Note, however, that we don't allow for x_{+1} > 10 because that is outside our pre-defined grid #
new_vals = min.(s_vals + y_vals,10);
InĀ [8]:
#### In the function below, we will do subpart (iii) and (iv) ####

#### Subpart (iii): Find associated value V(x_{+1}) from initial grid of guesses V(X) i.e. V0 ####
# Intuition: given that our grid is discrete, we need to make sure that each V(x_{+1}) lands on a point on the grid. 
# I.E. We need to map each "new_vals" to a point on the V0 grid. For this part, we will use the relationship between 
#the {index} and the value of {x} described above.

#### Subpart (iv): for each s (i.e. for each row) take the expectation of y to get V^C(s). See below matrix representation. ####
#### To calculate the expectation of \tilde{y}_{+1}, we use the probabilities Py ####
$$ V^C (s) = \begin{bmatrix} (0.01, \mathbb{E} \tilde{y}_{+1}) \\ (0.02, \mathbb{E} \tilde{y}_{+1}) \\ (0.03, \mathbb{E} \tilde{y}_{+1}) \\ \vdots \\ (10,\mathbb{E} \tilde{y}_{+1}) \end{bmatrix} $$
InĀ [9]:
## For a given value function, calculate VC ##
function VC_function(V0)
    V_vals = V0[Int.(round.(new_vals*100))] # subpart (iii)
    VC = V_vals*Py # subpart (iv)
    return VC
end
Out[9]:
VC_function (generic function with 1 method)

Part (b) : Solve for the value function V(x)Ā¶

InĀ [10]:
#### Subpart (i): Create two matricies for every pair of wealth $x$ and saving $s$ ####
# This results in two 1000 x 1000 matrices. #
$$ X = \begin{bmatrix} (0.01) & (0.01) & \dots & (0.01) \\ (0.02) & (0.02) & \dots & (0.02) \\ \vdots & \vdots & \ddots & \vdots \\ (10) & (10) & \dots & (10) \end{bmatrix} $$$$ s = \begin{bmatrix} (0.01) & (0.02) & \dots & (10) \\ (0.01) & (0.02) & \dots & (10) \\ \vdots & \vdots & \ddots & \vdots \\ (0.01) & (0.02) & \dots & (10) \end{bmatrix} $$

Note that each row of $X$ represents a given value of X (importantly remember that the row index divided by 100 gives the value of $X$ i.e. the first row is 0.01, second row is 0.02, etc.). Each column of $s$ represents a given value of $s$ (importantly remember that the column index divided by 100 gives the value of $s$ i.e. the first column is 0.01, second column is 0.02, etc.). These matricies give us every possible combination of $(x,s)$.

InĀ [11]:
X_vals_2 = repeat(X,1,1000)
s_vals = X_vals_2';
InĀ [12]:
#### Subpart (ii): [will come back to below!] ####
InĀ [13]:
#### Subpart (iii): Create flow utility matrix U for each (x,s) as log (\max{x-s, \bar{c})#####

c_bar = 1e-10 # set c_bar given in problem

# Gives us a 1000 x 1000 matrix for a given par (x,s) #
U_vals = log.(max.(X_vals_2-s_vals,c_bar));
InĀ [14]:
#### Subpart (ii), (iv), & (v) ####

function iterate_V(V;policy = false)
    # Subpart (ii): for a given initial value function V, calculate the continuation value using VC_function defined above #
    # Input: Value function V over every possible value of s i.e. it is a 1000x1 vector #
    # Output w/ repeat function: 1000 x 1000 matrix with each column corrresponding to a value of (s) 
    VC_vals = repeat(VC_function(V)',1000,1) 

    # Subpart (iv) : Add flow utility + continuation value #
    V1_vals = U_vals + Ī“*VC_vals

    # Subpart (v) : Over every wealth of x (i.e. across every row since each row corresponds with a value of x), find the max value function  
    # value V1 AND also find the index location of where the max value occurs (this is important for using being able to find the amount 
    # of savings s that maximizes the value function)
    V1, temp = findmax(V1_vals,dims=2)
    if policy == true
        # Remember the relationship between the column index and the value of saving $s$. Here, we extract the column index for each row where 
        # the max value function occurs. And then, we divide by 100 which gives us the value of saving $s$ that maximizes the value function.
        return [s[2] for s in temp]./100 
    end
    return V1
end
Out[14]:
iterate_V (generic function with 1 method)

Part 3: Iterate on the Value function until convergenceĀ¶

InĀ [20]:
#### Subpart (a) ####

diff = 1000 # intialize difference variable
tol = 1e-8 # set tolerance as defined in problem

## Keep iterating until the difference between value function iterations gets very small ##
while diff > tol
    V1 = iterate_V(V0)
    diff = norm(V1-V0)
    V0=V1
end

V_final = V0;
InĀ [25]:
#### Subpart (b) ####

plot(X,V_final,legend = false)
Out[25]:
InĀ [31]:
#### Subpart (c) ####

# Extract the policy function from V_final #
saving = iterate_V(V_final, policy = true);

# Manually set saving when x = 0.01 to 0 - otherwise negative consumption #
saving[1] = 0

plot(X,saving,legend = false)
Out[31]:
InĀ [53]:
#### Subpart (d): Start with different guess and confirm same solution ####

diff = 1000 # Re-initialize the diff variable (!!!!!!!) - Don't forget to do this or the iteration won't work

V0 = fill(0, 1000, 1)

while diff > tol
    V1 = iterate_V(V0)
    diff = norm(V1-V0)
    V0=V1
end

V_final_v2 = V0;


# Plot this value function with other value function #
plot(X, V_final, legend=true, label="V_final")
plot!(X, V_final_v2, label="V_final_v2")
Out[53]: