% Vermandel 2013


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

close all;

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

var 
y pi r 
yema yemb ysema ysemb ycove ycore 
phi_y
s_b s_p s_r;

varexo e_b e_p e_r;

parameters  sigmaC sigmaL chi rho phi_pi theta
			beta alph_a  alph_b  alph_c chi_min cay  delt chi_max kay
			% shocks
			rho_b rho_p rho_r u;

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

beta    	= 0.991;			% discount factor
sigmaC		= 2;				% risk aversion consumption
sigmaL		= 1;				% labor disutility
theta		= 3/4;				% new keynesian Philips Curve, forward term
epsilon		= 6;				% subsituability/mark-up on prices
rho     	= 0.01;				% MPR Smoothing
phi_pi		= 1.5;				% MPR Inflation
alph_a      = 7/100;			% estimator coefficient
alph_b		= 10/100;			% estimator coefficient
alph_c		= 50/100;			% estimator coefficient
cay         = 1;				% central MPR phi_dy value
delt		= 0.1;				% difference agressiveness high/low vol
chi_min		= 0.3;				% steady state volatility

chi_max     = 2;				% peak volatility
kay         = 100;				% linear/asymptotic response to vol

% shock processes
rho_b   = 0.9;
rho_p   = 0.9;
rho_r 	= 0.4;
u       = 0.1;

% steady states
R		= 1/beta;
H		= 1/3;
MC		= (epsilon-1)/epsilon;
W		= MC;
Y		= H;
chi		= W*Y^-sigmaC*H^-sigmaL;

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

model(linear); 
	% IS curve

	y = y(+1) - 1/sigmaC*(r-pi(+1)) + s_b;

	% AS curve

	pi = beta*pi(+1) + ((1-theta)*(1-beta*theta)/theta)*(sigmaC+sigmaL)*y + s_p;
	

    % Rolling volatility estimator

	yema = (1-alph_a)*yema(-1) + alph_a*y;
	yemb = (1-alph_b)*yemb(-1) + alph_b*y;
	ysema = (1-alph_a)*ysema(-1) + alph_a*((y-yema)*(y-yema(-1)));
	ysemb = (1-alph_b)*ysemb(-1) + alph_b*((y-yemb)*(y-yemb(-1)));
	ycove = (1-alph_c)*(ycove(-1) + alph_c*((y-yema(-1))*(y-yemb(-1))));
	ycore = ycove/(sqrt(ysemb)*sqrt(ysema));


	% MPR coefficient logistic function

	phi_y = ((-cay - delt)*(exp(kay*ycore + (kay*(chi_max))/2)) + 
	(cay - delt)*(exp(kay*ycore + (kay*(chi_min))/2)) - 
	(cay - delt)*(exp(kay*(chi_max) + (kay*(chi_min))/2)) + 
	(cay + delt)*(exp((kay*(chi_max))/2 + kay*(chi_min))))/
	((-(exp((kay*(chi_max))/2)) + exp((kay*(chi_min))/2))*
	(exp(kay*ycore) + exp((1/2)*kay*((chi_max) + (chi_min)))));


    	% Monetary Policy Rule

	r = rho*r(-1) +  (1-rho)*(phi_pi*pi + phi_y*y) + s_r  ;


	% shocks
	s_b = rho_b*s_b(-1)+e_b;
	s_p = rho_p*s_p(-1)+e_p-u*e_p(-1);
	s_r = rho_r*s_r(-1)+e_r;
end;


initval;
yema = 0.6;
yemb = 0.5;
ysema = 1;
ysemb = 1;
ycove = 0.5;
ycore = 0.5;
phi_y = 0.2;
end;

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



shocks;
var e_b;  stderr 1;
var e_p;  stderr 1;
var e_r;  stderr 1;
end;


check;

% GENERATING IRF
% benchmark
%set_param_value('phi_pi',1.5);
% inflation target
%set_param_value('phi_pi',1.7);
%set_param_value('phi_pi',1.5);
%set_param_value('cay',0.5);

stoch_simul(order=3,irf=30) y r pi phi_y ycore;

%IRF_to_txt( char('y','r','pi'),char('e_b','e_p','e_r'))

%set_param_value('theta',0.01);
set_param_value('delt',0.3);

stoch_simul(order=3,irf=30,periods=200) y r pi phi_y ycore;
dynasave(ma);