Roadmap

1. How do we think about solving dynamic economic models
2. Value function iteration
3. Fixed point iteration
4. Time iteration
5. VFI + discretization

Things to do

1. Install: LinearAlgebra, Optim, Plots, Roots
2. Keep in mind that for VFI and TI we will be using optimization/rootfinding packages

• This matters because these packages typically only let the functions they work on have one input: the guesses for the maximizing input or root
• We get around this by expressing the function as a closure
• i.e. declare the function inside of a wrapper function that does the maximization/rootfinding so it can access the parameters in the wrapper function

1. Keep in mind we will be working with about the simplest example possible,
   more complex problems will be more difficult to solve in many ways

Solutions to economic models

How do we solve economic models?

First, what do we want?

We want to be able to compute things like optimal policy trajectories, welfare, etc

There are generally two objects that can deliver what we want:

1. Value functions
2. Policy functions

The idea behind the most commonly used solution concepts is
to recover good approximations to one of these two functions

Solutions to economic models

We recover these functions by exploiting two things:

1. Dynamic equilibrium conditions incorporating these functions (e.g. Bellman equations, Euler equations)
2. Fixed points

First lets look at recovering the value function

Our general example

Consider the following problem we will be using for all of these solution methods: $$\begin{gather} \max_{\left\{c_t \right\}_{t=0}^\infty} \sum_{t=1}^\infty \beta^t u(c_t) \notag \\ \text{subject to:} \,\,\,\,\, k_{t+1} = f(k_t) - c_t \notag \end{gather}$$ where both consumption and time $t+1$ capital are positive,
$k(0) = k_0$, $\alpha > 0$, and $\beta \in (0,1)$

Our general example

Represent the growth model as a Bellman equation $$\begin{gather} V(k) = \max_{c} u(c) + \beta V(k') \notag \\ \text{subject to:} \,\,\,\,\, k' = f(k) - c \notag \end{gather}$$

we can reduce this to $$V(k) = \max_{c} u(c) + \beta V(f(k) - c)$$

Method 1: Value function iteration

In VFI we approximate the value function with some flexible functional form $\Gamma(k_t;b)$ where $b$ is a vector of coefficients

The algorithm has 6 steps

Method 1: Value function iteration

Step 1: Select the number of collocation points in each dimension and the domain of the approximation space

Step 2: Select an initial vector of coefficients $b_0$ with the same number of elements as the collocation grid, and initial guesses for consumption for the solver

Step 3: Select a rule for convergence

Step 4: Construct the grid and basis matrix

Method 1: Value function iteration

Step 5: While convergence criterion $>$ tolerance (outer loop)

• Start iteration $p$
• For each grid point (inner loop)
• Solve the right hand side of the Bellman equation at each grid point using the
  value function approximant $\Gamma(k_{t+1};b^{(p)})$ in place of $V(k_{t+1})$
• Recover the maximized values at each grid point, conditional on the approximant
• Fit the polynomial to the values and recover a new vector of coefficients $\hat{b}^{(p+1)}$.
• Compute the vector of coefficients $b^{(p+1)}$ for iteration $p+1$ by
  $b^{(p+1)} = (1-\gamma) b^{(p)} + \gamma \hat{b}^{(p+1)}$ where $\gamma \in (0,1).$ (damping)
• Use the optimal controls for this iteration as our initial guess for next iteration

Step 6: Error check your approximation

Functional forms and parameters

Functional forms

• $u(c_t) = c_t^{1-\eta}/(1-\eta)$
• $f(k_t) = k_t^\alpha$

Parameters

• $\alpha = 0.75$
• $\beta = 0.95$
• $\eta = 2$

Initial capital value for simulating

• $k_0 = (\alpha \beta)^{1/(1-\alpha)}/2$

Step 1: Select the number of points and domain

If $k_0 = (\alpha \beta)^{1/(1-\alpha)}/2$ what are a logical set of bounds for the capital state?

$k^0$ and the steady state level $(\alpha \beta)^{1/(1-\alpha)}$

Put everything in a named tuple to make passing things easier

using LinearAlgebra
using Optim
using Plots
params = (alpha = 0.75, 
beta = 0.95, 
eta = 2,
steady_state = (0.75*0.95)^(1/(1 - 0.75)), 
k_0 = (0.75*0.95)^(1/(1 - 0.75))/2,
capital_upper = (0.75*0.95)^(1/(1 - 0.75))*1.01, 
capital_lower = (0.75*0.95)^(1/(1 - 0.75))/2,
num_points = 7, 
tolerance = 0.0001)

## (alpha = 0.75, beta = 0.95, eta = 2, steady_state = 0.25771486816406236, k_0 = 0.12885743408203118, capital_upper = 0.26029201684570297, capital_lower = 0.12885743408203118, num_points = 7, tolerance = 0.0001)

Step 2: Select an initial vector of coefficients $b_0$

In some cases you might have a good guess (e.g. increasing and concave so you know the second value is positive, third value is negative, rest maybe set to zero)

Other cases you might not,
guessing zeros effectively turns the initial iteration into a static problem,
the second iteration into a 2 period problem, and so on

coefficients = zeros(params.num_points)

## 7-element Array{Float64,1}:
##  0.0
##  0.0
##  0.0
##  0.0
##  0.0
##  0.0
##  0.0

Step 3: Select a convergence rule

There's a lot of potential options here to determine convergence of the function

Relative or absolute change? Or both?

Change in the value function? Change in the policy function?

Which norm?

Our rule for class: convergence is when the maximum relative change in value on the grid is < 0.001%

Step 4: Construct the grid and basis matrix

The function cheb_nodes constructs the grid on $[-1,1]$

$$x_k = cos\left(\frac{2k-1}{2n}\pi\right),\,\, k = 1,...,n$$

cheb_nodes(n) = cos.(pi * (2*(1:n) .- 1)./(2n))

## cheb_nodes (generic function with 1 method)

grid = cheb_nodes(params.num_points) # [-1, 1] grid

## 7-element Array{Float64,1}:
##   0.9749279121818236   
##   0.7818314824680298   
##   0.4338837391175582   
##   6.123233995736766e-17
##  -0.43388373911755806  
##  -0.7818314824680297   
##  -0.9749279121818236

Step 4: Construct the grid and basis matrix

But we need to expand the grid from $[-1,1]$ to our actual capital domain

expand_grid(grid, params) = (1 .+ grid)*(params.capital_upper - params.capital_lower)/2 .+ params.capital_lower

## expand_grid (generic function with 1 method)

capital_grid = expand_grid(grid, params)

## 7-element Array{Float64,1}:
##  0.2586443471450049 
##  0.2459545728087113 
##  0.2230883895732961 
##  0.19457472546386706
##  0.16606106135443804
##  0.14319487811902284
##  0.13050510378272925

Step 4: Construct the grid and basis matrix

Use cheb_polys to construct the basis matrix

# Chebyshev polynomial function
function cheb_polys(x, n)
if n == 0
return 1                    # T_0(x) = 1
elseif n == 1
return x                    # T_1(x) = x
else
cheb_recursion(x, n) =
2x.*cheb_polys.(x, n - 1) .- cheb_polys.(x, n - 2)
return cheb_recursion(x, n) # T_n(x) = 2xT_{n-1}(x) - T_{n-2}(x)
end
end;

Step 4a: Pre-invert your basis matrix

Pro tip: you will be using the exact same basis matrix in each loop iteration to recover the coefficients: just pre-invert it to save time because inverting the same matrix every loop is costly (especially when large)

basis_matrix = [cheb_polys.(grid, n) for n = 0:params.num_points - 1];
basis_matrix = hcat(basis_matrix...)

## 7×7 Array{Float64,2}:
##  1.0   0.974928      0.900969  …   0.62349    0.433884      0.222521
##  1.0   0.781831      0.222521     -0.900969  -0.974928     -0.62349 
##  1.0   0.433884     -0.62349      -0.222521   0.781831      0.900969
##  1.0   6.12323e-17  -1.0           1.0        3.06162e-16  -1.0     
##  1.0  -0.433884     -0.62349      -0.222521  -0.781831      0.900969
##  1.0  -0.781831      0.222521  …  -0.900969   0.974928     -0.62349 
##  1.0  -0.974928      0.900969      0.62349   -0.433884      0.222521

basis_inverse = basis_matrix\I

## 7×7 Array{Float64,2}:
##  0.142857    0.142857    0.142857   …   0.142857    0.142857    0.142857 
##  0.278551    0.22338     0.123967      -0.123967   -0.22338    -0.278551 
##  0.25742     0.0635774  -0.17814       -0.17814     0.0635774   0.25742  
##  0.22338    -0.123967   -0.278551       0.278551    0.123967   -0.22338  
##  0.17814    -0.25742    -0.0635774     -0.0635774  -0.25742     0.17814  
##  0.123967   -0.278551    0.22338    …  -0.22338     0.278551   -0.123967 
##  0.0635774  -0.17814     0.25742        0.25742    -0.17814     0.0635774

Step 4b: Evaluate the continuation value

To loop and maximize the Bellman at each grid point we need a function eval_value_function(coefficients, capital, params) that lets us evaluate the continuation value given a vector of coefficients coefficients, a vector of capital nodes capital, and the parameters params to scale capital back into $[-1,1]$

It needs to:

1. Scale capital back into $[-1,1]$
2. Use the coefficients and Chebyshev polynomials to evaluate the value function

Step 4b: Evaluate the continuation value

Here's a simple way to do it:

shrink_grid(capital) = 2*(capital - params.capital_lower)/(params.capital_upper - params.capital_lower) - 1;
eval_value_function(coefficients, capital, params) = 
coefficients' * [cheb_polys.(shrink_grid(capital), n) for n = 0:params.num_points - 1];

The top function inherits params from the bottom function

Step 5: Inner loop over grid points

Construct a function that loops over the grid points
and solves the Bellman given $\Gamma(x;b^{(p)})$

Pseudocode:

for each grid point:
define the Bellman as a closure so it can take in parameters
maximize the Bellman by choosing consumption with the optimize function
store maximize value in a vector
end
return vector of maximized values

Step 5: Inner loop over grid points

function loop_grid(params, basis_inverse, basis_matrix, grid, capital_grid, coefficients)
max_value = similar(coefficients); # initialized max value vector
# Inner loop over grid points
for (iteration, capital) in enumerate(capital_grid)
# Define Bellman as a closure
function bellman(consumption)
capital_next = capital^params.alpha - consumption # Next period state
cont_value = eval_value_function(coefficients, capital_next, params) # Continuation value
value_out = (consumption)^(1-params.eta)/(1-params.eta) + params.beta*cont_value # Utility + continuation value
return -value_out
end;
results = optimize(bellman, 0.00*capital^params.alpha, 0.99*capital^params.alpha) # maximize Bellman
max_value[iteration] = -Optim.minimum(results) # Store max value in vector
end
return max_value
end

## loop_grid (generic function with 1 method)

Step 5: Outer loop iterating on Bellman

function solve_vfi(params, basis_inverse, basis_matrix, grid, capital_grid, coefficients)
iteration = 1
error = 1e10;
max_value = similar(coefficients);
value_prev = .1*ones(params.num_points);
coefficients_store = Vector{Vector}(undef, 1)
coefficients_store[1] = coefficients
while error > params.tolerance # Outer loop iterating on Bellman eq
max_value = loop_grid(params, basis_inverse, basis_matrix, grid, capital_grid, coefficients) # Inner loop
coefficients = basis_inverse*max_value # \Psi \ y, recover coefficients
error = maximum(abs.((max_value - value_prev)./(value_prev))) # compute error
value_prev = deepcopy(max_value) # save previous values
if mod(iteration, 5) == 0
println("Maximum Error of $(error) on iteration $(iteration).")
append!(coefficients_store, [coefficients])
end
iteration += 1
end
return coefficients, max_value, coefficients_store
end

## solve_vfi (generic function with 1 method)

Step 5: Outer loop iterating on Bellman

solution_coeffs, max_value, intermediate_coefficients =
solve_vfi(params, basis_inverse, basis_matrix, grid, capital_grid, coefficients)

## Maximum Error of 0.3301919884226087 on iteration 5.
## Maximum Error of 0.1080139919745119 on iteration 10.
## Maximum Error of 0.05647917855011529 on iteration 15.
## Maximum Error of 0.034833389245083224 on iteration 20.
## Maximum Error of 0.02328643376111153 on iteration 25.
## Maximum Error of 0.016301543092581427 on iteration 30.
## Maximum Error of 0.011747480470414 on iteration 35.
## Maximum Error of 0.008631245645920336 on iteration 40.
## Maximum Error of 0.006427690604127001 on iteration 45.
## Maximum Error of 0.0048330736842455225 on iteration 50.
## Maximum Error of 0.0036597148900624123 on iteration 55.
## Maximum Error of 0.002785692376976503 on iteration 60.
## Maximum Error of 0.0021286867102129677 on iteration 65.
## Maximum Error of 0.0016314249677370316 on iteration 70.
## Maximum Error of 0.0012531160571391219 on iteration 75.
## Maximum Error of 0.000964170879127405 on iteration 80.
## Maximum Error of 0.0007428166750767927 on iteration 85.
## Maximum Error of 0.0005728521498806172 on iteration 90.
## Maximum Error of 0.0004421161961630202 on iteration 95.
## Maximum Error of 0.00034141814527728653 on iteration 100.
## Maximum Error of 0.0002637753971325223 on iteration 105.
## Maximum Error of 0.00020386108209645033 on iteration 110.
## Maximum Error of 0.00015759845903308656 on iteration 115.
## Maximum Error of 0.00012185979301915734 on iteration 120.

## ([-200.6291758538632, 9.991472391067802, -1.2278992641150595, 0.17379460460092133, -0.0262119101944549, 0.003954006913225783, -0.0007409750421256689], [-191.87342361439286, -193.14594489252312, -195.68963652528421, -199.42674752490046, -204.02721963808617, -208.61071718644223, -211.6305415954132], Array{T,1} where T[[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [-36.8589490106075, 6.259861237095166, -0.8221318972861571, 0.11882022830526928, -0.01784945487786871, 0.00269203315390687, -0.00037520235134058666], [-74.87645132736965, 8.931829874526976, -1.1238227131298792, 0.1608047649285531, -0.024292798807126204, 0.0039027323348115317, -0.0006447123917387998], [-103.80959945589598, 9.68901311684359, -1.1994026067805663, 0.1705609554663461, -0.0255481405498017, 0.003969430583547151, -0.0007161344250450341], [-125.9200308753978, 9.904403319377032, -1.219552059022945, 0.17300404781584788, -0.026010058299306138, 0.0039595903769296115, -0.0007341995822489622], [-142.9397465762615, 9.966570886453908, -1.2253830534940988, 0.1735771626051843, -0.02615385493585265, 0.003955453304666179, -0.0007391283603670473], [-156.08276002781142, 9.984453747414598, -1.2271575579534755, 0.1737309439806367, -0.02619527771992125, 0.003954320282257839, -0.0007404618310982158], [-166.2448607828197, 9.989518053393965, -1.2276872050894383, 0.1737761442042398, -0.026207193540553675, 0.003954072102075656, -0.0007408314670129812], [-174.10590746285365, 9.990932110798063, -1.2278398648722586, 0.17378937493059254, -0.026210590152359714, 0.003954021287697174, -0.0007409350663660774], [-180.18802216175297, 9.991323564318797, -1.2278828039994707, 0.1737931465991167, -0.026211544200975823, 0.003954010369653815, -0.0007409639787887556]  …  [-196.52576254532732, 9.991472327014705, -1.2278992569992653, 0.17379460396753643, -0.02621191003608203, 0.0039540069145189705, -0.0007409750373632562], [-197.5358441489631, 9.991472373473798, -1.2278992621604772, 0.17379460442691652, -0.026211910150955475, 0.0039540069135881595, -0.0007409750408040594], [-198.3174260190019, 9.991472386235223, -1.2278992635781947, 0.17379460455315154, -0.026211910182489362, 0.00395400691330039, -0.0007409750417757266], [-198.92219916559213, 9.991472389740558, -1.2278992639676147, 0.17379460458784024, -0.02621191019120772, 0.003954006913229335, -0.0007409750420261929], [-199.39016109641344, 9.991472390703436, -1.2278992640745656, 0.17379460459731177, -0.026211910193559618, 0.0039540069132151245, -0.0007409750420777073], [-199.75226111754048, 9.991472390967878, -1.227899264103975, 0.17379460459993368, -0.02621191019424174, 0.003954006913200914, -0.0007409750421203398], [-200.0324472112422, 9.991472391040567, -1.2278992641120254, 0.17379460460063711, -0.026211910194398058, 0.003954006913229335, -0.0007409750421185635], [-200.24924986946957, 9.991472391060476, -1.227899264114292, 0.1737946046008929, -0.026211910194433585, 0.0039540069131902555, -0.0007409750420972472], [-200.41700763359654, 9.99147239106594, -1.2278992641148818, 0.1737946046008929, -0.026211910194390953, 0.00395400691322223, -0.0007409750421132344], [-200.54681539359328, 9.991472391067411, -1.2278992641150026, 0.17379460460090712, -0.0262119101944549, 0.003954006913208019, -0.0007409750421345507]])

Now lets plot our solutions

capital_levels = range(params.capital_lower, params.capital_upper, length = 100);
eval_points = shrink_grid.(capital_levels);
solution = similar(intermediate_coefficients);
# Compute optimal value at all capital grid points
for (iteration, coeffs) in enumerate(intermediate_coefficients)
solution[iteration] = [coeffs' * [cheb_polys.(capital, n) for n = 0:params.num_points - 1] for capital in eval_points];
end

Plot the value function iterations

Plot the value function iterations

Plot the value function iterations

Plot the value function iterations

Plot the value function iterations

Plot the value function iterations

Plot the value function iterations

Plot the value function iterations

Plot the value function iterations

Plot the final value function

Now lets try simulating

function simulate_model(params, solution_coeffs, time_horizon = 100)
capital_store = zeros(time_horizon + 1)
consumption_store = zeros(time_horizon)
capital_store[1] = params.k_0
for t = 1:time_horizon
capital = capital_store[t]
function bellman(consumption)
capital_next = capital^params.alpha - consumption
capital_next_scaled = shrink_grid(capital_next)
cont_value = solution_coeffs' * [cheb_polys.(capital_next_scaled, n) for n = 0:params.num_points - 1]
value_out = (consumption)^(1-params.eta)/(1-params.eta) + params.beta*cont_value
return -value_out
end;
results = optimize(bellman, 0.0, capital^params.alpha)
consumption_store[t] = Optim.minimizer(results)
capital_store[t+1] = capital^params.alpha - consumption_store[t]
end
return consumption_store, capital_store
end;

Now lets try simulating

The consumption policy function

capital_levels = range(params.capital_lower, params.capital_upper, length = 100);
consumption = similar(capital_levels);
# Compute optimal consumption at all capital grid points
for (iteration, capital) in enumerate(capital_levels)
function bellman(consumption)
capital_next = capital^params.alpha - consumption
capital_next_scaled = shrink_grid(capital_next)
cont_value = solution_coeffs' * [cheb_polys.(capital_next_scaled, n) for n = 0:params.num_points - 1]
value_out = (consumption)^(1-params.eta)/(1-params.eta) + params.beta*cont_value
return -value_out
end
results = optimize(bellman, 0., capital^params.alpha)
consumption[iteration] = Optim.minimizer(results)
end;

The consumption policy function

Method 2: Fixed point iteration

In FPI we generally approximate a policy function with some flexible functional form $\Gamma(k_t;b)$ where $b$ is a vector of coefficients

FPI re-casts equilibrium conditions of the model as a fixed point

We then perform multi-dimensional function iteration to solve for the fixed point

This ends up being very simple and it works on any dimension function

It is also does not bear a terrible computational cost and is derivative-free

The drawback is that it will not always converge and is generally unstable

This can be solved by damping

Eq condition: Euler equation

Often we will iterate on the Euler equation which for our problem is $$u'(c_t) = \beta u'(c_{t+1}) f'(k_{t+1})$$

We need to put this in a fixed point form in order to iterate on it $$c_t = u'^{(-1)}\left(\beta u'(c_{t+1}) f'(k_{t+1})\right)$$

How do we solve this?

Eq condition: Euler equation

$$c_t = u'^{(-1)}\left(\beta u'(c_{t+1}) f'(k_{t+1})\right)$$

We approximate the consumption policy function $c_{t} = C(k_t)$ with some flexible functional form $\Psi(k_t; b)$

We have defined $c_{t}$ in two ways, once as an outcome of the policy function, and once as an equilibrium condition

Now we can form our consumption policy function as a fixed point by substituting $C(k_t)$ into the the Euler fixed point as follows $$C(k_t) = u'^{(-1)}\left(\beta u'(C(k_{t+1})) f'(k_{t+1}(C(k_t),k_t))\right)$$

Method 2: Fixed point iteration

The algorithm has 6 steps, very similar to VFI

Method 2: Fixed point iteration

Step 1: Select the number of collocation points in each dimension and the domain of the approximation space

Step 2: Select an initial vector of coefficients $b_0$ with the same number of elements as the collocation grid

Step 3: Select a rule for convergence

Step 4: Construct the grid and basis matrix

Method 2: Fixed point iteration

Step 5: While convergence criterion $>$ tolerance (outer loop)

• Start iteration $p$
• For each grid point (inner loop)
  • Substitute $C(k_{t+1};b^{(p)})$ into the right hand side of the Euler fixed point
  • Recover the LHS values of consumption at each grid point
• Fit the polynomial to the values and recover a new vector of coefficients $\hat{b}^{(p+1)}$.
• Compute the vector of coefficients $b^{(p+1)}$ for iteration $p+1$ by $b^{(p+1)} = (1-\gamma) b^{(p)} + \gamma \hat{b}^{(p+1)}$ where $\gamma \in (0,1).$ (damping)

Step 6: Error check your approximation

Notice: we did not have to perform a maximization step anywhere,
this leads to big speed gains

Step 1: Select the number of points and domain

Put everything in a named tuple to make passing things easier

using LinearAlgebra
using Optim
using Plots
params_fpi = (alpha = 0.75, beta = 0.95, eta = 2, damp = 0.7,
                steady_state = (0.75*0.95)^(1/(1-0.75)), k_0 = (0.75*0.95)^(1/(1-0.75))*0.5,
                capital_upper = (0.75*0.95)^(1/(1-0.75))*1.5, capital_lower = (0.75*0.95)^(1/(1-0.75))*0.5,
                num_points = 5, tolerance = 0.00001)

## (alpha = 0.75, beta = 0.95, eta = 2, damp = 0.7, steady_state = 0.25771486816406236, k_0 = 0.12885743408203118, capital_upper = 0.3865723022460935, capital_lower = 0.12885743408203118, num_points = 5, tolerance = 1.0e-5)

shrink_grid(capital) = 2*(capital - params_fpi.capital_lower)/(params_fpi.capital_upper - params_fpi.capital_lower) - 1

## shrink_grid (generic function with 1 method)

Step 2: Select an initial vector of coefficients $b_0$

coefficients = zeros(params_fpi.num_points)

## 5-element Array{Float64,1}:
##  0.0
##  0.0
##  0.0
##  0.0
##  0.0

Step 3: Select a convergence rule

Rule: maximum change in consumption on the grid < 0.001%

Step 4: Construct the grid and basis matrix

The function cheb_nodes constructs the grid on $[-1,1]$

$$x_k = cos\left(\frac{2k-1}{2n}\pi\right),\,\, k = 1,...,n$$

cheb_nodes(n) = cos.(pi * (2*(1:n) .- 1)./(2n))

## cheb_nodes (generic function with 1 method)

grid = cheb_nodes(params_fpi.num_points) # [-1, 1] grid

## 5-element Array{Float64,1}:
##   0.9510565162951535   
##   0.5877852522924731   
##   6.123233995736766e-17
##  -0.587785252292473    
##  -0.9510565162951535

Step 4: Construct the grid and basis matrix

But we need to expand the grid from $[-1,1]$ to our actual capital domain

expand_grid(grid, params_fpi) = (1 .+ grid)*(params_fpi.capital_upper - params_fpi.capital_lower)/2 .+ params_fpi.capital_lower

## expand_grid (generic function with 1 method)

capital_grid = expand_grid(grid, params_fpi)

## 5-element Array{Float64,1}:
##  0.3802655705208513 
##  0.33345536756572974
##  0.2577148681640623 
##  0.18197436876239492
##  0.1351641658072734

Step 4: Construct the grid and basis matrix

Use cheb_polys to construct the basis matrix

# Chebyshev polynomial function
function cheb_polys(x, n)
    if n == 0
        return 1                    # T_0(x) = 1
    elseif n == 1
        return x                    # T_1(x) = x
    else
        cheb_recursion(x, n) =
        2x.*cheb_polys.(x, n - 1) .- cheb_polys.(x, n - 2)
        return cheb_recursion(x, n) # T_n(x) = 2xT_{n-1}(x) - T_{n-2}(x)
    end
end;

Step 4a: Pre-invert your basis matrix

Pro tip: you will be using the exact same basis matrix in each loop iteration to recover the coefficients: just pre-invert it to save time because inverting the same matrix every loop is costly (especially when large)

basis_matrix = [cheb_polys.(grid, n) for n = 0:params_fpi.num_points - 1];
basis_matrix = hcat(basis_matrix...)

## 5×5 Array{Float64,2}:
##  1.0   0.951057      0.809017   0.587785      0.309017
##  1.0   0.587785     -0.309017  -0.951057     -0.809017
##  1.0   6.12323e-17  -1.0       -1.83697e-16   1.0     
##  1.0  -0.587785     -0.309017   0.951057     -0.809017
##  1.0  -0.951057      0.809017  -0.587785      0.309017

basis_inverse = basis_matrix\I

## 5×5 Array{Float64,2}:
##  0.2        0.2        0.2           0.2        0.2     
##  0.380423   0.235114   5.55112e-17  -0.235114  -0.380423
##  0.323607  -0.123607  -0.4          -0.123607   0.323607
##  0.235114  -0.380423  -8.32667e-17   0.380423  -0.235114
##  0.123607  -0.323607   0.4          -0.323607   0.123607

Step 4b: Evaluate the policy function

To loop and maximize the Bellman at each grid point we need a function eval_policy_function(coefficients, capital, params_fpi) that lets us evaluate the policy function given a vector of coefficients coefficients, a vector of capital nodes capital, and the parameters params_fpi to scale capital back into $[-1,1]$

It needs to:

1. Scale capital back into $[-1,1]$
2. Use the coefficients and Chebyshev polynomials to evaluate the value function

Step 4b: Evaluate the policy function

Here's a simple way to do it:

shrink_grid(capital) = 2*(capital - params_fpi.capital_lower)/(params_fpi.capital_upper - params_fpi.capital_lower) - 1;
eval_policy_function(coefficients, capital, params_fpi) = 
    coefficients' * [cheb_polys.(shrink_grid(capital), n) for n = 0:params_fpi.num_points - 1];

The top function inherits params_fpi from the bottom function

Step 5: Loop

Construct the Euler fixed point function

function consumption_euler(params_fpi, capital, coefficients)
    # RHS: Current consumption given current capital
    consumption = eval_policy_function(coefficients, capital, params_fpi)
    # RHS: Next period's capital given current capital and consumption
    capital_next = capital^params_fpi.alpha - consumption
    # RHS: Next period's consumption given current capital and consumption
    consumption_next = eval_policy_function(coefficients, capital_next, params_fpi)
    consumption_next = max(1e-10, consumption_next)
    # LHS: Next period's consumption from Euler equation
    consumption_lhs = (
        params_fpi.beta *
        consumption_next^(-params_fpi.eta) *
        params_fpi.alpha*(capital_next).^(params_fpi.alpha-1)
        ).^(-1/params_fpi.eta)
    return consumption_lhs
end

## consumption_euler (generic function with 1 method)

Step 5: Loop

Construct a function that loops over the grid points and solves the Euler given $\Psi(x;b^{(p)})$

function loop_grid_fpi(params_fpi, basis_inverse, basis_matrix, grid, capital_grid, coefficients)
    consumption = similar(coefficients)
    # Compute next period's consumption from the Euler equation
    for (iteration, capital) in enumerate(capital_grid)
        consumption[iteration] = consumption_euler(params_fpi, capital, coefficients)
    end
    return consumption
end

## loop_grid_fpi (generic function with 1 method)

Step 5: Loop

function solve_fpi(params_fpi, basis_inverse, basis_matrix, grid, capital_grid, coefficients)
    error = 1e10
    iteration = 1
    consumption = similar(coefficients)
    consumption_prev, coefficients_prev = similar(coefficients), similar(coefficients)
    coefficients_store = Vector{Vector}(undef, 1)
    coefficients_store[1] = coefficients
    while error > params_fpi.tolerance
        consumption = loop_grid_fpi(params_fpi, basis_inverse, basis_matrix, grid, capital_grid, coefficients)
        if iteration > 1
            coefficients = params_fpi.damp*(basis_inverse*consumption) + (1 - params_fpi.damp)*coefficients_prev
        else
            coefficients = basis_inverse*consumption
        end
        error = maximum(abs.((consumption - consumption_prev)./(consumption_prev)))
        coefficients_prev, consumption_prev = deepcopy(coefficients), deepcopy(consumption)
        if mod(iteration, 5) == 0
            println("Maximum Error of $(error) on iteration $(iteration).")
            append!(coefficients_store, [coefficients])
        end
        iteration += 1
    end
    return coefficients, consumption, coefficients_store
end

## solve_fpi (generic function with 1 method)

Step 5: Loop

solution_coeffs, consumption, intermediate_coefficients =
    solve_fpi(params_fpi, basis_inverse, basis_matrix, grid, capital_grid, coefficients)

## Maximum Error of 0.06899533243794735 on iteration 5.
## Maximum Error of 1.481484091434361 on iteration 10.
## Maximum Error of 0.7160886352974187 on iteration 15.
## Maximum Error of 0.3030829029947649 on iteration 20.
## Maximum Error of 0.1747314280387973 on iteration 25.
## Maximum Error of 505.49415653346773 on iteration 30.
## Maximum Error of 1.193581540407262 on iteration 35.
## Maximum Error of 0.2847765742028598 on iteration 40.
## Maximum Error of 0.1953733398047281 on iteration 45.
## Maximum Error of 236164.34258207615 on iteration 50.
## Maximum Error of 1.4738022100415356 on iteration 55.
## Maximum Error of 0.38721294652038635 on iteration 60.
## Maximum Error of 0.21074834398476477 on iteration 65.
## Maximum Error of 0.23613685495903153 on iteration 70.
## Maximum Error of 4.722253097188453 on iteration 75.
## Maximum Error of 0.03764118069395616 on iteration 80.
## Maximum Error of 0.007970939090958102 on iteration 85.
## Maximum Error of 0.0025735627849945903 on iteration 90.
## Maximum Error of 0.0009352774342816211 on iteration 95.
## Maximum Error of 0.0003237453048973147 on iteration 100.
## Maximum Error of 0.00010972350926806329 on iteration 105.
## Maximum Error of 3.685905281539792e-5 on iteration 110.
## Maximum Error of 1.2336267929334763e-5 on iteration 115.

## ([0.10216143756493291, 0.030492616130876768, -0.0018615048468132258, 0.0002412084612514174, -3.673979487435926e-5], [0.12978669111539357, 0.12046064710918444, 0.10398662400247835, 0.08507298585014693, 0.07150233223860529], Array{T,1} where T[[0.0, 0.0, 0.0, 0.0, 0.0], [1.2763130221485964e-10, 7.183782174427739e-12, -5.1738097170530396e-12, -2.5864997862853625e-13, -8.831154889933433e-14], [5.317142229921049e-10, 5.527847778513712e-10, 1.6703532060316548e-10, 1.8141243875066817e-11, -1.447815296122337e-13], [8.316967352701557e-9, 5.1875398530251166e-9, 2.0308557486015783e-10, -1.1491354530189611e-11, 1.7095212937975535e-12], [3.386275892864933e-8, 1.2764008226603924e-8, -3.091246647800147e-10, 4.0250161382908916e-11, -1.1414762945221957e-11], [7.52422641928251e-8, 1.0570161284975366e-8, -4.436605014141635e-9, -1.174432549622199e-10, -4.947771112908894e-11], [3.5815259527248966e-8, -1.4472091279968585e-8, 1.834736552630338e-8, 5.228523646511384e-9, -1.5039924808248198e-10], [8.649742496094897e-7, 6.24025102132348e-7, 2.5809596059877958e-8, -5.853265390446539e-9, 3.605947196166565e-10], [3.557675838337089e-6, 1.2735910170908482e-6, -3.0722078466656094e-8, 1.0900582853788604e-8, -1.2927503422953689e-9], [8.396083635269678e-6, 1.5757179125057394e-6, -4.34579079208472e-7, -2.0463275140734446e-8, -5.01406374442219e-9]  …  [0.014880293346307229, -0.0012599250822342666, -0.0016020724342596764, -7.459264578786532e-5, -2.116862090455988e-5], [0.06056412987240832, 0.05723653749340376, -0.002784166148579988, -0.004137596592328633, 0.0038802212866255025], [0.09274367453966734, 0.02965908339159335, -0.0016114105549141255, 0.00024575125177950144, -3.150220775509663e-5], [0.09883950641207304, 0.02942554740820863, -0.001774882853544672, 0.00023321088845588143, -3.695484542590681e-5], [0.10100737679048473, 0.030003914490261078, -0.0018414718812655218, 0.00023999644919096, -3.7123219490591355e-5], [0.10177119243335409, 0.03030979782120775, -0.0018570023721880147, 0.00024112552978790023, -3.6907250737000245e-5], [0.10203276998407984, 0.03042983416898875, -0.0018604002640126983, 0.00024122771243408221, -3.679946836529483e-5], [0.10212100280888622, 0.03047255019091084, -0.0018612121452150179, 0.0002412205464725234, -3.675909294053858e-5], [0.10215057465288752, 0.03048718729579106, -0.001861432535354346, 0.00024121237383887737, -3.6745039179336284e-5], [0.10216046106906573, 0.030492126707819756, -0.0018614985832312249, 0.00024120883727619163, -3.674026850712608e-5]])

Now lets plot our solutions

capital_levels = range(params_fpi.capital_lower, params_fpi.capital_upper, length = 100);
eval_points = shrink_grid.(capital_levels);
solution = similar(intermediate_coefficients);
for (iteration, coeffs) in enumerate(intermediate_coefficients)
    solution[iteration] = [coeffs'*[cheb_polys.(capital, n) for n = 0:params_fpi.num_points - 1] for capital in eval_points];
end

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the final consumption policy function

Now lets try simulating

function simulate_model(params, solution_coeffs, time_horizon = 100)
    capital_store = zeros(time_horizon + 1)
    consumption_store = zeros(time_horizon)
    capital_store[1] = params.k_0
    for t = 1:time_horizon
        capital = capital_store[t]
        consumption_store[t] = consumption_euler(params, capital, solution_coeffs)
        capital_store[t+1] = capital^params.alpha - consumption_store[t]
    end
    return consumption_store, capital_store
end

## simulate_model (generic function with 2 methods)

Now lets try simulating

time_horizon = 100;
consumption, capital = simulate_model(params_fpi, solution_coeffs, time_horizon);
plot(1:time_horizon, consumption, color = :red, linewidth = 4.0, tickfontsize = 14, guidefontsize = 14, label = "Consumption", legend = :right, grid = false, size = (500, 300));
plot!(1:time_horizon, capital[1:end-1], color = :blue, linewidth = 4.0, linestyle = :dash, label = "Capital");
plot!(1:time_horizon, params_fpi.steady_state*ones(time_horizon), color = :purple, linewidth = 2.0, linestyle = :dot, label = "Steady State Capital")

Method 3: Time iteration

In time iteration we approximate the policy function with some flexible functional form $\Psi(k_t;b)$ where $b$ is a vector of coefficients

The difference vs FPI is we use root-finding techniques on our $n$ node collocation grid where we search for the scalar $c^{(p+1)}(k_t)$ that solves $$u'(c^{(p+1)}(k^j_t)) = \beta u'(C^{(p)}(f(k^j_t)-c^{(p+1)}(k^j_t))) f'(f(k^i_t)-c^{(p+1)}(k^j_t)) \,\,\,\, \text{for j } = 1,...,n$$

$C^{(p)}()$ is our current approximation to the policy function, and we are searching for a scalar $c^{(p+1)}(k^j_t)$, given our collocation node $k_t^j$, that solves the Euler equation root-finding problem

In the Euler equation $c^{(p+1)}$ corresponds to today's policy function while $C^{(p)}$ corresponds to tomorrow's policy function: we are searching for today's policy that satisfies the Euler equation

Method 3: Time iteration

As we iterate and $p$ increases, $C^{(p)}(k)$ should converge because of a monotonicity property

If $C'^{(p)}(k) > 0$, and $C^{(p)}(k) < C^{(p-1)}(k)$, then $C^{(p+1)}(k) < C^{(p)}(k)$ and $C'^{(p+1)}(k) > 0$

It preserves the (first-order) shape of the policy function so it is reliable and convergent

Unfortunately time iteration tends to be slow,
especially as the number of dimensions grows

Method 3: Time iteration

Step 1: Select the number of collocation points in each dimension and the domain of the approximation space

Step 2: Select an initial vector of coefficients $b_0$ with the same number of elements as the collocation grid, and initial guesses for consumption for the root-finder

Step 3: Select a rule for convergence

Step 4: Construct the grid and basis matrix

Method 3: Time iteration

Step 5: While convergence criterion $>$ tolerance (outer loop)

• Start iteration $p$
• For each grid point (inner loop)
  • Substitute $C(k^j_{t+1};b^{(p)})$ in for $t+1$ consumption
  • Recover the $c^{(p+1)}(k^j_t) \in \mathbb{R}$ scalar values that satisfy the equation
• Fit the polynomial to the values and recover a new vector of coefficients $\hat{b}^{(p+1)}$
• Compute the vector of coefficients $b^{(p+1)}$ for iteration $p+1$ by $b^{(p+1)} = (1-\gamma) b^{(p)} + \gamma \hat{b}^{(p+1)}$ where $\gamma \in (0,1)$ (damping)
• Error check your approximation

Step 1: Select the number of points and domain

Put everything in a named tuple to make passing things easier

using LinearAlgebra
using Optim
using Plots
using Roots
params_ti = (alpha = 0.75, beta = 0.95, eta = 2, damp = 0.7,
                steady_state = (0.75*0.95)^(1/(1-0.75)), k_0 = (0.75*0.95)^(1/(1-0.75))*0.5,
                capital_upper = (0.75*0.95)^(1/(1-0.75))*1.5, capital_lower = (0.75*0.95)^(1/(1-0.75))*0.5,
                num_points = 5, tolerance = 0.00001)

## (alpha = 0.75, beta = 0.95, eta = 2, damp = 0.7, steady_state = 0.25771486816406236, k_0 = 0.12885743408203118, capital_upper = 0.3865723022460935, capital_lower = 0.12885743408203118, num_points = 5, tolerance = 1.0e-5)

shrink_grid(capital) = 2*(capital - params_ti.capital_lower)/(params_ti.capital_upper - params_ti.capital_lower) - 1

## shrink_grid (generic function with 1 method)

Step 2: Select an initial vector of coefficients $b_0$

In some cases you might have a good guess (e.g. increasing and concave so you know the second value is positive, third value is negative, rest maybe set to zero)

Other cases you might not,
guessing zeros effectively turns the initial iteration into a static problem,
the second iteration into a 2 period problem, and so on

coefficients = zeros(params_ti.num_points)

## 5-element Array{Float64,1}:
##  0.0
##  0.0
##  0.0
##  0.0
##  0.0

Step 3: Select a convergence rule

There's a lot of potential options here to determine convergence of the function

Relative or absolute change? Or both?

Change in the value function? Change in the policy function?

Which norm?

Our rule for class: convergence is when the maximum relative change in value on the grid is < 0.001%

Step 4: Construct the grid and basis matrix

The function cheb_nodes constructs the grid on $[-1,1]$

$$x_k = cos\left(\frac{2k-1}{2n}\pi\right),\,\, k = 1,...,n$$

cheb_nodes(n) = cos.(pi * (2*(1:n) .- 1)./(2n))

## cheb_nodes (generic function with 1 method)

grid = cheb_nodes(params_ti.num_points) # [-1, 1] grid

## 5-element Array{Float64,1}:
##   0.9510565162951535   
##   0.5877852522924731   
##   6.123233995736766e-17
##  -0.587785252292473    
##  -0.9510565162951535

Step 4: Construct the grid and basis matrix

But we need to expand the grid from $[-1,1]$ to our actual capital domain

expand_grid(grid, params_ti) = (1 .+ grid)*(params_ti.capital_upper - params_ti.capital_lower)/2 .+ params_ti.capital_lower

## expand_grid (generic function with 1 method)

capital_grid = expand_grid(grid, params_ti)

## 5-element Array{Float64,1}:
##  0.3802655705208513 
##  0.33345536756572974
##  0.2577148681640623 
##  0.18197436876239492
##  0.1351641658072734

Step 4: Construct the grid and basis matrix

Use cheb_polys to construct the basis matrix

# Chebyshev polynomial function
function cheb_polys(x, n)
    if n == 0
        return 1                    # T_0(x) = 1
    elseif n == 1
        return x                    # T_1(x) = x
    else
        cheb_recursion(x, n) =
            2x.*cheb_polys.(x, n - 1) .- cheb_polys.(x, n - 2)
        return cheb_recursion(x, n) # T_n(x) = 2xT_{n-1}(x) - T_{n-2}(x)
    end
end;

Step 4a: Pre-invert your basis matrix

Pro tip: you will be using the exact same basis matrix in each loop iteration to recover the coefficients: just pre-invert it to save time because inverting the same matrix every loop is costly (especially when large)

basis_matrix = [cheb_polys.(grid, n) for n = 0:params_ti.num_points - 1];
basis_matrix = hcat(basis_matrix...)

## 5×5 Array{Float64,2}:
##  1.0   0.951057      0.809017   0.587785      0.309017
##  1.0   0.587785     -0.309017  -0.951057     -0.809017
##  1.0   6.12323e-17  -1.0       -1.83697e-16   1.0     
##  1.0  -0.587785     -0.309017   0.951057     -0.809017
##  1.0  -0.951057      0.809017  -0.587785      0.309017

basis_inverse = basis_matrix\I

## 5×5 Array{Float64,2}:
##  0.2        0.2        0.2           0.2        0.2     
##  0.380423   0.235114   5.55112e-17  -0.235114  -0.380423
##  0.323607  -0.123607  -0.4          -0.123607   0.323607
##  0.235114  -0.380423  -8.32667e-17   0.380423  -0.235114
##  0.123607  -0.323607   0.4          -0.323607   0.123607

Step 4b: Evaluate the policy function

To loop and maximize the Bellman at each grid point we need a function eval_policy_function(coefficients, capital, params_ti) that lets us evaluate the continuation value given a vector of coefficients coefficients, a vector of capital nodes capital, and the parameters params_ti to scale capital back into $[-1,1]$

It needs to:

1. Scale capital back into $[-1,1]$
2. Use the coefficients and Chebyshev polynomials to evaluate the value function

Step 4b: Evaluate the continuation value

Here's a simple way to do it:

shrink_grid(capital) = 2*(capital - params_ti.capital_lower)/(params_ti.capital_upper - params_ti.capital_lower) - 1;
eval_policy_function(coefficients, capital, params_ti) = 
    coefficients' * [cheb_polys.(shrink_grid(capital), n) for n = 0:params_ti.num_points - 1];

The top function inherits params_ti from the bottom function

Step 5: Loop

Construct a function that loops over the grid points and solves the Euler given $\Psi(x;b^{(p)})$

function loop_grid_ti(params_ti, basis_inverse, basis_matrix, grid, capital_grid, coefficients)
    consumption = similar(coefficients)
    for (iteration, capital) in enumerate(capital_grid)
        function consumption_euler(consumption_guess)
            capital_next = capital^params_ti.alpha - consumption_guess
            # Next period consumption based on policy approximant
            consumption_next = eval_policy_function(coefficients, capital_next, params_ti)
            consumption_next = max(1e-10, consumption_next)
            # Organize Euler so it's g(c,k) = 0
            euler_error = consumption_guess^(-params_ti.eta) /
                (params_ti.beta*consumption_next^(-params_ti.eta)*params_ti.alpha*(capital_next)^(params_ti.alpha - 1)) - 1
            return euler_error
        end
        # Search over consumption such that Euler = 0
        consumption[iteration] = fzero(consumption_euler, 0., capital)
    end
    return consumption
end

## loop_grid_ti (generic function with 1 method)

Step 5: Loop

function solve_ti(params_ti, basis_inverse, basis_matrix, grid, capital_grid, coefficients)
    error = 1e10
    iteration = 1
    consumption = similar(coefficients)
    consumption_prev, coefficients_prev = similar(coefficients), similar(coefficients)
    coefficients_store = Vector{Vector}(undef, 1)
    coefficients_store[1] = coefficients
    while error > params_ti.tolerance
        consumption = loop_grid_ti(params_ti, basis_inverse, basis_matrix, grid, capital_grid, coefficients)
        if iteration > 1
            coefficients = params_ti.damp*(basis_inverse*consumption) + (1 - params_ti.damp)*coefficients_prev
        else
            coefficients = basis_inverse*consumption
        end
        error = maximum(abs.((consumption - consumption_prev)./(consumption_prev)))
        consumption_prev, coefficients_prev = deepcopy(consumption), deepcopy(coefficients)
        if mod(iteration, 5) == 0
            println("Maximum Error of $(error) on iteration $(iteration).")
            append!(coefficients_store, [coefficients])
        end
        iteration += 1
    end
    return coefficients, consumption, coefficients_store
end

## solve_ti (generic function with 1 method)

Step 5: Loop

solution_coeffs, consumption, intermediate_coefficients =
    solve_ti(params_ti, basis_inverse, basis_matrix, grid, capital_grid, coefficients)

## Maximum Error of 0.06899533243252934 on iteration 5.
## Maximum Error of 1.481484067728511 on iteration 10.
## Maximum Error of 0.7160885861973345 on iteration 15.
## Maximum Error of 0.3030828724582769 on iteration 20.
## Maximum Error of 0.17473143140005473 on iteration 25.
## Maximum Error of 505.4805245124472 on iteration 30.
## Maximum Error of 1.193558398749873 on iteration 35.
## Maximum Error of 0.28476702156043676 on iteration 40.
## Maximum Error of 0.19536671877310902 on iteration 45.
## Maximum Error of 235825.18031971727 on iteration 50.
## Maximum Error of 1.4632905297050023 on iteration 55.
## Maximum Error of 0.3819128797728578 on iteration 60.
## Maximum Error of 0.2069777248265044 on iteration 65.
## Maximum Error of 0.2882977765088428 on iteration 70.
## Maximum Error of 1.5557188732619867 on iteration 75.
## Maximum Error of 0.20172230910005298 on iteration 80.
## Maximum Error of 0.06127838612896111 on iteration 85.
## Maximum Error of 0.028153777682640414 on iteration 90.
## Maximum Error of 0.014466543089840038 on iteration 95.
## Maximum Error of 0.006878960963739522 on iteration 100.
## Maximum Error of 0.003130222595497513 on iteration 105.
## Maximum Error of 0.001391461005386054 on iteration 110.
## Maximum Error of 0.0006111544900213002 on iteration 115.
## Maximum Error of 0.00026681831311659726 on iteration 120.
## Maximum Error of 0.00011614043376537828 on iteration 125.
## Maximum Error of 5.047881962395737e-5 on iteration 130.
## Maximum Error of 2.1923791129685817e-5 on iteration 135.
## Maximum Error of 9.518360650356241e-6 on iteration 140.

## ([0.10216025688175603, 0.030492240762756, -0.0018614730785602319, 0.0002412050526698486, -3.673940758425717e-5], [0.12978513558181665, 0.120459207876453, 0.10398539694181286, 0.08507201135884981, 0.07150154187804965], Array{T,1} where T[[0.0, 0.0, 0.0, 0.0, 0.0], [1.276313022082098e-10, 7.183782173526791e-12, -5.1738097180129495e-12, -2.5864997905351185e-13, -8.831154896937806e-14], [5.317142136502967e-10, 5.527847639739073e-10, 1.6703531496499355e-10, 1.8141242726804023e-11, -1.4478161076640903e-13], [8.316966371969167e-9, 5.18753915343869e-9, 2.0308553662549534e-10, -1.1491350952729583e-11, 1.7095210124630882e-12], [3.386275145017677e-8, 1.2764005373987535e-8, -3.0912462874207683e-10, 4.0250155043961724e-11, -1.1414759632199032e-11], [7.52422466539421e-8, 1.057016474990513e-8, -4.436603667835434e-9, -1.1744349926279e-10, -4.947772660433226e-11], [3.5814952980292754e-8, -1.4472521472050847e-8, 1.8347071102204086e-8, 5.2283868000856916e-9, -1.5042279183557698e-10], [8.649173370341168e-7, 6.23976511971795e-7, 2.580472512479873e-8, -5.8530576443756555e-9, 3.6058883046539203e-10], [3.5573111138957096e-6, 1.2734335635389533e-6, -3.0722754197916816e-8, 1.0899793501970388e-8, -1.292670762307699e-9], [8.394998791841883e-6, 1.5754988650395768e-6, -4.3452939029739223e-7, -2.0464224738503144e-8, -5.013974305135995e-9]  …  [0.0937722121997489, 0.027049729315754424, -0.0017387008487714786, 0.00022715074679682435, -3.8481841856618996e-5], [0.09828586835137634, 0.028819431842777576, -0.0018084381167883102, 0.0002370462002013028, -3.7730003683528356e-5], [0.10043008020063873, 0.029723588228166724, -0.0018394468427340885, 0.00023995653980080382, -3.7218854408817556e-5], [0.10140171061061844, 0.030149857542289698, -0.001852360245859877, 0.00024079574250233525, -3.6957483395931e-5], [0.10183187323362813, 0.030342720870865328, -0.0018576775097620402, 0.00024105811168629785, -3.683563041501696e-5], [0.1020202511316293, 0.0304281893808404, -0.0018598919520605756, 0.0002411489742540459, -3.678086910692431e-5], [0.10210233687095274, 0.03046567454426408, -0.0018608274951932573, 0.00024118324939464226, -3.6756664109979044e-5], [0.10213802583793859, 0.030482029928345013, -0.0018612270529681154, 0.0002411969602829308, -3.674605381205512e-5], [0.10215352705674494, 0.03048914745670296, -0.0018613988826236723, 0.00024120264584137242, -3.674142338781434e-5], [0.10216025688175603, 0.030492240762756, -0.0018614730785602319, 0.0002412050526698486, -3.673940758425717e-5]])

Now lets plot our solutions

capital_levels = range(params_ti.capital_lower, params_ti.capital_upper, length = 100);
eval_points = shrink_grid.(capital_levels);
solution = similar(intermediate_coefficients);
for (iteration, coeffs) in enumerate(intermediate_coefficients)
    solution[iteration] = [coeffs' * [cheb_polys.(capital, n) for n = 0:params_ti.num_points - 1] for capital in eval_points];
end

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the consumption policy function

Plot the final consumption policy function

Now lets try simulating

function simulate_model(params_ti, solution_coeffs, time_horizon = 100)
    capital_store = zeros(time_horizon + 1)
    consumption_store = zeros(time_horizon)
    capital_store[1] = params_ti.k_0
    for t = 1:time_horizon
        capital = capital_store[t]
        consumption_store[t] = eval_policy_function(solution_coeffs, capital, params_ti)
        capital_store[t+1] = capital^params_ti.alpha - consumption_store[t]
    end
    return consumption_store, capital_store
end;

Now lets try simulating

time_horizon = 100;
consumption, capital = simulate_model(params_ti, solution_coeffs, time_horizon);
plot(1:time_horizon, consumption, color = :red, linewidth = 4.0, tickfontsize = 14, guidefontsize = 14, label = "Consumption", legend = :right, grid = false, size = (500, 300));
plot!(1:time_horizon, capital[1:end-1], color = :blue, linewidth = 4.0, linestyle = :dash, label = "Capital");
plot!(1:time_horizon, params_ti.steady_state*ones(time_horizon), color = :purple, linewidth = 2.0, linestyle = :dot, label = "Steady State Capital")

A short overview of discretization + VFI

When we use discretization methods we create a grid on our state space, typically evenly spaced

This becomes our actual state space, not just collocation points

How does it work?

A short overview of discretization + VFI

The discretized state space implies a discretized control space

If there are only a finite number of states tomorrow conditional on the current state, then there is only a finite number of valid controls

This makes solving easy!

Search over all possible controls today until you find the one that yields the highest value of the RHS of the Bellman: just requires looping and a max operator

The maximized value is the new value of this discretized state

3 loops now: outer VFI loop, middle capital grid loop, inner consumption loop

Discretizing the state space

using LinearAlgebra
using Optim
using Plots
params_dis = (alpha = 0.75, beta = 0.95, eta = 2,
                steady_state = (0.75*0.95)^(1/(1 - 0.75)), k_0 = (0.75*0.95)^(1/(1 - 0.75))*.75,
                capital_upper = (0.75*0.95)^(1/(1 - 0.75))*1.5, capital_lower = (0.75*0.95)^(1/(1 - 0.75))/2,
                tolerance = 0.0001, max_iterations = 1000)

## (alpha = 0.75, beta = 0.95, eta = 2, steady_state = 0.25771486816406236, k_0 = 0.19328615112304676, capital_upper = 0.3865723022460935, capital_lower = 0.12885743408203118, tolerance = 0.0001, max_iterations = 1000)

Discretizing the state space

function iterate_value(grid, params) 
    grid_size = size(grid, 1)
    V, V_prev = zeros(grid_size, 1), zeros(grid_size, 1)
    V_store = Array{Float64}(undef, grid_size, params.max_iterations)
    max_diff = 1e10
    it = 1
    while max_diff > params.tolerance && it <= params.max_iterations # iterate on the value function
        for (iteration, grid_point) in enumerate(grid) # iterate across the capital grid
            # possible consumption values (output + remaining capital - capital next period)
            c_vec = grid_point.^params.alpha .- grid
            value_max = -Inf
            # find consumption that maximizes the right hand side of the Bellman, search over ones with positive consumption
            for (it_inner, consumption) in enumerate(c_vec[c_vec .> 0]) # iterate across possible consumption values
                value_temp = consumption^(1 - params.eta)/(1 - params.eta) + params.beta*V[it_inner]
                value_max = max(value_temp, value_max)
            end
            V[iteration] = value_max
        end
        max_diff = maximum(abs.(V .- V_prev))
        if mod(it,10) == 0
            println("Current maximum value difference at iteration $it is $max_diff.")
        end
        V_prev = copy(V)
        V_store[:,it] = V
        if it == params.max_iterations
            println("Hit maximum iterations")
            break
        end
        it += 1
    end
    V_store = V_store[:, 1:it-1]
    return V, V_store
end

## iterate_value (generic function with 1 method)

Discretizing the state space

        max_diff = maximum(abs.((V .- V_prev)./V_prev))
        if mod(it,10) == 0
            println("Current maximum value difference at iteration $it is $max_diff.")
        end
        V_prev = copy(V)
        V_store[:,it] = V
        if it == params.max_iterations
            println("Hit maximum iterations")
            break
        end
        it += 1
    end
    V_store = V_store[:, 1:it-1]
    return V, V_store
end

Discretizing the state space

grid_size = 3;
grid = collect(range(params_dis.capital_lower,
    stop = params_dis.capital_upper,
    length = grid_size))

## 3-element Array{Float64,1}:
##  0.12885743408203118
##  0.25771486816406236
##  0.3865723022460935

value, v_store = @time iterate_value(grid, params_dis)

## Current maximum value difference at iteration 10 is 7.310316889342374.
## Current maximum value difference at iteration 20 is 4.376956759187493.
## Current maximum value difference at iteration 30 is 2.6206456931746516.
## Current maximum value difference at iteration 40 is 1.5690773811596443.
## Current maximum value difference at iteration 50 is 0.9394645886236788.
## Current maximum value difference at iteration 60 is 0.5624921523154001.
## Current maximum value difference at iteration 70 is 0.33678482962292833.
## Current maximum value difference at iteration 80 is 0.20164551807033604.
## Current maximum value difference at iteration 90 is 0.12073262030057208.
## Current maximum value difference at iteration 100 is 0.07228707954499214.
## Current maximum value difference at iteration 110 is 0.043280944753263384.
## Current maximum value difference at iteration 120 is 0.0259139003889004.
## Current maximum value difference at iteration 130 is 0.015515609402569908.
## Current maximum value difference at iteration 140 is 0.009289768484109118.
## Current maximum value difference at iteration 150 is 0.005562127548415674.
## Current maximum value difference at iteration 160 is 0.0033302512239856696.
## Current maximum value difference at iteration 170 is 0.0019939444247540905.
## Current maximum value difference at iteration 180 is 0.0011938481818845048.
## Current maximum value difference at iteration 190 is 0.0007148010063247057.
## Current maximum value difference at iteration 200 is 0.00042797776669090126.
## Current maximum value difference at iteration 210 is 0.0002562460981039294.
## Current maximum value difference at iteration 220 is 0.00015342400445206295.
##   0.150962 seconds (455.00 k allocations: 25.278 MiB, 6.89% gc time)

## ([-231.9798759489783; -192.32427374317618; -187.00837177812517], [-11.599085658067757 -22.618217033232128 … -231.97977925359004 -231.9798759489783; -9.616289844742465 -18.751765197247806 … -192.32419357729867 -192.32427374317618; -9.644703799156582 -18.807172408355335 … -187.00829562054153 -187.00837177812517])

The value function: every 20 iterations

The value function: final

Discretizing the state space

grid_size = 100;
grid = collect(range(params_dis.capital_lower,
    stop = params_dis.capital_upper,
    length = grid_size));
value, v_store = @time iterate_value(grid, params_dis)

## Current maximum value difference at iteration 10 is 6.914720355073825.
## Current maximum value difference at iteration 20 is 3.8092197025250414.
## Current maximum value difference at iteration 30 is 2.221007019891772.
## Current maximum value difference at iteration 40 is 1.3164056274758593.
## Current maximum value difference at iteration 50 is 0.7840556792955624.
## Current maximum value difference at iteration 60 is 0.4679885486935973.
## Current maximum value difference at iteration 70 is 0.27994507576624983.
## Current maximum value difference at iteration 80 is 0.16751908240780722.
## Current maximum value difference at iteration 90 is 0.1002998626648548.
## Current maximum value difference at iteration 100 is 0.05997003214901042.
## Current maximum value difference at iteration 110 is 0.03590627349490205.
## Current maximum value difference at iteration 120 is 0.0214984122918338.
## Current maximum value difference at iteration 130 is 0.01287189357412899.
## Current maximum value difference at iteration 140 is 0.007706878160803399.
## Current maximum value difference at iteration 150 is 0.004614392641116183.
## Current maximum value difference at iteration 160 is 0.0027628073264054365.
## Current maximum value difference at iteration 170 is 0.0016541948023416353.
## Current maximum value difference at iteration 180 is 0.0009904275328835865.
## Current maximum value difference at iteration 190 is 0.0005930055496037312.
## Current maximum value difference at iteration 200 is 0.00035505432774129986.
## Current maximum value difference at iteration 210 is 0.00021258414147951044.
## Current maximum value difference at iteration 220 is 0.0001272819782229817.
##   0.151063 seconds (159.38 k allocations: 132.604 MiB, 19.39% gc time)

## ([-212.42908333245703; -211.7639543138398; … ; -183.072370381581; -182.92812836469787], [-11.599085658067757 -22.618217033232128 … -212.42898484408866 -212.42908333245703; -11.512646026583107 -22.44965975183706 … -211.76385582547144 -211.7639543138398; … ; -9.63308569160569 -18.784517098631092 … -183.07232481286312 -183.072370381581; -9.644703799156582 -18.807172408355335 … -182.92808279597998 -182.92812836469787])

The value function: every 20 iterations

The value function: final