%-------------------------------------------------------------------------
%Programa para obtener la dinámica del modelo de Solow con crecimiento de
%la tecnología
%-------------------------------------------------------------------------
%Borro todo antes
clear
clc
%Parámetros del modelo
A0=1;
a=0.01;
rho=0.9;
alpha=0.2;
n=0.005;
delta=0.1;
s=0.15;
kbar=((s)/(a+n+a*n+delta))^(1/(1-alpha));
k0=kbar;
k(1)=k0;
theta(1)=1;
%Loop para la dinámica del capital en el modelo empezando con k0=kbar
for j=1:40
    if j==2
    k(j+1)=(((1-delta)*k(j))/((1+n)*(1+a)))+((s*(theta(j)*k(j)^(alpha)))/((1+n)*(1+a)));
    y(j)=theta(j)*k(j)^(alpha);
    Y(j)=(1+a)^(j)*(1+n)^(j)*y(j);
    K(j)=(1+a)^(j)*(1+n)^(j)*k(j);
    epsilon(j)=0.2;
    theta(j+1)=theta(j)^(rho)*exp(epsilon(j));
    logtheta(j)=log(theta(j));
    else
    k(j+1)=(((1-delta)*k(j))/((1+n)*(1+a)))+((s*(theta(j)*k(j)^(alpha)))/((1+n)*(1+a)));
    y(j)=theta(j)*k(j)^(alpha);
    Y(j)=(1+a)^(j)*(1+n)^(j)*y(j);
    K(j)=(1+a)^(j)*(1+n)^(j)*k(j);
    epsilon(j)=0;
    theta(j+1)=theta(j)^(rho)*exp(epsilon(j));
    logtheta(j)=log(theta(j));
    end
end
subplot(3,1,1),plot(y,'g')
hold on
plot(k)
hold off
legend('y','k','\theta')
title('Impulso Respuesta')
subplot(3,1,2),plot(epsilon)
legend('\epsilon')
title('El shock inicial')
subplot(3,1,3),plot(logtheta)
legend('log(\theta)')
title('La variable del shock')