# Value function iteration savings problem not working as expected

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


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


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

• 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) Commented May 14 at 1:48
• 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! Commented May 14 at 2:05