% fig1.m

% Script for fig 1.

% Calls p0f, p1f, qsf.

N=1000; R0=1.17; alpha1=1; alpha2=1; nmax=160;

p0=p0f(N,R0,alpha1,alpha2);

p1=p1f(N,R0,alpha1,alpha2);

qs=qsf(N,R0,alpha1,alpha2);

n=1:nmax;

K1=(R0-1)*N/(alpha1*R0+alpha2);

sig1=sqrt((alpha1+alpha2)*R0*N)/(alpha1*R0+alpha2);

y1=(n-K1)/sig1;

p1a=phi(y1)/sig1;

plot(n,p0(n),'r',n,p1(n),'g',n,qs(n),'b',n,p1a,':g')

%w=axis; w(4)=0.04; axis(w), grid

% fig2.m

% Script for fig 2.

% Show the left tails of p1, p0, qs, and their approximations, for R0>1.

% I have two approximations for each of p0 and p1.

% The approximations that hold for n=o(sqrt(N)) are called p0a1 and p1a1.

% The better approximations that hold for n=O(sqrt(N)) are called p0a2 and p1a2.

% The approximation of qs that holds for n=o(sqrt(N)) is called qsa1.

% Calls p0f, p1f, qsf, beta1f, phi.

N=1000; R0=1.16; alpha1=1; alpha2=1; nmax=20;

p1=p1f(N,R0,alpha1,alpha2);

p0=p0f(N,R0,alpha1,alpha2);

qs=qsf(N,R0,alpha1,alpha2);

beta1=beta1f(N,R0,alpha1,alpha2);

mu2=log(R0)*N/(alpha1+alpha2); sig2=sqrt(N)/sqrt(alpha1+alpha2); beta2=mu2/sig2;

a=(alpha1*R0+alpha2)/sqrt(alpha1+alpha2);

n=1:nmax;

y2=(n-mu2)/sig2;

a1=a*phi(beta1)/(R0*sqrt(N));

p1a1=a1*R0.^n;

a2=a1/phi(beta2);

p1a2=a2*phi(y2);

b1=(R0-1)*sqrt(N)*phi(beta1)/a;

p0a1=b1*R0.^n./n;

b2=b1/phi(beta2);

p0a2=b2*phi(y2)./n;

qsa1=b1*(R0.^n-1)./n;

K1=(R0-1)*N/(alpha1*R0+alpha2);

sig1=sqrt(R0*N)/a;

y1=(n-K1)/sig1;

tail=phi(y1)/sig1;

semilogy(n,p0(n),'r',n,p1(n),'g',n,qs(n),'b',n,tail,'m',n,p0a2,':r',n,p1a2,':g',n,qsa1,':b')

w=axis; w(3)=0.00005; axis(w)

% fig3.m

% Script for plotting fig3: Distributions p0, p1, qs for R0<1.

% Calls p0f, p1f, qsf.

N=1000; R0=0.9; alpha1=1; alpha2=1; nmax=25;

p0=p0f(N,R0,alpha1,alpha2);

p1=p1f(N,R0,alpha1,alpha2);

qs=qsf(N,R0,alpha1,alpha2);

n=1:nmax;

p1a=(1-R0)*R0.^(n-1);

p0a=-(1/log(1-R0))*R0.^n./n;

plot(n,p0(n),'r',n,p1(n),'g',n,qs(n),'b',n,p0a,':r',n,p1a,':g')

w=axis; w(4)=0.15; axis(w)

% fig4.m

% Script for fig4. Plots H(rho), H0(rho) and H1(rho) for rho<0.

% Calls h0f, h1f, hf.

rhob=-3; alpha=1;

rho=-3:0.1:0;

h0=h0f(rho,rhob);

h1=h1f(rho);

h=hf(rho,alpha);

plot(rho,h0,rho,h1,rho,h)

% fig5.m

% Script for fig5. Plots H(rho), H0(rho) and H1(rho) for rho>0.

% Calls h0f, h1f, hf.

rhob=-3; alpha=0.2;

rho=0:0.1:3;

h0=h0f(rho,rhob);

h1=h1f(rho);

h=hf(rho,alpha);

semilogy(rho,h0,'r',rho,h1,'g',rho,h,'b')

w=axis; w(3)=0.5; w(4)=250; axis(w)

% fig6.m

% Script for fig6.

% Plots p0, p1, qs, and their approximations, in the transition region.

% Calls p0f, p1f, qsf, phi, capphi, h0f, hf.

rho=2; N=1000; R0=1+rho*sqrt(alpha1+alpha2)/sqrt(N); alpha1=1; alpha2=1;

nmax=50; rhob=-3;

p0=p0f(N,R0,alpha1,alpha2);

p1=p1f(N,R0,alpha1,alpha2);

qs=qsf(N,R0,alpha1,alpha2);

l=min([length(p0),length(p1),length(qs)]);

n=1:nmax;

sig3=sqrt(N/(alpha1+alpha2)); mu3=rho*sig3; y3=(n-mu3)/sig3;

p1a=phi(y3)/(sig3*capphi(rho));

sig1=sqrt((alpha1+alpha2)*R0*N)/(alpha1*R0+alpha2);

a=(log(sig3)+h0f(rho,rhob))*phi(rho);

p0a=(1/a)*phi(y3)./n;

R1=1+(rho*sqrt(alpha1+alpha2)+1/hf(rho,alpha1+alpha2))/sqrt(N);

b=phi(rho)*(1+rho*sqrt(alpha1+alpha2)*hf(rho,alpha1+alpha2));

qsa=(1/b)*(1-1./R1.^n).*phi(y3)./n;

plot(n,p0(n),'r',n,p1(n),'g',n,qs(n),'b',n,p0a,':r',n,p1a,':g',n,qsa,':b')

w=axis; w(4)=0.2; axis(w)

% fig7.m

% Script for fig7.

% Plots Expectations of X0, X1, XQ, and their approximations, in the transition region.

% Calls p0f, p1f, qsf, h0f, h1f, hf.

rho=-3:0.2:3; N=1000; alpha1=0.5; alpha2=0.5; rhob=-3;

l=length(rho); ex0=zeros(size(rho)); ex1=ex0; exq=ex0; k1=ex0;

ex0a=ex0; ex1a=ex0; exqa=ex0;

alpha=alpha1+alpha2; M=N/alpha;

for j=1:l

R0=1+rho(j)/sqrt(M);

p0=p0f(N,R0,alpha1,alpha2);

n=1:length(p0);

ex0(j)=sum(n.*p0);

p1=p1f(N,R0,alpha1,alpha2);

n=1:length(p1);

ex1(j)=sum(n.*p1);

qs=qsf(N,R0,alpha1,alpha2);

n=1:length(qs);

exq(j)=sum(n.*qs);

h1=h1f(rho(j)); h0=h0f(rho(j),rhob); h=hf(rho(j),alpha);

ex1a(j)=(rho(j)+1/h1)*sqrt(M);

ex0a(j)=h1*sqrt(M)/(0.5*log(M)+h0);

k=1/(h*sqrt(alpha));

exqa(j)=(h1-h1f(-k))*sqrt(M)/(1+rho(j)/k);

if j<(l+1)/2

k1(j)=0;

else

k1(j)=rho(j)*sqrt(M);

end

end

plot(rho,ex0,'r',rho,ex1,'g',rho,exq,'b',rho,ex0a,':r',rho,ex1a,':g')

hold on

plot(rho,exqa,':b',rho,k1,'b')

%plot(rho,ex0,rho,ex1,rho,exq,rho,ex0a,':',rho,ex1a,':',rho,exqa,':',rho,k1)

% fig8.m

% Plot T[tauQ] as function of rho.

% Calls qsf, beta1f, hf.

N=1000; alpha1=0.5; alpha2=0.5; alpha=alpha1+alpha2; M=N/alpha;

rho=-5:0.2:5; R0=1+rho/sqrt(M);

et=zeros(size(rho));

for k=1:length(rho)

qs=qsf(N,R0(k),alpha1,alpha2);

et(k)=1/qs(1);

end

rho1=3:0.2:5; R1=1+rho1/sqrt(M);

beta1=beta1f(N,R1,alpha1,alpha2);

et1=(alpha1*R1+alpha2)./(sqrt(alpha*N)*(R1-1).^2.*phi(beta1));

rho2=-5:3;

et2=sqrt(N)*hf(rho2,alpha);

semilogy(rho,et,'r',rho1,et1,'b:',rho2,et2,'g:')

w=axis; w(3)=5; w(4)=2e5; axis(w)

function y=beta1f(N,R0,alpha1,alpha2)

% R0 is allowed to be a vector.

if alpha1==0

y=sign(R0-1).*sqrt(2*N*(R0-1-log(R0))/alpha2);

elseif alpha2==0

y=sign(R0-1).*sqrt(2*N*(log(R0)-(R0-1)./R0)/alpha1);

else

alpha=alpha1+alpha2;

a=(alpha/alpha1)*log((alpha1*R0+alpha2)/alpha);

y=sign(R0-1).*sqrt(2*N*(a-log(R0))/alpha2);

end

function y=capphi(x)

% capphi.m Computes the normal cdf Phi(x).

% For x<=-7 I use asymptotic approximation.

% Phi(x)=[erf(x/sqrt(2))+1]/2;

% The density is denoted phi(x).

% x is allowed to be a vector.

% Calls phi.

y=x; ii0=find(x<=-7); ii1=find(x>-7);

if length(ii0)>0

z=x(ii0); kk=26;

x1=-1./z; x2=x1;

for k=2:kk

x1=-x1*(2*k-3)./z.^2;

x2=x2+x1;

end

y(ii0)=phi(x(ii0)).*x2;

end

if length(ii1>0)

y(ii1)=(erf(x(ii1)/sqrt(2))+1)/2;

end

function y=h0f(rho,rhob)

% Computes the function H0(rho); 1/p1^(0) is

% approximated by 0.5*log(N)+H0(rho) for rho=O(1).

% rho is allowed to be a vector.

% Calls haf and h1f.

mb=floor((rhob^2-1+sqrt((rhob^2+1)^2+4*rhob^2))/4);

harhob=haf(rhob,mb);

l=length(rho); y=zeros(size(rho));

for k=1:l

if rho(k)<=rhob

y(k)=haf(rho(k),mb);

else

y(k)=harhob+quad8('h1f',rhob,rho(k));

end

end

function y=h1f(x)

% Computes H1(x)=Phi(x)/phi(x). x is allowed to be a vector.

% Calls capphi and phi.

y=capphi(x)./phi(x);

function y=haf(rho,mb)

% Determines the function Ha(rho) with the sum below running from 2 to mb.

sum=0; a=-1;

for k=2:mb

a=-a*(2*k-3);

term=a/((2*k-2)*rho^(2*k-2));

sum=sum+term;

end

y=-log(abs(rho))-sum;

function y=hf(rho,alpha)

% Determine the function H(rho,alpha) that is important for approximating the

% quasi-stationary distribution in the transition region.

% rho is allowed to be a vector.

% fzero finds a zero of the function hhf. h1(rho(k)) is initial approximation.

% rho(k) and alpha are parameters passed to hhf.

options=optimset('Display','off'); % Sets all other fields equal to [].

for k=1:length(rho)

y(k)=fzero('hhf',[0.001,h1f(rho(k))],options,rho(k),alpha);

end

function y=hhf(x,rho,alpha)

% This function is 0 for x=H(rho,alpha).

% Is used for determining H(rho,alpha) in hf.

% The main job done here is to integrate H1. Uses h1f.

if x>1/30,

aa=quad8('h1f',-1/x,rho);

else

aa=-log(30*x)+quad8('h1f',-1/30,rho);

end

y=aa-1-rho*x*sqrt(alpha);

function y=phi(x)

% phi.m Computes the normal density phi(x).

% x ix allowed to be a vector.

y=zeros(size(x));

y=exp(-x.^2/2)/sqrt(2*pi);

function p=p0f(N,R0,alpha1,alpha2)

% Determine the stationary distribution p0(N,R0,alpha1,alpha2) of the

% approximation X0 for the Verhulst model.

p=zeros(1,N);

p(1)=1; sum=1;

for k=2:N

p(k)=p(k-1)*(1-1/k)*(1-alpha1*(k-1)/N)*R0/(1+alpha2*k/N);

sum=sum+p(k);

if p(k)<2*eps*sum; break; end

end

% Determine where p is positive!

ispos=find(p>0); p=p(ispos);

p=p/sum;

function p=p1f(N,R0,alpha1,alpha2)

% Determine the stationary distribution p1(N,R0,alpha1,alpha2) of the

% approximation X1 for the Verhulst model.

p=zeros(1,N);

p(1)=1; sum=1;

for k=2:N

p(k)=p(k-1)*(1-alpha1*(k-1)/N)*R0/(1+alpha2*(k-1)/N);

sum=sum+p(k);

if p(k)<2*eps*sum; break; end

end

%Determine where p is positive!

ispos=find(p>0); p=p(ispos);

p=p/sum;

function qs1=qsf(N,R0,alpha1,alpha2)

% Determine the qs distribution using iteration.

% Each step of iteration starts from qs0 and determines qs1.

% The initial distribution qs0 equals p1.

% steps is number of iteration steps.

% Each term is called term(k); its numerator is num(k).

% The sum of the first n terms is called summa(n).

% Calls p1f, p0f.

p1=p1f(N,R0,alpha1,alpha2); p0=p0f(N,R0,alpha1,alpha2);

% The components of these two vectors are positive!

% Run iteration only up to an n-value where both p0 and p1

% are positive!

l=min(length(p1),length(p0));

p1=p1(1:l); p0=p0(1:l);

qs0=p1;

% Declare the vectors qs1, num, term, summa, ett:

qs1=zeros(size(qs0)); num=qs1; term=qs1; summa=qs1; ett=ones(size(qs0));

for steps=1:100

qsum=cumsum(qs0);

num=ett-[0 qsum(1:l-1)];

term=p1(1)*num./p1;

summa=cumsum(term);

qs1=p0.*summa/p0(1);

qs1=qs1/sum(qs1);

x1=max(qs1./qs0); x2=max(qs0./qs1);

if max(x1,x2)-1<eps*1e9, break, end

qs0=qs1;

end