
%-------------------------------------------------------------------------------------------------------------------------------------
% One country RBC model with flexible labour (hours worked) and government spending.  
% 
% It is a version of
% Eric M. Leeper, Michael Plante, Nora Traum (2009), Dynamics of Fiscal financing in the United States,  Journal of Econometrics
%
% U(c,n) = BETA * u_b * [ (c(t) - hc(t-1)) ** (1-gamma)/(1-gamma) -  u_l* n**(1-KAPPA)/(1-KAPPA) ]
%
% The paper is modified to allow for either IAC or CAC
%  investment adjustment costs as in Christiano Eichenbaum and Evans (2005)
%  capital adjustment costs as in Backus, Routledge and Zin (2007)
%
% h subscript stands for home country
% UPPERCASE is a parameter or steady state value
% lowercase variables are in logarithms
%
% In this version I alter the original leeper file "leeperorig.mod" model equations:
% (A.4)and multiply the errors on the shocks by the variances
% --------------------------------------------------------------------------------------------------------------------------------------



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

var uc_h ul_h y_h c_h i_h k_h q_i_h r_k_h l_h w_h r_h b_h v_h delv_h g_h 
tk_h tl_h tc_h tau_k_h tau_l_h tau_c_h z_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 
i_h_obs c_h_obs l_h_obs tc_h_obs tl_h_obs z_h_obs tk_h_obs g_h_obs b_h_obs ; 

%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% 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 BETA GAMMA H KAPPA DELTA0 DELTA1 DELTA2 ALPHA  CHI ZETA UPS
G_Y_bar B_Y_bar TAU_K_h_bar TAU_L_h_bar TAU_C_h_bar
Y_h_bar B_h_bar Z_h_bar G_h_bar
RHO_B RHO_L RHO_I RHO_G RHO_A RHO_TK RHO_TL RHO_TC RHO_Z   
SIGMA_B_H SIGMA_L_H SIGMA_I_H SIGMA_A_H SIGMA_G_H SIGMA_TK_H  SIGMA_TL_H SIGMA_TC_H SIGMA_Z_H
EPSI_G EPSI_K EPSI_L EPSI_Z 
GAMMA_G GAMMA_K GAMMA_L GAMMA_Z PHI_KL PHI_KC PHI_LC ;



%--------------------------------------------------------------------------
% Calibration:
%--------------------------------------------------------------------------

%Discount factor
BETA       = 0.99;

DELTA0     = 0.025;

%Capital share
ALPHA      = 0.30;

%Steady state ratio of governemnt spending to GDP
G_Y_bar    = .0922;

%Steady state ratio of governemnt debt to GDP
B_Y_bar    = .3396;

%Steady state capital tax rate
TAU_K_h_bar = 0.184;

%Steady state labour tax rate
TAU_L_h_bar = 0.223;

// Steady state consumption tax rate
TAU_C_h_bar = 0.028;

%Parameter values based on Table 2 in Leeper et al (2009).  I use the mean of the posterior distribution under column titled "All Adjust".

%Risk aversion
GAMMA      = 2.7;

%Iinverse Frish elasticity
KAPPA      = 1.9;

%Habit formation
H          = 0.50;

%Investment adjustment cost parameters and is 0 <= CHI 
CHI      = 5.5;

%Capital adjustment cost parameters and 0 < UPS <= 1
UPS      = 0.500; 

%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     = 1.0; 

%Capital utilization
DELTA2      = 0.29;

%Government spending debt coefficient
GAMMA_G     = 0.23;

%Capital tax rate debt coefficient
GAMMA_K     = 0.39;

%Labour tax rate debt coefficient
GAMMA_L     = 0.049;

%Transfers debt coefficient
GAMMA_Z     = 0.5;

%Capital tax rate output coefficient
EPSI_K     = 1.7;

%Labour tax rate output coefficient
EPSI_L     = 0.36;

%Government spending output coefficient
EPSI_G     = 0.34;

%Transfers output coefficient
EPSI_Z     = 0.13;

%Capital-Labour co-term
PHI_KL     = 0.19;

%Capital-Consumption co-term
PHI_KC     = 0.024;

%Labour-Consumption co-term
PHI_LC     = -0.028;

%Parameters for the exogenous shocks 

RHO_A        = 0.96;
SIGMA_A_H    = 0.0062;
RHO_B        = 0.66;
SIGMA_B_H    = 0.07;
RHO_L        = 0.99;
SIGMA_L_H    = 0.028;
RHO_I        = 0.55;
SIGMA_I_H    = 0.064;
RHO_G        = 0.97;
SIGMA_G_H    = 0.031;
RHO_TK       = 0.93;
SIGMA_TK_H   = 0.044;
RHO_TL       = 0.97;
SIGMA_TL_H   = 0.03;
RHO_TC       = 0.93;
SIGMA_TC_H   = 0.04;
RHO_Z        = 0.94;
SIGMA_Z_H    = 0.034;

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

R_h_bar   = 1/BETA;
R_K_h_bar = (1/BETA + 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-G_Y_bar)-DELTA0))^(-GAMMA)*(1-TAU_L_h_bar)*(1-ALPHA)/(1+TAU_C_h_bar))/ (R_K_h_bar/ALPHA)^((KAPPA+ALPHA)/(1-ALPHA))  ) ^(1/(KAPPA+GAMMA));
A1   =    (((1-H)*((R_K_h_bar/ALPHA)*(1-G_Y_bar)-DELTA0))^(-GAMMA)*(1-TAU_L_h_bar)*(1-ALPHA)/(1+TAU_C_h_bar));
A2 = (R_K_h_bar/ALPHA)^((KAPPA+ALPHA)/(1-ALPHA));
A3 = (A1/A2)^(1/(KAPPA+GAMMA));
//disp('Steady state A1');      disp(A1);
//disp('Steady state A2');      disp(A2);
//disp('Steady state A3');      disp(A3);

Y_h_bar   =  (R_K_h_bar/ALPHA)                       *K_h_bar;

C_h_bar   = ((R_K_h_bar/ALPHA)*(1-G_Y_bar)-DELTA0)   *K_h_bar;
I_h_bar   = DELTA0                                   *K_h_bar;
L_h_bar   =  (R_K_h_bar/ALPHA)^(1/(1-ALPHA))         *K_h_bar;
G_h_bar   =  G_Y_bar   *   Y_h_bar;
B_h_bar   =  B_Y_bar   *   Y_h_bar;
Z_h_bar   =  TAU_K_h_bar*ALPHA*Y_h_bar + TAU_L_h_bar*(1-ALPHA)*Y_h_bar + TAU_C_h_bar*C_h_bar - Y_h_bar * (G_Y_bar + ((1-BETA)/BETA)*B_Y_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; 

Uc_h_bar = (C_h_bar-H*C_h_bar) ^(-GAMMA) ;
Ul_h_bar = -(L_h_bar ^ KAPPA) ;
W_h_bar   = (1-ALPHA)* Y_h_bar / L_h_bar ;
DELV_h_bar   = DELTA0 ;
end;

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 z_h');      disp(Z_h_bar);
disp('Steady state r_h');      disp(R_h_bar);
disp('Steady state r_k_h');    disp(R_K_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);

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

model;

%MU of consumption
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;

%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))   +  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 ) )  +   (1-exp(tau_k_h(+1)))*ALPHA * exp(y_h(+1))/exp(k_h)  ) / 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
exp(c_h) + exp(i_h) + exp(g_h)  = exp(y_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.8):  Government budget 
exp(b_h) + exp(tk_h) + exp(tl_h) + exp(tc_h)  = exp(r_h(-1))*exp(b_h(-1)) + exp(g_h) + exp(z_h) ;

%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

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

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

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

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

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

% Exogenous shocks

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

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

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

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

u_g_h  = RHO_G  * u_g_h(-1)  + SIGMA_G_H  * e_u_g_h;
u_tk_h = RHO_TK * u_tk_h(-1) + SIGMA_TK_H * e_u_tk_h;
u_tl_h = RHO_TL * u_tl_h(-1) + SIGMA_TL_H * e_u_tl_h;
u_tc_h = RHO_TC * u_tc_h(-1) + SIGMA_TC_H * e_u_tc_h;
u_z_h  = RHO_Z  * u_z_h(-1)  + SIGMA_Z_H  * e_u_z_h;

i_h_obs  = (i_h);  
//- log (0.4623);
c_h_obs  = (c_h); 
//- log (1.9448);
l_h_obs  = (l_h); 
//- log (1.1535);
tc_h_obs = (tc_h); 
//- log (0.0545);
tl_h_obs = (tl_h);
// - log (0.4139);
z_h_obs  = (z_h); 
// - log (0.3612);
tk_h_obs = (tk_h);
// - log (0.1464);
g_h_obs  = (g_h);
// - log (0.2445);
b_h_obs  = (b_h);
// - log (0.9005);

end;


%--------------------------------------------------------------------------
% Observables:
%--------------------------------------------------------------------------

varobs 
i_h_obs c_h_obs l_h_obs tc_h_obs tl_h_obs z_h_obs tk_h_obs g_h_obs b_h_obs; 

%shocks;
%var e_u_b_h  = 1^2;
%var e_u_l_h  = 1^2;
%var e_u_i_h  = 1^2;
%var e_u_a_h  = 1^2;
%var e_u_g_h  = 1^2;
%var e_u_z_h  = 1^2;
%var e_u_tk_h = 1^2;
%var e_u_tl_h = 1^2;
%var e_u_tc_h = 1^2;
%end;

%--------------------------------------------------------------------------
% Computation
%--------------------------------------------------------------------------

initval;
uc_h    = (Uc_h_bar);
ul_h    = (Ul_h_bar);
y_h     = log(Y_h_bar);
c_h     = log(C_h_bar);
i_h     = log(I_h_bar);
l_h     = log(L_h_bar);
k_h     = log(K_h_bar);
w_h     = log(W_h_bar);
g_h     = log(G_h_bar);
b_h     = log(B_h_bar);
z_h     = log(Z_h_bar);
r_h     = log(R_h_bar);
r_k_h   = log(R_K_h_bar);
tk_h    = log(TK_h_bar);
tl_h    = log(TL_h_bar);
tc_h    = log(TC_h_bar);
tau_k_h = log(TAU_K_h_bar);
tau_l_h = log(TAU_L_h_bar);
tau_c_h = log(TAU_C_h_bar);
v_h     = 0;
delv_h  = log(DELV_h_bar);
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;

%--------------------------------------------------------------------------
% Priors
%--------------------------------------------------------------------------

estimated_params; 

GAMMA,  gamma_pdf, 1.75, 0.5;
KAPPA,  gamma_pdf, 2, 0.5;
DELTA2, gamma_pdf, 0.7, 0.5; 
CHI,    gamma_pdf, 5, 0.25;
H,      beta_pdf, 0.5, 0.2; 

RHO_B,  beta_pdf, 0.7, 0.2;
RHO_L,  beta_pdf, 0.7, 0.2;
RHO_I,  beta_pdf, 0.7, 0.2;
RHO_G,  beta_pdf, 0.7, 0.2;
RHO_A,  beta_pdf, 0.7, 0.2;
RHO_TK, beta_pdf, 0.7, 0.2;
RHO_TL, beta_pdf, 0.7, 0.2;
RHO_TC, beta_pdf, 0.7, 0.2;
RHO_Z,  beta_pdf, 0.7, 0.2;

SIGMA_B_H,   inv_gamma_pdf, 1, 4; 
SIGMA_L_H,   inv_gamma_pdf, 1, 4; 
SIGMA_I_H,   inv_gamma_pdf, 1, 4; 
SIGMA_A_H,   inv_gamma_pdf, 1, 4; 
SIGMA_G_H,   inv_gamma_pdf, 1, 4; 
SIGMA_TK_H,  inv_gamma_pdf, 1, 4; 
SIGMA_TL_H,  inv_gamma_pdf, 1, 4; 
SIGMA_TC_H,  inv_gamma_pdf, 1, 4; 
SIGMA_Z_H,   inv_gamma_pdf, 1, 4; 

stderr e_u_b_h,  inv_gamma_pdf, 0.01, inf;
stderr e_u_l_h,  inv_gamma_pdf, 0.01, inf;
stderr e_u_i_h,  inv_gamma_pdf, 0.01, inf;
stderr e_u_a_h,  inv_gamma_pdf, 0.01, inf;
stderr e_u_g_h,  inv_gamma_pdf, 0.01, inf;
stderr e_u_tk_h, inv_gamma_pdf, 0.01, inf;
stderr e_u_tl_h, inv_gamma_pdf, 0.01, inf;
stderr e_u_tc_h, inv_gamma_pdf, 0.01, inf;
stderr e_u_z_h,  inv_gamma_pdf, 0.01, inf;

GAMMA_G,     gamma_pdf, 0.4, 0.20; 
GAMMA_K,     gamma_pdf, 0.4, 0.20; 
GAMMA_L,     gamma_pdf, 0.4, 0.20;  
GAMMA_Z,     gamma_pdf, 0.4, 0.20;  
EPSI_K,      gamma_pdf,   1, 0.30;
EPSI_L,      gamma_pdf, 0.5, 0.25;
EPSI_G,      gamma_pdf, 0.07,0.05; 
EPSI_Z,      gamma_pdf, 0.20,0.10; 

PHI_KL,     normal_pdf, 0.25, 0.10; 
PHI_KC,     normal_pdf, 0.05, 0.10;
PHI_LC,     normal_pdf, 0.05, 0.10;

end;

%--------------------------------------------------------------------------
% Estimation
%--------------------------------------------------------------------------


estimation(datafile =transf_d,nobs=193,first_obs=1, mh_replic=2000, mh_nblocks=2, mh_drop=0.45, mh_jscale=0.8, mode_compute=4);


