<p>I tried the code you provided.</p>
<ol>
<li>The problem seems to occur because, in those points, consumption turns negative and messes up with the utility function. So, in your Class <code>OptimalGrowth_VI</code> definition, I added the following piece for the <code>objective</code> method:</li>
</ol>
<pre><code>    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 &lt;= 0:
            u = -888 - 800 * abs(c)
        else:
            u = u(c)
        return u + β*v(kp)
</code></pre>
<p>Basically, if <code>c</code> is negative, it will return an extremely small utility, so the <code>brent_max</code> function will ignore it. (McCandless' treatment). After that, you will have the textbook curved policy function.</p>
<ol start="2">
<li>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 <code>grid_max= 50</code> and added a 45-degree line to check the convergence to the steady state.</li>
</ol>
<pre><code>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')
plt.xlabel('k(t)')
plt.ylabel('k(t+1)')
</code></pre>
<p><a href="https://i.sstatic.net/64Btf.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/64Btf.png" alt="policy function" /></a></p>

CC BY-SA 4.0
2023-09-28T06:27:55.977Z
2023-09-28T06:27:55.977Z
teddi
<p>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.</p>
<p>I'm practicing solving economic model using python programming. I refer to
<a href="https://python.quantecon.org/optgrowth.html" rel="nofollow noreferrer">https://python.quantecon.org/optgrowth.html</a> and
<a href="https://python.quantecon.org/optgrowth_fast.html" rel="nofollow noreferrer">https://python.quantecon.org/optgrowth_fast.html</a>.</p>
<p>I wanna solve Cass-Koopmans planning problem with bellman equation.</p>
<ul>
<li><p>There is only one agent who governs all resource allocations in the economy. She produces a (single) good through a production function:</p>
<p><span class="math-container">$$
y_{t} = f(k_{t}) \tag{1}
$$</span></p>
</li>
<li><p>The good can be used for consumption and investment.</p>
</li>
<li><p>In this exercise, we introduce a partial capital depreciation, i.e. capital depreciates at the rate <span class="math-container">$0&lt;\delta&lt;1$</span> each period. Hence the resource constraint for the agent is:</p>
<p>
<span class="math-container">$$
k_{t+1} + c_t = y_t + (1-\delta)k_t \tag{2}
$$</span></p>
</li>
<li><p>The agent wants to maximize</p>
<p>
<span class="math-container">$$
\mathbb E \left[ \sum_{t = 0}^{\infty} \beta^t u(c_t) \right] \tag{3}
$$</span>
subject to (2), where <span class="math-container">$ \beta \in (0, 1) $</span> is a discount factor.</p>
</li>
<li><p>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.</p>
</li>
</ul>
<p>We can write the value function for the agent's utility maximization problem in the form of <strong>Bellman equation</strong>. In this exercise we will use capital (<span class="math-container">$k$</span>) as state variable instead of output (<span class="math-container">$y$</span>). The value function is defined as follows:</p>
<p>
<span class="math-container">$$
v(k) = \max_{0 \leq c \leq y + (1-\delta)k}
\left\{ u(c) + \beta v(k') \right\} \tag{4}
$$</span>
subject to
<span class="math-container">$$
k' = f(k) + (1-\delta)k - c 
$$</span></p>
<p>This formulation takes consumption (<span class="math-container">$c$</span>) as control variable. For the computation below we will use the next period capital (<span class="math-container">$k'$</span>) as control variable. Then (4) can be rewritten as follows:</p>
<p>
<span class="math-container">$$
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}
$$</span></p>
<p>Essentially we converted the constrained maximization problem in (4) into the unconstrained maximization problem in (5).</p>
<p>So I made my python code:</p>
<pre><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
</code></pre>
<pre><code>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)
</code></pre>
<pre><code>@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
</code></pre>
<pre><code>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 &lt; max_iter and error &gt; 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&quot;Error at iteration {i} is {error}.&quot;)
        v = v_new
    
    if i == max_iter:
        print(&quot;Failed to converge!&quot;)

    if i &lt; max_iter:
        print(f&quot;\nConverged in {i} iterations.&quot;)

    return v_greedy, v_new
</code></pre>
<pre><code>og_VI = OptimalGrowth_VI()
v_greedy, v_solution = solve_model_VI(og_VI)
plt.plot(og_VI.grid, v_greedy)
</code></pre>
<p><a href="https://i.sstatic.net/A7D6A.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/A7D6A.png" alt="enter image description here" /></a>
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.</p>

CC BY-SA 4.0
2022-04-20T12:28:45.500Z
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
CC BY-SA 4.0
2022-04-21T08:28:36.410Z
It would be great to add the code in readable font
Papayapap
CC BY-SA 4.0
2022-04-21T09:26:14.420Z
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
CC BY-SA 4.0
2022-04-21T10:45:55.167Z
I made it runnable. Sorry for the mistake.
dueun
<macroeconomics><recursive-macroeconomics>
Solving cass-koopmans planning model with python value function iteration
dueun
# Solving cass-koopmans planning model with python value function iteration

(Asked by: username_0 on 2022-04-20)

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)

```

[![enter image description here](https://i.sstatic.net/A7D6A.png)](https://i.sstatic.net/A7D6A.png)
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.

**Tags:** macroeconomics, recursive-macroeconomics

> **Comments:**
> **username_1**: 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.  
> **username_2**: It would be great to add the code in readable font  
> **username_3**: 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  
> **username_0**: I made it runnable. Sorry for the mistake.  

---

## 1 answers

### Answer by username_4 on 2023-09-28. Score: 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)
        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.

2. 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')
plt.xlabel('k(t)')
plt.ylabel('k(t+1)')

```

[![policy function](https://i.sstatic.net/64Btf.png)](https://i.sstatic.net/64Btf.png)


---

