Lecture 9

class: center, middle, inverse, title-slide

.title[
# Lecture 9
]
.subtitle[
## Advanced Methods for Numerical Dynamic Models
]
.author[
### Ivan Rudik
]
.date[
### AEM 7130
]

---

exclude: true

```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. Regression
2. Endogenous grid method
3. Envelope condition method
4. Modified policy iteration

---

# Chebyshev regression

Chebyshev regression works just like normal regression

For a degree `$n$` polynomial approximation, we choose `$m > n+1$` grid points

We then build our basis function matrix `$\Psi$`, but instead of being `$n \times n$` it is `$m \times n$`

Finally we use the standard least-squares equation
`$$c = (\Psi'\Psi)^{-1}\Psi'y$$`

Where `$(\Psi'\Psi)^{-1}\Psi'$` is the Moore-Penrose matrix inverse

We can apply Chebyshev regression to even our regular tensor approaches, this has the advantage of dropping higher order terms which often oscillate due to error, giving us a smoother approximation

---

# Chebyshev regression: practice

Go back to our original VFI example and convert it to a regression approach

```julia
 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))*.75,
                 capital_upper = (0.75*0.95)^(1/(1 - 0.75))*1.5, capital_lower = (0.75*0.95)^(1/(1 - 0.75))/2,
                 num_basis = 7, num_points = 15, tolerance = 1e-6, fin_diff = 1e-6, mpi_start = 5);
 
 coefficients = zeros(params.num_basis);
 coefficients[1:2] = [100 5];
```

---

# Chebyshev regression: practice

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

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

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

```
## 15-element Vector{Float64}:
##   0.9945218953682733
##   0.9510565162951535
##   0.8660254037844387
##   0.7431448254773942
##   0.5877852522924731
##   0.4067366430758004
##   0.20791169081775945
##   2.83276944882399e-16
##  -0.20791169081775912
##  -0.4067366430758001
##  -0.587785252292473
##  -0.7431448254773941
##  -0.8660254037844387
##  -0.9510565162951535
##  -0.9945218953682733
```

```julia
 expand_grid(grid, params) = (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)
```

```
## 15-element Vector{Float64}:
##  0.38586640773961633
##  0.3802655705208513
##  0.3693086795455801
##  0.35347460352641824
##  0.33345536756572974
##  0.31012590833794895
##  0.28450583515849537
##  0.25771486816406236
##  0.2309239011696293
##  0.20530382799017577
##  0.18197436876239492
##  0.16195513280170648
##  0.1461210567825446
##  0.1351641658072734
##  0.1295633285885084
```

---

# Chebyshev regression: practice

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` will inherit `params` from wrapper functions

---

# Chebyshev regression: practice

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

---

# Chebyshev regression: practice

```julia
 construct_basis_matrix(grid, params) = hcat([cheb_polys.(shrink_grid.(grid), n) for n = 0:params.num_basis - 1]...);
 basis_matrix = construct_basis_matrix(capital_grid, params);
 basis_inverse = inv(basis_matrix'*basis_matrix)*(basis_matrix') # pre-compute pseudoinverse for regressions
```

```
## 7×15 Matrix{Float64}:
##  0.0666667   0.0666667     0.0666667    …   0.0666667     0.0666667
##  0.132603    0.126808      0.11547         -0.126808     -0.132603
##  0.13042     0.107869      0.0666667        0.107869      0.13042
##  0.126808    0.0783714     4.46986e-16     -0.0783714    -0.126808
##  0.121806    0.0412023    -0.0666667        0.0412023     0.121806
##  0.11547    -8.47892e-17  -0.11547      …  -1.85126e-16  -0.11547
##  0.107869   -0.0412023    -0.133333        -0.0412023     0.107869
```

---

# Chebyshev regression: practice

```julia
 eval_value_function(coefficients, grid, params) = construct_basis_matrix(grid, params) * coefficients;
```

---

# Chebyshev regression: practice

```julia
 function loop_grid_regress(params, capital_grid, coefficients)
     max_value = -.0*ones(params.num_points);
     consumption_store = -.0*ones(params.num_points);
 
     for (iteration, capital) in enumerate(capital_grid)
         function bellman(consumption)
             capital_next = capital^params.alpha - consumption
             cont_value = eval_value_function(coefficients, capital_next, params)[1]
             value_out = (consumption)^(1-params.eta)/(1-params.eta) + params.beta*cont_value
             return -value_out
         end;
 
         results = optimize(bellman, 0.00*capital^params.alpha, 0.99*capital^params.alpha)
         max_value[iteration] = -Optim.minimum(results)
         consumption_store[iteration] = Optim.minimizer(results)
     end
 
     return max_value, consumption_store
 end;
```

---

# Chebyshev regression: practice

```julia
 function solve_vfi_regress(params, basis_inverse, capital_grid, coefficients)
 
     max_value = -.0*ones(params.num_points);
     error = 1e10;
     value_prev = .1*ones(params.num_points);
     coefficients_store = Vector{Vector}(undef, 1)
     coefficients_store[1] = coefficients
     iteration = 1
 
     while error > params.tolerance
         max_value, consumption_store = loop_grid_regress(params, capital_grid, coefficients)
         coefficients = basis_inverse*max_value
         error = maximum(abs.((max_value - value_prev)./(value_prev)))
         value_prev = deepcopy(max_value)
         if mod(iteration, 25) == 0
             println("Maximum Error of $(error) on iteration $(iteration).")
             append!(coefficients_store, [coefficients])
         end
         iteration += 1
     end
 
     return coefficients, max_value, coefficients_store
 end;
```
---

# Chebyshev regression: practice

```julia
 @time solution_coeffs, max_value, intermediate_coefficients =
     solve_vfi_regress(params, basis_inverse, capital_grid, coefficients)
```

```
## Maximum Error of 0.04560791678414923 on iteration 25.
## Maximum Error of 0.007635436575597669 on iteration 50.
## Maximum Error of 0.0019075236037348097 on iteration 75.
## Maximum Error of 0.0005149316099488123 on iteration 100.
## Maximum Error of 0.00014178153569906261 on iteration 125.
## Maximum Error of 3.92482976918936e-5 on iteration 150.
## Maximum Error of 1.0880896630296866e-5 on iteration 175.
## Maximum Error of 3.0177727176252443e-6 on iteration 200.
##   0.847524 seconds (22.49 M allocations: 515.604 MiB, 5.23% gc time, 0.25% compilation time)
```

```
## ([-194.85536958622183, 14.142104524187651, -2.664424683176605, 0.5749549884000286, -0.1333725115671519, 0.03457002344598274, -0.008458351978988204], [-182.9478828542588, -183.26236062704822, -183.89615953570603, -184.85865884335664, -186.1638852575174, -187.83024816516803, -189.87710126487326, -192.31615989435315, -195.13870272444962, -198.2999931535903, -201.69819337761513, -205.14355420741794, -208.32841615695403, -210.83615415431385, -212.233073576706], Vector[[100.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0], [-111.84446088811816, 14.101018923457382, -2.6579296540031114, 0.5738072224183246, -0.13323996387170098, 0.03452677955052085, -0.008431681338534515], [-171.83327647230453, 14.142022326296566, -2.664411611854907, 0.574952651575576, -0.13337223722756203, 0.03456993567525585, -0.008458298282616294], [-188.4717369787631, 14.142104360071517, -2.66442465707678, 0.5749549837338038, -0.13337251101925873, 0.034570023270718855, -0.008458351871781953], [-193.08706876568482, 14.14210452386, -2.6644246831245004, 0.5749549883907289, -0.13337251156606242, 0.03457002344563234, -0.008458351978754048], [-194.36731367254603, 14.142104524186939, -2.664424683176515, 0.574954988400002, -0.13337251156715654, 0.03457002344598024, -0.008458351978990098], [-194.72244026073702, 14.142104524187669, -2.6644246831765925, 0.5749549883999996, -0.13337251156716917, 0.03457002344598161, -0.008458351978986358], [-194.82094867343935, 14.14210452418768, -2.664424683176591, 0.5749549883999837, -0.13337251156715232, 0.034570023445998875, -0.008458351978979827], [-194.8482738799878, 14.14210452418769, -2.664424683176578, 0.5749549884000195, -0.1333725115671582, 0.03457002344596549, -0.00845835197898225]])
```

---

# Chebyshev regression: practice

```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
             cont_value = eval_value_function(solution_coeffs, capital_next, params)[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;
```

---

# Chebyshev regression: practice

---

# Endogenous grid method (Carroll, 2006)

Suppose now we are working with a model with an inelastic labor supply with logarithmic utility `$\eta=1$`, and capital that fully depreciates

Leisure does not enter the utility function nor does labor enter the production function, i.e. `$B = 0, l = 1$`

This yields closed form solutions to the model
`\begin{align}
 	k_{t+1} &= \beta \alpha \theta_t k_t^\alpha \notag\\
 	c_t &= (1-\beta\alpha)\theta_t k_t^\alpha
\end{align}`

The endogenous grid method was introduced by Carroll (2006) for value function iteration

---

# Endogenous grid method (Carroll, 2006)

The idea behind EGM is .hi-blue[super simple]

instead of constructing a grid on the current states, construct the grid on .hi-red[future states] (making current states endogenous)

This works to our advantage because typically it is easier to solve for `$k$` given `$k'$` than the reverse

Let's see how this works

---

# Endogenous grid method

1. Choose a grid `$\left\{k_m',\theta_m\right\}_{m=1,...,M}$` on which the value function is approximated
2. Choose nodes `$\epsilon_j$` and weights `$\omega_j$`, `$j=1,...,J$` for approximating integrals.
3. Compute next period productivity, `$\theta'_{m,j} = \theta_m^\rho exp(\epsilon_j)$`.
4. Solve for `$b$` and `$\left\{ c_m,k_m \right\}$` such that
  + (inner loop) The quantities `$\left\{c_m,k_m	\right\}$` solve the following given `$V(k'_m,\theta'_m)$`:
      - `$u'(c_m) = \beta E\left[ V_k(k'_m,\theta'_{m,j}) \right]$`,
      - `$c_m + k'_m = \theta_m f(k_m) + (1-\delta)k_m$`
  + (outer loop) The value function `$\hat{V}(k,\theta;b)$` solves the following given `$\{ c_m,k'_m \}$`:
      - `$\hat{V}(k_m,\theta_m;b) = u(c_m) + \beta \sum_{j=1}^J \omega_j \left[\hat{V}(k'_m,\theta'_{m,j};b)\right]$`

---

# Endogenous grid method

.hi-blue[Focus the inner loop of VFI:]  
+ (inner loop) The quantities `$\left\{c_m,k_m	\right\}$` solve the following given `$V(k'_m,\theta'_m)$`:
    - `$u'(c_m) = \beta E\left[ V_k(k'_m,\theta'_{m,j}) \right]$`,
    - `$c_m + k'_m = \theta_m f(k_m) + (1-\delta)k_m$`

Notice that the values of `$k'$` are fixed since they are grid points

This means that we can pre-compute the expectations of the value function and value function derivatives and let `$W(k',\theta) = E[V(k',\theta';b)]$`

We can then use the consumption FOC to solve for consumption, `$c = [\beta W_k(k',\theta)]^{-1/\gamma}$` and then rewrite the resource constraint as,
`$$(1-\delta)k + \theta k^\alpha = [\beta W_k(k',\theta)]^{-1/\gamma} + k'$$`

---

# Endogenous grid method

This is easier to solve than the necessary conditions we would get out of standard value function iteration
`$$(k'-(1-\delta)k - \theta k^\alpha)^{-\gamma} = \beta W_k(k',\theta')$$`

Why?

We do not need to do any interpolation ( `$k'$` is on our grid)

We do not need to approximate a conditional expectation (already did it before hand and can do it with very high accuracy since it is a one time cost)

Can we make the algorithm better?

---

# Endogenous grid method: turbo speed

Let's make a change of variables
`$$Y \equiv (1-\delta)k + \theta k^\alpha = c + k'$$`

so we can rewrite the Bellman as
`\begin{gather}
	V(Y,\theta) = \max_{k'} \left\{ \frac{c^{1-\gamma}-1}{1-\gamma} + \beta E\left[ V(Y',\theta') \right] \right\} \notag\\
	\text{s.t.} \,\,\, c = Y - k' \notag\\
	Y' = (1-\delta)k' + \theta'(k')^\alpha \notag
\end{gather}`

---

# Endogenous grid method: turbo speed

This yields the FOC
`$$u'(c) = \beta E\left[ V_Y(Y',\theta') (1-\delta + \alpha\theta'(k')^{\alpha-1}) \right]$$`

`$Y'$` is a simple function of `$k'$` (our grid points) so we can compute it, and the entire conditional expectation on the RHS, directly from the endogenous grid points

---

# Endogenous grid method: turbo speed

`$$u'(c) = \beta E\left[ V_Y(Y',\theta') (1-\delta + \alpha\theta'(k')^{\alpha-1}) \right]$$`

This allows us to compute `$c$` from the FOC

Then from `$c$` we can compute `$Y = c + k'$` and then `$V(Y,\theta)$` from the Bellman

At no point did we need to use a numerical solver

Once we have converged on some `$\hat{V}^*$` we then solve for `$k$` via `$Y = (1-\delta)k + \theta k^\alpha$` which does require a solver, but only once and after we have recovered our value function approximant

---

# Endogenous grid method: practice

Let's solve our previous basic growth model using EGM

```julia
 coefficients = zeros(params.num_basis);
 coefficients[1:2] = [100 5];
```

---

# Endogenous grid method: practice

```julia
 function loop_grid_egm(params, capital_grid, coefficients)
 
     max_value = similar(capital_grid)
     capital_store = similar(capital_grid)
 
     for (iteration, capital_next) in enumerate(capital_grid)
 
         function bellman(consumption)
             cont_value = eval_value_function(coefficients, capital_next, params)[1]
             value_out = (consumption)^(1-params.eta)/(1-params.eta) + params.beta*cont_value
             return value_out
         end;
         value_deriv = (eval_value_function(coefficients, capital_next + params.fin_diff, params)[1] -
             eval_value_function(coefficients, capital_next - params.fin_diff, params)[1])/(2params.fin_diff)
         consumption = (params.beta*value_deriv)^(-1/params.eta)
         max_value[iteration] = bellman(consumption)
         capital_store[iteration] = (capital_next + consumption)^(1/params.alpha)
     end
 
     grid = shrink_grid.(capital_store)
     basis_matrix = [cheb_polys.(grid, n) for n = 0:params.num_basis - 1];
     basis_matrix = hcat(basis_matrix...)
     return basis_matrix, capital_store, max_value
 
 end
```

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

---

# Endogenous grid method: practice

```julia
 function solve_egm(params, capital_grid, coefficients)
     iteration = 1
     error = 1e10;
     max_value = -.0*ones(params.num_points);
     value_prev = .1*ones(params.num_points);
     coefficients_store = Vector{Vector}(undef, 1)
     coefficients_store[1] = coefficients
     while error > params.tolerance
         coefficients_prev = deepcopy(coefficients)
         current_poly, current_capital, max_value =
             loop_grid_egm(params, capital_grid, coefficients)
         coefficients = current_poly\max_value
         error = maximum(abs.((max_value - value_prev)./(value_prev)))
         value_prev = deepcopy(max_value)
         if mod(iteration, 25) == 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_egm (generic function with 1 method)
```
No need to pass basis function matrices or `$[-1,1]$` grid since we will be constructing an endogenous grid on time `$t$` capital

---

# Endogenous grid method: practice

```julia
 @time solution_coeffs, max_value, intermediate_coefficients = solve_egm(params, capital_grid, coefficients)
```

```
## Maximum Error of 0.04656312802519048 on iteration 25.
## Maximum Error of 0.0077279558439712175 on iteration 50.
## Maximum Error of 0.0019283105651684272 on iteration 75.
## Maximum Error of 0.0005203908554999235 on iteration 100.
## Maximum Error of 0.00014327347852089944 on iteration 125.
## Maximum Error of 3.966043939472439e-5 on iteration 150.
## Maximum Error of 1.0995091700787616e-5 on iteration 175.
## Maximum Error of 3.0494442960141223e-6 on iteration 200.
##   0.229561 seconds (5.69 M allocations: 144.159 MiB, 9.64% gc time, 0.76% compilation time)
```

```
## ([-194.86588167567055, 14.166854450284145, -2.659830643535021, 0.5619970720353987, -0.1363231862642804, 0.042584891304797305, -0.008520257414629136], [-181.06657322568805, -181.43461128908683, -182.1747231895469, -183.304987257002, -184.86176848996226, -186.88566388459972, -189.39801023315908, -192.3909157214092, -195.8381047584532, -199.70210567344222, -203.90039227218548, -208.22572768575577, -212.27795990787445, -215.49012260822093, -217.28206934422852], Vector[[100.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0], [-111.84269630136266, 14.118165997839721, -2.650785400463155, 0.5609056021882664, -0.1364992921461478, 0.04248754631691099, -0.008438668682153542], [-171.8414355773418, 14.166739903739217, -2.659808413026581, 0.5619942072650489, -0.13632352557703392, 0.04258464026804016, -0.008520061862729576], [-188.48159997738063, 14.16685418293722, -2.659830591621627, 0.5619970653438966, -0.13632318706504137, 0.042584890719581464, -0.008520256956436444], [-193.0974010625483, 14.166854449659024, -2.6598306434084686, 0.5619970720132431, -0.13632318627523834, 0.04258489130363603, -0.008520257410556667], [-194.3777761387717, 14.166854450281805, -2.6598306435297983, 0.56199707203115, -0.13632318627409784, 0.0425848913055752, -0.008520257413112204], [-194.73293883456404, 14.16685445029134, -2.659830643529197, 0.5619970720300576, -0.13632318627469658, 0.042584891306481895, -0.008520257412638278], [-194.83145726313631, 14.166854450277802, -2.6598306435420427, 0.5619970720345356, -0.1363231862644243, 0.042584891304662566, -0.008520257414191656], [-194.85878524798287, 14.166854450290375, -2.659830643529872, 0.5619970720300339, -0.1363231862764906, 0.042584891305316105, -0.00852025741245583]])
```

---

# Endogenous grid method: practice

```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
             cont_value = eval_value_function(solution_coeffs, capital_next, params)[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;
```

---

# Endogenous grid method: practice

---

# Envelope condition method

We can simplify rootfinding in an alternative way than an endogenous grid for infinite horizon problems

The idea here is that we want to use the envelope conditions instead of FOCs to construct policy functions

These will end up being easier to solve and sometimes we can solve them in closed form

---

# Envelope condition method

For our old basic growth model problem (fully depreciating capital, no tech) the envelope condition (combined with the consumption FOC) is given by
`$$V_k(k) = u'(c)f'(k)$$`

Notice that the envelope condition is an intratemporal condition,  
it only depends on time `$t$` variables

We can use it to solve for `$c$` as a function of current variables
`$$c = \left( \frac{V_k(k)}{\alpha k^{\alpha-1}} \right)^{-1/\eta}$$`

We can then recover `$k'$` from the budget constraint given our current state

We never need to use a solver at any point in time!

---

# Envelope condition method

The algorithm is
1. Choose a grid `$\left\{k_m\right\}_{m=1,...,M}$` on which the value function is approximated
2. Solve for `$b$` and `$\left\{ c_m,k'_m \right\}$` such that
  - (inner loop) The quantities `$\left\{c_m,k'_m	\right\}$` solve the following given `$V(k_m)$`:
    + `$V_k(k_m) = u'(c_m)f'(k_m)$`,
    + `$c_m + k'_m = f(k_m)$`
  - (outer loop) The value function `$\hat{V}(k;b)$` solves the following given `$\{ c_m,k_m \}$`:
    + `$\hat{V}(k_m;b) = u(c_m) + \beta \sum_{j=1}^J \omega_j \left[\hat{V}(k'_m;b)\right]$`

---

# Envelope condition method

In more complex settings (e.g. elastic labor supply) we will not necessarily be able to solve for policies without a solver

However we will generally be able to solve a system of conditions via function iteration to recover the optimal controls as a function of current states and future states that are perfectly known at the current time

Thus at no point in time during the value function approximation algorithm do we need to interpolate off the grid or approximate expectations: this yields large speed and accuracy gains

---

# Envelope condition method: practice

```julia
 function loop_grid_ecm(params, capital_grid, coefficients)
 
     max_value = similar(capital_grid);
 
     for (iteration, capital) in enumerate(capital_grid)
 
         function bellman(consumption)
             capital_next = capital^params.alpha - consumption
             cont_value = eval_value_function(coefficients, capital_next, params)[1]
             value_out = (consumption)^(1-params.eta)/(1-params.eta) + params.beta*cont_value
             return value_out
         end;
 
         value_deriv = (eval_value_function(coefficients, capital + params.fin_diff, params)[1] -
             eval_value_function(coefficients, capital - params.fin_diff, params)[1])/(2params.fin_diff)
         consumption = (value_deriv/(params.alpha*capital^(params.alpha-1)))^(-1/params.eta)
         consumption = min(consumption, capital^params.alpha)
         max_value[iteration] = bellman(consumption)
 
     end
 
     return max_value
 end
```

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

---

# Envelope condition method: practice

```julia
 function solve_ecm(params, basis_inverse, capital_grid, coefficients)
     iteration = 1
     error = 1e10;
     max_value = similar(capital_grid);
     value_prev = .1*ones(params.num_points);
     coefficients_store = Vector{Vector}(undef, 1)
     coefficients_store[1] = coefficients
     while error > params.tolerance
         coefficients_prev = deepcopy(coefficients)
         max_value = loop_grid_ecm(params, capital_grid, coefficients)
         coefficients = basis_inverse*max_value
         error = maximum(abs.((max_value - value_prev)./(value_prev)))
         value_prev = deepcopy(max_value)
         if mod(iteration, 25) == 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_ecm (generic function with 1 method)
```

---

# Envelope condition method: practice

```julia
 @time solution_coeffs, max_value, intermediate_coefficients =
     solve_ecm(params, basis_inverse, capital_grid, coefficients)
```

```
## Maximum Error of 0.0453525403650495 on iteration 25.
## Maximum Error of 0.0076079815408730215 on iteration 50.
## Maximum Error of 0.0019013403845842475 on iteration 75.
## Maximum Error of 0.0005133070957501089 on iteration 100.
## Maximum Error of 0.00014133753500579239 on iteration 125.
## Maximum Error of 3.912563883912878e-5 on iteration 150.
## Maximum Error of 1.0846910981672597e-5 on iteration 175.
## Maximum Error of 3.0083484348532563e-6 on iteration 200.
##   0.174326 seconds (4.53 M allocations: 106.518 MiB, 5.97% gc time, 0.99% compilation time)
```

```
## ([-194.85531932176127, 14.142062593106905, -2.6644837015279976, 0.5749531960546624, -0.13337430101896322, 0.034551695371124104, -0.008484748971169142], [-182.94799968004708, -183.26235258561604, -183.89612956524388, -184.85864174862832, -186.16385889769955, -187.83020022777666, -189.87700857244968, -192.3160354679106, -195.1385722506377, -198.29988788664537, -201.69811571031212, -205.1434915535787, -208.3283534443464, -210.8360937491323, -212.2330484863955], Vector[[100.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0], [-112.10480287950833, 14.107931378395325, -2.6591791476193576, 0.57400531884664, -0.13325275180176394, 0.03451091903707708, -0.00847435508444089], [-171.90510776971806, 14.14199420909911, -2.6644730344230134, 0.5749512770627394, -0.13337405299096824, 0.03455161324815789, -0.008484727914318257], [-188.4916175567322, 14.142062456260323, -2.664483680184769, 0.5749531922147479, -0.1333743005216747, 0.0345516952063283, -0.00848474892937651], [-193.09253945121353, 14.142062592834266, -2.664483701486141, 0.5749531960469292, -0.13337430101847123, 0.034551695370051864, -0.008484748972061315], [-194.36878720553952, 14.142062593105248, -2.664483701529717, 0.5749531960525941, -0.13337430102116343, 0.03455169536876924, -0.008484748973553577], [-194.72280502529446, 14.142062593107418, -2.6644837015275153, 0.5749531960567261, -0.13337430101720807, 0.03455169537105062, -0.008484748971776859], [-194.821005877197, 14.142062593107045, -2.6644837015291127, 0.5749531960527132, -0.13337430102148554, 0.03455169537076611, -0.008484748971882054], [-194.84824576958684, 14.142062593105138, -2.6644837015309104, 0.574953196053043, -0.1333743010207937, 0.03455169536880788, -0.008484748973840613]])
```

---

# Envelope condition method: practice

```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
             cont_value = eval_value_function(solution_coeffs, capital_next, params)[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;
```
---

# Envelope condition method: practice

```julia
 time_horizon = 100;
 consumption, capital = simulate_model(params, solution_coeffs, time_horizon);
 plot(1:time_horizon, consumption, color = :red, linewidth = 4.0, label = "Consumption", legend = :right, 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.steady_state*ones(time_horizon), color = :purple, linewidth = 2.0, linestyle = :dot, label = "Analytic Steady State")
```

---

# Modified policy iteration

When doing VFI what is the most expensive part of the algorithm?

The maximization step!

If we can reduce how often we need to maximize the Bellman we can significantly improve speed

It turns out that between VFI iterations, the optimal policy does not change all that much

This means that we may be able to skip the maximization step and re-use our old policy function to get new values for polynomial fitting

This is called **modified policy iteration**

---

# Modified policy iteration

It only changes step 5 of VFI:

While convergence criterion `$>$` tolerance
  + Start iteration `$p$`
  + Solve the right hand side of the Bellman equation
  + Recover the maximized values, conditional on `$\Gamma(k_{t+1};b^{(p)})$`
  + Fit the polynomial to the values and recover new coefficients `$\hat{b}^{(p+1)}$`
  + Compute `$b^{(p+1)} = (1-\gamma) b^{(p)} + \gamma \hat{b}^{(p+1)}$` where `$\gamma \in (0,1)$`
  + While MPI stop criterion `$>$` tolerance
      - Use policies from last VFI iteration to re-fit the polynomial (no maximizing!)
      - Compute `$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)$`

---

# Modified policy iteration

Stop criteron can be a few things:
1. Fixed number of iterations
2. Stop when change in value function is sufficient small, QuantEcon suggests stopping MPI when
`$$\max(V_p(x;c) - V_{p-1}(x;c)) -  \min(V_p(x;c) - V_{p-1}(x;c)) < \epsilon(1-\beta)\beta$$`
where the max and min are over the values on the grid

Only MPI after a few VFI iterations unless you have a good initial guess, if your early policy functions are bad then starting MPI too early will blow up your problem

---

# Modified policy iteration

```julia
 function solve_vfi_regress_mpi(params, basis_inverse, basis_matrix, grid, capital_grid, coefficients)
     max_value = -.0*ones(params.num_points);
     error = 1e10;
     value_prev = .1*ones(params.num_points);
     value_prev_outer = .1*ones(params.num_points);
     coefficients_store = Vector{Vector}(undef, 1)
     coefficients_store[1] = coefficients
     iteration = 1
     while error > params.tolerance
         max_value, consumption_store =
             loop_grid_regress(params, capital_grid, coefficients)
         coefficients = basis_inverse*max_value
         if iteration > params.mpi_start # modified policy iteration loop
             mpi_iteration = 1
             while maximum(abs.(max_value - value_prev)) - 
                     minimum(abs.(max_value - value_prev)) > 
                     (1 - params.beta)/params.beta*params.tolerance
                 value_prev = deepcopy(max_value)
```

---

# Modified policy iteration

```julia
               function bellman(consumption, capital)
                     capital_next = capital^params.alpha - consumption
                     cont_value = eval_value_function(coefficients, capital_next, params)[1]
                     value_out = (consumption)^(1-params.eta)/(1-params.eta) + params.beta*cont_value
                     return value_out
                 end
                 max_value = bellman.(consumption_store, capital_grid) # greedy policy
                 coefficients = basis_inverse*max_value
                 if mod(mpi_iteration, 5) == 0
                     println("MPI iteration $mpi_iteration on VFI iteration $iteration.")
                 end
                 mpi_iteration += 1
             end
         end
         error = maximum(abs.((max_value .- value_prev_outer)./(value_prev_outer)))
         value_prev_outer = deepcopy(max_value)
 
         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;
```

---

# Modified policy iteration

---

# Modified policy iteration

```julia
 @time solution_coeffs, max_value, intermediate_coefficients =
     solve_vfi_regress_mpi(params, basis_inverse, basis_matrix, grid, capital_grid, coefficients)
```

```
## Maximum Error of 0.33871304913135464 on iteration 5.
## MPI iteration 25 on VFI iteration 6.
## MPI iteration 50 on VFI iteration 6.
## MPI iteration 75 on VFI iteration 6.
## MPI iteration 25 on VFI iteration 7.
## MPI iteration 50 on VFI iteration 7.
## MPI iteration 25 on VFI iteration 8.
## Maximum Error of 1.141822859504868e-5 on iteration 10.
## Maximum Error of 4.0892905314530945e-6 on iteration 15.
## Maximum Error of 3.0060755201005212e-6 on iteration 20.
## Maximum Error of 2.3260768798925272e-6 on iteration 25.
## Maximum Error of 1.7998932693973723e-6 on iteration 30.
## Maximum Error of 1.3927345521584257e-6 on iteration 35.
## Maximum Error of 1.0776782693136323e-6 on iteration 40.
##   0.223322 seconds (5.89 M allocations: 136.588 MiB, 4.52% gc time, 0.88% compilation time)
```

```
## ([-194.8621441678187, 14.1421045241982, -2.6644246831782934, 0.5749549884003013, -0.1333725115671613, 0.03457002344599215, -0.008458351978991155], [-182.95465743584654, -183.26913520863621, -183.9029341172946, -184.86543342494596, -186.17065983910774, -187.8370227467597, -189.88387584646654, -192.32293447594827, -195.14547730604696, -198.30676773519, -201.70496795921736, -205.15032878902258, -208.3351907385608, -210.84292873592221, -212.23984815831525], Vector[[100.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0], [36.80824126072954, 8.655027968115757, -1.4832910460371644, 0.3154836954815494, -0.07041869937297221, 0.015492390858873, -0.003125302805521049], [-194.87810217521672, 14.142104573873516, -2.664424691399574, 0.5749549899065862, -0.13337251176021203, 0.034570023507038646, -0.008458352013435776], [-194.87297913249904, 14.142104535322538, -2.6644246849689233, 0.5749549887223999, -0.13337251160585856, 0.03457002345832334, -0.008458351986372807], [-194.86921322872692, 14.142104526689847, -2.6644246835760215, 0.5749549884715544, -0.13337251157560837, 0.034570023448690965, -0.008458351980620617], [-194.86684924389925, 14.142104524909199, -2.6644246832915175, 0.5749549884205752, -0.13337251156957403, 0.03457002344675323, -0.008458351979455069], [-194.8650200376443, 14.14210452439574, -2.664424683209724, 0.5749549884059388, -0.13337251156784122, 0.03457002344620518, -0.008458351979124092], [-194.86360463275415, 14.142104524247642, -2.6644246831861524, 0.5749549884017184, -0.13337251156734098, 0.034570023446056024, -0.008458351979021731], [-194.86250941944303, 14.142104524204983, -2.6644246831793468, 0.5749549884005244, -0.13337251156720537, 0.034570023446011434, -0.00845835197899945]])
```

---

# Modified policy iteration

---

# Solve times

```
## The solve times are:
```

```
## Regression: 0.867 seconds
```

```
## Endogenous Grid + Regression: 0.216 seconds
```

```
## Envelope Condition + Regression: 0.172 seconds
```

```
## Modified Policy Iteration + Regression: 0.257 seconds
```

Individually they give 3-6 times speed up