
% Parameters
betta=0.98;
% Model 
No_controls=1;
No_states=1;
f_ps=linspace(0,3,100);
gamas=linspace(0,0.1,100);


% Construct matrices to store the feedback values
f_p_unique = [];
gama_unique = [];

f_p_unstable = [];
gama_unstable = [];

f_p_indet = [];
gama_indet = [];

tic;
%parfor jj = 1:length(f_ps)
for jj = 1:length(f_ps)
    f_p = f_ps(jj);
    
%    parfor kk = 1:length(gamas)
for kk = 1:length(gamas)
        gama = gamas(kk);
        [jj kk]
        % A Matrix
        A = [f_p * betta, 0, 0; 
            -(1/betta), (1/betta - gama), f_p; 
            0, 0, 1];
        [V D]=eig(A);
        roots=diag(D);

        J=D;
        P=V;
% 
%         %Check Jordan decomposition
% 
%         Check_J=norm(A - P*J*P^(-1));

%         if max(Check_J)<10^(-10)
%             display('Jordan decomposition checked')
% 
%         else
%             display('Jordan decomposition accuracy failed')
%         end


        % % Apply the Jordan Decomposition and check it
%         [P, J] = jordan(A); 
%         
%         check = A - P * J * inv(P); 
%         if norm(check) > 1e-6
%             warning('Jordan decomposition accuracy issue');
%         end
% 
%         roots = diag(J); 

        % Count unstable and stable roots
        No_unstable_roots = sum(abs(roots) > 1);
        No_stable_roots = sum(abs(roots) < 1);
   
        if No_controls == No_unstable_roots
             f_p_unique = [ f_p_unique; f_p];
            gama_unique = [gama_unique; gama];
            %disp('BK conditions are satisfied');
        elseif No_controls > No_unstable_roots
            %disp('No stable solution');
            f_p_unstable = [ f_p_unstable; f_p];
            gama_unstable = [gama_unstable; gama];
        elseif No_controls < No_unstable_roots
            %disp('Infinite solutions');
            f_p_indet = [ f_p_indet; f_p];
            gama_indet = [gama_indet; gama];
        end
    end
end

time_elapsed=toc
%%
scatter(f_p_unique,gama_unique,'b', 'filled')
hold on
scatter(f_p_unstable,gama_unstable,'g', 'filled')
hold on
scatter(f_p_indet,gama_indet,'r', 'filled')
title('Parameter space','Interpreter', 'latex','FontSize',20)
ylabel('Fiscal Policy, $\gamma$', 'Interpreter', 'latex','FontSize',20);
xlabel('Monetary Policy, $\phi_{\pi}$', 'Interpreter', 'latex','FontSize',20);
legend('Uniqueness','Instability','Indeterminacy','Interpreter', 'latex','FontSize',12)

saveas(gcf, 'BK_Leeper.pdf');
