Bell Polynomial Approach for Time-Inhomogeneous Linear Birth–Death Process with Immigration

We considered the time-inhomogeneous linear birth–death processes with immigration. For these processes closed form expressions for the transition probabilities were obtained in terms of the complete Bell polynomials. The conditional mean and the conditional variance were explicitly evaluated. Several time-inhomogeneous processes were studied in detail in view of their potential applications in population growth models and in queuing systems. A time-inhomogeneous linear birth–death processes with finite state-space was also taken into account. Special attention was devoted to the cases of periodic immigration intensity functions that play an important role in the description of the evolution of dynamic systems influenced by seasonal immigration or other regular environmental cycles. Various numerical computations were performed for periodic immigration intensity functions.


Introduction
Birth-death processes are continuous-time Markov chains on the state space of non-negative integers, in which only transitions to adjacent states are allowed. These processes have been used as models in populations growth and in queuing systems and in many other fields of both theoretical and applied interest (cf., for instance, Bailey [1], Conolly [2], Feldman [3], Iosifescu and Tautu [4], Medhi [5], Ricciardi [6] and Thieme [7]).
Very often, in growth models, immigration's effects may occur, due to the circumstance that the population is not isolated. In these cases, it is necessary take into account birth-death processes with a reflecting condition in the zero state (see, for instance, Di Crescenzo et al. [8], Crawford and Suchard [9], Giorno and Nobile [10], Lenin et al. [11] and Tavaré [12]). These processes provide interesting applications in queuing models in which a reflecting boundary must be imposed to describe the number of customers in the system (cf., for instance, Crawford et al. [13], Di Crescenzo et al. [14] and Giorno et al. [15]).
In some instances, birth-death processes have also been studied under the effect of catastrophes of various types, interpretable as failures of the system, that produce transitions from the current state to the zero state from which the process can start again (see, Di Crescenzo et al. [16], Dharmaraja et al. [17], Giorno and Nobile [18], Economou and Fakinos [19] and Kapodistria et al. [20]).
Proposition 1. For t ≥ t 0 and j ∈ N 0 , the probability generating function of the NHBDI process N(t) is: Proof. The proof is given in Appendix A.
The expression of the probability generating function (5) allows us to determine the conditional mean and the conditional variance of the NHBDI process N(t). Indeed, for t ≥ t 0 and j ∈ N 0 one has:  From (5), for t ≥ t 0 we note that where G 0 (z, t) is the probability generating function of the process N(t) for j = 0 and where is the probability generating function of a linear time-inhomogeneous birth-death (NHBD) process {M(t), t ≥ t 0 }, whose intensity functions are λ n (t) = n λ(t) and µ n (t) = n µ(t) for n ∈ N. Therefore, to determine the probabilities p j,n (t|t 0 ) for j, n ∈ N 0 of the NHBDI process N(t), we proceed as follows: (1) we determine the transition probabilities p 0,n (t|t 0 ) for n ∈ N 0 and t ≥ t 0 ; (2) we calculate the probabilities p j,n (t|t 0 ) as a convolution between p 0n (t|t 0 ) and the transition probabilities of the process M(t) for j ∈ N, n ∈ N 0 and t ≥ t 0 .
We now recall some results on the NHBD process M(t), which will be useful in the following. Denoting by f j,n (t|t 0 ) = P{M(t) = n|M(t 0 ) = j}, j ∈ N, n ∈ N 0 , t ≥ t 0 the transition probabilities of M(t), by expanding (9) in powers series of z, for t ≥ t 0 one obtains (cf. Bailey [1]): where we have set:

Determination of the Transition Probabilities Starting from the Zero State
We obtain the transition probabilities for the NHBDI process N(t) when the process moves starting from the zero state.

Proposition 2.
For t ≥ t 0 and j = 0, the transition probabilities of the NHBDI process N(t) are: where A(t|t 0 ) is given in (6) and B n (d 1 , d 2 , . . . , d n ) are the complete Bell polynomials recurrently defined as follows: Proof. The proof is given in Appendix B.
Remark 1. For the NHBDI process N(t), one has: where we have set: with A(t|t 0 ) given in (6).

Determination of the Transition Probabilities Starting from j ∈ N
We determine the closed form expressions for the transition probabilities p j,n (t|t 0 ) of the NHBDI process N(t) when j ∈ N and n ∈ N 0 . Proposition 3. For t ≥ t 0 and j ∈ N, the transition probabilities of the NHBDI process N(t) are: with α(t|t 0 ), β(t|t 0 ) defined in (11).
Proof. The proof is given in Appendix C.
In Figure 1, the conditional mean and the conditional variance, given in (7), of the NHBDI process with constant birth-death rates and periodic immigration intensity function (25) are plotted as function of t for j = 15, a = 0.9, Q = 1.0 and some choices of λ, µ and ν. The dashed lines indicate the asymptotic mean E(N) and the asymptotic variance Var(N), given in (24), related to HBDI process in the case λ < µ and ν > 0.   λ n (t) = λ n + ν(t) and µ n (t) = µ n, with ν(t) given in (25), the conditional mean and the conditional variance (7) are plotted as function of t for j = 15, a = 0.9, Q = 1.0 and for some choices of λ, µ, ν. The dashed lines indicate the asymptotic mean and the asymptotic variance of the HBDI process, given in (24).
In Figure 2, making use of Remarks 1 and 2, the transition probabilities p j,n (t|0) are plotted as function of t for j = 0, 1, 2, 15 and n = 0, 1, 2 for some fixed choices of parameters. Since λ < µ and ν > 0, we note that the transition probabilities admit a periodic asymptotic behavior, strongly influenced by ν(t). Moreover, the asymptotic behavior of p j,n (t|0) oscillate around the asymptotic probabilities q n , given in (23), of the HBDI process.

Special Discrete Processes
We derive the expressions for the transition probabilities for some known and not known processes making use of the general results obtained in Section 2. We consider the processes indicated in Table 1. Table 1. Special time-inhomogeneous linear discrete processes conditioned to start from the state j.

Time-Inhomogeneous Poisson Process
We consider the time-inhomogeneous Poisson (NHP) process having state-space {j, j + 1, . . .}, conditioned to start from j ∈ N 0 at time t 0 . We assume that the birth intensity function is λ n (t) = ν(t) for n = j, j + 1, . . ., with ν(t) positive, bounded and continuous function for t ≥ t 0 . Many applications of the NHP process can be found in reliability growth models, in risk analysis, in financial problems and in queuing models. For instance, in queuing systems the NHP process is often used to describe the arrival process to a queue in which the come of customers varies according to the time of day (cf., for instance, Medhi [5], Konno [34]).
In the general NHBDI process, we set λ(t) = µ(t) = 0 and λ n (t) = ν(t) for n = j, j + 1, . . ., so that from (4) and (6) one has Λ(t|t 0 ) = M(t|t 0 ) = A(t|t 0 ) = 0. Hence, from (5), we obtain the well-known expression of the probability generating function for the NHP process: Since from (14) we have . ., the Bell complete polynomials follow from (13): For the NHP process, from (11) one has α(t|t 0 ) = β(t|t 0 ) = 0. Then, from (12) and from (17) with k = n − j and r = j, for j ∈ N 0 and t ≥ t 0 it follows: that is a time-inhomogeneous Poisson distribution. From (7), we obtain the conditional mean and the conditional variance for the NHP process N(t): Example 2. We consider the NHP process N(t), with t 0 = 0, in which the intensity function ν(t) is the periodic function given in (25). In Figure 3, for the NHP process with ν(t) as in (25), the transition probabilities, given in (26), and the conditional mean, given in (27), are plotted as function of t for some choices of parameters. We note that the conditional mean increases linearly with a periodic modulation.  (25), the probabilities p j,n (t|0) (a) and the conditional mean (b) are plotted as function of t for a = 0.9, Q = 1.0.

Time-Inhomogeneous Linear Birth Process
We consider the time-inhomogeneous linear birth (NHB) process having state-space {j, j + 1, . . .}, conditioned to start from j ∈ N at time t 0 . We assume that the birth intensity function is λ n (t) = n λ(t) for n = j, j + 1, . . ., with λ(t) positive, bounded and continuous function for t ≥ t 0 . This process is also called "generalized inhomogeneous Yule-Furry process" (cf., for instance, Konno [34]). The NHB process can be used to modeling the growth of a population of unicellular organism, such as bacteria, taking into account that whenever two new organisms born, the reproducing individual ceases to exist. The NHB process can also describe population models in which the parent organism coexists with the newly generated individual. In both the cases, the population size increases exactly by one unit as a result of a single reproduction (cf., for instance, Kendall [35], Van Den Broek and Heesterbeek [36]).

Time-Inhomogeneous Linear Death Process
We consider the time-inhomogeneous linear death (NHD) process having state-space {0, 1, . . . , j}, conditioned to start from j ∈ N at time t 0 . We assume that the death intensity function is µ n (t) = n µ(t) for n = 0, 1, . . . , j, with µ(t) positive, bounded and continuous function for t ≥ t 0 . For the NHD process, the population size decreases over time (cf., for instance, Ricciardi [6], Van Den Broek and Heesterbeek [36]).
Then, when lim t→+∞ M(t|t 0 ) = +∞, one has lim t→+∞ p j,0 (t|t 0 ) = 1, implying that the asymptotic extinction of the population is a sure event. Making use of (7), the conditional mean and the conditional variance for the NHD process N(t) are: Finally, we note that the conditional mean decreases with t and approaches to zero when lim t→+∞ M(t|t 0 ) = +∞ according to the fact that the population is doomed to extinction.

Time-Inhomogeneous Linear Birth-Death Process
The linear birth-death process (NHBD) process is obtained by combining the assumptions underlying the NHB and the NHD processes (cf., for instance, Bailey [1], Kendall [35] and Tavaré [37]). We now consider the NHBD process having space-state N 0 , conditioned to start from j ∈ N at time t 0 . We assume that the birth and death intensity functions are λ n (t) = n λ(t) and µ n (t) = n µ(t), respectively, with λ(t) and µ(t) positive, bounded and continuous functions for t ≥ t 0 .
In the general NHBDI process, we set ν(t) = 0, λ n (t) = n λ(t) and µ n (t) = n µ(t) for n ∈ N. Hence, from (5) we derive the probability generating function for the NHBD process: where Λ(t|t 0 ), M(t|t 0 ) and A(t|t 0 ) are defined in (4) and (6), respectively. Recalling (13) and (14), we have: d r = 0 for r ∈ N, B 0 = 1 and B n (d 1 , d 2 , . . . , d n ) = 0 for n ∈ N. Furthermore, by choosing k = 0 in (17), for j ∈ N and t ≥ t 0 one obtains: where α(t|t 0 ) and β(t|t 0 ) are defined in (11). For the NHBD process, the use of (7) for j ∈ N and t ≥ t 0 leads to: Moreover, by setting n = 0 in (30), for j ∈ N one has: so that for the NHBD process the probability of ultimate extinction tends to unity as t → +∞ if and only if

Time-Homogeneous Linear Birth-Death Process
A special case is the time-homogeneous linear birth-death (HBD) process, characterized by birth and death rates λ n (t) = n λ and µ n (t) = n µ for n ∈ N 0 , with λ and µ positive real numbers (cf., for instance, Bailey [1], Ricciardi [6] and Crawford and Suchard [9]). For the HBD process, by assuming that ν → 0 in (20) and in (21), we obtain the transition probabilities. For λ = µ and j ∈ N one has: with (t|t 0 ) defined in (19), whereas, for λ = µ and j ∈ N one obtains: Moreover, from (33) and (34), for j ∈ N and t ≥ t 0 it follows: Hence, for the HBD process the probability of ultimate extinction tends to unity as t → +∞ if and only if λ ≤ µ. Furthermore, for ν → 0 in (22), one derives the conditional mean and conditional variance for the HBD process: Then, for the HBD process the mean population size exponentially increases for λ > µ, exponentially decreases if λ < µ and remains constant if λ = µ.
In Figure 4, the conditional mean and the conditional variance (35) of the HBD process are plotted as function of t for some choices of λ and µ.

Time-Inhomogeneous Linear Death Process with Immigration
We consider the time-inhomogeneous linear death with immigration (NHDI) process having state-space N 0 , conditioned to start from j ∈ N 0 at time t 0 . We assume that the birth and death intensity functions are λ n (t) = ν(t) and µ n (t) = n µ(t) for n ∈ N 0 , with ν(t) and µ(t) positive, bounded and continuous functions for t ≥ t 0 (cf., for instance, Ohkubo [38]). The NHDI process also describes the multi-server queuing systems M(t)/M(t)/∞ in which the arrivals are governed by a time-inhomogeneous Poisson process with intensity function ν(t), there are infinitely many servers in parallel and the service times have an non-homogeneous exponential density with intensity function µ(t). In the M(t)/M(t)/∞ queue, there are always sufficient servers such that every arriving job is served immediately (cf., for instance, Di Crescenzo et al. [8], Giorno et al. [23]).

Time-Homogeneous Linear Death Process with Immigration
A special case is the time-homogeneous linear death-immigration (HDI) process, characterized by birth and death rates λ n (t) = ν and µ n (t) = n µ for n ∈ N 0 , with ν and µ positive real numbers. The HDI process can be used to describe the M/M/∞ queuing system with infinitely many servers in parallel, exponential interarrival and service times with mean 1/ν and 1/µ, respectively (cf., for instance, Medhi [5], Di Crescenzo et al. [8] and Giorno et al. [23]). For the HDI process, by setting λ = 0 in (20), for t ≥ t 0 we obtain the transition probabilities: j ∈ N, n ∈ N 0 and for λ → 0 in (22), for j ∈ N 0 and t ≥ t 0 we obtain the conditional mean and the conditional variance of the HDI process: Furthermore, from (38) it follows that the HDI process always admits the Poisson steady-state distribution: so that the asymptotic mean and variance are E(N) = Var(N) = ν/µ.

Example 3.
We consider the NHDI process N(t), with t 0 = 0, characterized by λ n (t) = ν(t) and µ n (t) = µ n, with µ > 0, subject to periodic immigration phenomena that occur with intensity function (25). This process can describe a queuing system M(t)/M/∞, in which the customers arrive according to a time-inhomogeneous Poisson process with the periodic intensity function (25). In Figure 5, the transition probabilities p j,n (t|0), given in (36), for the NHDI process with constant death rate and periodic immigration intensity function (25) are plotted as function of t for j = 2, 15 and n = 0, 1, 2 for some fixed choices of parameters. We note that, for large times, the transition probabilities p j,n (t|0) oscillate around the probabilities q n , given in (39), related to the HDI process.  For the NHDI process, having λ n (t) = ν(t) and µ n (t) = µ n, with ν(t) given in (25), the probabilities p j,n (t|0) are plotted as function of t for µ = 0.8, ν = 2.0, a = 0.9, Q = 1.0 with n = 0, 1, 2. The dashed lines indicate the probabilities q 0 , q 1 , q 2 , given in (39).
In Figure 6, the conditional mean and the conditional variance, given in (37), of the NHDI process of Figure 5 are plotted as function of t for j = 15, a = 0.9, Q = 1.0 and some choices of µ and ν. The dashed lines indicate E(N) = Var(N) = ν/µ, related to HDI process. In the M(t)/M/∞ queue this means that, for large times, the number of the customers exhibits a periodic behavior, which strongly depends on the periodicity of the arrivals in the system.

Time Inhomogeneous Linear Birth Process with Immigration
We consider the time-inhomogeneous linear birth with immigration (NHBI) process having state-space {j, j + 1, . . .}, conditioned to start to j ∈ N 0 at time t 0 . We assume that the birth intensity function is λ n (t) = λ(t) n + ν(t) for n = j, j + 1, . . ., with λ(t) and ν(t) positive, bounded and continuous functions for t ≥ t 0 . In the NHBI process a new individual can be regarded as arising from two sources, one involving the time-inhomogeneous Poisson process with intensity function ν(t), the other being a linear time-inhomogeneous birth process with intensity function λ(t) for individual.

Time-Homogeneous Linear Birth Immigration Process
A special case is the time-homogeneous linear birth immigration (HBI) process, characterized by birth and death rates λ n (t) = n λ + ν for n = j, j + 1, . . . and j ∈ N 0 , with λ and ν positive real numbers (cf., for instance, Tavaré [12]). For the HBI process, by setting µ = 0 in (20), for t ≥ t 0 we obtain the transition probabilities: Furthermore, for µ → 0 in (22), for j ∈ N 0 and t ≥ t 0 we have the conditional mean and the conditional variance of the HBI process:
In Figure 8, the conditional mean and the conditional variance, given in (44), of the NHBI process of Figure 7 are plotted as function of t for j = 0, a = 0.9, Q = 1.0 and some choices of λ and ν. We note that the shapes of the conditional mean and of the conditional variance are little influenced by the periodicity of ν(t).
For the GPy process, recalling that ν(t) = ν λ(t), from (40) one has: Moreover, from (41) we have d r = ν (r − 1)! 1 − e −Λ(t|t 0 ) r for r ∈ N, so that from (13) we have B 0 = 1 and The solution of the recurrence equation (49) is: obtained by using the following identities: with Γ(x) denoting the Euler's gamma function. By virtue of (50), from (42) and (43) for t ≥ t 0 one obtains the transition probabilities of the GPy process: where the last identity follows by noting that Finally, from (44), for j ∈ N 0 and t ≥ t 0 the conditional mean and the conditional variance of the GPy process can be derived: Example 5. We consider the GPy process N(t), with t 0 = 0, characterized by λ n (t) = λ(t) (n + ν), with ν > 0 and with λ > 0 and 0 ≤ a < 1.
In Figure 10, the conditional mean and the conditional variance, given in (53), of the GPy process of Figure 9 are plotted as function of t for j = 0, λ = 1.0, a = 0.9, Q = 1.0 and some choices of ν.   Figure 9, the conditional mean and the conditional variance (53) are plotted as function of t for j = 0, λ = 1.0, a = 0.9, Q = 1.0 and for some choices of ν.

Generalized Polya-Death Process
We consider the generalized Polya-death (GPyD) process {N(t), t ≥ t 0 }, obtained by setting ν(t) = ν λ(t) in the NHBDI process. Then, in the GPyD process we assume that the birth and death intensity functions are λ n (t) = λ(t) (n + ν) for n ∈ N 0 and µ n (t) = n µ(t) for n ∈ N, with ν positive real number and λ(t), µ(t) bounded and continuous functions for t ≥ t 0 (cf., for instance, Ohkubo [39]). A special case of GPyD process has been considered in Giorno et al. [23] and Di Crescenzo and Nobile [40] to describe a time-inhomogeneous adaptive queuing system, by assuming that λ n (t) = (λ n + α) k(t) for n ∈ N 0 and µ n (t) = n µ k(t) for n ∈ N. Moreover, a time-homogeneous GPyD process, characterized by λ n (t) = λ (n + α) and µ n (t) = n µ, has been used to describe an adaptive queuing system, known as "Model D" with panic-buying and compensatory reaction of service (cf., for instance, Conolly [2], Lenin et al. [11] and Giorno and Nobile [18]).

Proposition 4.
For t ≥ t 0 , the transition probabilities of the GPyD process are: Proof. Thew proof is given in Appendix D.
Finally, from (7) for j ∈ N 0 we obtain the conditional mean and the conditional variance for the GPyD process:
In Figure 12, the conditional mean, given in (61), of the GPyD process of Figure 11 is plotted as function of t for j = 0, λ = 1.0, a = 0.9, Q = 1.0 and some choices of ν. The dashed lines in Figure 12b give the asymptotic averages E(N) = λν/(µ − λ) of the time-homogeneous Polya-death process. We note that if λ < µ, for large times the conditional mean oscillates around the asymptotic mean E(N) of the time-homogeneous Polya-death process.
In Figure 13, the transition probabilities p j,n (t|0), given in (68), of the time-inhomogeneous Prendiville process, with λ(t) given in (54), are plotted as function of t with j = 15, K = 30, a = 0.9, Q = 1.0 and for some choices of λ and µ. The dashed line indicate the steady-state probabilities q 0 , q 1 , q 2 on the left and the steady-state probabilities q 28 , q 29 , q 30 on the right, given in (71), of the time-homogeneous Prendiville process having λ n (t) = (K − n) λ and µ n (t) = n µ.  Figure 13. For the Prendiville process, having λ n (t) = λ(t) (n + ν) and µ n (t) = µ n, with λ(t) given in (54), the probabilities p j,n (t|0) are plotted as function of t with j = 15, K = 30, a = 0.9, Q = 1.0 and for some choices of λ and µ.
In Figure 14, the conditional mean and the conditional variance are plotted as function of t for j = 15, K = 30, a = 0.9, Q = 1.0 and some choices of λ, µ. The dashed line indicate the asymptotic mean and the asymptotic variance, given in (72), of the time-homogeneous Prendiville process.

Conclusions
In this paper, we have presented a detailed analysis of the time-inhomogeneous linear birth-death processes with immigration. The transition probabilities, the conditional averages and the conditional variances are determined in closed form. Specifically, in Section 2 the transition probabilities are obtained in terms of the complete Bell polynomials. Special time-inhomogeneous processes of interest in population growth models and in queuing systems are carefully analyzed in Section 3. A time-inhomogeneous linear birth-death processes with finite state-space is also taken into account in Section 4. Special attention is devoted to the cases of periodic immigration intensity functions, that play an important role in the description of the evolution of dynamic systems in various applied fields. For instance, in population dynamics they express the existence of fluctuation in the growth due to seasonal immigration or other regular environmental cycles. Various numerical computations with MATHEMATICA are performed in the cases of periodic immigration intensity functions.
We conclude by mentioning that future research will concern the analysis of the first-passage time problem for time-inhomogeneous birth-death processes with immigration through specific boundaries interpretable as the extinction level or the saturation level (carrying capacity) in population dynamics. Moreover, our investigation will be extended to the continuous approximations of the time-inhomogeneous birth-death processes with immigration. It is expected to obtain new theoretical results on stochastic time-inhomogeneous diffusion processes, restricted to the interval (0, +∞), where the state zero is a reflecting or absorbing boundary. These studies will be used to model biological, physical and chemical systems. The knowledge of the transition distributions for discrete processes and their continuous approximations will allow the determination of the first-passage densities both through analytical methods and through numerical and simulation techniques.