%% General New Keynesian DSGE model
%
% Last edited by Ben on 12 August 2016

% Parameters as chosen in Section 19.4 in FoMM3

clear
clf
clc
close all

% Structural parameters
delta_a       = 0.10;                 % Annual depreciation rate for capital
PAR.delta     = 1-(1-delta_a)^(1/4);  % Quarterly depreciation rate for capital
rho_a         = 0.04;                 % Annual rate of time preference
PAR.rho       = (1+rho_a)^(1/4)-1;    % Quarterly depreciation rate for capital
PAR.alpha     = 1/3;                  % Efficiency parameter capital
PAR.phi       = 0.10;                 % Steady-state employment share of overhead workers
PAR.Omega_0   = 0.00; 
PAR.theta     = 0.00;
PAR.epsilon_c = 0.00;
PAR.epsilon_l = 0.00;
PAR.epsilon_m = 0.00;
PAR.zeta      = 0.75;
PAR.xi_Z      = 0.95;


% Exogenous variables
EXO.M_0       = 1.00;                 % Nominal money supply
EXO.G         = 0.20;                 % Government consumption
EXO.time      = linspace(0,200,2001); % Post-shock time (for plotting) 
EXO.Z0        = 1;

% Target variables
Y0            = 1.00;                     % Baseline per capita output per quarter 
P0            = 1.00;
L0            = 1/3;
PAR.Lbar      = PAR.phi * L0;
R0            = PAR.rho;
MARKUP0       = 1.30;
PAR.theta     = MARKUP0 / (MARKUP0 - 1);
mc0           = 1 / MARKUP0;
KYratio0      = PAR.alpha * mc0 / (PAR.rho + PAR.delta);
K0            = KYratio0 * Y0;
I0            = PAR.delta * K0;
C0            = Y0 - EXO.G - I0;
w0            = (1-PAR.alpha) * mc0 * Y0 / ((1-PAR.phi) * L0);
rk0           = PAR.alpha * mc0 * Y0 / K0 ;
eta_L         = w0 *(1 - L0) / C0;
eta_M         = R0 * EXO.M_0 / ((1+R0) * P0 * C0) ;
PAR.Omega_0   = Y0 * (K0)^(-PAR.alpha) * ((1-PAR.phi)*L0)^(PAR.alpha-1);
PAR.epsilon_c = 1 / (1 + eta_L + eta_M) ;
PAR.epsilon_l = eta_L / (1 + eta_L + eta_M) ;
PAR.epsilon_m = eta_M / (1 + eta_L + eta_M) ;
XI_D0         = ((1+PAR.rho) / (1+PAR.rho - PAR.zeta)) *  P0^(PAR.theta-1) * Y0 / C0;
XI_N0         = ((1+PAR.rho) / (1+PAR.rho - PAR.zeta)) *  P0^PAR.theta * Y0 * mc0 / C0;



tlin      = linspace(0,1.0,1001);

% Compute the steady-state equilibrium

x = zeros(8,1);
[x(:,1),fval,exitflag,output] = fsolve(@Function19_01_SS,[Y0,K0,L0,C0,w0,R0,mc0,P0]'...
    ,optimset('MaxFunEvals',1e12,'MaxIter',1000,'TolX',1.0e-8,'TolFun',1.0e-8),PAR,EXO);
                                     
fval;
output;
x;


VAR.Y0     = x(1,1);
VAR.K0     = x(2,1);
VAR.L0     = x(3,1);
VAR.C0     = x(4,1);
VAR.w0     = x(5,1);
VAR.R0     = x(6,1);
VAR.mc0    = x(7,1);
VAR.P0     = x(8,1);
VAR.rk0    = PAR.alpha * VAR.mc0 * VAR.Y0 / VAR.K0 ;
VAR.RP0    = (1 - VAR.mc0) * VAR.Y0 - VAR.w0 * PAR.Lbar ;  


shock = 3;
    
% Impulse Response Functions

Nperiods = 200;
PStime = linspace(0,Nperiods,Nperiods+1);

switch shock ;
    
    case 1    
%% Productivity shock

% ************* *
% Sticky prices *
% ************* *

zetaval = 0.75;
dynare MBC_Zshock noclearall

% For jumping variables:
% first observation is the initial steady state value:
Yinitial = oo_.endo_simul(1,1);
Linitial = oo_.endo_simul(4,1);
Zinitial = oo_.endo_simul(15,1);

% second observation is the impact value:
Ypath      = oo_.endo_simul(1,2:Nperiods+2);
Yapath     = oo_.endo_simul(2,2:Nperiods+2);
Lpath      = oo_.endo_simul(4,2:Nperiods+2);
Cpath      = oo_.endo_simul(5,2:Nperiods+2);
wpath      = oo_.endo_simul(6,2:Nperiods+2);
rkpath     = oo_.endo_simul(7,2:Nperiods+2);
Rpath      = oo_.endo_simul(8,2:Nperiods+2);
mcpath     = oo_.endo_simul(9,2:Nperiods+2);
Ppath      = oo_.endo_simul(10,2:Nperiods+2);
Papath     = oo_.endo_simul(11,2:Nperiods+2);
Pnpath     = oo_.endo_simul(12,2:Nperiods+2);
XI_Npath   = oo_.endo_simul(13,2:Nperiods+2);
XI_Dpath   = oo_.endo_simul(14,2:Nperiods+2);
Ztildepath = oo_.endo_simul(15,2:Nperiods+2);
Zpath      = exp(Ztildepath);

% For predetermined variables
% first observation is the initial steady state value
Kpath    = oo_.endo_simul(3,1:Nperiods+1);


% *************** *
% Flexible prices *
% *************** *

zetaval = 0.00;
dynare MBC_Zshock noclearall

% For jumping variables:
% first observation is the initial steady state value:
YinitialFP = oo_.endo_simul(1,1);
LinitialFP = oo_.endo_simul(4,1);
ZinitialFP = oo_.endo_simul(15,1);

% second observation is the impact value:
YpathFP      = oo_.endo_simul(1,2:Nperiods+2);
YapathFP     = oo_.endo_simul(2,2:Nperiods+2);
LpathFP      = oo_.endo_simul(4,2:Nperiods+2);
CpathFP      = oo_.endo_simul(5,2:Nperiods+2);
wpathFP      = oo_.endo_simul(6,2:Nperiods+2);
rkpathFP     = oo_.endo_simul(7,2:Nperiods+2);
RpathFP      = oo_.endo_simul(8,2:Nperiods+2);
mcpathFP     = oo_.endo_simul(9,2:Nperiods+2);
PpathFP      = oo_.endo_simul(10,2:Nperiods+2);
PapathFP     = oo_.endo_simul(11,2:Nperiods+2);
PnpathFP     = oo_.endo_simul(12,2:Nperiods+2);
XI_NpathFP   = oo_.endo_simul(13,2:Nperiods+2);
XI_DpathFP   = oo_.endo_simul(14,2:Nperiods+2);
ZtildepathFP = oo_.endo_simul(15,2:Nperiods+2);
ZpathFP      = exp(ZtildepathFP); 

% For predetermined variables
% first observation is the initial steady state value
KpathFP    = oo_.endo_simul(3,1:Nperiods+1);

% Some graphs 

figure(1)
hold on;
plot(PStime,100*(Zpath - EXO.Z0)/EXO.Z0,'.','MarkerSize',18)
xlim([0,50])
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{Z}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Productivity_Zshock')
pause(1)

figure(2)
hold on;
plot(PStime,100*(Ypath - VAR.Y0)/VAR.Y0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(YpathFP - VAR.Y0)/VAR.Y0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{Y}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Output_Zshock')
pause(1)

figure(3)
hold on;
plot(PStime,100*(Kpath - VAR.K0)/VAR.K0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(KpathFP - VAR.K0)/VAR.K0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{K}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Capital_Zshock')
pause(1)

figure(4)
hold on;
plot(PStime,100*(Lpath - VAR.L0)/VAR.L0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(LpathFP - VAR.L0)/VAR.L0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{L}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Employment_Zshock')
pause(1)

figure(5)
hold on;
plot(PStime,100*(Cpath - VAR.C0)/VAR.C0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(CpathFP - VAR.C0)/VAR.C0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{C}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Consumption_Zshock')
pause(1)

figure(6)
hold on;
plot(PStime,100*(wpath - VAR.w0)/VAR.w0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(wpathFP - VAR.w0)/VAR.w0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{w}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Wage_Zshock')
pause(1)

figure(7)
hold on;
plot(PStime,100*(rkpath - VAR.rk0),'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(rkpathFP - VAR.rk0),'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\left[r^{K}_{t} - \left(r^{K}\right)^{*}\right] \times 100$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Rentalcapital_Zshock')
pause(1)

figure(8)
hold on;
plot(PStime,100*(mcpath - VAR.mc0)/VAR.mc0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(mcpathFP - VAR.mc0)/VAR.mc0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{\mathit{mc}}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','mc_Zshock')
pause(1)

figure(9)
hold on;
plot(PStime,100*(Ppath - VAR.P0)/VAR.P0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(PpathFP - VAR.P0)/VAR.P0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{P}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Pricelevel_Zshock')
pause(1)

figure(10)
hold on;
plot(PStime,100*(Rpath - PAR.rho)/(1+PAR.rho),'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(RpathFP - PAR.rho)/(1+PAR.rho),'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{R}_{t} \times 100 $','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Interestrate_Zshock')
pause(1)

figure(11)
hold on;
plot(PStime,100*(Pnpath - VAR.P0)/VAR.P0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(PnpathFP - VAR.P0)/VAR.P0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{P}^{n}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Newprice_Zshock')
pause(1)


    case 2
%% Money supply shock

% ************* *
% Sticky prices *
% ************* *

zetaval = 0.75;
dynare MBC_Mshock noclearall

% For jumping variables:
% first observation is the initial steady state value:
Yinitial = oo_.endo_simul(1,1);
Linitial = oo_.endo_simul(4,1);
Zinitial = oo_.endo_simul(15,1);

% second observation is the impact value:
Ypath      = oo_.endo_simul(1,2:Nperiods+2);
Yapath     = oo_.endo_simul(2,2:Nperiods+2);
Lpath      = oo_.endo_simul(4,2:Nperiods+2);
Cpath      = oo_.endo_simul(5,2:Nperiods+2);
wpath      = oo_.endo_simul(6,2:Nperiods+2);
rkpath     = oo_.endo_simul(7,2:Nperiods+2);
Rpath      = oo_.endo_simul(8,2:Nperiods+2);
mcpath     = oo_.endo_simul(9,2:Nperiods+2);
Ppath      = oo_.endo_simul(10,2:Nperiods+2);
Papath     = oo_.endo_simul(11,2:Nperiods+2);
Pnpath     = oo_.endo_simul(12,2:Nperiods+2);
XI_Npath   = oo_.endo_simul(13,2:Nperiods+2);
XI_Dpath   = oo_.endo_simul(14,2:Nperiods+2);
Mtildepath = oo_.endo_simul(15,2:Nperiods+2);
Mpath      = exp(Mtildepath);

% For predetermined variables
% first observation is the initial steady state value
Kpath    = oo_.endo_simul(3,1:Nperiods+1);

% *************** *
% Flexible prices *
% *************** *

zetaval = 0.00;
dynare MBC_Mshock noclearall

% For jumping variables:
% first observation is the initial steady state value:
YinitialFP = oo_.endo_simul(1,1);
LinitialFP = oo_.endo_simul(4,1);
MinitialFP = oo_.endo_simul(15,1);

% second observation is the impact value:
YpathFP      = oo_.endo_simul(1,2:Nperiods+2);
YapathFP     = oo_.endo_simul(2,2:Nperiods+2);
LpathFP      = oo_.endo_simul(4,2:Nperiods+2);
CpathFP      = oo_.endo_simul(5,2:Nperiods+2);
wpathFP      = oo_.endo_simul(6,2:Nperiods+2);
rkpathFP     = oo_.endo_simul(7,2:Nperiods+2);
RpathFP      = oo_.endo_simul(8,2:Nperiods+2);
mcpathFP     = oo_.endo_simul(9,2:Nperiods+2);
PpathFP      = oo_.endo_simul(10,2:Nperiods+2);
PapathFP     = oo_.endo_simul(11,2:Nperiods+2);
PnpathFP     = oo_.endo_simul(12,2:Nperiods+2);
XI_NpathFP   = oo_.endo_simul(13,2:Nperiods+2);
XI_DpathFP   = oo_.endo_simul(14,2:Nperiods+2);
MtildepathFP = oo_.endo_simul(15,2:Nperiods+2);
MpathFP      = exp(MtildepathFP);


% For predetermined variables
% first observation is the initial steady state value
KpathFP    = oo_.endo_simul(3,1:Nperiods+1);

% Some graphs 


figure(1)
hold on;
plot(PStime,100*(Mpath - EXO.M_0)/EXO.M_0,'.','MarkerSize',18)
xlim([0,50])
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\Delta M_{t+1} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Productivity_Mshock')
pause(1)

figure(2)
hold on;
plot(PStime,100*(Ypath - VAR.Y0)/VAR.Y0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(YpathFP - VAR.Y0)/VAR.Y0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{Y}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Output_Mshock')
pause(1)

figure(3)
hold on;
plot(PStime,100*(Kpath - VAR.K0)/VAR.K0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(KpathFP - VAR.K0)/VAR.K0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{K}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Capital_Mshock')
pause(1)

figure(4)
hold on;
plot(PStime,100*(Lpath - VAR.L0)/VAR.L0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(LpathFP - VAR.L0)/VAR.L0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{L}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Employment_Mshock')
pause(1)

figure(5)
hold on;
plot(PStime,100*(Cpath - VAR.C0)/VAR.C0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(CpathFP - VAR.C0)/VAR.C0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{C}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Consumption_Mshock')
pause(1)

figure(6)
hold on;
plot(PStime,100*(wpath - VAR.w0)/VAR.w0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(wpathFP - VAR.w0)/VAR.w0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{w}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Wage_Mshock')
pause(1)

figure(7)
hold on;
plot(PStime,100*(rkpath - VAR.rk0),'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(rkpathFP - VAR.rk0),'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\left[r^{K}_{t} - \left(r^{K}\right)^{*}\right] \times 100$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Rentalcapital_Mshock')
pause(1)

figure(8)
hold on;
plot(PStime,100*(mcpath - VAR.mc0)/VAR.mc0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(mcpathFP - VAR.mc0)/VAR.mc0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{\mathit{mc}}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','mc_Mshock')
pause(1)

figure(9)
hold on;
plot(PStime,100*(Ppath - VAR.P0)/VAR.P0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(PpathFP - VAR.P0)/VAR.P0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{P}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Pricelevel_Mshock')
pause(1)

figure(10)
hold on;
plot(PStime,100*(Rpath - PAR.rho)/(1+PAR.rho),'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(RpathFP - PAR.rho)/(1+PAR.rho),'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{R}_{t} \times 100 $','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Interestrate_Mshock')
pause(1)

figure(11)
hold on;
plot(PStime,100*(Pnpath - VAR.P0)/VAR.P0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(PnpathFP - VAR.P0)/VAR.P0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{P}^{n}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Newprice_Mshock')
pause(1)



    case 3
%% Government consumption shock

% ************* *
% Sticky prices *
% ************* *

zetaval = 0.75;
dynare MBC_Gshock noclearall

% For jumping variables:
% first observation is the initial steady state value:
Yinitial = oo_.endo_simul(1,1);
Linitial = oo_.endo_simul(4,1);
Zinitial = oo_.endo_simul(15,1);

% second observation is the impact value:
Ypath    = oo_.endo_simul(1,2:Nperiods+2);
Yapath   = oo_.endo_simul(2,2:Nperiods+2);
Lpath    = oo_.endo_simul(4,2:Nperiods+2);
Cpath    = oo_.endo_simul(5,2:Nperiods+2);
wpath    = oo_.endo_simul(6,2:Nperiods+2);
rkpath   = oo_.endo_simul(7,2:Nperiods+2);
Rpath    = oo_.endo_simul(8,2:Nperiods+2);
mcpath   = oo_.endo_simul(9,2:Nperiods+2);
Ppath    = oo_.endo_simul(10,2:Nperiods+2);
Papath   = oo_.endo_simul(11,2:Nperiods+2);
Pnpath   = oo_.endo_simul(12,2:Nperiods+2);
XI_Npath = oo_.endo_simul(13,2:Nperiods+2);
XI_Dpath = oo_.endo_simul(14,2:Nperiods+2);
Gtildepath = oo_.endo_simul(15,2:Nperiods+2);
Gpath    = 0.2 * exp(Gtildepath);

% For predetermined variables
% first observation is the initial steady state value
Kpath    = oo_.endo_simul(3,1:Nperiods+1);


% *************** *
% Flexible prices *
% *************** *

zetaval = 0.00;
dynare MBC_Gshock noclearall

% For jumping variables:
% first observation is the initial steady state value:
YinitialFP  = oo_.endo_simul(1,1);
LinitialFP  = oo_.endo_simul(4,1);
GinitialFP  = oo_.endo_simul(15,1);

% second observation is the impact value:
YpathFP      = oo_.endo_simul(1,2:Nperiods+2);
YapathFP     = oo_.endo_simul(2,2:Nperiods+2);
LpathFP      = oo_.endo_simul(4,2:Nperiods+2);
CpathFP      = oo_.endo_simul(5,2:Nperiods+2);
wpathFP      = oo_.endo_simul(6,2:Nperiods+2);
rkpathFP     = oo_.endo_simul(7,2:Nperiods+2);
RpathFP      = oo_.endo_simul(8,2:Nperiods+2);
mcpathFP     = oo_.endo_simul(9,2:Nperiods+2);
PpathFP      = oo_.endo_simul(10,2:Nperiods+2);
PapathFP     = oo_.endo_simul(11,2:Nperiods+2);
PnpathFP     = oo_.endo_simul(12,2:Nperiods+2);
XI_NpathFP   = oo_.endo_simul(13,2:Nperiods+2);
XI_DpathFP   = oo_.endo_simul(14,2:Nperiods+2);
GtildepathFP = oo_.endo_simul(15,2:Nperiods+2);
GpathFP      = 0.2 * exp(GtildepathFP);

% For predetermined variables
% first observation is the initial steady state value
KpathFP    = oo_.endo_simul(3,1:Nperiods+1);

% Some graphs 


figure(1)
hold on;
plot(PStime,100*(Gpath - EXO.G)/EXO.G,'.','MarkerSize',18)
xlim([0,50])
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{G}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Governmentconsumption_Gshock')
pause(1)

figure(2)
hold on;
plot(PStime,100*(Ypath - VAR.Y0)/VAR.Y0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(YpathFP - VAR.Y0)/VAR.Y0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{Y}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Output_Gshock')
pause(1)

figure(3)
hold on;
plot(PStime,100*(Kpath - VAR.K0)/VAR.K0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(KpathFP - VAR.K0)/VAR.K0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{K}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Capital_Gshock')
pause(1)

figure(4)
hold on;
plot(PStime,100*(Lpath - VAR.L0)/VAR.L0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(LpathFP - VAR.L0)/VAR.L0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{L}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Employment_Gshock')
pause(1)

figure(5)
hold on;
plot(PStime,100*(Cpath - VAR.C0)/VAR.C0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(CpathFP - VAR.C0)/VAR.C0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{C}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Consumption_Gshock')
pause(1)

figure(6)
hold on;
plot(PStime,100*(wpath - VAR.w0)/VAR.w0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(wpathFP - VAR.w0)/VAR.w0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{w}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Wage_Gshock')
pause(1)

figure(7)
hold on;
plot(PStime,100*(rkpath - VAR.rk0),'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(rkpathFP - VAR.rk0),'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\left[r^{K}_{t} - \left(r^{K}\right)^{*}\right] \times 100$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Rentalcapital_Gshock')
pause(1)

figure(8)
hold on;
plot(PStime,100*(mcpath - VAR.mc0)/VAR.mc0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(mcpathFP - VAR.mc0)/VAR.mc0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{\mathit{mc}}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','mc_Gshock')
pause(1)

figure(9)
hold on;
plot(PStime,100*(Ppath - VAR.P0)/VAR.P0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(PpathFP - VAR.P0)/VAR.P0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{P}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Pricelevel_Gshock')
pause(1)

figure(10)
hold on;
plot(PStime,100*(Rpath - PAR.rho)/(1+PAR.rho),'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(RpathFP - PAR.rho)/(1+PAR.rho),'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{R}_{t} \times 100 $','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Interestrate_Gshock')
pause(1)

figure(11)
hold on;
plot(PStime,100*(Pnpath - VAR.P0)/VAR.P0,'.','MarkerSize',18)
xlim([0,50])
plot(PStime,100*(PnpathFP - VAR.P0)/VAR.P0,'o','MarkerSize',5)
hold off;
xlabel('post-shock time $t$','interpreter','latex','FontSize',14)
ylabel('$\tilde{P}^{n}_{t} \times 100 \%$','interpreter','latex','FontSize',14)
box off
print('-deps2','-loose','Newprice_Gshock')
pause(1)

    case 4
        
 % Check roots of the linearized model in the B-K format
 
 Gamma     = eye(6,6);
 Deltastar = eye(6,6);
 alpha     = PAR.alpha;
 delta     = PAR.delta;
 rho       = PAR.rho;
 phi       = PAR.phi;
 zeta      = PAR.zeta;
 xi_Z      = PAR.xi_Z;
 
 omega_C    = VAR.C0/VAR.Y0;
 omega_I    = delta * VAR.K0/VAR.Y0;
 omega_G    = 1 - omega_C - omega_I;
 omega_LL   = (1-VAR.L0)/VAR.L0;
 omega_Y    = (1-phi+alpha*omega_LL)/((1-alpha) * omega_LL);
 omega_K    = (1-phi+omega_LL)/((1-alpha) * omega_LL);
 gamma      = (1-zeta) * (1-zeta+rho) / zeta;
 A          = (delta+rho) / (1+rho);
 
 Gamma(3,1) = -rho;
 Gamma(3,5) = 1;
 Gamma(4,2) =  A * (1 + alpha*omega_K);
 Gamma(4,3) =  (1-delta)/(1+rho);
 Gamma(4,4) = -A * omega_K;
 Gamma(4,6) = A * (1-alpha)*omega_K;

 Gamma;
 
 Gaminv     = inv(Gamma);
 
 Deltastar(1,5) = 1;
 Deltastar(2,2) = 1-delta; 
 Deltastar(2,3) = -(delta * omega_C / omega_I) ; 
 Deltastar(2,4) = delta/omega_I ; 
 Deltastar(3,3) = 1+rho ; 
 Deltastar(4,3) = 1 ; 
 Deltastar(4,4) = 0 ;
 Deltastar(5,2) = alpha * gamma * omega_K;
 Deltastar(5,3) = -gamma;
 Deltastar(5,4) = -gamma * (omega_K-1);
 Deltastar(5,5) = 1+rho ; 
 Deltastar(5,6) = (1- alpha) * gamma * omega_K;
 Deltastar(6,6) = xi_Z;
 
 Deltastar;

 Delta = Gaminv * Deltastar;
 Ben = eig(Delta);
 
 Analytic = eye(6,6);
 Analytic(1,5) = 1;
 Analytic(2,2) = 1 - delta;
 Analytic(2,3) = -(delta * omega_C / omega_I);
 Analytic(2,4) = delta/omega_I;
 Analytic(3,1) = rho;
 Analytic(3,2) = -alpha * gamma * omega_K;
 Analytic(3,3) = 1 + rho + gamma;
 Analytic(3,3) = 1 + rho + gamma;
 Analytic(3,4) = gamma * omega_Y ;
 Analytic(3,5) = -1;
 Analytic(3,6) = -(1-alpha)*gamma*omega_K;
 Analytic(4,1) = rho * ((1-A)/A) /(1+omega_Y)  ;
 Analytic(4,2) = (alpha*gamma*omega_K* ((A-1)/A) + (1-delta)*(alpha*omega_K +1) ) ...
  /(1+omega_Y);
 Analytic(4,3) = ((gamma+rho+1)*(1-A) - 1)/(A*(1+omega_Y)) ...
     - (delta * omega_C / omega_I) * (1+alpha*omega_K) / (1+omega_Y) ;
 Analytic(4,4) = (delta / omega_I) * (1+alpha*omega_K) / (1+omega_Y) ...
     - gamma * omega_Y * ((A-1)/A) / (1+omega_Y) ;
 Analytic(4,5) = ((A-1)/A)/(1+omega_Y) ;
 Analytic(4,6) = gamma*omega_K*(1-alpha)*((A-1)/A) /(1+omega_Y) ...
     + (1-alpha)*omega_K*xi_Z / (1+omega_Y) ;
 Analytic(5,2) = alpha * gamma * omega_K; 
 Analytic(5,3) = -gamma;
 Analytic(5,4) = -gamma*omega_Y;
 Analytic(5,5) = 1 + rho;
 Analytic(5,6) = gamma*omega_K *(1-alpha);
 Analytic(6,6) = xi_Z;
 
 Jitka = eig(Analytic)
 
 Analytic
 
 Check1 = (1-delta) * ( ... 
 (Delta(4,1)-Delta(4,5))*gamma*(1-omega_K) + Delta(4,4)*(1+rho)) ...
 + delta/omega_I * ((Delta(4,5)-Delta(4,1))*alpha*gamma*omega_K ...
 - Delta(4,2)*(1+rho)) ;
 
phi_a2 = -(2 + rho + gamma - delta + Delta(4,4))
phi_a1 = -gamma*(1-omega_K)*(Delta(4,5)-Delta(4,3)) ...
    + (1 + rho + gamma) *Delta(4,4) ...
    + (1-delta)*(1 + rho + gamma + Delta(4,4)) ...
    - (delta/omega_I) *(alpha*gamma*omega_K*omega_C + Delta(4,2))
phi_a0 = (1-delta)* (gamma*(1-omega_K)*(Delta(4,5)-Delta(4,3)) ...
    - (1 + rho + gamma) * Delta(4,4)) ...
    + (delta*omega_C/omega_I) * (alpha*gamma*omega_K*Delta(4,4) ...
    - gamma*(1-omega_K)*Delta(4,2)) ...
    + (delta/omega_I) * ((1 + rho + gamma)*Delta(4,2) ...
    - alpha*gamma*omega_K*(Delta(4,5)-Delta(4,3)))

phi_b2 = 2 + rho - delta + Delta(4,4)
phi_b1 = gamma*(1-omega_K)*(Delta(4,5)-Delta(4,1)) ...
    -  (1 + rho) *Delta(4,4) - (1-delta)*(1 + rho + Delta(4,4)) ...
    +  (delta/omega_I) * Delta(4,2) 

rootprod = phi_a0 + phi_b1
rootprod + Jitka(1,1) * Jitka(4,1) * Jitka(5,1)

p = [1 (phi_a2-1) (phi_a1+phi_b2) (phi_a0+phi_b1)]
Pietje = roots(p)
Pietje(1,1)* Pietje(2,1) * Pietje(3,1);

del_41 = (1-delta)*rho / ((delta+rho)*omega_K);
del_42 = (1-delta)*(alpha*omega_K*(delta+rho-gamma) + (delta+rho))  / ((delta+rho)*omega_K);
del_43 = ((1-delta)*gamma - delta*(1+rho)) /  ((delta+rho)*omega_K) ...
    - (delta*omega_C/omega_I) * (1+alpha*omega_K) / omega_K;
del_44 = (delta/omega_I)*(1+alpha*omega_K)/omega_K + gamma*(1-delta)*(omega_K-1)/ ((delta+rho)*omega_K);
del_45 = -(1-delta) / ((delta+rho)*omega_K);
del_46 = (1-alpha)*(xi_Z - gamma*(1-delta)/(delta+rho));

fie_b1 = -(delta/(omega_I*(rho+delta)*omega_K)) ...
    * ( alpha*omega_K * ( gamma*(1-delta) + (1+rho)*(rho+delta)) ...
    + (delta+rho)*(1+rho) ) ...
    - (1-delta) * ((gamma*(1-delta) + (1+rho)*(rho+delta))*(omega_K - 1) + (1+rho)*(rho+delta) ) ...
    /((rho+delta)*omega_K) ;

fie_b2 = (delta/omega_I) * (1+alpha*omega_K)/omega_K ...
    + (1/((rho+delta)*omega_K)) *( omega_K*(2+rho-delta)*(delta+rho) ...
    + gamma*(1-delta)*(omega_K-1) ) ;

fie_a0 = -alpha*gamma*delta/omega_I - gamma*(1-delta)*(omega_K-1)/omega_K;

fie_a1 = (delta/(omega_I*(rho+delta)*omega_K)) * ...
( alpha*omega_K * ...
( gamma*(1-delta) + (delta+rho)*(1+rho+gamma*(1-omega_C)))  ... 
+(delta+rho)*(1+rho+gamma*(1+(omega_K-1)*omega_C)) ) ...
+ (1/((rho+delta)*omega_K)) * ...
( (omega_K-1)*(gamma*(1+rho*(2-delta)) ...
+ (1-delta)*(1+rho)*(delta+rho) ) + (1-delta)*(1+rho+gamma)*(delta+rho) );

fie_a2 = -(delta/omega_I) * (1+alpha*omega_K)/omega_K ...
    - (1/((rho+delta)*omega_K)) *( omega_K*(2+rho+gamma-delta)*(delta+rho) ...
    + gamma*(1-delta)*(omega_K-1) ) ;

Sofie = [fie_a0 fie_a1 fie_a2 fie_b1 fie_b2 ;
    phi_a0 phi_a1 phi_a2 phi_b1 phi_b2]'
end;
    