Branching Chains

Generically, suppose that we have a system of particles that can generate or split into other particles of the same type. Here are some typical examples:

We assume that the lifetime of each particle is exponentially distributed with parameter

α \in (0, \infty)

, and at the end of its life, is replaced by a random number of new particles that we will refer to as children of the original particle. The number of children

N

of a particle has probability density function

f

N

. The particles act independently, so in addition to being identically distributed, the lifetimes and the number of children are independent from particle to particle. Finally, we assume that

f (1) = 0

, so that a particle cannot simply die and be replaced by a single new particle. Let

μ

and

σ^{2}

denote the mean and variance of the number of offspring of a single particle. So

μ = E (N) = \sum_{n = 0}^{\infty} n f (n), σ^{2} = var (N) = \sum_{n = 0}^{\infty} (n - μ)^{2} f (n)

We assume that

μ

is finite and so

σ^{2}

makes sense. In our study of discrete-time Markov chains, we studied branching chains in terms of generational time. Here we want to study the model in real time.

Let $X_{t}$ denote the number of particles at time $t \in [0, \infty)$ . Then $X = {X_{t} : t \in [0, \infty)}$ is a continuous-time Markov chain on $N$ , known as a branching chain. The exponential parameter function $λ$ and jump transition matrix $Q$ are given by

$λ (x) = α x$ for $x \in N$
$Q (x, x + k - 1) = f (k)$ for $x \in N_{+}$ and $k \in N$ .

Details:

That $X$ is a continuous-time Markov chain follows from the assumptions and the basic structure of continuous-time Markov chains. In turns out that the assumption that $μ < \infty$ implies that $X$ is regular, so that $τ_{n} \to \infty$ as $n \to \infty$ , where $τ_{n}$ is the time of the $n$ th jump for $n \in N_{+}$ .

Starting with $x$ particles, the time of the first state change is the minimum of $x$ independent variables, each exponentially distributed with parameter $α$ . As we know, this minimum is also exponentially distributed with parameter $α x$ .
Starting in state $x \in N_{+}$ , the next state will be $x + k - 1$ for $k \in N$ , if the particle dies and leaves $k$ children in her place. This happens with probability $f (k)$ .

Of course 0 is an absorbing state, since this state means extinction with no particles. (Note that

λ (0) = 0

and so by default,

Q (0, 0) = 1

.) So with a branching chain, there are essentially two types of behavior: population extinction or population explosion.

For the branching chain $X = {X_{t} : t \in [0, \infty)}$ one of the following events occurs with probability 1:

Extinction: $X_{t} = 0$ for some $t \in [0, \infty)$ and hence $X_{s} = 0$ for all $s \geq t$ .
Explosion: $X_{t} \to \infty$ as $t \to \infty$ .

Details:

If $f (0) > 0$ then all states lead to the absorbing state 0 and hence the set of positive staties $N_{+}$ is transient. With probability 1, the jump chain $Y$ visits a transient state only finitely many times, so with probability 1 either $Y_{n} = 0$ for some $n \in N$ or $Y_{n} \to \infty$ as $n \to \infty$ . If $f (0) = 0$ then $Y_{n}$ is strictly increasing in $n$ , since $f (1) = 0$ by assumption. Hence with probability 1, $Y_{n} \to \infty$ as $n \to \infty$ .

Without the assumption that

μ < \infty

, explosion can actually occur in finite time. On the other hand, the assumption that

f (1) = 0

is for convenience. Without this assumption,

X

would still be a continuous-time Markov chain, but as discussed in the Introduction, the exponential parameter function would be

λ (x) = α f (1) x

for

x \in N

and the jump transition matrix would be

Q (x, x + k - 1) = \frac{f (k)}{1 - f (1)}, x \in N_{+}, k \in {0, 2, 3, \dots}

Because all particles act identically and independently, the branching chain starting with

x \in N_{+}

particles is essentially

x

independent copies of the branching chain starting with 1 particle. In many ways, this is the fundamental insight into branching chains, and in particular, means that we can often condition on

X (0) = 1

Generator and Transition Matrices

As usual, we will let

P = {P_{t} : t \in [0, \infty)}

denote the semigroup of transition matrices of

X

, so that

P_{t} (x, y) = P (X_{t} = y ∣ X = x)

for

(x, y) \in N^{2}

. Similarly,

G

denotes the infinitesimal generator matrix of

X

The infinitesimal generator $G$ is given by $\begin{aligned} G (x, x) & = - α x, x \in N \\ G (x, x + k - 1) & = α x f (k), x \in N_{+}, k \in N \end{aligned}$

Details:

This follows immediately from the exponential parameter function and the jump transition matrix in [1].

The Kolmogorov backward equation is $\frac{d}{d t} P_{t} (x, y) = - α x P_{t} (x, x) + α x \sum_{k = 0}^{\infty} f (k) P_{t} (x + k - 1, y), (x, y) \in N^{2}$

Details:

The backward equation is $\frac{d}{d t} P_{t} = G P_{t}$ , so the result follows from [3].

Unlike some of our other continuous-time models, the jump chain

Y

governed by

Q

is not the discrete-time version of the model. That is,

Y

is not a discrete-time branching chain, since in discrete time, the index

n

represents the

n

th generation, whereas here it represent the

n

th time that a particle reproduces. However, there are lots of discrete-time branching chain embedded in the continuous-time chain.

Fix $t \in (0, \infty)$ and define $Z_{t} = {X_{n t} : n \in N}$ . Then $Z_{t}$ is a discrete-time branching chain with offspring probability density function $f_{t}$ given by $f_{t} (x) = P_{t} (1, x)$ for $x \in N$ .

Details:

In general, we know that sampling a (homogeneous) continuous-time Markov chain at multiples of a fixed $t \in (0, \infty)$ , results in a (homogeneous) discrete-time Markov chain. For $Z_{t}$ to be a branching chain, we just need to note that $P_{t} (x, y) = f_{t}^{* x} (y), (x, y) \in N^{2}$ where $f_{t}^{* x}$ is the convolution power of $f_{t}$ of order $x$ . This is a consequence of the fundamental fact that $X_{t}$ given $X_{0} = x$ has the same distribution as the sum of $x$ independent copies of $X_{t}$ given $X_{0} = 1$ . Recall that the PDF of a sum of independent variables is the convolution of the individual PDFs.

Probability Generating Functions

For $t \in [0, \infty)$ let $Φ_{t}$ denote the probability generating function of $X_{t}$ given $X_{0} = 1$ $Φ_{t} (r) = E (r^{X_{t}} ∣ X_{0} = 1) = \sum_{x = 0}^{\infty} r^{x} P_{t} (1, x)$ Let $Ψ$ denote the probability generating function of $N$ $Ψ (r) = E (r^{N}) = \sum_{n = 0}^{\infty} r^{n} f (n)$ The generating functions are defined (the series are absolutely convergent) at least for $r \in (- 1, 1]$ .

The collection of generating functions

Φ = {Φ_{t} : t \in [0, \infty)}

gives the same information as the collection of probability density functions

{P_{t} (1, \cdot) : t \in [0, \infty)}

. With the fundamental insight that the branching process starting with one particle determines the branching process in general,

Φ

actually determines the transition semigroup

P = {P_{t} : t \in [0, \infty)}

For $t \in [0, \infty)$ and $x \in N$ , the probability generating function of $X_{t}$ given $X_{0} = x$ is $Φ_{t}^{x}$ : $\sum_{y = 0}^{\infty} r^{y} P_{t} (x, y) = [Φ_{t} (r)]^{x}$

Details:

Again, given $X_{0} = x$ , the number of particles $X_{t}$ at time $t$ has the same distribution as the sum of $x$ independent copies of $X_{t}$ given $X_{0} = 1$ . Recall that the PGF of a sum of independent variables is the product of the PGFs of the variables.

Note that

Φ_{t}

is the generating function of the offspring distribution for the embedded discrete-time branching chain

Z_{t} = {X_{n t} : n \in N}

for

t \in (0, \infty)

. On the other hand,

Ψ

is the generating function of the offspring distribution for the continuous-time chain. So our main goal in this discussion is to see how

Φ

is built from

Ψ

. Because

P

is a semigroup under matrix multiplication, and because the particles act identically and independently,

Φ

is a semigroup under composition.

$Φ_{s + t} = Φ_{s} \circ Φ_{t}$ for $s, t \in [0, \infty)$ .

Details:

Using the semigroup property (the Chapman-Kolmogorov equations) and [7] we have $\begin{aligned} Φ_{s + t} (r) & = \sum_{y = 0}^{\infty} r^{y} P_{s + t} (1, y) = \sum_{y = 0}^{\infty} r^{y} \sum_{x = 0}^{\infty} P_{s} (1, x) P_{t} (x, y) = \sum_{x = 0}^{\infty} P_{s} (1, x) \sum_{y = 0}^{\infty} r^{y} P_{t} (x, y) \\ = \sum_{x = 0}^{\infty} P_{s} (1, x) [Φ_{t} (r)]^{x} = Φ_{s} [Φ_{t} (r)] \end{aligned}$

Note also that

Φ_{0} (r) = E (r^{X_{0}} ∣ X_{0} = 1) = r

for all

r \in R

. This also follows from the semigroup property:

Φ_{0} = Φ_{0} \circ Φ_{0}

. The fundamental relationship between the collection of generating functions

Φ

and the generating function

Ψ

is given in the following theorem:

The mapping $t \mapsto Φ_{t}$ satisfies the differential equation $\frac{d}{d t} Φ_{t} = α (Ψ \circ Φ_{t} - Φ_{t})$

Details:

Using the Kolmogorov backward equation in [4] we have $\frac{d}{d t} Φ_{t} (r) = \sum_{x = 0}^{\infty} r^{x} \frac{d}{d t} P_{t} (1, x) = \sum_{x = 0}^{\infty} r^{x} G P_{t} (1, x)$ Using the generator in [2], $G P_{t} (1, x) = \sum_{y = 0}^{\infty} G (1, y) P_{t} (y, x) = - α P_{t} (1, x) + \sum_{k = 0}^{\infty} α f (k) P_{t} (k, x), x \in N$ Substituting and using [7] gives $\begin{aligned} \frac{d}{d t} Φ_{t} (r) & = \sum_{x = 0}^{\infty} r^{x} [- α P_{t} (1, x) + \sum_{k = 0}^{\infty} α f (k) P_{t} (k, x)] = - α \sum_{x = 0}^{\infty} r^{x} P_{t} (1, x) + α \sum_{x = 0}^{\infty} \sum_{k = 0}^{\infty} r^{x} f (k) P_{t} (k, x) \\ = - α Φ_{t} (r) + α \sum_{k = 0}^{\infty} f (k) \sum_{x = 0}^{\infty} r^{x} P_{t} (k, x) = - α Ψ_{t} (r) + α \sum_{k = 0}^{\infty} f (k) [Φ_{t} (r)]^{k} = - α Φ_{t} (r) + α Ψ [Φ_{t} (r)] \end{aligned}$

This differential equation, along with the initial condition

Φ_{0} (r) = r

for all

r \in R

determines the collection of generating functions

Φ

. In fact, an implicit solution for

Φ_{t} (r)

is given by the integral equation

\int_{r}^{Φ_{t} (r)} \frac{1}{Ψ (u) - u} d u = α t

Another relationship is given in the following theorem. Here,

Φ_{t}^{'}

refers to the derivative of the generating function

Φ_{t}

with respect to its argument, of course (so

r

, not

t

For $t \in [0, \infty)$ , $Φ_{t}^{'} = \frac{Ψ \circ Φ_{t} - Φ_{t}}{Ψ - Φ_{0}}$

Details:

From the semigroup property, we have $Φ_{t + s} (r) = Φ_{t} [Φ_{s} (r)]$ for $s, t \in [0, \infty)$ . Differentiating with respect to $s$ and using the chain rule along with [9] gives $\frac{d}{d s} Φ_{t + s} (r) = Φ_{t}^{'} [Φ_{s} (r)] \frac{d}{d s} Φ_{s} (r) = Φ_{t}^{'} [Φ_{s} (r)] α [Ψ (Φ_{s} (r)) - Φ_{s} (r)]$ Evaluating at $s = 0$ and using the condition $Φ_{0} (r) = r$ we have $\frac{d}{d t} Φ_{t} (r) = Φ_{t}^{'} (r) α [Ψ (r) - r]$ Using the [9] once again gives $α [Ψ (Φ_{t} (r)) - Φ_{t} (r)] = Φ_{t}^{'} (r) α [Ψ (r) - r]$ Solving for $Φ_{t}^{'} (r)$ gives the result.

Moments

In this discussion, we wil study the mean and variance of the number of particles at time

t \in [0, \infty)

. Let

m_{t} = E (X_{t} ∣ X_{0} = 1), v_{t} = var (X_{t} ∣ X_{0} = 1), t \in [0, \infty)

so that

m_{t}

and

v_{t}

are the mean and variance, starting with a single particle. As always with a branching process, it suffices to consider a single particle:

For $t \in [0, \infty)$ and $x \in N$ ,

$E (X_{t} ∣ X_{0} = x) = x m_{t}$
$var (X_{t} ∣ X_{0} = x) = x v_{t}$

Details:

Once again, the distribution of $X_{t}$ given $X_{0} = x$ is the same as the distribution of the sum of $x$ independent copies of $X_{t}$ given $X_{0} = 1$ . Recall that the mean of a sum of variables is the sum of the individual means, and the variance of the sum of independent variables is the sum of the individual variances.

Recall also that

μ

and

σ^{2}

are the the mean and variance of the number of offspring of a particle. Here is the connection between the means:

$m_{t} = e^{α (μ - 1) t}$ for $t \in [0, \infty)$ .

If $μ < 1$ then $m_{t} \to 0$ as $t \to \infty$ . This is extinction in the mean.
If $μ > 1$ then $m_{t} \to \infty$ as $t \to \infty$ . This is explosion in the mean.
If $μ = 1$ then $m_{t} = 1$ for all $t \in [0, \infty)$ . This is stability in the mean.

Details:

From the proof of [10], $\frac{d}{d t} Φ_{t} (r) = α Φ_{t}^{'} (r) [Ψ (r) - r]$ Differentiating with respect to $r$ , interchanging the order of integration on the left, and using the product rule on the right gives $\frac{d}{d t} Φ_{t}^{'} (r) = α Φ_{t}^{''} (r) [Ψ (r) - r] + α Φ_{t}^{'} (r) [Ψ^{'} (r) - 1]$ Now let $r = 1$ and recall that $Φ (1) = 1$ . We get $\frac{d}{d t} Φ_{t}^{'} (1) = α Φ_{t}^{'} (1) [Ψ^{'} (1) - 1]$ From the basic theory of probability generating functions, $m_{t} = Φ_{t}^{'} (1)$ and similarly, $μ = Ψ^{'} (1)$ . Hence we have $\frac{d}{d t} m_{t} = α (μ - 1) m_{t}$ Of course we have the initial condition $m_{0} = 1$ .

This result is intuitively very appealing. As a function of time, the expected number of particles either grows or decays exponentially, depending on whether the expected number of offspring of a particle is greater or less than one. The connection between the variances is more complicated. We assume that

σ^{2} < \infty

If $μ \neq 1$ then $v_{t} = [\frac{σ^{2}}{μ - 1} + (μ - 1)] [e^{2 α (μ - 1) t} - e^{α (μ - 1) t}], t \in [0, \infty)$ If $μ = 1$ then $v_{t} = α σ^{2} t$ .

If $μ < 1$ then $v_{t} \to 0$ as $t \to \infty$
If $μ \geq 1$ then $v_{t} \to \infty$ as $t \to \infty$

Details:

Probability generating functions are naturally connected to factorial moments, so it's best to work with these. Thus, let $w_{t} = E [X_{t} (X_{t} - 1) ∣ X_{0} = 1]$ for $t \in [0, \infty)$ and let $δ = E [N (N - 1)]$ . These are the factorial moments of order 2. In the proof of the last theorem we showed that $\frac{d}{d t} Φ_{t}^{'} (r) = α Φ_{t}^{''} (r) [Ψ (r) - r] + α Φ_{t}^{'} (r) [Ψ^{'} (r) - 1]$ Differentiating with respect to $r$ again gives $\frac{d}{d t} Φ_{t}^{''} (r) = α Φ_{t}^{'''} (r) [Ψ (r) - r] + 2 α Φ_{t}^{''} (r) [Ψ^{'} (r) - 1] + α Φ_{t}^{'} (r) Ψ^{''} (r)$ Now substitute $r = 1$ . Recall that $Φ_{t}^{''} (1) = w_{t}$ , $Φ_{t}^{'} (1) = m_{t} = e^{α (μ - 1) t}$ , $Ψ^{''} (1) = δ$ , $Ψ^{'} (1) = μ$ , and $Ψ (1) = 1$ . We get the differential equation $\frac{d}{d t} w_{t} = 2 α (μ - 1) w_{t} + α δ e^{α (μ - 1) t}$ with the initial condition $w_{0} = 0$ .

Suppose that $μ \neq 1$ . Then using standard methods for a linear, first order differential equations with constant coefficients and an exponential forcing function, the solution is $w_{t} = \frac{δ}{μ - 1} [e^{2 α (μ - 1)} - e^{α (μ - 1) t}]$ But $δ = σ^{2} + μ^{2} - μ$ , and similarly $w_{t} = v_{t} + m_{t}^{2} - m_{t}$ with $m_{t} = e^{α (μ - 1) t}$ . Substitution and some algebra then gives the result.

Suppose now that $μ = 1$ . Then also $m_{t} = 1$ for all $t \in [0, \infty)$ and so $δ = σ^{2}$ and $v_{t} = w_{t}$ . The differential equation above reduces simply to $\frac{d}{d t} v_{t} = α σ^{2}$ with initial condition $v_{0} = 0$ so trivially $v_{t} = α σ^{2} t$ . Finally, in the context of part (b), note that if $μ = 1$ we must have $σ^{2} > 0$ since we have assumed that $f (1) = 0$ .

μ < 1

so that

m_{t} \to 0

t \to \infty

and we have extinction in the mean, then

v_{t} \to 0

t \to \infty

also. If

μ > 1

so that

m_{t} \to \infty

t \to \infty

and we have explosion in the mean, then

v_{t} \to \infty

t \to \infty

also. We would expect these results. On the other hand, if

μ = 1

so that

m_{t} = 1

for all

t \in [0, \infty)

and we have stability in the mean, then

v_{t}

grows linearly in

t

. This gives some insight into what to expect next when we consider the probability of extinction.

The Probability of Extinction

As shown in [2], there are two types of behavior for a branching process, either population extinction or population explosion. In this discussion, we study the extinction probability, starting as usual with a single particle:

q = P (X_{t} = 0 for some t \in (0, \infty) ∣ X_{0} = 1) = lim_{t \to \infty} P (X_{t} = 0 ∣ X_{0} = 1)

Need we say it? The extinction probability starting with an arbitrary number of particles is easy.

For $x \in N$ , $P (X_{t} = 0 for some t \in (0, \infty) ∣ X_{0} = x) = lim_{t \to \infty} P (X_{t} = 0 ∣ X_{0} = x) = q^{x}$

Details:

Given $X_{0} = x$ , extinction has occurred by time $t$ if and only if extinction has occurred by time $t$ for each of the $x$ independent branching chains formed from the descendents of the $x$ initial particles.

We can easily relate extinction for the continuous-time branching chain

X

to extinction for any of the embedded discrete-time branching chains.

If extinction occurs for $X$ then extinction occurs for $Z_{t}$ for every $t \in (0, \infty)$ . Conversely if extinction occurs for $Z_{t}$ for some $t \in (0, \infty)$ then extinction occurs for $Z_{t}$ for every $t \in (0, \infty)$ and extinction occurs for $X$ . Hence $q$ is the minimum solution in $(0, 1]$ of the equation $Φ_{t} (r) = r$ for every $t \in (0, \infty)$ .

Details:

The statements about the extinction event follow immediately from the fact that $0$ is absorbing, so that if $X_{t} = 0$ for some $t \in (0, \infty)$ then $X_{s} = 0$ for every $s \in [t, \infty)$ . The result for the extinction probability $q$ follows from the theory of discrete-time branching chains.

So whether or not extinction is certain depends on the critical parameter

μ

The extinction probability $q$ and the mean of the offspring distribution $μ$ are related as follows:

If $μ \leq 1$ then $q = 1$ , so extinction is certain.
If $μ > 1$ then $0 < q < 1$ , so there is a positive probability of extinction and a positive probability of explosion.

Details:

These results follow from the corresponding results for discrete-time branching chains. Fix $t \in (0, \infty)$ and recall that $m_{t}$ is the mean of the offspring distribution for the discrete-time chain $Z_{t} = {X_{n t} : n \in N}$ . From [12],

If $μ \leq 1$ then $m_{t} \leq 1$ .
If $μ > 1$ then $m_{t} > 1$ .

It would be nice to have an equation for

q

in terms of the offspring probability generating function

Ψ

. This is also easy

The probability of extinction $q$ is the minimum solution in $(0, 1]$ of the equation $Ψ (r) = r$ .

Details:

From [16], $Φ_{t} (q) = 1$ for every $t \in (0, \infty)$ . Substituting $r = q$ in exercose [9], we have $\frac{d}{d t} Φ_{t} (q) = 0$ and hence $Ψ (q) = q$ . As in the theory of discrete branching chains, the equation $Ψ (r) = r$ has only the solution 1 in (0, 1] if $μ = Ψ^{'} (1) \leq 1$ or there are two solutions $q \in (0, 1)$ and $1$ if $μ > 0$ . In both cases, $q$ is the smaller solution.

Special Models

We now turn our attention to a number of special branching chains that are important in applications or lead to interesting insights. We will use the notation established above, so that

α

is the parameter of the exponential lifetime of a particle,

Q

is the transition matrix of the jump chain,

G

is the infinitesimal generator matrix, and

P_{t}

is the transition matrix at time

t \in [0, \infty)

. Similarly,

m_{t} = E (X_{t} ∣ X_{0} = x)

v_{t} = var (X_{t} ∣ X_{0} = x)

, and

Φ_{t}

are the mean, variance, and generating function of the number of particles at time

t \in [0, \infty)

, starting with a single particle. As always, be sure to try these exercises yourself before looking at the proofs and solutions.

The Pure Death Branching Chain

First we consider the branching chain in which each particle simply dies without offspring. Sadly for these particles, extinction is inevitable, but this case is still a good place to start because the analysis is simple and lead to explicit formulas. Thus, suppose that

X = {X_{t} : t \in [0, \infty)}

is a branching process with lifetime parameter

α \in (0, \infty)

and offspring probability density function

f

with

f (0) = 1

The transition matrix of the jump chain and the generator matrix are given by

$Q (0, 0) = 1$ and $Q (x, x - 1) = 1$ for $x \in N_{+}$
$G (x, x) = - α x$ for $x \in N$ and $G (x, x - 1) = α x$ for $x \in N_{+}$

Let $t \in [0, \infty)$ . Then

$m_{t} = e^{- α t}$
$v_{t} = e^{- α t} - e^{- 2 α t}$
$Φ_{t} (r) = 1 - (1 - r) e^{- α t}$ for $r \in R$
Given $X_{0} = x$ the distribution of $X_{t}$ is binomial with trial parameter $x$ and success parameter $e^{- α t}$ . $P_{t} (x, y) = (\binom{x}{y}) e^{- α t y} (1 - e^{- α t})^{x - y}, x \in N, y \in {0, 1, \dots, x}$

Details:

All of these results follow from the general methods above, with $μ = σ = 0$ and $Ψ (r) = 1$ for $r \in R$ . But it's helpful to give direct proofs. Given $X_{0} = 1$ , let $τ$ be the time until the first transition, which is simply the lifetime of the particle. So $τ$ has the exponential distribution with parameter $α$ . For $t \in [0, \infty)$ , $X_{t}$ is an indicator random variable (taking just values 0 and 1) with $P (X_{t} = 1 ∣ X_{0} = 1) = P (τ > t ∣ X_{0} = 1) = e^{- α t}$ Part (a), (b), and (c) are standard results for an indicator variable. For part (d), given $X_{0} = x$ , each of the $x$ particles, independently, is still alive at time $t$ with probability $e^{- α t}$ . Hence the number of particles still alive has the binomial distribution with parameters $x$ and $e^{- α t}$ .

In particular, note that

P_{t} (x, 0) = (1 - e^{- α t})^{x} \to 1

t \to \infty

. that is, the probability of extinction by time

t

increases to 1 exponentially fast. Since we have an explicit formula for the transition matrices, we can find an explicit formula for the potential matrices as well. The result uses the beta function

B

For $β \in (0, \infty)$ the potential matrix $U_{β}$ is given by $U_{β} (x, y) = \frac{1}{α} (\binom{x}{y}) B (y + β / α, x - y + 1), x \in N, y \in {0, 1, \dots, x}$ For $β = 0$ , the potential matrix $U$ is given by

$U (x, 0) = \infty$ for $x \in N$
$U (x, y) = 1 / α y$ for $x, y \in N_{+}$ and $x \leq y$ .

Details:

Suppose that $β > 0$ and that $x, y \in N$ with $x \leq y$ . By definition $U_{β} (x, y) = \int_{0}^{\infty} e^{- β t} P_{t} (x, y) d t = \int_{0}^{\infty} e^{- β t} (\binom{x}{y}) e^{- α t y} (1 - e^{- α t})^{x - y} d t$ Substitute $u = e^{- α t}$ so that $d u = - α e^{- α t} d t$ or equivalently $d t = - d u / α u$ . After some algebra, the result is $U_{β} (x, y) = \frac{1}{α} (\binom{x}{y}) \int_{0}^{1} u^{y + β / α - 1} (1 - u)^{x - y} d u$ By definition, the last integral is $B (y + β / α, x - y + 1)$ .

For $x \in N$ , $U (x, 0) = \int_{0}^{\infty} (1 - e^{- α t})^{x} d d = \infty$
For $x, y \in N_{+}$ with $x \leq y$ , the derivation above and properties of the beta function give $U (x, y) = \frac{1}{α} (\binom{x}{y}) B (y, x - y + 1) = \frac{1}{α} (\binom{x}{y}) \frac{(y - 1)! (x - y)!}{x!} = \frac{1}{α y}$

We could argue the results for the potential

U

directly. Recall that

U (x, y)

is the expected time spent in state

y

starting in state

x

. Since 0 is absorbing and all states lead to 0,

U (x, 0) = \infty

for

x \in N

. If

x, y \in N_{+}

and

x \leq y

, then

x

leads to

y

with probability 1. Once in state

y

the time spent in

y

has an exponential distribution with parameter

λ (y) = α y

, and so the mean is

1 / α y

. Of course, when the chain leaves

y

, it never returns.

Recall that

β U_{β}

is a transition probability matrix for

β > 0

, and in fact

β U_{β} (x, \cdot)

is the probability density function of

X_{T}

given

X_{0} = x

where

T

is independent of

X

has the exponential distribution with parameter

β

. For the next result, recall the ascending power notation

a^{[k]} = a (a + 1) \dots (a + k - 1), a \in R, k \in N

For $β > 0$ and $x \in N_{+}$ , the function $β U_{β} (x, \cdot)$ is the beta-binomial probability density function with parameters $x$ , $β / α$ , and 1. $β U_{β} (x, y) = (\binom{x}{y}) \frac{(β / α)^{[y]} 1^{[x - y]}}{(1 + β / α)^{[x]}}, x \in N, y \in {0, 1, \dots x}$

Details:

From [20] and properties of the beta function. $β U_{β} (x, y) = \frac{β}{α} (\binom{x}{y}) B (y + β / α, x - y + 1), x \in N, y \in {0, 1, \dots, x}$ But from properties of the beta function, $B (y + β / α, x - y + 1) = B (β / α, 1) \frac{(β / α)^{[y]} 1^{[x - y]}}{(1 + β / α)^{[x]}} = \frac{α}{β} \frac{(β / α)^{[y]} 1^{[x - y]}}{(1 + β / α)^{[x]}}$ Substituting gives the result

The Yule Process

Next we consider the pure birth branching chain in which each particle, at the end of its life, is replaced by 2 new particles. Equivalently, we can think of particles that never die, but each particle gives birth to a new particle at a constant rate. This chain could serve as the model for an unconstrained nuclear reaction, and is known as the Yule process, named for George Yule. So specifically, let

X = {X_{t} : t \in [0, \infty)}

be the branching chain with exponential parameter

α \in (0, \infty)

and offspring probability density function given by

f (2) = 1

. Explosion is inevitable, starting with at least one particle, but other properties of the Yule process are interesting. in particular, there are fascinating parallels with the pure death branching chain in . Since 0 is an isolated, absorbing state, we will sometimes restrict our attention to positive states.

The transition matrix of the jump chain and the generator matrix are given by

$Q (0, 0) = 1$ and $Q (x, x + 1) = 1$ for $x \in N_{+}$
$G (x, x) = - α x$ for $x \in N$ and $G (x, x + 1) = α x$ for $x \in N_{+}$

Since the Yule process is a pure birth process and the birth rate in state

x \in N

α x

, the process is also called the linear birth chain. As with the pure death process, we can give the distribution of

X_{t}

specifically.

Let $t \in [0, \infty)$ . Then

$m_{t} = e^{α t}$
$v_{t} = e^{2 α t} - e^{α t}$
$Φ_{t} (r) = \frac{r e^{- α t}}{1 - r + r e^{- α t}}$ for $| r | < \frac{1}{1 - e^{- α t}}$
Given $X_{0} = x$ , $X_{t}$ has the negative binomial distribution on $N_{+}$ with stopping parameter $x$ and success parameter $e^{- α t}$ . $P_{t} (x, y) = (\binom{y - 1}{x - 1}) e^{- x α t} (1 - e^{- α t})^{y - x}, x \in N_{+}, y \in {x, x + 1, \dots}$

Details:

Proof from the general results: Parts (a) and (b) follow from the general moment results above, with $μ = 2$ and $σ^{2} = 0$ . For part (c), note that $Ψ (r) = r^{2}$ for $r \in R$ , so the integral equation for $Φ_{t}$ is $\int_{r}^{Φ_{t} (r)} \frac{1}{u^{2} - u} = α t$ From partial fractions, $\frac{1}{u^{2} - u} = \frac{1}{u - 1} - \frac{1}{u}$ , so the result follows by standard integration and algebra. We recognize $Φ_{t}$ as the probability generating function of the geometric distribution on $N_{+}$ with success parameter $e^{- α t}$ , so for part (d) we use our standard argument. Given $X_{0} = x \in N_{+}$ , $X_{t}$ has the same distribution as the sum of $x$ independent copies of $X_{t}$ given $X_{0} = 1$ , and so this is the distribution of the sum of $x$ independent variables each with the geometric distribution on $N_{+}$ with parameter $e^{- α t}$ . But this is the negative binomial distribution on $N_{+}$ with parameters $x$ and $e^{- α t}$ .

Direct proof: As usual, let $τ_{0} = 0$ and let $τ_{n}$ denote the time of the $n$ th transition (birth) for $n \in N_{+}$ . Given $X_{0} = 1$ , the population is $n$ at time $τ_{n - 1}$ . So the random interval $τ_{n} - τ_{n - 1}$ (the time until the next birth) has the exponential distribution with parameter $α n$ and these intervals are independent as $n$ varies. From a result in the section on the exponential distribution, it follows that $τ_{n} = \sum_{k = 1}^{n} (τ_{k} - τ_{k - 1})$ has distribution function given by $P (τ_{n} \leq t ∣ X_{0} = 1) = (1 - e^{- α t})^{n}, t \in [0, \infty)$ Curiously, this is also the distribution function of the maximum of $n$ independent variables, each with the exponential distribution with rate $α$ . Hence $P (X_{t} \geq n ∣ X_{0} = 1) = P (τ_{n - 1} \leq t ∣ X_{0} = 1) = (1 - e^{- α t})^{n - 1}, n \in N_{+}$ and therefore $P (X_{t} = n ∣ X_{0} = 1) = P (X_{t} \geq n ∣ X_{0} = 1) - P (X_{t} \geq n + 1 ∣ X_{0} = 1) = (1 - e^{- α t})^{n - 1} e^{- α t}, n \in N_{+}$ So given $X_{0} = 1$ , $X_{t}$ has the geometric distribution with parameter $e^{- α t}$ . The other results then follow easily.

Recall that the negative binomial distribution with parameters

k \in N_{+}

and

p \in (0, 1)

governs the trial number of the

k

th success in a sequence of Bernoulli Trials with success parameter

p

. So the occurrence of this distribution in the Yule process suggests such an interpretation. However this interpretation is not nearly as obvious as with the binomial distribution in the pure death branching chain. Next we give the potential matrices.

For $β \in [0, \infty)$ the potential matrix $U_{β}$ is given by $U_{β} (x, y) = \frac{1}{α} (\binom{y - 1}{x - 1}) B (x + β / α, y - x + 1), x \in N_{+}, y \in {x, x + 1, \dots}$ If $β > 0$ , the function $β U_{β} (x, \cdot)$ is the beta-negative binomial probability density function with parameters $x$ , $β / α$ , and 1: $β U_{β} (x, y) = (\binom{y - 1}{x - 1}) \frac{(β / α)^{[x]} 1^{[y - x]}}{(1 + β / α)^{[y]}}, x \in N, y \in {x, x + 1, \dots}$

Details:

The proof is very similar to the one in [23]. Suppose that $β \geq 0$ and that $x, y \in N_{+}$ with $y \geq x$ . By definition $U_{β} (x, y) = \int_{0}^{\infty} e^{- β t} P_{t} (x, y) d t = \int_{0}^{\infty} e^{- β t} (\binom{y - 1}{x - 1}) e^{- α t x} (1 - e^{- α t})^{y - x} d t$ Substitute $u = e^{- α t}$ so that $d u = - α e^{- α t} d t$ or equivalently $d t = - d u / α u$ . After some algebra, the result is $U_{β} (x, y) = \frac{1}{α} (\binom{y - 1}{x - 1}) \int_{0}^{1} u^{x + β / α - 1} (1 - u)^{y - x} d u$ By definition, the last integral is $B (x + β / α, y - x + 1)$ .

If we think of the Yule process in terms of particles that never die, but each particle gives birth to a new particle at rate

α

, then we can study the age of the particles at a given time. As usual, we can start with a single, new particle at time 0. So to set up the notation, let

X = {X_{t} : t \in [0, \infty)}

be the Yule branching chain with birth rate

α \in (0, \infty)

, and assume that

X_{0} = 1

. Let

τ_{0} = 0

and for

n \in N_{+}

, let

τ_{n}

denote the time of the

n

th transition (birth).

For $t \in [0, \infty)$ , let $A_{t}$ denote the total age of the particles at time $t$ . Then $A_{t} = \sum_{n = 0}^{X_{t} - 1} (t - τ_{n}), t \in [0, \infty)$ The random process $A = {A_{t} : t \in [0, \infty)}$ is the age process.

Details:

Note that there have been $X_{t} - 1$ births in the interval $[0, t]$ . For $n \in {0, 1, \dots, X_{t} - 1}$ , the age at time $t$ of the particle born at time $τ_{n}$ is $t - τ_{n}$ .

Again, let $A = {A_{t} : t \in [0, \infty)}$ be the age process for the Yule chain starting with a single particle. Then $A_{t} = \int_{0}^{t} X_{s} d s, t \in [0, \infty)$

Details:

Suppose that $X_{t} = k + 1$ where $k \in N$ , so that $τ_{k} \leq t < τ_{k + 1}$ . Note that $X_{s} = n$ for $τ_{n - 1} \leq s < τ_{n}$ and $n \in {1, 2, \dots, k}$ , while $X_{s} = k + 1$ for $τ_{k} \leq s \leq t$ . Hence $\int_{0}^{t} X_{s} d s = \sum_{n = 1}^{k} n (τ_{n} - τ_{n - 1}) + (k + 1) (t - τ_{k}) = (k + 1) t - \sum_{n = 0}^{k} τ_{n}$ From [25], $A_{t} = \sum_{n = 0}^{k} (t - τ_{n}) = (k + 1) t - \sum_{n = 0}^{k} τ_{n}$

With the last representation, we can easily find the expected total age at time

t

Again, let $A = {A_{t} : t \in [0, \infty)}$ be the age process for the Yule chain starting with a single particle. Then $E (A_{t}) = \frac{e^{α t} - 1}{α}, t \in [0, \infty)$

Details:

We can interchange the expected value and the integral by Fubini's theorem. So using [23], $E (A_{t}) = E (\int_{0}^{t} X_{s} d s) = \int_{0}^{t} E (X_{s}) d s = \int_{0}^{t} e^{α s} d s = \frac{e^{α t} - 1}{α}$

The General Birth-Death Branching Chain

Next we consider the continuous-time branching chain in which each particle, at the end of its life, leaves either no children or two children. At each transition, the number of particles either increases by 1 or decreases by 1, and so such a branching chain is also a continuous-time birth-death chain. Specifically, let

X = {X_{t} : t \in [0, \infty)}

be a continuous-time branching chain with lifetime parameter

α \in (0, \infty)

and offspring probability density function

f

given by

f (0) = 1 - p

f (2) = p

, where

p \in [0, 1]

. When

p = 0

we have the pure death chain in [18], and when

p = 1

we have the Yule process in [22]. We have already studied these, so the interesting case is when

p \in (0, 1)

so that both extinction and explosion are possible.

The transition matrix of the jump chain and the generator matrix are given by

$Q (0, 0) = 1$ , and $Q (x, x - 1) = 1 - p$ , $Q (x, x + 1) = p$ for $x \in N_{+}$
$G (x, x) = - α x$ for $x \in N$ , and $G (x, x - 1) = α (1 - p) x$ , $G (x, x + 1) = α p x$ for $x \in N_{+}$

As mentioned earlier,

X

is also a continuous-time birth-death chain on

N

, with 0 absorbing. In state

x \in N_{+}

, the birth rate is

α p x

and the death rate is

α (1 - p) x

. The moment functions are given next.

For $t \in [0, \infty)$ ,

$m_{t} = e^{α (2 p - 1) t}$
If $p \neq \frac{1}{2}$ , $v_{t} = [\frac{4 p (1 - p)}{2 p - 1} + (2 p - 1)] [e^{2 α (2 p - 1) t} - e^{α (2 p - 1) t}]$ If $p = \frac{1}{2}$ , $v_{t} = 4 α p (1 - p) t$ .

Details:

These results follow from the general formulas above for $m_{t}$ and $v_{t}$ , since $μ = 2 p$ and $σ^{2} = 4 p (1 - p)$ .

The next result gives the generating function of the offspring distribution and the extinction probability.

For the birth-death branching chain,

$Ψ (r) = p r^{2} + (1 - p)$ for $r \in R$ .
$q = 1$ if $0 < p \leq \frac{1}{2}$ and $q = \frac{1 - p}{p}$ if $\frac{1}{2} < p < 1$ .

Details:

By definition, $Ψ (r) = (1 - p) r^{0} + p r^{2}$ .
The equation $Ψ (r) = r$ factors into $(r - 1) [p r - (1 - p)] = 0$ . The extinction probability $q$ is the smaller solution in $(0, 1]$ .

For $t \in [0, \infty)$ , the generating function $Φ_{t}$ is given by $\begin{aligned} Φ_{t} (r) & = \frac{p r - (1 - p) + (1 - p) (1 - r) e^{α (2 p - 1) t}}{p r - (1 - p) + p (1 - r) e^{α (2 p - 1) t}}, if p \neq 1 / 2 \\ Φ_{t} (r) & = \frac{2 r + (1 - r) α t}{2 + (1 - r) α t}, if p = \frac{1}{2} \end{aligned}$

Details:

The integral equation for $Φ_{t}$ is $\int_{r}^{Φ_{t} (r)} \frac{d u}{p u^{2} + (1 - p) - u} = α t$ The denominator in the integral factors into $(u - 1) [p u - (1 - p)]$ . If $p \neq \frac{1}{2}$ , use partial fractions, standard integration, and some algebra. If $p = \frac{1}{2}$ the factoring is $\frac{1}{2} (u - 1)^{2}$ and partial fractions is not necessary. Again, use standard integration and algebra.

23. Continuous-Time Branching Chains

Basic Theory

Introduction

Generator and Transition Matrices

Probability Generating Functions

Moments

The Probability of Extinction

Special Models

The Pure Death Branching Chain

The Yule Process

The General Birth-Death Branching Chain