var c l k b w r tc e p z te y tb;

varexo x;

parameters mu r1 beta teta gama phi alpha pe prho rho sigma delta;

mu = 1/3;
delta = 0.028;
r1 = 0.044+delta;
beta = 1/(1+r1); 
teta = 1/3;
gama = 1/3;
phi = 0.15;
alpha = 0.4;
pe = 0.04;
prho = 0.09;
rho = 0.96;
sigma = 0.3;

model;


y = z*(k^(teta))*(l^gama)*(e^(1-teta-gama));  % eq 1

%  Shock de productividad

log(z) = rho*log(z(+1))+ x;   % eq 2

c + k(+1) - (1-delta)*k + tb + 0.15*y = y;      % eq 3

%  Shock de precios 

log(p) = pe + prho*log(p(+1))+ x;   % eq 4


tb = p*e + b(+1) - (1+r1)*b;      % eq 5


(c/(1-l)) = ((1-mu)/mu)*(w/(1+tc));   % eq 6 


((c/(1-l))^(-mu))*((1+phi*y)^(-alpha))/(1-tc) = beta*((((c(+1)/(1-l(+1)))^(-mu)))*((1+phi*y(+1))^(-alpha))*(1-delta+r(+1))/(1+tc(+1)));  % eq 7


r(+1) - delta = r1; % eq 8


(1+tc)*c + k(+1) - (1-delta)*k + b(+1) = w*l + r*k + (1+r1)*b;  % eq 9


w = gama*(y/l) ; % eq 10


r = teta*(y/k) ; % eq 11


(1+te)*p = (1-teta-gama)*(y/e); % eq 12


0.15*y = tc*c + te*p*e;      % eq 13

end;

initval;

c = 1.4;
l = 0.5;
k = 17.4;
e = 1.15;
z = 1;
y = z*(k^(teta))*(l^gama)*(e^(1-teta-gama));
b = 21;
w = gama*(y/l);
r = 0.04;
tc = 0.1 ;
p = 0.04;
te = (-tc*c)/(p*e)+0.15*y;
tb = -1;

end;

steady;

check;

shocks;

var x = sigma^2;

end;

stoch_simul(periods=1000);