% Author: Olivier Loisel (CREST). Date: February 2023.

% This .m file produces Figures 1-4 in the paper "New Principles For Stabilization Policy" (O. Loisel, February 2023).
% In so doing, it also checks the analytical results of Appendices A.1-A.2 with some numerical examples.
% It also produces some additional figures (more specifically, the counterparts of Figures 2 and 4 for phi<0). 

% To run this file, one needs to have Matlab's Optimization Toolbox (in order to use "fminbnd").
% On my computer, the running time is approximately 2 minutes.

% 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. Gabaix (2020)
% 3. Bilbiie (2008)
% 4. 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
% In the paper, Model 1 corresponds to Woodford's calibration, but the results under Gali's calibration are also briefly discussed

% RULES (indexed by rr): setting the interest rate as a function of
% 1. Inflation: rule i_t = phi*E_t{pi_{t+h}}
% 2. Real output: rule i_t = phi*E_t{y_{t+h}}

% SIGN OF PHI IN THE RULE (indexed by ss):
% 1. Positive: phi>0
% 2. Negative: phi<0
sign = 'Phi positive, Phi negative';

% 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} = [-15, -15;
           -10, -10;
           -10, -10;
           -5, -5];
hmin{2} = [-15, -15;
           -40, -40;
           -45, -30;
           -5, -5];
hmax{1} = [5, 5;
           10, 10;
           10, 10;
           15, 15];
hmax{2} = [5, 5;
           10, 10;
           15, 10;
           15, 15];
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(4,2); % minimum of phi in the figures
phimax = NaN(4,2); % maximum of phi in the figures
philowerbar = NaN(4,2);
phiupperbar = NaN(4,2);
phiw = NaN(4,2);

%%%%% MODELS 1 TO 3 %%%%%%
% 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.85; 1];
beta = [0.99; 0.792; 0.99];
sigma = [0.16; 5; -0.11];
kappa = [0.022; 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:3
    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,x+i*sqrt(1-x^2))/polyval(R,x+i*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)]);
        phiw(mm,rr) = -polyval(Q,1)/polyval(R,1);
    end
end

% Check of the phi thresholds using their analytical expression
phiupperbarpia = NaN(3,1);
for mm=1:3
    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(3,1);
for mm=1:3
    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(3,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:3
    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:3,1:2);
checkphiupperbar = phiupperbara - phiupperbar(1:3,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:3],[1 2]) = exp((band*log(phiupperbar([1:3],[1 2]))-(1-band)*log(philowerbar([1:3],[1 2])))/(2*band-1));
phimax([1:3],[1 2]) = exp((band*log(philowerbar([1:3],[1 2]))-(1-band)*log(phiupperbar([1:3],[1 2])))/(2*band-1));

% Determination of the determinacy status
for mm=1:3
    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 4 %%%%%%
% 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 = 4;

% 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,x+i*sqrt(1-x^2))/polyval(R,x+i*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)]);
    phiw(mm,rr) = -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,x+i*sqrt(1-x^2))/polyval(R,x+i*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 AND 4 IN THE PAPER %%%%%%

width = 11;
height = 4.33;
width2 = 11;
height2 = 13;
ss=1;

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;
for rr=1:2
    subplot(1,2,rr);
    contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
    colormap(gray);
    if (rr==1)
        title('Model 1, Rule 1','FontSize',6);
    else
        title('Model 1, Rule 2','FontSize',6);
    end
    caxis([-1, 1]);
    xlabel('$h$','FontSize',6);
    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);
    if (rr==1)
        ylabel('$\phi$ (log scale)','FontSize',6);
    end
    xline(-hmin{ss}(mm,rr)*pix,':k');
    xline((-hmin{ss}(mm,rr)+1)*pix,':k');
    if (rr==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);
        set(gca,'YTickLabel',{round(phimin(mm,rr),2,'significant'); ['$\underline{\phi}=\phi_W=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
            ['$\bar{\phi}=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]; round(phimax(mm,rr),2,'significant')},'YTick',...
            [1, 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)))), nphi],'FontSize',6);
    else
        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);
        set(gca,'YTickLabel',{round(phimin(mm,rr),2,'significant'); ['$\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
            ['$\bar{\phi}=\phi_W=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]; round(phimax(mm,rr),2,'significant')},'YTick',...
            [1, 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)))), nphi],'FontSize',6);
    end
    axis square;
end
print(f,'-r1000','-dpdf','Fig2');


f=figure('Name','Figure 4','Units','centimeters');
set(gcf,'PaperUnits','centimeters','Position',[0 0 width2 height2],'PaperSize',[width2 height2],'PaperPosition',[0 0 width2 height2]);
ax = gca;
ax.TitleFontSizeMultiplier = 1;
for mm=[2 3]
    for rr=1:2
        subplot(3,2,2*(mm-2)+rr);
        contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
        colormap(gray);
        if (rr==1)
            if (mm==2)
                title('Model 2, Rule 1','FontSize',6);
            else
                title('Model 3, Rule 1','FontSize',6);
            end
        else
            if (mm==2)
                title('Model 2, Rule 2','FontSize',6);
            else
                title('Model 3, Rule 2','FontSize',6);
            end
        end
        caxis([-1, 1]);
        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);
        if (rr==1)
            ylabel('$\phi$ (log scale)','FontSize',6);
        end
        xline(-hmin{ss}(mm,rr)*pix,':k');
        xline((-hmin{ss}(mm,rr)+1)*pix,':k');
        if (mm==2)
            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);
            set(gca,'YTickLabel',{round(phimin(mm,rr),2,'significant'); ['$\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
                ['$\bar{\phi}=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]; round(phimax(mm,rr),2,'significant')},'YTick',...
                [1, 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)))), nphi],'FontSize',6);
        else
            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);
            set(gca,'YTickLabel',{round(phimin(mm,rr),2,'significant'); ['$\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
                ['$\bar{\phi}=\phi_W=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]; round(phimax(mm,rr),2,'significant')},'YTick',...
                [1, 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)))), nphi],'FontSize',6);
        end
    end
end
mm=4;
for rr=1:2
    subplot(3,2,4+rr);
    contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
    colormap(gray);
    if (rr==1)
        title('Model 4, Rule 1','FontSize',6);
    else
        title('Model 4, Rule 2','FontSize',6);
    end
    caxis([-1, 1]);
    if (rr==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);
    else
        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);
    end
    xlabel('$h$','FontSize',6);
    xline((-hmin{ss}(mm,rr)+2)*pix,':k');
    xline((-hmin{ss}(mm,rr)+3)*pix,':k');
    if (rr==1)
        ylabel('$\phi$ (log scale)','FontSize',6);
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$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);
        set(gca,'YTickLabel',{round(phimin(mm,rr),2,'significant'); ['$\underline{\phi}=\phi_W=$ ' num2str(round(philowerbar(mm,rr),2,'significant'))];...
            ['$\bar{\phi}=$ ' num2str(round(phiupperbar(mm,rr),2,'significant'))]; round(phimax(mm,rr),2,'significant')},'YTick',...
            [1, 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)))), nphi],'FontSize',6);
    else
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$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);
        set(gca,'YTickLabel',{round(phimin(mm,rr),2,'significant'); ['$\underline{\phi}=$ ' num2str(round(philowerbar(mm,rr),...
            2,'significant'))]; round(phimax(mm,rr),2,'significant')},'YTick',[1, round(1+(nphi-1)*...
            (log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))), nphi],'FontSize',6);
    end
end
axesHandles = findobj(get(gcf,'Children'),'flat','Type','axes');
axis(axesHandles,'square')
print(f,'-r1000','-dpdf','Fig4');



%%%%%% FIGURES 1 AND 3 IN THE PAPER %%%%%%

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); 

width = 11;
height = 4.3;

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);
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}$'; '$\bar{\phi}$'; '$+\infty$'},'YTick',...
    [1, round(nphi*band), round(nphi*(1-band)), nphi],'FontSize',6);
axis square;
print(f,'-r1000','-dpdf','Fig1');


f=figure('Name','Figure 3','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('Models with $S_{peg}=M$','FontSize',5);
    elseif (tt==2)
        contourf(Nv1{3},[0, 0.4, 0.7, 1]);
        title('Models with $S_{peg}=D$','FontSize',5);
    else
        contourf(Nv1{2},[0, 0.4, 0.7, 1]);
        title('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','Fig3');



%%%%%% COUNTERPARTS OF FIGURES 2 AND 4 FOR PHI NEGATIVE %%%%%%

ss=2;

% Figures for Models 1-3

for mm=1:3
    f=figure('Name',['Model',' ',num2str(mm),', Rules 1-2,',' ',sign(14*ss-13:14*ss-2)],'PaperUnits','inches','PaperPosition',[-0.67 0 9.6 3],'PaperSize',[8.1 3]);
    for rr=1:2
        subplot(1,2,rr);
        contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
        colormap(gray);
        title(['Model',' ',num2str(mm),', Rule',' ',num2str(rr)]);            
        caxis([-1, 1]);
        xlabel('$h$');
        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)]);
        ylabel('$\phi$ (log scale)');
        if (mm==2)
            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);
            set(gca,'YTickLabel',{-round(phimax(mm,rr),2,'significant'); ['-$\bar{\phi}=$ ' num2str(-round(phiupperbar(mm,rr),2,'significant'))];...
                ['-$\underline{\phi}=\phi_W=$ ' num2str(-round(philowerbar(mm,rr),2,'significant'))]; -round(phimin(mm,rr),2,'significant')},'YTick',[1, ...
                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)))), nphi]);
        else
            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);
            set(gca,'YTickLabel',{-round(phimax(mm,rr),2,'significant'); ['-$\bar{\phi}=$ ' num2str(-round(phiupperbar(mm,rr),2,'significant'))];...
                ['-$\underline{\phi}=$ ' num2str(-round(philowerbar(mm,rr),2,'significant'))]; -round(phimin(mm,rr),2,'significant')},'YTick',[1, ...
                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)))), nphi]);
        end
        axis square;
    end
    print(f,'-r1000','-dpdf',['Model',' ',num2str(mm),', Rules 1-2,',' ',sign(14*ss-13:14*ss-2)]);
end

% Figures for Model 4

mm=4;

f=figure('Name',['Model',' ',num2str(mm),', Rules 1-2,',' ',sign(14*ss-13:14*ss-2)],'PaperUnits','inches','PaperPosition',[-0.67 0 9.6 3],'PaperSize',[8.1 3]);
for rr=1:2
    subplot(1,2,rr);
    contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
    colormap(gray);
    title(['Model',' ',num2str(mm),', Rule',' ',num2str(rr)]);
    caxis([-1, 1]);
    if (rr==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);
    else
        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);
    end
    xlabel('$h$');
    ylabel('$\phi$ (log scale)');
    if (rr==1)
        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)]);
        set(gca,'YTickLabel',{-round(phimax(mm,rr),2,'significant'); ['-$\bar{\phi}=$ ' num2str(-round(phiupperbar(mm,rr),2,'significant'))];...
            ['-$\underline{\phi}=$ ' num2str(-round(philowerbar(mm,rr),2,'significant'))]; -round(phimin(mm,rr),2,'significant')},'YTick',[1, ...
            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))))]);
    else
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$1$';hmax{ss}(mm,rr)},'XTick',...
        [round(pix/2), (1-hmin{ss}(mm,rr))*pix+round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)]);
        set(gca,'YTickLabel',{-round(phimax(mm,rr),2,'significant'); ['-$\underline{\phi}=$ ' num2str(-round(philowerbar(mm,rr),...
            2,'significant'))]; -round(phimin(mm,rr),2,'significant')},'YTick',[1, nphi+1-round(1+(nphi-1)*...
            (log(philowerbar(mm,rr))-log(phimin(mm,rr)))/(log(phimax(mm,rr))-log(phimin(mm,rr)))), nphi]);
    end
    axis square;
end
print(f,'-r1000','-dpdf',['Model',' ',num2str(mm),', Rules 1-2,',' ',sign(14*ss-13:14*ss-2)]);
