%
% By Nikolas Kontogiannis
% Department of Economics
% University of Leicester
%
% This is a two country New Keynesian DSGE model with labour market frictions. Below I solve the model to display the IRF
% 

%variables
var 
u_a %unemployment in country A
u_b %unemployment in country B

pi_ad %domestic inflation country A
pi_bd %domestic inflation country B 

s_ba %terms of trade

cpi_a %cpi inflation in country A
cpi_b %cpi inflation in country B

r_a %real interest rate in country A
r_b %real interest rate in country B

q_ab %short run nominal interest rate common for both countries

pi_wd %union wide inflation
u_w %union wide unemployment

z_a %productivity in country A
z_b; %productivity in country B

varexo 
e_a %stochastic technology shock in country A
e_b; %stochastic technology shcok in country B

%parameters

parameters 
%Preferences parameters
zeta %country size, openness 
beta %discount factor
sigma %coefficient of relative risk aversion
gamma %elasticity of substitution among differentiated final goods
varphi %Frisch elasticity of labour
UCC %marginal utility of consumption times consumption. useful for the calculation of the optimal weight of the objective function
UCN %marginal utility of consumption times labour. useful for the calculation of the optimal weight of the objective function

%Final good firms parameters
omega_a %nominal price rigidity (probability at period t), country A
omega_b %nominal price rigidity (probability at period t), country B
delta_ap %elasticity of inflation with respect to the real marginal cost, country A
delta_bp %elasticity of inflation with respect to the real marginal cost, country B

%Labour market parameters (intermediate good firms)
q_theta %probability of firms to fill a vacancy
delta %exogenous separation rate
kappa %elasticity of matching function with respect to vacancies
xi_a %bargaining power of firms, country A
xi_b %bargaining power of firms, country B
mu_a %degree of real wage rigidity, country A
mu_b %degree of real wage rigidity, country B

%Labour market parameters (after linearization)
alpha_0 %parameter3 for optimal weight in unemployment
alpha_a1 alpha_b1
alpha_a2 alpha_b2
alpha_a3 alpha_b3
gamma_a0 gamma_b0
gamma_a1 gamma_b1
gamma_a2 gamma_b2
delta_a1 %efficient steady state parameter1
delta_2 %parameter1 for optimal weight in unemployment(eliminated for 2nd order) and efficient steady state parameter2
delta_3 %parameter2 for optimal weight in unemployment
delta_u

%Steady-state levels. Because I linearize by hand, in matlab code, I treat to these as parameters
Z %steady state of productivity (level)
S %steady state of terms of trade (level)
phi %steady state of real marginal cost (level)
u %steady state of unemployment (level)
C_a %steady state value of consumption (level), country A
C_b %steady state value of consumption (level), country B
W %steady state of real wage
N %steady state value of employment (level)
Y %steady state value of output (level)
v %steady state value of vacancies (level)
theta %steady state of labour market tightness (level)
psi %cost of posting vacancies per vacancy
d_a %disutility of labour constant coefficient, country A
d_b %disutility of labour constant coefficient, country B

%Interest rate rule parameters (weights for lagged values are given as well
%to satisfy the equilibrium determinacy
omega_r %interest rate smoothing parameter in Taylor
omega_p %weight in lagged inflation in Taylor
omega_s %weight in the lagged terms of trade in Taylor
omega_u %weight in lagged unemployment in Taylor
omega_p1 %weight in current inflation
omega_s1 %weight in current terms of trade
omega_u1 %weight in current unemployment

%Parameters of the constraints
%parameters in dynamic IS in country A
phi_a1 phi_a2 phi_a3 phi_a4 

%parameters in dynamic IS in country B
phi_b1 phi_b2 phi_b3 phi_b4

%parameters in NKPC in country A
rho_a0 rho_a1 rho_a2 rho_a3 rho_a4

%parameters in NKPC in country B
rho_b0 rho_b1 rho_b2 rho_b3 rho_b4

%Autocorrelation parameters of AR(1) productivity
rho_a rho_b

%weights on the objective function
omega_pia omega_pib omega_ua omega_ub omega_sss;

% parameters values
beta=0.99;
zeta=0.5;
sigma=1;
varphi=0;
kappa=0.5;
xi_a=0.5;
omega_a=0.75;
q_theta=0.97;
delta=0.10;
Z=1;
S=1;
phi=0.83;
u=0.10;
mu_a=0.5;
rho_a=0.95;
rho_b=0.95;
gamma=6;
delta_p=0.0858;

%parameter values different for country B
xi_b=0.5;
omega_b=0.75;
mu_b=0.5;

%interest rate rule parameters (satisfy determinacy of equilibrium)
omega_r = 0.85;
omega_s=1.5;
omega_p=1.5;
omega_u=1.5;
omega_p1=0.5;
omega_u1=0.5;
omega_s1=0.5;

% steady state values (levels)
N=(1-u)/(1-delta);
Y=Z*N;
v=(delta*N)/q_theta;
theta=v/u;
psi=(0.01*Y)/v;
C_a=Y*S^(zeta-1) - psi*v;
W=phi +(1-delta)*beta*(psi/q_theta) - psi/q_theta;
d_a=(W-(1-xi_a)*(phi + (1-delta)*beta*psi*theta))/( xi_a*(N^varphi)*C_a^(1/sigma));

% steady state values (levels) different for country B 
C_b=Y*S^(zeta) - psi*v;
d_b=(W-(1-xi_b)*(phi + (1-delta)*beta*psi*theta))/( xi_b*(N^varphi)*C_b^(1/sigma));

%parameters after linearization, useful for the dynamic IS, country A
alpha_0=(1-delta)*(N/u);
alpha_a2=(psi*v)/(alpha_0*delta*kappa*C_a);
alpha_a1=Y/((S^(1-zeta))*C_a);
alpha_a3=alpha_a1/alpha_0;
delta_a1=1/(S^(1-zeta)) - (psi*v)/delta*kappa*N - d_a*(N^varphi)*(C_a^(1/sigma));
delta_2=(1-delta)*((psi*v)/kappa)*(1/(delta*N) - ((1-kappa)/u));
delta_3=(N/(sigma*C_a))*S^(2*(1-zeta));
delta_u=((1-delta)-alpha_0*delta);
gamma_a0=(1-mu_a)*(1-xi_a);
gamma_a1=(1-gamma_a0)*phi;
gamma_a2=(d_a)*(N^varphi)*C_a^(1/sigma);
delta_ap=((1-omega_a*beta)*(1-omega_a)/(omega_a));


%parameters after linearization, useful for the dynamic IS, country B
alpha_b1=Y/((S^(-zeta))*C_b);
alpha_b2=(psi*v)/(alpha_0*delta*kappa*C_b);
alpha_b3=alpha_b1/alpha_0;
delta_b1=1/(S^(-zeta)) - (psi*v)/delta*kappa*N - d_b*(N^varphi)*(C_b^(1/sigma));
gamma_b0=(1-mu_b)*(1-xi_b);
gamma_b1=(1-gamma_b0)*phi;
gamma_b2=(d_b)*(N^varphi)*C_b^(1/sigma);
delta_bp=((1-omega_b*beta)*(1-omega_b)/(omega_b));


% rho's parameters for NKPC, country A
rho_a0=(delta_u/(alpha_0*delta*kappa*gamma_a1))*((1-kappa)*psi*theta^(1-kappa)-(((1-mu_a)*xi_a*psi*v*gamma_a2)/(sigma*C_a))) - ((1-mu_a)*xi_a*psi*v*gamma_a2)/(sigma*C_a);
rho_a1=(1/(alpha_0*delta*kappa*gamma_a1))*(delta_u*(1-delta)*beta*theta*psi*(gamma_a0-(1-kappa))-(1-kappa)*psi*theta^(1-kappa)-(1-mu_a)*xi_a*gamma_a2*((((alpha_a1)*delta*kappa)/sigma)+varphi*delta*kappa - ((psi*v)/(sigma*C_a))));
rho_a2=(((1-delta)*beta*theta*psi)/(alpha_0*delta*kappa*gamma_a1))*(gamma_a0-(1-kappa)*theta^(-kappa));
rho_a3=1 - ((1-mu_a)*xi_a*((gamma_a2*alpha_a1)/(sigma*gamma_a1))); 
rho_a4=((((1-delta)*beta*theta)*psi)/gamma_a1)*(theta^(-kappa)-gamma_a0);

% rho's parameters for NKPC country B
rho_b0=(delta_u/(alpha_0*delta*kappa*gamma_b1))*((1-kappa)*psi*theta^(1-kappa)-(((1-mu_b)*xi_b*psi*v*gamma_b2)/(sigma*C_b))) - ((1-mu_b)*xi_b*psi*v*gamma_b2)/(sigma*C_b);
rho_b1=(1/(alpha_0*delta*kappa*gamma_b1))*(delta_u*(1-delta)*beta*theta*psi*(gamma_b0-(1-kappa))-(1-kappa)*psi*theta^(1-kappa)-(1-mu_b)*xi_b*gamma_b2*((((alpha_b1)*delta*kappa)/sigma)+varphi*delta*kappa - ((psi*v)/(sigma*C_b))));
rho_b2=(((1-delta)*beta*theta*psi)/(alpha_0*delta*kappa*gamma_b1))*(gamma_b0-(1-kappa)*theta^(-kappa));
rho_b3=(((1-mu_b)*xi_b*((gamma_b2*alpha_b1)/(sigma*gamma_b1))) + ((zeta-1)/zeta));
rho_b4=((((1-delta)*beta*theta)*psi)/gamma_b1)*(theta^(-kappa)-gamma_b0);

%eta parameters, useful for IS 
eta_a1=((psi*v)/C_a)*(((1-delta)/(alpha_0*delta*kappa))*(1-((delta*N)/u))+1);
eta_a2=(1/alpha_0*C_a)*(psi*v/(delta*kappa) - Y/(S^(1-zeta)));
eta_b1=((psi*v)/C_b)*(((1-delta)/(alpha_0*delta*kappa))*(1-((delta*N)/u))+1);
eta_b2=(1/alpha_0*C_b)*(psi*v/(delta*kappa) - Y/(S^(-zeta)));


% phi's parameters, IS, country A
phi_a1=eta_a1/(eta_a1 + eta_a2);
phi_a2=eta_a2/(eta_a1 + eta_a2);
phi_a3=sigma/(eta_a1 + eta_a2);
phi_a4=alpha_a1/(eta_a1 + eta_a2);

% phi's parameters, IS, country B
phi_b1=eta_b1/(eta_b1 + eta_b2);
phi_b2=eta_b2/(eta_b1 + eta_b2);
phi_b3=sigma/(eta_b1 + eta_b2);
phi_b4=alpha_b1/(eta_b1 + eta_b2);

UCC = C_a*(C_a)^(-1/sigma);
UCN = N*(C_a)^(-1/sigma);

% weights on the policy objectives
omega_pia=UCC*(zeta*gamma/(2*delta_ap));
omega_piab=UCC*((1-zeta)*gamma/(2*delta_bp));
omega_ua=UCN*((zeta*delta_3)/(2*(alpha_0)^2));
omega_ub=UCN*(((1-zeta)*delta_3)/(2*(alpha_0)^2));
omega_sss=UCC*zeta*(1-zeta)*((1+sigma)/(2*sigma));


model (linear);
%dynamic IS expressed in terms of unemployment
u_a = phi_a1*u_a(-1) + phi_a2*u_a(+1) - phi_a3*r_a - phi_a4*(1-zeta)*s_ba(+1) + phi_a4*(1-zeta)*s_ba  - phi_a4*(1-rho_a)*z_a;
u_b = phi_b1*u_b(-1) + phi_b2*u_b(+1) - phi_b3*r_b + phi_b4*zeta*s_ba(+1) - phi_b4*zeta*s_ba  - phi_b4*(1-rho_b)*z_b;

%NKPC with unemployment. When the mu_a=mu_b then the rho coefficients are identical
pi_ad = beta*pi_ad(+1) + delta_ap*rho_a0*u_a(-1) + delta_ap*rho_a1*u_a - delta_ap*rho_a2*u_a(+1) + delta_ap*rho_a3*(1-zeta)*s_ba + delta_ap*rho_a4*r_a - delta_ap*rho_a3*z_a; 
pi_bd = beta*pi_bd(+1) + delta_bp*rho_b0*u_b(-1) + delta_bp*rho_b1*u_b - delta_bp*rho_b2*u_b(+1) + delta_bp*zeta*rho_b3*s_ba + delta_bp*rho_b4*r_b + delta_bp*rho_b3*z_b; 

%Link terms of trade with domestic unemployment
s_ba = s_ba(-1) + pi_bd - pi_ad;

%union wide inflation
pi_wd = zeta*pi_ad + (1-zeta)*pi_bd;

%union wide unemployment
u_w = zeta*u_a + (1-zeta)*u_b;

%Fisher equations for country A and B
r_a = q_ab - cpi_a(+1);
r_b = q_ab - cpi_b(+1);

%Taylor rule with domestic inflation and unemployment
q_ab = omega_r*q_ab(-1) + (1-omega_r)*(omega_p*(zeta*pi_ad(-1) + (1-zeta)*pi_bd(-1)) + omega_s*s_ba(-1) + omega_u*(zeta*u_a(-1) + (1-zeta)*u_b(-1))) + omega_p1*(zeta*pi_ad + (1-zeta)*pi_bd) + omega_s1*s_ba + omega_u1*(zeta*u_a + (1-zeta)*u_b);

%link of cpi with domestic inflation and terms of trade
cpi_a = pi_ad + (1-zeta)*(s_ba - s_ba(-1));
cpi_b = pi_bd + (zeta-1)*(s_ba - s_ba(-1));

%AR(1) technology in countries A and B
z_a = rho_a*z_a(-1) + e_a;
z_b = rho_b*z_b(-1) + e_b; 
end;


initval; 
u_a=0;
u_b=0;
pi_ad=0;
pi_bd=0;
pi_wd=0;
u_w=0;
cpi_a=0;
cpi_b=0;
s_ba=0;
r_a=0;
r_b=0;
q_ab=0;
z_a=0;
z_b=0;
end;

resid(1); 

shocks;
var e_a;
stderr 0.624;
var e_b;
stderr 0;
var e_a, e_b= 0*0.624*0.624;


end;

stoch_simul(order=2, periods=2100, irf=40);