% modelo de chumacero y fuente
%
% con 2 schocks z y p
% 

%----------------------------------------------------------------
% 0. Housekeeping
%----------------------------------------------------------------

%close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

var y c k i l  z1 z2 z3 z4;
varexo e1 e2 e3 e4;

parameters beta psi delta alpha rho1 rho2 rho3 rho4 sigma1 sigma2 sigma3 sigma4 yexp;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

alpha   = 1/3;
beta    = 0.98;
delta   = 0.06;
psi     =0.27; %0.27 expk 2.65
yexp    =2; %mejor valor 2 
rho1    = 0.73; 
rho2    = 0.844; 
rho3    = 0.895;
rho4    = 0.892;
sigma1  = 0.04;
sigma2  = 0.1;
sigma3   = 0.024;
sigma4  = 0.14;
%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model; 
  (1/exp(c)) = beta*(1/exp(c(+1)))*(1/exp(z2))*((alpha*(exp(k)^(alpha-1))*exp(exp(z1(+1)))*(exp(l(+1))^(1-alpha)))+((1-delta)*exp(z2(+1))));
  ((1-psi)/psi)*exp(c)/(1-exp(l)) = (1-alpha)*(exp(k(-1))^alpha)*exp(exp(z1))*(exp(l)^(-alpha));
  exp(c)+exp(i)+exp(z3)= exp(y) + exp(z4)*(yexp);
  exp(y) = (exp(k(-1))^alpha)*exp(exp(z1))*(exp(l)^(1-alpha));
  exp(k)*exp(z2) = (1-delta)*exp(k(-1))*exp(z2)+exp((i));
  z1 = rho1*z1(-1)+e1;
  z2 = rho2*z2(-1)+e2;
  z3= rho3*z3(-1)+e3; 
  z4= rho4*z4(-1)+e4;
end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

initval;
  k  = 4;
  c  = 0.76;
  l  = 0.3;
  z1 = 0; 
  z2 = 0;
  z3=  0;
  z4= 0;
  
end;

shocks;
var e1 = sigma1^2;
var e2 = sigma2^2;
var e3 = sigma3^2;
var e4 = sigma4^2;
end;

steady;
%h = waitbar(0,'Por favor espere');
%for tt=0:100
%waitbar(tt/120)
%end;
stoch_simul(irf=72 , order = 1, periods=5000);
%stoch_simul(hp_filter = 1600, irf=72 , order = 1);
%graficos
%----------------------------------------------------------------
% 5. Some Results
%----------------------------------------------------------------

%waitbar(1)

%statistic1 = 100*sqrt(diag(oo_.var(1:6,1:6)))./oo_.mean(1:6);
graficos1
graficosarr
open('C:\dynare\4.0.3\matlab\rbc_investmentplog.log')

%table('Relative standard deviations in %',strvcat('VARIABLE','REL. S.D.'),lgy_(1:6,:),statistic1,10,8,4)