% Real Business Cycle Model

close all;

estimation_dataset= input('Which dataset? Simulated rbc_simul = 1; Simulated rbc_estimmodel = 2   '); % Select the Dataset

if estimation_dataset==1
          %  use simulated data generated with rbc_simul.mod
load simuldata
xsim  = x;
ysim  = y;
save simuldata1 ysim xsim;
elseif estimation_dataset==2
          %  use simulated data generated with rbc_estimmodel.mod
load simuldata_estimmodel
xsim  = x;
ysim  = y;
save simuldata1 ysim xsim;
end

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

var y c k l x z g eps_z eps_g ysim xsim;

varexo eps_zg u_z u_g;
parameters gn gz betta delta psi sigma theta rho_z sig_z rho_g sig_g q12 zbar gbar;


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

gn    = 0.00243918275778010;
gz    = 0.00499789993972673;
betta = 0.9722^(1/4)/(1+gz); 
delta = 1-(1-0.0464)^(1/4);
psi   = 2.24;
sigma = 1.000001;
theta = 0.35;

zbar  = 0.0023;  
gbar  = -0.0382;

rho_z = 0.8;
sig_z = 0.0086;
rho_g = 0.8;
sig_g = 0.0248;

q12 = 0.0001;



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


model;
% Intratemporal Optimality Condition
(psi*c)/(1-l) = (1-theta)*(k(-1)^theta)*(l^(-theta))*(exp(z)^(1-theta)); 
% Intertemporal Optimality Condition
(c^(-sigma))*((1-l)^(psi*(1-sigma))) = betta*((c(+1))^(-sigma))*((1-l(+1))^(psi*(1-sigma)))*(theta*(k^(theta-1))*((exp(z(+1))*l(+1))^(1-theta))+(1-delta));
% Aggregate Resource Constraint
y = c + exp(g) + x;
% Capital Accumulation Law
(1+gz)*(1+gn)*k = (1-delta)*k(-1) + x;
% Production Function
y = (k(-1)^theta)*((exp(z)*l)^(1-theta));

% Stochastic Processes
z = zbar + rho_z*z(-1) + eps_z;
g = gbar + rho_g*g(-1) + eps_g;
eps_z =     eps_zg + u_z; 
eps_g = q12*eps_zg + u_g;

ysim = y;
xsim = x;

end;


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


shocks;
var u_z ; stderr sig_z;
var u_g ; stderr sig_g;
var eps_zg = 1; 
end;


initval;
z = 0;
k = (betta*theta)^(1/(1-theta))*exp(zbar/(1-rho_z));
c = (1-betta*theta)/(betta*theta)*k;
y = exp(zbar/(1-rho_z))^(1-theta)*k^theta;
l = k/(((1-betta*(1-delta))/(betta*theta*(exp(zbar/(1-rho_z))^(1-theta))))^(1/(theta-1)));
end;



estimated_params;
stderr u_z,   sig_z, 0.00001,1;
stderr u_g,   sig_g, 0.00001,1;
q12,          q12,        -1,1;
end;

varobs ysim xsim;


estimation(datafile=simuldata1,mode_compute=5, mode_check) ysim xsim;