%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BKK(94) economy with single international bond
% Kyriacos LAMBRIAS
% This version: 17 May 2016

% KPR preferences
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@#define countries = [ "1", "2" ]
% Number of lags for which news shocks as in play
@#define q = 0
 q= @{q}

@#for co in countries
var c@{co} n@{co} k@{co} y@{co} i@{co} z@{co} qa@{co} qb@{co} a@{co} b@{co} lamda@{co} 
G@{co} B@{co} r@{co} w@{co}  ; 

varexo u@{co} ;
@#endfor
var p Q P2 c_rel; %nx1_y1_nipa y_rel
 
@#for co in countries
parameters alpha_@{co} beta_@{co} delta_@{co} gamma_@{co} theta_@{co} omega_@{co} is_@{co}  
n_ss_@{co} qa@{co}_ss qb@{co}_ss A_1@{co} A_2@{co} sigma_@{co} su_@{co} phi_@{co} B_ss_@{co} tau_@{co}; %PSI_@{co} gammaJR_@{co} ni_@{co}
@#endfor
parameters corru y1_k1 c1_y1 i1_y1;

@#for co in countries
alpha_@{co} = 0.36 ;
beta_@{co} = 0.99;
delta_@{co} = 0.025;
gamma_@{co} = -1;
theta_@{co} = 1.5;
phi_@{co} = 0.00001;
@#endfor
B_ss_1 = 0 ;
B_ss_2 = 0 ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           DETERMINATION OF THE STEADY-STATE                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Note that I do not loop over countries. Assuming symmetry, all params are the same

is_1 = 0.15 ; % import share to GDP
r1_ss = 1/beta_1 - 1 + delta_1 ;
n_ss_1 = 0.3072 ; 

omega_1 = 1/((is_1/(1-is_1))^(1/theta_1)+1);
G1_y1 = (omega_1*(1-is_1)^((theta_1-1)/theta_1) + (1-omega_1)*is_1^((theta_1-1)/theta_1))^(theta_1/(theta_1-1));
qa1_ss = omega_1*(G1_y1/(1-is_1))^(1/theta_1); 
qb1_ss = qa1_ss ;
y1_k1 = r1_ss/(alpha_1*qa1_ss) ;
i1_y1 = delta_1/y1_k1 ; 
c1_y1 = G1_y1-i1_y1;
w1_y1 = (1-alpha_1)*qa1_ss/n_ss_1 ;
% Getting TAU right:
tau_1 = 1/( w1_y1*(1-n_ss_1)/c1_y1 + 1) ;
tau_2 = tau_1 ;

y1_ss = (y1_k1^(-alpha_1/(1-alpha_1)))*n_ss_1 ;
c1_ss = c1_y1*y1_ss ;
i1_ss = i1_y1*y1_ss ;
k1_ss = (1/y1_k1)*y1_ss ;
w1_ss = w1_y1*y1_ss ;
lamda1_ss = tau_1*c1_ss^(gamma_1*tau_1-1)*(1-n_ss_1)^(gamma_1*(1-tau_1));
G1_ss = G1_y1*y1_ss;
check2 = lamda1_ss*w1_ss - (1-tau_1)*c1_ss^(gamma_1*tau_1)*(1-n_ss_1)^(gamma_1*(1-tau_1)-1);

% Symmetry in ctry 2
omega_2 = omega_1 ;
is_2 = is_1;
n_ss_2 = n_ss_1 ;
qa2_ss = qa1_ss ;
qb2_ss = qa1_ss ;
P2_ss = tau_1*c1_ss^(gamma_1*tau_1-1)*(1-n_ss_1)^(gamma_1*(1-tau_1))/lamda1_ss; 
P2_ss1 = ((1/is_2)^(1/theta_2)*(1-omega_2)*( (1-omega_2)*is_2^((theta_2-1)/theta_2) + omega_2*(1-is_2)^((theta_2-1)/theta_2) )^(1/(theta_2-1)))^(-1)*qa1_ss; 

// Elements of the VAR(1) matrix of technology as in BKK (1994)
@#for co in countries
sigma_@{co} = 0.00852; 
su_@{co} = 0.00852 ;
@#endfor
A_11 = 0.906;  
A_12 = 0.088;
A_22 = A_11 ;
A_21 = A_12 ;
corru = 0.258;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                         MODEL                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

model; 
% Households

% FOC consumption
exp(lamda1) = tau_1*exp(c1)^(gamma_1*tau_1-1)*(1-exp(n1))^(gamma_1*(1-tau_1));
exp(lamda2)*exp(P2) = tau_2*exp(c2)^(gamma_2*tau_2-1)*(1-exp(n2))^(gamma_2*(1-tau_2));

% FOC leisure
@#for co in countries
exp(lamda@{co})*exp(w@{co}) = (1-tau_@{co})*exp(c@{co})^(gamma_@{co}*tau_@{co})*(1-exp(n@{co}))^(gamma_@{co}*(1-tau_@{co}) - 1);
@#endfor

% EE - Capital
exp(lamda1) = beta_1*exp(lamda1(+1))*(exp(r1(+1)) + 1 - delta_1);
exp(lamda2)*exp(P2) = beta_2*exp(lamda2(+1))*(exp(r2(+1)) + (1 - delta_2)*exp(P2(+1)));

% EE - Bond
@#for co in countries
exp(lamda@{co})*exp(Q) = beta_@{co}*exp(lamda@{co}(+1))*(1-phi_@{co}*(B@{co} - B_ss_@{co}));
@#endfor

% Final-good producer
exp(a1) = (omega_1/exp(qa1))^theta_1*exp(G1) ;
exp(b1) = ((1-omega_1)/exp(qb1))^theta_1*exp(G1) ; 
exp(a2) = ((1-omega_2)/(exp(qa2)/exp(P2)))^theta_2*exp(G2) ;
exp(b2) = (omega_2/(exp(qb2)/exp(P2)))^theta_2*exp(G2) ;

% Intermediate-good producer
exp(r1) = alpha_1*exp(qa1)*exp(y1)/exp(k1(-1)) ; 
exp(w1) = (1-alpha_1)*exp(qa1)*exp(y1)/exp(n1) ; 
exp(r2) = alpha_2*exp(qb2)*exp(y2)/exp(k2(-1)) ; 
exp(w2) = (1-alpha_2)*exp(qb2)*exp(y2)/exp(n2) ; 

% Resource constraints

exp(G1) = (omega_1*exp(a1)^((theta_1-1)/theta_1) + (1-omega_1)*exp(b1)^((theta_1-1)/theta_1))^(theta_1/(theta_1-1));
exp(G2) = (omega_2*exp(b2)^((theta_2-1)/theta_2) + (1-omega_2)*exp(a2)^((theta_2-1)/theta_2))^(theta_2/(theta_2-1));
//exp(G1) = exp(c1) + exp(i1);
exp(G2) = exp(c2) + exp(i2);
exp(y1) = exp(a1) + exp(a2); 
exp(y2) = exp(b1) + exp(b2);

@#for co in countries
% CD PF
exp(y@{co}) = exp(z@{co})*exp(k@{co}(-1))^alpha_@{co}*exp(n@{co})^(1-alpha_@{co}) ;
% Law of motion of capital 
exp(k@{co}) = (1-delta_@{co})*exp(k@{co}(-1)) + exp(i@{co}); 
@#endfor

% Budget constraints
exp(c1) + exp(i1) + exp(Q)*B1 + phi_1*(B1(-1) - B_ss_1)^2/2 = exp(w1)*exp(n1) + exp(r1)*exp(k1(-1)) + B1(-1) ;
exp(P2)*(exp(c2) + exp(i2)) + exp(Q)*B2 + phi_2*(B2(-1) - B_ss_2)^2/2 = exp(w2)*exp(n2) + exp(r2)*exp(k2(-1)) + B2(-1) ;

% Equilibium in the international bond market
B1(-1) + B2(-1) = 0 ;

% Prices
exp(p) = exp(qb1)/exp(qa1) ;
exp(qa1) = exp(qa2) ;
exp(qb1) = exp(qb2) ;

// Stochastic processes 
z1 = A_11*z1(-1) + A_12*z2(-1) + su_1*u1(-@{q}) ;
z2 = A_21*z1(-1) + A_22*z2(-1) + su_2*u2(-@{q}) ;

% Definitions
exp(c_rel) = exp(c1)/exp(c2) ; // relative consumptions
//exp(n_rel) = exp(n1)/exp(n2) ;
//exp(y_rel) = exp(y1)/exp(y2) ;
//nx1_y1_nipa = (qa1_ss*exp(a2) - qb1_ss*exp(b1))/exp(y1);

end;

initval;
lamda1 = log(lamda1_ss) ;
p = log(1);
c1= log(c1_ss);
n1= log(n_ss_1);
k1= log(k1_ss);
y1= log(y1_ss);
G1 = log(G1_ss);
i1= log(i1_ss);
z1= log(1);
qa1= log(qa1_ss);
qb1= log(qa1_ss);
a1= log((1-is_1)*y1_ss);
b1= log(is_1*y1_ss);
B1 = B_ss_1;
r1 = log(r1_ss);
w1 = log(w1_ss);
Q = log(beta_1); 

% country 2
lamda2 = log(lamda1_ss) ;
c2= log(c1_ss);
n2= log(n_ss_1);
k2= log(k1_ss);
y2=  log(y1_ss);
G2 = log(G1_ss);
i2= log(i1_ss);
z2= log(1);
qa2= log(qa1_ss);
qb2= log(qb1_ss);
a2= log(is_1*y1_ss);
b2= log((1-is_1)*y1_ss);
P2= log(P2_ss) ;

B2 = B_ss_2 ;
r2 = log(r1_ss);
w2 = log(w1_ss);

end;

steady;
resid(1);
check ;

//write_latex_dynamic_model;

shocks ;
// That would adjust the var-cov matrix even when the sd of the shocks are not symmetric. be careful though...
@#for co in countries
var u@{co} = (su_@{co}^2)/(sigma_@{co}*su_@{co}) ;
@#endfor
var u1, u2 = (corru) ;
end ;

stoch_simul(order=1, hp_filter = 1600, ar=0, nograph, irf=10) P2 c_rel c1 c2 y1 y2 i1 i2 n1 n2 Q B1 B2 ;