MODELLING QUASISTATIONARITY IN EVANESCENT PROCESSES

by

Phil Pollett

Department of Mathematics
The University of Queensland

Figure 1. Quasistationarity of an explosive population process
[Simulation | java code]

QUASISTATIONARITY

Processes appear to be stationary over any reasonable time scale, yet they are evanescent.

Examples:

• Absorption (eg population extinction)
• Positive drift (eg population explosion)
• Tunnelling between stable states (eg bistability in communication networks)

AN AUTO-CATALYTIC REACTION

Consider the following reaction scheme:

where X is a catalyst. Suppose that there are two stages, namely

Let number of X molecules at time t.

Suppose that the concentration of A is held constant; let a be the number of molecules of A. The state space is and the transition rates, , are given by

Figure 2. Auto-catalytic reaction simulation
[Simulation | java code]

AN EPIDEMIC MODEL
(Ridler-Rowe)

The state at time t:

- number of susceptibles

- number of infectives

State space:

Transition rates :

if , then

Figure 3. State transition diagram for
Ridler-Rowe epidemic model.

Clearly,

is an irreducible transient class.

Ridler-Rowe (1967): Q is regular (nonexplosive) and absorption occurs with probability 1.

Figure 4. Simulation of the Ridler-Rowe epidemic model
[Simulation | java code]

Figure 5. Simple contour plot of the usual diffusion
approximation for the Ridler-Rowe epidemic model
(to be imagined as superimposed on Figure 4)

Figure 6. Surface plot of the usual diffusion
approximation for the Ridler-Rowe epidemic model

QUASISTATIONARY DISTRIBUTIONS

Notation. Let be the state of the system at time t and suppose that , where S, the state space, is a discrete set. Let be the matrix of transition rates (q-matrix), so that , for , represents the transition rate from state i to state j and , where represents the transition rate out of state i.

Assumptions. Take 0 to be the sole absorbing state, that is, , and suppose that , where is an irreducible class. In order that there be a positive probability of ever reaching 0, given that the process starts in C, suppose that for at least one .

Let , .

Definition. Let be a probability distribution over C and let

Then, is a quasistationary distribution if

That is, if the process has as its initial distribution, then is a quasistationary distribution if the state probabilities at time t, conditional on the process being in C at t, are the same for all t.

If is a quasistationary distribution then, under certain conditions,

no matter what the value of .

A CHARACTERIZATION

If is a quasistationary distribution, then clearly

where . It is easy to show that g satisfies:

and . Thus, , for some . Since the converse is also true, we have the following result:

Proposition. A probability distribution, , over C is a quasistationary distribution if and only if, for some , is a -invariant measure, that is

HOW DO WE DETERMINE from Q?

Rewrite (1) as

and remember that is the right-hand derivative of near 0. So, on dividing by t and letting , we get (formally)

or, equivalently,

Accordingly, we shall say that is a -invariant measure for Q whenever (2) holds.

Proposition. If is a quasistationary distribution then, for some , is a -invariant measure for Q.

IS THE CONVERSE TRUE?

Suppose that, for some , is a -invariant measure for Q, that is

Is a quasistationary distribution?

Sum this equation over : we get (formally)

Theorem. Let be a probability distribution over C and suppose that is a -invariant measure for Q. Then, with equality if and only if is a quasistationary distribution.

AN EXAMPLE

The birth-death-catastrophe process. Let and suppose that

where , a>0, for at least one value of and . Thus, jumps occur at a constant ``per capita'' rate and, at a jump time, a birth occurs with probability a, or otherwise a catastrophe occurs, the size of which is determined by the probabilities , .

Clearly, 0 is an absorbing state and is an irreducible class.

So, does the process admit a quasistationary distribution?

Substitute the transition rates in (2):

and, for ,

If we try a solution of the form , the first equation tells us that , where

and, on substituting both of these quantities in the second equation, we find that . This latter equation has a unique solution, , on . Thus, by setting we can obtain a positive solution, , to (2), such that whenever . Further, it is easy to show that .

The condition is satisfied only in the subcritical case, that is, when the drift, D, given by

is strictly negative; incidentally, this also guarantees that absorption occurs with probability 1.

Proposition. The subcritical birth-death-catastrophe process has a geometric quasistationary distribution , , given by

where is the unique solution to on the interval .

COMPUTATIONAL METHODS

Finite S : Mandl (1960) showed that the restriction of Q to C has eigenvalues with negative real parts and the one with maximal real part (called above) is real and has multiplicity 1, and, the corresponding left eigenvector, , has positive entries; this is, of course, a -invariant measure for Q (unique up to constant multiples). Since S is finite, the quasistationary distribution, , exists and is given by

Infinite S : Truncate the restriction of Q to an matrix, , and construct a sequence, , of eigenvectors and hope that this converges to a -invariant measure, m, for Q, et cetera.

HOW SHOULD WE EVALUATE m?

Consider the Ridler-Rowe epidemic model.

First truncate C to

and restrict Q to . Use the transformation to convert the restricted q-matrix into an matrix, , where .

Evaluation of the eigenvectors of Q is not completely trivial. For example, if (as well shall assume) N=100, that is , Q needs 400 Mbytes of storage.

However, for N large, this matrix is large and sparse. Does this help?

THE ARNOLDI METHOD

We need to solve , where .

Using an initial estimate of x, the basic Arnoldi method produces an (upper-Hessenberg) matrix and an matrix with

and such that if is an eigenvector of , then, for m large, is close to an eigenvector of A.

Solve for using standard (dense-matrix) methods.

AN ITERATIVE ARNOLDI METHOD

Take m small (we found that m=20 worked best). Then, using an initial estimate, , of the eigenvector x, apply the basic Arnoldi method (to obtain and ) and set to be the dominant eigenvalue of if this is real, or set equal to zero otherwise.

Now solve

with z chosen at random and repeat the procedure with a new initial estimate, given by

Continue until the residual, , is sufficiently small.

Figure 7. Surface plot of the quasistationary distribution
for the Ridler-Rowe epidemic model

Figure 8. Contour plot of the quasistationary distribution
for the Ridler-Rowe epidemic model

Figure 8. Comparison between the quasistationary distribution
and the diffusion approximation (contour plot) for the
Ridler-Rowe epidemic model