% RBC Model with an underground sector. 
%
% Busato and Chiarini Model
%
% I borrow their calibration.
%
% Catalina Granda-Carvajal
% Storrs, July 2010.

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

close all;

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

var c k l i r w g y_m y_u y_tot z_m z_u tau_m tau_y;
varexo e_m e_u e_tm e_ty;

parameters beta h f eta gamma phi psi alpha delta rho_1 rho_2 rho_3 rho_4 sigmae_m sigmae_u sigmae_tm sigmae_ty;

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

% Parameters taken directly from BC
beta  = 0.98;         
alpha = 0.33;
h     = 0.55;
f     = 1.99;    
eta   = 0.40;
gamma = 3.00;
phi   = 0.03;
psi   = 1.30;
delta = 0.025;
rho_1 = 0.90;
rho_2 = 0.90;
rho_3 = 0.90;
rho_4 = 0.90;
sigmae_m = 0.003;
sigmae_u = 0.003;
sigmae_tm = 0.003;
sigmae_ty = 0.003;

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

model; 
  1/c = beta*(1/c(+1))*(1+r(+1)-delta);
  (w*tau_y)/c = -h*(l^gamma)+h*((2+gamma)/(1+gamma))*(l^(1+gamma))+f*((1-l)^(-gamma));
  r = alpha*(1-tau_m)*z_m*(k^(alpha-1))*(l^(1-alpha));
  w*tau_m = (1-tau_m)*(1-alpha)*z_m*(k^(alpha))*(l^(-alpha))-(1-phi*psi*tau_m)*z_u;
  g = w*tau_y*l+phi*psi*tau_m*y_u+tau_m*y_m;
  c+i+g = y_tot;
  y_m = z_m*(k^(alpha))*(l^(1-alpha));
  y_u = z_u*(1-l);
  y_tot = y_m+y_u;
  i = k-(1-delta)*k(-1);
  ln(z_m) = rho_1*ln(z_m(-1))+e_m;
  ln(z_u) = rho_2*ln(z_u(-1))+e_u;
  ln(tau_m) = rho_3*ln(tau_m(-1))+e_tm;
  ln(tau_y) = rho_4*ln(tau_y(-1))+e_ty;
end;

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

initval;
  c     = 0.04050;
  k     = 8.77967;
  i     = 0.87796;
  l     = 0.735;
  r     = 0.04540;
  w     = 0.40718;
  g     = 0.56134;
  y_m   = 1.66632;
  y_u   = 0.265;
  y_tot = 1.93132;
  z_m   = 1;
  z_u   = 1;
  tau_m = 0.275;
  tau_y = 0.335;	
  e_m   = 0;
  e_u   = 0;
  e_tm  = 0;
  e_ty  = 0;
end;

shocks;
  var e_m  = sigmae_m^2;
  var e_u  = sigmae_u^2;
  var e_tm = sigmae_tm^2;
  var e_ty = sigmae_ty^2;
end;

steady;

solve_algo = 0

stoch_simul(hp_filter = 1600, order = 1);

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

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