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

% 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 = del_pi_matrix(i,j);  
    del_y  = del_y_matrix(i,j);
    % system matrices specifying the HE model, equation (23) in the paper.
    F = [1 0; -lamb 1];
    B = [alp-del_y/sig (alp-del_pi)/sig; 0 bet*alp];
    C = [(1-alp)*thet^2 ((1-alp)*thet^2)/sig; 0 bet*(1-alp)*sig^2];
    % matirces M, eq. 24 in the paper 
    M = [inv(B)*F -inv(B)*C; eye(2,2) zeros(2,2)];
    lambda = eig(M);
    abs_value_EV = abs(lambda);
    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!')     
      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 of the next round
%Results = NaN(size(del_pi_matrix));
end