% 
% 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

%%%%%%%%%%%%%%%%%%%%%%%%% BEGINNING OF THE CONVENTIONAL .mod FILE %%%%%%%%%%%%%%%%%%%%%%%%%

var y c1 c2 i pai w n k rl rd re l d e r eta lcb a ut vp ui x rdcb rlcb rdaux rlaux;

varexo eps_a, eps_ut, eps_vp,eps_ui ;

%%%%%%%%%%%%%%%%%%%%%>>> BEWARE! <<<%%%%%%%%%%%%%%%%%%%%%%%
% Shadow price shocks need to be the last shocks defined! %
% ------------------------------------------------------- %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

@#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

parameters gammac1 gammac2 gammai gammal gammad gammadss gammae beta1 beta2 
sigma1 sigma2 phi theta alpha miu delta psi phipai phiy phir phil 
phid philcb philpai phily kappa phietaeta phietapai tao bigphi 
lambda cl ce cd rr_ss re_ss rd_ss rl_ss d_ss e_ss r_ss l_ss  rhoa rhout rhovp rhoui phirl phird;

gammac1=0.2;
gammac2=0.3;
gammai=0.5;
gammal=0.4;
gammad=0.3;
gammadss=0.5;
gammae=0.01;
beta1 = 0.995;
beta2 = 0.985;
sigma1 = 1;
sigma2 = 1;
phi = 1/3;
theta = 0.75;
alpha = 0.5;
miu = 0.1;
delta = 0.03;
psi = 2;
phipai = 1.01;
phiy = 0;
phir = 0.5;
phil=0.7;
phid=0.7;
philcb=0.2;
philpai=50;
phily=5;
kappa=50;
phietaeta=0.6;
phietapai=10;
tao = 0.2;
bigphi =(theta+tao*(1-theta*(1-theta*(beta1)))^(-1));
lambda=(1-tao)*(1-theta)*(1-beta1*theta)*bigphi;
cl = 1;
ce = 1;
cd = 2;
re_ss =1.009;
rr_ss =1.009;
rd_ss=beta1^(-1);
rl_ss=beta2^(-1);
l_ss=(rl_ss-r_ss)/cl;
d_ss=(eta*rr_ss+(1-eta)*r_ss-rd_ss)/((1-eta)*(1-eta)*cd);
e_ss=(re_ss-r_ss)/ce;
r_ss=(cd*ce*rl_ss+cd*cl*re_ss)/(cl*ce+cd*ce+cd*cl)+(cl*ce*(rd_ss-eta*rr_ss))/(1-eta)*(cl*ce+cd*ce+cd*cl);
rhoa=0.7;
rhout=0.7;
rhovp=0.7;
rhoui=0.7;
phirl=0.2;
phird=0.2;

model( linear );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% '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
;
//%---------------------------------------------
//% model equations
//%---------------------------------------------

//% 1
y=gammac2*c2+gammac1*c1+gammai*i;

//% 2
c1=c1(+1)-sigma1^(-1)*(rd-pai(+1)-(1-rhovp)*vp);

//% 3
w=sigma1*c2+phi*n;

//% 4
c2=c2(+1)-sigma2*(rl*pai(+1));
//% 5
gammac2*c2=y/(1+miu)+gammal*l-(1-alpha)*(w+n)/(1+miu)-rl_ss*gammal*(rl(-1)-pai-l(-1))-gammai*i;

//% 6
k=delta*i+(1-delta)*k(-1);

//% 7
i=ui+k(-1)+beta2*(i(+1)-ui(+1)-k)+(1-beta2*(1-delta))*(y(+1)-k)/psi+sigma2*(c2-c2(+1));

//% 8
y=a+alpha*k+(1-alpha)*n;

//% 9
w=y-n;

//% 10
pai=bigphi*(theta*beta1*pai(+1))-lambda*x+ut;

//% 11 

rl_ss*rl=r_ss*r-kappa*l_ss*lcb+(cl+kappa)*l_ss*l+ shadow_price;

//% 12 
rd_ss*rd=(rr_ss-r_ss+2*cd*d_ss*(1-eta))*eta+eta*rr_ss*0.01+(1-eta)*r_ss*r-cd*(1-eta)^2*d_ss*d;

//% 13 
re_ss*re=r_ss*r+ce*e_ss*e;

//% 14 
gammal*l=(1-eta)*gammad*d-gammadss*eta-gammae*e;

//% 15 
r=phipai*pai+phiy*y+phir*r(-1);

//% 16
lcb=-(1-philcb)*(philpai*pai+phily*y)+philcb*lcb(-1);

//% 17
eta=(1-phietaeta)*phietapai*pai+phietaeta*eta(-1);

//% 18
a=rhoa*a(-1)*eps_a;

//% 19
ut=rhoa*ut(-1)*eps_ut;

//% 20
vp=rhovp*vp(-1)*eps_vp;

//% 21
ui=rhoui*ui(-1)*eps_ui;

//% 22
x=y-n-pai;

//% 23
rdcb=(1-phird)*r+phird*rdcb(-1);

//% 24
rlcb=(1-phirl)*r+phirl*rlcb(-1);

//% 25
rlaux=rl-rlcb+shadow_price;

//% 26
rdaux=rdcb-rd+shadow_price;
end;


initval;
y = 0;
pai = 0;
end;

check;
steady;

shocks;

var eps_a; stderr 1;
var eps_ut; stderr 1;
var eps_vp; stderr 1;
var eps_ui; 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 pai i;

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