Macroeconomía I

Método de Perturbación basado en la Aproximación de las Funciones de Política

(El Modelo Neoclásico de Crecimiento)

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}\left[\frac{c_{t}^{1-\sigma}-1}{1-\sigma}\right]\\ s.a & & c_{t}+i_{t}=Ak_{t}^{\alpha}\\ & & k_{t+1}=i_{t}+(1-\delta)k_{t}\\ & & k_{0}\,dado. \end{eqnarray*}

Alternativamente:

\begin{eqnarray*} \max U_{0} & = & \sum_{t=0}^{\infty}\beta^{t}\left[\frac{c_{t}^{1-\sigma}-1}{1-\sigma}\right]\\ s.a & & c_{t}+k_{t+1}=Ak_{t}^{\alpha}+(1-\delta)k_{t}\\ & & k_{0}\,dado. \end{eqnarray*}

De las CPO tenemos el siguiente sistema de ecuaciones en diferencias (no lineal):

\begin{eqnarray*} c_{t}^{-\sigma} & = & \beta c_{t+1}^{-\sigma}\left[\alpha Ak_{t+1}^{\alpha-1}+(1-\delta)\right]\\ c_{t}+k_{t+1} & = & Ak_{t}^{\alpha}+(1-\delta)k_{t} \end{eqnarray*}

Este sistema se puede escribir como:

$$F(c_{t+1}, c_{t}, k_{t+1}, k_{t}) = 0$$

Estado Estacionario

El estado estacionario resuelve: \begin{eqnarray*} 1 & = & \beta\left[\alpha Ak^{*\alpha-1}+(1-\delta)\right]\\ c^{*} & = & Ak^{*\alpha}-\delta k^{*} \end{eqnarray*}

La solución puede ser hallada algebraicamente:

\begin{eqnarray*} k^{*} & = & \left[\frac{1-\beta(1-\delta)}{\beta\alpha A}\right]^{\frac{1}{\alpha-1}}\\ c^{*} & = & A\left[\frac{1-\beta(1-\delta)}{\beta\alpha A}\right]^{\frac{\alpha}{\alpha-1}}-\delta\left[\frac{1-\beta(1-\delta)}{\beta\alpha A}\right]^{\frac{1}{\alpha-1}} \end{eqnarray*}

In [1]:
# Parametrización

α = 0.3
β = 0.9
σ = 0.5
A = 2
δ = 0.25

params = (α, β, σ, A, δ);
In [2]:
# Definimos el número de variables de estado y de control

n = 1  # Variables de estado
m = 1;  # Variables de control

Calculamos el estado estacionario usando un solver de julia.

In [3]:
# Definimos primero es sistema de ecuaciones a resolver.

function systemss(x,pr)
    # Nota: los parámetros dentro la función son considerados variables locales 
    α, β, σ, A, δ = pr
    
    ks = x[1]
    cs = x[2]
    
    fout = zeros(2)

    fout[1] = ks - ((1-β*(1-δ))/(α*β*A))^(1/(α-1))
    fout[2] = cs - A*(ks^α) + δ*ks
    
    return fout
end
Out[3]:
systemss (generic function with 1 method)
In [5]:
# Pkg.add("NLsolve")
using NLsolve

xss = nlsolve(x -> systemss(x,params), [0.5; 0.5], inplace=false)

ks = xss.zero[1]
cs = xss.zero[2]
ys = A*(ks^α)
is = ks - (1-δ)*ks

println("ks = $ks")
println("cs = $cs")
println("ys = $ys")
println("is = $is")
ks = 2.065450805481485
cs = 1.9698280830054893
ys = 2.486190784375861
is = 0.5163627013703711

Aproximación Lineal de las Funciones de Política

Recordemos que el sistema de ecuaciones que caracteriza el óptimo era:

$$F(c_{t+1}, c_{t}, k_{t+1}, k_{t}) = 0$$

La solución de este sistema, son dos funciones (denominadas funciones de política) que permiten recuperar la secuencia de consumo y capital:

\begin{eqnarray*} c_t & = & g(k_t)\\ k_{t+1} & = & h(k_t) \end{eqnarray*}

Ahora podemos aproximar linealmente, alrededor del estado estacionario, las funciones de política:

\begin{eqnarray*} c_{t} &=& g(k^*) + g_k(k^*)(k_t - k^*) \\ k_{t+1} &=& h(k^*) + h_k(k^*)(k_t - k^*) \end{eqnarray*}

El problema radica entonces en hallar las derivadas $g_k(k^*)$ y $h_k(k^*)$.

El método general:

El sistema de CPOs es:

$$F(y_{t+1},y_{t},x_{t+1},x_{t}) = 0$$

donde $y$ es el vector de variables de control (de tamaño $m \times 1$) y $x$ es el vector de variables de estado (de tamaño $n \times 1$).

La aproximación de primer orden de las funciones de política es:

\begin{eqnarray*} g(x) &=& g(x^*) + g_x(x^*)(x-x^*) \\ h(x) &=& h(x^*) + h_x(x^*)(x-x^*) \end{eqnarray*}

Por definición de función de política, en estado estacionario tenemos que: \begin{eqnarray*} g(x^*) &=& y^* \\ h(x^*) &=& x^* \end{eqnarray*}

Usando el sistema descrito en $F(\cdot)$ y las funciones de política:

$$F(g(h(x_{t})),g(x_{t}),h(x_{t}),x_{t}) = F(x_t) = 0$$

Derivando el sistema de ecuaciones de las CPOs respecto de $x$ y evaluando en el estado estacionario tenemos:

$$F_x(x^*) = f_{y'} g_x h_x + f_y g_x + f_{x'}h_x + f_x = 0$$

En notación matricial:

\begin{align*} \left[\begin{array}{cc} f_{x'} & f_{y'}\end{array}\right]\left[\begin{array}{c} I\\ g_{x} \end{array}\right]h_{x} & =-\left[\begin{array}{cc} f_{x} & f_{y}\end{array}\right]\left[\begin{array}{c} I\\ g_{x} \end{array}\right]\\ A\left[\begin{array}{c} I\\ g_{x} \end{array}\right]h_{x} & =B\left[\begin{array}{c} I\\ g_{x} \end{array}\right] \end{align*}

Por otro lado, sea $P$ la matrix de autovectores de la matriz $h_x$ de tal manera que:

$$h_x P = P \Lambda$$

Si definimos que $Z=\left[\begin{array}{c} I\\ g_{x} \end{array}\right]P$ entonces la ecuación anterior se puede escribir como:

$$A Z \Lambda = B Z$$

Este problema se puede expresar como un problema de autovalores generalizado. Para matrices $A$ y $B$ dadas, existe una matriz V y una matriz diagonal $D$ tal que:

$$A\left[\begin{array}{cc} V_{1} & V_{2}\end{array}\right]\left[\begin{array}{cc} D_{11} & 0\\ 0 & D_{22} \end{array}\right]=B\left[\begin{array}{cc} V_{1} & V_{2}\end{array}\right]$$

Si todos los autovalores en $D_{11}$ están dentro el circulo unitario, tenemos que:

\begin{align*} \Lambda & =D_{11}\\ \left[\begin{array}{c} I\\ g_{x} \end{array}\right]P & =V_{1}=\left[\begin{array}{c} V_{11}\\ V_{12} \end{array}\right] \end{align*}

con lo cual:

\begin{align*} h_{x} & =V_{11}D_{11}V_{11}^{-1}\\ g_{x} & =V_{12}V_{11}^{-1} \end{align*}

Implementación

Definimos primero el sistema de ecuaciones en diferencia caracterizado por las CPOs del problema:

In [6]:
function systemcpo(x, pr)
    α, β, σ, A, δ = pr

    kf = x[1]
    cf = x[2]
    k  = x[3]
    c  = x[4]
    
    fout = similar(x, 2)

    fout[1] = c^-σ - β*(cf^-σ)*(α*A*(kf^(α-1))+(1-δ))
    fout[2] = c + kf - A*(k^α) - (1-δ)*k
    
    return fout
end
Out[6]:
systemcpo (generic function with 1 method)

Verificamos si el sistema asociado a las CPOs esta bien escrito. En estado estacionario debiera cumplirse exáctamente con cero.

In [7]:
systemcpo([ks; cs; ks; cs], params)
Out[7]:
2-element Array{Float64,1}:
 0.0
 0.0

Usamos el paquete ForwardDiff para calcular la matriz jacobiana del sistema (esto es $f_{k'}$, $f_{k}$, $f_{c'}$ y $f_{c}$) evaluadas en el estado estacionario:

In [8]:
# Pkg.add("ForwardDiff")
using ForwardDiff

jac_systemcpo(y) = ForwardDiff.jacobian(x -> systemcpo(x,params), y)
xvars = [ks; cs; ks; cs]
F = jac_systemcpo(xvars)
Out[8]:
2×4 Array{Float64,2}:
 0.0784788  0.180854   0.0      -0.180854
 1.0        0.0       -1.11111   1.0     

A partir de este punto el procedimiento es estándar. Por tanto, el código puede ser reutilizado.

In [29]:
using LinearAlgebra

A = F[:,1:n+m]     
B = -F[:,n+m+1:end] 

D = eigvals(B,A);
V = eigvecs(B,A);

DD = Diagonal(D)
D11 = DD[1:n,1:n]
D22 = DD[n+1:n+m, n+1:n+m]

V1 = V[:,1:n]
V11 = V1[1:n,:]
V12 = V1[n+1:n+m,:];

Las funciones de política se calculan como sigue:

In [30]:
hx = V11*D11*inv(V11)
gx = V12*inv(V11)

println("Solucion:")
println("c(t) = $(gx[1]) * k(t)")
println("k(t+1) = $(hx[1]) * k(t)")
Solucion:
c(t) = 0.5514722813918117 * k(t)
k(t+1) = 0.5596388297192996 * k(t)

Simulando Transiciones al Estado Estacionario

In [31]:
k0 = 0.2*ks
println("Stock de Capital Inicial es $k0")
Stock de Capital Inicial es 0.413090161096297
In [32]:
using Plots
gr()