Continuous State DP Continued
We solve the same model as the previous lecture
Once again, we use fitted value iteration
However, we now use linear interpolation instead of step functions
In one dimensional concave models, this can give a better fit
We can now present our approximate Bellman operator
As before, corresponds to the model
## Filename fvi.py
## Author: John Stachurski
from scipy import linspace, mean, exp, randn
from scipy.optimize import fminbound
from lininterp import LinInterp # From listing 6.4
theta, alpha, rho = 0.5, 0.8, 0.9 # Parameters
def U(c): return 1 - exp(- theta * c) # Utility
def f(k, z): return (k**alpha) * z # Production
W = exp(randn(1000)) # Draws of shock
gridmax, gridsize = 8, 150
grid = linspace(0, gridmax**1e-1, gridsize)**10
def maximum(h, a, b):
return h(fminbound(lambda x: -h(x), a, b))
def bellman(w):
"""The approximate Bellman operator.
Parameters: w is a vectorized function (i.e., a callable object
which acts pointwise on arrays).
Returns: An instance of LinInterp.
"""
vals = []
for y in grid:
h = lambda k: U(y - k) + rho * mean(w(f(k,W)))
vals.append(maximum(h, 0, y))
return LinInterp(grid, vals)
This is listing 6.5 in EDTC
- See the discussion there for detailed explanation
- Note that the grid is constructed so that there are many points are close to zero
- More curvature near zero, so need more grid points
- Integral is computed via Monte Carlo, but quadrature (
quad, etc.) also works fine
Problem:
Generate iterates of the approximate Bellman operator and plot them
Start from initial condition w = U
Iterate until deviation between current and previous iterate is less than 0.005
Here, deviation between functions v and w is measured by
When these functions are instances of
LinInterp, the supremum equals the supremum over the grid points alone
You should produce the following graph
Solution:
from fvi import *
from scipy import absolute as abs
import matplotlib.pyplot as plt
v = U # Initial condition
tol = 0.005 # Error tolerance
while 1:
plt.plot(grid, v(grid), 'k-')
new_v = bellman(v)
if max(abs(new_v(grid) - v(grid))) < tol:
break
v = new_v
plt.show()
Problem:
Compute an approximate optimal policy from the last of these iterates
- Compute the value of the policy at the grid points
- Plot the result
Your figure should look like this:
Note that, for our parameters, all income is saved when income is low
Solution:
from fvi import *
from scipy import absolute as abs
import matplotlib.pyplot as plt
v = U # Initial condition
tol = 0.005 # Error tolerance
while 1:
new_v = bellman(v)
if max(abs(new_v(grid) - v(grid))) < tol:
break
v = new_v
def maximizer(h, a, b):
g = lambda x: - h(x) # Negate
return fminbound(g, a, b) # and minimize
policy = []
for y in grid:
h = lambda k: U(y-k) + rho * mean(new_v(f(k,W)))
policy.append(maximizer(h, 0, y))
plt.plot(grid, grid, 'k--')
plt.plot(grid, policy, 'g-')
plt.show()
0 comments:
Post a Comment