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
    iter += 1
    V0 = V1 # reassign VF to be iterated again if loop does not break

plt.plot(ES, V1)

enter image description here

enter image description here


1 Answer 1


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

  • $\begingroup$ 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) $\endgroup$
    – manifold
    Commented May 14 at 1:48
  • 1
    $\begingroup$ 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! $\endgroup$
    – teddi
    Commented May 14 at 2:05

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.