Lecture 7

.title[
# Lecture 7
]
.subtitle[
## Solution methods for discrete time dynamic models
]
.author[
### Ivan Rudik
]
.date[
### AEM 7130
]

---

```r
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  xaringanExtra, JuliaCall
)
options(htmltools.dir.version = FALSE)
library(knitr)
opts_chunk$set(
  prompt = T, ## See hook below. I basically want a "$" prompt for every bash command in this lecture.
  fig.align = "center", fig.width=10, fig.height=6, 
  out.width="748px", out.length="520.75px",
  dpi = 300, #fig.path='Figs/',
  cache = F#, echo=F, warning=F, message=F
  )
## Next hook based on this SO answer: https://stackoverflow.com/a/39025054
knit_hooks$set(
  prompt = function(before, options, envir) {
    options(
      prompt = if (options$engine %in% c('sh','bash')) '$ ' else ' ',
      continue = if (options$engine %in% c('sh','bash')) '$ ' else ' '
      )
})
julia_setup(JULIA_HOME = "/Applications/Julia-1.8.app/Contents/Resources/julia/bin")
```

---

# Roadmap

1. Intuition for solving dynamic 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

---

# Things to do

3. 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

---

# How do we solve dynamic models?

---

# 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

---

# Solutions to economic models

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
  - Bellman
  - Euler
2. Fixed points

---

# Our general example

Consider the following problem we will be using for all the 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_t) = \max_{c_t}  u(c_t) + \beta V(k_{t+1}) \notag \\
\text{subject to:} \,\,\,\,\, k_{t+1}' = f(k_t) - c_t \notag 
\end{gather}`

We can then express the value function in terms of itself, the current state, and the current consumption choice:
`$$V(k_t) = \max_{c_t}  u(c_t) + \beta V(f(k_t) - c_t)$$`

---

# Our general example

`$$V(k_t) = \max_{c_t}  u(c_t) + \beta V(f(k_t) - c_t)$$`

How do we solve this?

Main idea:

1. Guess `$V(k_t)$`
2. Given guess, do the maximization on the right hand side at some set of states `$\mathbf{k^i_t}$`
3. Maximized right hand side gives us new values of `$V(k_t)$`
4. Repeat

---

# Our general example

Another equilibrium condition is the .hi[Euler equation]

For our problem it is
`$$u'(c_t) = \beta u'(c_{t+1}) f'(k_{t+1})$$`

Plug in the policy function `$c_t = C(k_t)$`:

`$$u'(C(k_t)) = \beta u'(C(k_{t+1})) f'(k_{t+1})$$`

---

# Our general example

Recognize `$k_{t+1}$` is a function of the current policy and state `$(C(k_t),k_t)$`:

`$$k_{t+1} = f(k_t) - C(k_t)$$`

Use this to express the Euler equation in terms of `$k_t$` and `$C$`

`$$C(k_t) = u'^{(-1)}\left[\beta u'\left[C(k_{t+1}(C(k_t),k_t))\right] f'\left[k_{t+1}(C(k_t),k_t)\right]\right]$$`

---

# Our general example

`$$C(k_t) = u'^{(-1)}\left[\beta u'\left[C(k_{t+1}(C(k_t),k_t))\right] f'\left[k_{t+1}(C(k_t),k_t)\right]\right]$$`

How do we solve this?

Main idea:

1. Guess `$C(k_t)$`
2. Given guess, evaluate the right hand side at some set of states `$\mathbf{k^i_t}$`
3. Evaluated right hand side gives us new values of `$C(k_t)$`
4. Repeat

---

# Value function iteration

---

# Method 1: Value function iteration

In VFI we approximate the .hi[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

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

.hi[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

---

# Method 1: Value function iteration

.hi[Step 5:] While convergence criterion `$>$` tolerance **(outer loop [fixed point])**
+ Start iteration `$p$`
+ For each grid point **(inner loop [right hand side maximization])**
    - Maximize the right hand side of the Bellman equation at each grid point using the approximating value function `$\Gamma(k_{t+1};b^{(p)})$` in place of `$V(k_{t+1})$`
    - Recover the maximized values `$V^{(p)}$` at each grid point, conditional on `$\Gamma(k_{t+1};b^{(p)})$`

...

---

# Method 1: Value function iteration

.hi[Step 5:] While convergence criterion `$>$` tolerance **(outer loop, continued)**
+ Fit the polynomial to the maximized values `$V^{(p)}$` 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)

---

# Functions we will code up

We will need to code up six key functions for all of the algorithms

- `cheb_nodes(n)`: construct degree `$n$` collocation grid 
- `cheb_polys(x, n)`: evaluate degree `$n$` Chebyshev polynomials
- `construct_basis_matrix(grid, params)`: construct full `$n\times n$` Chebyshev basis matrix
- `eval_approx_function(coefficients, grid, params)`: evaluate approximating function at grid points `grid`
- `loop_grid(params, capital_grid, coefficients)`: loop over the collocation grid
- `solve_algorithm(params, basis_inverse, capital_grid, coefficients)`: iterate on the fixed point

---

# 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: `$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)}$`

---

# Step 1: Select the number of points and domain

Put everything in a **named tuple** to make passing things easier

```julia
 using LinearAlgebra
 using Optim
 using Plots
 params = (alpha = 0.75, # capital share
           beta = 0.95, # discount
           eta = 2, # EMUC
           steady_state = (0.75*0.95)^(1/(1 - 0.75)),
           k_0 = (0.75*0.95)^(1/(1 - 0.75))/2, # initial state
           capital_upper = (0.75*0.95)^(1/(1 - 0.75))*1.01, # upper bound
           capital_lower = (0.75*0.95)^(1/(1 - 0.75))/2, # lower bound
           num_points = 7, # number of grid points
           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

```julia
 coefficients = zeros(params.num_points) # # coeffs = # grid points in collocation
```

```
## 7-element Vector{Float64}:
##  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?

.hi[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` from last lecture constructs the grid on `$[-1,1]$`

Recall:

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

---

# Step 4: Construct the grid and basis matrix

```julia
 cheb_nodes(n) = cos.(pi * (2*(1:n) .- 1)./(2n));
 grid = cheb_nodes(params.num_points) # [-1, 1] grid with n points
```

```
## 7-element Vector{Float64}:
##   0.9749279121818236
##   0.7818314824680298
##   0.4338837391175582
##   6.123233995736766e-17
##  -0.43388373911755806
##  -0.7818314824680297
##  -0.9749279121818236
```

Our actual capital domain isn't on `$[-1,1]$`, we need to expand the grid to some arbitrary `$[a,b]$` using a function `expand_grid(grid, params)`

---

# Step 4: Construct the grid and basis matrix

```julia
 expand_grid(grid, params) = # function that expands [-1,1] to [a,b]
    (1 .+ grid)*(params.capital_upper - params.capital_lower)/2 .+ params.capital_lower
```

```
## expand_grid (generic function with 1 method)
```

```julia
 capital_grid = expand_grid(grid, params)
```

```
## 7-element Vector{Float64}:
##  0.2586443471450049
##  0.2459545728087113
##  0.2230883895732961
##  0.19457472546386706
##  0.16606106135443804
##  0.14319487811902284
##  0.13050510378272925
```

---

# Step 4: Construct the grid and basis matrix

Make the inverse function to shrink from capital to Chebyshev space `shrink_grid(capital)`

```julia
 shrink_grid(capital) = 
   2*(capital - params.capital_lower)/(params.capital_upper - params.capital_lower) - 1;
 shrink_grid.(capital_grid)
```

```
## 7-element Vector{Float64}:
##   0.9749279121818237
##   0.7818314824680297
##   0.43388373911755806
##  -2.220446049250313e-16
##  -0.43388373911755806
##  -0.7818314824680297
##  -0.9749279121818236
```

`shrink_grid` will inherit `params` from wrapper functions

---

# Step 4: Construct the grid and basis matrix

`cheb_polys(x, n)` from last lecture gives us the nth degree Chebyshev polynomial at point x

```julia
 # 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 4: Construct the grid and basis matrix

`cheb_polys.(grid, n)` gives us the nth degree Chebyshev polynomial at all points on our grid

```julia
 cheb_polys.(grid, 2) # 2nd degree Cheb poly at each grid point
```

```
## 7-element Vector{Float64}:
##   0.9009688679024193
##   0.22252093395631434
##  -0.6234898018587334
##  -1.0
##  -0.6234898018587336
##   0.22252093395631412
##   0.9009688679024193
```

---

# Step 4: Construct the grid and basis matrix

In our basis matrix, rows are grid points, columns are basis functions, make a function `construct_basis_matrix(grid, params)` that makes the basis matrix for some arbitrary grid of points

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

```
## 7×7 Matrix{Float64}:
##  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  -2.22045e-16  -1.0           1.0       -1.11022e-15  -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
```

---

# Step 4: Pre-invert your basis matrix

.hi[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)

```julia
 basis_inverse = basis_matrix \ I # pre-invert
```

```
## 7×7 Matrix{Float64}:
##  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
```

---

# Pre-Step 5: Evaluate the continuation value

To maximize the Bellman at each grid point we need to evaluate the value function

We need to make 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 model parameters `params`

It needs to:

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

---

# Pre-Step 5: Evaluate the continuation value

```julia
 # evaluates V on the [-1,1]-equivalent grid
 eval_value_function(coefficients, grid, params) = construct_basis_matrix(grid, params) * coefficients;
```

---

# Step 5: Inner loop over grid points

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

Pseudocode:
```
for each grid point i:
define the Bellman as a closure so it can take in parameters
maximize the Bellman by choosing consumption
store maximized value in a vector v[i]
end

return vector of maximized values v
```

---

# Step 5: Inner loop over grid points

```julia
 function loop_grid(params, 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)[1] # 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

Construct a function `solve_vfi(params, basis_inverse, capital_grid, coefficients)` that iterates on `loop_grid` and solves for the coefficient vector `$b$` until the maximized values on the grid converge

Pseudocode:
```
while error > tolerance
call loop_grid to get maximized values
use maximized values and basis matrix to get new coefficients
error is maximum relative difference between current and previous maximized values
end

return vector of maximized values v
```

---

# Step 5: Outer loop iterating on Bellman

```julia
 function solve_vfi(params, basis_inverse, 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, 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

```julia
 solution_coeffs, max_value, intermediate_coefficients =
     solve_vfi(params, basis_inverse, capital_grid, coefficients)
```

```
## Maximum Error of 0.3301919884226089 on iteration 5.
## Maximum Error of 0.10801399197451178 on iteration 10.
## Maximum Error of 0.05647917855011511 on iteration 15.
## Maximum Error of 0.034833389245083446 on iteration 20.
## Maximum Error of 0.023286433761111308 on iteration 25.
## Maximum Error of 0.016301543092581618 on iteration 30.
## Maximum Error of 0.011747480470413624 on iteration 35.
## Maximum Error of 0.008631245645920323 on iteration 40.
## Maximum Error of 0.006427690604127157 on iteration 45.
## Maximum Error of 0.004833073684245354 on iteration 50.
## Maximum Error of 0.003659714890062249 on iteration 55.
## Maximum Error of 0.0027856923769765014 on iteration 60.
## Maximum Error of 0.0021286867102128133 on iteration 65.
## Maximum Error of 0.0016314249677370307 on iteration 70.
## Maximum Error of 0.0012531160571389703 on iteration 75.
## Maximum Error of 0.0009641708791272537 on iteration 80.
## Maximum Error of 0.0007428166750767924 on iteration 85.
## Maximum Error of 0.000572852149880767 on iteration 90.
## Maximum Error of 0.0004421161961630204 on iteration 95.
## Maximum Error of 0.0003414181452772865 on iteration 100.
## Maximum Error of 0.0002637753971325223 on iteration 105.
## Maximum Error of 0.00020386108209659878 on iteration 110.
## Maximum Error of 0.00015759845903308656 on iteration 115.
## Maximum Error of 0.00012185979301930555 on iteration 120.
```

```
## ([-200.6291758538633, 9.991472391067827, -1.22789926411501, 0.17379460460100873, -0.02621191019442686, 0.00395400691320583, -0.0007409750421374391], [-191.87342361439286, -193.14594489252323, -195.68963652528447, -199.42674752490058, -204.02721963808625, -208.61071718644237, -211.63054159541332], Vector[[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [-36.858949010607404, 6.259861237095008, -0.8221318972860758, 0.11882022830523278, -0.01784945487785583, 0.00269203315390105, -0.0003752023513341855], [-74.87645132736962, 8.931829874526924, -1.1238227131298495, 0.16080476492855306, -0.024292798807138135, 0.003902732334813939, -0.0006447123917418033], [-103.80959945589599, 9.68901311684356, -1.1994026067805432, 0.17056095546632696, -0.02554814054980478, 0.0039694305835632155, -0.0007161344250393939], [-125.9200308753979, 9.904403319376973, -1.2195520590229452, 0.17300404781585907, -0.026010058299272096, 0.003959590376948559, -0.0007341995822562963], [-142.93974657626163, 9.9665708864539, -1.2253830534940031, 0.17357716260519426, -0.02615385493582288, 0.003955453304668797, -0.0007391283603728321], [-156.08276002781162, 9.98445374741454, -1.2271575579533587, 0.1737309439806476, -0.026195277719921334, 0.003954320282277471, -0.000740461831093684], [-166.24486078281993, 9.989518053393983, -1.2276872050893808, 0.17377614420429102, -0.02620719354052565, 0.003954072102106873, -0.0007408314670065023], [-174.10590746285402, 9.9909321107981, -1.227839864872191, 0.17378937493063698, -0.02621059015233894, 0.003954021287676794, -0.0007409350663560629], [-180.18802216175337, 9.991323564318826, -1.227882803999439, 0.17379314659915696, -0.02621154420094639, 0.003954010369678741, -0.0007409639787900078]  …  [-196.52576254532752, 9.991472327014671, -1.2278992569991485, 0.1737946039676496, -0.026211910036028768, 0.003954006914552045, -0.0007409750373485068], [-197.53584414896338, 9.991472373473814, -1.2278992621603635, 0.1737946044269843, -0.026211910150932702, 0.003954006913575448, -0.0007409750408210317], [-198.31742601900214, 9.99147238623527, -1.227899263578116, 0.17379460455317255, -0.02621191018249914, 0.003954006913323826, -0.0007409750417566008], [-198.92219916559227, 9.99147238974058, -1.2278992639675093, 0.17379460458788432, -0.026211910191119247, 0.003954006913259123, -0.000740975042017718], [-199.3901610964136, 9.991472390703427, -1.2278992640744515, 0.1737946045974085, -0.026211910193502778, 0.003954006913246129, -0.0007409750420960553], [-199.7522611175407, 9.991472390967925, -1.2278992641038726, 0.17379460460000581, -0.026211910194188358, 0.003954006913236502, -0.0007409750421330943], [-200.03244721124233, 9.99147239104054, -1.2278992641119348, 0.17379460460071294, -0.026211910194353222, 0.003954006913220467, -0.0007409750420998224], [-200.24924986946974, 9.991472391060487, -1.2278992641141895, 0.17379460460092855, -0.026211910194371166, 0.003954006913207492, -0.0007409750421181931], [-200.4170076335967, 9.991472391065942, -1.2278992641147288, 0.17379460460096674, -0.02621191019440478, 0.003954006913239182, -0.0007409750421174219], [-200.54681539359345, 9.991472391067477, -1.2278992641149715, 0.1737946046009521, -0.02621191019442265, 0.003954006913247727, -0.0007409750421140032]])
```

---

# Now lets plot our solutions

```julia
 
 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

```julia
 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

```julia
 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

---

# Fixed point iteration

---

# 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

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

The drawback is that it will not always converge and may be unstable

This can be solved by .hi[damping]

---

# Eq condition: Euler equation

Standard procedure is to iterate on the Euler equation

`$$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

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

.hi[Step 2:] Select an initial vector of coefficients `$b_0$` with the same number of elements as the collocation grid

---

# Method 2: Fixed point iteration

.hi[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)

Notice: we did not have to perform a maximization step .hi-red[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

```julia
 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)
```

```julia
 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$`

```julia
 coefficients = zeros(params_fpi.num_points)
```

```
## 5-element Vector{Float64}:
##  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` from last lecture constructs the grid on `$[-1,1]$`

Recall:

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

---

# Step 4: Construct the grid and basis matrix

```julia
 cheb_nodes(n) = cos.(pi * (2*(1:n) .- 1)./(2n));
 grid = cheb_nodes(params_fpi.num_points) # [-1, 1] grid with n points
```

```
## 5-element Vector{Float64}:
##   0.9510565162951535
##   0.5877852522924731
##   6.123233995736766e-17
##  -0.587785252292473
##  -0.9510565162951535
```

Our actual capital domain isn't on `$[-1,1]$`, we need to expand the grid to some arbitrary `$[a,b]$`

---

# Step 4: Construct the grid and basis matrix

```julia
 expand_grid(grid, params_fpi) = # function that expands [-1,1] to [a,b]
    (1 .+ grid)*(params_fpi.capital_upper - params_fpi.capital_lower)/2 .+ params_fpi.capital_lower
```

```
## expand_grid (generic function with 1 method)
```

```julia
 capital_grid = expand_grid(grid, params_fpi)
```

```
## 5-element Vector{Float64}:
##  0.3802655705208513
##  0.33345536756572974
##  0.2577148681640623
##  0.18197436876239492
##  0.1351641658072734
```

---

# Step 4: Construct the grid and basis matrix

Make the inverse function to shrink from capital to Chebyshev space `shrink_grid(capital, params)`

```julia
 shrink_grid(capital) = 
   2*(capital - params_fpi.capital_lower)/(params_fpi.capital_upper - params_fpi.capital_lower) - 1;
 shrink_grid.(capital_grid)
```

```
## 5-element Vector{Float64}:
##   0.9510565162951536
##   0.5877852522924731
##  -2.220446049250313e-16
##  -0.5877852522924731
##  -0.9510565162951536
```

`shrink_grid` will inherit `params_fpi` from wrapper functions

---

# Step 4: Construct the grid and basis matrix

`cheb_polys(x, n)` from last lecture gives us the nth degree Chebyshev polynomial at point x

---

# Step 4: Construct the grid and basis matrix

`cheb_polys.(grid, n)` gives us the nth degree Chebyshev polynomial at all points on our grid

```julia
 cheb_polys.(grid, 2) # 2nd degree Cheb poly at each grid point
```

```
## 5-element Vector{Float64}:
##   0.8090169943749472
##  -0.30901699437494745
##  -1.0
##  -0.3090169943749477
##   0.8090169943749472
```

---

# Step 4: Construct the grid and basis matrix

In our basis matrix, rows are grid points, columns are basis functions, make a function `construct_basis_matrix(grid, params)` that makes the basis matrix for some arbitrary grid of points

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

```
## 5×5 Matrix{Float64}:
##  1.0   0.951057      0.809017   0.587785      0.309017
##  1.0   0.587785     -0.309017  -0.951057     -0.809017
##  1.0  -2.22045e-16  -1.0        6.66134e-16   1.0
##  1.0  -0.587785     -0.309017   0.951057     -0.809017
##  1.0  -0.951057      0.809017  -0.587785      0.309017
```

---

# Step 4: Pre-invert your basis matrix

```julia
 basis_inverse = basis_matrix \ I # pre-invert
```

```
## 5×5 Matrix{Float64}:
##  0.2        0.2        0.2           0.2        0.2
##  0.380423   0.235114  -4.69836e-17  -0.235114  -0.380423
##  0.323607  -0.123607  -0.4          -0.123607   0.323607
##  0.235114  -0.380423   2.66269e-17   0.380423  -0.235114
##  0.123607  -0.323607   0.4          -0.323607   0.123607
```

---

# Pre-Step 5: Evaluate the policy function

We need to make 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 model parameters `params_fpi`

It needs to:

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

---

# Pre-Step 5: Evaluate the policy function

```julia
 # evaluates V on the [-1,1]-equivalent grid
 eval_policy_function(coefficients, capital, params_fpi) = 
     construct_basis_matrix(capital, params_fpi) * coefficients;
```

---

# Step 5: Loop

Construct a function `consumption_euler(params_fpi, capital, coefficients)` that evaluates the RHS of the Euler

```julia
 function consumption_euler(params_fpi, capital, coefficients)
     # RHS: Current consumption given current capital
     consumption = eval_policy_function(coefficients, capital, params_fpi)[1]
     # 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)[1]
     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 `loop_grid_fpi(params_fpi, capital_grid, coefficients)` that loops over the grid points and evaluates the RHS of the Euler given `$\Psi(x;b^{(p)})$`

```julia
 function loop_grid_fpi(params_fpi, 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

Construct a function `solve_fpi(params_fpi, basis_inverse, capital_grid, coefficients)` that iterates on `loop_grid_fpi` and solves for the coefficient vector `$b$` until the consumption values on the grid converge

```julia
 function solve_fpi(params_fpi, basis_inverse, 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, 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

```julia
 solution_coeffs, consumption, intermediate_coefficients =
     solve_fpi(params_fpi, basis_inverse, capital_grid, coefficients)
```

```
## Maximum Error of 0.06899533243794273 on iteration 5.
## Maximum Error of 1.4814840914335052 on iteration 10.
## Maximum Error of 0.7160886352973411 on iteration 15.
## Maximum Error of 0.30308290299478047 on iteration 20.
## Maximum Error of 0.17473142803887126 on iteration 25.
## Maximum Error of 505.4941565432374 on iteration 30.
## Maximum Error of 1.1935815404160586 on iteration 35.
## Maximum Error of 0.2847765742094213 on iteration 40.
## Maximum Error of 0.19537333980794477 on iteration 45.
## Maximum Error of 236164.34266168138 on iteration 50.
## Maximum Error of 1.4738022100603876 on iteration 55.
## Maximum Error of 0.38721294651838284 on iteration 60.
## Maximum Error of 0.21074834397654027 on iteration 65.
## Maximum Error of 0.2361368548338292 on iteration 70.
## Maximum Error of 4.722253098129778 on iteration 75.
## Maximum Error of 0.03764118064354074 on iteration 80.
## Maximum Error of 0.007970939078890096 on iteration 85.
## Maximum Error of 0.002573562781229235 on iteration 90.
## Maximum Error of 0.0009352774328684446 on iteration 95.
## Maximum Error of 0.0003237453044014674 on iteration 100.
## Maximum Error of 0.00010972350909872122 on iteration 105.
## Maximum Error of 3.685905275827949e-5 on iteration 110.
## Maximum Error of 1.2336267910941606e-5 on iteration 115.
```

```
## ([0.10216143756493908, 0.030492616130879852, -0.0018615048468132854, 0.00024120846125141784, -3.673979487433395e-5], [0.12978669111540167, 0.12046064710919156, 0.10398662400248394, 0.08507298585015081, 0.07150233223860815], Vector[[0.0, 0.0, 0.0, 0.0, 0.0], [1.2763130221485729e-10, 7.183782174425037e-12, -5.173809717053596e-12, -2.5864997862855907e-13, -8.831154889936558e-14], [5.317142229917236e-10, 5.527847778508575e-10, 1.670353206030028e-10, 1.814124387504668e-11, -1.4478152961250365e-13], [8.316967352693103e-9, 5.187539853019581e-9, 2.03085574859909e-10, -1.1491354530168337e-11, 1.7095212937976074e-12], [3.386275892861658e-8, 1.2764008226593374e-8, -3.0912466477914574e-10, 4.025016138289332e-11, -1.141476294521091e-11], [7.52422641927768e-8, 1.0570161284982647e-8, -4.4366050141360056e-9, -1.1744325496207776e-10, -4.947771112903972e-11], [3.581525952746443e-8, -1.447209127964556e-8, 1.834736552650865e-8, 5.2285236466049244e-9, -1.5039924806733335e-10], [8.649742496411903e-7, 6.240251021601133e-7, 2.5809596063174637e-8, -5.8532653904643626e-9, 3.605947196171803e-10], [3.557675838547508e-6, 1.2735910171892387e-6, -3.0722078466502715e-8, 1.0900582854042943e-8, -1.2927503423060941e-9], [8.396083635910527e-6, 1.5757179126951788e-6, -4.3457907922500314e-7, -2.0463275140970205e-8, -5.014063744733601e-9]  …  [0.014880293352421442, -0.0012599250812935312, -0.0016020724346032572, -7.459264580808092e-5, -2.1168620911705394e-5], [0.06056412990309252, 0.05723653751221666, -0.0027841661603526958, -0.004137596590460779, 0.003880221293456899], [0.09274367455252915, 0.029659083392742768, -0.001611410555591043, 0.0002457512519012199, -3.150220774416544e-5], [0.09883950641692525, 0.0294255474096157, -0.001774882853701155, 0.00023321088847494772, -3.6954845425772695e-5], [0.10100737679223325, 0.030003914490970344, -0.0018414718813008382, 0.00023999644919362001, -3.7123219490068424e-5], [0.10177119243395825, 0.0303097978214859, -0.0018570023721958258, 0.0002411255297881449, -3.690725073674978e-5], [0.10203276998428429, 0.030429834169087837, -0.0018604002640145794, 0.00024122771243407208, -3.679946836519153e-5], [0.10212100280895481, 0.030472550190944813, -0.0018612121452155315, 0.000241220546472494, -3.675909294050314e-5], [0.10215057465291039, 0.03048718729580253, -0.0018614325353545004, 0.00024121237383886463, -3.6745039179331276e-5], [0.1021604610690734, 0.030492126707823614, -0.0018614985832312806, 0.0002412088372761912, -3.674026850711039e-5]])
```

---

# Now lets plot our solutions

```julia
 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

```
## simulate_model (generic function with 2 methods)
```

---

# Now lets try simulating

```julia
 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")
```

---

# Now lets try simulating

---

# Time iteration

---

# 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 .hi-red[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))$$`

for `$j=1,\dots,N$`

---

# Method 3: Time iteration

`$$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))$$`

`$C^{(p)}()$` is our current approximation to the policy function, and we are searching for a .hi-red[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

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

.hi[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

---

# Method 3: Time iteration

.hi[Step 5:] While convergence criterion `$>$` tolerance **(outer loop [fixed point])**
  + Start iteration `$p$`
  + For each grid point **(inner loop [rootfinding])**
    + 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)
  
.hi[Step 6:] Error check your approximation

---

# Step 1: Select the number of points and domain

Put everything in a **named tuple** to make passing things easier

```julia
 using LinearAlgebra, Optim, Plots, 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 = 6, tolerance = 0.00001)
```

---

# 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

```julia
 coefficients = zeros(params_ti.num_points)
```

```
## 6-element Vector{Float64}:
##  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?

.hi[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]$`

Recall:

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

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

```
## cheb_nodes (generic function with 1 method)
```

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

```
## 6-element Vector{Float64}:
##   0.9659258262890683
##   0.7071067811865476
##   0.25881904510252096
##  -0.25881904510252063
##  -0.7071067811865475
##  -0.9659258262890682
```
---

# Step 4: Construct the grid and basis matrix

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

```julia
 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)
```

```julia
 capital_grid = expand_grid(grid, params_ti)
```

```
## 6-element Vector{Float64}:
##  0.3821815916532374
##  0.3488308336097651
##  0.2910656262075347
##  0.22436411012059004
##  0.16659890271835956
##  0.13324814467488724
```

---

# Step 4: Construct the grid and basis matrix

Make the inverse function to shrink from capital to Chebyshev space `shrink_grid(capital)`

```julia
 shrink_grid(capital) = 
   2*(capital - params_ti.capital_lower)/(params_ti.capital_upper - params_ti.capital_lower) - 1;
 shrink_grid.(capital_grid)
```

```
## 6-element Vector{Float64}:
##   0.9659258262890678
##   0.7071067811865472
##   0.25881904510252096
##  -0.2588190451025205
##  -0.7071067811865475
##  -0.9659258262890683
```

`shrink_grid` will inherit `params_ti` from wrapper functions

---

# Step 4: Construct the grid and basis matrix

Use `cheb_polys` to construct the basis matrix

---

# Step 4: Construct the grid and basis matrix

In our basis matrix, rows are grid points, columns are basis functions, make a function `construct_basis_matrix(grid, params)` that makes the basis matrix for some arbitrary grid of points

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

```
## 6×6 Matrix{Float64}:
##  1.0   0.965926   0.866025      0.707107   0.5   0.258819
##  1.0   0.707107  -7.77156e-16  -0.707107  -1.0  -0.707107
##  1.0   0.258819  -0.866025     -0.707107   0.5   0.965926
##  1.0  -0.258819  -0.866025      0.707107   0.5  -0.965926
##  1.0  -0.707107  -2.22045e-16   0.707107  -1.0   0.707107
##  1.0  -0.965926   0.866025     -0.707107   0.5  -0.258819
```

---

# Step 4: Pre-invert your basis matrix

```julia
 basis_inverse = basis_matrix \ I # pre-invert
```

```
## 6×6 Matrix{Float64}:
##  0.166667   0.166667      0.166667   0.166667   0.166667      0.166667
##  0.321975   0.235702      0.086273  -0.086273  -0.235702     -0.321975
##  0.288675  -1.52278e-15  -0.288675  -0.288675   2.62176e-16   0.288675
##  0.235702  -0.235702     -0.235702   0.235702   0.235702     -0.235702
##  0.166667  -0.333333      0.166667   0.166667  -0.333333      0.166667
##  0.086273  -0.235702      0.321975  -0.321975   0.235702     -0.086273
```

---

# Pre-Step 5: Evaluate the policy function

We need to make a function `eval_policy_function(coefficients, capital, params_ti)` that lets us evaluate the policy function given a vector of coefficients `coefficients`, a vector of capital nodes `capital`, and the model parameters `params_ti`

It needs to:

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

---

# Pre-Step 5: Evaluate the policy function

```julia
 # evaluates V on the [-1,1]-equivalent grid
 eval_policy_function(coefficients, capital, params_ti) = 
     construct_basis_matrix(capital, params_ti) * coefficients;
```

---

# Step 5: Loop

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

---

# Step 5: Loop

```julia
 function loop_grid_ti(params_ti, 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 approximating policy function
             consumption_next = eval_policy_function(coefficients, capital_next, params_ti)[1]
             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

Construct a function `solve_ti(params_fpi, basis_inverse, capital_grid, coefficients)` that iterates on `loop_grid_ti` and solves for the coefficient vector `$b$` until the scalar `$c$` values on the grid converge

---

# Step 5: Loop

```julia
 function solve_ti(params_ti, basis_inverse, 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, 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

```julia
 solution_coeffs, consumption, intermediate_coefficients =
     solve_ti(params_ti, basis_inverse, capital_grid, coefficients)
```

```
## Maximum Error of 0.25321486014602274 on iteration 5.
## Maximum Error of 0.5080049627820759 on iteration 10.
## Maximum Error of 0.437979808520151 on iteration 15.
## Maximum Error of 0.3028938891677158 on iteration 20.
## Maximum Error of 0.37071932036765914 on iteration 25.
## Maximum Error of 0.38905423220148744 on iteration 30.
## Maximum Error of 0.3603910536945966 on iteration 35.
## Maximum Error of 0.3351549146984757 on iteration 40.
## Maximum Error of 0.3591346757111148 on iteration 45.
## Maximum Error of 0.3585490158260728 on iteration 50.
## Maximum Error of 0.34330698497556356 on iteration 55.
## Maximum Error of 0.3276713699071389 on iteration 60.
## Maximum Error of 0.29460121133680833 on iteration 65.
## Maximum Error of 0.20990288276886732 on iteration 70.
## Maximum Error of 0.11884185187600733 on iteration 75.
## Maximum Error of 0.060469194796897 on iteration 80.
## Maximum Error of 0.02995709287567328 on iteration 85.
## Maximum Error of 0.015429040567671607 on iteration 90.
## Maximum Error of 0.0074537390569381365 on iteration 95.
## Maximum Error of 0.003416503884818834 on iteration 100.
## Maximum Error of 0.0015221065099616587 on iteration 105.
## Maximum Error of 0.0006686232223713311 on iteration 110.
## Maximum Error of 0.0002917121153220349 on iteration 115.
## Maximum Error of 0.00012685261438477697 on iteration 120.
## Maximum Error of 5.507475726633629e-5 on iteration 125.
## Maximum Error of 2.3892828871781487e-5 on iteration 130.
## Maximum Error of 1.036132860800212e-5 on iteration 135.
```

```
## ([0.10216203236737176, 0.030498670725803596, -0.001857789759114553, 0.00023368783860645652, -3.8642977764246074e-5, 6.661170022615049e-6], [0.13016076813641286, 0.12359702900223878, 0.11148684905718577, 0.0960171043079617, 0.08080506089306848, 0.07090760036622795], Vector[[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [1.4567242688844085e-10, 3.169878115427886e-11, 3.1662133457247602e-12, 1.2113284161164136e-12, 4.4881896990345055e-14, 1.520496566234212e-14], [6.84917675662042e-10, 4.293100807728819e-10, 4.883706342750169e-11, 3.2439271308940757e-12, 1.2813191339134166e-13, 2.2283417796619725e-14], [4.153191672017407e-9, 2.1556300433675913e-9, 8.056510564979214e-11, 5.15973072450304e-12, -4.7855685412401687e-14, 7.654617276973873e-14], [1.6259446055882043e-8, 6.701488218837707e-9, 1.2334234184283768e-10, 5.5300885900205965e-11, 6.408323340941311e-14, 8.858141904396351e-13], [6.269198404363711e-8, 3.0026030536563036e-8, 2.1500002143837277e-9, 3.1674774666817147e-10, 6.42939703591455e-12, 3.60952586191918e-12], [3.104034507901032e-7, 1.631171301935438e-7, 1.0496290192130644e-8, 1.0101053465935438e-9, 2.0769493005477152e-11, 1.0770405404829365e-11], [1.4511219360840434e-6, 6.965645852407328e-7, 2.94282064343977e-8, 4.132068130163973e-9, 3.598395703039279e-11, 5.622111067660242e-11], [6.068482008258223e-6, 2.8224098604480987e-6, 1.3713361071277905e-7, 2.2808122410119338e-8, 2.913990796638209e-10, 2.9672560573567834e-10], [2.662075146568554e-5, 1.3098893818433901e-5, 7.800633216848286e-7, 1.0098427249360026e-7, 1.7440842600487954e-9, 1.1915374548022634e-9]  …  [0.09306282153869429, 0.026853070483444617, -0.0017114664341436594, 0.00021942138015794825, -3.6341032630253375e-5, 6.961272080433971e-6], [0.09793771923786879, 0.028697855510694176, -0.0018021672596051639, 0.00022801152587639448, -3.802377628125088e-5, 6.845623895345248e-6], [0.100271754109348, 0.029664899861432014, -0.0018359977141039117, 0.00023150892114705274, -3.8478684487564766e-5, 6.752804240512059e-6], [0.10133272672202615, 0.030126138248820176, -0.0018489619631233446, 0.0002328347618073334, -3.859634172244187e-5, 6.7033441984452135e-6], [0.10180274792901248, 0.030335700641499843, -0.0018541354555431323, 0.00023334374699115378, -3.862834026773767e-5, 6.679896994251582e-6], [0.10200849176805281, 0.030428664461198504, -0.001856267058332237, 0.00023354669493784306, -3.863792523557689e-5, 6.6692802876339084e-6], [0.10209806674284169, 0.030469425675171274, -0.0018571640749662317, 0.00023363034114399621, -3.864112729565336e-5, 6.664573542458802e-6], [0.10213697123386364, 0.030487196141929306, -0.0018575464296204961, 0.0002336655792271789, -3.864230349249091e-5, 6.662508815736765e-6], [0.10215385030838046, 0.030494921606974985, -0.0018577106263749452, 0.00023368061656785626, -3.864276602943436e-5, 6.661608056503669e-6], [0.1021611699833687, 0.030498275415941793, -0.001857781435494125, 0.00023368707984655295, -3.8642955892818527e-5, 6.66121624338682e-6]])
```

---

# Now lets plot our solutions

```julia
 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

```julia
 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)[1]
         capital_store[t+1] = capital^params_ti.alpha - consumption_store[t]
     end
 
     return consumption_store, capital_store
 end;
```

---

# Now lets try simulating

```julia
 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")
```

---

# Now lets try simulating

---

# Discretization

---

# 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 .hi[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!

---

# A short overview of discretization + VFI

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

```julia
 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

```julia
 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

```julia
         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

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

```
## 3-element Vector{Float64}:
##  0.12885743408203118
##  0.25771486816406236
##  0.3865723022460935
```

```julia
 
 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.087516 seconds (396.82 k allocations: 18.709 MiB, 98.93% compilation 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

```julia
 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.8092197025250982.
## Current maximum value difference at iteration 30 is 2.221007019891772.
## Current maximum value difference at iteration 40 is 1.316405627475831.
## 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.081975 seconds (91.57 k allocations: 38.907 MiB, 26.47% gc time, 2.54% compilation 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