# Solving cass-koopmans planning model with python value function iteration

edit: sorry, I used stack exchange first time, so I didn't know I can type python code. you can now copy and run the code.

I'm practicing solving economic model using python programming. I refer to https://python.quantecon.org/optgrowth.html and https://python.quantecon.org/optgrowth_fast.html.

I wanna solve Cass-Koopmans planning problem with bellman equation.

• There is only one agent who governs all resource allocations in the economy. She produces a (single) good through a production function:

$$y_{t} = f(k_{t}) \tag{1}$$

• The good can be used for consumption and investment.

• In this exercise, we introduce a partial capital depreciation, i.e. capital depreciates at the rate $$0<\delta<1$$ each period. Hence the resource constraint for the agent is:

$$k_{t+1} + c_t = y_t + (1-\delta)k_t \tag{2}$$

• The agent wants to maximize

$$\mathbb E \left[ \sum_{t = 0}^{\infty} \beta^t u(c_t) \right] \tag{3}$$ subject to (2), where $$\beta \in (0, 1)$$ is a discount factor.

• The model is almost identical to Cass-Koopmans planning model. In order to solve this model, we will apply value function iteration algorithm. Hence, we first reformulate the maximization problem into a Bellman equation.

We can write the value function for the agent's utility maximization problem in the form of Bellman equation. In this exercise we will use capital ($$k$$) as state variable instead of output ($$y$$). The value function is defined as follows:

$$v(k) = \max_{0 \leq c \leq y + (1-\delta)k} \left\{ u(c) + \beta v(k') \right\} \tag{4}$$ subject to $$k' = f(k) + (1-\delta)k - c$$

This formulation takes consumption ($$c$$) as control variable. For the computation below we will use the next period capital ($$k'$$) as control variable. Then (4) can be rewritten as follows:

$$v(k) = \max_{0 \leq k' \leq y + (1-\delta)k} \left\{ u\left(f(k) + (1-\delta)k - k'\right) + \beta v(k') \right\} \tag{5}$$

Essentially we converted the constrained maximization problem in (4) into the unconstrained maximization problem in (5).

So I made my python code:

import numpy as np
import matplotlib.pyplot as plt
from numba import njit, float64
from numba.experimental import jitclass
from quantecon.optimize.scalar_maximization import brent_max
from interpolation import interp

opt_growth_data = [('α', float64),
('β', float64),
('γ', float64),
('δ', float64),
('grid',float64[:])]

@jitclass(opt_growth_data)
class OptimalGrowth_VI:
def __init__(self, α=0.4, β=0.96, γ=2.0, δ=0.1, grid_max=10, grid_size=500):
self.α, self.β, self.γ, self.δ = α, β, γ, δ
self.grid = np.linspace(0.1, grid_max, grid_size)

def f(self, k):
return k**self.α

def u(self, c):
return c**(1 - self.γ) / (1 - self.γ)

def objective(self, k, kp, v_array):
f, u, β, δ = self.f, self.u, self.β, self.δ
v = lambda x: interp(self.grid, v_array, x)

return u(f(k)+(1-δ)*k-kp) + β*v(kp)

@njit
def T(v, og_VI):
v_greedy = np.empty_like(v)
v_new = np.empty_like(v)

for i in range(len(og_VI.grid)):
k = og_VI.grid[i]
lower = 1e-10
upper = og_VI.f(k) + (1-og_VI.δ)*k
result = brent_max(og_VI.objective, lower, upper, args=(k,v))
v_greedy[i], v_new[i] = result[0], result[1]

return v_greedy, v_new

def solve_model_VI(og_VI, tol=1e-4, max_iter=1000, print_skip=20):
v = og_VI.grid
error = tol+1
i=0

while i < max_iter and error > tol:
v_greedy, v_new = T(v, og_VI)
error = np.max(np.abs(v - v_new))
i += 1
if i % print_skip == 0:
print(f"Error at iteration {i} is {error}.")
v = v_new

if i == max_iter:
print("Failed to converge!")

if i < max_iter:
print(f"\nConverged in {i} iterations.")

return v_greedy, v_new

og_VI = OptimalGrowth_VI()
v_greedy, v_solution = solve_model_VI(og_VI)
plt.plot(og_VI.grid, v_greedy)


I thought that I made my code correctly, but it produced weird optimal k's graph. Could you review my code and tell me what's wrong with my implementation? It's a nuisance, but any help would be really appreciated.

• Try fitting a curve instead of a line to your graph and I suspect you should be able to get the usual path. If you can post your code, I' ll take a look at it. You might also want to look into stack-overflow for this.
– Rumi
Commented Apr 20, 2022 at 12:28
• It would be great to add the code in readable font Commented Apr 21, 2022 at 8:28
• Please instead of posting pictures of code copy code as code you can do that by using two  around the code . I think most people who would be willing to help won’t want to waste time copying code from pictures to try to run it to see where the problem is
– 1muflon1
Commented Apr 21, 2022 at 9:26
• I made it runnable. Sorry for the mistake. Commented Apr 21, 2022 at 10:45

I tried the code you provided.

1. The problem seems to occur because, in those points, consumption turns negative and messes up with the utility function. So, in your Class OptimalGrowth_VI definition, I added the following piece for the objective method:
    def objective(self, k, kp, v_array):
f, u, β, δ = self.f, self.u, self.β, self.δ
v = lambda x: interp(self.grid, v_array, x)
c = f(k)+(1-δ)*k-kp
if c <= 0:
u = -888 - 800 * abs(c)
else:
u = u(c)
return u + β*v(kp)


Basically, if c is negative, it will return an extremely small utility, so the brent_max function will ignore it. (McCandless' treatment). After that, you will have the textbook curved policy function.

1. Another thing I noticed is that the policy function, according to the grid and parameters you defined, does not touch the 45-degree line. Perhaps the predefined grid does not contain the steady-state value. I increased your grid_max= 50 and added a 45-degree line to check the convergence to the steady state.
og_VI = OptimalGrowth_VI()
v_greedy, v_solution = solve_model_VI(og_VI)
plt.plot(og_VI.grid, v_greedy)
`