% Basic RBC Model 
%
% Jesus Fernandez-Villaverde
% Philadelphia, March 3, 2005

%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------

close all;

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

var y m c R n r pi P im ex q w b g t Pf cf Rf bf z; %20 variables 
varexo e1 e2 e3 e4 e5 e6 eta epsilon ;

parameters beta gamma delta rho1 rho2 rho3 rho4 rho5 theta omega1 omega2 alpha sigma1 sigma2 mu phi_t phi_g;

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

beta    = 0.97;
gamma   = 0.02;
delta   = 1.00;
rho1    = 0.95;  
rho2    = 0.95;
rho3    = 0.95;
rho4    = 0.95;
rho5    = 0.95;
theta   = 1.20;
omega1  = 0.7;
omega2  = 0.7;
alpha   = 0.8;
sigma1  = 2; 
sigma2  = 2;
mu      = 0.5;
phi_t   = -1;
phi_g   = 1;

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model; 
 
  1/R = beta*((c(+1)/c)^-theta)*(1+1/R(+1))*(P/P(+1)); %equation 1
  m^(-1)*c^(-theta)=mu*((1+1/R(+1)-1/R)/(1+1/R(+1))); %equation 2
  y = z*(n^alpha); %equation 3
  n = (alpha*y)/w; %equation 4 labour demand
  c+g+im-ex = y; %equation 5
  exp(im)=sigma1*exp(1-omega1)-sigma1*exp((omega1^sigma1)/(omega1^sigma1+(1-omega1)^sigma1))*exp(q)+ exp(c); %equation 6
  exp(ex)=sigma2*exp(1-omega2)+sigma2*exp((omega2^sigma1)/(omega2^sigma1+(1-omega2)^sigma1))*exp(q)+ exp(cf); %equation 7
  1-n = ((mu*c^(-theta)*exp(exp(w)-(omega1^sigma1/(omega1^sigma1+(1-omega1)^sigma1))*exp(q)+ exp(e6)))/(1-mu))^(-1/delta);
  Rf/R = ((1+1/Rf(+1))/(1+1/R))*q(+1)/q*Pf/Pf(+1); %equation 9
  R*y(-1)*(t-g)= b*(1+pi)*(r-gamma); %equation 10
  R = r+pi(+1); %equation 11
  b*(1+pi)/R(-1)=(g(-1)-t(-1))*y(-1)+b(-1)+(b(-1)/R(-1));  %equation 12
  pi = exp(P)-exp(P(-1));  %equation 13
  t = phi_t*t(-1)+ eta; %equation 14
  g = phi_g*g(-1)+ epsilon; %equation 15
  z = rho1*z(-1)+e1; %equation 16
  bf = rho2*bf(-1)+e2; %equation 17
  cf = rho3*cf(-1)+e3; %equation 18
  Rf = rho4*Rf(-1)+e4; %equation 19
  Pf = rho5*Pf(-1)+e5; %equation 20
  
end;

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

%y c R n r pi P im ex q w b g t Pf cf Rf bf z

initval;
  y = 1.08;
  c = 0.76;
  %m = 0.25;
  R = 0.05;
  n = 0.33;
  r = 0.01;
  pi= 0.56;
  P = 1;
  im= 0.3;
  ex= 0.31;
  q = 0.17; 
  w = 1.13;
  b = 20;
  g = 0.2;
  t = 1;
  Pf= 1;
  cf= 25;
  Rf= 0.05;
  bf= -0.82;
  z = 0; 
  e1 = 0;
  e2 = 0;
  e3 = 0;
  e4 = 0;
  e5 = 0;
  e6 = 0;
  eta = 0;
  epsilon = 0;
end;

shocks;
var e1 = 0.000045;
var e2 = 0.00003;
var e3 = 0.00005;
var e4 = 0.00004;
var e5 = 0.000036;
var eta = 0.000026;
var epsilon = 0.00002;
end;

steady;

check;

stoch_simul(hp_filter = 1600, order = 1);



%----------------------------------------------------------------
% 5. Some Results
%----------------------------------------------------------------

statistic1 = 100*sqrt(diag(oo_.var(1:6,1:6)))./oo_.mean(1:6);
table('Relative standard deviations in %',strvcat('VARIABLE','REL. S.D.'),lgy_(1:6,:),statistic1,10,8,4);