% Dynare Code for "Risk Matters: The Real Effects of 
%Volatility Shocks" Fernandez-Villalverde, Guerron-Quintana, Rubio-Ramirez, Uribe 2009.
%By Mario Gonzalez

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

close all;

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

var c d k h i r  x y lambda varphi eps_tb eps_r sigma_tb sigma_r;
predetermined_variables k d;
varexo e_ux  e_utb e_ur e_usigmatb e_usigmar;

parameters p_beta p_omega p_nu p_delta p_phid p_D p_phi p_alpha p_sigmax p_r p_rhox p_rhotb p_rhor rho_sigmatb p_sigmatb eeta_tb rho_sigmar p_sigmar eeta_r;

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

p_beta	    = 0.98;
p_omega	    = 1.6;
p_nu	    = 2;
p_delta	    = 0.014; 
p_phid	    = 1.4/10000;
p_D	        = 27;
p_phi	    = 280; 
p_alpha	    = 0.32;
p_sigmax	= 0.0075;
p_r	        = 0.02;
p_rhox	    = 0.95;
p_rhotb	    = 0.95;
p_rhor	    = 0.97;
rho_sigmatb	= 0.94;
p_sigmatb	=-8.05;
eeta_tb	    = 0.13;
rho_sigmar	= 0.94;
p_sigmar	=-5.71;
eeta_r	    = 0.46;


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

model; 
  y           =  k^p_alpha*(exp(x)*h)^(1-p_alpha);
  c           = (c-h^p_omega/p_omega)^-p_nu;
  lambda/(1+r)= lambda*p_phid*(d(+1)-p_D)+p_beta*lambda(+1);
  varphi      = p_beta*((1-p_delta)*varphi(+1)+p_alpha*lambda(+1)*y(+1)/k(+1));
  h^p_omega   = (1-p_alpha)*y;
  lambda      = varphi*(1-0.5*p_phi*((i-i(-1))/i(-1))^2-p_phi*i/i(-1)*((i-i(-1))/i(-1)))+p_beta*(varphi(+1)*p_phi*(i(+1)/i)^2*((i-i(-1))/i(-1)));
  k(+1)       = (1-p_delta)*k+(1-0.5*p_phi*(i/i(-1)-1)^2)*i;
  x           = p_rhox*x(-1)+exp(p_sigmax)*e_ux;
  y           = c+i+d-d(+1)/(1+r)+0.5*p_phid*(d(+1)-d)^2;
  r           = p_r + eps_tb + eps_r;
  eps_tb      = p_rhotb*eps_tb(-1)+exp(sigma_tb)*e_utb;
  eps_r       = p_rhor*eps_r(-1)+exp(sigma_r)*e_ur;
  sigma_tb    = (1-rho_sigmatb)*p_sigmatb+rho_sigmatb*sigma_tb(-1)+eeta_tb*e_usigmatb;
  sigma_r     = (1-rho_sigmar)*p_sigmar+rho_sigmar*sigma_r(-1)+eeta_r*e_usigmar;
  %w           = (1-p_alpha)*exp(x)*y/h;
  %R           = p_alpha*y/k;
  %d(+1)/(1+r) = d-w*h-R*k+c+i+0.5*p_phid*(d(+1)-d)^2;


 % c_obs      = c-c(-1) +(g_gamma+log(c_g_gamma));
 % i_obs      = i-i(-1)+(g_gamma+log(c_g_gamma));
 % n_obs      = n+(log(0.25));
 % ps_obs     = ps-ps(-1)+(g_gamma+log(c_g_gamma));
 % p_obs      = p-p(-1)+(-g_z-log(c_lambda_z_bar));
end;


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

initval;
	c           =	25.05289233	;
	d           =	27	;
	k           =	273.5329007	;
	h           =	60.03954594	;
	i           =	3.82946061	;
	r           =	0.02	;
	x           =	0	;
	y           =	29.41176471	;
	lambda      =	5.86693E-06	;
	varphi      =	5.86693E-06	;
	eps_tb      =	0	;
	eps_r       =	0	;
	sigma_tb	=	-8.05	;
	sigma_r     =	-5.71	;
end;


shocks;
var 	e_ux        =	p_sigmax^2;
var 	e_utb       =	1;
var 	e_ur        =	1;
var 	e_usigmatb	=	1;
var 	e_usigmar	=	1;	
end;

steady_state_model;
    y       = (p_r/(1*p_r))*(p_D/0.018);
    k       = -p_alpha*y/((1-p_delta)*p_beta^-1);
    i       = p_delta*k;
    c       = 0.982*y-i;
    h       = (y/k)^(1/(1-p_alpha))*k;
    lambda  = (c-h^p_omega/p_omega)^-p_nu;
    varphi  = lambda;
    eps_tb	= 0;
	eps_r	= 0;
	sigma_tb= -8.05	;
	sigma_r	= -5.71	;
end;




%----------------------------------------------------------------
% 5. Bayesian Estimation
%----------------------------------------------------------------


%estimated_params;
%c_h, beta_pdf, 0.33, 0.235;
%c_Omega, gamma_pdf, 2, 2; 
%c_delta_delta, gamma_pdf, 1, 1;
%c_xi_bar, beta_pdf, 0.2, 0.2;
%c_mu, gamma_pdf, 2, 2;
%rho_a, beta_pdf, 0.5, 0.2;
%rho_am, beta_pdf, 0.5, 0.2;      
%rho_z, beta_pdf, 0.5, 0.2;       
%rho_zm, beta_pdf, 0.5, 0.2;      
%rho_theta, beta_pdf, 0.5, 0.2;	
%rho_psi, beta_pdf, 0.5, 0.2;     
%rho_xi, beta_pdf, 0.5, 0.2;      
%stderr e_a, inv_gamma_pdf, 0.01, inf; 
%stderr e_am, inv_gamma_pdf, 0.01, inf;         
%stderr e_z, inv_gamma_pdf, 0.01, inf;          
%stderr e_zm, inv_gamma_pdf, 0.01, inf;         
%stderr e_psi, inv_gamma_pdf, 0.01, inf;        
%stderr e_xi, inv_gamma_pdf, 0.01, inf;         
%stderr e_theta, inv_gamma_pdf, 0.01, inf;      
%end;



stoch_simul(solve_algo=1,periods=2096, irf=40, order = 3, pruning  );
































%----------------------------------------------------------------
% 6. 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);