Mauricio M. Tejada
ILADES - Universidad Alberto Hurtado
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\} $$
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.
Algoritmo:
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)$.
# Discretizamos la variable de estado. Definimos primero el número de puntos en el grid de k
n = 500
500
Parametros del modelo:
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:
crit = 1e-6 # valor critico de convergencia
1.0e-6
Definimos el grid de puntos de $k$ como $k_1,...,k_n \in [0.6k^*,1.4k^*]$
# 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] $$
# Definir valores inciales
val0 = zeros(n) # Conjetura Inicial
valiter = val0; # Vector donde guardamos las iteraciones
using Optim
using Interpolations
Usamos interpolación lineal para evaluar la función entre los puntos definidos en el grid del capital.
v_interpolate = LinearInterpolation(kgrid, val0, extrapolation_bc = Line())
v_interpolate(kss)
0.0
Resolvemos usando el método de iteración de la función valor:
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
Para este problema particular conocemos la solución alebraica.
# 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
using Plots
gr()
Plots.GRBackend()
Función de política:
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])
Función valor:
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])
Proceso de iteración de la función valor:
plot(kgrid, valiter, title = "Iteración de la Función Valor",
xlabel = "k(t)", ylabel = "v(k(t))", legend = false, linewidth = 2,
grid = true)