%{

%%% General presentation %%%

Author: Olivier Loisel (CREST).
Date: October 2024.
Description: This program produces the figures that are not about MMB models in the paper "New Principles For Stabilization
             Policy" (O. Loisel, October 2024): Figures 1-6 (in the main text) and Figures A.2-A.3 (in the Online Appendix).
Compatibility: The program works with Matlab R2023b. It requires Matlab's Optimization Toolbox (in order to use "fminbnd").
Running time: approximately 45 seconds with my computer.

%}

% CALIBRATED MODELS (indexed by mm):
% 1. Basic New Keynesian model with the calibration of Woodford (2003, Chapter 4) (if M1=1) or Gali (2015, Chapter 3) (if M1=2)
% 2. McKay, Nakamura, and Steinsson (2017)
% 3. Gabaix (2020)
% 4. Bilbiie (2008)
% 5. Ball (1999) and Svensson (1997) with Ball's (1999) calibration

M1 = 1; % please change to M1=2 to consider Gali's calibration of the basic New Keynesian model, instead of Woodford's calibration

% RULES (indexed by rr):
% 1. i_t = phi*E_t{pi_{t+h}} (Rule 1 in the paper)
% 2. i_t = phi*E_t{y_{t+h}} (Rule (9) in the paper)
% 3. i_t = rho*i_{t-1} + (1-rho)*phi*E_t{pi_{t+h}} (Rule 2 in the paper)

% SIGN OF PHI IN THE RULE (indexed by ss):
% 1. Positive: phi>0
% 2. Negative: phi<0

% General settings for optimization
options = optimset('TolX',eps);
tolerance = 1.e-12;

% General settings for the look of the figures
pix = 100; % number of pixels per horizon
nphi = 1000; % number of values for phi
band = 3.5/16; % for the rectangle defined by philowerbar<phi<phiupperbar and hmin<h<hmax to have a 16:9 aspect ratio
set(0,'DefaultAxesTickLabelInterpreter','latex');
set(0,'defaulttextInterpreter','latex');
set(0,'defaultAxesFontSize',6);
set(0,'DefaultAxesTitleFontWeight','normal');

% Setting the bounds of h in the figures
% The minimum hmin{ss}(mm,rr) should be lower than minus one
% The maximum hmax{ss}(mm,rr) should be higher than two
hmin{1} = [-10, -10;
           -10, -10;
           -10, -10;
           -10, -10;
           -8, -8];
hmin{2} = [-10, -10;
           -10, -10;
           -10, -10;
           -10, -10;
           -8, -8];
hmax{1} = [10, 10;
           10, 10;
           10, 10;
           10, 10;
           12, 12];
hmax{2} = [10, 10;
           10, 10;
           10, 10;
           10, 10;
           12, 12];
for ss=1:2
    nh{ss} = hmax{ss}-hmin{ss}+1; % number of horizons
end

% Initial NaN(mm,rr) settings for the bounds and thresholds of phi
phimin = NaN(5,4); % minimum of phi in the figures
phimax = NaN(5,4); % maximum of phi in the figures
philowerbar = NaN(5,4);
phiupperbar = NaN(5,4);
phi1 = NaN(5,4);
phiminus1 = NaN(5,4);


%%%%% MODELS 1 TO 4 %%%%%%
% IS equation: y_t = alpha*E_t{y_{t+1}} - (1/sigma)*(i_t-E_t{pi_{t+1}}), with 0<alpha<=1
% Phillips curve: pi_t = beta*E_t{pi_{t+1}} + kappa*y_t, with 0<beta<1 and kappa>0

% Calibrations:
alpha = [1; 0.97; 0.85; 1];
beta = [0.99; 0.99; 0.792; 0.99];
sigma = [0.16; 2.67; 5; -0.11];
kappa = [0.022; 0.02; 0.11; 0.228];
if M1==2
    sigma(1) = 1;
    kappa(1) = 0.125;
end

% Computation of the phi thresholds using Matlab's Optimization Toolbox
for mm=1:4
    for rr=[1 2]
        Q = [1, -(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm)), alpha(mm)*beta(mm)];
        if (rr==1)
            R = [kappa(mm)/sigma(mm)];
        else
            R = [1/sigma(mm), -beta(mm)/sigma(mm)];
        end
        f = @(x) abs(polyval(Q,complex(x,sqrt(1-x^2)))/polyval(R,complex(x,sqrt(1-x^2))));
        [~,mini] = fminbnd(f,-1,1,options);
        philowerbar(mm,rr) = min([mini f(1) f(-1)]);
        nf = @(x) -f(x);
        [~,mini] = fminbnd(nf,-1,1,options);
        phiupperbar(mm,rr) = -min([mini nf(1) nf(-1)]);
        phi1(mm,rr) = -polyval(Q,1)/polyval(R,1);
        phiminus1(mm,rr) = abs(polyval(Q,-1)/polyval(R,-1));
    end
end

% Check of the phi thresholds using their analytical expression
phiupperbarpia = NaN(4,1);
for mm=1:4
    if (sigma(mm)>0)
        phiupperbarpia(mm) = 1+(1+alpha(mm))*(1+beta(mm))*sigma(mm)/kappa(mm); % phiupperbar for inflation
    else
        phiupperbarpia(mm) = max(1,2*(1+beta(mm))*abs(sigma(mm))/kappa(mm)-1);
    end
end
phiupperbarya = max(abs((1+alpha).*sigma+kappa./(1+beta)),abs((1-alpha).*sigma-kappa./(1-beta))); % phiupperbar for output
phiupperbara = [phiupperbarpia, phiupperbarya];
philowerbarpia = NaN(4,1);
for mm=1:4
    if (sigma(mm)>0)
        philowerbarpia(mm) = abs(1-(1-alpha(mm))*(1-beta(mm))*sigma(mm)/kappa(mm)); % philowerbar for inflation
    elseif (abs(1+beta(mm)+kappa(mm)/sigma(mm))<4*beta(mm)/(1+beta(mm)))
        philowerbarpia(mm) = (1-beta(mm))*abs(sigma(mm))/kappa(mm)*sqrt(1-(1+beta(mm)+kappa(mm)/sigma(mm))^2/(4*beta(mm)));
    else
        philowerbarpia(mm) = min(1,abs(1+2*(1+beta(mm))*sigma(mm)/kappa(mm)));
    end
end    
philowerbarya = NaN(4,1);
eta = beta.*((1-alpha.*beta).^2+(alpha+beta+kappa./sigma).^2) - (1+beta.^2).*(1+alpha.*beta).*(alpha+beta+kappa./sigma);
for mm=1:4
    if abs(eta(mm)-4*alpha(mm)*beta(mm)^2) < 4*alpha(mm)*beta(mm)*(1+beta(mm)^2)
        philowerbarya(mm) = abs(sigma(mm)) / sqrt(beta(mm)) * sqrt((1+alpha(mm)*beta(mm))*(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm)) ...
                           + 2*sqrt(alpha(mm)^2*(1+beta(mm)^2)^2+alpha(mm)*eta(mm)) - 2*alpha(mm)*(1+beta(mm)^2));
    else
        philowerbarya(mm) = min(abs((1+alpha(mm))*sigma(mm)+kappa(mm)/(1+beta(mm))),abs((1-alpha(mm))*sigma(mm)-kappa(mm)/(1-beta(mm))));
    end
end
philowerbara = [philowerbarpia, philowerbarya];
checkphilowerbar = philowerbara - philowerbar(1:4,1:2);
checkphiupperbar = phiupperbara - phiupperbar(1:4,1:2);
if any([(abs(checkphilowerbar)>tolerance),(abs(checkphiupperbar)>tolerance)],'all')
    disp('Mistake in the computation of some phi thresholds.');
end

% Setting the bounds of phi in the figures
phimin([1:4],[1 2]) = exp((band*log(phiupperbar([1:4],[1 2]))-(1-band)*log(philowerbar([1:4],[1 2])))/(2*band-1));
phimax([1:4],[1 2]) = exp((band*log(philowerbar([1:4],[1 2]))-(1-band)*log(phiupperbar([1:4],[1 2])))/(2*band-1));

% Determination of the determinacy status
for mm=1:4
    for ss=1:2
        for rr=1:2
            M=ones(nphi,nh{ss}(mm,rr));
            for hh=1:nh{ss}(mm,rr)
                for pp=1:nphi
                    P = zeros(1,nh{ss}(mm,rr)+rr-1);
                    phi = (3-2*ss)*exp(log(phimin(mm,rr))+(pp-1)*(log(phimax(mm,rr))-log(phimin(mm,rr)))/(nphi-1)); % log scale on y axis
                    if (rr==1)
                        P(1,hh) = kappa(mm)*phi/sigma(mm);
                    else
                        P(1,hh) = phi/sigma(mm);
                        P(1,hh+1) = -beta(mm)*phi/sigma(mm);
                    end
                    P = P + [zeros(1,-hmin{ss}(mm,rr)) 1 -(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm)) alpha(mm)*beta(mm) zeros(1,hmax{ss}(mm,rr)+rr-3)];
                    Roo = sort(abs(roots(P)));
                    Roo = Roo(Roo~=0);
                    nnpv = max(2,hh+hmin{ss}(mm,rr)+rr-2); % number of non-predetermined variables
                    if (Roo(nnpv)>1)
                        M((ss-1)*(nphi+1)+(3-2*ss)*pp,hh)=0.7; % In this case, there is multiplicity
                    elseif ((hh+hmin{ss}(mm,rr)<=0) & (Roo(nnpv+1)<1))
                        M((ss-1)*(nphi+1)+(3-2*ss)*pp,hh)=0.4; % In this case, there is explosiveness
                    end
                end
            end
            N{mm,rr,ss}=kron(M,ones(1,pix));
        end
    end
end


%%%%%% MODEL 5 %%%%%%
% IS equation: y_t = lambda*y_{t-1} - mu*(i_{t-1}-pi_{t-1}), with 0<lambda<1 and mu>0
% Phillips curve: pi_t = pi_{t-1} + chi*y_{t-1}, with chi>0
mm = 5;

% Calibration
chi = 0.4;
mu = 1;
lambda = 0.8;

% Computation of the phi thresholds using Matlab's optimization toolbox
for rr=1
    Q = [lambda-chi*mu -(1+lambda) 1];
    R = [chi*mu];
    f = @(x) abs(polyval(Q,complex(x,sqrt(1-x^2)))/polyval(R,complex(x,sqrt(1-x^2))));
    [~,mini] = fminbnd(f,-1,1);
    philowerbar(mm,rr) = min([mini f(1) f(-1)]);
    nf = @(x) -f(x);
    [~,mini] = fminbnd(nf,-1,1);
    phiupperbar(mm,rr) = -min([mini nf(1) nf(-1)]);
    phi1(mm,rr) = -polyval(Q,1)/polyval(R,1);
    phiminus1(mm,rr) = abs(polyval(Q,-1)/polyval(R,-1));
end
for rr=2
    Q = [lambda-chi*mu -(1+lambda) 1];
    R = [-mu mu];
    f = @(x) abs(polyval(Q,complex(x,sqrt(1-x^2)))/polyval(R,complex(x,sqrt(1-x^2))));
    [~,mini] = fminbnd(f,-1,1);
    philowerbar(mm,rr) = min([mini f(1) f(-1)]);
end
    
% Check of the phi thresholds using their analytical expression
philowerbara(mm,1) = min(1,abs(2*(1+lambda)/(chi*mu)-1));
if (lambda<chi*mu*(1+chi*mu/16))
    philowerbara(mm,2) = abs((1+lambda)/mu-chi/2);
else
    philowerbara(mm,2) = (1/mu)*sqrt((1-lambda)^2+chi*mu*(3-lambda+2*sqrt(lambda-chi*mu)));
end
if (chi*mu<(lambda^2+6*lambda+1)/(lambda+5))
    phiupperbara(mm,1) = 2*(1+lambda)/(chi*mu)-1;
else
    phiupperbara(mm,1) = (1+chi*mu-lambda)/(chi*mu)*sqrt(1+(1+lambda)^2/(4*(chi*mu-lambda)));
end
checkphilowerbar = philowerbara(mm,1:2) - philowerbar(mm,1:2);
checkphiupperbar = phiupperbara(mm,1) - phiupperbar(mm,1);
if any([(abs(checkphilowerbar)>tolerance),(abs(checkphiupperbar)>tolerance)],'all')
    disp('Mistake in the computation of some phi thresholds.');
end

% Setting the bounds of phi in the figures
phimin(mm,1) = exp((band*log(phiupperbar(mm,1))-(1-band)*log(philowerbar(mm,1)))/(2*band-1));
phimax(mm,1) = exp((band*log(philowerbar(mm,1))-(1-band)*log(phiupperbar(mm,1)))/(2*band-1));
phimax(mm,2) = 2;
phimin(mm,2) = exp((log(philowerbar(mm,2))-band*log(phimax(mm,2)))/(1-band));

% Determination of the determinacy status
for ss=1:2
    for rr=1:2
        M=ones(nphi,nh{ss}(mm,rr));
        for hh=1:nh{ss}(mm,rr)
            for pp=1:nphi
                P = zeros(1,nh{ss}(mm,rr)+rr-1);
                phi = (3-2*ss)*exp(log(phimin(mm,rr))+(pp-1)*(log(phimax(mm,rr))-log(phimin(mm,rr)))/(nphi-1)); % log scale on y axis
                if (rr==1)
                    P(1,hh) = chi*mu*phi;
                else
                    P(1,hh) = -mu*phi;
                    P(1,hh+1) = mu*phi;
                end
                P = P + [zeros(1,-hmin{ss}(mm,rr)) lambda-chi*mu -(1+lambda) 1 zeros(1,hmax{ss}(mm,rr)+rr-3)];
                Roo = sort(abs(roots(P)));
                Roo = Roo(Roo~=0);
                nnpv = max(0,hh+hmin{ss}(mm,rr)+rr-4); % number of non-predetermined variables
                if (Roo(nnpv+1)<1)
                    M((ss-1)*(nphi+1)+(3-2*ss)*pp,hh)=0.4; % In this case, there is explosiveness
                elseif ((nnpv>=1) & (Roo(nnpv)>1))
                    M((ss-1)*(nphi+1)+(3-2*ss)*pp,hh)=0.7; % In this case, there is multiplicity
                end
            end
        end
        N{mm,rr,ss}=kron(M,ones(1,pix));
    end
end


%%%%%% FIGURES 2, 5 AND 6 %%%%%%

width = 11;
height = 4.33;
height2 = 13;
height3 = 8.66;

f=figure('Name','Figure 2','Units','centimeters');
set(gcf,'PaperUnits','centimeters','Position',[0 0 width height],'PaperSize',[width height],'PaperPosition',[0 0 width height]);
ax = gca;
ax.TitleFontSizeMultiplier = 1;
mm=1;
rr=1;
for ss=1:2
    subplot(1,2,ss);
    contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
    colormap(gray);
    if (ss==1)
        title('(a) $\phi>0$','FontSize',6);
    else
        title('(b) $\phi<0$','FontSize',6);
    end
    caxis([-1, 1]);
    xlabel('$h$','FontSize',6);
    if (ss==1)
        ylabel('$\phi$ (log scale)','FontSize',6);
        xline(-hmin{ss}(mm,rr)*pix,':k');
        xline((-hmin{ss}(mm,rr)+1)*pix,':k');
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$0$';'$\lfloor h_1\rfloor =1$';hmax{ss}(mm,rr)},'XTick',...
            [round(pix/2), -hmin{ss}(mm,rr)*pix+round(pix/4), (-hmin{ss}(mm,rr)+1)*pix+round(3*pix/4), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
        ylowerbar = round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))));
        line([(-hmin{ss}(mm,rr)+2)*pix (-hmin{ss}(mm,rr)+2)*pix],[0 ylowerbar],'LineStyle',':','Color',[.25 .25 .25]);
        set(gca,'YTickLabel',{['$\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
            ['$\bar{\phi}=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]},'YTick',...
            [round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    else
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$0$';hmax{ss}(mm,rr)},'XTick',...
            [round(pix/2), -hmin{ss}(mm,rr)*pix+round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
        set(gca,'YTickLabel',{['-$\bar{\phi}=$ ' num2str(-round(phiupperbar(mm,rr),2,'significant'))];...
            ['-$\underline{\phi}=$ ' num2str(-round(philowerbar(mm,rr),2,'significant'))]},'YTick', ...
            [nphi+1-round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            nphi+1-round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    end
    yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/...
        (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
    yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/...
        (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
    axis square;
end
print(f,'-r1000','-dpdf','Fig2');



f=figure('Name','Figure 5','Units','centimeters');
set(gcf,'PaperUnits','centimeters','Position',[0 0 width height3],'PaperSize',[width height3],'PaperPosition',[0 0 width height3]);
ax = gca;
ax.TitleFontSizeMultiplier = 1;
rr=1;
ss=1;
for mm=2:5
    subplot(2,2,mm-1);
    contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
    colormap(gray);
    caxis([-1, 1]);
    yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/...
        (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
    yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/...
        (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
    if mm~=5
        xline(-hmin{ss}(mm,rr)*pix,':k');
        xline((-hmin{ss}(mm,rr)+1)*pix,':k');
    end
    if (mm==2)
        title('(a) Model 2','FontSize',6);
        ylabel('$\phi$ (log scale)','FontSize',6);
        ylowerbar = round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))));
        line([(-hmin{ss}(mm,rr)+6)*pix (-hmin{ss}(mm,rr)+6)*pix],[0 ylowerbar],'LineStyle',':','Color',[.25 .25 .25]);
        line([(-hmin{ss}(mm,rr)+7)*pix (-hmin{ss}(mm,rr)+7)*pix],[0 ylowerbar],'LineStyle',':','Color',[.25 .25 .25]);
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$h^\ast=0$';'$\lfloor h_1\rfloor =6$';hmax{ss}(mm,rr)},'XTick',...
            [round(pix/2), -hmin{ss}(mm,rr)*pix+round(pix/2), (-hmin{ss}(mm,rr)+6)*pix+round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
        set(gca,'YTickLabel',{['$\phi_1=\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
            ['$\phi_{-1}=\bar{\phi}=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]},'YTick',...
            [round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    elseif (mm==3)
        title('(b) Model 3','FontSize',6);
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$h^\ast=0$';hmax{ss}(mm,rr)},'XTick',...
            [round(pix/2), -hmin{ss}(mm,rr)*pix+round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
        xtickangle(30);
        set(gca,'YTickLabel',{['$\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
            ['$\phi_{-1}=\bar{\phi}=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]},'YTick',...
            [round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    elseif (mm==4)
        title('(c) Model 4','FontSize',6);
        xlabel('$h$','FontSize',6);
        ylabel('$\phi$ (log scale)','FontSize',6);
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$h^\ast=0$';hmax{ss}(mm,rr)},'XTick',...
            [round(pix/2), -hmin{ss}(mm,rr)*pix+round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
        xtickangle(30);
        yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phiminus1(mm,rr))-log(phimin(mm,rr)))/...
            (log(phimax(mm,rr))-log(phimin(mm,rr)))),':k','LineWidth',1);
        annotation(gcf,'arrow', [0.156 0.170],[0.344 0.372],'LineWidth',0.25,'HeadStyle','vback3','HeadLength',3,'HeadWidth',3);
        set(gca,'YTickLabel',{['$\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
            ['$\phi_{-1}=$ ' num2str(round(phiminus1(mm,rr),2,'significant')) ' '];...
            ['$\phi_1=\bar{\phi}=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]},'YTick',...
            [round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            round(1+(nphi-1)*(log(0.3*phiminus1(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    elseif (mm==5)
        title('(d) Model 5','FontSize',6);
        xlabel('$h$','FontSize',6);
        xline((-hmin{ss}(mm,rr)+2)*pix,':k');
        xline((-hmin{ss}(mm,rr)+3)*pix,':k');
        ylowerbar = round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))));
        line([(-hmin{ss}(mm,rr))*pix (-hmin{ss}(mm,rr))*pix],[0 ylowerbar],'LineStyle',':','Color',[.25 .25 .25]);
        line([(-hmin{ss}(mm,rr)+1)*pix (-hmin{ss}(mm,rr)+1)*pix],[0 ylowerbar],'LineStyle',':','Color',[.25 .25 .25]);
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$\lceil h_1\rceil =0$';'$h^\ast=2$';hmax{ss}(mm,rr)},'XTick',...
            [round(pix/2), (-hmin{ss}(mm,rr))*pix+round(pix/4), (2-hmin{ss}(mm,rr))*pix+round(3*pix/4), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
        set(gca,'YTickLabel',{['$\phi_1=\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
            ['$\phi_{-1}=\bar{\phi}=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]},'YTick',...
            [round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    end
end
axesHandles = findobj(get(gcf,'Children'),'flat','Type','axes');
axis(axesHandles,'square')
print(f,'-r1000','-dpdf','Fig5');



f=figure('Name','Figure 6','Units','centimeters');
set(gcf,'PaperUnits','centimeters','Position',[0 0 width height],'PaperSize',[width height],'PaperPosition',[0 0 width height]);
ax = gca;
ax.TitleFontSizeMultiplier = 1;
rr=2;
ss=1;
for mm=1:2
    subplot(1,2,mm);
    contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
    colormap(gray);
    if (mm==1)
        title('(a) Model 1','FontSize',6);
    else
        title('(b) Model 2','FontSize',6);
    end
    caxis([-1, 1]);
    xlabel('$h$','FontSize',6);
    set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$h^\ast=0$';hmax{ss}(mm,rr)},'XTick',...
        [round(pix/2), -hmin{ss}(mm,rr)*pix+round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
    xtickangle(30);
    xline(-hmin{ss}(mm,rr)*pix,':k');
    xline((-hmin{ss}(mm,rr)+1)*pix,':k');
    if (mm==1)
        ylabel('$\phi$ (log scale)','FontSize',6);
        yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phiminus1(mm,rr))-log(phimin(mm,rr)))/...
            (log(phimax(mm,rr))-log(phimin(mm,rr)))),':k','LineWidth',1);
        set(gca,'YTickLabel',{['$\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
            ['$\phi_{-1}=$ ' num2str(round(phiminus1(mm,rr),2,'significant'))];...
            ['$\phi_1=\bar{\phi}=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]},'YTick',...
            [round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            round(1+(nphi-1)*(log(phiminus1(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    else
        yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phi1(mm,rr))-log(phimin(mm,rr)))/...
            (log(phimax(mm,rr))-log(phimin(mm,rr)))),':k','LineWidth',1);
        set(gca,'YTickLabel',{['$\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
            ['$\phi_1=$ ' num2str(round(phi1(mm,rr),2,'significant'))];...
            ['$\phi_{-1}=\bar{\phi}=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]},'YTick',...
            [round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            round(1+(nphi-1)*(log(phi1(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    end
    yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/...
        (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
    yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/...
        (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
    axis square;
end
print(f,'-r1000','-dpdf','Fig6');



%%%%%% FIGURE 3 %%%%%%

mm = 1;
rr = 3;
ss = 1;
rho = 0.8;

Q = [-rho, 1+rho*(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm)), -(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm))-rho*alpha(mm)*beta(mm), alpha(mm)*beta(mm)]/(1-rho);
R = [kappa(mm)/sigma(mm)];
f = @(x) abs(polyval(Q,complex(x,sqrt(1-x^2)))/polyval(R,complex(x,sqrt(1-x^2))));
[~,mini] = fminbnd(f,-1,1,options);
philowerbar(mm,rr) = min([mini f(1) f(-1)]);
nf = @(x) -f(x);
[~,mini] = fminbnd(nf,-1,1,options);
phiupperbar(mm,rr) = -min([mini nf(1) nf(-1)]);
phi1(mm,rr) = -polyval(Q,1)/polyval(R,1);
phiminus1(mm,rr) = abs(polyval(Q,-1)/polyval(R,-1));

phimin(1,3) = exp((band*log(phiupperbar(1,3))-(1-band)*log(philowerbar(1,3)))/(2*band-1));
phimax(1,3) = exp((band*log(philowerbar(1,3))-(1-band)*log(phiupperbar(1,3)))/(2*band-1));

M=0.999*ones(nphi,21);
for hh=1:21
    for pp=1:nphi
        P = zeros(1,21);
        phi = (3-2*ss)*exp(log(phimin(mm,rr))+(pp-1)*(log(phimax(mm,rr))-log(phimin(mm,rr)))/(nphi-1)); % log scale on y axis
        P(1,hh) = kappa(mm)*phi/sigma(mm);
        P = P + [zeros(1,9) -rho 1+rho*(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm)) ...
            -(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm))-rho*alpha(mm)*beta(mm) alpha(mm)*beta(mm) zeros(1,8)]/(1-rho);
        Roo = sort(abs(roots(P)));
        Roo = Roo(Roo~=0);
        nnpv = max(2,hh-11); % number of non-predetermined variables
        if (Roo(nnpv)>1)
            M((ss-1)*(nphi+1)+(3-2*ss)*pp,hh)=0.7; % In this case, there is multiplicity
        elseif ((hh-10<=0) & (Roo(nnpv+1)<1))
            M((ss-1)*(nphi+1)+(3-2*ss)*pp,hh)=0.4; % In this case, there is explosiveness
        end
    end
end
      
for hh=1:21
    for pp=1:nphi
        P = zeros(1,21);
        phi = exp(log(phimin(mm,rr))+(pp-1)*(log(phimax(mm,rr))-log(phimin(mm,rr)))/(nphi-1)); % log scale on y axis
        P(1,hh) = kappa(mm)*phi/sigma(mm);
        P = P + [zeros(1,10) 1 -(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm)) alpha(mm)*beta(mm) zeros(1,8)];
        Roo = sort(abs(roots(P)));
        Roo = Roo(Roo~=0);
        nnpv = max(2,hh-11);
        if (Roo(nnpv)>1)
        elseif ((hh-10<=0) & (Roo(nnpv+1)<1))
        else
            M(pp,hh)=1; % This will serve to draw the boundaries of the D region under Rule 1
        end
    end
end

N{mm,rr,ss}=kron(M,ones(1,pix));

f=figure('Name','Figure 3','Units','centimeters');
set(gcf,'PaperUnits','centimeters','Position',[0 0 width height],'PaperSize',[width height],'PaperPosition',[0 0 width height]);
ax = gca;
ax.TitleFontSizeMultiplier = 1;
contourf(N{mm,rr,1},[0.4, 0.7, 0.999],'-k');
colormap(gray);
hold on
contourf(N{mm,rr,1},[1, 1],'-.k');
colormap(gray);
hold off
axis equal;
caxis([-1, 1]);
xlabel('$h$','FontSize',6);
set(gca,'Ticklength', [0 0],'XTickLabel',{-10;'0';'$\lfloor h_1\rfloor =5$';10},'XTick',...
    [round(pix/2), 10*pix+round(pix/2), 15*pix+round(pix/2), 20*pix+round(pix/2)],'FontSize',6);
ylabel('$\phi$ (log scale)','FontSize',6);
yupperbar = round(1+(nphi-1)*(log(phiupperbar(mm,1))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))));
line([10*pix 10*pix],[0 yupperbar],'LineStyle',':','Color',[.25 .25 .25]);
line([11*pix 11*pix],[0 yupperbar],'LineStyle',':','Color',[.25 .25 .25]);
ylowerbar = round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))));
line([15*pix 15*pix],[0 ylowerbar],'LineStyle',':','Color',[.25 .25 .25]);
line([16*pix 16*pix],[0 ylowerbar],'LineStyle',':','Color',[.25 .25 .25]);
set(gca,'YTickLabel',{['$\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
    ['$\bar{\phi}=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]},'YTick',...
    [round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
    round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/...
    (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/...
    (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
axis square;
print(f,'-r1000','-dpdf','Fig3');



%%%%%% FIGURES A.2 AND A.3 %%%%%%

f=figure('Name','Figure A.2','Units','centimeters');
set(gcf,'PaperUnits','centimeters','Position',[0 0 width height3],'PaperSize',[width height3],'PaperPosition',[0 0 width height3]);
ax = gca;
ax.TitleFontSizeMultiplier = 1;
rr=1;
ss=2;
for mm=2:5
    subplot(2,2,mm-1);
    contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
    colormap(gray);
    caxis([-1, 1]);
    yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/...
        (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
    yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/...
        (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
    if mm~=5
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$h^\ast=0$';hmax{ss}(mm,rr)},'XTick',...
            [round(pix/2), -hmin{ss}(mm,rr)*pix+round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
        xtickangle(30);
    end
    if (mm==2)
        title('(a) Model 2 ','FontSize',6);
        ylabel('$\phi$ (log scale)','FontSize',6,'Position',[-12*pix 5*pix]);
        set(gca,'YTickLabel',{['-$\phi_{-1}=$ -$\bar{\phi}=$ ' num2str(-round(phiupperbar(mm,rr),2,'significant'))];...
            ['-$\underline{\phi}=$ ' num2str(-round(philowerbar(mm,rr),2,'significant'))]},'YTick', ...
            [nphi+1-round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            nphi+1-round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    elseif (mm==3)
        title('(b) Model 3 ','FontSize',6);
        set(gca,'YTickLabel',{['-$\phi_{-1}=$ -$\bar{\phi}=$ ' num2str(-round(phiupperbar(mm,rr),2,'significant'))];...
            ['$\phi_1=$ -$\underline{\phi}=$ ' num2str(-round(philowerbar(mm,rr),2,'significant'))]},'YTick', ...
            [nphi+1-round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            nphi+1-round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    elseif (mm==4)
        title('(c) Model 4 ','FontSize',6);
        xlabel('$h$','FontSize',6);
        ylabel('$\phi$ (log scale)','FontSize',6);
        yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phiminus1(mm,rr))-log(phimin(mm,rr)))/...
            (log(phimax(mm,rr))-log(phimin(mm,rr)))),':k','LineWidth',1);
        annotation(gcf,'arrow', [0.163 0.177],[0.222 0.191],'LineWidth',0.25,'HeadStyle','vback3','HeadLength',3,'HeadWidth',3);
        set(gca,'YTickLabel',{['-$\bar{\phi}=$ ' num2str(-round(phiupperbar(mm,rr),2,'significant'))];...
            ['-$\phi_{-1}=$ ' num2str(-round(phiminus1(mm,rr),2,'significant')) ' '];...
            ['-$\underline{\phi}=$ ' num2str(-round(philowerbar(mm,rr),2,'significant'))]},'YTick', ...
            [nphi+1-round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            nphi+1-round(1+(nphi-1)*(log(0.3*phiminus1(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            nphi+1-round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    elseif (mm==5)
        title('(d) Model 5 ','FontSize',6);
        xlabel('$h$','FontSize',6);
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$h^\ast=2$';hmax{ss}(mm,rr)},'XTick',...
            [round(pix/2), (2-hmin{ss}(mm,rr))*pix+round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
        xtickangle(30);
        set(gca,'YTickLabel',{['-$\phi_{-1}=$ -$\bar{\phi}=$ ' num2str(-round(phiupperbar(mm,rr),2,'significant'))];...
            ['-$\underline{\phi}=$ ' num2str(-round(philowerbar(mm,rr),2,'significant'))]},'YTick', ...
            [nphi+1-round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            nphi+1-round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    end
end
axesHandles = findobj(get(gcf,'Children'),'flat','Type','axes');
axis(axesHandles,'square')
print(f,'-r1000','-dpdf','FigA2');


f=figure('Name','Figure A.3','Units','centimeters');
set(gcf,'PaperUnits','centimeters','Position',[0 0 width height],'PaperSize',[width height],'PaperPosition',[0 0 width height]);
ax = gca;
ax.TitleFontSizeMultiplier = 1;
rr=2;
ss=2;
for mm=1:2
    subplot(1,2,mm);
    contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
    colormap(gray);
    if (mm==1)
        title('(a) Model 1','FontSize',6);
    else
        title('(b) Model 2','FontSize',6);
    end
    caxis([-1, 1]);
    xlabel('$h$','FontSize',6);
    set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$h^\ast=0$';hmax{ss}(mm,rr)},'XTick',...
        [round(pix/2), -hmin{ss}(mm,rr)*pix+round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
    xtickangle(30);
    if (mm==1)
        ylabel('$\phi$ (log scale)','FontSize',6);
        yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phiminus1(mm,rr))-log(phimin(mm,rr)))/...
            (log(phimax(mm,rr))-log(phimin(mm,rr)))),':k','LineWidth',1);
        set(gca,'YTickLabel',{['-$\bar{\phi}=$ ' num2str(-round(phiupperbar(mm,rr),2,'significant'))];...
            ['-$\phi_{-1}=$ ' num2str(-round(phiminus1(mm,rr),2,'significant'))];...
            ['-$\underline{\phi}=$ ' num2str(-round(philowerbar(mm,rr),2,'significant'))]},'YTick', ...
            [nphi+1-round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            nphi+1-round(1+(nphi-1)*(log(phiminus1(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            nphi+1-round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    else
        set(gca,'YTickLabel',{['-$\phi_{-1}=$ -$\bar{\phi}=$ ' num2str(-round(phiupperbar(mm,rr),2,'significant'))];...
            ['-$\underline{\phi}=$ ' num2str(-round(philowerbar(mm,rr),2,'significant'))]},'YTick', ...
            [nphi+1-round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))),...
            nphi+1-round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr))))],'FontSize',6);
    end
    yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(philowerbar(mm,rr))-log(phimin(mm,rr)))/...
        (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
    yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phiupperbar(mm,rr))-log(phimin(mm,rr)))/...
        (log(phimax(mm,rr))-log(phimin(mm,rr)))),'--k','LineWidth',1);
    axis square;
end
print(f,'-r1000','-dpdf','FigA3');


%%%%%% FIGURES 1 AND 4 %%%%%%

hmint = -10;
hmaxt = 10;
nht = hmaxt-hmint+1;
factor1 = 2.471;
factor2 = 2.62;

for tt=1:2
    Nv1{tt} = ones(nphi,nht);
    Nv1{tt}(1:round(nphi*band),:) = 1-0.3*tt;
    Nv1{tt}(round(nphi*band)+5:round(nphi*(1-band))-1,-19+20*tt) = 0.1+0.3*tt;
    Nv1{tt}(round(nphi*band)+round(factor1)+5:round(nphi*(1-band))-1,-16+18*tt) = 0.1+0.3*tt;
    Nv1{tt}(round(nphi*band)+round(factor1^2)+5:round(nphi*(1-band))-1,-13+16*tt) = 0.1+0.3*tt;
    Nv1{tt}(round(nphi*band)+round(factor1^3)+5:round(nphi*(1-band))-1,-10+14*tt) = 0.1+0.3*tt;
    Nv1{tt}(round(nphi*band)+round(factor1^4)+5:round(nphi*(1-band))-1,-7+12*tt) = 0.1+0.3*tt;
    Nv1{tt}(round(nphi*band)+round(factor1^5)+5:round(nphi*(1-band))-1,-4+10*tt) = 0.1+0.3*tt;
    Nv1{tt}(round(nphi*band)+round(factor1^6)+5:round(nphi*(1-band))-1,-1+8*tt) = 0.1+0.3*tt;
    Nv1{tt}(round(nphi*band)+1:round(nphi*(1-band))-1,29-14*tt:35-14*tt) = 1-0.3*tt;
    Nv1{tt}(round(nphi*(1-band)):nphi,1:10) = 0.4;
    Nv1{tt}(round(nphi*(1-band)):nphi,12:nht) = 0.7;
    Nv1{tt} = kron(Nv1{tt},ones(1,pix));
    for uu=1:2
        for ww=1
            Nv1{tt}(round(nphi*band)+5*ww,(-3+5*tt)*pix+mod(ww+uu-1,10)*10:100:(9+5*tt)*pix-10) = 0.1;
        end
        for ww=2:3
            Nv1{tt}(round(nphi*band)+5*ww,(-1+4*tt)*pix+mod(ww+uu-1,10)*10:100:(10+4*tt)*pix-10) = 0.1;
        end
        for ww=4:7
            Nv1{tt}(round(nphi*band)+5*ww,(1+3*tt)*pix+mod(ww+uu-1,10)*10:100:(11+3*tt)*pix-10) = 0.1;
        end
        for ww=8:19
            Nv1{tt}(round(nphi*band)+5*ww,(3+2*tt)*pix+mod(ww+uu-1,10)*10:100:(12+2*tt)*pix-10) = 0.1;
        end
        for ww=20:46
            Nv1{tt}(round(nphi*band)+5*ww,(5+1*tt)*pix+mod(ww+uu-1,10)*10:100:(13+1*tt)*pix-10) = 0.1;
        end
        for ww=47:round(nphi*(1-2*band)/5)-1
            Nv1{tt}(round(nphi*band)+5*ww,7*pix+mod(ww+uu-1,10)*10:100:(14+0*tt)*pix-10) = 0.1;
        end
    end
end
Nv1{3}(1:round(nphi*band),:) = 1;
Nv1{3}(round(nphi*band)+1:nphi,1:10.5*pix) = Nv1{1}(round(nphi*band)+1:nphi,1:10.5*pix);  
Nv1{3}(round(nphi*band)+1:nphi,10.5*pix+1:21*pix) = Nv1{2}(round(nphi*band)+1:nphi,10.5*pix+1:21*pix); 


f=figure('Name','Figure 1','Units','centimeters');
set(gcf,'PaperUnits','centimeters','Position',[0 0 width height],'PaperSize',[width height],'PaperPosition',[0 0 width height]);
contourf(Nv1{1},[0, 0.4, 0.7, 1]);
colormap(gray);
caxis([-1, 1]);
xlabel('$h$','FontSize',6);
set(gca,'Ticklength', [0 0],'XTickLabel',{'$-\infty$';'$0$';'$+\infty$'},'XTick',...
    [round(pix/2), round(nht*pix/2), (nht-1)*pix+round(pix/2)],'FontSize',6);
ylabel('$\left|\phi\right|$','FontSize',6,'Position',[-8*pix 5*pix]);
xline(round((nht-1)*pix/2),':k','LineWidth',0.5);
xline(round((nht+1)*pix/2),':k','LineWidth',0.5);
yline(round(nphi*band),'--k','LineWidth',0.75);
yline(round(nphi*(1-band)),'--k','LineWidth',0.75);
set(gca,'YTickLabel',{'$0$'; '$\underline{\phi}=1$'; '$\bar{\phi}=1+2(1+\beta)\frac{\sigma}{\kappa}$'; '$+\infty$'},'YTick',...
    [1, round(nphi*band), round(nphi*(1-band)), nphi],'FontSize',6);
axis square;
print(f,'-r1000','-dpdf','Fig1');


f=figure('Name','Figure 4','Units','centimeters');
set(gcf,'PaperUnits','centimeters','Position',[0 0 width height],'PaperSize',[width height],'PaperPosition',[0 0 width height]);
for tt=1:3
    subplot(1,3,tt);
    colormap(gray);
    if (tt==1)
        contourf(Nv1{tt},[0, 0.4, 0.7, 1]);
        title('(a) Models with $S_{peg}=M$','FontSize',5);
    elseif (tt==2)
        contourf(Nv1{3},[0, 0.4, 0.7, 1]);
        title('(b) Models with $S_{peg}=D$','FontSize',5);
    else
        contourf(Nv1{2},[0, 0.4, 0.7, 1]);
        title('(c) Models with $S_{peg}=E$','FontSize',5);
    end
    caxis([-1, 1]);
    xlabel('$h$','FontSize',5);
    set(gca,'Ticklength', [0 0],'XTickLabel',{'$-\infty$';'$h^\ast$';'$+\infty$'},'XTick',...
        [round(pix/2), round(nht*pix/2), (nht-1)*pix+round(pix/2)],'FontSize',5);
    if (tt==1)
        ylabel('$\left|\phi\right|$','FontSize',5);
        set(gca,'YTickLabel',{'$0$'; '$\underline{\phi}$'; '$\bar{\phi}$'; '$+\infty$'},'YTick',...
            [1, round(nphi*band), round(nphi*(1-band)), nphi],'FontSize',5);
    else
        set(gca,'YTickLabel',[],'FontSize',5);
    end
    xline(round((nht-1)*pix/2),':k','LineWidth',0.5);
    xline(round((nht+1)*pix/2),':k','LineWidth',0.5);
    yline(round(nphi*band),'--k','LineWidth',0.75);
    yline(round(nphi*(1-band)),'--k','LineWidth',0.75);
    axis square;
end
print(f,'-r1000','-dpdf','Fig4');

