//DYNARE CODE

% Implementation of the ZLB
% Based on Holden, Thomas and Michael Paetz ( 2012 ), "Efficient Simulation of DSGE Models with Inequality Constraints"
% 
% DEFINITIONS:
% ------------
% TShadow:  Number of periods you believe the ZLB binds
% TIRF:     Number of periods for IRF graphs
% TInitIRF: Number of periods for IRF simulation
% TSimul:   Number of periods for complete simulation
% TDrop:    Number of simulation periods to be dropped before plotting

@#define TShadow = 30  
@#define TIRF = 15    
@#define TSimul = 300   
@#define TDrop = 100    

@#if TShadow > TIRF
	@#define TInitIRF = TShadow
@#else
	@#define TInitIRF = TIRF
@#endif



//VARIABLE DEACLARATION

var x pi i r rbar ybar y n c wp mcr mu shock_i shock_v shock_g shock_a;
varexo epsilon_v epsilon_a epsilon_g epsilon_i;

//Implementation of the ZLB  Based on Holden, Thomas and Michael Paetz ( 2012 ), "Efficient Simulation of DSGE Models with Inequality Constraints"
@#for index in 0:( TShadow - 1 )
	varexo epsilon_shadow_@{index};
	var LAG0_shadow_@{index};
	@#for offset in 1:index
		var LAG@{offset}_shadow_@{index};
	@#endfor
	@#define offset = 10000000000
@#endfor
@#define index = 10000000000

//MODEL PARAMETERS

parameters rho_g rho_v rho_a rho_i sigma phi beta lambda kappa gamma_pi gamma_x epsilon mue;

rho_a=0.9;
rho_g=0.9;
rho_v=0.5;
rho_i=0.7;
sigma=1;
phi=1;
beta=0.99;
theta=0.01;
lambda=(1-theta)*(1-beta*theta)/theta;
kappa=lambda*(sigma+phi);
gamma_pi=1.5;
gamma_x=0.5/4;
epsilon=6;
mue=log(epsilon/(epsilon-1));

//MODEL
model (linear);

//Implementation of the ZLB  Based on Holden, Thomas and Michael Paetz ( 2012 ), "Efficient Simulation of DSGE Models with Inequality Constraints"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'shadow_price' is defined so that the shadow price shocks hit the equation of the constrained variable in consecutive periods: %
% ------------------------------------------------------------------------------------------------------------------------------ %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

@#for index in 0:( TShadow - 1 )
    LAG0_shadow_@{index} = epsilon_shadow_@{index};
	@#for offset in 1:index
		LAG@{offset}_shadow_@{index}( 0 ) = LAG@{offset-1}_shadow_@{index}( -1 );
	@#endfor
	@#define offset = 10000000000
@#endfor
@#define index = 10000000000
	
#shadow_price = 0
@#for index in 0:( TShadow - 1 )
	+ LAG@{index}_shadow_@{index}( 0 )
@#endfor
@#define index = 10000000000
;

x=-(1/sigma)*(i-pi(+1)-rbar)+x(+1);
pi=beta*pi(+1)+kappa*x;
i=gamma_pi*pi+gamma_x*x+shock_i + shadow_price;
r=i-pi(+1);
rbar=-sigma*((1+phi)/(phi+sigma))*(1-rho_a)*shock_a+sigma*(phi/(phi+sigma))*(1-rho_g)*shock_g;
ybar=((1+phi)/(phi+sigma))*shock_a+(sigma/(phi+sigma))*shock_g;
y=x+ybar;
y=n+shock_a;
y=c+shock_g;
wp=phi*n+sigma*c;
mcr=wp-shock_a;
mu=shock_a-wp;

shock_g=rho_g*shock_g(-1)+epsilon_g;
shock_a=rho_a*shock_a(-1)+epsilon_a;
shock_v=rho_v*shock_v(-1)+epsilon_v;
shock_i = rho_i * shock_i( -1 ) -  0.1 * epsilon_i;
end;
initval;
x=0; 
pi=0; 
i=0; 
r=0; 
rbar=0; 
ybar=0; 
y=0; 
n=0; 
c=0; 
wp=0; 
mcr=0; 
mu=0;
end;

check;
steady;


//SHOCKS

shocks;

var epsilon_v;stderr 1;
var epsilon_g;stderr 1;
var epsilon_a;stderr 1;

//Implementation of the ZLB  Based on Holden, Thomas and Michael Paetz ( 2012 ), "Efficient Simulation of DSGE Models with Inequality Constraints"
var epsilon_i; stderr 1;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Standard deviation of shadow price shocks is set to one: %
% -------------------------------------------------------- %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

@#for index in 0:( TShadow - 1 )
	var epsilon_shadow_@{index} = 1;
@#endfor
@#define index = 10000000000

end;
t0 = tic;

stoch_simul( order = 1, irf = @{TInitIRF}, noprint, nograph, nomoments, periods=0 ) y ybar i mcr mu;

ZLBData = INITIAL_CHECKS( @{TShadow}, @{TInitIRF}, 'i', var_list_ );
display( sprintf( '\n%s', 'Computation time for set-up and initial feasibility checks:' ) ); toc( t0 ); 
t1 = tic;
[ IRF_, IRF_CONS_ ] = STEADYSTATE_IRF_CONS( ZLBData );
display( sprintf( '\n%s', 'Additional computation time for IRFs:' ) ); toc( t1 ); 
t2 = tic;
[ SIMU_, SIMU_CONS_ ] = SIMU_CONS(  ZLBData, @{TSimul} );
display( sprintf( '\n%s', 'Additional computation time for simulations:' ) ); toc( t2 );

%%%%%%%%%%%%%%%%%%%%%
% Generate Figures: %
% ----------------- %
%%%%%%%%%%%%%%%%%%%%%

exovars( 1, 1:( size( M_.exo_names, 1 ) - @{TShadow} ) ) = cellstr( M_.exo_names( 1:( size( M_.exo_names, 1 ) - @{TShadow} ), : ) );
vars( 1, 1:size( var_list_ ) ) = cellstr( var_list_( 1:size( var_list_ ), : ) );

xx =( 0:@{TIRF}-1 )';
for exovar = exovars
    @#define number = 0
    for var = vars
        @#define number = number +1
        figure( 'Name', cell2mat( [ 'IRF: ' var ', Shock: ' exovar ] ) );
        hold on;
        plot( xx, IRF_.( cell2mat( [ var '_' exovar ] ) )( 1:@{TIRF} ), '-k', 'MarkerSize', 4, 'LineWidth', 1.5 );
        plot( xx, IRF_CONS_.( cell2mat( [ var '_' exovar ] ) )( 1:@{TIRF} ), '--k', 'MarkerSize', 4, 'LineWidth', 1.5 );
        plot( xx, zeros( @{TIRF} )', '-k', 'MarkerSize', 2, 'LineWidth', 0.5 );
        set( gca, 'PlotBoxAspectRatio', [ 1 2 1 ] );  
        set( gca, 'fontsize', 25 ); 
        title(var);
        axis tight;
        hold off;
    end
end

yy =( 0:@{TSimul}-1-@{TDrop} )';
@#define number = 0
for var = vars
	@#define number = number +1
	figure( 'Name', cell2mat( [ 'Simulation: ' var ] ) );
	hold on;
	plot( yy, SIMU_.( cell2mat( [ var  ] ) )( 1:( @{TSimul}-@{TDrop} ) ), '-k', 'MarkerSize', 4, 'LineWidth', 1.5 );
	plot( yy, SIMU_CONS_.( cell2mat( [ var ] ) )( 1:( @{TSimul}-@{TDrop} ) ), '--k', 'MarkerSize', 4, 'LineWidth', 1.5 );
	plot( yy, zeros( @{TSimul}-@{TDrop} ), '-k', 'MarkerSize', 2, 'LineWidth', .5 );
	set( gca, 'PlotBoxAspectRatio', [ 2.5 1 1 ] );  
	set( gca, 'fontsize', 25 ); 
    title( var);
	axis tight;
	hold off;
end