% Author: Olivier Loisel (CREST). Date: August 2022.

% This .m file produces Figures 1-10 in the (longer) working-paper version of the paper "New 
% Principles For Stabilization Policy" (CREST Working Paper No. 2022-16, O. Loisel, August 2022).
% In so doing, it also checks the analytical results of Appendices A.2 and A.3 with some numerical examples.
% Finally, it produces additional figures (altogether, all the figures corresponding to Models 1-6,
% Rules 1-6, and phi>0 or phi<0), a few of which are briefly discussed but not shown in the paper. 

% 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 8 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. McKay, Nakamura, and Steinsson (2017)
% 3. Gabaix (2020)
% 4. Bilbiie (2008)
% 5. Ball (1999) and Svensson (1997) with Ball's (1999) calibration
% 6. Rudebusch and Svensson (1999)

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}}
% 3. Real-output growth: rule i_t = phi*E_t{y_{t+h}-y_{t+h-1}}
% 4. Price level: rule i_t = phi*E_t{p_{t+h}}
% 5. Nominal output: rule i_t = phi*E_t{p_{t+h}+y_{t+h}}
% 6. Nominal-output growth: rule i_t = phi*E_t{pi_{t+h}+y_{t+h}-y_{t+h-1}}

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

%%%%% 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 6]
        Q = [1, -(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm)), alpha(mm)*beta(mm)];
        if (rr==1)
            R = [kappa(mm)/sigma(mm)];
        elseif (rr==2)
            R = [1/sigma(mm), -beta(mm)/sigma(mm)];
        else
            R = [-1/sigma(mm), (1+beta(mm)+kappa(mm))/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
    for rr=3
        Q = [1, -(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm)), alpha(mm)*beta(mm)];
        R = [-1/sigma(mm) (1+beta(mm))/sigma(mm) -beta(mm)/sigma(mm)];
        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)]);
    end
    for rr=[4 5]
        Q = [-1, 1+alpha(mm)+beta(mm)+kappa(mm)/sigma(mm), -(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm)+alpha(mm)*beta(mm)), alpha(mm)*beta(mm)];
        if (rr==4)
            R = [kappa(mm)/sigma(mm)];
        else
            R = [-1/sigma(mm), (1+beta(mm)+kappa(mm))/sigma(mm), -beta(mm)/sigma(mm)];
        end
        nf = @(x) -abs(polyval(Q,x+i*sqrt(1-x^2))/polyval(R,x+i*sqrt(1-x^2)));
        [~,mini] = fminbnd(nf,-1,1,options);
        phiupperbar(mm,rr) = -min([mini nf(1) nf(-1)]);
        phiw(mm,rr) = 0;
    end
end

% Check of the phi thresholds for Rules 1 and 2 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 6]) = exp((band*log(phiupperbar([1:4],[1 2 6]))-(1-band)*log(philowerbar([1:4],[1 2 6])))/(2*band-1));
phimax([1:4],[1 2 6]) = exp((band*log(philowerbar([1:4],[1 2 6]))-(1-band)*log(phiupperbar([1:4],[1 2 6])))/(2*band-1));
phimax([1:4],3) = [25; 35; 7; 1000];
phimin([1:4],3) = exp((log(philowerbar([1:4],3))-band*log(phimax([1:4],3)))/(1-band));
phimin([1:4],[4 5]) = 0;
phimax([1:4],[4 5]) = exp(log(1+phiupperbar([1:4],[4 5]))/(1-band))-1;
if M1==2
    phimax(1,3) = 100;
    phimin(1,6) = 0.1;
    phimax(1,6) = 10;
end

% 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
        for rr=[3 6]
            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)+2);
                    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) = -phi/sigma(mm);
                    P(1,hh+1) = (((rr-3)/3)*kappa(mm)+1+beta(mm))*phi/sigma(mm);
                    P(1,hh+2) = -beta(mm)*phi/sigma(mm);
                    P = P + [zeros(1,-hmin{ss}(mm,rr)+1) 1 -(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm)) alpha(mm)*beta(mm) zeros(1,hmax{ss}(mm,rr)-1)];
                    Roo = sort(abs(roots(P)));
                    Roo = Roo(Roo~=0);
                    nnpv = max(2,hh+hmin{ss}(mm,rr)); % 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)<=1) & (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
        for rr=4:5
            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)+2*rr-8);
                    phi = (3-2*ss)*(exp(log(1+phimin(mm,rr))+(pp-1)*(log(1+phimax(mm,rr))-log(1+phimin(mm,rr)))/(nphi-1))-1); % pseudo-log scale on y axis
                    if (rr==4)
                        P(1,hh) = kappa(mm)*phi/sigma(mm);
                    else
                        P(1,hh) = -phi/sigma(mm);
                        P(1,hh+1) = (kappa(mm)+1+beta(mm))*phi/sigma(mm);
                        P(1,hh+2) = -beta(mm)*phi/sigma(mm);
                    end
                    P = P + [zeros(1,-hmin{ss}(mm,rr)+rr-5) -1 1+alpha(mm)+beta(mm)+kappa(mm)/sigma(mm) ...
                        -(alpha(mm)+beta(mm)+kappa(mm)/sigma(mm)+alpha(mm)*beta(mm)) alpha(mm)*beta(mm) zeros(1,hmax{ss}(mm,rr)+rr-6)];
                    Roo = sort(abs(roots(P)));
                    Roo = Roo(Roo~=0);
                    nnpv = (phi~=0)*max(2,hh+hmin{ss}(mm,rr)+rr-5) + (phi==0)*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 (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,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 3 6]
    Q = [lambda-chi*mu -(1+lambda) 1];
    if (rr==2)
        R = [-mu mu];
    elseif (rr==3)
        R = [mu -2*mu mu];
    else
        R = [mu mu*(chi-2) mu];
        phiw(mm,rr) = -polyval(Q,1)/polyval(R,1);
    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);
    philowerbar(mm,rr) = min([mini f(1) f(-1)]);
end
for rr=4
    Q = [chi*mu-lambda 1+2*lambda-chi*mu -(2+lambda) 1];
    R = [chi*mu];
    nf = @(x) -abs(polyval(Q,x+i*sqrt(1-x^2))/polyval(R,x+i*sqrt(1-x^2)));
    [~,mini] = fminbnd(nf,-1,1);
    phiupperbar(mm,rr) = -min([mini nf(1) nf(-1)]);
    phiw(mm,rr) = 0;
end
for rr=5
    Q = [chi*mu-lambda 1+2*lambda-chi*mu -(2+lambda) 1];
    R = [mu mu*(chi-2) mu];
    f = @(x) abs(polyval(Q,x+i*sqrt(1-x^2))/polyval(R,x+i*sqrt(1-x^2)));
    phiw(mm,rr) = 0;
end
    
% Check of the phi thresholds for Rules 1 and 2 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 3 6]) = [2, 2, 1.7];
phimin(mm,[2 3 6]) = exp((log(philowerbar(mm,[2 3 6]))-band*log(phimax(mm,[2 3 6])))/(1-band));
phimin(mm,4) = 0;
phimax(mm,4) = exp(log(1+phiupperbar(mm,4))/(1-band))-1;
phimin(mm,5) = 0;
phimax(mm,5) = 3;

% 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
    for rr=[3 6]
        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)+2);
                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) = mu*phi;
                P(1,hh+1) = (((rr-3)/3)*chi-2)*mu*phi;
                P(1,hh+2) = mu*phi;
                P = P + [zeros(1,-hmin{ss}(mm,rr)+1) lambda-chi*mu -(1+lambda) 1 zeros(1,hmax{ss}(mm,rr)-1)];
                Roo = sort(abs(roots(P)));
                Roo = Roo(Roo~=0);
                nnpv = max(0,hh+hmin{ss}(mm,rr)-2); % 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
    for rr=4:5
        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)+2*rr-8);
                phi = (3-2*ss)*(exp(log(1+phimin(mm,rr))+(pp-1)*(log(1+phimax(mm,rr))-log(1+phimin(mm,rr)))/(nphi-1))-1); % pseudo-log scale on y axis
                if (rr==4)
                    P(1,hh) = chi*mu*phi;
                else
                    P(1,hh) = mu*phi;
                    P(1,hh+1) = (chi-2)*mu*phi;
                    P(1,hh+2) = mu*phi;
                end
                P = P + [zeros(1,-hmin{ss}(mm,rr)+rr-5) chi*mu-lambda 1+2*lambda-chi*mu -(2+lambda) 1 zeros(1,hmax{ss}(mm,rr)+rr-6)];
                Roo = sort(abs(roots(P)));
                Roo = Roo(Roo~=0);
                nnpv = (phi~=0)*max(0,hh+hmin{ss}(mm,rr)+rr-7); % 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
         

%%%%%% MODEL 6 %%%%%%
% IS equation: y_t = betay1*y_{t-1} + betay2*y_{t-2} - (betar/4)*(i_{t-1}+i_{t-2}+i_{t-3}+i_{t-4}-pi_{t-1}-pi_{t-2}-pi_{t-3}-pi_{t-4})
% Phillips curve: pi_t = alphapi1*pi_{t-1} + alphapi2*pi_{t-2} + alphapi3*pi_{t-3} + alphapi4*pi_{t-4} + alphay*y_{t-1}
mm = 6;

% Estimation
betay1 = 1.16;
betay2 = -0.25;
betar = 0.10;
alphapi1 = 0.70;
alphapi2 = -0.10;
alphapi3 = 0.28;
alphapi4 = 0.12;
alphay = 0.14;

% Computation of the phi thresholds using Matlab's optimization toolbox
for rr=[1 2 3 6]
    Q = [alphapi4*betay2, alphapi4*betay1+alphapi3*betay2-alphay*betar/4, -alphapi4+alphapi3*betay1+alphapi2*betay2-alphay*betar/4, ...
        -alphapi3+alphapi2*betay1+alphapi1*betay2-alphay*betar/4, -alphapi2+alphapi1*betay1-betay2-alphay*betar/4, -alphapi1-betay1, 1];
    if (rr==1)
        R = alphay*betar/4*[1, 1, 1, 1];
    elseif (rr==2)
        R = betar/4*[-alphapi4, -alphapi4-alphapi3, -alphapi4-alphapi3-alphapi2, -alphapi4-alphapi3-alphapi2-alphapi1, ...
        -alphapi3-alphapi2-alphapi1+1, -alphapi2-alphapi1+1, -alphapi1+1, 1];
    elseif (rr==3)
        R = betar/4*[alphapi4, alphapi3, alphapi2, alphapi1, -alphapi4-1, -alphapi3, -alphapi2, -alphapi1, 1];
    else
        R = betar/4*[alphapi4, alphapi3, alphapi2, alphapi1, alphay-alphapi4-1, alphay-alphapi3, alphay-alphapi2, alphay-alphapi1, 1];
    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);
    philowerbar(mm,rr) = min([mini f(1) f(-1)]);
    if (rr==1)|(rr==6)
        phiw(mm,rr) = -polyval(Q,1)/polyval(R,1);
    end
end
for rr=[4 5]
    phiw(mm,rr) = 0;
end
            
% Setting the bounds of phi in the figures
phimax(mm,:) = [10000, 1000, 100, 50000, 1000, 200];
phimin(mm,[1 2 3 6]) = exp((log(philowerbar(mm,[1 2 3 6]))-band*log(phimax(mm,[1 2 3 6])))/(1-band));
phimin(mm,[4 5]) = 0;

% 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
                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 = zeros(1,nh{ss}(mm,rr)+4*rr-1);
                if (rr==1)
                    P(1,hh:hh+3) = alphay*betar*phi/4;
                else
                    P(1,hh:hh+7) = betar*phi/4*[-alphapi4, -alphapi4-alphapi3, -alphapi4-alphapi3-alphapi2, ...
                        -alphapi4-alphapi3-alphapi2-alphapi1, -alphapi3-alphapi2-alphapi1+1, -alphapi2-alphapi1+1, -alphapi1+1, 1];
                end
                P = P + [zeros(1,-hmin{ss}(mm,rr)+3*rr-4), alphapi4*betay2, alphapi4*betay1+alphapi3*betay2-alphay*betar/4, ...
                    -alphapi4+alphapi3*betay1+alphapi2*betay2-alphay*betar/4, -alphapi3+alphapi2*betay1+alphapi1*betay2-alphay*betar/4, ...
                    -alphapi2+alphapi1*betay1-betay2-alphay*betar/4, -alphapi1-betay1, 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
    for rr=4:5
        M=ones(nphi,nh{ss}(mm,rr));
        for hh=1:nh{ss}(mm,rr)
            for pp=1:nphi
                phi = (3-2*ss)*(exp(log(1+phimin(mm,rr))+(pp-1)*(log(1+phimax(mm,rr))-log(1+phimin(mm,rr)))/(nphi-1))-1); % pseudo-log scale on y axis
                P = zeros(1,nh{ss}(mm,rr)+5*rr-17);
                if (rr==4)
                    P(1,hh:hh+3) = alphay*betar*phi/4;
                else
                    P(1,hh:hh+8) = betar*phi/4*[alphapi4, alphapi3, alphapi2, alphapi1, ...
                        alphay-1-alphapi4, alphay-alphapi3, alphay-alphapi2, alphay-alphapi1, 1];
                end
                P = P + [zeros(1,-hmin{ss}(mm,rr)+4*rr-18), -alphapi4*betay2, ...
                    alphapi4*betay2-alphapi4*betay1-alphapi3*betay2+alphay*betar/4, ...
                    alphapi4*betay1+alphapi3*betay2+alphapi4-alphapi3*betay1-alphapi2*betay2, ...
                    -alphapi4+alphapi3*betay1+alphapi2*betay2+alphapi3-alphapi2*betay1-alphapi1*betay2, ...
                    -alphapi3+alphapi2*betay1+alphapi1*betay2+alphapi2-alphapi1*betay1+betay2, ...
                    -alphapi2+alphapi1*betay1-betay2-alphay*betar/4+alphapi1+betay1, ...
                    -alphapi1-betay1-1, 1, zeros(1,hmax{ss}(mm,rr)+rr-6)];
                Roo = sort(abs(roots(P)));
                Roo = Roo(Roo~=0);
                nnpv = (phi~=0)*max(0,hh+hmin{ss}(mm,rr)+rr-7); % 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
    for rr=[3 6]
        M=ones(nphi,nh{ss}(mm,rr));
        for hh=1:nh{ss}(mm,rr)
            for pp=1:nphi
                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 = zeros(1,nh{ss}(mm,rr)+8);
                if (rr==3)
                    P(1,hh:hh+8) = betar*phi/4*[alphapi4, alphapi3, alphapi2, alphapi1, ...
                        -1-alphapi4, -alphapi3, -alphapi2, -alphapi1, 1];
                else
                    P(1,hh:hh+8) = betar*phi/4*[alphapi4, alphapi3, alphapi2, alphapi1, ...
                        alphay-1-alphapi4, alphay-alphapi3, alphay-alphapi2, alphay-alphapi1, 1];
                end
                P = P + [zeros(1,-hmin{ss}(mm,rr)+3), alphapi4*betay2, alphapi4*betay1+alphapi3*betay2-alphay*betar/4, ...
                    -alphapi4+alphapi3*betay1+alphapi2*betay2-alphay*betar/4, -alphapi3+alphapi2*betay1+alphapi1*betay2-alphay*betar/4, ...
                    -alphapi2+alphapi1*betay1-betay2-alphay*betar/4, -alphapi1-betay1, 1, zeros(1,hmax{ss}(mm,rr)-1)];
                Roo = sort(abs(roots(P)));
                Roo = Roo(Roo~=0);
                nnpv = max(0,hh+hmin{ss}(mm,rr)-2); % 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, 4, 6, 8, AND 10 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=[3 4]
    for rr=1:2
        subplot(3,2,2*(mm-3)+rr);
        contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
        colormap(gray);
        if (rr==1)
            if (mm==3)
                title('Model 3, Rule 1','FontSize',6);
            else
                title('Model 4, Rule 1','FontSize',6);
            end
        else
            if (mm==3)
                title('Model 3, Rule 2','FontSize',6);
            else
                title('Model 4, 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==3)
            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=5;
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 5, Rule 1','FontSize',6);
    else
        title('Model 5, 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');


f=figure('Name','Figure 6','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=1:6
    for rr=3
        subplot(3,2,mm);
        contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
        colormap(gray);
        if (mm==1)
            title('Model 1, Rule 3','FontSize',6);
        elseif (mm==2)
            title('Model 2, Rule 3','FontSize',6);
        elseif (mm==3)
            title('Model 3, Rule 3','FontSize',6);
        elseif (mm==4)
            title('Model 4, Rule 3','FontSize',6);
        elseif (mm==5)
            title('Model 5, Rule 3','FontSize',6);
        elseif (mm==6)
            title('Model 6, Rule 3','FontSize',6);
        end
        caxis([-1, 1]);
        if ((mm==5)|(mm==6))
            xlabel('$h$','FontSize',6);
        end
        if ((mm==1)|(mm==2)|(mm==4))
            set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$1$';hmax{ss}(mm,rr)},'XTick',...
                [round(pix/2), -hmin{ss}(mm,rr)*pix+round(pix/2)+pix, (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
            xline((-hmin{ss}(mm,rr)+1)*pix,':k');
            xline((-hmin{ss}(mm,rr)+2)*pix,':k');
        elseif (mm==3)
            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);
            xline(-hmin{ss}(mm,rr)*pix,':k');
            xline((-hmin{ss}(mm,rr)+1)*pix,':k');
        elseif (mm==5)
            set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$2$';hmax{ss}(mm,rr)},'XTick',...
                [round(pix/2), -hmin{ss}(mm,rr)*pix+round(pix/2)+2*pix, (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
            xline((-hmin{ss}(mm,rr)+2)*pix,':k');
            xline((-hmin{ss}(mm,rr)+3)*pix,':k');
        else
            set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);hmax{ss}(mm,rr)},'XTick',...
                [round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
        end
        if ((mm==1)|(mm==3)|(mm==5))
            ylabel('$\phi$ (log scale)','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);
        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','Fig6');


f=figure('Name','Figure 8','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=[1 4 5]
    for rr=4:5
        if (mm==1)
            subplot(3,2,rr-3);
        else
            subplot(3,2,2*mm+rr-9);
        end
        contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
        colormap(gray);
        if ((mm==1) & (rr==4))
            title('Model 1, Rule 4','FontSize',6);
        elseif ((mm==1) & (rr==5))
            title('Model 1, Rule 5','FontSize',6);
        elseif ((mm==4) & (rr==4))
            title('Model 4, Rule 4','FontSize',6);
        elseif ((mm==4) & (rr==5))
            title('Model 4, Rule 5','FontSize',6);
        elseif ((mm==5) & (rr==4))
            title('Model 5, Rule 4','FontSize',6);
        elseif ((mm==5) & (rr==5))
            title('Model 5, Rule 5','FontSize',6);
        end
        caxis([-1, 1]);
        if (mm==5)
            xlabel('$h$','FontSize',6);
        end
        if ((mm==1)|(mm==4))
            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);
            xline((-hmin{ss}(mm,rr))*pix,':k');
            xline((-hmin{ss}(mm,rr)+1)*pix,':k');
        elseif (rr==4)
            set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$2$';hmax{ss}(mm,rr)},'XTick',...
                [round(pix/2), (-hmin{ss}(mm,rr)+2)*pix+round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
            xline((-hmin{ss}(mm,rr)+2)*pix,':k');
            xline((-hmin{ss}(mm,rr)+3)*pix,':k');
        else
            set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);hmax{ss}(mm,rr)},'XTick',...
                [round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
        end
        if (rr==4)
            ylabel('$\phi$ (pseudo-log scale)','FontSize',6);
        end
        if ((mm==5) & (rr==5))
            set(gca,'YTickLabel',{['$\underline{\phi}=\phi_W=$ ' num2str(round(phimin(mm,rr),2,'significant'))]; ...
                round(phimax(mm,rr),2,'significant')},'YTick',[1, nphi],'FontSize',6);
        else
            yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(1+phiupperbar(mm,rr))-log(1+phimin(mm,rr)))/...
                (log(1+phimax(mm,rr))-log(1+phimin(mm,rr)))),'--k','LineWidth',1);
            set(gca,'YTickLabel',{['$\underline{\phi}=\phi_W=$ ' num2str(round(phimin(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(1+phiupperbar(mm,rr))-log(1+phimin(mm,rr)))/(log(1+phimax(mm,rr))-log(1+phimin(mm,rr)))), nphi],'FontSize',6);
        end
    end
end
axesHandles = findobj(get(gcf,'Children'),'flat','Type','axes');
axis(axesHandles,'square')
print(f,'-r1000','-dpdf','Fig8');


f=figure('Name','Figure 10','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=6;
for rr=4:5
    subplot(1,2,rr-3);
    contourf(N{mm,rr,ss},[0.4, 0.7, 1]);
    colormap(gray);
    if (rr==4)
        title('Model 6, Rule 4','FontSize',6);
    else
        title('Model 6, Rule 5','FontSize',6);
    end
    caxis([-1, 1]);
    xlabel('$h$','FontSize',6);
    if (rr==4)
        ylabel('$\phi$ (pseudo-log scale)','FontSize',6);
        xline((-hmin{ss}(mm,rr)+3)*pix,':k');
        xline((-hmin{ss}(mm,rr)+4)*pix,':k');
        xline((-hmin{ss}(mm,rr)+5)*pix,':k');
        xline((-hmin{ss}(mm,rr)+6)*pix,':k');
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$3$';'$5$';hmax{ss}(mm,rr)},'XTick',...
            [round(pix/2), (3-hmin{ss}(mm,rr))*pix+round(pix/2), (5-hmin{ss}(mm,rr))*pix+round(pix/2), ...
            (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
    else
        xline((-hmin{ss}(mm,rr)+2)*pix,':k');
        xline((-hmin{ss}(mm,rr)+3)*pix,':k');
        xline((-hmin{ss}(mm,rr)+4)*pix,':k');
        xline((-hmin{ss}(mm,rr)+5)*pix,':k');
        set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$2$';'$4$';hmax{ss}(mm,rr)},'XTick',...
            [round(pix/2), (2-hmin{ss}(mm,rr))*pix+round(pix/2), (4-hmin{ss}(mm,rr))*pix+round(pix/2), ...
            (nh{ss}(mm,rr)-1)*pix+round(pix/2)],'FontSize',6);
    end
    set(gca,'YTickLabel',{'$\underline{\phi}=\phi_W=0$'; round(phimax(mm,rr),2,'significant')},'YTick',[1, nphi],'FontSize',6);
    axis square;
end
print(f,'-r1000','-dpdf','Fig10');



%%%%%% FIGURES 1, 3, 5, 7, AND 9 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); 
    
for tt=1:2
    Nv2{tt} = ones(nphi,nht);
    Nv2{tt}(1:round(nphi*band),:) = 1-0.3*tt;
    Nv2{tt}(round(nphi*band)+5:round(nphi*(1-band))-1,-19+20*tt) = 0.1+0.3*tt;
    Nv2{tt}(round(nphi*band)+round(factor1)+5:round(nphi*(1-band))-1,-16+18*tt) = 0.1+0.3*tt;
    Nv2{tt}(round(nphi*band)+round(factor1^2)+5:round(nphi*(1-band))-1,-13+16*tt) = 0.1+0.3*tt;
    Nv2{tt}(round(nphi*band)+round(factor1^3)+5:round(nphi*(1-band))-1,-10+14*tt) = 0.1+0.3*tt;
    Nv2{tt}(round(nphi*band)+round(factor1^4)+5:round(nphi*(1-band))-1,-7+12*tt) = 0.1+0.3*tt;
    Nv2{tt}(round(nphi*band)+round(factor1^5)+5:round(nphi*(1-band))-1,-4+10*tt) = 0.1+0.3*tt;
    Nv2{tt}(round(nphi*band)+round(factor1^6)+5:round(nphi*(1-band))-1,-1+8*tt) = 0.1+0.3*tt;
    Nv2{tt}(round(nphi*band)+1:round(nphi*(1-band))-1,29-14*tt:35-14*tt) = 1-0.3*tt;
    Nv2{tt}(round(nphi*(1-band)):nphi,1:9) = 0.4;
    Nv2{tt}(round(nphi*(1-band)):nphi,13:nht) = 0.7;
    Nv2{tt} = kron(Nv2{tt},ones(1,pix));
    for uu=1:2
        for ww=1
            Nv2{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
            Nv2{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
            Nv2{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
            Nv2{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
            Nv2{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
            Nv2{tt}(round(nphi*band)+5*ww,7*pix+mod(ww+uu-1,10)*10:100:(14+0*tt)*pix-10) = 0.1;
        end
    end
    for uu=10:16
        for vv=1:12
            for ww=round(nphi*(1-2*band)/30):round(nphi*(1-band)/30-1)
                Nv2{tt}(round(nphi*band)+30*ww+2+uu,9*pix+vv+44+mod(ww,2)*100:200:(12+0*tt)*pix) = 0.1;
            end
        end
    end
end
Nv2{3}(1:round(nphi*band),:) = 1;
Nv2{3}(round(nphi*band)+1:nphi,1:10.5*pix) = Nv2{1}(round(nphi*band)+1:nphi,1:10.5*pix);  
Nv2{3}(round(nphi*band)+1:nphi,10.5*pix+1:21*pix) = Nv2{2}(round(nphi*band)+1:nphi,10.5*pix+1:21*pix); 

for tt=1:2
    Nv3{tt} = ones(nphi,nht);
    Nv3{tt}(1+5:round(nphi*(1-band))-1,-19+20*tt) = 0.1+0.3*tt;
    Nv3{tt}(1+round(factor2)+5:round(nphi*(1-band))-1,-16+18*tt) = 0.1+0.3*tt;
    Nv3{tt}(1+round(factor2^2)+5:round(nphi*(1-band))-1,-13+16*tt) = 0.1+0.3*tt;
    Nv3{tt}(1+round(factor2^3)+5:round(nphi*(1-band))-1,-10+14*tt) = 0.1+0.3*tt;
    Nv3{tt}(1+round(factor2^4)+5:round(nphi*(1-band))-1,-7+12*tt) = 0.1+0.3*tt;
    Nv3{tt}(1+round(factor2^5)+5:round(nphi*(1-band))-1,-4+10*tt) = 0.1+0.3*tt;
    Nv3{tt}(1+round(factor2^6)+5:round(nphi*(1-band))-1,-1+8*tt) = 0.1+0.3*tt;
    Nv3{tt}(1+1:round(nphi*(1-band))-1,29-14*tt:35-14*tt) = 1-0.3*tt;
    Nv3{tt}(round(nphi*(1-band)):nphi,1:10) = 0.4;
    Nv3{tt}(round(nphi*(1-band)):nphi,12:nht) = 0.7;
    Nv3{tt} = kron(Nv3{tt},ones(1,pix));
    for uu=1:2
        for ww=14:25
            Nv3{tt}(5*ww,(3+2*tt)*pix+mod(ww+uu-1,10)*10:100:(12+2*tt)*pix-10) = 0.1;
        end
        for ww=26:65
            Nv3{tt}(5*ww,(5+1*tt)*pix+mod(ww+uu-1,10)*10:100:(13+1*tt)*pix-10) = 0.1;
        end
        for ww=66:round(nphi*(1-band)/5)-1
            Nv3{tt}(5*ww,7*pix+mod(ww+uu-1,10)*10:100:(14+0*tt)*pix-10) = 0.1;
        end
    end
end
Nv3{3}(1:nphi,1:10.5*pix) = Nv3{1}(1:nphi,1:10.5*pix);  
Nv3{3}(1:nphi,10.5*pix+1:21*pix) = Nv3{2}(1:nphi,10.5*pix+1:21*pix); 
for uu=10:16
    for vv=15:26
        Nv3{3}(uu,[3*pix+30+vv,17*pix+30+vv]) = 0.1;
        Nv3{3}(uu+5,5*pix+30+vv:200:16*pix-10) = 0.1;
        Nv3{3}(uu+30,[4*pix+30+vv,16*pix+30+vv]) = 0.1;
        Nv3{3}(uu+35,6*pix+30+vv:200:15*pix-10) = 0.1;
    end
end
    
for tt=1:2
    Nv4{tt} = ones(nphi,nht);
    Nv4{tt}(1+5:round(nphi*(1-band))-1,-19+20*tt) = 0.1+0.3*tt;
    Nv4{tt}(1+round(factor2)+5:round(nphi*(1-band))-1,-16+18*tt) = 0.1+0.3*tt;
    Nv4{tt}(1+round(factor2^2)+5:round(nphi*(1-band))-1,-13+16*tt) = 0.1+0.3*tt;
    Nv4{tt}(1+round(factor2^3)+5:round(nphi*(1-band))-1,-10+14*tt) = 0.1+0.3*tt;
    Nv4{tt}(1+round(factor2^4)+5:round(nphi*(1-band))-1,-7+12*tt) = 0.1+0.3*tt;
    Nv4{tt}(1+round(factor2^5)+5:round(nphi*(1-band))-1,-4+10*tt) = 0.1+0.3*tt;
    Nv4{tt}(1+round(factor2^6)+5:round(nphi*(1-band))-1,-1+8*tt) = 0.1+0.3*tt;
    Nv4{tt}(1+1:round(nphi*(1-band))-1,29-14*tt:35-14*tt) = 1-0.3*tt;
    Nv4{tt}(round(nphi*(1-band)):nphi,1:9) = 0.4;
    Nv4{tt}(round(nphi*(1-band)):nphi,13:nht) = 0.7;
    Nv4{tt} = kron(Nv4{tt},ones(1,pix));
    for uu=1:2
        for ww=14:25
            Nv4{tt}(5*ww,(3+2*tt)*pix+mod(ww+uu-1,10)*10:100:(12+2*tt)*pix-10) = 0.1;
        end
        for ww=26:65
            Nv4{tt}(5*ww,(5+1*tt)*pix+mod(ww+uu-1,10)*10:100:(13+1*tt)*pix-10) = 0.1;
        end
        for ww=66:round(nphi*(1-band)/5)-1
            Nv4{tt}(5*ww,7*pix+mod(ww+uu-1,10)*10:100:(14+0*tt)*pix-10) = 0.1;
        end
    end
    for uu=10:16
        for vv=1:12
            for ww=round(nphi*(1-2*band)/30):round(nphi*(1-band)/30-1)
                Nv4{tt}(round(nphi*band)+30*ww+2+uu,9*pix+vv+44+mod(ww,2)*100:200:(12+0*tt)*pix) = 0.1;
            end
        end
    end
end
Nv4{3}(1+1:nphi,1:10.5*pix) = Nv4{1}(1+1:nphi,1:10.5*pix);  
Nv4{3}(1+1:nphi,10.5*pix+1:21*pix) = Nv4{2}(1+1:nphi,10.5*pix+1:21*pix);  
for uu=10:16
    for vv=15:26
        Nv4{3}(uu,[3*pix+30+vv,17*pix+30+vv]) = 0.1;
        Nv4{3}(uu+5,5*pix+30+vv:200:16*pix-10) = 0.1;
        Nv4{3}(uu+30,[4*pix+30+vv,16*pix+30+vv]) = 0.1;
        Nv4{3}(uu+35,6*pix+30+vv:200:15*pix-10) = 0.1;
    end
end
   

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 $q\leq \delta-1$','FontSize',5);
    elseif (tt==2)
        contourf(Nv1{3},[0, 0.4, 0.7, 1]);
        title('Models with $q=\delta$','FontSize',5);
    else
        contourf(Nv1{2},[0, 0.4, 0.7, 1]);
        title('Models with $q\geq \delta+1$','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');


f=figure('Name','Figure 5','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(Nv2{tt},[0, 0.4, 0.7, 1]);
        title('Models with $q\leq \delta-1$','FontSize',5);
    elseif (tt==2)
        contourf(Nv2{3},[0, 0.4, 0.7, 1]);
        title('Models with $q=\delta$','FontSize',5);
    else
        contourf(Nv2{2},[0, 0.4, 0.7, 1]);
        title('Models with $q\geq \delta+1$','FontSize',5);
    end
    caxis([-1, 1]);
    xlabel('$h$','FontSize',5);
    set(gca,'Ticklength', [0 0],'XTickLabel',{'$-\infty$';'$\underline{h}^\ast$';'$\bar{h}^\ast$';'$+\infty$'},'XTick',...
        [round(pix/2), round(nht*pix/2-pix), round(nht*pix/2+pix), (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-pix),':k','LineWidth',0.5);
    xline(round((nht-1)*pix/2),':k','LineWidth',0.5);
    xline(round((nht+1)*pix/2),':k','LineWidth',0.5);
    xline(round((nht+1)*pix/2+pix),':k','LineWidth',0.5);
    yline(round(nphi*band),'--k','LineWidth',0.75);
    yline(round(nphi*(1-band)),':k','LineWidth',0.75);
    line([900 1200], [round(nphi*(1-band)) round(nphi*(1-band))],'Color','k','LineWidth',0.5);
    axis square;
end
print(f,'-r1000','-dpdf','Fig5');


f=figure('Name','Figure 7','Units','centimeters');
set(gcf,'PaperUnits','centimeters','Position',[0 0 width height],'PaperSize',[width height],'PaperPosition',[0 0 width height]);
contourf(Nv3{3},[0, 0.4, 0.7, 1]);
colormap(gray);
caxis([-1, 1]);
xlabel('$h$','FontSize',6);
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',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*(1-band)),'--k','LineWidth',0.75);
line([500 1600], [1+round(factor2^4)+20 1+round(factor2^4)+20],'Color','k','LineWidth',0.5);
set(gca,'YTickLabel',{'$\underline{\phi}=0$'; '$\bar{\phi}$'; '$+\infty$'},'YTick',...
    [1, round(nphi*(1-band)), nphi],'FontSize',6);
axis square;
print(f,'-r1000','-dpdf','Fig7');


f=figure('Name','Figure 9','Units','centimeters');
set(gcf,'PaperUnits','centimeters','Position',[0 0 width height],'PaperSize',[width height],'PaperPosition',[0 0 width height]);
contourf(Nv4{3},[0, 0.4, 0.7, 1]);
colormap(gray);
caxis([-1, 1]);
xlabel('$h$','FontSize',6);
set(gca,'Ticklength', [0 0],'XTickLabel',{'$-\infty$';'$\underline{h}^\ast$';'$\bar{h}^\ast$';'$+\infty$'},'XTick',...
    [round(pix/2), round(nht*pix/2-pix), round(nht*pix/2+pix), (nht-1)*pix+round(pix/2)],'FontSize',6);
ylabel('$\left|\phi\right|$','FontSize',6);
xline(round((nht-1)*pix/2-pix),':k','LineWidth',0.5);
xline(round((nht-1)*pix/2),':k','LineWidth',0.5);
xline(round((nht+1)*pix/2),':k','LineWidth',0.5);
xline(round((nht+1)*pix/2+pix),':k','LineWidth',0.5);
yline(round(nphi*(1-band)),':k','LineWidth',0.75);
line([900 1200], [round(nphi*(1-band)) round(nphi*(1-band))],'Color','k','LineWidth',0.5);
line([500 1600], [1+round(factor2^4)+20 1+round(factor2^4)+20],'Color','k','LineWidth',0.5);
set(gca,'YTickLabel',{'$\underline{\phi}=0$'; '$\bar{\phi}$'; '$+\infty$'},'YTick',...
    [1, round(nphi*(1-band)), nphi],'FontSize',6);
axis square;
print(f,'-r1000','-dpdf','Fig9');



%%%%%% ALL THE FIGURES FOR MODELS 1-6, RULES 1-6, PHI POSITIVE OR NEGATIVE %%%%%%

% Figures for Models 1-4
for ss=1:2  
    for mm=1: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]);
            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 (ss==1)
                xline(-hmin{ss}(mm,rr)*pix,':k');
                xline((-hmin{ss}(mm,rr)+1)*pix,':k');
                if ((mm<=2) & (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]);
                elseif ((mm==4)|((mm==1)&(rr==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}=\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]);
                elseif ((mm==2) & (rr==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);
                    yline((ss-1)*(nphi+1)+(3-2*ss)*round(1+(nphi-1)*(log(phiw(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'))];...
                        ['$\phi_W=$ ' num2str(round(phiw(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(phiw(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]);
                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}=$ ' 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]);
                end
            else
                if (mm==3)
                    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
            end
            axis square;
        end
        print(f,'-r1000','-dpdf',['Model',' ',num2str(mm),', Rules 1-2,',' ',sign(14*ss-13:14*ss-2)]);
    end
end

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

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

% Figures for Model 5

mm=5;

for ss=1:2
    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) & (ss==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);
        elseif ((rr==1) & (ss==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);
        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$');
        if (ss==1)
            xline((-hmin{ss}(mm,rr)+2)*pix,':k');
            xline((-hmin{ss}(mm,rr)+3)*pix,':k');
        end
        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)]);
            if (ss==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]);
            else
                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))))]);
            end
        else
            if (ss==1)
                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)]);
                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]);
            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
        end
        axis square;
    end
    print(f,'-r1000','-dpdf',['Model',' ',num2str(mm),', Rules 1-2,',' ',sign(14*ss-13:14*ss-2)]);
end

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

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

% Figures for Model 6

mm=6;

for ss=1:2
    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$');
        ylabel('$\phi$ (log scale)');
        if ((rr==1) & (ss==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);
        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
        if (rr==1)
            if (ss==1)
                xline((-hmin{ss}(mm,rr)+3)*pix,':k');
                xline((-hmin{ss}(mm,rr)+4)*pix,':k');
                xline((-hmin{ss}(mm,rr)+5)*pix,':k');
                xline((-hmin{ss}(mm,rr)+6)*pix,':k');
                set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);'$3$';'$5$';hmax{ss}(mm,rr)},'XTick',...
                    [round(pix/2), (3-hmin{ss}(mm,rr))*pix+round(pix/2), (5-hmin{ss}(mm,rr))*pix+round(pix/2), ...
                    (nh{ss}(mm,rr)-1)*pix+round(pix/2)]);
            else
                set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);hmax{ss}(mm,rr)},'XTick',...
                    [round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)]);
            end
        else
            set(gca,'Ticklength', [0 0],'XTickLabel',{hmin{ss}(mm,rr);hmax{ss}(mm,rr)},'XTick',...
                [round(pix/2), (nh{ss}(mm,rr)-1)*pix+round(pix/2)]);
        end
        if ((rr==1) & (ss==1))
            set(gca,'YTickLabel',{round(phimin(mm,rr),2,'significant'); ['$\underline{\phi}=\phi_W=$ ' 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]);
        elseif ((rr==2) & (ss==1))
            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]);
        else
            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)]);
end

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

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