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

```

[![enter image description here](https://i.sstatic.net/QOEVm9nZ.png)](https://i.sstatic.net/QOEVm9nZ.png)

[![enter image description here](https://i.sstatic.net/Tp6D3GtJ.png)](https://i.sstatic.net/Tp6D3GtJ.png)

**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.

[![enter image description here](https://i.sstatic.net/JpcfEN72.png)](https://i.sstatic.net/JpcfEN72.png)

> **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!  

---

