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),

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)
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
    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)

enter image description here 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.

  • 1
    $\begingroup$ 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. $\endgroup$
    – Rumi
    Apr 20, 2022 at 12:28
  • 2
    $\begingroup$ It would be great to add the code in readable font $\endgroup$
    – Papayapap
    Apr 21, 2022 at 8:28
  • 1
    $\begingroup$ 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 $\endgroup$
    – 1muflon1
    Apr 21, 2022 at 9:26
  • $\begingroup$ I made it runnable. Sorry for the mistake. $\endgroup$
    – dueun
    Apr 21, 2022 at 10:45

1 Answer 1


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)
            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)
# add 45 degree line
plt.plot(og_VI.grid, og_VI.grid, ls=':', label='45 degree line')

policy function


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.