La iteración de la función de valor en el problema de ahorro no está funcionando como se esperaba

Question

Estoy codificando un programa en Python para encontrar una aproximación numérica de la función de valor $v:\mathbb{R}\to\mathbb{R}$ que resuelve la ecuación de Bellman:

$$v(a)=\max\limits_{a'}\left\{\ln\left(w+(1+r)a-a'\right)+\beta v(a')\right\}$$

con $w, r>0$ y $\beta\in(0,1)$ dados. Cuando codifico un algoritmo de iteración de función de valor en Matlab para resolver este problema, obtengo el resultado esperado, una función $v$ concava y suave (1ª imagen adjunta). Sin embargo, al codificar el exactamente mismo algoritmo en Python usando numpy y scipy, obtengo resultados extraños. En particular, el algoritmo 'converge' en la primera iteración y la VF resultante claramente se distancia demasiado de una buena aproximación (2ª imagen adjunta).

Agradecería ayuda para encontrar por qué el código de Python no logra los mismos resultados que Matlab utilizando el mismo algoritmo. Gracias.

PD: Mi entorno es Python 3.10.14, scipy 1.13.0.

Código de Python utilizado:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fminbound

# Parámetros
bbeta = 0.98
w = 1
r = 0.05
tol = 1e-10 # tolerancia de distancia de la función de valor
iters = 2000 # número máximo de iteraciones
N = 1000 # número de puntos en la cuadrícula
amin = -6 # nivel mínimo de activos en la cuadrícula
amax = 30 # nivel máximo de activos en la cuadrícula

ES = np.linspace(amin, amax, N) # espacio de estado

def VF(a, a0, V0):
    """
    Esta función devuelve la función de valor interpolada en un dado a0
    """
    # interpolar valor
    g = np.interp(a, ES, V0)
    # consumo
    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) # suposición inicial de la función de valor, este vector se reciclará
a1 = np.zeros(N) # vector que almacena los optimizadores, es decir, la función de política aproximada
V1 = np.zeros(N) # vector para almacenar la función de valor iterada
iter = 1
maxMet = 100

while (iter < iters) and (maxMet>tol):
    for j in range(N):
        # fijar el valor del punto actual en el espacio de estado para iterar la VF
        a = ES[j]
        # aplicar el operador de Bellman a V0
        res = fminbound(lambda x: VF(x, a, V0), amin, amax, xtol=1e-10,
                          full_output = True)
        a1[j], fval, ierr, numfunc = res # optimizador
        V1[j] = -VF(a1[j], a, V0) # valor optimizado

    maxMet = np.max(np.abs(V1-V0)) # calcular la distancia de Chebyshev entre V1 y V0
    print(maxMet)
    iter += 1
    V0 = V1 # reasignar VF para ser iterada nuevamente si el bucle no se rompe

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

Answer 1

En la última línea del bucle externo, estás actualizando V0 asignándole V1. Sin embargo, al hacerlo, el cálculo de maxMet fallará justo después del primer bucle. Es mejor copiar solo los valores.

V0 = V1.copy()

Comentarios menores:

• La función de valor debería estar escrita de esta manera

  def VF(a, a0, V0, bbeta, w, r): """ Esta función devuelve la función de valor interpolada en un determinado a0. """

  Interpolación de valor

  g = np.interp(a, ES, V0)
  # Consumo
  c = w + (1 + r) * a0 - a
  if c <= 0:
      V = -1e10 - 1e10 * abs(c)
  else:
      V = np.log(c) + bbeta * g
  return -V

• La conversión parece ser muy lenta. Sugiero hacer que la tol sea alrededor de 1e-5. Dado que estás interpolando de todos modos, elegir N=1000 parece un poco excesivo. Creo que 100 es suficiente.