% Dynare Code for "Bayesian DSGE Model of Stock Market Bubbles 
%      and Business Cycles" Miao, Wang, XU, 2012.
%By Mario Gonzalez

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

close all;

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

var y c i b q x p n big_lambda k u m lambda_a a lambda_z zm psi xi theta epsilon g_gamma g_z ;
predetermined_variables k;
varexo e_a e_am e_z e_zm e_psi e_xi e_theta ;

parameters c_betti c_alpha c_delta1 c_delta_e c_N c_g_gamma c_lambda_z_bar c_u_ss c_i_y c_k0k c_theta_bar c_omega c_h c_Omega c_delta_delta c_xi_bar 
c_mu rho_a rho_am rho_z rho_zm rho_theta rho_psi rho_xi sigma_a sigma_am sigma_z sigma_zm sigma_theta sigma_psi sigma_xi c_y b_y phi_x phi_k phi_g c_delta_d1 c_1_Phie G;

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

c_betti     = 0.99;
c_alpha     = 0.3;
c_delta1	= 0.025;
c_delta_e   = 0.020;
c_N         = 0.25;
c_g_gamma	= 1.0042;
c_lambda_z_bar	= 1.0121;
c_u_ss      = 1;
c_i_y       = 0.2;
c_k0k       = 0.2;
c_theta_bar	= 1;
c_omega     = 0.5;
c_h         = 0.33;
c_Omega     = 2;
c_delta_delta	= 1;
c_xi_bar	= 0.2;
c_mu        = 2;
rho_a       = 0.5;
rho_am      = 0.5;
rho_z       = 0.5;
rho_zm      = 0.5;
rho_theta	= 0.5;
rho_psi     = 0.5;
rho_xi      = 0.5;
sigma_a     = 0.01;
sigma_am	= 0.01;
sigma_z     = 0.01;
sigma_zm	= 0.01;
sigma_theta	= 0.1;
sigma_psi	= 0.01;
sigma_xi	= 0.01;
  phi_x      = c_alpha/(c_lambda_z_bar*c_g_gamma*c_theta_bar-(1-c_delta1)*c_betti*(1-c_delta_e)*c_theta_bar-c_xi_bar*(1-c_betti*(1-c_delta_e)*c_theta_bar));
% c_mu       =  
  c_i_y      = 0.2; % (c_1_Phie*(phi_k-(1-c_delta1))*phi_x)/(G+c_1_Phie);
  phi_k      = ((1-c_delta_e)/(c_lambda_z_bar*c_g_gamma)+c_delta_e*c_k0k)^(-1);
  c_delta_d1 = (c_alpha/(c_betti*(1-c_delta_e)*c_theta_bar*phi_x));
  c_1_Phie   = (1/(c_betti*(1-c_delta_e)*c_theta_bar)-1)/((c_i_y^(-1))*(phi_k-(1-c_delta1))*phi_x-1);
  G          = 1/(c_betti*(1-c_delta_e)*c_theta_bar)-1; 
  phi_g      = -c_1_Phie/G-1;
  c_y        = 1-c_i_y;
  b_y        = c_i_y*(1/c_1_Phie)-c_alpha-c_xi_bar*phi_x;

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

model(linear); 
  y          = c_y*c+c_i_y*i ;
  i          = (c_alpha/(c_alpha+c_xi_bar*phi_x+b_y))*y+(c_xi_bar*phi_x/(c_alpha+c_xi_bar*phi_x+b_y))*(xi+q+x)+(b_y/(c_alpha+c_xi_bar*phi_x+b_y))*b-c_mu*epsilon-p;
  y          = c_alpha*(u+x)+(1-c_alpha)*n;
  psi        = big_lambda+y-n;
  k(+1)      = -(c_delta_d1/phi_k)*u+((1-c_delta1)/phi_k)*x+(1-((1-c_delta1)/phi_k))*(i-(c_mu/phi_g)*epsilon);
  y-x+(1-c_betti*(1-c_delta_e)*c_theta_bar)*phi_g*epsilon = q + (1+c_delta_delta)*u;
  q          = (big_lambda(+1)-big_lambda)+(q(+1)-g_z(+1)-g_gamma(+1))+((c_betti*(1-c_delta_e)*c_delta_d1)/(c_lambda_z_bar*c_g_gamma))*c_delta_delta*u(+1)+((c_xi_bar*c_betti*(1-c_delta_e)*G)/(c_lambda_z_bar*c_g_gamma))*(xi(+1)+phi_g*epsilon(+1));
  x          = ((1-c_delta_e)/(c_lambda_z_bar*c_g_gamma))*phi_k*(k-g_z-g_gamma);
  p          = (1+c_betti)*c_Omega*(c_g_gamma^2)*(c_lambda_z_bar^2)*i+c_Omega*(c_g_gamma^2)*(c_lambda_z_bar^2)*(g_gamma+g_z)-c_Omega*(c_g_gamma^2)*(c_lambda_z_bar^2)*i(-1)-c_betti*c_Omega*(c_g_gamma^2)*(c_lambda_z_bar^2)*(i(+1)+g_z(+1)+g_gamma(+1));
  b          = (big_lambda(+1)-big_lambda+b(+1))+(1+c_betti*(1-c_delta_e)*c_theta_bar)*phi_g*epsilon(+1)+((1-(1-c_delta_e)*c_theta_bar)/((1-c_delta_e)*c_theta_bar))*m(+1);
  m          = (1-c_delta_e)*c_theta_bar*m(-1)+(1-c_delta_e)*c_theta_bar*theta(-1);
  big_lambda = (c_g_gamma/(c_g_gamma-c_betti*c_h))*(-(c_g_gamma/(c_g_gamma-c_h))*c+(c_h/(c_g_gamma-c_h))*(c(-1)-g_gamma))-((c_betti*c_h)/(c_g_gamma-c_betti*c_h))*(-(c_g_gamma/(c_g_gamma-c_h))*(c(+1)+g_gamma(+1))+(c_h/(c_g_gamma-c_h))*c);
  g_gamma    = (c_alpha/(1-c_alpha))*(lambda_z+zm-zm(-1))+(lambda_a+a-a(-1));
  g_z        = lambda_z+zm-zm(-1);   
  lambda_a   = rho_a*lambda_a(-1)+e_a;
  a          = rho_am(-1)*a+e_am;
  lambda_z   = rho_z*lambda_z(-1)+e_z;
  zm         = rho_zm*zm(-1)+e_zm;
  psi        = rho_psi*psi(-1)+e_psi;
  xi         = rho_xi*xi(-1)+e_xi;
  theta      = rho_theta*theta(-1)+e_theta;
  epsilon    = p-q;
end;

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

initval;
y           =	0;
c           =	0;
i           =	0;
b           =	0;
q           =	0;
x           =	0;
p           =	0;
n           =	0;
big_lambda	=	0;
k           =	0;
u           =	0;
m           =	0;
lambda_a	=	0;
a           =	0;
lambda_z	=	0;
zm          =	0;
psi         =	0;
xi          =	0;
theta       =	0;
e_a         =	0;
e_am        =	0;
e_z         =	0;
e_zm        =	0;
e_psi       =	0;
e_xi        =	0;
e_theta     =	0;
end;

shocks;
var e_a  =sigma_a^2;     
var e_am  =sigma_am^2;	
var e_z  =sigma_z^2;     
var e_zm  =sigma_zm^2;	
var e_theta  =sigma_theta^2;	
var e_psi  =sigma_psi^2;	
var e_xi  =sigma_xi^2;	
end;

steady;
check;

stoch_simul(periods=1000);

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

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