MODELLING QUASISTATIONARITY IN EVANESCENT PROCESSES

by

Phil Pollett

Department of Mathematics

The University of Queensland

(Download the following image (as postscript)):

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

(Download the following image (as postscript)):

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

(Download the following image (as postscript)):

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.

(Download the following image (as postscript)):

Figure 4. Simulation of the Ridler-Rowe epidemic model

[Simulation |
java code]

(Download the following image (as postscript)):

Figure 5. Simple contour plot of the usual diffusion

approximation for the Ridler-Rowe epidemic model

(to be imagined as superimposed on Figure 4)

(Download the following image (as postscript)):

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

** 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.

(Download the following image (as postscript)):

Figure 7. Surface plot of the quasistationary distribution

for the Ridler-Rowe epidemic model

(Download the following image (as postscript)):

Figure 8. Contour plot of the quasistationary distribution

for the Ridler-Rowe epidemic model

(Download the following image (as postscript)):

Figure 8. Comparison between the quasistationary distribution

and the diffusion approximation (contour plot) for the

Ridler-Rowe epidemic model

(Download the following image (as postscript)):

Figure 10. Graph of the residual norm versus the

number of outer iterations for the

Iterative Arnoldi Method

(Download the following image (as postscript)):

Figure 11. Graph of the residual norm versus the

number of outer iterations for various

values of m, the number of inner iterations