2. MODEL AND ANALYSIS

Consider a linear birth–death process with immigration and emigration. The state of the system is the population size. When the state of the system is n,n = 0,1,…, the immigration rate is ν, the emigration rate is θ, the birth rate is nλ, and the death rate is nμ.

Let X(t) be the population size at time t. Then, {X(t); t ≥ 0} is a continuous-time Markov process with transition rates

We are concerned with the probability distribution of the population size at any time t ≥ 0 [i.e., of X(t)] . Assume that the initial population distribution is

and let

Note that if the system is initially empty, then h(z) ≡ 1.

Let

then the Kolmogorov forward differential equation for P_n(t) is (see, e.g., Ross [24] and Bailey [1])

The initial condition for these differential equations is P_n(0) = h_n, n = 0,1,….

Let

be the z-transform of P_n(t),n ≥ 0, then it follows from (2) that P(t,z) satisfies the partial differential equation (PDE)

with initial-boundary conditions

Thus, the problem boils down to solving the PDE (3) with initial condition (4). We remark that the appearance of the boundary term P(t,0) in the equation makes the problem unusual from the view point of standard PDEs. Thus, a key step is to determine P(t,0) ≡ P₀(t), which is the probability that the population size at time t is zero.

To solve the differential equation, we first study its characteristics. The characteristic equation is (see, e.g., John [9])

We note that this characteristic equation is the same as that studied in Chao and Zheng [3]. Solving this simple ordinary differential equation gives the characteristics of PDE (3):

See Chao and Zheng [3, Figs. 1–5] for the different scenarios of the characteristic curve. Some characteristic curves intersect the t-axis if and only if μ > 0. We will assume μ > 0 in this article.

Along any characteristic curve, the PDE becomes an ordinary differential equation (ODE):

For convenience, we let

We now introduce the integration factor

where τ > 0 and z₀ ∈ [0,1] are two independent variables and the function z(ζ) depends on z₀. Then, the solution of the ODE can be expressed implicitly as

where the parameter z₀ is an arbitrary number in [0,1] and z(0) = z₀.

To obtain P(t,z) at a given point (t,z), we let z(t) = z in (5) and solve for z₀ to obtain

for λ ≠ μ and

for λ = μ. Thus, the z-transform P(t,z) is given by

where z₀ is given by (8) for λ ≠ μ and by (9) for λ = μ, and z(τ) is given by (5) (where t is replaced by τ), with z₀ as just defined. The remaining key issue is to obtain a(t).

If we let z₀ > 0 be less than both μ/λ and 1, then the characteristic curve starting at z₀ must intersect the t axis at a finite time. See Figure 1. Denote this finite time by t*. By setting z(t) = 0 in the characteristic solution (5), we find that

as t → t*. Thus, the integration factor (6) satisfies

as t → t* (assuming θ > 0). Since P(t,z(t)) is expected to approach P(t*,0) as t → t*, we conclude that the second factor in (10) must approach zero for z₀ ∈ [0,min{1,μ/λ}); that is,

This is the equation that determines a(t), which is a key observation in our derivation.

Characteristics connecting the two axes.

The following is the main result of this section. The reader is referred to [5], [22], and [23] for the concept of Volterra integral equations (IEs).

Theorem 1: Suppose μ > 0 and θ > 0. Then, the following hold:

Proof:

It follows from Theorem 1 that the solution P(t, z) is obtained once we solve the singular Volterra integral equation (12). From the theory of integral equations, we know that the singular integral equation (12) can be reduced to a regular equation and this first-kind equation can be brought to the second kind, for which the iteration method can be used to construct an infinite-series solution. The details will not be given here and the reader is referred to [5], [22], or [23].

In the following sections, we consider the asymptotic behavior of the solution to this equation and we also consider several cases for which we are able to find closed-form solutions to this equation.

Remark 1: The probability P₀(t) = a(t) is determined by (11). Once P₀(t) is obtained, recursive formulas can be developed for P_n(t), n = 1,2,…, using (2), which can then be used for calculating explicit numerical solutions. However, for many other purposes, we prefer to obtain a closed-form explicit formula for P_n(t), which, among other things, shows more explicitly its dependencies on the systems parameters. Indeed, very few stochastic systems possess closed-form solutions for transient probabilities, and in Section 4, we focus on finding closed-form solutions for the process.

3. CHANCE OF EXTINCTION

The probability that the population size is zero at time t, i.e., P₀(t), is referred to as the chance of extinction (Kendall [13]). In this section we apply Theorem 1 to analyze the limiting behavior of P₀(t) as time tends to infinity.

Theorem 2: As time tends to infinity, we have the following:

Proof: We notice that (12) can be written as

For the case λ > μ, we let y → 1/(λ − μ) in (21) to find

If ν/λ ≥ θ/μ, the term (1 − s)^{−ν/λ+θ/μ−1} is not integrable on [0,1] . Thus, A(s/(λ − μ)) must vanish as s → 1; hence, A(1/(λ − μ)) = 0. In this case, we have

This proves the first half of part (a).

If λ = μ, we let y → ∞ in (21) and use h(1) = 1 to find

This implies A(∞) = 0 if ν ≥ θ and A(∞) = (θ − ν)/θ if ν < θ. This proves part (b) and completes part (a).

If λ < μ, we similarly let y → ∞ in (21) and use h(1) = 1 to find

Thus,

If ν = 0, we find A(∞) = 1. Otherwise, it is a number between zero and one since the factor (1 − λs/μ)^−ν/λ is greater than one. This proves parts (c) and (d). The proof of Theorem 2 is complete. █

Remark 2: The limit of P₀(t) for the case λ > μ and ν/λ < θ/μ is generally dependent on the initial condition h(z). For example, a(t) → 1 if h(z) ≡ 1 for the case λ > μ, ν = 0, any θ > 0, and μ > 0. However, a(t) → 1 − ((λ − μ)/λ)² if h(z) = z for the case λ > μ, ν = 0, and θ = μ. These results can be obtained from the explicit formulas of Corollaries 3 and 4 of the next section.

Remark 3: It follows from parts (b) and (d) of Theorem 2 that if λ ≤ μ and ν = 0, then lim_t→∞ P₀(t) = 1; that is, the population is “almost certain” to die out for linear birth–death process with emigration when λ ≤ μ. This result was proved in Kendall [13] for the special case ν = θ = 0.

4. CLOSED-FORM SOLUTIONS

Birth and death processes have been extensively studied in the literature and some cases have been solved explicitly a long time ago. For example, when λ = μ = 0, we obtain an M/M/1 queue and its transient solution can be found, for example, in Kleinrock [15]. The case λ = 0 and θ = 0 is equivalent to an M/M/∞ queue and its transient solution can also be found in many texts on stochastic processes (e.g., Example 4.5 of Cox and Miller [4] and Massey and Whitt [19]). The case ν = 0 and θ = 0 is the typical linear birth–death process and the case θ = 0 is a linear birth–death process with immigration, also referred to as the Kendall process (Kendall [13]). The transient solutions to all of these cases have been solved in closed forms (see, e.g., Bailey [1]).

Based on the analysis in Section 2, we present in this section the closed-form transient solutions for four cases. Our first case is covered by Ismail et al. [8], but we choose to include it here since our presentation is simple and more elementary. The transient solutions to the other cases appear to be new.

We first consider the case λ = μ. In this case, the integral equation (12) becomes

We recognize that the right-hand side of this equation is a convolution. In what follows, we use f * g to represent the convolution of f and g. Let

then,

By using the Laplace transform ℒ, we have

This analysis yields the following result.

Corollary 1: The integral equation (12) has closed-form solutions if λ = μ. In this case, a(t) = A(t) is given by (27) via Laplace transform, and the solution is

We next consider the case ν = 0, resulting in a linear birth–death process with emigration. In this case, the integral equation (12) becomes

where 0 ≤ y < 1/(λ − μ)⁺. We note that since ν = 0, we can assume that (28) holds for all y ≥ 0 (in this case, the term (1 + μy − λs)^−ν/λ is not present).

Note that the right-hand side of the IE (28) is also a convolution. Let

Then

Using the Laplace transform ℒ, we have

More specifically, we find

where Γ is the gamma function. Using αΓ(α) = Γ(α + 1), we have

If θ/μ = m is a positive integer, then

since H⁽ⁿ⁾(0) = 0 for all n = 0,1,2,…,m − 1. If θ/μ = m + α, where m is a nonnegative integer and α ∈ (0,1), then

since H⁽ⁿ⁾(0) = 0 for all n = 0,1,2,…,m. Using the convolution formula, we find

We note that (31) is valid even for α = 0; thus, (31) includes (30) as a special case.

We comment in passing that the equation for A(s) for the case θ/μ = α ∈ (0,1) is an Abel's equation; see [5], [22], [23], or [28]. In particular, if θ/μ = ½, then we have

Corollary 2: If ν = 0, μ = λ, and θ = mμ for some positive integer m, we have

where

and

with

We remark that the function G(y) is very similar to H(y). We prefer to use G(y) here. Also, the notation G^(k)(λt) denotes the kth derivative of G(y) with respect to the generic variable y at λt. One can see that λ and t appear as a product in the formula, which is consistent with the equation.

Proof: The function a(t) follows from (30) and a simple change of variable. From (15), we obtain that

Through Taylor expansion and the beta function, we find

That simplifies to

Using the expansion

we find for all n ≥ 1 that

where the sum is over all k ≥ 1 such that l + k = n, and it simplifies to yield the result stated in Corollary 2. The expressions in (33) can be verified easily through either integration by parts on (34) or inspection of their Taylor expansion. █

Corollary 3: For the case ν = 0, μ ≠ λ, and θ = mμ for some positive integer m, we have

where G(y) = h(y/(1 + y))y^m is the same as in (32), and

where s = s(t) is defined in (14), with

Proof: Once we make use of the variable s(t), the proof will be parallel to that of Corollary 2. We omit the details. █

Corollary 4: Suppose ν = 0, θ/μ = m + α > 0 for some nonnegative integer m and 0 ≤ α < 1, and μ > 0, λ > 0, we have

where s = s(t) is given in (14). The corresponding P_n(t) for all n ≥ 1 are

Proof: It is a simple consequence of Theorem 1 and (31). █

We next present an alternative method. As used in (21) in the proof of Theorem 2, we notice that (12) can always be written as

where we have defined

When ν = 0, we have

It follows from (39) that the nth derivative of F at zero is

for all nonnegative integers n. Thus, for all n, we obtain A⁽ⁿ⁾(0) as

where β is the beta function. Then, we find

Since a(τ(s)) = A(s), we have

This method can be used to deal with the general case. We can use Taylor expansion on F(y), A(sy), and (1 + (μ − λs)y)^−ν/λ at y = 0 in (37) to derive a recursive formula for A⁽ⁿ⁾(0):

where α = −ν/λ. In this way we can find an expression for a(t); however, since it is very complicated, we choose not to present its formula.