Solución del Modelo Neoclásico de Crecimiento usando Programación Dinámica e Iteración de la Función Valor

Mauricio M. Tejada

ILADES - Universidad Alberto Hurtado


El Problema de Opmización

El problema de optimización del planificador central es el siguiente: \begin{eqnarray*} \max U_{0} & = & \sum_{t=0}^{\infty}\beta^{t}\ln c_{t}\\ s.a & & c_{t}+i_{t}=Ak_{t}^{\alpha}\\ & & k_{t+1}=i_{t}\\ & & k_{0}\,dado. \end{eqnarray*}

Alternativamente: \begin{eqnarray*} \max U_{0} & = & \sum_{t=0}^{\infty}\beta^{t}\ln c_{t}\\ s.a & & c_{t}+k_{t+1}=Ak_{t}^{\alpha}\\ & & k_{0}\,dado. \end{eqnarray*}

La ecuación de Bellman es: $$ v(k_{t})=\max_{k_{t+1}}\left\{ \log\left(Ak_{t}^{\alpha}-k_{t+1}\right)+\beta v(k_{t+1})\right\} $$

Solución Algebraica

La solución del problema funcional está dada por la función valor y la función de política:

$$ v(k_{k})\,\,y\,\,k_{t+1}=g(k_{t}) $$

Conocemos la solución de este problema:

\begin{eqnarray*} k_{t+1} & = & \frac{\beta F}{1+\beta F}Ak_{t}^{\alpha}\\ v(k_{t}) & = & E+F\ln(k_{t}) \end{eqnarray*}

Donde: \begin{eqnarray*} E & = & \frac{1}{1-\beta}\left[\ln\left(A\left(1-\alpha\beta\right)\right)+\frac{\alpha\beta}{1-\alpha\beta}\ln\left(Aa\beta\right)\right]\\ F & = & \frac{\alpha}{1-\alpha\beta} \end{eqnarray*}

La idea es implementar la solución numérica en Julia y compararla con la obtenida algebraicamente.

Implementación

Algoritmo:

  1. Definir un grid de puntos $k$ (discretización) y realizar una conjetura para $v(k)$ en cada posible punto del grid.
  2. Usando la conjetura evaluar el operador $Tv$ para cada punto en el grid (esto implica resolver el $\max$).
  3. Evaluar si $v=Tv$ (usando algún grado de tolerancia), en tal caso terminar el proceso. Caso contrario volver al punto 2 usando la función valor resultante.

En términos de la ecuación de Bellman, la iteración sería: $$ v_{j+1}(k)=\max_{k'}\left\{ \ln\left(Ak^{\alpha}-k'\right)+\beta v_{j}(k')\right\} $$

Partir con algún $v_{0}(k)$.

In [1]:
# Discretizamos la variable de estado. Definimos primero el número de puntos en el grid de k
n = 500
Out[1]:
500

Parametros del modelo:

In [2]:
A = 1      # tecnologia
α = 0.36   #  participacion del capital   
β = 0.9;   #  factor de descuento

Definimos el grado de tolerancia de la iteración de la función valor:

In [3]:
crit = 1e-6   # valor critico de convergencia
Out[3]:
1.0e-6

Definimos el grid de puntos de $k$ como $k_1,...,k_n \in [0.6k^*,1.4k^*]$

In [4]:
# Calcular kss y definir kgrid  
kss   = (α*β)^(1/(1-α))                       # capital en EE    
kgrid = kss * collect(range(0.6,stop=1.4,length=n));  # n x 1  
#println(kgrid)

Note que mientras más puntos tengamos en el grid, más precisión tendremos.

El siguiente paso es iterar la función valor partiendo de la siguiente conjetura: $$ v_{0}(k)=\left[\begin{array}{c} 0\\ \vdots \\ 0 \end{array}\right] $$

In [5]:
# Definir valores inciales

val0 = zeros(n)  # Conjetura Inicial  
valiter = val0;  # Vector donde guardamos las iteraciones
In [14]:
using Optim
using Interpolations

Usamos interpolación lineal para evaluar la función entre los puntos definidos en el grid del capital.

In [15]:
v_interpolate = LinearInterpolation(kgrid, val0, extrapolation_bc = Line())
v_interpolate(kss)
Out[15]:
0.0

Resolvemos usando el método de iteración de la función valor:

In [16]:
diff = 1

val1  = zeros(n)  # Vector para guardar la función valor
kdeci = zeros(n)  # Vector para guardar la función de política 

while diff > crit
    
    global val0, val1, valiter, kdeci, diff
    
    v_interpolate = LinearInterpolation(kgrid, val0, extrapolation_bc = Line())                      
    Tv(k,kf) = log(A*k^α-kf) + β*v_interpolate(kf) 
    
    for i in 1:n
        res = optimize(x -> -Tv(kgrid[i],x), 1e-10, A*kgrid[i]^α)
        kdeci[i] = res.minimizer
        val1[i] = -res.minimum
    end
    
    diff = abs(maximum(val1-val0))      
    #println(diff)
    valiter = [valiter val1]
    val0 = copy(val1) 
end

println("Convergencia = $diff");
Convergencia = 9.774898082071104e-7

Solución Algebraica

Para este problema particular conocemos la solución alebraica.

In [17]:
# Computar los coeficiente E y F 
E1 = log(A*(1 - α*β))
E2 = α*β/(1 - α*β) * log(A*α*β)
E = 1/(1 - β)*(E1 + E2)
F = α/(1 - α*β)

# Evaluar en la solucion algebraica la FV y la FP en kgrid
val_cs = E*ones(n) + F * log.(kgrid)            # n x 1
kdeci_cs = β*F/(1 + β*F)*A*kgrid.^α;            # n x 1

Gráficos

In [10]:
using Plots
gr()
Out[10]:
Plots.GRBackend()

Función de política:

In [18]:
plot(kgrid, [kdeci kdeci_cs], title = "Función de Política", 
     xlabel = "k(t)", ylabel = "k(t+1)", label=["Sol. Num." "Sol. Alg."], 
     linewidth = 2, grid = true, line=[:solid :dot], color=[:red :blue])
Out[18]:

Función valor:

In [19]:
plot(kgrid, [val1 val_cs], title = "Función Valor", 
     xlabel = "k(t)", ylabel = "v(k(t))", label=["Sol. Num." "Sol. Alg."], 
     linewidth = 2, grid = true, line=[:solid :dot], color=[:red :blue])
Out[19]:

Proceso de iteración de la función valor:

In [20]:
plot(kgrid, valiter, title = "Iteración de la Función Valor", 
     xlabel = "k(t)", ylabel = "v(k(t))", legend = false, linewidth = 2, 
     grid = true)
Out[20]: