% 3-equation New Keynesian model

% Daniel Geissmann, based on Nicolas Cuche-Curti
% University of St.Gallen, Master Thesis

%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------

close all;


%----------------------------------------------------------------
% 1. Defining variables and parameters
%----------------------------------------------------------------

var 		mc x y z c pi_h pi_f phi_f pi q s piw yw i iw ygap yflex xflex veps_pif veps_z veps_i veps_piw veps_iw veps_yw veps_ad veps_q;
varexo 		eps_pif eps_q eps_i eps_z eps_yw eps_piw eps_iw eps_ad;
parameters 	gamma phi sigma h delta_h beta kappa_h delta_f kappa_f theta_h theta_f eta rho_i phi_pi phi_a phi_s phi_yy A B phi_xz rho_z rho_piw rho_iw omega_s rho_y rho_yw rho_epif rho_epiw rho_eiw rho_eyw rho_ead rho_eq rho_ei rho_ez;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

beta = 0.9961;
gamma = 0.25;
sigma = 1;
phi = 0.75;
h = 0.85;
eta = 0.06;
theta_h = 0.97;
theta_f = 0.96;
rho_z = 0.73;
rho_i = 0.90;
phi_pi = 1.12;
phi_a = 0.38;
phi_s = 0.30;
delta_h = 0.78;
delta_f = 0.5;
phi_yy = -0.82;
phi_xz = 0.35;
rho_iw = 0.1;
rho_piw = 0.5;
rho_yw = 0.8;
rho_y = 0.8;

A = h*sigma*(phi*gamma*eta*(2-gamma)+1)/(sigma*(phi*gamma*eta*(2-gamma)+1)+(1-h)*((1-gamma)^2)*phi);
B = sigma*(1+phi)/(sigma*(phi*gamma*eta*(2-gamma)+1) + (1-h)*((1-gamma)^2)*phi);
kappa_h = (1-theta_h)*(1-beta*theta_h)/theta_h;
kappa_f = (1-theta_f)*(1-beta*theta_f)/theta_f;
omega_s = 1+gamma*(2-gamma)*(sigma*eta-1);

rho_epif = 0.7; % Share of Persistence of UIP shock
rho_eq   = 0.7; % Share of Persistence of Monetary policy shock
rho_ei   = 0.7; % Share of Persistence of world Demand shock
rho_ez   = 0.7; % Share of Persistence of domestic cost push shock
rho_eyw  = 0.7; % Share of Persistence of mark-up shock of imported inflation
rho_eiw  = 0.7; % Share of Persistence of world interest rate shock
rho_ead  = 0.7; % Share of Persistence of domestic demand shock
rho_epiw = 0.7;

% Variance of the Shocks
sigma_pih = 0.01;
sigma_pif = 0.01;
sigma_q   = 0.01;
sigma_i   = 0.01;
sigma_z   = 0.01;
sigma_yw  = 0.01;
sigma_piw = 0.01;
sigma_iw  = 0.01;
sigma_ad  = 0.01;
%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model(linear);
    //eqn 1 marginal costs
	mc = gamma * x + phi * y - (1+phi) * z + sigma * 1/(1-h) * (c - h*c(-1));
    
    //eqn 2 Domestic inflation
    pi_h = delta_h * pi_h(-1) + beta*(pi_h(+1) - delta_h * pi_h) + kappa_h * mc;

    //eqn 3 foreign inflation
    pi_f = delta_f * pi_f(-1) + beta*(pi_f(+1) - delta_f * pi_f) + kappa_f * phi_f + veps_pif;

    //eqn 4 CPI inflation
    pi = pi_h + gamma * (x-x(-1));

    //eqn 5 Real exchange rate
    q = (1-gamma) * x + phi_f;

    //eqn 6 LOP gap
    phi_f = phi_f(-1) + s - s(-1) + piw  - pi_f;

    //eqn 7 Terms of trade
    x = x(-1)+ pi_f - pi_h;

    //eqn 8 International Risk sharing
    c = h*c(-1) + yw - h*yw(-1) + 1/sigma *(1-h) *((1-gamma)*x + phi_f);
    
    // Aggregate Demand
    //y = (1-rho_y) * y(+1) + rho_y * y(-1) + (1-omega_s/(1-gamma)) * (yw(+1) - yw) + gamma*eta / (1-gamma) * (phi_f(+1) - phi_f) - omega_s/(sigma*(1-gamma)) * (i - pi(+1)) + veps_ad;

    //eqn 9 UIP
    i = pi(+1) + iw - piw(+1) + q(+1) - q + veps_q;
    
    //eqn 10 Aggregate Demand
    c = 1/(1-gamma) *(y - gamma*eta*(2-gamma)*x - gamma * eta* phi_f - gamma *yw);

    //eqn 11 MPR
    i = rho_i*i(-1) + (1-rho_i)*(phi_pi*pi + phi_a* ygap + phi_s*(s-s(-1))) + veps_i;
    
    //eqn 12 output gap
    ygap = y - yflex;
    
    //eqn 13 potential output
    yflex = (1+phi)/phi*z - xflex/phi+phi_yy*yw;
   
    //eqn 14 potential output
    xflex = A*xflex(-1) + B*(z-h*z(-1)) - phi_xz*yw;
    
    //eqn 15 Technology process
    z = rho_z*z(-1) + veps_z;

    //eqn 16 foreign processes
    yw = rho_yw * yw(-1) + veps_yw;
    piw = rho_piw * piw(-1) + veps_piw;
    iw = rho_iw * iw(-1) + veps_iw;

    // shock equations:
    veps_piw = rho_epiw * veps_piw(-1) + eps_piw;
    veps_q   = rho_eq   * veps_q(-1)   + eps_q;
    veps_i   = rho_ei   * veps_i(-1)   + eps_i;
    veps_yw  = rho_eyw  * veps_yw(-1)  + eps_yw;
    veps_z   = rho_ez   * veps_z(-1)   + eps_z;
    veps_pif = rho_epif * veps_pif(-1) + eps_pif;
    veps_iw  = rho_eiw  * veps_iw(-1)  + eps_iw;
    veps_ad  = rho_ead  * veps_ad(-1)  + eps_ad;
end;


%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

initval; mc=0; x=0; y=0; z=0; c=0; pi_h=0; pi_f=0; phi_f=0; pi=0; q=0; s=0; piw=0; yw=0; iw=0; ygap=0; yflex=0; xflex=0; veps_q=0; veps_pif=0; veps_piw=0; veps_yw=0; veps_iw=0; veps_i=0; veps_z=0; veps_ad=0; end;

shocks;
var eps_pif = sigma_pif^2;
var eps_q   = sigma_q^2;
var eps_i   = sigma_i^2;
var eps_z   = sigma_z^2;
var eps_yw  = sigma_yw^2;
var eps_piw = sigma_piw^2;
var eps_iw  = sigma_iw^2;
var eps_ad  = sigma_ad^2;
end;

check;
//steady;
stoch_simul(irf=32,order=1,periods=10000,drop=200,nograph);