ANALYTICAL AND COMPUTATIONAL

METHODS FOR MODELLING THE

LONG-TERM BEHAVIOUR OF EVANESCENT RANDOM PROCESSES

by

Phil Pollett

Department of Mathematics The University of Queensland

AN AUTO-CATALYTIC REACTION

Consider the following reaction scheme:

A →^{X} B,

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

A + X →^{k_{1}} 2X and 2X →^{k_{2}} B.

Let X(t) = 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 S = {0, 1, 2, …} and the transition rates, Q = (q_{xy}, x, y ∈ S), are given by

q_{xy} = \{ \begin{matrix}k_{1}ax & if y = x + 1, \\ k_{2} \binom{x}{2} & if y = x − 2, \\ 0 & otherwise . \end{matrix}

AN EPIDEMIC MODEL (Ridler-Rowe)

The state at time t:

X(t) - number of susceptibles Y (t) - number of infectives

State space:

S = {(x, y) : x, y = 0, 1, 2, …} Transition rates Q = (q_{ij}, i, j ∈ S):

if i = (x, y), then

q_{ij} = \{ \begin{matrix}βxy & if j = (x − 1, y + 1), \\ γy & if j = (x, y − 1), \\ α & if j = (x + 1, y), \\ 0 & otherwise . \end{matrix}

Clearly,

C = {(x, y) : x = 0, 1, …; y = 1, 2, …} is an irreducible transient class. Ridler-Rowe (1967): Q is regular (nonexplosive) and absorption occurs with probability 1.

QUASISTATIONARY DISTRIBUTIONS

Notation. Let X(t) be the state of the system at time t and suppose that X(t) ∈ S, where S, the state space, is a discrete set. Let Q = (q_{ij}, i, j ∈ S) be the matrix of transition rates (q-matrix), so that q_{ij} (≥ 0), for j ≠ i, represents the transition rate from state i to state j and q_{ii} = −q_{i}, where q_{i} = ∑_{j≠i} q_{ij} (< ∞) represents the transition rate out of state i.

Assumptions. Take 0 to be the sole absorbing state, that is, q_{0} = 0, and suppose that S = {0} ∪ C, where C = {1, 2, …} 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 q_{i0} > 0 for at least one i ∈ C.

Let p_{ij}(t) = Pr(X(t) = j|X(0) = i), t ≥ 0.

Definition. Let π = (π_{j}, j ∈ C) be a probability distribution over C and let

p_{j}(t) = _{i∈C}^{∑} π_{i}p_{ij}(t), j ∈ S, t ≥ 0.

Then, π is a quasistationary distribution if

\frac{p_{j}(t)}{∑_{i∈C} p_{i}(t)} = π_{j}, j ∈ C, t ≥ 0.

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,

lim_{t→∞} Pr(X(t) = j|X(0) = i, X(t) ∈ C) = π_{j},

no matter what the value of i(∈ C).

A CHARACTERIZATION

If π is a quasistationary distribution, then clearly

_{i∈C}^{∑} π_{i}p_{ij}(t) = g(t)π_{j}, j ∈ C, t > 0,

where g(t) = ∑_{i∈C} p_{i}(t). It is easy to show that g satisfies:

g(s + t) = g(s)g(t), s, t ≥ 0,

and 0 < g(t) < 1. Thus, g(t) = e^{−μt}, for some μ > 0. Since the converse is also true, we have the following result:

Proposition. A probability distribution, π = (π_{j}, j ∈ C), over C is a quasistationary distribution if and only if, for some μ > 0, π is a μ-invariant measure, that is

_{i∈C}^{∑} π_{i}p_{ij}(t) = e^{−μt}π_{j}, j ∈ C, t ≥ 0.  (1)

HOW DO WE DETERMINE π from Q?

Rewrite (1) as

_{i∈C: i≠j}^{∑} π_{i}p_{ij}(t) = ((1 − p_{jj}(t)) − (1 − e^{−μt})) π_{j}

and remember that q_{ij} is the right-hand derivative of p_{ij}(·) near 0. So, on dividing by t and letting t ↓ 0, we get (formally)

_{i∈C: i≠j}^{∑} π_{i}q_{ij} = (q_{j} − μ)π_{j}, j ∈ C,

or, equivalently,

_{i∈C}^{∑} π_{i}q_{ij} = −μπ_{j}, j ∈ C.  (2)

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

Proposition. If π is a quasistationary distribution then, for some μ > 0, π is a μ-invariant measure for Q.

IS THE CONVERSE TRUE?

Suppose that, for some μ > 0, π is a μ-invariant measure for Q, that is

_{i∈C}^{∑} π_{i}q_{ij} = −μπ_{j}, j ∈ C.

Is π a quasistationary distribution?

Sum this equation over j ∈ C: we get (formally)

_{i∈C}^{∑} π_{i}q_{i0} = − _{i∈C}^{∑} π_{i} _{j∈C}^{∑} q_{ij} = − _{j∈C}^{∑} _{i∈C}^{∑} π_{i}q_{ij}

= μ _{j∈C}^{∑} π_{j} = μ.

Theorem. Let π = (π_{j}, j ∈ C) be a probability distribution over C and suppose that π is a μ-invariant measure for Q. Then, μ ≥ ∑_{j∈C} π_{j}q_{j0} with equality if and only if π is a quasistationary distribution.

AN EXAMPLE

The birth-death-catastrophe process. Let S = {0, 1, 2, …} and suppose that

\begin{matrix}q_{i,i+1} = aρi, & i ≥ 0, \\ q_{i,i} = −ρi, & i ≥ 0, \\ q_{i,i−k} = ρib_{k}, & i ≥ 2, k = 1, 2 … i − 1, \\ q_{i,0} = ρi ∑_{k=i}^{∞} b_{k}, & i ≥ 1, \end{matrix}

where ρ > 0, a > 0, b_{i} > 0 for at least one value of i ≥ 1 and a+∑_{i=1}^{∞} b_{i} = 1. 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 b_{i}, i ≥ 1.

Clearly, 0 is an absorbing state and C = {1, 2, …} is an irreducible class.

So, does the process admit a quasistationary distribution?

Substitute the transition rates in (2):

−(ρ − μ)π_{1} + ∞ ; ∑ ; k=2 kρb_{k−1}π_{k} = 0,

and, for j ≥ 2,

(j −1)ρaπ_{j−1} −(jρ−μ)π_{j} + ∞ ; ∑ ; k=j+1 kρb_{k−j}π_{k} = 0.

If we try a solution of the form π_{j} = t^{j}, the first equation tells us that μ = −ρ(f′(t) − 1), where

f(s) = a + _{i∈C}^{∑} b_{i}s^{i+1} , |s| ≤ 1,

and, on substituting both of these quantities in the second equation, we find that f(t) = t. This latter equation has a unique solution, σ, on [0, 1]. Thus, by setting t = σ we can obtain a positive solution, π = (π_{j}, j ∈ C), to (2), such that ∑_{j∈C} π_{j} = 1 whenever σ < 1. Further, it is easy to show that ∑_{i∈C} π_{i}q_{i0} = μ.

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

D = a − _{i∈C}^{∑} ib_{i} = 1 − _{i∈C}^{∑}(i + 1)b_{i},

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, π = (π_{j}, j ∈ C), given by

π_{j} = (1 − σ)σ^{j−1}, j ∈ C,

where σ is the unique solution to f(t) = t on the interval [0, 1].

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, m = (m_{j}, j ∈ C), has positive entries; this is, of course, a μ-invariant measure for Q (unique up to constant multiples). Since S is finite, the quasistationary distribution, π = (π_{i}, i ∈ C), exists and is given by

π_{j} = m_{j}/ _{k∈C}^{∑} m_{k}, j ∈ C.

Infinite S : Truncate the restriction of Q to an n × n matrix, Q^{(n)}, and construct a sequence, {m^{(n)}}, 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

C_{N} = {(x, y) : x = 0, …, N − 1; y = 1, …, N}

and restrict Q to C_{N}. Use the transformation i = x + N(y − 1) to convert the restricted qmatrix into an n × n matrix, Q = (q_{ij}, i, j = 0, 1, …, n − 1), where n = N^{2}.

Evaluation of the eigenvectors of Q is not completely trivial. For example, if (as well shall assume) N = 100, that is n = 10^{4}, 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 Ax = λx, where A = Q^{T}.

Using an initial estimate of x, the basic Arnoldi method produces an m×m (upper-Hessenberg) matrix H_{m} and an n × m matrix V_{m} with

V_{m}^{T}AV_{m} = H_{m},

and such that if z_{m} is an eigenvector of H_{m}, then, for m large, V_{m}z_{m} is close to an eigenvector of A.

Solve for z_{m} using standard (dense-matrix) methods.

AN ITERATIVE ARNOLDI METHOD

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

Now solve

(H_{m} − λ̂I)u_{1} = z

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

v_{1} = V_{m}u_{1}/∥V_{m}u_{1}∥_{2}.

Continue until the residual, ∥Av_{1} − λ̂v_{1}∥_{2}, is sufficiently small.