clear;
% Parameters of the model (without policy parameters)
alp  = 1;  % share of rational agents
bet  = 0.99;
sig  = 0.157;
kap  = 0.024;
thet = 1.1;  % adaptive or extrapolative expectations

% policy parameters
del_pi_1 = 0;
del_pi_0 = 0;
del_pi_m1= 0;
del_y_1  = 0;
del_y_0  = 0;
del_y_m1 = 0;

% Grid of policy parameters
del_pi_space = linspace(0,2,200);
del_y_space  = linspace(0,2,200);
[del_pi_matrix,del_y_matrix] = meshgrid(del_pi_space,del_y_space);

% Number of forward-looking variables
Number_of_FLV = 2;

% Save the results: wheter we have determinacy = 1, indeterminacy = 0, insability = -1.
Results = NaN(size(del_pi_matrix));

alp_values = [1 0.99 0.9 0.8];

for a = 1:length(alp_values)
  alp = alp_values(a);

for i = 1:length(del_y_space)
  for j = 1:length(del_pi_space)
    del_pi_1 = del_pi_matrix(i,j);  
    del_y_1  = del_y_matrix(i,j);
    % matirces of normal form AEx(+1)=Bx
    A_22 = [(alp*sig-del_y_1)/sig (alp-del_pi_1)/sig; 0 bet*alp];
    B_21 = [-(((1-alp)*(thet^2)*sig-del_y_m1)/sig) -(((1-alp)*(thet^2)-del_pi_m1)/sig); 0 -bet*(1-alp)*(thet^2)];
    B_22 = [(sig+del_y_0)/sig del_pi_0/sig; -kap 1];
    A = [eye(2,2) zeros(2,2); zeros(2,2) A_22];
    B = [zeros(2,2) eye(2,2); B_21 B_22];
    %compute the eigenvalues of A^-1*B with eig()
    EV = eig(inv(A)*B);
    abs_value_EV = abs(EV);
    Number_of_unstable_EV = sum(abs_value_EV > 1);
      if Number_of_FLV==Number_of_unstable_EV
        Results(i,j) = 1; % determinacy = 1
      elseif Number_of_FLV < Number_of_unstable_EV
        Results(i,j) = -1; % instability = -1
      elseif (Number_of_FLV-Number_of_unstable_EV)==1
          Results(i,j) = 0.5; % order one indeterminacy = 0.5
      elseif (Number_of_FLV-Number_of_unstable_EV)==2
          Results(i,j) = 0; % oder two indeterminacy = 0
      else
          error('Something is wrong! Too many forward-looking variables!') 
      end  
  end
end
subplot(2,2,a)
contourf(del_pi_matrix,del_y_matrix,Results,2);
colormap(gray);
xlabel('\phi_\pi')
ylabel('\phi_y')
title(['\alpha = ' num2str(alp)]);
% new matrix to save results 
Results2(a,:) = [sum(sum(Results==1)) sum(sum(Results==-1)) sum(sum(Results==0.5)) sum(sum(Results==0))];
end
disp('Results based on eig() command:')
disp(['For alpha = 1: (det,instab,Or1ind,Or2ind) = ' num2str(Results2(1,:))])
disp(['For alpha = 0.99: (det,instab,Or1ind,Or2ind) = ' num2str(Results2(2,:))])
disp(['For alpha = 0.9: (det,instab,Or1ind,Or2ind) = ' num2str(Results2(3,:))])
disp(['For alpha = 0.8: (det,instab,Or1ind,Or2ind) = ' num2str(Results2(4,:))])

% Save the results: wheter we have determinacy = 1, indeterminacy = 0, insability = -1.
Results = NaN(size(del_pi_matrix));
figure
for a = 1:length(alp_values)
  alp = alp_values(a);

for i = 1:length(del_y_space)
  for j = 1:length(del_pi_space)
    del_pi_1 = del_pi_matrix(i,j);  
    del_y_1  = del_y_matrix(i,j);
    % matirces of normal form AEx(+1)=Bx
    A_22 = [(alp*sig-del_y_1)/sig (alp-del_pi_1)/sig; 0 bet*alp];
    B_21 = [-(((1-alp)*(thet^2)*sig-del_y_m1)/sig) -(((1-alp)*(thet^2)-del_pi_m1)/sig); 0 -bet*(1-alp)*(thet^2)];
    B_22 = [(sig+del_y_0)/sig del_pi_0/sig; -kap 1];
    A = [eye(2,2) zeros(2,2); zeros(2,2) A_22];
    B = [zeros(2,2) eye(2,2); B_21 B_22];
    %compute the eigenvalues with generalized Schur
    [S,T,Q,Z] = qz(A,B);
    t_ii = diag(T);
    s_ii = diag(S);
    abs_value_EV = abs(t_ii./s_ii);
    Number_of_unstable_EV = sum(abs_value_EV > 1);
      if Number_of_FLV==Number_of_unstable_EV
        Results(i,j) = 1; % determinacy = 1
      elseif Number_of_FLV < Number_of_unstable_EV
        Results(i,j) = -1; % instability = -1
      elseif (Number_of_FLV-Number_of_unstable_EV)==1
          Results(i,j) = 0.5; % order one indeterminacy = 0.5
      elseif (Number_of_FLV-Number_of_unstable_EV)==2
          Results(i,j) = 0; % oder two indeterminacy = 0
      else
          error('Something is wrong! Too many forward-looking variables!') 
      end  
  end
end
subplot(2,2,a)
contourf(del_pi_matrix,del_y_matrix,Results,2);
colormap(gray);
xlabel('\phi_\pi')
ylabel('\phi_y')
title(['\alpha = ' num2str(alp)]);
% new matrix to save results 
Results3(a,:) = [sum(sum(Results==1)) sum(sum(Results==-1)) sum(sum(Results==0.5)) sum(sum(Results==0))];
end
disp('')
disp('Results based on generalized Schur:')
disp(['For alpha = 1: (det,instab,Or1ind,Or2ind) = ' num2str(Results3(1,:))])
disp(['For alpha = 0.99: (det,instab,Or1ind,Or2ind) = ' num2str(Results3(2,:))])
disp(['For alpha = 0.9: (det,instab,Or1ind,Or2ind) = ' num2str(Results3(3,:))])
disp(['For alpha = 0.8: (det,instab,Or1ind,Or2ind) = ' num2str(Results3(4,:))])