%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Shocks and Optimal Policy Response in High-Frequency New-Keynesian Models%
%%%%%%%%%%%%%by Stephen Sacht & Hans-Werner Wohltmann (2013)%%%%%%%%%%%%%%%
%%%Matlab Source Code (c) by Jan-Niklas Brenneisen & Stephen Sacht (2013)%%
%%%%%%%%%%%%%%%Representations of the Maximum Value of h %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%under Variation of the Shock Persistence%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%(Adjusted for Skip Sampling Aggregation Scheme)%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Structural Parameters
for j=1:99
phi=0.01*j;
h=colon(1/90,1/180,1).'; %[1;1/2;1/3;1/12;1/90];
alpha_2=0.05; %[0.5,0.01]
r_bar=0.0101;
theta=0.6667;
sigma=1; %[1,2]
eta=sigma;
epsilon_0=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Composite Parameters
for i=1:179
K_tilde(i,1)=(sigma+eta)*(1-theta)*h(i,1)^2*(1-theta+r_bar)/((1-(1-theta)*h(i,1))*(1+h(i,1)*r_bar));
alpha_h(i,1)=1-h(i,1)*(1-phi);
beta_h(i,1)=1/(1+h(i,1)*r_bar);
Psi(i,1)=(1-alpha_h(i,1)*beta_h(i,1))*alpha_2+K_tilde(i,1)^2/h(i,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Regime Commitment Parameters
lambda_h(i,1)=1/beta_h(i,1)+beta_h(i,1)^(1-h(i,1))+K_tilde(i,1)^2/(h(i,1)*alpha_2*beta_h(i,1));
C([2*i-1,2*i],[1,2])=[lambda_h(i,1),-beta_h(i,1)^-h(i,1);1,0];
E(i,[1,2])=eig(C([2*i-1,2*i],[1,2]));
r_1(i,1)=E(i,1);
r_2(i,1)=E(i,2);
phi_1(i,1)=alpha_2+h(i,1)*alpha_2^2/(K_tilde(i,1)^2);
phi_2(i,1)=h(i,1)*(alpha_2^2/K_tilde(i,1)^2)*beta_h(i,1)^(1-h(i,1));
phi_3(i,1)=h(i,1)*(alpha_2^2/K_tilde(i,1)^2)*beta_h(i,1)^(2*(1-h(i,1)));
gamma_1(i,1)=phi_1(i,1)*r_2(i,1)^2-2*phi_2(i,1)*r_2(i,1)+phi_3(i,1);
gamma_2(i,1)=phi_1(i,1)*alpha_h(i,1)*r_2(i,1)-phi_2(i,1)*(alpha_h(i,1)+r_2(i,1))+phi_3(i,1);
gamma_3(i,1)=phi_1(i,1)*alpha_h(i,1)^2-2*phi_2(i,1)*alpha_h(i,1)+phi_3(i,1);
W([2*i-1,2*i],[1,2])=[gamma_1(i,1)/(1-beta_h(i,1)*r_2(i,1)^(2/h(i,1))), gamma_2(i,1)/(1-beta_h(i,1)*(r_2(i,1)*alpha_h(i,1))^(1/h(i,1)));gamma_2(i,1)/(1-beta_h(i,1)*(r_2(i,1)*alpha_h(i,1))^(1/h(i,1))), gamma_3(i,1)/(1-beta_h(i,1)*alpha_h(i,1)^(2/h(i,1)))];
h_33_inv(i,1)=1/(alpha_h(i,1)^2-lambda_h(i,1)*alpha_h(i,1)+beta_h(i,1)^-h(i,1))*K_tilde(i,1)/(alpha_2*beta_h(i,1));
H_22_inv([2*i-1,2*i],[1,2])=[1, -h_33_inv(i,1);0,h_33_inv(i,1)];
q_tilde([2*i-1,2*i],1)=H_22_inv([2*i-1,2*i],[1,2])*[0;1];
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loss Functions
for i=1:179
V_D(i,1)=(alpha_2*(h(i,1)*alpha_2+K_tilde(i,1)^2))/Psi(i,1)^2*1/(1-beta_h(i,1)*alpha_h(i,1)^(2/h(i,1)))*epsilon_0^2;
V_OSR(i,1)=h(i,1)*alpha_2/(K_tilde(i,1)^2/h(i,1)+alpha_2*(1-alpha_h(i,1)*beta_h(i,1))^2)*1/(1-beta_h(i,1)*alpha_h(i,1)^(2/h(i,1)))*epsilon_0^2;
V_C(i,1)=trace(W([2*i-1,2*i],[1,2])*q_tilde([2*i-1,2*i],1)*q_tilde([2*i-1,2*i],1).')*epsilon_0^2;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loss Function Relations
for i=1:179
V_D_C(i,1)=V_D(i,1)-V_C(i,1);
V_D_OSR(i,1)=V_D(i,1)-V_OSR(i,1);
V_OSR_C(i,1)=V_OSR(i,1)-V_C(i,1);
V_D_C_rel(i,1)=(V_D(i,1)-V_C(i,1))/V_C(i,1);
V_OSR_C_rel(i,1)=(V_OSR(i,1)-V_C(i,1))/V_C(i,1);
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Display Maximum Values of h
[A,B(j,1)]=max(V_OSR_C_rel);
max_V_OSR_C_rel(j,1)=h(B(j,1),1);
[C,D(j,1)]=max(V_D_C_rel);
max_V_D_C_rel(j,1)=h(D(j,1),1);
[E,F(j,1)]=min(V_OSR_C_rel);
min_V_OSR_C_rel(j,1)=h(F(j,1),1);
[G,H(j,1)]=min(V_D_C_rel);
min_V_D_C_rel(j,1)=h(H(j,1),1);
phi_st=colon(0,0.01,0.98).';
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
plot(phi_st,max_V_OSR_C_rel)
legend('(V^{OSR}-V^C)/V^C');
figure(2);
plot(phi_st,max_V_D_C_rel)
legend('(V^D-V^C)/V^C');
figure(3);
plot(phi_st,min_V_OSR_C_rel)
legend('(V^{OSR}-V^C)/V^C');
figure(4);
plot(phi_st,min_V_D_C_rel)
legend('(V^D-V^C)/V^C');