% RBC Model with Household Production. 
%
% Benhabib, Rogerson, and Wright Model
%
% I borrow their calibration.
%
% Jesus Fernandez-Villaverde
% Philadelphia, March 4th, 2005.

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

close all;

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

var c k n_f n_i l y_f y_i i z_f;
varexo e_f e_i;

parameters beta tau rho psi Time alpha theta_f theta_i gamma delta rho_1 sigmae_f;

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

% Parameters taken directly from BRW
beta  = 0.96;
tau   = 0.25;
rho   = 0.2;
psi   = 0.12;
Time  = 1000;
alpha = 0.36;
theta_f = 1;
theta_i = 23;         
gamma = 0.425;
delta = 0.08;
rho_1 = 0.99;
sigmae_f = 0.001;
%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model; 
 c/c(-1) = beta*((1-tau)*theta_f*(z_f)*alpha*((k)^(alpha-1))*(((n_f))^(1-alpha))+1-delta);
 c*psi/(Time-n_f-n_i)=(1-tau)*z_f*theta_f*(1-alpha)*(k^alpha)*(n_f)^(-alpha);
 c*psi/(Time-n_f-n_i)=(1-rho*tau)*theta_i*gamma*(n_i)^(gamma-1);
 Time-n_f-n_i = l;
 c+k = (1-tau)*exp(z_f)*theta_f*(k^alpha)*(n_f^(1-alpha))+(1-delta)*k(-1)+(1-rho*tau)*theta_i*(n_i^(gamma));
  y_f = exp(z_f)*theta_f*(k^alpha)*(n_f^(1-alpha));
  y_i = theta_i*(n_i^(gamma));
  i = k-(1-delta)*k(-1);
  z_f = rho_1*z_f(-1)+e_f;
end;

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

initval;
  c    = 802.2348;
  k    = 2566.727;
  i    = 205.3381;
  n_f  = 738.6766;
  n_i  = 79.2343;
  l    = 182.0919;
  y_f  = 1156.611;
  y_i  = 147.48;
  z_f  = 0;
  e_f  = 0;
end;

shocks;
  var e_f = sigmae_f^2;
end;

steady;

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);