%-------------------------------------------------------------------------
%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.02;
theta=0.36;
n=0.05;
delta=0.1;
s=0.25;
k0=0.5;
k(1)=k0;
v0=5;
v(1)=v0;
h0=1.8102;
h(1)=h0;
%Calculo el estado estacionario para luego verificar
kbar=((s)/(a+n+a*n+delta))^(1/(1-theta))
%Loop para la dinámica del capital en el modelo empezando con k0=0.5
for j=1:50
    k(j+1)=(((1-delta)*k(j))/((1+n)*(1+a)))+((s*(k(j)^(theta)))/((1+n)*(1+a)));
    %Quiero ver acá también como es la tasa de crecimiento del capital
    %agregado a la que voy a llamar gammaK
    gammaK(j)=(k(j+1)*(1+a)*(1+n))/(k(j));
    %Acá planteamos la ecuación de la tasa de crecimiento del salario salario y su evolución en el tiempo en
    %el modelo
    w(j)=(1-theta)*((1+a)^(j))*A0*(k(j)^(theta));
end
%Loop para la dinámica del capital en el modelo empezando con k0=5
for j=1:50
    v(j+1)=(((1-delta)*v(j))/((1+n)*(1+a)))+((s*(v(j)^(theta)))/((1+n)*(1+a)));
end
%Loop para la dinámica del capital en el modelo empezando con k0=1.8102
for j=1:50
    h(j+1)=(((1-delta)*h(j))/((1+n)*(1+a)))+((s*(h(j)^(theta)))/((1+n)*(1+a)));
end
%Ahora quiero usar otra forma de ver el estado estacionario
%Recordemos que (delta+n)k=sAk^(theta)
%Genero un vector para k que voy a nombrar x
x=[0:0.01:4];
e=[1:1:50];
%Genero un vector y que va a reportarme el lado izquierdo de la ecuación
y=((1+a)*(1+n)-(1-delta))*x;
%Genero un vector z que me va a reportar el lado derecho
z=(s*(x.^(theta)));
%Ahora queremos ver la tasa de crecimiento de capital per cápita efectivo y hacer estática
%comparativa. Para eso usamos el vector x que habíamos usado antes.
%Fijarse que gamma hace referencia al vector de la tasa de crecimiento del
%capital y vamos a ver como se mueve ante diferentes niveles de capital per
%cápita efectivo
gamma=(((1-delta))/((1+n)*(1+a)))+((s*(x.^(theta-1)))./((1+n)*(1+a)));
%Gráfico de la autopista
subplot(2,2,1), plot(k,'b')
hold on
plot(v,'r')
hold on
plot(h,'g')
title('Turnpike')
hold off
%Gráfico del equilibrio para llegar al estado estacionario
subplot(2,2,2), plot(x,y,'r')
hold on
plot(x,z,'g')
legend('[(1+a)(1+n)-(1-\delta)]k','sk^{\theta}')
title('Estado Estacionario')
hold off
%Gráfico de la estática comparativa para el capital per cápita efectivo
subplot(2,2,3), plot(x,gamma)
legend('\gamma')
title('Estática comparativa')
%Gráfico de la tasa de crecimiento de capital agregado de la economía y del
%salario en la economía
subplot(2,2,4), plot(gammaK,'b')
hold on
plot(w,'r')
legend('\gamma_{K}','w')
title('Tasa de crecimiento de capital y salarios')
hold off