% Running the open economy with open economy parameter estimates 

%--------------------------------------------------------------------------
% 40 Endogenous variables
%--------------------------------------------------------------------------

var uc_h ul_h y_h c_h i_h k_h q_i_h r_k_h R_h b_h l_h w_h v_h delv_h  
g_h z_h tk_h tl_h tc_h     tau_k_h tau_l_h tau_c_h 
u_b_h u_l_h u_i_h u_a_h u_g_h u_tk_h u_tl_h u_tc_h u_z_h 
b_f  b  by_h by_f by    tb_h tby_h   cay_h risk_prem U welfare;

%--------------------------------------------------------------------------
% Endogenous variable definitions
%--------------------------------------------------------------------------

% uc_h            : marginal utility of consumption
% ul_h            : marginal utility of hours worked
% y_h             : log of home output
% c_h             : log of home consumption
% i_h             : log of home investment
% k_h             : log of home capital stock
% q_i_h           : log of home investment price
% l_h             : log of home labour
% w_h             : log of the home real wage
% b_h             : log of the number of home bonds
% v_h             : log of the utilization level
% delv_h          : log of the depreciation rate

% g_h             : log of home government consumption
% tk_h            : log of home total capital tax revenue
% tl_h            : log of home total labour tax revenue
% tc_h            : log of home total consumption tax revenue
% tau_k_h         : log of home capital tax rate
% tau_l_h         : log of home labour tax rate
% tau_c_h         : log of home consumption tax rate
% z_h             : log of home government lump sum transfers
% b_h			  : log of home debt 
% b_f		 	  : log of foreign debt 
% by_h            : log of home debt to output ratio
% by_f            : log of foreign debt to output ratio 
% by			  : log of total debt to output ratio
% b 			  : log of total debt 
% R_h			  : log of debt elastic interest rate 
% tb_h            : log of home trade balance 
% tby_h   		  : log of home trade balance to output ratio
% cay_h 		  : log of home current account to output ratio 

% u_b_h           : log of home general preference shock
% u_l_h        	  : log of home labour supply shock
% u_i_h        	  : log of home investment specific shock
% u_a_h           : log of home technology
% u_g_h        	  : log of home government spending shock
% u_tk_h          : log of home capital tax rate shock
% u_tl_h       	  : log of home labour tax rate shock
% u_tc_h       	  : log of home consumption tax rate shock
% u_z_h        	  : log of home government lump sum transfers shock

%--------------------------------------------------------------------------
% 9 Exogenous shocks (error terms)
%--------------------------------------------------------------------------

varexo e_u_b_h e_u_l_h e_u_i_h e_u_a_h e_u_g_h e_u_tk_h e_u_tl_h e_u_tc_h 
e_u_z_h;

%--------------------------------------------------------------------------
% Parameters
%--------------------------------------------------------------------------

parameters BETA GAMMA H DELTA0 DELTA1 DELTA2 ALPHA  CHI ZETA UPS PSI2 R_w KAPPA
TK_h_bar TL_h_bar TC_h_bar L_h_bar GY_h_bar BY_h_bar BY_f_bar  
G_h_bar  Y_h_bar B_bar B_h_bar B_f_bar Z_h_bar TAU_K_h_bar TAU_L_h_bar TAU_C_h_bar 

EPSI_G EPSI_K EPSI_L EPSI_Z GAMMA_G GAMMA_K GAMMA_L GAMMA_Z PHI_KL PHI_KC PHI_LC    

RHO_B RHO_L RHO_I RHO_G RHO_A RHO_TK RHO_TL RHO_TC RHO_Z   
SIGMA_A_H SIGMA_B_H SIGMA_L_H SIGMA_I_H SIGMA_G_H SIGMA_Z_H SIGMA_TK_H SIGMA_TL_H SIGMA_TC_H ;

%--------------------------------------------------------------------------
% Calibration
%--------------------------------------------------------------------------

% We set BETA. This implies a value for steady state R.  
% From here we can work out steady state Rk, K, Y, I and finally C (in that order).

% Parameter values are from Letendre (2007) and SGU (2003)

BETA 			= 0.993;
R_w  			= 1/BETA;
PSI2 			= 0.021;
% 0.001 Letendre (2007)
KAPPA     		= 0.57;
DELTA0    	 	= 0.025;
ALPHA     	    = 0.32;
GAMMA      		= 1.27;
H          		= 0.6;
CHI        		= 5.09;
UPS        		= 0.500; 
ZETA       		= 1.0; 
DELTA2    	 	= 0.83;

GY_h_bar    	= .23;
BY_h_bar    	=0.56;
BY_f_bar    	= .10;

% Parameter values come from the mean of the respective series 
TAU_K_h_bar 	= 0.22;
TAU_L_h_bar 	= 0.11;
TAU_C_h_bar 	= 0.099;

% Parameter values from Leeper (2010)
GAMMA_G   	 	= 0.30;
GAMMA_K    		= 0.02;
GAMMA_L    		= 0.03;
GAMMA_Z    		= 0.35;
EPSI_K     		= 1.93;
EPSI_L     		= 0.54;
EPSI_G     		= 0.05;
EPSI_Z     		= 0.15;
PHI_KL     		= 0.11;
PHI_KC     		= -0.05;
PHI_LC     		= -0.04;

%--------------------------------------------------------------------------
% Definitions of parameters 
%--------------------------------------------------------------------------

% BETA          : discount factor 
% R_w           : world interest rate 
% PSI2          : coefficient for debt elastic interest rate 
% KAPPA      	: inverse Frish elasticity 
% DELTA0     	:
% ALPHA 		: capital share 
% GY_h_bar 		: SS ratio of government spending to GDP
% BY_h_bar    	: SS debt, Ratio of government (public held at home) debt to GDP 
% BY_f_bar 		: Ratio of government (public held abroad) debt to GDP 
% TAU_K_h_bar	: Steady state capital tax rate
% TAU_L_h_bar	: Steady state labor tax rate
% TAU_C_h_bar 	: Steady state consumption tax rate 

% GAMMA         : Risk aversion
% H    		    : Habit formation 
% CHI 		    : Investment adjustment cost parameters and is 0 <= CHI 
% UPS 		    : Capital adjustment cost parameters and 0 < UPS <= 1

% Capital and/or Investment adjustment cost parameters and 0 <= ZETA <= 1
% ZETA =1 Pure investment adjustment costs as in Christiano and Eichenbaum %(2005)
% ZETA =0 Pure capital adjustment costs as in Backus, Routledge and Zin (2007)

% ZETA     		: Capital and/or Investment adjustment cost parameters 
% DELTA2   		: Capital utilization 
% GAMMA_G     	: Government spending debt coefficient
% GAMMA_K     	: Capital tax rate debt coefficient 
% GAMMA_L 		: Labor tax rate debt coefficient
% GAMMA_Z     	: Transfers debt coefficient 
% EPSI_K   		: Capital tax rate output coefficient
% EPSI_L		: Labor tax rate output coefficient
% EPSI_G     	: Government spending output coefficient	
% EPSI_Z     	: Transfers output coefficient 
% PHI_KL     	: Capital-Labor co-term
% PHI_KC     	: Capital-Consumption co-term 
% PHI_LC     	: Labor-Consumption co-term 

%--------------------------------------------------------------------------
% Parameters for the exogenous shocks
%--------------------------------------------------------------------------

RHO_A        = 0.99;
SIGMA_A_H    = 0.59;

RHO_B        = 0.91;
SIGMA_B_H    = 0.76;

RHO_L        = 0.99;
SIGMA_L_H    = 2.15;

RHO_I        = 0.51;
SIGMA_I_H    = 1.61;

RHO_G        = 0.98;
SIGMA_G_H    = 2.12;

RHO_TK       = 0.98;
SIGMA_TK_H   = 2.23;

RHO_TL       = 0.95;
SIGMA_TL_H   = 1.95;

RHO_TC       = 0.86;
SIGMA_TC_H   = 2.78;

RHO_Z        = 0.99;
SIGMA_Z_H    = 3.05;

%--------------------------------------------------------------------------
% Steady states
%--------------------------------------------------------------------------

R_h_bar     = 1/BETA;
R_K_h_bar   = (R_h_bar + DELTA0 - 1) / (1-TAU_K_h_bar);
DELTA1      = R_K_h_bar * (1-TAU_K_h_bar);

K_h_bar     = ( ( ((1-H)*( (R_K_h_bar/ALPHA)*(1-GY_h_bar-BY_f_bar*(R_h_bar-1))-DELTA0) ))^(-GAMMA)* ((1-TAU_L_h_bar)*(1-ALPHA)/(1+TAU_C_h_bar))*(R_K_h_bar/ALPHA)^((ALPHA+KAPPA)/(ALPHA-1)) ) ^ (1/(GAMMA+KAPPA)) ;                      

L_h_bar     = (R_K_h_bar/ALPHA)^(1/(1-ALPHA)) * K_h_bar;
Y_h_bar     = K_h_bar^ALPHA  * L_h_bar^(1-ALPHA);
I_h_bar     = DELTA0 * K_h_bar;
G_h_bar     = GY_h_bar   *   Y_h_bar;
B_h_bar     = BY_h_bar   *   Y_h_bar;
W_h_bar     = (1-ALPHA)* Y_h_bar / L_h_bar ;

B_f_bar     = BY_f_bar   *   Y_h_bar;
C_h_bar     = Y_h_bar - I_h_bar - G_h_bar - (R_h_bar-1)*B_f_bar;

TB_h_bar    = Y_h_bar - (C_h_bar+I_h_bar+G_h_bar);
TBY_h_bar   = TB_h_bar/Y_h_bar;
CAY_h_bar   = (-(R_h_bar-1) * B_f_bar + TB_h_bar) / Y_h_bar;
BY_h_bar    =  B_h_bar/Y_h_bar;
BY_f_bar    =  B_f_bar/Y_h_bar;
B_bar       =  B_h_bar + B_f_bar;
BY_bar      =  BY_h_bar + BY_f_bar;

TK_h_bar    =  TAU_K_h_bar *  ALPHA    *   Y_h_bar;
TL_h_bar    =  TAU_L_h_bar * (1-ALPHA) *   Y_h_bar; 
TC_h_bar    =  TAU_C_h_bar *               C_h_bar; 
Z_h_bar     =  TK_h_bar + TL_h_bar + TC_h_bar - G_h_bar - (R_h_bar-1)*B_bar;

Uc_h_bar    = (C_h_bar-H*C_h_bar) ^(-GAMMA) ;
Ul_h_bar    = -(L_h_bar ^ KAPPA) ;

disp('Steady state DELTA0');     disp(DELTA0);
disp('Steady state DELTA1');     disp(DELTA1);
disp('Steady state KAPPA');      disp(KAPPA);

disp('=====================');    
disp('LEVELS');
disp('Steady state R_h_bar');    disp(R_h_bar);
disp('Steady state R_K_h_bar');  disp(R_K_h_bar);
disp('Steady state Uc_h_bar');   disp(Uc_h_bar);
disp('Steady state Ul_h_bar');   disp(Ul_h_bar);
disp('Steady state y_h');       disp(Y_h_bar);
disp('Steady state c_h');       disp(C_h_bar);
disp('Steady state i_h');       disp(I_h_bar);
disp('Steady state l_h');       disp(L_h_bar);
disp('Steady state w_h');       disp(W_h_bar);
disp('Steady state k_h');       disp(K_h_bar);
disp('Steady state g_h');       disp(G_h_bar);
 
disp('Steady state b_h');       disp(B_h_bar);
disp('Steady state b_f');       disp(B_f_bar);
disp('Steady state by_h');      disp(BY_h_bar);
disp('Steady state by_f');      disp(BY_f_bar);
disp('Steady state b');         disp(B_bar);
disp('Steady state by');        disp(BY_bar);

disp('Steady state z_h');       disp(Z_h_bar);
disp('Steady state tb_h');      disp(TB_h_bar);
disp('Steady state tby_h');     disp(TBY_h_bar);
disp('Steady state cay_h');     disp(CAY_h_bar);
disp('Steady state tk_h');      disp(TK_h_bar);
disp('Steady state tl_h');      disp(TL_h_bar);
disp('Steady state tc_h');      disp(TC_h_bar);
disp('Steady state tau_k_h');   disp(TAU_K_h_bar);
disp('Steady state tau_l_h');   disp(TAU_L_h_bar);
disp('Steady state tau_c_h');   disp(TAU_C_h_bar);
disp('=====================');    
disp('  ');    

disp('=====================');    
disp('Convert to logs');
R_K_h_bar   = log(R_K_h_bar);
Y_h_bar     = log(Y_h_bar);
C_h_bar     = log(C_h_bar);
I_h_bar     = log(I_h_bar);
L_h_bar     = log(L_h_bar);
K_h_bar     = log(K_h_bar);
W_h_bar     = log(W_h_bar);
G_h_bar     = log(G_h_bar);
TK_h_bar    = log(TK_h_bar);
TL_h_bar    = log(TL_h_bar);
TC_h_bar    = log(TC_h_bar);
TAU_K_h_bar = log(TAU_K_h_bar);
TAU_L_h_bar = log(TAU_L_h_bar);
TAU_C_h_bar = log(TAU_C_h_bar);

disp('Steady state R_K_h_bar');  disp(R_K_h_bar);
disp('Steady state y_h');        disp(Y_h_bar);
disp('Steady state c_h');        disp(C_h_bar);
disp('Steady state i_h');        disp(I_h_bar);
disp('Steady state l_h');        disp(L_h_bar);
disp('Steady state w_h');        disp(W_h_bar);
disp('Steady state k_h');        disp(K_h_bar);
disp('Steady state g_h');        disp(G_h_bar);
disp('Steady state tk_h');       disp(TK_h_bar);
disp('Steady state tl_h');       disp(TL_h_bar);
disp('Steady state tc_h');       disp(TC_h_bar);
disp('Steady state tau_k_h');    disp(TAU_K_h_bar);
disp('Steady state tau_l_h');    disp(TAU_L_h_bar);
disp('Steady state tau_c_h');    disp(TAU_C_h_bar);
disp('=====================');   

%--------------------------------------------------------------------------
% Model
%--------------------------------------------------------------------------

model;

U = (u_b_h) * ( ( c_h - H*c_h(-1) )^(1-GAMMA)/(1-GAMMA) - (u_l_h)* l_h^(1+KAPPA)/(1+ KAPPA));

welfare = U + BETA*welfare(+1);

risk_prem = PSI2*(exp(b_f - B_f_bar) - 1);

exp(R_h)  = R_w + risk_prem ;

by_h   = b_h/exp(y_h);
by_f   = b_f/exp(y_h);
b      = b_h + b_f;
by     = by_h + by_f;

tb_h   = exp(y_h) - exp(c_h) - exp(i_h) - exp(g_h)  ;
tby_h  = tb_h/exp(y_h);
cay_h  = tby_h - (exp(R_h(-1))-1)*b_f(-1)/exp(y_h);

% MU of consumption (Habit is not internalized)
uc_h = exp(u_b_h)*( exp(c_h) - H*exp(c_h(-1)) ) ^(-GAMMA) ;

% MU of hours worked
ul_h = -exp(u_b_h)*exp(u_l_h)*exp(l_h) ^ KAPPA ;

% Equation 6:  Endogenous depreciation rate of the capital stock
exp(delv_h)  = DELTA0 + DELTA1*(exp(v_h) - 1) + (DELTA2/2)*(exp(v_h) - 1)^2;

% Equation 9: marginal product of labour
exp(w_h) = (1-ALPHA) * exp(y_h)/exp(l_h)  ;

% Marginal product of capital
exp(r_k_h) = ALPHA   * exp(y_h)/(exp(v_h)*exp(k_h(-1))) ;

% Equation (A.1):  Euler equation - intertemporal condition
uc_h/(1+exp(tau_c_h))   = BETA*exp(R_h)*uc_h(+1)/(1+exp(tau_c_h(+1)));

% Equation (A.2):  Intratemporal condition
ul_h + uc_h * exp(w_h) * (1-exp(tau_l_h))/(1+exp(tau_c_h)) = 0;

% Equation (A.3):  marginal conditions for capital 
exp(q_i_h) = ( exp(q_i_h(+1))*(1-exp(delv_h(+1)))  +   (1-exp(tau_k_h(+1)))*ALPHA * exp(y_h(+1))/exp(k_h) 
+  exp(q_i_h(+1))*( ( (1-ZETA)*(1-UPS)*DELTA0/UPS)*(  ( exp(u_i_h(+1))*exp(i_h(+1))/exp(k_h))^UPS * DELTA0^(-UPS) - 1 ) )  ) /exp(R_h) ;

% Equation (A.4):  marginal conditions for intensity of capital
exp(r_k_h)*(1-exp(tau_k_h))   = exp(q_i_h)*( DELTA1 + DELTA2*(exp(v_h) - 1) ) ;

% Equation (A.5):  marginal conditions for investment
1.0 = exp(q_i_h)* ( ZETA*( 1.0 - 0.5*CHI - 1.5*CHI*(exp(u_i_h)*exp(i_h)/exp(i_h(-1)))^2 + 2.0*CHI*(exp(u_i_h)*exp(i_h)/exp(i_h(-1))) ) 
+ (1-ZETA)* ( (exp(u_i_h)*exp(i_h)/exp(k_h(-1)))^(UPS-1) * DELTA0^(1-UPS) )  ) 
+ exp(q_i_h(+1)) *ZETA *CHI *( (exp(u_i_h(+1))*exp(i_h(+1))/exp(i_h))^3 -(exp(u_i_h(+1))*exp(i_h(+1))/exp(i_h))^2 )  / exp(R_h) ;


% Equation (A.6):  final goods market equilibrium LIKE SGU BUT ADJUSTMENT COSTS % ARE PUT IN CAPITAL ACCUMULATION PROCESS 
b_f =   exp(R_h(-1)) * b_f(-1) - exp(y_h) + exp(c_h) + exp(i_h) +  exp(g_h);

%Equation :  Household resource constraint 
(1+exp(tau_c_h))*exp(c_h) + exp(i_h) + b_h = exp(l_h)*exp(w_h)*(1-exp(tau_l_h)) + exp(r_k_h)*(1-exp(tau_k_h))*exp(k_h(-1)) + exp(R_h(-1))*b_h(-1) + z_h ;


% Equation (A.7):  Capital stock accumulation (adjusted to allow for either IAC % and/or CAC)
% Note if ZETA=1 it is the same as Leeper et al (2009).

exp(k_h)   = (1-exp(delv_h))*exp(k_h(-1)) 
+    ZETA *( (1 - 0.5 * CHI * (exp(u_i_h)*exp(i_h)/exp(i_h(-1))-1.0)^2) * exp(i_h) ) 
+ (1-ZETA)*( (  (exp(u_i_h)*exp(i_h)/exp(k_h(-1)))^UPS * DELTA0^(1-UPS) - (1-UPS)* DELTA0 )*(exp(k_h(-1))/UPS) );

% Equation (A.9):  Total capital tax revenue 
exp(tk_h)   = exp(tau_k_h) * ALPHA * exp(y_h) ;

% Equation (A.10):  Total labour tax revenue 
exp(tl_h)   = exp(tau_l_h) * (1-ALPHA) * exp(y_h) ;

% Equation (A.11):  Total consumption tax revenue 
exp(tc_h)   = exp(tau_c_h) * exp(c_h) ;

% Equation (A.12):  Production function
exp(y_h) = exp(u_a_h) * (exp(v_h)*exp(k_h(-1)))^ALPHA * exp(l_h)^(1-ALPHA);


%--------------------------------------------------------------------------
% Fiscal rules where instruments respond to foreign debt
%--------------------------------------------------------------------------


%Equation 11: Government spending policy rule
g_h - G_h_bar         = - EPSI_G*(y_h-Y_h_bar) - GAMMA_G*(b(-1)-B_bar) + (u_g_h)  ;

%Equation 12: Capital tax policy rule
tau_k_h - TAU_K_h_bar =   EPSI_K*(y_h-Y_h_bar) + GAMMA_K*(b(-1)-B_bar) + PHI_KL*u_tl_h + PHI_KC*u_tc_h + (u_tk_h)  ;

%Equation 13: Labour tax policy rule
tau_l_h - TAU_L_h_bar =   EPSI_L*(y_h-Y_h_bar) + GAMMA_L*(b(-1)-B_bar)  + PHI_KL*u_tk_h + PHI_LC*u_tc_h + (u_tl_h)  ;

%Equation 14: Consumption tax policy rule
tau_c_h - TAU_C_h_bar =   PHI_KC*(u_tk_h) + PHI_LC*(u_tl_h) + (u_tc_h)  ;

%Equation 15: Transfers tax policy rule
z_h - Z_h_bar         = - EPSI_Z*(y_h-Y_h_bar) - GAMMA_Z*(b(-1)-B_bar) + (u_z_h)  ;


%--------------------------------------------------------------------------
% Exogenous shocks
%--------------------------------------------------------------------------

% Equation 1: Taste shock
u_b_h = RHO_B * u_b_h(-1) +  e_u_b_h;

% Equation 2: Leisure shock
u_l_h = RHO_L * u_l_h(-1) +  e_u_l_h;

% Equation 5: Utilization shock
u_i_h = RHO_I * u_i_h(-1) +  e_u_i_h;

% Equation 8: Technology shock
u_a_h = RHO_A * u_a_h(-1) +  e_u_a_h;

u_g_h  = RHO_G  * u_g_h(-1)  + e_u_g_h;

u_tk_h = RHO_TK * u_tk_h(-1) + e_u_tk_h;

u_tl_h = RHO_TL * u_tl_h(-1) + e_u_tl_h;

u_tc_h = RHO_TC * u_tc_h(-1) + e_u_tc_h;

u_z_h  = RHO_Z  * u_z_h(-1)  + e_u_z_h;


end;

%--------------------------------------------------------------------------
% Computation/Initial Values
%--------------------------------------------------------------------------

initval;
U       = 0;
welfare = 0;
R_h     = log(R_h_bar);
uc_h    = Uc_h_bar;
ul_h    = Ul_h_bar;
y_h     = Y_h_bar;
c_h     = C_h_bar;
i_h     = I_h_bar;
l_h     = L_h_bar;
k_h     = K_h_bar;
w_h     = W_h_bar;
g_h     = G_h_bar;
b       = B_bar;
b_h     = B_h_bar;
b_f     = B_f_bar;
by      = BY_bar;
by_h    = BY_h_bar;
by_f    = BY_f_bar;
tb_h    = TB_h_bar;
tby_h   = TBY_h_bar;
cay_h   = CAY_h_bar;
risk_prem = 0;
z_h     = Z_h_bar;
r_k_h   = R_K_h_bar;
tk_h    = TK_h_bar;
tl_h    = TL_h_bar;
tc_h    = TC_h_bar;
tau_k_h = TAU_K_h_bar;
tau_l_h = TAU_L_h_bar;
tau_c_h = TAU_C_h_bar;
v_h     = 0;
delv_h  = log(DELTA0);
q_i_h   = 0;
u_b_h   = 0;
u_l_h   = 0;
u_i_h   = 0;
u_a_h   = 0;
u_g_h   = 0;
u_tk_h  = 0;
u_tl_h  = 0;
u_tc_h  = 0;
u_z_h   = 0;

end;

%--------------------------------------------------------------------------
% shocks block
%--------------------------------------------------------------------------

shocks;

var e_u_b_h  = SIGMA_B_H ^2;
var e_u_l_h  = SIGMA_L_H ^2;
var e_u_i_h  = SIGMA_I_H ^2;
var e_u_a_h  = SIGMA_A_H ^2;
var e_u_g_h  = SIGMA_G_H ^2;
var e_u_z_h  = SIGMA_Z_H ^2;
var e_u_tk_h = SIGMA_TK_H ^2;
var e_u_tl_h = SIGMA_TL_H ^2;
var e_u_tc_h = SIGMA_TC_H ^2;

end;

%--------------------------------------------------------------------------
% Simulation
%--------------------------------------------------------------------------

steady(solve_algo=2);
resid;
check;
%model_info;

stoch_simul(order=2,pruning,periods=100000,nograph);


%stoch_simul(order=1,hp_filter=1600,irf=40, conditional_variance_decomposition=1) y_h c_h i_h l_h b b_h b_f w_h r_k_h R_h risk_prem cay_h tby_h 
%g_h z_h tau_k_h tau_l_h tau_c_h tk_h tl_h tc_h;


