Resolviendo el modelo de planificación de cass-koopmans con la iteración de función de valor en Python

Question

Editado: lo siento, utilicé stack exchange por primera vez, así que no sabía que podía escribir código python. ahora puedes copiar y ejecutar el código.

Estoy practicando resolviendo modelos económicos utilizando programación en python. Me refiero a https://python.quantecon.org/optgrowth.html y https://python.quantecon.org/optgrowth_fast.html.

Quiero resolver el problema de planificación de Cass-Koopmans con la ecuación de Bellman.

• Hay solo un agente que gobierna todas las asignaciones de recursos en la economía. Ella produce un (único) bien a través de una función de producción:

  $$ y_{t} = f(k_{t}) \tag{1} $$

• El bien puede ser usado para consumo e inversión.

• En este ejercicio, introducimos una depreciación parcial del capital, es decir, el capital se deprecia a la tasa $0 < \delta < 1$ cada período. Por lo tanto, la restricción de recursos para el agente es:

  $$ k_{t+1} + c_t = y_t + (1-\delta)k_t \tag{2} $$

• El agente quiere maximizar

  $$ \mathbb E \left[ \sum_{t = 0}^{\infty} \beta^t u(c_t) \right] \tag{3} $$ sujeto a (2), donde $ \beta \in (0, 1) $ es un factor de descuento.

• El modelo es casi idéntico al modelo de planificación Cass-Koopmans. Para resolver este modelo, aplicaremos el algoritmo de iteración de función de valor. Por lo tanto, primero reformulamos el problema de maximización en una ecuación de Bellman.

Podemos escribir la función de valor para el problema de maximización de utilidad del agente en forma de ecuación de Bellman. En este ejercicio utilizaremos el capital ($k$) como variable de estado en lugar de la producción ($y$). La función de valor se define de la siguiente manera:

$$ v(k) = \max_{0 \leq c \leq y + (1-\delta)k} \left\{ u(c) + \beta v(k') \right\} \tag{4} $$ sujeto a $$ k' = f(k) + (1-\delta)k - c $$

Esta formulación toma el consumo ($c$) como variable de control. Para los cálculos a continuación, utilizaremos el capital del próximo período ($k'$) como variable de control. Entonces (4) puede ser reescrito de la siguiente manera:

$$ v(k) = \max_{0 \leq k' \leq y + (1-\delta)k} \left\{ u\left(f(k) + (1-\delta)k - k'\right) + \beta v(k') \right\} \tag{5} $$

Básicamente convertimos el problema de maximización restringida en (4) en un problema de maximización irrestricta en (5).

Así que hice mi código en python:

import numpy as np
import matplotlib.pyplot as plt
from numba import njit, float64
from numba.experimental import jitclass
from quantecon.optimize.scalar_maximization import brent_max
from interpolation import interp

datos_crecimiento_opt = [('', float64),
                  ('', float64),
                  ('', float64),
                  ('', float64),
                  ('grid',float64[:])]

@jitclass(datos_crecimiento_opt)
class CrecimientoOptimo_VI:
    def __init__(self, =0.4, =0.96, =2.0, =0.1, grid_max=10, grid_size=500):
        self., self., self., self. = , , , 
        self.grid = np.linspace(0.1, grid_max, grid_size)

    def f(self, k):
        return k**self.

    def u(self, c):
        return c**(1 - self.) / (1 - self.)

    def objetivo(self, k, kp, v_array):
        f, u, ,  = self.f, self.u, self., self.
        v = lambda x: interp(self.grid, v_array, x)

        return u(f(k)+(1-)*k-kp) + *v(kp)

@njit
def T(v, og_VI):
    v_greedy = np.empty_like(v)
    v_new = np.empty_like(v)

    for i in range(len(og_VI.grid)):
        k = og_VI.grid[i]
        lower = 1e-10
        upper = og_VI.f(k) + (1-og_VI.)*k
        result = brent_max(og_VI.objetivo, lower, upper, args=(k,v))
        v_greedy[i], v_new[i] = result[0], result[1]

    return v_greedy, v_new

def resolver_modelo_VI(og_VI, tol=1e-4, max_iter=1000, print_skip=20):
    v = og_VI.grid
    error = tol+1
    i=0

    while i < max_iter and error > tol:
        v_greedy, v_new = T(v, og_VI)
        error = np.max(np.abs(v - v_new))
        i += 1
        if i % print_skip == 0:
            print(f"Error en la iteración {i} es {error}.")
        v = v_new

    if i == max_iter:
        print("¡No se pudo converger!")

    if i < max_iter:
        print(f"\nConvergido en {i} iteraciones.")

    return v_greedy, v_new

og_VI = CrecimientoOptimo_VI()
v_greedy, v_solucion = resolver_modelo_VI(og_VI)
plt.plot(og_VI.grid, v_greedy)

Pensé que hice mi código correctamente, pero produjo un gráfico de k óptimo extraño. ¿Podrías revisar mi código y decirme qué hay de malo en mi implementación? Es una molestia, pero cualquier ayuda sería muy apreciada.

Answer 1

Intenté el código que proporcionaste.

1. El problema parece ocurrir porque, en esos puntos, el consumo se vuelve negativo y se desordena con la función de utilidad. Entonces, en la definición de tu Clase OptimalGrowth_VI, agregué la siguiente pieza para el método objective:

   def objective(self, k, kp, v_array):
       f, u, ,  = self.f, self.u, self., self.
       v = lambda x: interp(self.grid, v_array, x)
       c = f(k)+(1-)*k-kp
       if c <= 0:
           u = -888 - 800 * abs(c)
       else:
           u = u(c)
       return u + *v(kp)

Básicamente, si c es negativo, devolverá una utilidad extremadamente pequeña, por lo que la función brent_max la ignorará. (Tratamiento de McCandless). Después de eso, tendrás la función de política curvada del libro de texto.

2. Otra cosa que noté es que la función de política, de acuerdo con la cuadrícula y los parámetros que definiste, no toca la línea de 45 grados. Tal vez la cuadrícula predefinida no contiene el valor de estado estacionario. Aumenté tu grid_max= 50 y agregué una línea de 45 grados para verificar la convergencia al estado estacionario.

   og_VI = OptimalGrowth_VI() v_greedy, v_solution = solve_model_VI(og_VI) plt.plot(og_VI.grid, v_greedy)

   agregar línea de 45 grados

   plt.plot(og_VI.grid, og_VI.grid, ls=':', label='línea de 45 grados') plt.xlabel('k(t)') plt.ylabel('k(t+1)')