<p>In the last line of the outer loop, you are updating <code>V0</code> by assigning <code>V1</code> to it. However, doing so will make the <code>maxMet</code> calculation fail right after the first loop. It is better to copy only the values.</p>
<pre><code>V0 = V1.copy()
</code></pre>
<p>Minor comments:</p>
<ul>
<li>The value function should be written like this</li>
</ul>
<pre><code>def VF(a, a0, V0, bbeta, w, r):
    &quot;&quot;&quot;
    This function returns the interpolated value function at a given a0.
    &quot;&quot;&quot;
    # Interpolate value
    g = np.interp(a, ES, V0)
    # Consumption
    c = w + (1 + r) * a0 - a
    if c &lt;= 0:
        V = -1e10 - 1e10 * abs(c)
    else:
        V = np.log(c) + bbeta * g
    return -V
</code></pre>
<ul>
<li>Conversion seems very slow. I suggest making the <code>tol</code> around <code>1e-5</code>. Since you are interpolating anyway, choosing <code>N=1000</code> seems a bit excessive. I think <code>100</code> is enough.</li>
</ul>
<p><a href="https://i.sstatic.net/JpcfEN72.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/JpcfEN72.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2024-05-14T01:48:55.363Z
Thanks! Can you elaborate on the technicality of why this happens? What is the difference between both `V0=V1` and `V0=V1.copy()`? (thanks for the recommendations too)
manifold
CC BY-SA 4.0
2024-05-14T02:05:26.173Z
After the first loop, if you directly assign `V0=V1`, both `V0` and `V1` now refer to the same array, which means any change to `V1` will now change `V0`. Then, in the second loop, when you calculate `V1[j]`, you are actually changing `V0` and `V1` simultaneously. This is why `maxMet` in the second loop is exactly 0, so the iteration stops!


By using the `.copy()` method, you make sure that `V0` remains unchanged during the next inner loop. Since the optimized value from the last iteration is preserved until the entire array of `V1` is computed, iteration continues. Hope it helps!
teddi
CC BY-SA 4.0
2024-05-14T00:53:19.333Z
2024-05-14T01:05:09.000Z
2024-05-14T01:05:09.000Z
teddi
<p>I'm coding a program in Python to find a numerical approximation of the value function <span class="math-container">$v:\mathbb{R}\to\mathbb{R}$</span> that solves the Bellman equation:</p>
<p><span class="math-container">$$v(a)=\max\limits_{a'}\left\{\ln\left(w+(1+r)a-a'\right)+\beta v(a')\right\}$$</span></p>
<p>with <span class="math-container">$w, r&gt;0$</span> and <span class="math-container">$\beta\in(0,1)$</span> given. When I code a value function iteration algorithm in Matlab to solve this problem I get the expected result a concave and smooth <span class="math-container">$v$</span> (1st image attached). Nevertheless, coding the exact same algorithm in Python using numpy and scipy I get weird results. In particular the algorithm 'converges' on the first iteration and the resulting VF is clearly too far from a good approximation (2nd image attached).</p>
<p>I would appreciate help in finding why the Python code does not achieve the same results as the Matlab while using the same algorithm. Thank you.</p>
<p>PD: My environment is Python 3.10.14, scipy 1.13.0.</p>
<p>Python code used:</p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fminbound

# Parameters
bbeta = 0.98
w = 1
r = 0.05
tol = 1e-10 # val. fun. distance tolerante
iters = 2000 # maximum number of iterations
N = 1000 # number of point in grid
amin = -6 # minimum level of assets in grid
amax = 30 # maximum level of assets in grid

ES = np.linspace(amin, amax, N) # state space

def VF(a, a0, V0):
    &quot;&quot;&quot;
    This function returns the interpolated value function at a given a0
    &quot;&quot;&quot;
    # interpolate value
    g = np.interp(a, ES, V0)
    # consumption
    c = w + (1+r)*a0 - a
    if c&lt;=0:
        V = 1e10 + 1e10*np.abs(c)
        return V
    V = np.log(c) + bbeta*g
    return -V

V0 = np.zeros(N) # initial guess value function, this vector will be recycled
a1 = np.zeros(N) # vector that stores optimizers, i.e. approx. policy function
V1 = np.zeros(N) # vector that to store the iterated value function
iter = 1
maxMet = 100

while (iter &lt; iters) and (maxMet&gt;tol):
    for j in range(N):
        # fix the value of current point in the state space to iterate the VF
        a = ES[j]
        # apply the Bellman operator to V0
        res = fminbound(lambda x: VF(x, a, V0), amin, amax, xtol=1e-10,
                          full_output = True)
        a1[j], fval, ierr, numfunc = res # optimizer
        V1[j] = -VF(a1[j], a, V0) # optimized value

    maxMet = np.max(np.abs(V1-V0)) # compute Chebyshev dist. between V1 and V0
    print(maxMet)
    iter += 1
    V0 = V1 # reassign VF to be iterated again if loop does not break

plt.plot(ES, V1)
plt.xlabel('<span class="math-container">$a$</span>')
plt.ylabel('<span class="math-container">$v(a)$</span>')
plt.show()
</code></pre>
<p><a href="https://i.sstatic.net/QOEVm9nZ.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/QOEVm9nZ.png" alt="enter image description here" /></a></p>
<p><a href="https://i.sstatic.net/Tp6D3GtJ.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Tp6D3GtJ.png" alt="enter image description here" /></a></p>

<computation><python>
Value function iteration savings problem not working as expected
manifold
# Value function iteration savings problem not working as expected

(Asked by: username_0 on 2024-05-14)

I'm coding a program in Python to find a numerical approximation of the value function $v:\mathbb{R}\to\mathbb{R}$ that solves the Bellman equation:
$$v(a)=\max\limits_{a'}\left\{\ln\left(w+(1+r)a-a'\right)+\beta v(a')\right\}$$
with $w, r>0$ and $\beta\in(0,1)$ given. When I code a value function iteration algorithm in Matlab to solve this problem I get the expected result a concave and smooth $v$ (1st image attached). Nevertheless, coding the exact same algorithm in Python using numpy and scipy I get weird results. In particular the algorithm 'converges' on the first iteration and the resulting VF is clearly too far from a good approximation (2nd image attached).
I would appreciate help in finding why the Python code does not achieve the same results as the Matlab while using the same algorithm. Thank you.
PD: My environment is Python 3.10.14, scipy 1.13.0.
Python code used:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fminbound

# Parameters
bbeta = 0.98
w = 1
r = 0.05
tol = 1e-10 # val. fun. distance tolerante
iters = 2000 # maximum number of iterations
N = 1000 # number of point in grid
amin = -6 # minimum level of assets in grid
amax = 30 # maximum level of assets in grid

ES = np.linspace(amin, amax, N) # state space

def VF(a, a0, V0):
    """
    This function returns the interpolated value function at a given a0
    """
    # interpolate value
    g = np.interp(a, ES, V0)
    # consumption
    c = w + (1+r)*a0 - a
    if c<=0:
        V = 1e10 + 1e10*np.abs(c)
        return V
    V = np.log(c) + bbeta*g
    return -V

V0 = np.zeros(N) # initial guess value function, this vector will be recycled
a1 = np.zeros(N) # vector that stores optimizers, i.e. approx. policy function
V1 = np.zeros(N) # vector that to store the iterated value function
iter = 1
maxMet = 100

while (iter < iters) and (maxMet>tol):
    for j in range(N):
        # fix the value of current point in the state space to iterate the VF
        a = ES[j]
        # apply the Bellman operator to V0
        res = fminbound(lambda x: VF(x, a, V0), amin, amax, xtol=1e-10,
                          full_output = True)
        a1[j], fval, ierr, numfunc = res # optimizer
        V1[j] = -VF(a1[j], a, V0) # optimized value

    maxMet = np.max(np.abs(V1-V0)) # compute Chebyshev dist. between V1 and V0
    print(maxMet)
    iter += 1
    V0 = V1 # reassign VF to be iterated again if loop does not break

plt.plot(ES, V1)
plt.xlabel('$a$')
plt.ylabel('$v(a)$')
plt.show()

**Tags:** computation, python


---

## 1 answers

### Answer by username_1 on 2024-05-14. Score: 2 | Accepted answer

In the last line of the outer loop, you are updating V0 by assigning V1 to it. However, doing so will make the maxMet calculation fail right after the first loop. It is better to copy only the values.
V0 = V1.copy()

Minor comments:

The value function should be written like this

def VF(a, a0, V0, bbeta, w, r):
    """
    This function returns the interpolated value function at a given a0.
    """
    # Interpolate value
    g = np.interp(a, ES, V0)
    # Consumption
    c = w + (1 + r) * a0 - a
    if c <= 0:
        V = -1e10 - 1e10 * abs(c)
    else:
        V = np.log(c) + bbeta * g
    return -V


Conversion seems very slow. I suggest making the tol around 1e-5. Since you are interpolating anyway, choosing N=1000 seems a bit excessive. I think 100 is enough.

> **Comments:**
> **username_0**: Thanks! Can you elaborate on the technicality of why this happens? What is the difference between both `V0=V1` and `V0=V1.copy()`? (thanks for the recommendations too)  
> **username_1**: After the first loop, if you directly assign `V0=V1`, both `V0` and `V1` now refer to the same array, which means any change to `V1` will now change `V0`. Then, in the second loop, when you calculate `V1[j]`, you are actually changing `V0` and `V1` simultaneously. This is why `maxMet` in the second loop is exactly 0, so the iteration stops!   By using the `.copy()` method, you make sure that `V0` remains unchanged during the next inner loop. Since the optimized value from the last iteration is preserved until the entire array of `V1` is computed, iteration continues. Hope it helps!  

---

