%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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)%%
%%%%%%%%%%%%%%%Computation of Loss Relations under Variation of%%%%%%%%%%%%
%%%%%%%%%%%%%%%%the Shock Persistence and the Calvo Parameter%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%(Adjusted for Skip Sampling Aggregation Scheme)%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Structural Parameters
h=colon(1/90,1/180,1).'; %[1;1/2;1/3;1/12;1/90];
alpha_2=1; %[0.5,0.01]
r_bar=1;
theta=colon(0.5,0.004,0.9);
sigma=1; %[1,2]
eta=sigma;
phi=0.5; %[0, 0.5, 0.8]
epsilon_0=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Composite parameters
for i=1:179
for j=1:101
K_tilde(i,j)=(sigma+eta)*(1-theta(1,j))*h(i,1)^2*(1-theta(1,j)+r_bar)/((1-(1-theta(1,j))*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,j)=(1-alpha_h(i,1)*beta_h(i,1))*alpha_2+K_tilde(i,j)^2/h(i,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Regime Commitment Parameters
lambda_h(i,j)=1/beta_h(i,1)+beta_h(i,1)^(1-h(i,1))+K_tilde(i,j)^2/(h(i,1)*alpha_2*beta_h(i,1));
C([2*i-1,2*i],[2*j-1,2*j])=[lambda_h(i,j),-beta_h(i,1)^-h(i,1);1,0];
E(i,[2*j-1,2*j])=eig(C([2*i-1,2*i],[2*j-1,2*j]));
r_1(i,j)=E(i,2*j-1);
r_2(i,j)=E(i,2*j);
phi_1(i,j)=alpha_2+h(i,1)*alpha_2^2/K_tilde(i,j)^2;
phi_2(i,j)=h(i,1)*alpha_2^2/K_tilde(i,j)^2*beta_h(i,1)^(1-h(i,1));
phi_3(i,j)=h(i,1)*alpha_2^2/K_tilde(i,j)^2*beta_h(i,1)^(2*(1-h(i,1)));
gamma_1(i,j)=phi_1(i,j)*r_2(i,j)^2-2*phi_2(i,j)*r_2(i,j)+phi_3(i,j);
gamma_2(i,j)=phi_1(i,j)*alpha_h(i,1)*r_2(i,j)-phi_2(i,j)*(alpha_h(i,1)+r_2(i,j))+phi_3(i,j);
gamma_3(i,j)=phi_1(i,j)*alpha_h(i,1)^2-2*phi_2(i,j)*alpha_h(i,1)+phi_3(i,j);
W([2*i-1,2*i],[2*j-1,2*j])=[gamma_1(i,j)/(1-beta_h(i,1)*r_2(i,j)^(2/h(i,1))), gamma_2(i,j)/(1-beta_h(i,1)*(r_2(i,j)*alpha_h(i,1))^(1/h(i,1)));gamma_2(i,j)/(1-beta_h(i,1)*(r_2(i,j)*alpha_h(i,1))^(1/h(i,1))), gamma_3(i,j)/(1-beta_h(i,1)*alpha_h(i,1)^(2/h(i,1)))];
h_33_inv(i,j)=1/(alpha_h(i,1)^2-lambda_h(i,j)*alpha_h(i,1)+beta_h(i,1)^-h(i,1))*K_tilde(i,j)/(alpha_2*beta_h(i,1));
H_22_inv([2*i-1,2*i],[2*j-1,2*j])=[1, -h_33_inv(i,j);0,h_33_inv(i,j)];
q_tilde([2*i-1,2*i],j)=H_22_inv([2*i-1,2*i],[2*j-1,2*j])*[0;1];
%W_1([2*i-1,2*i],[2*j-1,2*j])=[(r_2(i,j)^2-r_2(i,j)^(2/h(i,1)))*gamma_1(i,j)/((1-beta_h(i,1)*r_2(i,j)^(2/h(i,1)))*(1-r_2(i,j)^2)),(alpha_h(i,1)*r_2(i,j)-(alpha_h(i,1)*r_2(i,j))^(1/h(i,1)))*gamma_2(i,j)/((1-beta_h(i,1)*(alpha_h(i,1)*r_2(i,j))^(1/h(i,1)))*(1-alpha_h(i,1)*r_2(i,j)));(alpha_h(i,1)*r_2(i,j)-(alpha_h(i,1)*r_2(i,j))^(1/h(i,1)))*gamma_2(i,j)/((1-beta_h(i,1)*(alpha_h(i,1)*r_2(i,j))^(1/h(i,1)))*(1-alpha_h(i,1)*r_2(i,j))),(alpha_h(i,1)^2-alpha_h(i,1)^(2/h(i,1)))*gamma_3(i,j)/((1-beta_h(i,1)*alpha_h(i,1)^(2/h(i,1)))*(1-alpha_h(i,1)^2))];
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loss Functions
for i=1:179
for j=1:101
V_D(i,j)=(alpha_2*(h(i,1)*alpha_2+K_tilde(i,j)^2)/Psi(i,j)^2)*1/(1-beta_h(i,1)*alpha_h(i,1)^(2/h(i,1)))*epsilon_0^2;
V_OSR(i,j)=h(i,1)*alpha_2/(K_tilde(i,j)^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,j)=trace(W([2*i-1,2*i],[2*j-1,2*j])*q_tilde([2*i-1,2*i],j)*q_tilde([2*i-1,2*i],j).')*epsilon_0^2;
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loss Function Relations
for i=1:179
for j=1:101
V_D_C(i,j)=V_D(i,j)-V_C(i,j);
V_D_OSR(i,j)=V_D(i,j)-V_OSR(i,j);
V_OSR_C(i,j)=V_OSR(i,j)-V_C(i,j);
V_D_C_rel(i,j)=(V_D(i,j)-V_C(i,j))/V_C(i,j);
V_OSR_C_rel(i,j)=(V_OSR(i,j)-V_C(i,j))/V_C(i,j);
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:101
h_mat(:,i)=h(:,1);
end;
for j=1:179
theta_mat(j,:)=theta(1,:);
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plots der Loss Function Relations (3D-Plots)
figure(1);
%subplot(3,2,1);
plot3(h_mat,theta_mat,V_D_C);
legend('\alpha_2=0.5, \sigma=\eta=1, \phi=0.8');
title('V^D-V^C');
xlim([0.2 1]);
xlabel('h')
ylabel('\theta')
zlabel('V^D-V^C')
%subplot(3,2,2);
figure(2);
plot3(h_mat,theta_mat,V_D_OSR);
legend('\alpha_2=0.5, \sigma=\eta=1, \phi=0.8');
title('V^D-V^{OSR}');
xlim([0.2 1]);
xlabel('h')
ylabel('\theta')
zlabel('V^D-V^{OSR}')
%subplot(3,2,3);
figure(3);
plot3(h_mat,theta_mat,V_OSR_C);
legend('\alpha_2=0.5, \sigma=\eta=1, \phi=0.8');
title('V^{OSR}-V^C');
xlim([0.2 1]);
xlabel('h')
ylabel('\theta')
zlabel('V^{OSR}-V^C')
%subplot(3,2,4);
figure(4);
plot3(h_mat,theta_mat,V_D_C_rel);
legend('\alpha_2=0.5, \sigma=\eta=1, \phi=0.8');
title('(V^D-V^C)/V^C');
xlabel('h')
ylabel('\theta')
zlabel('(V^D-V^C)/V^C')
%subplot(3,2,5);
figure(5);
plot3(h_mat,theta_mat,V_OSR_C_rel);
legend('\alpha_2=0.5, \sigma=\eta=1, \phi=0.8');
title('(V^{OSR}-V^C)/V^C');
xlabel('h')
ylabel('\theta')
zlabel('(V^{OSR}-V^C)/V^C')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Display Maximum Values of h
[A,B]=max(V_OSR_C_rel);
max_V_OSR_C_rel=h(B,1)
[C,D]=max(V_D_C_rel);
max_V_D_C_rel=h(D,1)