% 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; 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) + (1-delta_h)*beta*pi_h(+1) + kappa_h * mc;

    //eqn 3 foreign inflation
    pi_f = delta_f * pi_f(-1) + (1-delta_f)*beta*(pi_f(+1)) + 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 = yw + 1/sigma *((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)) + 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) 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;

% Path and Parameters
path = 'C:\Dani\Bildung\Uni\Master\MasterThesis\DSGE\';
filename = 'Input_allHP10000_20111019.xls';
n_endo = 8;
n_exo  = 8;
n = n_endo + n_exo;

% Forecast with zero shocks of DSGE Model

% VAR - Matrix
/*
A = oo_.dr.ghx([7 8 9 20 19 21 23 24],[4 5 6 17 16 18 20 21]);
B = oo_.dr.ghu([7 8 9 20 19 21 23 24],:);
C = eye(8); C(1,1) = rho_epif; C(2,2)=rho_eq; C(3,3)=rho_ei; C(4,4)=rho_ez; C(5,5)=rho_eyw; C(6,6) = rho_epiw; C(7,7) = rho_eiw; C(8,8) = rho_ead;
save A B C;
*/
/*
Af = oo_.dr.ghx([7 8 9 20 19 21 23 24],[4 5 6 17 16 18 20 21]);
Bf = oo_.dr.ghu([7 8 9 20 19 21 23 24],:);
Cf = eye(8); C(1,1) = rho_epif; C(2,2)=rho_eq; C(3,3)=rho_ei; C(4,4)=rho_ez; C(5,5)=rho_eyw; C(6,6) = rho_epiw; C(7,7) = rho_eiw; C(8,8) = rho_ead;
save Af Bf Cf;
A = Af;
B = Bf;
C = Cf;
*/
% *** Forecast ***
% Load Initial Values, e.g. 1983 Q1
/*
[data_init, vnam] = xlsread(strcat(path,filename));
date = vnam(15:end-1,1);
data = [data_init(14:14+rows(date)-1,[26 12 13 15 18 16 20 21])];
*/
/*
sh = zeros(rows(date),n_exo);
vsh = zeros(rows(date),n_exo);
vsh(1,:) = inv(B) * data(1,:)';
sh(1,:)  = vsh(1,:);
for t = 2:rows(date)
    vsh(t,:) = (inv(B) * (data(t,:) - (A * data(t-1,:)')')')'; 
    sh(t,:)  = vsh(t,:) - (C * vsh(t-1,:)')';
end
save sh;
save vsh;

sh1 = zeros(112,n_endo); sh2 = zeros(112,n_endo); sh3 = zeros(112,n_endo); sh4 = zeros(112,n_endo); sh5 = zeros(112,n_endo); sh6 = zeros(112,n_endo); sh7 = zeros(112,n_endo);  sh8 = zeros(112,n_endo);
sh1(:,1) = sh(:,1); sh2(:,2) = sh(:,2); sh3(:,3) = sh(:,3); sh4(:,4) = sh(:,4); sh5(:,5) = sh(:,5); sh6(:,6) = sh(:,6); sh7(:,7) = sh(:,7); sh8(:,8) = sh(:,8);
save sh1; save sh2; save sh3; save sh4; save sh5; save sh6; save sh7; save sh8;
*/
/*
load sh; load sh1; load sh2; load sh3; load sh4; load sh5; load sh6; load sh7; load sh8;
load A B C;
*/
% *** Rebuild Y and others
/*
forc_y  = zeros(rows(date),n_endo);
vsh = zeros(112,n_endo);
vsh(1,:) = sh(1,:);
forc_y(1,:) = B * vsh(1,:)';
for t = 2:rows(date)
    vsh(t,:) = C * vsh(t-1,:)' + sh(t,:)';
    forc_y(t,:) = (A * forc_y(t-1,:)')' + (B * vsh(t,:)')';
end

forc_y1  = zeros(rows(date),n_endo);
vsh1 = zeros(112,n_endo);
vsh1(1,:) = sh1(1,:);
forc_y1(1,:) = B * vsh1(1,:)';
for t = 2:rows(date)
    vsh1(t,:) = C * vsh1(t-1,:)' + sh1(t,:)';
    forc_y1(t,:) = (A * forc_y1(t-1,:)')' + (B * vsh1(t,:)')';
end

forc_y2  = zeros(rows(date),n_endo);
vsh2 = zeros(112,n_endo);
vsh2(1,:) = sh2(1,:);
forc_y2(1,:) = B * vsh2(1,:)';
for t = 2:rows(date)
    vsh2(t,:) = C * vsh2(t-1,:)' + sh2(t,:)';
    forc_y2(t,:) = (A * forc_y2(t-1,:)')' + (B * vsh2(t,:)')';
end

forc_y3  = zeros(rows(date),n_endo);
vsh3 = zeros(112,n_endo);
vsh3(1,:) = sh3(1,:);
forc_y3(1,:) = B * vsh3(1,:)';
for t = 2:rows(date)
    vsh3(t,:) = C * vsh3(t-1,:)' + sh3(t,:)';
    forc_y3(t,:) = (A * forc_y3(t-1,:)')' + (B * vsh3(t,:)')';
end

forc_y4  = zeros(rows(date),n_endo);
vsh4 = zeros(112,n_endo);
vsh4(1,:) = sh4(1,:);
forc_y4(1,:) = B * vsh4(1,:)';
for t = 2:rows(date)
    vsh4(t,:) = C * vsh4(t-1,:)' + sh4(t,:)';
    forc_y4(t,:) = (A * forc_y4(t-1,:)')' + (B * vsh4(t,:)')';
end

forc_y5  = zeros(rows(date),n_endo);
vsh5 = zeros(112,n_endo);
vsh5(1,:) = sh5(1,:);
forc_y5(1,:) = B * vsh5(1,:)';
for t = 2:rows(date)
    vsh5(t,:) = C * vsh5(t-1,:)' + sh5(t,:)';
    forc_y5(t,:) = (A * forc_y5(t-1,:)')' + (B * vsh5(t,:)')';
end

forc_y6  = zeros(rows(date),n_endo);
vsh6 = zeros(112,n_endo);
vsh6(1,:) = sh6(1,:);
forc_y6(1,:) = B * vsh6(1,:)';
for t = 2:rows(date)
    vsh6(t,:) = C * vsh6(t-1,:)' + sh6(t,:)';
    forc_y6(t,:) = (A * forc_y6(t-1,:)')' + (B * vsh6(t,:)')';
end

forc_y7  = zeros(rows(date),n_endo);
vsh7 = zeros(112,n_endo);
vsh7(1,:) = sh7(1,:);
forc_y7(1,:) = B * vsh7(1,:)';
for t = 2:rows(date)
    vsh7(t,:) = C * vsh7(t-1,:)' + sh7(t,:)';
    forc_y7(t,:) = (A * forc_y7(t-1,:)')' + (B * vsh7(t,:)')';
end

forc_y8  = zeros(rows(date),n_endo);
vsh8 = zeros(112,n_endo);
vsh8(1,:) = sh8(1,:);
forc_y8(1,:) = B * vsh8(1,:)';
for t = 2:rows(date)
    vsh8(t,:) = C * vsh8(t-1,:)' + sh8(t,:)';
    forc_y8(t,:) = (A * forc_y8(t-1,:)')' + (B * vsh8(t,:)')';
end
*/
/*
Ag = oo_.dr.ghx([2],[4 5 6 17 16 18 20 21]);
Bg = oo_.dr.ghu([2],:);

forc_ygap  = zeros(rows(date),1);
forc_ygap(1,:) = Bg * vsh(1,:)';
for t = 2:rows(date)
    forc_ygap(t,:) = (Ag * forc_y(t-1,:)')' + (Bg * vsh(t,:)')';
end

forc_ygap1  = zeros(rows(date),1);
forc_ygap1(1,:) = Bg * vsh1(1,:)';
for t = 2:rows(date)
    forc_ygap1(t,:) = (Ag * forc_y1(t-1,:)')' + (Bg * vsh1(t,:)')';
end

forc_ygap2  = zeros(rows(date),1);
forc_ygap2(1,:) = Bg * vsh2(1,:)';
for t = 2:rows(date)
    forc_ygap2(t,:) = (Ag * forc_y2(t-1,:)')' + (Bg * vsh2(t,:)')';
end

forc_ygap3  = zeros(rows(date),1);
forc_ygap3(1,:) = Bg * vsh3(1,:)';
for t = 2:rows(date)
    forc_ygap3(t,:) = (Ag * forc_y3(t-1,:)')' + (Bg * vsh3(t,:)')';
end

forc_ygap4  = zeros(rows(date),1);
forc_ygap4(1,:) = Bg * vsh4(1,:)';
for t = 2:rows(date)
    forc_ygap4(t,:) = (Ag * forc_y4(t-1,:)')' + (Bg * vsh4(t,:)')';
end

forc_ygap5  = zeros(rows(date),1);
forc_ygap5(1,:) = Bg * vsh5(1,:)';
for t = 2:rows(date)
    forc_ygap5(t,:) = (Ag * forc_y5(t-1,:)')' + (Bg * vsh5(t,:)')';
end

forc_ygap6  = zeros(rows(date),1);
forc_ygap6(1,:) = Bg * vsh6(1,:)';
for t = 2:rows(date)
    forc_ygap6(t,:) = (Ag * forc_y6(t-1,:)')' + (Bg * vsh6(t,:)')';
end

forc_ygap7  = zeros(rows(date),1);
forc_ygap7(1,:) = Bg * vsh7(1,:)';
for t = 2:rows(date)
    forc_ygap7(t,:) = (Ag * forc_y7(t-1,:)')' + (Bg * vsh7(t,:)')';
end

forc_ygap8  = zeros(rows(date),1);
forc_ygap8(1,:) = Bg * vsh8(1,:)';
for t = 2:rows(date)
    forc_ygap8(t,:) = (Ag * forc_y8(t-1,:)')' + (Bg * vsh8(t,:)')';
end
*/
% save forc_yflex; save forc_yflex1; save forc_yflex2; save forc_yflex3; save forc_yflex4; save forc_yflex5; save forc_yflex6;
% save forc_y; save forc_y1; save forc_y2; save forc_y3; save forc_y4; save forc_y5; save forc_y6; save forc_y7; save forc_y8;
% load y; load y1; load y2; load y3; load y4; load y5; load y6; load y7; load y8;
% load forc_y; load forc_y1; load forc_y2; load forc_y3; load forc_y4; load forc_y5; load forc_y6; load forc_y7; load forc_y8;
% load forc_yflex forc_yflex1 forc_yflex2 forc_yflex3 forc_yflex4 forc_yflex5 forc_yflex6 forc_yflex7 forc_yflex8;
% [out1,out2]=xlswrite(strcat(path,'output\Out&PotOut.xls'),[forc_ygap forc_ygap1 forc_ygap2 forc_ygap3 forc_ygap4 forc_ygap5 forc_ygap6 forc_ygap7 forc_ygap8]);