Блог

В процессе обучения я пишу задачу итерации функции значения двумя разными способами.

1. Весь код внутри файла python

 
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import rutinas as rt

# Tenemos que definir lo parametros 
alpha = .33
delta = .07
sigma = 1.5
gamma = 0     # Vamos a considerar que los hogares no valoran el ocio
beta  = .97
tfp   = 1.0

# Definicion de estados estaionarios
sss = alpha*delta/(1.0/beta-1.0 delta)
lss = (1.0-alpha)/(1.0-alpha gamma*(1-sss))
kss = (alpha*tfp/(1.0/beta-1.0 delta))**(alpha/(1.0-alpha))
yss = tfp*kss**alpha*lss**(1.0-alpha)
css = (1.0-sss)*yss

# Definimos un grid para la variable de estado de la economía
nk = 500

kmin = 1.0e-5
kmax = 4*kss
kgrid = np.linspace(kmin, kmax, nk)
# Definimos una matriz de consumos. Esto puede hacerse de dos formas:

# 1. Preasignando y llenado espacios
matrizC = np.empty((nk, nk))
matrizU = np.empty((nk, nk))
for i1 in range(0,nk):
    for i2 in range(0,nk):
        matrizC[i1, i2] = tfp*kgrid[i1]**alpha
                          (1.0-delta)*kgrid[i1]
                        - kgrid[i2]

# Un indicador de posiciones de consumo positivas
indC = matrizC <= 0
matrizU = matrizC**(1.0-sigma)/(1.0-sigma)
matrizU[indC] = -1.e 6

# A partir de aquí podemos inicializar el algoritmo de programacion dinamica
valorInit = np.zeros(nk)
valorObj  = np.zeros((nk, nk))

# Criterio de convergencia del algoritmo
tol = 1e-6

# Inicializamos el bucle while y un error arbitrario inicial
cont = 0
error = 10

while error > tol:
    for i2 in range(0, nk):
        valorObj[:, i2] = matrizU[:, i2]   beta*valorInit[i2]
    pos   = valorObj.argmax(axis = 1) 
    valor = valorObj.max(axis = 1) 
    error = np.max(np.abs(valor-valorInit))
    valorInit = valor
    cont  = 1
    print(error)

  

2. Внутри модуля с другими подпрограммами

• Алгоритм внутри подпрограммы представляет собой копирование-вставку предыдущего файла

 def mnc_bellman(param):
    import numpy as np
    # Declaramos los argumentos del problema de programacion dinamica
    # tol, kmin, kmax, tam, valorInit = args
    alpha, delta, sigma, gamma, beta, tfp = param

    kss = (alpha * tfp / (1.0 / beta - 1.0   delta))**(alpha / (1.0 - alpha))
    tol = 1.0e-6
    kmin = 1.0e-4
    kmax = 4*kss
    tam = 500
    valorInit = np.zeros((tam))
    # Inicializamos objetos que vamos a necesitar
    matrizC = np.empty((tam, tam))
    valorObj = np.zeros((tam, tam))
    kgrid = np.linspace(kmin, kmax, tam)

    for i1 in range(0, tam):
        for i2 in range(0, tam):
            matrizC[i1, i2] = tfp * kgrid[i1] ** alpha 
                                (1.0 - delta) * kgrid[i1] 
                              - kgrid[i2]

    # Un indicador de posiciones de consumo positivas
    indC = matrizC <= 0
    matrizU = matrizC ** (1.0 - sigma) / (1.0 - sigma)
    matrizU[indC] = -1.e 6

    # A partir de aqui podemos inicializar el algoritmo ITERACION de la
    # funcion de VALOR

    # Esta variable va a guardar cuantas veces iteramos hasta cumplir con los
    # criterios de convergencia
    cont = 0

    # Tambien tenemos que establecer un error inicial que no se satisfaga
    # trivialmente en el primer paso del bucle
    error = 10

    while error > tol:
        for i2 in range(0, tam):
            valorObj[:, i2] = matrizU[:, i2]   beta * valorInit[i2]
            pos = valorObj.argmax(axis=1)
            valor = valorObj.max(axis=1)
            error = np.max(np.abs(valor - valorInit))
            valorInit = valor
            cont  = 1
        print(error)

    # Ahora reconstruimos las reglas de decision
    reglaK = kgrid[pos]
    reglaC = tfp * kgrid ** alpha   (1.0 - delta) * kgrid - reglaK

    return valor, reglaK, reglaC

  

Который я позже вызываю из основной программы

 import numpy as np
import numpy.matlib
#import matplotlib.pyplot as plt
#from scipy.optimize import fsolve
import rutinas as rt

# Definicion de parametros del modelo (posteriormente intentare meterlos en un
# diccionario)
alpha = .36
delta = .08
sigma = 1.5
gamma = 0.5
beta = .97
tfp = 1.0

# Definicion de estado estacionario correspondiente a los parametros del
# modelo
sss = alpha * delta / (1.0 / beta - 1.0   delta)
lss = (1.0 - alpha) / (1.0 - alpha   gamma * (1.0 - sss))
kss = (alpha * tfp / (1.0 / beta - 1.0   delta))**(alpha / (1.0 - alpha))
yss = tfp * kss**alpha * lss**(1.0 - alpha)
css = (1.0 - sss) * yss

# Parametros relacionados con el algoritmo de iteracion de la funcion de valor
tol = 1.0e-6
kmin = 1.0e-4
kmax = 3.0*kss
tam = 500
valorInit = np.zeros((tam))

param = [alpha, delta, sigma, gamma, beta, tfp]

valor, reglaK, reglaC = rt.mnc_bellman(param)
  

Оба кода выполняются и выдают один и тот же результат. Однако второй вариант на порядки медленнее.

Я надеюсь, что кто-нибудь может указать мне, как решить это несоответствие. Я думал об использовании Jit, но сначала я хочу понять проблему.

Спасибо