%--------------------------------------------------------------------------
% ROT model with distortionary taxes and debt
% All adjust where ROT agents pay no transfer tax 
% Open economy 
% 2 wage rates and labour tax rates 
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% 32 Endogenous variables  
%--------------------------------------------------------------------------
var     y_h c_h c_h_o c_h_r l_h l_h_o l_h_r i_h k_h b_h  w_h_o w_h_r r_h r_k_h
                g_h tau_k_h tau_c_h tau_l_h  z_h 
        u_a_h u_g_h u_tk_h  u_tc_h  u_tl_h  u_z_h      
        tk_h   tc_h   tl_h  b_f; 

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

varexo  e_u_a_h e_u_g_h e_u_tk_h e_u_tc_h e_u_tl_h e_u_z_h;

%--------------------------------------------------------------------------
% 45 parameters
%--------------------------------------------------------------------------
parameters BETA DELTA ALPHA GAMMA KAPPA H THETA PHI 

Y_h_bar G_h_bar B_h_bar Z_h_bar 
TAU_K_h_bar TAU_C_h_bar TAU_L_h_bar 
TK_h_bar TC_h_bar TL_h_bar 
RHO_A RHO_G RHO_TK RHO_TC RHO_TL  RHO_Z 
SIGMA_A_H SIGMA_G_H SIGMA_TK_H SIGMA_TC_H SIGMA_TL_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 R_w PSI2 B_f_bar B_Y_h_bar;

%--------------------------------------------------------------------------
% parameter values
%--------------------------------------------------------------------------

BETA   = 0.99;
DELTA  = 0.025;
ALPHA  = 0.30;
GAMMA  = 5.0;
KAPPA  = -3.050;
H      = 0.50;
THETA  = 0.50;
PHI    = 10.0;

R_w  			= 1/BETA;
PSI2 			= 0.021;

Y_h_bar  = 0.956;
B_h_bar  = 0.324;
B_f_bar  = 0.10;
B_bar  = B_h_bar + B_f_bar;
B_Y_h_bar= B_bar * Y_h_bar;
G_h_bar  = 0.088;
Z_h_bar  = 0.132;

// Steady state capital tax rate
TAU_K_h_bar = 0.184;

// Steady state consumption tax rate
TAU_C_h_bar = 0.028;

// Steady state ricardian labour tax rate
TAU_L_h_bar = 0.223;

TK_h_bar  = 0.058;
TL_h_bar = 0.071;
TC_h_bar  = 0.019;

// 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.9609;
SIGMA_A_H        = 0.6241;

RHO_G            = 0.9691;
SIGMA_G_H        = 3.0582;

RHO_TK           = 0.9349;
SIGMA_TK_H       = 4.3819;

RHO_TC           = 0.9321;
SIGMA_TC_H       = 4.0445;

RHO_TL          = 0.9725;
SIGMA_TL_H      = 3.0029;

RHO_Z            = 0.9437;
SIGMA_Z_H        = 3.3809;


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

model;

%#R_K_h_ss   = (1/BETA + DELTA - 1) / (1- exp(TAU_K_h_bar));

%#K_h_ss    = ( ( ((1-H)*( exp(R_K_h_ss/ALPHA)*(1-GY_h_bar-BY_f_bar*(exp(1/BETA)-1))-DELTA) ))^(-GAMMA)* ((1-exp(TAU_L_h_bar))*(1-ALPHA)/(1+exp(TAU_C_h_bar)))*(exp(R_K_h_ss)/ALPHA)^((ALPHA+KAPPA)/(ALPHA-1)) ) ^ (1/(GAMMA+KAPPA)) ;                

%#L_h_ss     =  (exp(R_K_h_ss)/ALPHA)^(1/(1-ALPHA)) * exp(K_h_ss);     

%#Y_h_ss     = exp(K_h_ss)^ALPHA  * exp(L_h_ss)^(1-ALPHA);

%#G_h_ss     = GY_h_bar   *   exp(Y_h_ss);

%#B_h_ss     = BY_h_bar   *   exp(Y_h_ss);

%#I_h_ss     = DELTA * K_h_ss;

%#C_h_ss     = exp(Y_h_ss) - exp(I_h_ss) - G_h_bar - (1/BETA -1)*B_f_bar;

%#Z_h_ss     = exp(TK_h_bar) + exp(TL_h_bar) + exp(TC_h_bar) - G_h_bar - exp(1/BETA -1)*B_bar;

//Aggregation
exp(c_h)     = THETA * exp(c_h_r) + (1-THETA) * exp(c_h_o)  ;
exp(l_h)     = THETA * exp(l_h_r) + (1-THETA) * exp(l_h_o)  ;

// Production function
exp(y_h) = exp(u_a_h)*(exp(k_h(-1)))^ALPHA * exp(l_h)^(1-ALPHA);

//Law of motion
exp(k_h)   = (1-DELTA)*exp(k_h(-1)) + exp(i_h);

//   MPL and MPK
exp(w_h_r) =      THETA *(1-ALPHA)* exp(y_h)/exp(l_h)  ;
exp(w_h_o) =   (1-THETA)*(1-ALPHA)* exp(y_h)/exp(l_h)  ;
exp(r_k_h) =                ALPHA * exp(y_h)/(exp(k_h(-1))) ;


// intertemporal conditions
(exp(c_h_o) - H*exp(c_h_o(-1)) )^(-GAMMA)/(1+exp(tau_c_h)) = BETA * (exp(c_h_o(+1))- H*exp(c_h_o) )^(-GAMMA) / (1+exp(tau_c_h(+1))) *  exp(r_h(+1));
(exp(c_h_o) - H*exp(c_h_o(-1)) )^(-GAMMA)/(1+exp(tau_c_h)) = BETA * (exp(c_h_o(+1))- H*exp(c_h_o) )^(-GAMMA) / (1+exp(tau_c_h(+1))) *((1-exp(tau_k_h(+1))) * exp(r_k_h(+1)) + 1 - DELTA);

//intratemporal conditions
exp(w_h_o) = PHI*(1-exp(l_h_o))^ KAPPA /(1-exp(tau_l_h)) * (1+exp(tau_c_h))/ (exp(c_h_o)     - H*exp(c_h_o(-1)) )^(-GAMMA); 
exp(w_h_r) = PHI*(1-exp(l_h_r))^ KAPPA /(1-exp(tau_l_h)) * (1+exp(tau_c_h))/ (exp(c_h_r)     - H*exp(c_h_r(-1)) )^(-GAMMA); 

//GBC
exp(b_h) + exp(b_f) + exp(tk_h) + exp(tc_h) + exp(tl_h)  = exp(r_h(-1)) * exp(b_h(-1)+b_f(-1)) + exp(g_h) + exp(z_h);


// Optimizing agents BC
(1-exp(tau_l_h)) * exp(w_h_o) * exp(l_h_o) + (1-exp(tau_k_h)) * exp(r_k_h) * exp(k_h(-1))  +  exp(r_h) * exp( b_h(-1)) +   (1-THETA) *exp(z_h)  = 
 (1+exp(tau_c_h))*exp(c_h_o) +  exp(i_h) + exp(b_h);


//ROT constraint
(1+exp(tau_c_h)) * exp(c_h_r) = (1 - exp(tau_l_h)) * exp(l_h_r)* exp(w_h_r) + THETA * exp(z_h) ;

//Total capital tax revenue 
exp(tk_h)  = exp(tau_k_h) * exp(r_k_h)  *  exp(k_h(-1));

//Total  labor tax revenue 
exp(tl_h) = exp(tau_l_h) * exp(w_h_o) * exp(l_h_o);

// Total consumption tax revenue 
exp(tc_h) = exp(tau_c_h) * exp(c_h);

//Interest rate equation
exp(r_h)  = R_w + PSI2*(exp(b_f - B_f_bar) - 1) ;


%--------------------------------------------------------------------------
% Fiscal rules
%--------------------------------------------------------------------------

//Government spending policy rule
g_h = - EPSI_G*y_h - GAMMA_G*(b_h(-1)+b_f(-1)) + u_g_h  + ( EPSI_G*log(Y_h_bar) + GAMMA_G*log(B_h_bar) + log(G_h_bar) ) ;

//Capital tax policy rule
tau_k_h =  EPSI_K*y_h + GAMMA_K*(b_h(-1)+b_f(-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) ) ;

// Labour tax policy rule
tau_l_h =  EPSI_L*y_h + GAMMA_L*(b_h(-1)+b_f(-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) ) ;

//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)  ;

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

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_tc_h  = RHO_TC  * u_tc_h(-1)  +  e_u_tc_h;
u_tl_h  = RHO_TL  * u_tl_h(-1)  +  e_u_tl_h;
u_z_h   = RHO_Z   * u_z_h(-1)   +  e_u_z_h;

end;

%--------------------------------------------------------------------------
% Initial values
%--------------------------------------------------------------------------

initval;
y_h     = -0.046;
c_h     = -0.366;
c_h_o   = -0.009;
c_h_r   = -0.926;
l_h     = -1.109;
l_h_o   = -2.124;
l_h_r   = -0.615;
i_h     = -1.057;
k_h     = 2.806;
b_h     = -1.126;
w_h_o   = -0.030;
w_h_r   = -0.030;
r_h     = 0.010;
r_k_h   = -3.267;

u_a_h   = 0;
u_g_h   = 0;
u_tk_h  = 0;
u_tc_h  = 0;
u_tl_h  = 0;
u_z_h   = 0;

tk_h    = -2.8;
tc_h    = -3.9;
tl_h    = -2.6;

tau_k_h  = log(TAU_K_h_bar);
tau_l_h  = log(TAU_L_h_bar);
tau_c_h  = log(TAU_C_h_bar);
g_h      = log(G_h_bar);
z_h      = log(Z_h_bar);

end;

%--------------------------------------------------------------------------
% shocks
%--------------------------------------------------------------------------

shocks;
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;

%--------------------------------------------------------------------------
% Compute 
%--------------------------------------------------------------------------

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


stoch_simul(order=1,hp_filter=1600,irf=0);


%--------------------------------------------------------------------------
% Observables:
%--------------------------------------------------------------------------
 
varobs 
c_h g_h tl_h tk_h b_h z_h;
 
%--------------------------------------------------------------------------
% Priors: 28 parameters to estimate
%--------------------------------------------------------------------------


estimated_params;

GAMMA,  gamma_pdf, 1, 0.4;
KAPPA,  normal_pdf, 1, 0.3; 
H,      beta_pdf, 0.5, 0.1; 
THETA,  beta_pdf, 0.5, 0.1;
PSI2,   inv_gamma_pdf, 0.001, inf;

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;

stderr e_u_a_h,  inv_gamma_pdf, 0.1, 2;
stderr e_u_g_h,  inv_gamma_pdf, 0.1, 2;
stderr e_u_tk_h, inv_gamma_pdf, 0.1, 2;
stderr e_u_tl_h, inv_gamma_pdf, 0.1, 2;
stderr e_u_tc_h, inv_gamma_pdf, 0.1, 2;
stderr e_u_z_h,  inv_gamma_pdf, 0.1, 2;

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;
estimated_params_init(use_calibration);


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

estimation (datafile = soe_data_obs_2,mode_compute=6,mh_replic=0);

%estimation(datafile = soe_data_obs_2,mode_compute=0, mode_file = leeper_soe_est_mode,nobs=193,first_obs=1, mh_replic=2000, mh_nblocks=2, mh_drop=0.005, mh_jscale=0.80, mode_check);

 
%--------------------------------------------------------------------------
% Options 
%--------------------------------------------------------------------------
% mh_replic : the number of draws for MH algorithm, default = 2000
 
% mh_nblocks: the number of parallel chains for MH algorithm, default = 2
 
% mh_drop: fraction of intially generated parameter vectors to be dropped as
% a burn in before doing posterior simulations, default = 0.5
 
% mh_jscale: scale parameter of the jumping distribution's covariance matrix
% Default (0.2) is rarely satisfactory, tune for an acceptance ratio of 25%-33%. 
% Increase if acceptance ratio is too high,decrease if too low. 
% Consider parameter specific values for this scale parameter, this can be done in the [estimated params], page 49 block. 
 
% -------------------------------------------------------------------------
% Algorithms
%--------------------------------------------------------------------------
% If mh_nblocks = 1; convergence diagnostics of Geweke (1992, 1999) is computed
% If mh_nblocks > 1; convergence diagnostics of Brooks and Gelman (1998) is used instead
% If mh_replic > 2000; MCMC diagnostics are generated