%Hybrid NKM under commitment
%% Step 1 define the parameters
%same Parameters as from the script
%Parameters
gam_f = 0.5
gam_b =0.5
beta= 1
gam_x = 0.2
lambda = 0.5
phi = 0.5                       % IS Parameter for interest rate
sig = 1                         % IS Parameter for interest rate
%AR_par = 0.8
AR_par1 = [0.8 0];
figure
for n = 1:2
      AR_par = AR_par1(n);

%% Step 2 
% System is: (w;v)(+1) = A*[w;v] + [ 1;0;0;0;0]*eps
A11= [AR_par 0 0;0 0 0; 0 0 0];
A12= [ 0 0; 1 0; 0 1];
A21 =[ 0 -gam_f/(beta^2*gam_b) 0; -1/(beta*gam_f) 0 -gam_b/gam_f];   %Klammern nicht vergessen
A22 =[ 1/(beta^2*gam_b) gam_x/(beta^2*lambda*gam_b); -gam_x/(beta*gam_f) 1/(beta*gam_f)];

A = [ [A11 A12] ; [A21 A22] ]

%% Step 3
% using the Schur Decomposition to solve the state equations
% solve the system
disp('Schur decomposition')
 [Z, T] = schur(A, 'complex')
  disp('reorder eigenvalues in increasing order along the principal diagonal')
 [Z T] = ordschur(Z,T, 1:5)

 if abs(sum(sum(Z*T*Z'-A))) > 0.0001 && sum(sum(Z'*Z-eye(lenght(Z)))) > 0.0001 
disp('Error in Schur decomposition')
end

 disp('check Blanchard-Kahn')
 abs_diag_T = abs(diag(T))';
 

%% Calculating the solution time path for nu, x and pi using the following law of motion:
% z(+1) = E[z(+1)] + Z_11^-1 * [1;0;0] * eps
T_11 = T(1:3,1:3)
Z_11 = Z(1:3,1:3)
Z_21 = Z(4:5,1:3)
T=1000000;
z_solution= zeros(3,T);                     % zeros because we have variables in  i with t-1
w_solution= zeros(3,T);                     % zeros because we have variables in i with t-1
w_solution(:,1)=[ 1; 0; 0];    %initial jump
z_solution(:,1)=inv(Z_11)* w_solution(: ,1);%initial jump
v_solution= zeros(2,T);
i_solution= zeros(1,T);           % nominal interest rate: IS umgestellt nach der Variable i.Hier liegt anscheinend das Problem.
for t= 2:T
z_solution(:,t)= T_11* z_solution(: ,t-1 );
w_solution(:,t)= Z_11 * z_solution(:,t);
v_solution(:,t)= Z_21 *inv(Z_11)* w_solution(:,t);
end
for t= 1:T-1
    i_solution(:,t) =((1- phi)*v_solution(1,t+1)+phi*w_solution(2,t)-w_solution(2,t+1))*sig+ v_solution(2,t+1); % Jump in 1 anstatt 2. umgestellte IS-Kurve.
end;

%% plots
t=0: 20;
w_solution = real(w_solution);
v_solution = real(v_solution);
i_solution = real(i_solution);

lineStyle = {'k->', 'k-o'};
subplot(2,2,1);  hold on;
plot(t,  w_solution(2,t+2), lineStyle{n}, 'LineWidth',0.5)
xlabel(' t '); ylabel (' x '); title(' Output gap ')
axis([0 20 -4 2])
set(gca,'XTickLabel',{'0','5','10','15','20'})
box on

%figure
%lineStyle = {'k->', 'k-o'};
%subplot(2,2,1);  hold on;
%plot(t,  w_solution(2,t+2), lineStyle{n}, 'LineWidth',2)
%xlabel(' t '); ylabel (' x '); title(' Output gap ')
%axis([0 20 -4 2])
%set(gca,'XTickLabel',{'0','5','10','15','20'})
%box on

subplot(2,2,2); hold on;
plot(t, w_solution(3,t+2),lineStyle{n},  'LineWidth',0.5)
xlabel(' t ');  ylabel(' \pi '); title('Inflation')
axis([0 20 -0.5 1.5])
set(gca,'XTickLabel',{'0','5','10','15','20'})
box on

subplot(2,2,3); hold on;
plot(t, i_solution(t+1),lineStyle{n},'Linewidth',0.5)
xlabel('t'); ylabel('i'); title('Interest rate')
axis([0 20 -1 3])
set(gca,'XTickLabel',{'0','5','10','15','20'})
box on

subplot(2,2,4); hold on;
plot(t, w_solution(1,t+1),lineStyle{n},'Linewidth',0.5)
xlabel('t'); ylabel('nu'); title('Shock')
legend('\rho = 0','\rho = 0.8') 
box on
end
 
%data = struct('v_solution',[1,T;2,T]); 
%T=1000000;
%w_pi = 0.5 ;   
%w_x = 0.25; 
%[LossVal,vol_pi, vol_x ] = Loss_fun( data,w_pi, w_x )

w_pi = 0.5 ;   
w_x = 0.25; 
LossVal = w_pi *var(v_solution(2,:)) + w_x *var(v_solution(1,:))
%variance.x=oo_.var(x_pos,x_pos);
%variance.pi=oo_.var(pi_pos,pi_pos);

%loss_val = sum(alp1*ppi_epsu.^2 + alp2*x_epsu.^2 + alp3*int_epsu.^2)