by Fedor Iskhakov, ANU
Description: Cake eating problem setup. Solution “on the grid”.
Let the flow utility be given by
$$ u(c_{t})=\log(c_t) $$Overall goal is to maximize the discounted expected utility
$$ \max_{\{c_{t}\}_{0}^{\infty}}\sum_{t=0}^{\infty}\beta^{t}u(c_{t}) \longrightarrow \max $$Value function $ V(W_t) $ = the maximum attainable value given the size of cake $ W_t $ (in period $ t $)
F.O.C. for $ c $
$$ \frac{1}{c} - \frac{\beta B}{W - c} = 0, \quad c = \frac {W} {1 + \beta B}, W - c = \frac {\beta B W} {1 + \beta B} $$Then we have
$$ A + B\log W = \log W + \log\frac{1}{1+\beta B} + \beta A + \beta B \log W + \beta B \log\frac{\beta B}{1+\beta B} $$$$ \begin{eqnarray*} A &=& \beta A + \log\frac{1}{1+\beta B} + \beta B \log\frac{\beta B}{1+\beta B} \\ B &=& 1 + \beta B \end{eqnarray*} $$After some algebra
$$ c^{\star}(W) = \frac {W} {1 + \beta B} = \frac {W} {1 + \frac{\beta}{1-\beta}} = (1-\beta)W $$$$ V(W) = \frac{\log(W)}{1-\beta} + \frac{\log(1-\beta)}{1-\beta} + \frac{\beta \log(\beta)}{(1-\beta)^2} $$The Bellman equation becomes operator in functional space
$$ T(V)(W) \equiv \max_{0 \le c \le W} \big[u(c)+\beta V(W-c)\big] $$The Bellman equations is then $ V(W) = T({V})(W) $, with the solution given by the fixed point (solution to $ T({V}) = V $)
Can you spot the potential problem?
import numpy as np
class cake_ongrid():
'''Simple class to implement cake eating problem on the grid'''
def __init__(self,beta=.9, Wbar=10, ngrid=50):
'''Initializer'''
self.beta = beta # Discount factor
self.Wbar = Wbar # Upper bound on cake size
self.ngrid = ngrid # Number of grid points
self.epsilon = np.finfo(float).eps # smallest positive float number
self.grid = np.linspace(self.epsilon,Wbar,ngrid) # grid for both state and decision space
def bellman(self,V0):
'''Bellman operator, V0 is one-dim vector of values on grid'''
c = self.grid - self.grid[:,np.newaxis] # current state in columns and choices in rows
c[c==0] = self.epsilon # add small quantity to avoid log(0)
mask = c>0 # mask off infeasible choices
matV1 = np.full((self.ngrid,self.ngrid),-np.inf) # init V with -inf
matV0 = np.repeat(V0.reshape(self.ngrid,1),self.ngrid,1) #current value function repeated in columns
matV1[mask] = np.log(c[mask])+self.beta*matV0[mask] # maximand of the Bellman equation
V1 = np.amax(matV1,axis=0,keepdims=False) # maximum in every column
c1 = self.grid - self.grid[np.argmax(matV1,axis=0)] # consumption (index of maximum in every column)
return V1, c1
def solve(self, maxiter=1000, tol=1e-4, callback=None):
'''Solves the model using VFI (successive approximations)'''
V0=np.log(self.grid) # on first iteration assume consuming everything
for iter in range(maxiter):
V1,c1=self.bellman(V0)
if callback: callback(iter,self.grid,V1,c1) # callback for making plots
if np.all(abs(V1-V0) < tol):
break
V0=V1
else: # when i went up to maxiter
print('No convergence: maximum number of iterations achieved!')
return V1,c1
model = cake_ongrid(beta=0.92,Wbar=10,ngrid=50)
V,c = model.solve()
print(c)
[0. 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.20408163 0.40816327 0.40816327]
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['axes.autolimit_mode'] = 'round_numbers'
plt.rcParams['axes.xmargin'] = 0
plt.rcParams['axes.ymargin'] = 0
plt.rcParams['patch.force_edgecolor'] = True
from cycler import cycler
plt.rcParams['axes.prop_cycle'] = cycler(color='bgrcmyk')
fig1, ax1 = plt.subplots(figsize=(12,8))
plt.grid(b=True, which='both', color='0.65', linestyle='-')
ax1.set_title('Value function convergence with VFI')
ax1.set_xlabel('Cake size, W')
ax1.set_ylabel('Value function')
def callback(iter,grid,v,c):
'''Callback function for DP solver'''
if iter<5 or iter%10==0:
ax1.plot(grid[1:],v[1:],label='iter={:3d}'.format(iter),linewidth=1.5)
V,c = model.solve(callback=callback)
plt.legend(loc=4)
# plt.savefig('cake1value.eps', format='eps', dpi=300)
plt.show()
We can find some derived theoretical property of the model and check if it holds in the computed numerical solution
fun = lambda w: np.log(w)/(1 - model.beta) + np.log(1 - model.beta)/(1 - model.beta) + model.beta*np.log(model.beta)/((1 - model.beta)**2)
grid = model.grid
plt.plot(grid[1:],V[1:],linewidth=1.5)
plt.plot(grid[1:],fun(grid[1:]),linewidth=1.5)
[<matplotlib.lines.Line2D at 0x7f8bb8ec28d0>]
apol = lambda w: (1 - model.beta) * w
grid = model.grid
plt.plot(grid[1:],c[1:],linewidth=1.5)
plt.plot(grid[1:],apol(grid[1:]),linewidth=1.5)
[<matplotlib.lines.Line2D at 0x7f8bb8f8fc90>]
m = cake_ongrid(beta=0.92,Wbar=10,ngrid=250)
V,c = m.solve()
plt.plot(m.grid[1:],c[1:],linewidth=1.5)
plt.plot(m.grid[1:],apol(m.grid[1:]),linewidth=1.5)
[<matplotlib.lines.Line2D at 0x7f8b90298510>]
Solving “on the grid” allows to avoid interpolation of the value function, but leads to huge inaccuracies for low levels of wealth!