// With Nominal GDP indexed bond
// With bestfit of RS2012

var	A			% Productivity 
	c			% Consumption
	Disp		% Price dispersion
	DZ  		% Trend growth (Zt+1/Zt)
	G			% Government spending
	i			% Short term nominal rate
	ir			% Ex post real short term rate
	L			% Labour
	Mc			% Maginal cost
	p0			% part of price rule
	pi			% inflation
	piavg 		% Average inflation
	pistar		% Target inflation 
	pricebond 	% Price of nominal bond (with maturity n)
	pricebondrn	% Price of nominal bond (risk neutral)  
	termprem	% Term premium 
	V			% Value function
	Valphaexp 	% part of certainty equivalence of value function
	Vkp			% part of certainty equivalence of value function
	wreal		% Real wage 
	Y 			% Income/Production
	ytm 		% YTM of nominal bond
	ytmrn		% YTM of risk neutral bond
	zd 			% part of price function
	zn 			% part of price function
	ehpr		% excess holding period return
	slope		% slope
	sdf			% SDF
%--------------------------
	Q			% Price of typical bond
	QR			%
	QG			% Price of NGDP-indexed bond
	b			% Amount of typical debt
	bG			% Amount of NGDP debt
	T			% real tax
	RG			% real expected return on NGDP bond
	RR			% real return on pi-indexed bond
	GP			% growth risk premium
	IP			% inflation risk premium
	TD			% Total gov. liability
	PSoY		% PS over output
	DoY			% liability over output
;

varexo epsa epsg epsi epsz epspistar;

parameters 
	betatilde	% discount factor (such that beta~ = beta*DZBar^(-phi) = 0.99)
	delta 			% depreciation
	LMax	% From utility kernel (on labour)
	theta			% From CES for final goods 
	xi				% prob. of price no change
	eta				% Labour share of production function 
	ABar piBar 
	DZBar			% gamma in the paper 
	taylrho taylpi tayly	% From monetary function
	rhoa rhog rhoinflavg rhoz rhopistar 
	consoldelta				% From RS(2008b), consol decaying rate
	gssload					% Weight on inflation gap for pistar
	sigmaa sigmag sigmai sigmaz  sigmapistar %sigmap Sigma 
	IES Frisch CRRA			% Structural parameters
%-----------------------------------------------------------------------------
	zeta	omegaG	DoYss
	Gamma1	Gamma2 ;

%% Given Parameter values (Table 1)

	DZBar 		= 1.0025;
	eta 		= 2/3;	
	theta 		= 0.2;	
	xi			= 0.75;	% Baseline = 0.75;	Best = 0.76
	IES			= 0.5;	% Baseline = 0.5; 	Best = 0.09 *****
	betatilde	= 0.99;	
	delta 		= 0.02;	
	LMax		= 3;	
	Frisch		= 2/3;	% Baseline = 2/3;	Best = 0.28
	CRRA		= 75;	% Baseline = 75;	Best = 110
	taylrho		= 0.73; % baseline = 0.73; % rhoi, g_pi, g_y in Table 1
	taylpi 		= 0.53;% baseline = 0.53;
	tayly 		= 0.93;% baseline = 0.93;
	piBar		= 0;	% pi* 
	rhoa 		= 0.95;
	rhog 		= 0.95;
	rhoinflavg	= 0.7;
	rhoz		= 0.90;
	rhopistar	= 0.99;	%0.995;	% Baseline = 0.99	Best = 0.995
	
	sigmaa		= 0.005;	%0.005;
	sigmag		= 0.004;	%0.004;
	sigmai		= 0.003;	%0.003;
	sigmaz		= 0;		%0;		% without prod. trend shock 
	sigmapistar = 0.0005;	% Baseline = 5bp;	Best = 0.0007

	gssload		= 0.003;	% Baseline = 0.01;	Best = 0.003

	zeta		= 1.00;
	omegaG		= 0;
	DoYss		= 1;
	Gamma1		= 0.2;
	Gamma2		= 5;
%% Using Steady State equations (See my derivations)
	ABar		= 1;
	consoldelta	= 0.9848;	% From Rudebusch and Swanson (2008b)
%---------------------------------------------------------------------------


model;

# phi		= (1/IES);
# beta 		= (betatilde/DZBar^(-(1/IES)));
# alpha		= (CRRA - 1/(IES + 1/((LMax-(LMax/3))/Frisch)*(LMax-(LMax/3)))) * (1/(1-1/IES) + 1/(1-((LMax-(LMax/3))/Frisch))*(LMax-(LMax/3)));
# chi0		= (((1/(1+theta))*eta*10^((1-eta)/eta))*((10^((1-eta)/eta)*(LMax/3))-(0.17 * (10^((1-eta)/eta)*(LMax/3)))-((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3)))))^(-1/IES)*(LMax-(LMax/3))^((LMax-(LMax/3))/Frisch));
# YBar		= (10^((1-eta)/eta)*(LMax/3));
# KBar		= (10 * (10^((1-eta)/eta)*(LMax/3)));
# IBar		= ((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3))));
# lss		= (LMax/3);
# GBar		= (0.17 * (10^((1-eta)/eta)*(LMax/3)));
# chi		= ((LMax-(LMax/3))/Frisch);
# css		= ((10^((1-eta)/eta)*(LMax/3))-(0.17 * (10^((1-eta)/eta)*(LMax/3)))-((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3)))));
# mcss		= (1/(1+theta));
# wrealss	= ((1/(1+theta))*eta*10^((1-eta)/eta));
# iss		= (log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES)))));

# vss		= ((((10^((1-eta)/eta)*(LMax/3))-(0.17 * (10^((1-eta)/eta)*(LMax/3)))-((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3)))))^(1-(1/IES))/(1-(1/IES))+((((1/(1+theta))*eta*10^((1-eta)/eta))*((10^((1-eta)/eta)*(LMax/3))-(0.17 * (10^((1-eta)/eta)*(LMax/3)))-((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3)))))^(-1/IES)*(LMax-(LMax/3))^((LMax-(LMax/3))/Frisch)))*(LMax-(LMax/3))^(1-((LMax-(LMax/3))/Frisch))/(1-((LMax-(LMax/3))/Frisch))) / (1-(betatilde/DZBar^(-(1/IES)))*DZBar^(1-(1/IES))));

# vkpss		= (((((10^((1-eta)/eta)*(LMax/3))-(0.17 * (10^((1-eta)/eta)*(LMax/3)))-((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3)))))^(1-(1/IES))/(1-(1/IES))+((((1/(1+theta))*eta*10^((1-eta)/eta))*((10^((1-eta)/eta)*(LMax/3))-(0.17 * (10^((1-eta)/eta)*(LMax/3)))-((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3)))))^(-1/IES)*(LMax-(LMax/3))^((LMax-(LMax/3))/Frisch)))*(LMax-(LMax/3))^(1-((LMax-(LMax/3))/Frisch))/(1-((LMax-(LMax/3))/Frisch))) / (1-(betatilde/DZBar^(-(1/IES)))*DZBar^(1-(1/IES))))*DZBar^(1-(1/IES)));

#znss  		= ((10^((1-eta)/eta)*(LMax/3)) / (1-(betatilde/DZBar^(-(1/IES)))*xi*DZBar^-(1/IES)));
#zdss  		= ((10^((1-eta)/eta)*(LMax/3)) / (1-(betatilde/DZBar^(-(1/IES)))*xi*DZBar^-(1/IES)));
#pricebondss	= (1/(1-(betatilde/DZBar^(-(1/IES)))*consoldelta*DZBar^-(1/IES)));
#pricebondrnss	= (1 /(1-consoldelta*exp(-(log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES))))))));

#ytmss		= (log(consoldelta*(1/(1-(betatilde/DZBar^(-(1/IES)))*consoldelta*DZBar^-(1/IES)))/((1/(1-(betatilde/DZBar^(-(1/IES)))*consoldelta*DZBar^-(1/IES)))-1)) *400);  % yield in annualized percent *)

#ytmrnss		= (log(consoldelta*(1 /(1-consoldelta*exp(-(log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES))))))))/((1 /(1-consoldelta*exp(-(log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES))))))))-1)) *400); 

#Qss		= ((betatilde/DZBar^(-(1/IES)))*DZBar^(-(1/IES)));
#QGss		= ((betatilde/DZBar^(-(1/IES)))/zeta*DZBar^(1-(1/IES)));
#RGBar		= ((betatilde/DZBar^(-(1/IES)))^-1 * DZBar^(1/IES) * zeta);
#TDss		= (DoYss*(10^((1-eta)/eta)*(LMax/3)));
#bBar		= ((DoYss*(10^((1-eta)/eta)*(LMax/3)))*(DZBar^-1 + omegaG*zeta/DZBar/(1-omegaG))^-1);
#bGBar		= (omegaG*((betatilde/DZBar^(-(1/IES)))*DZBar^(-(1/IES)))/(1-omegaG)/((betatilde/DZBar^(-(1/IES)))/zeta*DZBar^(1-(1/IES))));

#TBar 		= ((DoYss*(10^((1-eta)/eta)*(LMax/3))) + (0.17 * (10^((1-eta)/eta)*(LMax/3))) - ((betatilde/DZBar^(-(1/IES)))*DZBar^(-(1/IES)))*((DoYss*(10^((1-eta)/eta)*(LMax/3)))*(DZBar^-1 + omegaG*zeta/DZBar/(1-omegaG))^-1)*(1-omegaG)^-1);

%%	1) Value functions and consumer;s Euler Equation
	V			=	c^(1-phi)/(1-phi) + chi0*(LMax-L)^(1-chi)/(1-chi) + beta*Vkp ; %EQ(1.1)

%%	2) 3) certainty equivalent
	Valphaexp	=	(V(+1)*DZ(+1)^(1-phi)/vss/DZBar^(1-phi))^(1-alpha); %EQ(1.3)
	Vkp			=	(vss*DZBar^(1-phi)) * Valphaexp^(1/(1-alpha));	%EQ(1.2)

%% 4) 5) 6) 7) Price-setting equations
	zn			=	(1+theta)*Mc*Y + beta*xi*(c(+1)/c)^-phi*DZ(+1)^-phi * 
					(V(+1)*DZ(+1)^(1-phi)/Vkp)^-alpha * pi(+1)^((1+theta)/theta/eta) * zn(+1); %EQ(1.5)
	zd			=	Y + beta*xi*(c(+1)/c)^-phi*DZ(+1)^-phi * 
					(V(+1)*DZ(+1)^(1-phi)/Vkp)^-alpha * pi(+1)^(1/theta) * zd(+1); %EQ(1.6)
	p0^(1+(1+theta)/theta*(1-eta)/eta) = zn/zd ; %EQ(1.7)
	pi^(-1/theta) =	(1-xi)*(p0*pi)^(-1/theta) + xi; %EQ(1.8)

%% 8) 9) Marginal cost and real wage
	Mc			=	wreal*(1/eta)*Y^((1-eta)/eta) / A^(1/eta)/KBar^((1-eta)/eta); %EQ(1.9)
	chi0*(LMax-L)^-chi/c^-phi = wreal ; %EQ(1.10)

%% 10) 11) 12) Output equations
	Y			=	A*KBar^(1-eta)*L^eta / Disp; %EQ(1.11)
	Disp^(1/eta)=	(1-xi)*p0^(-(1+theta)/theta/eta) + xi*pi^((1+theta)/theta/eta)*Disp(-1)^(1/eta); %EQ(1.12)
	c			=	Y - G - IBar; %EQ(1.13)

%% 13) 14) Monetary Policy Rule
	log(piavg)	=	rhoinflavg*log(piavg(-1)) + (1-rhoinflavg)*log(pi); %EQ(1.14)
	4*i			=	taylrho * 4*i(-1) + 
					(1-taylrho) * ( 4*log(1/beta*DZBar^phi) + 4*log(piavg) + taylpi*(4*log(piavg)-pistar) + 
					tayly*(Y-YBar)/YBar ) + epsi;  %EQ(1.15)
					
%% 15) 16) 17) 18) Exogenous Shocks (EQ (1.16)~(1.19))
	log(A/ABar)	=	rhoa*log(A(-1)/ABar) + epsa ;
	log(DZ/DZBar)=	rhoz*log(DZ(-1)/DZBar)+epsz ;	
	log(G/GBar)	=	rhog*log(G(-1)/GBar) + epsg ;
	pistar		= 	(1-rhopistar)*piBar + rhopistar*pistar(-1) + gssload*(4*log(piavg)-pistar) + epspistar;

%% 19) - 26) Long-Term Bond Price; Yield; and Term Premium *)
	ir			=	log(exp(i(-1))/pi); %EQ(1.20)
	pricebond	=	1 + beta*consoldelta*(c(+1)/c)^-phi*DZ(+1)^-phi * 
					(V(+1)*DZ(+1)^(1-phi)/Vkp)^-alpha / pi(+1) * pricebond(+1); %EQ(1.21)
	pricebondrn	=	1 + consoldelta*exp(-i)*pricebondrn(+1); %EQ(1.22)
	ytm			=	log(consoldelta*pricebond/(pricebond-1)) * 400; %EQ(1.23)
	ytmrn  		=	log(consoldelta*pricebondrn/(pricebondrn-1)) * 400; %EQ(1.24)
	termprem 	=	100 * (ytm  - ytmrn ); %EQ(1.25)
	ehpr  		= 	((consoldelta*pricebond + exp(i(-1)))/pricebond(-1) - exp(i(-1))) * 400; %EQ(1.26)
	slope		=	ytm -i*400;	%EQ(1.27)
	
%% 27) SDF from t-1 to t
	sdf			= 	beta*(c/c(-1))^-phi*DZ^-phi * 
					(V*DZ^(1-phi)/Vkp(-1))^-alpha;	% real Stochastic discount factor

%---------------------------------------------------------------------------
% New Equations
	Q			=	beta * (V(+1)*DZ(+1)^(1-phi)/Vkp)^(-alpha) * (c(+1)*DZ(+1)/c)^(-phi) * (1/pi(+1));
	Q			=	1/(1+i);
	QR			=	beta * (V(+1)*DZ(+1)^(1-phi)/Vkp)^(-alpha) * (c(+1)*DZ(+1)/c)^(-phi);
	QG			=	(beta/zeta) * (V(+1)*DZ(+1)^(1-phi)/Vkp)^(-alpha) * (c(+1)*DZ(+1)/c)^(-phi) * (Y(+1)/Y*DZ(+1));

	RR			=	1/QR;
	RG 			=	(Y(+1)/Y*DZ(+1))/QG;		%  Expected Return of NGDP
	GP			=	log(RG/RR);
	IP			=	log(exp(i)/RR/pi(+1));
	
	Q*b + QG*bG + T-G = b(-1)/DZ(-1)/pi + (Y/Y(-1))*bG(-1);
	TD			=	 b(-1)/DZ(-1)/pi + (Y/Y(-1))*bG(-1);
	log(T*YBar/TBar/Y)	=	Gamma1*(2/3.141592)*atan(Gamma2*log(TD*YBar/(DoYss*(10^((1-eta)/eta)*(LMax/3)))/Y)) + log(G*YBar/GBar/Y);
	bG			=	omegaG/(1-omegaG) * (Q/QG) * b;
	PSoY		=	(T-G)/Y;
	DoY			=	TD/Y;
end;


initval;
	A			= ABar;
	c			= ((10^((1-eta)/eta)*(LMax/3))-(0.17 * (10^((1-eta)/eta)*(LMax/3)))-((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3)))));
	Disp		= 1;
	DZ			= DZBar;
	G			= (0.17 * (10^((1-eta)/eta)*(LMax/3)));
	i			= (log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES)))));
	ir			= (log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES)))));
	L			= (LMax/3);
	Mc			= (1/(1+theta));
	p0			= 1;
	pi			= 1;
	piavg		= 1;
	pricebond	= (1/(1-(betatilde/DZBar^(-(1/IES)))*consoldelta*DZBar^-(1/IES)));
	pricebondrn	= (1 /(1-consoldelta*exp(-(log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES))))))));
	V			= ((((10^((1-eta)/eta)*(LMax/3))-(0.17 * (10^((1-eta)/eta)*(LMax/3)))-((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3)))))^(1-(1/IES))/(1-(1/IES))+((((1/(1+theta))*eta*10^((1-eta)/eta))*((10^((1-eta)/eta)*(LMax/3))-(0.17 * (10^((1-eta)/eta)*(LMax/3)))-((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3)))))^(-1/IES)*(LMax-(LMax/3))^((LMax-(LMax/3))/Frisch)))*(LMax-(LMax/3))^(1-((LMax-(LMax/3))/Frisch))/(1-((LMax-(LMax/3))/Frisch))) / (1-(betatilde/DZBar^(-(1/IES)))*DZBar^(1-(1/IES))));
	Valphaexp	= 1;
	Vkp			= (((((10^((1-eta)/eta)*(LMax/3))-(0.17 * (10^((1-eta)/eta)*(LMax/3)))-((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3)))))^(1-(1/IES))/(1-(1/IES))+((((1/(1+theta))*eta*10^((1-eta)/eta))*((10^((1-eta)/eta)*(LMax/3))-(0.17 * (10^((1-eta)/eta)*(LMax/3)))-((delta + DZBar - 1) * (10 * (10^((1-eta)/eta)*(LMax/3)))))^(-1/IES)*(LMax-(LMax/3))^((LMax-(LMax/3))/Frisch)))*(LMax-(LMax/3))^(1-((LMax-(LMax/3))/Frisch))/(1-((LMax-(LMax/3))/Frisch))) / (1-(betatilde/DZBar^(-(1/IES)))*DZBar^(1-(1/IES))))*DZBar^(1-(1/IES)));
	wreal		= ((1/(1+theta))*eta*10^((1-eta)/eta));
	Y			= (10^((1-eta)/eta)*(LMax/3));
	ytm			= (log(consoldelta*(1/(1-(betatilde/DZBar^(-(1/IES)))*consoldelta*DZBar^-(1/IES)))/((1/(1-(betatilde/DZBar^(-(1/IES)))*consoldelta*DZBar^-(1/IES)))-1)) *400);
	ytmrn		= (log(consoldelta*(1 /(1-consoldelta*exp(-(log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES))))))))/((1 /(1-consoldelta*exp(-(log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES))))))))-1)) *400);
	zd			= ((10^((1-eta)/eta)*(LMax/3)) / (1-(betatilde/DZBar^(-(1/IES)))*xi*DZBar^-(1/IES)));
	zn			= ((10^((1-eta)/eta)*(LMax/3)) / (1-(betatilde/DZBar^(-(1/IES)))*xi*DZBar^-(1/IES)));
	termprem  	= 100 * ((log(consoldelta*(1/(1-(betatilde/DZBar^(-(1/IES)))*consoldelta*DZBar^-(1/IES)))/((1/(1-(betatilde/DZBar^(-(1/IES)))*consoldelta*DZBar^-(1/IES)))-1)) *400)  - (log(consoldelta*(1 /(1-consoldelta*exp(-(log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES))))))))/((1 /(1-consoldelta*exp(-(log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES))))))))-1)) *400) );
	slope		= (log(consoldelta*(1/(1-(betatilde/DZBar^(-(1/IES)))*consoldelta*DZBar^-(1/IES)))/((1/(1-(betatilde/DZBar^(-(1/IES)))*consoldelta*DZBar^-(1/IES)))-1)) *400) - 400*(log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES)))));
	sdf			= (betatilde/DZBar^(-(1/IES)))*DZBar^-(1/IES);
%-----------------------------------------------------
	Q			= ((betatilde/DZBar^(-(1/IES)))*DZBar^(-(1/IES)));
	QR			= (betatilde/DZBar^(-(1/IES)))*DZBar^(-1/IES);
	QG			= ((betatilde/DZBar^(-(1/IES)))/zeta*DZBar^(1-(1/IES)));
	b			= ((DoYss*(10^((1-eta)/eta)*(LMax/3)))*(DZBar^-1 + omegaG*zeta/DZBar/(1-omegaG))^-1);
	bG			= (omegaG*((betatilde/DZBar^(-(1/IES)))*DZBar^(-(1/IES)))/(1-omegaG)/((betatilde/DZBar^(-(1/IES)))/zeta*DZBar^(1-(1/IES))));
	T			= ((DoYss*(10^((1-eta)/eta)*(LMax/3))) + (0.17 * (10^((1-eta)/eta)*(LMax/3))) - ((betatilde/DZBar^(-(1/IES)))*DZBar^(-(1/IES)))*((DoYss*(10^((1-eta)/eta)*(LMax/3)))*(DZBar^-1 + omegaG*zeta/DZBar/(1-omegaG))^-1)*(1-omegaG)^-1);
	TD			= (DoYss*(10^((1-eta)/eta)*(LMax/3)));
	RR			= DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES)));
	RG			= ((betatilde/DZBar^(-(1/IES)))^-1 * DZBar^(1/IES) * zeta);
	GP			= (((betatilde/DZBar^(-(1/IES)))^-1 * DZBar^(1/IES) * zeta) - 1) - (log(DZBar^(1/IES)/(betatilde/DZBar^(-(1/IES)))));
	IP			= 0;
	
	PSoY		= (((DoYss*(10^((1-eta)/eta)*(LMax/3))) + (0.17 * (10^((1-eta)/eta)*(LMax/3))) - ((betatilde/DZBar^(-(1/IES)))*DZBar^(-(1/IES)))*((DoYss*(10^((1-eta)/eta)*(LMax/3)))*(DZBar^-1 + omegaG*zeta/DZBar/(1-omegaG))^-1)*(1-omegaG)^-1)-(0.17 * (10^((1-eta)/eta)*(LMax/3))))/Y;
	DoY			= (DoYss*(10^((1-eta)/eta)*(LMax/3)))/(10^((1-eta)/eta)*(LMax/3));
end;

steady (solve_algo=2);                      
model_diagnostics ;

shocks;
	var epsa	= sigmaa^2;
	var epsg	= sigmag^2;
	var epsi	= sigmai^2;
	var epsz	= sigmaz^2;
	var epspistar	= sigmapistar^2;

end;

@#define ABC = IRF
@#if IRF
	stoch_simul(order = 1, pruning, nograph, noprint, irf=0);
@#else
	stoch_simul(order = 2, pruning, nograph, irf=80);
@#endif