# El Modelo de Ciclos Económicos Reales con Agentes Ricardianos y No Ricardianos¶

In [ ]:
using NLsolve, LaTeXStrings, Plots


## El Modelo¶

### Familias Ricardianas¶

Las familias ricardianas maximizan:

$$\sum_{t=0}^{\infty} \beta^{t}\left[\gamma \log C_{i, t}+(1-\gamma) \log \left(1-L_{i, t}\right)\right]$$

sujeto a:

$$C_{i, t}+K_{i, t+1}=W_{t} L_{i, t}+\left(R_{t}+1-\delta\right) K_{i, t}$$

### Familias No Recardianas¶

Las familias no ricardianas maximizan:

$$\sum_{t=0}^{\infty} \beta^{t}\left[\gamma \log C_{j, t}+(1-\gamma) \log \left(1-L_{j, t}\right)\right]$$

sujeto a:

$$C_{j, t}=W_{t} L_{j, t}$$

### Empresas¶

Las empresas maximizan:

$$\max _{\left(K_{t}, L_{t}\right)} \Pi_{t}=A_{t} K_{t}^{\alpha} L_{t}^{1-\alpha}-R_{t} K_{t}-W_{t} L_{t}$$

con la tecnología que evoluciona de acuerdo a:

$$\ln A_{t}=\left(1-\rho_{A}\right) \ln \bar{A}+\rho_{A} \ln A_{t-1}+\varepsilon_{t}$$

### Agregación¶

Consumo:

$$C_{t}=\omega C_{i, t}+(1-\omega) C_{j, t}$$

$$L_{t}=\omega L_{i, t}+(1-\omega) L_{j, t}$$

\begin{aligned} K_{t} &=\omega K_{i, t} \\ I_{t} &=\omega I_{i, t} \end{aligned}

## Condiciones de Equilibrio¶

$$\begin{array}{c} \frac{1-\gamma}{\gamma} \frac{C_{i, t}}{1-L_{i, t}}=W_{t} \\ \frac{1-\gamma}{\gamma} \frac{C_{j, t}}{1-L_{j, t}}=W_{t} \\ \frac{C_{i, t+1}}{C_{i, t}}=\beta\left[R_{t+1}+1-\delta\right] \\ C_{j, t}=W_{t} L_{j, t} \\ C_{t}=\omega C_{i, t}+(1-\omega) C_{j, t} \\ L_{t}=\omega L_{i, t}+(1-\omega) L_{j, t} \\ Y_{t}=C_{t}+I_{t} \\ K_{t}=\omega K_{i, t} \\ I_{t}=\omega I_{i, t} \\ Y_{t}=A_{t} K_{t}^{\alpha} L_{t}^{1-\alpha} \\ W_{t}=(1-\alpha) Y_{t} / L_{t} \\ R_{t}=\alpha Y_{t} / K_{t} \\ K_{i, t+1}=(1-\delta) K_{i, t}+I_{i, t} \\ \ln A_{t}=\left(1-\rho_{A}\right) \ln \bar{A}+\rho_{A} \ln A_{t-1}+\varepsilon_{t} \end{array}$$
In [ ]:
# Parametros del modelo

α  = 0.35
β  = 0.97
γ  = 0.40
δ  = 0.06
ω  = 0.50
ρA = 0.95
σA = 0.01
Ā  = 1.0;
η  = [0, 1]

params = [α, β, γ, δ, ω, ρA, σA, Ā];


In [ ]:
# Sistema de estado estacionario

function sistema_ee(x,p)
y, ci, cj, c, i, ii, k, ki, l, li, lj, w, r = x
α, β, γ, δ, ω, ρA, σA, Ā = p

fval = similar(x)

fval[1]  = ((1-γ)/(γ))*(ci/(1-li)) - w
fval[2]  = ((1-γ)/(γ))*(cj/(1-lj)) - w
fval[3]  = 1 - β*(r + 1 - δ)
fval[4]  = cj - w*lj
fval[5]  = c - ω*ci - (1-ω)*cj
fval[6]  = l - ω*li - (1-ω)*lj
fval[7]  = y - c - i
fval[8]  = k - ω*ki
fval[9]  = i - ω*ii
fval[10] = y - Ā*(k^α)*(l^(1-α))
fval[11] = w - (1-α)*y/l
fval[12] = r - α*y/k
fval[13] = ki - ii - (1-δ)*ki

return fval

end

In [ ]:
x0 = [1.5, 0.6, 0.5, 0.5, 0.2, 0.2, 2.0, 5.0, 0.3, 0.3, 0.45, 1.0, 0.1]
solee = nlsolve(x -> sistema_ee(x, params), x0)
yee, ciee, cjee, cee, iee, iiee, kee, kiee, lee, liee, ljee, wee, ree =  solee.zero;

In [ ]:
solee.zero


## Método de Perturbación¶

In [ ]:
# Incluímos funciones necesarias para aplicar el método de perturbación
include("second_order_approx.jl");

In [ ]:
# Número de variables y shocks
nx = 2; # variables de estado ki, a
ny = 12; # variables de control k, y, ci, cj, c, i, ii, l, li, lj, w, r
ne = 1; # shock eps

# Definimos el modelo. Orden de variables
# orden: ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r
# las variables serán definidas en logs

function eq1(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return ((1-γ)/(γ))*(ci/(1-li)) - w
end

function eq2(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return ((1-γ)/(γ))*(cj/(1-lj)) - w
end

function eq3(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return cil/ci - β*(rl + 1 - δ)
end

function eq4(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return cj - w*lj
end

function eq5(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return c - ω*ci - (1-ω)*cj
end

function eq6(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return l - ω*li - (1-ω)*lj
end

function eq7(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return y - c - i
end

function eq8(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return k - ω*ki
end

function eq9(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return i - ω*ii
end

function eq10(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return y - a*(k^α)*(l^(1-α))
end

function eq11(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return w - (1-α)*y/l
end

function eq12(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return r - α*y/k
end

function eq13(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return kil - ii - (1-δ)*ki
end

function eq14(x)
kil, al, kl, yl, cil, cjl, cl, il, iil, ll, lil, ljl, wl, rl,
ki, a, k, y, ci, cj, c, i, ii, l, li, lj, w, r = exp.(x)

return al - (1-ρA)*log(Ā) - ρA*a
end

eqs = (eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12,eq13,eq14)

In [ ]:
xee = repeat([kiee, Ā, kee, yee, ciee, cjee, cee, iee, iiee, lee, liee, ljee, wee, ree],2)
lxee = log.(xee);

In [ ]:
# Calculando las matrices jacobiana y hesiana del sistema

In [ ]:
# Solución usando una aproximación de primer orden
Gs,Hs,Gx,Hx = solve_first_order_approx(lxee,J,nx,ny,ne);


## Análisis de Impulso Respuesta¶

In [ ]:
horizonte = 50
xx_a, yy_a = fir_first_order(Hx,Gx,horizonte,1,η,σA);

In [ ]:
nombre_vars = ["Capital" "PIB" "Consumo" "Trabajo"]

plot(1:horizonte, [xx_a[:,1] yy_a[:,2] yy_a[:,5] yy_a[:,8]],
xlabel = "horizonte", ylabel = "Desvío Porcentual del EE", label=nombre_vars,
title = "Funciones Impulso Respuesta: Asignaciones")

In [ ]:
plot(1:horizonte, [yy_a[:,3] yy_a[:,4]],
xlabel = "horizonte", ylabel = "Desvío Porcentual del EE", label=[L"C_R" L"C_{NR}"],
title = "Funciones Impulso Respuesta: Asignaciones")

In [ ]:
plot(1:horizonte, yy_a[:,11:12],
xlabel = "horizonte", ylabel = "Desvío Porcentual del EE", label=["Salarios" "Rental del Capital"],
title = "Funciones Impulso Respuesta: Precios")