Instantly share code, notes, and snippets.

Embed
What would you like to do?
Endo-Labour
import numpy as np
from scipy.optimize import fminbound, brentq
from scipy import interp
class ConsumerProblem:
def __init__(self, eta=0.5, r=0.01, beta=0.96, Pi=((0.6, 0.4), (0.05, 0.95)),
z_vals=(0.5, 1.0), phi=0, grid_max=16, grid_size=50):
self.eta = eta
self.r, self.R = r, 1 + r
self.beta, self.phi = beta, phi
self.Pi, self.z_vals = np.array(Pi), tuple(z_vals)
self.asset_grid = np.linspace(-phi, grid_max, grid_size)
def coleman_operator(self, c):
# === simplify names, set up arrays === #
eta, R, Pi, beta, phi = self.eta, self.R, self.Pi, self.beta, self.phi
asset_grid, z_vals = self.asset_grid, self.z_vals
z_size = len(z_vals)
gamma = R * beta
vals = np.empty(z_size)
du=lambda x: 1/x
# === linear interpolation to get consumption function === #
def cf(a):
for i in range(z_size):
vals[i] = interp(a, asset_grid, c[:, i])
return vals
# === solve for root to get Kc === #
Kap = np.empty(c.shape)
Kc, Kl = np.empty(c.shape),np.empty(c.shape)
for i_a, a in enumerate(asset_grid):
for i_z, z in enumerate(z_vals):
#print "asset: %.2f z: %.1f" %(a,z)
def l_func(l,t):
#print "l: %.4f t: %.4f" %(l,t)
c=z*l+R*a-t
o=1-l
#print "1: %.4f 2: %.4f" %(o*z, eta*c)
return o*z - eta*c
# Find ap
def h(t):
l=brentq(l_func, 0, 1, args=(t))
#print "select l: %.4f for a': %.4f" %(l, t)
expectation = np.dot(du(cf(t)), Pi[i_z, :])
return du(z*l + R*a - t) - max(gamma * expectation, du(z*l + R*a + phi))
l_max = 1
Kap[i_a, i_z] = brentq(h, 0, R*a + z*l_max + phi)
#print "ap_star: %.4f when a%d z%d" %(Kap[i_a,i_z], i_a, i_z)
Kl[i_a, i_z] = brentq(l_func, 0, 1, args=(Kap[i_a, i_z]))
Kc[i_a, i_z] = z*Kl[i_a, i_z] + R*a - Kap[i_a, i_z]
return Kc, Kap, Kl
def initialize(self):
# === Simplify names, set up arrays === #
R, beta, phi = self.R, self.beta, self.phi
asset_grid, z_vals = self.asset_grid, self.z_vals
shape = len(asset_grid), len(z_vals)
V, c = np.empty(shape), np.empty(shape)
for i_a, a in enumerate(asset_grid):
for i_z, z in enumerate(z_vals):
l_max = 1
c_max = R*a + z*l_max + phi
c[i_a, i_z] = c_max
return c
from quantecon import mc_sample_path
import matplotlib.pyplot as plt
def compute_fixed_point(T, c, error_tol=1e-3, max_iter=50, verbose=True):
iterate = 0
error = error_tol + 1
while iterate < max_iter and error > error_tol:
new_c, new_ap, new_l = T(c)
iterate += 1
error = np.max(np.abs(new_c - c))
if verbose:
print("Computed iterate %d with error %f" % (iterate, error))
try:
c[:] = new_c
ap, l = new_ap, new_l
except TypeError:
c = new_c
ap, l = new_ap, new_l
return c, ap, l
def compute_asset_series(cp, T=500000, verbose=False):
Pi, z_vals, R = cp.Pi, cp.z_vals, cp.R # Simplify names
c_init = cp.initialize()
c, ap, l = compute_fixed_point(cp.bellman_operator, c_init, verbose=verbose)
cf = lambda a, i_z: interp(a, cp.asset_grid, c[:, i_z])
lf = lambda a, i_z: interp(a, cp.asset_grid, l[:, i_z])
a = np.zeros(T+1)
z_seq = mc_sample_path(Pi, sample_size=T)
for t in range(T):
i_z = z_seq[t]
a[t+1] = R * a[t] + z_vals[i_z]*lf(a[t], i_z) - cf(a[t], i_z)
return a
# This works
cp = ConsumerProblem(r=0.01, eta=0.1, grid_max=4)
a = compute_asset_series(cp, T=100000, verbose=True)
print a
# This doesn't
cp = ConsumerProblem(r=0.01, eta=0.2, grid_max=4)
a = compute_asset_series(cp, T=100000, verbose=True)
print a
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment