%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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 and the Calvo Parameter%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%(Adjusted for Skip Sampling Aggregation Scheme)%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Structural Parameters
for j=1:99
phi=0.01*j;
theta=colon(0.5,0.004,0.895);
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=0.0101;
sigma=1; %[1,2]
eta=sigma;
epsilon_0=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Composite Parameters
for t=1:99
for i=1:179
K_tilde(i,t)=(sigma+eta)*(1-theta(1,t))*h(i,1)^2*(1-theta(1,t)+r_bar)/((1-(1-theta(1,t))*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,t)=(1-alpha_h(i,1)*beta_h(i,1))*alpha_2+K_tilde(i,t)^2/h(i,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Regime Commitment Parameters
lambda_h(i,t)=1/beta_h(i,1)+beta_h(i,1)^(1-h(i,1))+K_tilde(i,t)^2/(h(i,1)*alpha_2*beta_h(i,1));
C([2*i-1,2*i],[2*t-1,2*t])=[lambda_h(i,t),-beta_h(i,1)^-h(i,1);1,0];
E(i,[2*t-1,2*t])=eig(C([2*i-1,2*i],[2*t-1,2*t]));
r_1(i,t)=E(i,2*t-1);
r_2(i,t)=E(i,2*t);
phi_1(i,t)=alpha_2+h(i,1)*alpha_2^2/K_tilde(i,t)^2;
phi_2(i,t)=h(i,1)*alpha_2^2/K_tilde(i,t)^2*beta_h(i,1)^(1-h(i,1));
phi_3(i,t)=h(i,1)*alpha_2^2/K_tilde(i,t)^2*beta_h(i,1)^(2*(1-h(i,1)));
gamma_1(i,t)=phi_1(i,t)*r_2(i,t)^2-2*phi_2(i,t)*r_2(i,t)+phi_3(i,t);
gamma_2(i,t)=phi_1(i,t)*alpha_h(i,1)*r_2(i,t)-phi_2(i,t)*(alpha_h(i,1)+r_2(i,t))+phi_3(i,t);
gamma_3(i,t)=phi_1(i,t)*alpha_h(i,1)^2-2*phi_2(i,t)*alpha_h(i,1)+phi_3(i,t);
W([2*i-1,2*i],[2*t-1,2*t])=[gamma_1(i,t)/(1-beta_h(i,1)*r_2(i,t)^(2/h(i,1))), gamma_2(i,t)/(1-beta_h(i,1)*(r_2(i,t)*alpha_h(i,1))^(1/h(i,1)));gamma_2(i,t)/(1-beta_h(i,1)*(r_2(i,t)*alpha_h(i,1))^(1/h(i,1))), gamma_3(i,t)/(1-beta_h(i,1)*alpha_h(i,1)^(2/h(i,1)))];
h_33_inv(i,t)=1/(alpha_h(i,1)^2-lambda_h(i,t)*alpha_h(i,1)+beta_h(i,1)^-h(i,1))*K_tilde(i,t)/(alpha_2*beta_h(i,1));
H_22_inv([2*i-1,2*i],[2*t-1,2*t])=[1, -h_33_inv(i,t);0,h_33_inv(i,t)];
q_tilde([2*i-1,2*i],t)=H_22_inv([2*i-1,2*i],[2*t-1,2*t])*[0;1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loss Functions
V_D(i,t)=alpha_2*(h(i,1)*alpha_2+K_tilde(i,t)^2)/Psi(i,t)^2*(1/(1-beta_h(i,1)*alpha_h(i,1)^(2/h(i,1))))*epsilon_0^2;
V_OSR(i,t)=h(i,1)*alpha_2/(K_tilde(i,t)^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,t)=trace(W([2*i-1,2*i],[2*t-1,2*t])*q_tilde([2*i-1,2*i],t)*q_tilde([2*i-1,2*i],t).')*epsilon_0^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loss Function Relations
V_D_C(i,t)=V_D(i,t)-V_C(i,t);
V_D_OSR(i,t)=V_D(i,t)-V_OSR(i,t);
V_OSR_C(i,t)=V_OSR(i,t)-V_C(i,t);
V_D_C_rel(i,t)=(V_D(i,t)-V_C(i,t))/V_C(i,t);
V_OSR_C_rel(i,t)=(V_OSR(i,t)-V_C(i,t))/V_C(i,t);
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Display Maximum Values of h
[A,B(j,:)]=max(V_OSR_C_rel);
max_V_OSR_C_rel(j,:)=h(B(j,:),1);
[C,D(j,:)]=max(V_D_C_rel);
max_V_D_C_rel(j,:)=h(D(j,:),1);
phi_st=colon(0,0.01,0.98).';
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for a=1:99
phi_mat(:,a)=phi_st(:,1);
end;
for b=1:99
theta_mat(b,:)=theta(1,:);
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%3D-Plots
figure(1);
plot3(phi_mat,theta_mat,max_V_OSR_C_rel)
legend('(V^{OSR}-V^C)/V^C');
xlabel('\phi')
ylabel('\theta')
zlabel('(V^{OSR}-V^C)/V^C')
figure(2);
plot3(phi_mat,theta_mat,max_V_D_C_rel)
legend('(V^D-V^C)/V^C');
xlabel('\phi')
ylabel('\theta')
zlabel('(V^D-V^C)/V^C')