% Code for computing and plotting indeterminacy plot for baseline model

dynare Basis_model_sim_bagud % Run the model one first time

lbound=0;
ubound=1;
theta_vec=linspace(lbound,ubound,50);
omega_vec=linspace(0,1,50);
[theta_mat,omega_mat]=meshgrid(theta_vec,omega_vec);
Z_plot_total=zeros(size(theta_mat));
info_mat=NaN(size(theta_mat));

for j=1:length(theta_vec)
    for k=1:length(omega_vec)
        set_param_value('theta',theta_mat(j,k));
        set_param_value('omega',omega_mat(j,k));
        [dr,info]=resol(0,M_,options_,oo_);
        if info(1)==0
            stoch_simul pihat;
            if pihat_epsilon_R(1) > -0.4    % Check sign of cumulative response
                info_mat(j,k)=1;
            else info_mat(j,k)=2;
            end
        else info_mat(j,k)=3;
        end
    end
end

%%

figure(1);
contourf(theta_mat,omega_mat,info_mat,2); 
%colormap(gray);
mymap=[0.7 0.7 0.7
       1 1 1  ];
colormap(mymap);

%hold on
%plot(1.122,4.925,'o','MarkerSize',10,'MarkerEdgeColor','Black','MarkerFaceColor','Black')
xlabel('\theta')
ylabel('\omega')
axis([lbound ubound 0 1]);

%psname  = 'tauplot_cumul_boot.ps';
%print('-depsc2',psname); 

