On the Generalized Birth-Death Process and its Linear Versions

P. Vishwakarma Pradeep Vishwakarma, Department of Mathematics, Indian Institute of Technology Bhilai, Durg 491001, INDIA. [email protected] and K.k. Kataria Kuldeep Kumar Kataria, Department of Mathematics, Indian Institute of Technology Bhilai, Durg 491001, INDIA. [email protected]

Abstract.

In this paper, we consider a generalized birth-death process (GBDP) and examined its linear versions. Using its transition probabilities, we obtain the system of differential equations that governs its state probabilities. The distribution function of its waiting-time in state $s$ given that it starts in state $s$ is obtained. For a linear version of it, namely, the generalized linear birth-death process (GLBDP), we obtain the probability generating function, mean, variance and the probability of ultimate extinction of population. Also, we obtain the maximum likelihood estimate of one of its parameter. The differential equations that govern the joint cumulant generating functions of the population size with cumulative births and cumulative deaths are derived. In the case of constant birth and death rates in GBDP, the explicit forms of the state probabilities, joint probability mass functions of population size with cumulative births and cumulative deaths, and their marginal probability mass functions are obtained. It is shown that the Laplace transform of a stochastic integral of GBDP satisfies its Kolmogorov backward equation with certain scaled parameters. Also, the first two moments of the stochastic path integral of GLBDP are obtained. Later, we consider the immigration effect in GLBDP for two different cases. An application of a linear version of GBDP and its stochastic path integral to vehicles parking management system is discussed.

Key words and phrases:

Extinction probability; birth-death process; generalized birth-death process; generalized linear birth-death process; generalized linear birth-death process with immigration.

2010 Mathematics Subject Classification:

Primary : 60J27; Secondary: 60J20

1. Introduction

The birth-death process (BDP) is a continuous-time and discrete state-space Markov process which is used to model the growth of population with time. If $\lambda_{n}>0$ for all $n\geq 0$ are the birth rates and $\mu_{n}>0$ for all $n\geq 1$ are the death rates then, in an infinitesimal time interval of length $h$ , the probability of a birth equals $\lambda_{n}h+o(h)$ , probability of a death equals $\mu_{n}h+o(h)$ and probability of any other event is $o(h)$ . Here, $n$ denotes the population size at a given time $t\geq 0$ . In BDP, the transition rates can depend on the state of the process. If the transition rates in BDP are linear, that is, $\lambda_{n}=n\lambda$ and $\mu_{n}=n\mu$ then it is called the linear birth-death process (LBDP). The BDP has various applications in the fields of engineering, queuing theory and biological sciences etc. For example, it can be used to model the number of customers in queue at a service center or the growth of bacteria over time. However, it has certain limitation, for example, it is not a suitable process to model the situations involving multiple births or multiple deaths in an infinitesimal time interval.

Doubleday (1973) introduced and studied a linear birth-death process with multiple births and single death. Here, we consider a generalized version of the BDP where in an infinitesimal time interval there is a possibility of multiple but finitely many births or deaths with the assumption that the chance of simultaneous birth and death is negligible. We call this process as the generalized birth-death process (GBDP) and denote it by $\{\mathcal{N}(t)\}_{t\geq 0}$ . It is shown that the state probabilities $q(n,t)=\mathrm{Pr}\{\mathcal{N}(t)=n\}$ , $n\geq 0$ of GBDP solve the following system of differential equations:

\frac{\mathrm{d}}{\mathrm{d}t}q(n,t)=\begin{cases}-\sum_{i=1}^{k_{1}}\lambda_{(0)_{i}}q(0,t)+\sum_{j=1}^{k_{2}}\mu_{(j)_{j}}q(j,t),\ n=0,\vspace{0.1cm}\\ -\big{(}\sum_{i=1}^{k_{1}}\lambda_{(n)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(n)_{j}}\big{)}q(n,t)\\ \hskip 51.21504pt+\sum_{i=1}^{\mathrm{min}\{n,k_{1}\}}\lambda_{(n-i)_{i}}q(n-i,t)+\sum_{j=1}^{k_{2}}\mu_{(n+j)_{j}}q(n+j,t),\ n\geq 1\end{cases}

(1.1)

with initial conditions $q(1,0)=1$ and $q(n,0)=0$ for all $n\neq 1$ . Here, for each $n\geq 0$ , $\lambda_{(n)_{i}}$ is the rate of birth of size $i\in\{1,2,\dots,k_{1}\}$ and for each $n\geq 1$ , $\mu_{(n)_{j}}$ is the rate of death of size $j\in\{1,2,\dots,k_{2}\}$ .

If birth and death rates are linear in GBDP, that is, $\lambda_{(n)_{i}}=n\lambda_{i}$ and $\mu_{(n)_{j}}=n\mu_{j}$ for all $n\geq 1$ , then we call this process as the generalized linear birth-death process (GLBDP). It is denoted by $\{N(t)\}_{t\geq 0}$ . The state probabilities $p(n,t)=\mathrm{Pr}\{N(t)=n\}$ of GLBDP are the solution of the following system of differential equations:

\frac{\mathrm{d}}{\mathrm{d}t}p(n,t)=\begin{cases}\sum_{j=1}^{k_{2}}j\mu_{j}p(j,t),\ n=0,\vspace{0.1cm}\\ -\big{(}\sum_{i=1}^{k_{1}}n\lambda_{i}+\sum_{j=1}^{k_{2}}n\mu_{j}\big{)}p(n,t)\\ \hskip 54.06006pt+\sum_{i=1}^{k_{1}}(n-i)\lambda_{i}p(n-i,t)+\sum_{j=1}^{k_{2}}(n+j)\mu_{j}p(n+j,t),\ n\geq 1,\end{cases}

where $p(n,t)=0$ for all $n<0$ .

It is shown that the cumulative distribution function $W_{s}(t)=\mathrm{Pr}\{T_{s}\leq t\}$ is given by

W_{s}(t)=\exp\Bigg{(}-\sum_{i=1}^{k_{1}}\lambda_{(n)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(n)_{j}}\Bigg{)}t,\ t\geq 0,

where $T_{s}$ is the waiting time of GBDP in state $s$ given that it starts in state $s$ . On using the distributions of these waiting times, we obtain the maximum likelihood estimate of the parameter $\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}$ in GLBDP.

In Section 3, for GLBDP, we obtain some results for the cumulative births $B(t)$ , that is, the total number of births by time $t$ and for the cumulative deaths $D(t)$ , that is, the total number of deaths by time $t$ . We study a particular case of GBDP, where the birth and death rates are constant. We denote this process by $\{N^{*}(t)\}_{t\geq 0}$ . Let $B^{*}(t)$ and $D^{*}(t)$ be the total number of births and deaths in $\{N^{*}(t)\}_{t\geq 0}$ up to time $t$ , respectively. We obtain the explicit forms of probability mass functions (pmfs) of $N^{*}(t)$ , $B^{*}(t)$ and $D^{*}(t)$ and their joint distributions.

In Section 4, we consider a stochastic integral of GBDP and show that its Laplace transform satisfies the backward equation for GBDP with some scaled parameters. We study the joint distribution of GLBDP and its stochastic path integral. A similar study is done for GBDP with constant birth and death rates.

In Section 5, we study the effect of immigration in GLBDP for two different cases. Let $\nu>0$ is the rate of immigration from outside environment. First, we consider the immigration effect only if the population vanishes at any time $t\geq 0$ . In this case, the birth rates are $\lambda_{(0)_{i}}=\nu$ , $\lambda_{(n)_{i}}=n\lambda_{i},\ n\geq 1$ for all $i\in\{1,2,\dots,k_{1}\}$ and the death rates are $\mu_{(n)_{j}}=n\mu_{j}$ , $n\geq 1$ for all $j\in\{1,2,\dots,k_{2}\}$ . Then, we consider the effect of immigration at every state of the GLBDP. In this case, the birth and death rates are $\lambda_{(n)_{i}}=\nu+n\lambda_{i}$ and $\mu_{(n)_{j}}=n\mu_{j}$ , $n\geq 0$ , respectively.

In the last section, we discuss an application of a linear case of GBDP to vehicle parking management system.

2. Generalized birth-death process

First, we consider a generalized birth-death process (GBDP) where in an infinitesimal time interval of length $h$ there is a possibility of either finitely many births $i\in\{1,2,\ldots,k_{1}\}$ with positive rates $\lambda_{(n)_{i}}$ or finitely many deaths $j\in\{1,2,\ldots,k_{2}\}$ with positive rates $\mu_{(n)_{j}}$ . It is important to note that the birth and death rates may depend on the number of individuals $n$ present in the population at time $t\geq 0$ . It is assumed that the chances of simultaneous occurrence of births and deaths in such small intervals are negligible.

We denote the GBDP by $\{\mathcal{N}(t)\}_{t\geq 0}$ . With the above assumptions, its transition probabilities are given by

\mathrm{Pr}\{\mathcal{N}(t+h)=n+j|\mathcal{N}(t)=n\}=\begin{cases}1-\sum_{i=1}^{k_{1}}\lambda_{(n)_{i}}h-\sum_{i=1}^{k_{2}}\mu_{(n)_{i}}h+o(h),\ j=0,\vspace{0.1cm}\\ \lambda_{(n)_{j}}h+o(h),\ j=1,2,\ldots,k_{1},\vspace{0.1cm}\\ \mu_{(n)_{-j}}h+o(h),\ j=-1,-2,\ldots,-k_{2},\vspace{0.1cm}\\ o(h),\ \mathrm{otherwise},\end{cases}

(2.1)

where $o(h)/h\to 0$ as $h\to 0$ .

Let $q(n,t)=\mathrm{Pr}\{\mathcal{N}(t)=n\}$ , $n\geq 0$ be the state probabilities of GBDP. Then, we have

	$\displaystyle q(n,t+h)$	$\displaystyle=\mathrm{Pr}\{\mathcal{N}(t+h)=n\|\mathcal{N}(t)=n\}q(n,t)+\sum_{i=1}^{\mathrm{min}\{n,k_{1}\}}\mathrm{Pr}\{\mathcal{N}(t+h)=n\|\mathcal{N}(t)=n-i\}q(n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}\mathrm{Pr}\{\mathcal{N}(t+h)=n\|\mathcal{N}(t)=n+j\}q(n+j,t).$

Using (2.1), we get

	$\displaystyle q(n,t+h)$	$\displaystyle=\bigg{(}1-\sum_{i=1}^{k_{1}}\lambda_{(n)_{i}}h-\sum_{j=1}^{k_{2}}\mu_{(n)_{j}}h\bigg{)}q(n,t)+\sum_{i=1}^{\mathrm{min}\{n,k_{1}\}}\lambda_{(n-i)_{i}}hq(n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}\mu_{(n+j)_{j}}hq(n+j,t)+o(h).$

Equivalently,

	$\displaystyle\frac{q(n,t+h)-q(n,t)}{h}$	$\displaystyle=-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{(n)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(n)_{j}}\bigg{)}q(n,t)+\sum_{i=1}^{\mathrm{min}\{n,k_{1}\}}\lambda_{(n-i)_{i}}q(n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}\mu_{(n+j)_{j}}q(n+j,t)+\frac{o(h)}{h}.$

On letting $h\to 0$ , we get

\frac{\mathrm{d}}{\mathrm{d}t}q(n,t)=-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{(n)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(n)_{j}}\bigg{)}q(n,t)+\sum_{i=1}^{\mathrm{min}\{n,k_{1}\}}\lambda_{(n-i)_{i}}q(n-i,t)+\sum_{j=1}^{k_{2}}\mu_{(n+j)_{j}}q(n+j,t),\ n\geq 1.

Similarly, for $n=0$ , we have

q(0,t+h)=\bigg{(}1-\sum_{i=1}^{k_{1}}\lambda_{(0)_{i}}h\bigg{)}q(0,t)+\sum_{j=1}^{k_{2}}\mu_{(j)_{j}}q(j,t)h+o(h).

Thus,

\frac{\mathrm{d}}{\mathrm{d}t}q(0,t)=-\sum_{i=1}^{k_{1}}\lambda_{(0)_{i}}q(0,t)+\sum_{j=1}^{k_{2}}\mu_{(j)_{j}}q(j,t).

Let us assume that there is one progenitor at time $t=0$ . Thus, the state probabilities of GBDP solves the following system of differential equations:

\frac{\mathrm{d}}{\mathrm{d}t}q(n,t)=-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{(n)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(n)_{j}}\bigg{)}q(n,t)+\sum_{i=1}^{k_{1}}\lambda_{(n-i)_{i}}q(n-i,t)+\sum_{j=1}^{k_{2}}\mu_{(n+j)_{j}}q(n+j,t),\ n\geq 0

(2.2)

with initial conditions $q(1,0)=1$ and $q(n,0)=0$ for all $n\neq 1$ . Here, $q(n-i,t)=0$ for all $i>n$ and $\mu_{(0)_{j}}=0$ for all $j=1,2,\ldots,k_{2}$ .

Remark 2.1.

For $k_{1}=k_{2}=1$ , GBDP reduces to BDP. If $\mu_{(n)_{j}}=0$ for all $n$ and for all $j\in\{1,2,\ldots,k_{2}\}$ , then the GBDP reduces to the generalized pure birth process.

Let $T_{s}$ be the waiting time of GBDP in the state $s$ , that is, $T_{s}$ is the total time before the process leave the state $s$ given that it start from $s$ .

Proposition 2.1.

Let $W_{s}(t)=\mathrm{Pr}\{T_{s}\geq t\}.$ Then, $W_{s}(t)$ solves the following linear differential equation:

\frac{\mathrm{d}}{\mathrm{d}t}W_{s}(t)=-\left(\sum_{i=1}^{k_{1}}\lambda_{(s)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(s)_{j}}\right)W_{s}(t)

(2.3)

with initial condition $W_{s}(0)=1.$

Proof.

Consider an infinitesimal time interval of length $h$ such that $o(h)/h\to 0$ as $h\to 0$ . On using (2.1), we have

\displaystyle W_{s}(t+h)

\displaystyle=W_{s}(t)\mathrm{Pr}\{\mathcal{N}(t+h)=s|\mathcal{N}(t)=s\}=W_{s}(t)\bigg{(}1-\sum_{i=1}^{k_{1}}\lambda_{(s)_{i}}h-\sum_{j=1}^{k_{2}}\mu_{(s)_{j}}h\bigg{)}+o(h).

So,

\frac{W_{s}(t+h)-W_{s}(t)}{h}=-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{(s)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(s)_{j}}\bigg{)}W_{s}(t)+\frac{o(h)}{h}.

On letting $h\to 0$ , we get the required result. ∎

Remark 2.2.

On solving (2.3), we get

W_{s}(t)=e^{-\left(\sum_{i=1}^{k_{1}}\lambda_{(s)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(s)_{j}}\right)t},\ t\geq 0.

Thus, the waiting time $T_{s}$ of GBDP is exponentially distributed with parameter $\sum_{i=1}^{k_{1}}\lambda_{(s)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(s)_{j}}$ . For $k_{1}=k_{2}=1$ , it follows exponential distribution with parameter $\lambda_{s}+\mu_{s}$ (see Karlin and Taylor (1975), p. 133).

2.1. A linear case of GBDP

When the birth and death rates are linear, that is, $\lambda_{(n)_{i}}=n\lambda_{i}$ and $\mu_{(n)_{j}}=n\mu_{j}$ for all $i$ , $j$ and $n$ , we call the GBDP as the generalized linear birth-death process (GLBDP) and denote it by $\{N(t)\}_{t\geq 0}$ . From (2.2), it follows that the state probabilities $p(n,t)=\mathrm{Pr}\{N(t)=n\}$ of GLBDP solve the following system of differential equations:

\frac{\mathrm{d}}{\mathrm{d}t}p(n,t)=-\left(\sum_{i=1}^{k_{1}}n\lambda_{i}+\sum_{j=1}^{k_{2}}n\mu_{j}\right)p(n,t)+\sum_{i=1}^{k_{1}}(n-i)\lambda_{i}p(n-i,t)+\sum_{j=1}^{k_{2}}(n+j)\mu_{j}p(n+j,t),\ n\geq 0,

(2.4)

with initial conditions

p(n,0)=\begin{cases}1,\ n=1,\vspace{0.1cm}\\ 0,\ n\neq 1,\end{cases}

(2.5)

where $p(n-i,t)=0$ for all $i>n$ .

Remark 2.3.

If $\mu_{j}=0$ for all $j\in\{1,2,\ldots,k_{2}\}$ , then the GLBDP reduces to generalized linear pure birth process.

Remark 2.4.

In the case of GLBDP, the waiting time in state $s$ has the exponential distribution with parameter $s\big{(}\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\big{)}$ .

Next, we obtain the system of differential equations that governs the probability generating function (pgf) of GLBDP.

Proposition 2.2.

Let $H(u,t)=\mathbb{E}(u^{N(t)}),\ |u|\leq 1$ be the pgf of GLBDP. It solves the following partial differential equation:

\frac{\partial}{\partial t}H(u,t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}u(u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}u(u^{-j}-1)\right)\frac{\partial}{\partial u}H(u,t),\ t\geq 0,

(2.6)

with initial condition $H(u,0)=u$ .

Proof.

By using (2.5), we have $H(u,0)=\sum_{n=0}^{\infty}u^{n}p(n,0)=u.$ On multiplying $u^{n}$ on both sides of (2.4) and by adjusting the terms, we get

	$\displaystyle\frac{\partial}{\partial t}u^{n}p(n,t)$	$\displaystyle=-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\bigg{)}u\frac{\partial}{\partial u}u^{n}p(n,t)+\sum_{i=1}^{k_{1}}\lambda_{i}u^{i+1}\frac{\partial}{\partial u}u^{n-i}p(n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}\mu_{j}u^{-j+1}\frac{\partial}{\partial u}u^{n+j}p(n+j,t).$		(2.7)

Summing over $n=-k_{2},-k_{2}+1,\ldots$ on both sides of (2.7) gives the required result. Here, the change of sum and derivative is justified because

	$\displaystyle\left\|\sum_{n=0}^{\infty}\frac{\partial}{\partial u}u^{n}p(n,t)\right\|$	$\displaystyle=\left\|\sum_{n=0}^{\infty}nu^{n-1}p(n,t)\right\|$
		$\displaystyle\leq\sum_{n=0}^{\infty}np(n,t)=\mathbb{E}(N(t))<\infty.$

This completes the proof. ∎

Remark 2.5.

For $k_{1}=k_{2}=1$ , the GLBDP reduces to LBDP. Its pgf solves (see Bailey (1964))

\frac{\partial}{\partial t}{H}(u,t)=(\lambda_{1}u-\mu_{1})(u-1)\frac{\partial}{\partial u}{H}(u,t)

with ${H}(u,0)=u$ .

Remark 2.6.

The Laplace transform $\tilde{H}(u,z)=\int_{0}^{\infty}e^{-zt}H(u,t)\,\mathrm{d}t$ of the pgf of GLBDP is the solution of following differential equation:

-u+z\tilde{H}(u,z)=\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}u(u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}u(u^{-j}-1)\bigg{)}\frac{\partial}{\partial u}\tilde{H}(u,z).

Next, we discuss the solution of (2.6). First, we rewrite (2.6) as follows:

\frac{\partial}{\partial t}H(u,t)=\phi(u)\frac{\partial}{\partial u}H(u,t),

(2.8)

where $\phi(u)=\sum_{i=1}^{k_{1}}\lambda_{i}u(u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}u(u^{-j}-1).$ We apply the method of characteristics to solve (2.8), and thus we get the following subsidiary equations:

\mathrm{d}t=-\frac{u^{k_{2}-1}}{\psi(u)}\mathrm{d}u

(2.9)

and

\frac{\partial}{\partial t}H(u,t)=0,

(2.10)

where $\psi(u)=u^{k_{2}-1}\phi(u)$ . Moreover, we assume that the roots $r_{i},\ i=1,2,\ldots,k_{1}+k_{2}$ of $\psi(u)$ are distinct which can be ensured by a small perturbation in $\lambda_{i}$ and $\mu_{j}$ . So, (2.9) can be rewritten as

\mathrm{d}t=-\left(\sum_{i=1}^{k_{1}+k_{2}}\frac{c_{i}}{(u-r_{i})}\right)\mathrm{d}u,

(2.11)

where $c_{i}=\lambda_{k_{1}}^{-1}{r_{i}^{k_{2}-1}}/{\prod_{j\neq i}(r_{i}-r_{j})}.$

On solving (2.10) and (2.11), we get $H(u,t)=\beta_{1}$ and $e^{-t}g(u)=\beta_{2}$ , respectively. Here, $\beta_{1}$ and $\beta_{2}$ are real constants and

g(u)=\prod_{i=1}^{k_{1}+k_{2}}(u-r_{i})^{-c_{i}}.

(2.12)

Thus, the arbitrary solution of (2.6) is obtained in the following form:

H(u,t)=f\left(e^{-t}g(u)\right),

(2.13)

where $f$ is some real-valued function which can be determined by using the given initial condition, that is, $H(u,0)=u=f\left(g(u)\right)$ . Note that there always exist a $\delta>0$ such that $|r_{i}|>\delta$ and $\psi(u)$ has same sign for all $u\in[0,\delta]$ . So, $g(u)$ is a monotone function in $[0,\delta]$ and its inverse exist. Hence, (2.13) holds with $f(u)=g^{-1}(u)$ for all $u\in[0,\delta]$ and by using the analytic continuation it can be extended to all $u$ such that $|u|\leq 1$ .

Remark 2.7.

Let $\left(a_{k}:k\in\{0,1,2,\ldots\}\right)$ be the initial distribution of GLBDP. Then, its pgf has the following form:

H(u,t)=\sum_{k=0}^{\infty}a_{k}\left(f\left(e^{-t}\prod_{i=1}^{k_{1}+k_{2}}(u-r_{i})^{-c_{i}}\right)\right)^{k}.

In particular, if $N(0)=n_{0}>1$ then

H(u,t)=\left(f\left(e^{-t}\prod_{i=1}^{k_{1}+k_{2}}(u-r_{i})^{-c_{i}}\right)\right)^{n_{0}}.

Its proof follows along the similar lines to that of (2.13).

Next, we obtain the probability of ultimate extinction of GLBDP.

Let $p(0)=\lim_{t\to\infty}p(0,t)$ be the probability of ultimate extinction in GLBDP. We use the inverse function representation of $H(u,t)$ to approximate the extinction probability $p(0,t)$ . Let ${g}^{\prime}(u)=\frac{\mathrm{d}g(u)}{\mathrm{d}u}$ , where $g$ is given in (2.12). If $\epsilon$ is the least positive root of $\psi(u)$ then $g(\cdot)$ is a monotone function in $[0,\epsilon]$ and $g(\epsilon)=0$ . So, $f(u)=g^{-1}(u)$ for all $u\in[0,\epsilon]$ and $\ f(0)=\epsilon.$ As

	$\displaystyle f(u)$	$\displaystyle=\epsilon+\int_{0}^{u}f^{\prime}(x)\,\mathrm{d}x$
		$\displaystyle=\epsilon+\int_{0}^{u}\frac{\mathrm{d}x}{{g}^{\prime}(f(x))}\ \ \ \text{for}\ u\in[0,\delta],$

where $\delta$ is some positive real number. So, we have

\displaystyle H(u,t)=\epsilon+\int_{0}^{e^{-t}g(u)}\frac{\mathrm{d}x}{{g}^{\prime}(f(x))}.

On taking $u=0$ , we get

p(0,t)=H(0,t)=\epsilon+\int_{0}^{e^{-t}g(0)}\frac{\mathrm{d}x}{{g}^{\prime}(f(x))}.

Thus, if the process run for infinitely long time then the probability of ultimate extinction converges exponentially to $\epsilon$ , that is,

p(0)=\lim_{t\to\infty}p(0,t)=\epsilon.

(2.14)

Remark 2.8.

In the case of LBDP, the above results reduce to the results given in Bailey (1964).

For $k_{1}=k_{2}=1$ , $\lambda_{1}=\lambda$ and $\mu_{1}=\mu$ such that $\lambda\neq\mu$ , we get the two distinct roots of $\psi(u)$ , that is, $r_{1}=\mu/\lambda$ and $r_{2}=1$ . So, $c_{1}=1/(\mu-\lambda)$ and $c_{2}=1/(\lambda-\mu)$ . From (2.12), we have

g(u)=\left(\frac{u-\mu/\lambda}{u-1}\right)^{1/(\lambda-\mu)}.

Thus,

f(u)=g^{-1}(u)=\frac{\mu-\lambda u^{\lambda-\mu}}{\lambda(1-u^{\lambda-\mu})}.

So, from (2.13), the pgf of LBDP is

H(u,t)=\frac{\mu(u-1)-e^{-t(\lambda-\mu)}(\lambda u-\mu)}{\lambda(u-1)-e^{-t(\lambda-\mu)}(\lambda u-\mu)}.

Its extinction probability is given by (see Bailey (1964), p. 95)

p(0,t)=H(0,t)=\frac{\mu-\mu e^{-t(\lambda-\mu)}}{\lambda-\mu e^{-t(\lambda-\mu)}}

and its state probabilities are (see Bailey (1964), p. 95)

p(n,t)=(\lambda-\mu)^{2}e^{-t(\lambda-\mu)}\frac{(\lambda-\lambda e^{-t(\lambda-\mu)})^{n-1}}{(\lambda-\mu e^{-t(\lambda-\mu)})^{n+1}},\ n\geq 1.

From (2.14), the probability of ultimate extinction $p(0)$ of LBDP is given by (see Bailey (1964), p. 96)

p(0)=\begin{cases}\frac{\mu}{\lambda},\ \lambda>\mu,\vspace{0.1cm}\\ 1,\ \lambda<\mu.\end{cases}

Remark 2.9.

The moment generating function (mgf) $M(u,t)\coloneqq H(e^{u},t)$ of GLBDP solves

\frac{\partial}{\partial t}M(u,t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}(e^{iu}-1)+\sum_{j=1}^{k_{2}}\mu_{j}(e^{-ju}-1)\right)\frac{\partial}{\partial u}M(u,t),

with $M(u,0)=e^{u}$ .

Next, we obtain the mean and variance of GLBDP.

Theorem 2.1.

Let $\eta=\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}$ . For $t\geq 0$ , the mean $m(t)=\mathbb{E}(N(t))$ and the variance $v(t)=\mathrm{Var}(N(t))$ of GLBDP are given by

m(t)=e^{\eta t}\ \ \text{and}\ \ v(t)=\frac{\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}}{\eta}\left(e^{2\eta t}-e^{\eta t}\right),

(2.15)

respectively.

Proof.

From (2.5), we get $m(0)=\sum_{n=1}^{\infty}np(n,0)=1$ . On multiplying $n$ on both sides of (2.4) and substituting $n(n-i)=(n-i)^{2}+i(n-i)$ and $n(n+j)=(n+j)^{2}-j(n+j)$ yields

\frac{\mathrm{d}}{\mathrm{d}t}m(t)=-\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\right)n^{2}p(n,t)+\sum_{i=1}^{k_{1}}\lambda_{i}(n-i)^{2}p(n-i,t)+\sum_{j=1}^{k_{2}}\mu_{j}(n+j)^{2}p(n+j,t)\\ +\sum_{i=1}^{k_{1}}\lambda_{i}i(n-i)p(n-i,t)-\sum_{j=1}^{k_{2}}\mu_{j}j(n+j)p(n+j,t).

Summing over $n=-k_{2},-k_{2}+1,\ldots$ gives

\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}m(t)=\left(\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}\right)m(t).

(2.16)

On solving (2.16) with initial condition $m(0)=1$ , we get $m(t)=\exp{\big{(}\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}\big{)}t}$ .

Let $m_{2}(t)=\sum_{n=0}^{\infty}n^{2}p(n,t)$ . In order to obtain the variance, we multiply $n^{2}$ on both sides of (2.4) and substitute $n^{2}(n-i)=(n-i)^{3}+2i(n-i)^{2}+i^{2}(n-i)$ and $n^{2}(n+j)=(n+j)^{3}-2j(n+j)^{2}+j^{2}(n+j)$ , to get

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}m_{2}(t)$	$\displaystyle=-\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\right)n^{3}p(n,t)+\sum_{i=1}^{k_{1}}\lambda_{i}(n-i)^{3}p(n-i,t)+2\sum_{i=1}^{k_{1}}\lambda_{i}i(n-i)^{2}p(n-i,t)$
		$\displaystyle\quad+\sum_{i=1}^{k_{1}}\lambda_{i}i^{2}(n-i)p(n-i,t)+\sum_{j=1}^{k_{2}}\mu_{j}(n+j)^{3}p(n+j,t)-2\sum_{j=1}^{k_{2}}\mu_{j}j(n+j)^{2}p(n+j,t)$
		$\displaystyle\quad+\sum_{j=1}^{k_{2}}\mu_{j}j^{2}(n+j)p(n+j,t)$
		$\displaystyle=2\sum_{i=1}^{k_{1}}\lambda_{i}i(n-i)^{2}p(n-i,t)+\sum_{i=1}^{k_{1}}\lambda_{i}i^{2}(n-i)p(n-i,t)-2\sum_{j=1}^{k_{2}}\mu_{j}j(n+j)^{2}p(n+j,t)$
		$\displaystyle\quad+\sum_{j=1}^{k_{2}}\mu_{j}j^{2}(n+j)p(n+j,t).$

By summing over $n=-k_{2},-k_{2}+1,\ldots$ , we get

\frac{\mathrm{d}}{\mathrm{d}t}m_{2}(t)=2\left(\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}\right)m_{2}(t)+\left(\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}\right)m(t).

(2.17)

On solving (2.17) with initial condition $m_{2}(0)=1$ , we get

m_{2}(t)=\left(1+\frac{\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}}{\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}}\right)m(t)^{2}-\frac{\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}}{\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}}m(t).

Hence,

\displaystyle v(t)

\displaystyle=m_{2}(t)-m(t)^{2}=\left(\frac{\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}}{\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}}\right)(m(t)^{2}-m(t)).

This completes the proof. ∎

2.1.1. Parameter estimation in GLBDP

Moran (1951) and Moran (1953) used the inter-event time to estimate the parameters in evolutive processes. Doubleday (1973) used similar technique to obtain the maximum likelihood estimation (MLE) of the parameters in multiple birth-death process. Here, we use the same method to obtain the MLE of the parameter $\Lambda=\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{i=1}^{k_{2}}\mu_{i}$ in GLBDP.

Let $\mathcal{E}$ be the total number of transitions made by $N(t)$ upto some fixed time $\mathcal{T}$ . Suppose $t_{1},t_{2},\ldots,t_{\mathcal{E}}$ are epochs when these transitions have occurred. We define $\tau_{k}=t_{k}-t_{k-1}\ \text{for all}\ 1\leq k\leq\mathcal{E}\ \text{with}\ t_{0}=0,$ where $\tau_{1},\tau_{2},\ldots,\tau_{\mathcal{E}}$ denote the inter-event times. Note that $\tau_{k}$ ’s are independent and exponentially distributed with parameter $\Lambda N(t_{k-1})$ , $k=1,2,\ldots,\mathcal{E}$ . The likelihood function is given by

L(\Lambda;\tau_{1},\tau_{2},\ldots,\tau_{\mathcal{E}})=\Lambda^{n}\prod_{k=1}^{\mathcal{E}}N(t_{k-1})e^{-\Lambda\sum_{k=1}^{\mathcal{E}}N(t_{k-1})\tau_{k}}.

(2.18)

On taking the derivative of log-likelihood function $\ln L(\Lambda;\tau_{1},\tau_{2},\ldots,\tau_{\mathcal{E}})$ with respect to $\Lambda$ , we get

\frac{\mathrm{d}}{\mathrm{d}\Lambda}\ln L(\Lambda;\tau_{1},\tau_{2},\ldots,\tau_{\mathcal{E}})=\frac{\mathcal{E}}{\Lambda}-\sum_{k=1}^{\mathcal{E}}N(t_{k-1})\tau_{k}.

(2.19)

On equating (2.19) to zero, we get the MLE of $\Lambda$ as follows:

\hat{\Lambda}_{mle}=\frac{\mathcal{E}}{\sum_{k=1}^{\mathcal{E}}N(t_{k-1})\tau_{k}}.

Note that $2\Lambda N(t_{k-1})\tau_{k}$ follows $\chi^{2}$ distribution with two degrees of freedom. So,

2\Lambda\sum_{k=1}^{\mathcal{E}}N(t_{k-1})\tau_{k}=2\Lambda\int_{0}^{\mathcal{T}}N(s)\,\mathrm{d}s

has $\chi^{2}$ distribution with $2\mathcal{E}$ degrees of freedom. Thus, $\mathrm{E}\big{(}2\Lambda\int_{0}^{\mathcal{T}}N(s)\,\mathrm{d}s\big{)}=2\mathcal{E}.$ Hence, $\tilde{\Lambda}={\mathcal{E}}^{-1}\int_{0}^{\mathcal{T}}N(s)\,\mathrm{d}s$ is an unbiased estimator of $1/\Lambda$ whose distribution is known. From (2.18), it follows by using the Fisher factorization theorem that $\tilde{\Lambda}={\mathcal{E}}^{-1}\sum_{k=1}^{\mathcal{E}}N(t_{k-1})\tau_{k}$ is a sufficient statistics for $\Lambda$ . Indeed, it is a minimal sufficient statistics for $\Lambda$ .

2.2. GBDP with constant birth and death rates

Here, we obtain the explicit form of the state probabilities of GBDP in the case of constant birth and death rates, that is, $\lambda_{(n)_{i}}=\lambda_{i}>0$ for all $n\geq 0$ and $\mu_{(n)_{j}}=\mu_{j}>0$ for all $n\geq 1$ . In this case, we denote the process by $\{N^{*}(t)\}_{t\geq 0}$ .

From (2.2), the state probabilities $q^{*}(n,t)=\mathrm{Pr}\{N^{*}(t)=n\}$ solve the following system of differential equations:

\frac{\mathrm{d}}{\mathrm{d}t}q^{*}(n,t)=-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\bigg{)}q^{*}(n,t)+\sum_{i=1}^{k_{1}}\lambda_{i}q^{*}(n-i,t)+\sum_{j=1}^{k_{2}}\mu_{j}q^{*}(n+j,t),\ n\geq 0

(2.20)

with initial conditions $q^{*}(1,0)=1$ and $q^{*}(n,0)=0$ for all $n\neq 1$ .

Remark 2.10.

If there is no progenitor at $t=0$ and $\lambda_{(n)_{i}}=\lambda_{i}>0$ , $\mu_{(n)_{j}}=0$ for all $i$ , $j$ and $n\geq 0$ then the GBDP reduces to the generalized counting process (GCP) (see Crescenzo et. al. (2016)).

Remark 2.11.

The pgf $H^{*}(u,t)=\mathbb{E}\left(u^{N^{*}(t)}\right),\ |u|\leq 1$ is the solution of the following differential equation:

\frac{\partial}{\partial t}H^{*}(u,t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}\left(u^{i}-1\right)+\sum_{j=1}^{k_{2}}\mu_{j}\left(u^{-j}-1\right)\right)H^{*}(u,t)

with $H^{*}(u,0)=u$ . It is given by

H^{*}(u,t)=u\exp\left(\sum_{i=1}^{k_{1}}\lambda_{i}\left(u^{i}-1\right)+\sum_{j=1}^{k_{2}}\mu_{j}\left(u^{-j}-1\right)\right)t.

(2.21)

Theorem 2.2.

The state probabilities $q^{*}(n,t)$ of the process $\{N^{*}(t)\}_{t\geq 0}$ are given by

q^{*}(n,t)=\sum_{\Omega^{*}(k_{1},k_{2},n)}\prod_{i=1}^{k_{1}}\prod_{j=1}^{k_{2}}\frac{\lambda_{i}^{x_{i}}\mu_{j}^{y_{j}}t^{x_{i}+y_{j}}}{x_{i}!y_{j}!}e^{-\Lambda t},\ n\geq 0,

where $\Omega^{*}(k_{1},k_{2},n)=\bigg{\{}(x_{1},x_{2},\ldots,x_{k_{1}},y_{1},y_{2},\ldots,y_{k_{2}})\in\mathbb{N}_{0}^{k_{1}+k_{2}}:\sum_{i=1}^{k_{1}}ix_{i}-\sum_{j=1}^{k_{2}}jy_{j}=n-1\bigg{\}}$ and $\Lambda=\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\right)$ .

Proof.

From (2.21), we have

	$\displaystyle H^{*}(u,t)$	$\displaystyle=ue^{-\Lambda t}\sum_{l=0}^{\infty}\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}u^{i}+\sum_{j=1}^{k_{2}}\mu_{j}u^{-j}\bigg{)}^{l}\frac{t^{l}}{l!}$
		$\displaystyle=e^{-\Lambda t}\sum_{l=0}^{\infty}\sum_{S(k_{1},k_{2},l)}\prod_{i=1}^{k_{1}}\prod_{j=1}^{k_{2}}\frac{\lambda_{i}^{x_{i}}\mu_{j}^{y_{j}}t^{x_{i}+y_{j}}}{x_{i}!y_{j}!}u^{ix_{i}-jy_{j}+1},$		(2.22)

where

S(k_{1},k_{2},l)=\bigg{\{}(x_{1},x_{2},\ldots,x_{k_{1}},y_{1},y_{2},\ldots,y_{k_{2}}):x_{i},y_{j}\in\{0,1,\ldots,l\}\ \text{and}\ \sum_{i=1}^{k_{1}}x_{i}+\sum_{j=1}^{k_{2}}y_{j}=l\bigg{\}}.

On rearranging the terms and extracting the coefficients of $u^{n}$ in (2.22), we get the required result. ∎

Theorem 2.3.

For $t\geq 0$ , the mean and the variance of $N^{*}(t)$ are given by

\mathbb{E}(N^{*}(t))=1+\eta t\ \ \text{and}\ \ \mathrm{Var}(N^{*}(t))=\Bigg{(}\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}\Bigg{)}t,

respectively.

Proof.

On taking the derivatives with respect to $u$ on both sides of (2.21), we get

\frac{\partial}{\partial u}H^{*}(u,t)=e^{-\Lambda t}e^{\left(\sum_{i=1}^{k_{1}}\lambda_{i}u^{i}+\sum_{j=1}^{k_{2}}\mu_{j}u^{-j}\right)t}\left(1+\bigg{(}\sum_{i=1}^{k_{1}}i\lambda_{i}u^{i}-\sum_{j=1}^{k_{2}}j\mu_{j}u^{-j}\bigg{)}t\right)

and

	$\displaystyle\frac{\partial^{2}}{\partial u^{2}}H^{*}(u,t)$	$\displaystyle=e^{-\Lambda t}e^{\left(\sum_{i=1}^{k_{1}}\lambda_{i}u^{i}+\sum_{j=1}^{k_{2}}\mu_{j}u^{-j}\right)t}\left(\sum_{i=1}^{k_{1}}i\lambda_{i}u^{i-1}-\sum_{j=1}^{k_{2}}j\mu_{j}u^{-j-1}\right)t$
		$\displaystyle\ \ +\left(\sum_{i=1}^{k_{1}}i(i-1)\lambda_{i}u^{i-2}+\sum_{j=1}^{k_{2}}j(j+1)\mu_{j}u^{-j-2}\right)tH^{*}(u,t)$
		$\displaystyle\ \ +\left(\sum_{i=1}^{k_{1}}i\lambda_{i}u^{i-1}-\sum_{j=1}^{k_{2}}j\mu_{j}u^{-j-1}\right)t\frac{\partial}{\partial u}H^{*}(u,t).$

So, $\mathbb{E}(N^{*}(t))=\frac{\partial}{\partial u}H^{*}(u,t)\bigg{|}_{u=1}=1+\eta t$ and

\mathrm{Var}(N^{*}(t))=\frac{\partial^{2}}{\partial u^{2}}H^{*}(u,t)\bigg{|}_{u=1}+\mathbb{E}(N^{*}(t))-(\mathbb{E}(N^{*}(t)))^{2}=\left(\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}\right)t.

This completes the proof. ∎

3. Cumulative births and deaths in GLBDP

Let $B(t)$ denotes the total number of individuals born by time $t$ .

Proposition 3.1.

The joint pmf $p(b,n,t)=\mathrm{Pr}\{B(t)=b,{N}(t)=n\}$ solves the following system of differential equations:

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p(b,n,t)$	$\displaystyle=-\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\right)np(b,n,t)+\sum_{i=1}^{k_{1}}(n-i)\lambda_{i}p(b-i,n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}(n+j)\mu_{j}p(b,n+j,t),\ b\geq 1,\,n\geq 0,$		(3.1)

with initial conditions

p(b,n,0)=\begin{cases}1,\ b=n=1,\\ 0,\ \text{otherwise},\end{cases}

(3.2)

where $p(b,n,t)=0$ for all $b<0$ or $n<0$ .

Proof.

In an infinitesimal time interval of length $h$ such that $o(h)/h\to 0$ as $h\to 0$ , we have

\mathrm{Pr}\{B(t+h)=b+j,{N}(t+h)=n+k|\,B(t)=b,{N}(t)=n\}\\ =\begin{cases}1-\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{i=1}^{k_{2}}\mu_{i}\right)nh+o(h),\ j=k=0,\vspace{0.1cm}\\ n\lambda_{j}h+o(h),\ j=k=1,2,\ldots,k_{1},\vspace{0.1cm}\\ n\mu_{-k}h+o(h),\ j=0,\,-k=1,2,\ldots,k_{2},\vspace{0.1cm}\\ o(h),\ \mathrm{otherwise}.\end{cases}

So,

	$\displaystyle p(b,n,t+h)$	$\displaystyle=\left(1-\Bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\Bigg{)}nh\right)p(b,n,t)+\sum_{i=1}^{k_{1}}(n-i)\lambda_{i}hp(b-i,n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}(n+j)\mu_{j}hp(b,n+j,t)+o(h).$

Equivalently,

	$\displaystyle\frac{p(b,n,t+h)-p(b,n,t)}{h}$	$\displaystyle=-\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\right)np(b,n,t)+\sum_{i=1}^{k_{1}}(n-i)\lambda_{i}p(b-i,n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}(n+j)\mu_{j}p(b,n+j,t)+\frac{o(h)}{h}.$

On letting $h\to 0$ , we get the required result. ∎

Proposition 3.2.

The pgf $G(u,v,t)=\mathbb{E}\left(u^{B(t)}v^{{N}(t)}\right)$ , $|u|\leq 1$ , $|v|\leq 1$ solves

\frac{\partial}{\partial t}G(u,v,t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}v((uv)^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}v(v^{-j}-1)\right)\frac{\partial}{\partial v}G(u,v,t)

(3.3)

with $G(u,v,0)=uv.$

Proof.

Using (3.2), we get $G(u,v,0)=\sum_{b=0}^{\infty}\sum_{n=0}^{\infty}p(b,n,0)u^{b}v^{n}=uv$ . On multiplying with $u^{b}v^{n}$ on both sides of (3.1) and by adjusting the terms, we get

	$\displaystyle\frac{\partial}{\partial t}G(u,v,t)$	$\displaystyle=-\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\right)v\frac{\partial}{\partial v}u^{b}v^{n}p(b,n,t)+\sum_{i=1}^{k_{1}}\lambda_{i}u^{i}v^{i+1}\frac{\partial}{\partial v}u^{b-i}v^{n-i}p(b-i,n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}\mu_{j}v^{-j+1}\frac{\partial}{\partial v}u^{b}v^{n+j}p(b,n+j,t).$

On summing over $b=0,1,2,\dots$ and $n=-k_{2},-k_{2}+1,\dots$ , we get the required result. ∎

Let $\theta_{u}\leq 0$ and $\theta_{v}\leq 0$ . On taking $u=e^{\theta_{u}}$ and $v=e^{\theta_{v}}$ in (3.3), we get the cumulant generating function (cgf) $K(\theta_{u},\theta_{v},t)=\ln G(e^{\theta_{u}},e^{\theta_{v}},t)$ as a solution of the following differential equation:

\frac{\partial}{\partial t}K(\theta_{u},\theta_{v},t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}\left(e^{i(\theta_{u}+\theta_{v})}-1\right)+\sum_{j=1}^{k_{2}}\mu_{j}\left(e^{-j\theta_{v}}-1\right)\right)\frac{\partial}{\partial\theta_{v}}K(\theta_{u},\theta_{v},t)

(3.4)

with $K(\theta_{u},\theta_{v},0)={\theta_{u}+\theta_{v}}.$

In (3.4), we use the following expanded form of the cgf (see Kendall (1948), Eq. (35)):

K(\theta_{u},\theta_{v},t)=\theta_{u}\mathbb{E}(B(t))+\theta_{v}\mathbb{E}(N(t))+\theta_{u}^{2}/2\mathrm{Var}(B(t))+\theta_{v}^{2}/2\mathrm{Var}(N(t))+\theta_{u}\theta_{v}\mathrm{Cov}(B(t),N(t))+\dots.

(3.5)

to obtain

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\bigg{(}\theta_{u}\mathbb{E}(B(t))$	$\displaystyle+\theta_{v}\mathbb{E}({N}(t))+\frac{\theta_{u}^{2}}{2}\mathrm{Var}(B(t))+\frac{\theta_{v}^{2}}{2}\mathrm{Var}({N}(t))+\theta_{u}\theta_{v}\mathrm{Cov}(B(t),{N}(t))+\cdots\bigg{)}$
	$\displaystyle=\sum_{l=1}^{\infty}\left(\sum_{i=1}^{k_{1}}\lambda_{i}\frac{((\theta_{u}+\theta_{v})i)^{l}}{l!}+\sum_{j=1}^{k_{2}}\mu_{j}\frac{(-\theta_{v}j)^{l}}{l!}\right)\big{(}\mathbb{E}({N}(t))+\theta_{v}\mathrm{Var}({N}(t))$
	$\displaystyle\ \ +\theta_{u}\mathrm{Cov}(B(t),{N}(t))+\cdots\big{)}.$	(3.6)

Proposition 3.3.

For $t\geq 0$ , we have
(i) $\mathbb{E}({N}(t))=\begin{cases}e^{\eta t},\ \eta\neq 0,\vspace{0.1cm}\\ 1,\ \eta=0,\end{cases}$
(ii) $\mathrm{Var}({N}(t))=\begin{cases}\frac{\zeta(e^{2\eta t}-e^{\eta t})}{\eta},\ \eta\neq 0,\vspace{0.1cm}\\ \zeta t,\ \eta=0,\end{cases}$
(iii) $\mathrm{Cov}(B(t),{N}(t))=\begin{cases}\frac{\zeta\sum_{i=1}^{k_{1}}i\lambda_{i}}{\eta^{2}}\left(e^{2\eta t}-e^{\eta t}\right)-\frac{\xi}{\eta}te^{\eta t},\ \eta\neq 0,\vspace{0.1cm}\\ \sum_{i=1}^{k_{1}}i^{2}\lambda_{i}t+\zeta\sum_{i=1}^{k_{1}}i\lambda_{i}\frac{t^{2}}{2},\ \eta=0,\end{cases}$

(iv) $\mathbb{E}(B(t))=\begin{cases}\frac{\sum_{i=1}^{k_{1}}i\lambda_{i}}{\eta}e^{\eta t}-\frac{\sum_{j=1}^{k_{2}}j\mu_{j}}{\eta},\ \eta\neq 0,\vspace{0.1cm}\\ \sum_{i=1}^{k_{1}}i\lambda_{i}t+1,\ \eta=0,\end{cases}$

(v) $\mathrm{Var}(B(t))=\begin{cases}\frac{\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}}{\eta}(e^{\eta t}-1)+2\sum_{i=1}^{k_{1}}i\lambda_{i}\bigg{(}\frac{\zeta\sum_{i=1}^{k_{1}}i\lambda_{i}}{2\eta^{3}}\left(e^{\eta t}-1\right)^{2}-\frac{\xi}{\eta^{3}}\left(\eta te^{\eta t}-e^{\eta t}+1\right)\bigg{)},\ \eta\neq 0,\vspace{0.1cm}\\ \sum_{i=1}^{k_{1}}i^{2}\lambda_{i}t+\sum_{i=1}^{k_{1}}i\lambda_{i}\left(\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}{t^{2}}+\zeta\sum_{i=1}^{k_{1}}i\lambda_{i}\frac{t^{3}}{3}\right),\ \eta=0,\end{cases}$ where $\eta=\sum_{i=1}^{k_{1}}i\lambda_{1}-\sum_{j=1}^{k_{2}}j\mu_{j}$ , $\zeta=\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}$ and $\xi=\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}\sum_{j=1}^{k_{2}}j\mu_{j}+\sum_{i=1}^{k_{1}}i\lambda_{i}\sum_{j=1}^{k_{2}}j^{2}\mu_{j}$ .

Proof.

On comparing the coefficients on the both sides of the (3.6), we get

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}({N}(t))$	$\displaystyle=\begin{cases}\eta\mathbb{E}({N}(t)),\ \eta\neq 0,\vspace{0.1cm}\\ 0,\ \eta=0.\end{cases}$
	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{Var}({N}(t))$	$\displaystyle=\begin{cases}2\eta\mathrm{Var}({N}(t))+\zeta\mathbb{E}({N}(t)),\ \eta\neq 0,\vspace{0.1cm}\\ \zeta,\ \eta=0.\end{cases}$
	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{Cov}(B(t),{N}(t))$	$\displaystyle=\begin{cases}\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}\mathbb{E}({N}(t))+\sum_{i=1}^{k_{1}}i\lambda_{i}\mathrm{Var}({N}(t))+\eta\mathrm{Cov}(B(t),{N}(t)),\ \eta\neq 0,\vspace{0.1cm}\\ \sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{i=1}^{k_{1}}i\lambda_{i}\zeta t,\ \eta=0.\end{cases}$
	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}(B(t))$	$\displaystyle=\begin{cases}\sum_{i=1}^{k_{1}}i\lambda_{i}\mathbb{E}({N}(t)),\ \eta\neq 0,\vspace{0.1cm}\\ \sum_{i=1}^{k_{1}}i\lambda_{i},\ \eta=0,\end{cases}$

and

\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{Var}(B(t))=\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}\mathbb{E}({N}(t))+2\sum_{i=1}^{k_{1}}i\lambda_{i}\mathrm{Cov}(B(t),{N}(t)).

On solving the above differential equations with initial conditions $\mathbb{E}({N}(0))=1$ , $\mathrm{Var}({N}(0))=0$ , $\mathrm{Cov}(B(0),{N}(0))=0$ , $\mathbb{E}(B(0))=1$ and $\mathrm{Var}(B(0))=0$ , respectively, we get the required results. ∎

On using (ii), (iii) and (v) of Proposition 3.3, we can obtain the correlation coefficient between $B(t)$ and $N(t)$ as follows:

R_{B}(t)=\mathrm{Corr}(B(t),N(t))=\frac{\mathrm{Cov}(B(t),N(t))}{\sqrt{\mathrm{Var}(B(t))\mathrm{Var}(N(t))}}.

Remark 3.1.

The mean and variance of GLBDP obtained in Proposition 3.2 coincide with (2.15). From the mean of $N(t)$ and $B(t)$ , we get the following limiting results:

\lim_{t\to\infty}\mathbb{E}({N}(t))=\begin{cases}0,\ \eta<0,\vspace{0.1cm}\\ 1,\ \eta=0,\vspace{0.1cm}\\ \infty,\ \eta>0,\end{cases}

and

\lim_{t\to\infty}\mathbb{E}(B(t))=\begin{cases}\infty,\ \eta\geq 0,\\ -\frac{\sum_{j=1}^{k_{2}}j\mu_{j}}{\eta},\ \eta<0,\end{cases}

respectively. Thus, in the long run, if $\eta>0$ then population is expected to grow beyond any bound, for $\eta<0$ the population is expected to grow up to some finite number then it is expected go extinct, and for $\eta=0$ it is expected that infinitely many people will be born however all of them will die out except one.

In the case of $k_{1}=k_{2}=2$ , some numerical examples of cumulative births are given in Figure 1.

Refer to caption — Figure 1. Cumulative births in GLBDP for $k_{1}=k_{2}=2$ at $t=1$ .

For different values of $\eta$ , the variation of correlation coefficient $R_{B}(t)$ with respect to $t$ is shown in Figure 2. It can be observed that if $\eta>0$ then correlation is high for the larger value of $\eta$ and for significantly large value of $t$ , we have perfect positive correlation between $N(t)$ and $B(t)$ . For $\eta<0$ , $R_{B}(t)$ decreases with time and there is no linear correlation between $N(t)$ and $B(t)$ for the very large value of $t$ .

Let $D(t)$ denotes the total number of deaths by time $t$ .

Proposition 3.4.

The joint pmf $q(d,n,t)=\mathrm{Pr}\{D(t)=d,{N}(t)=n\}$ is the solution of the following system of differential equations:

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}q(d,n,t)$	$\displaystyle=-n\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\right)q(d,n,t)+\sum_{i=1}^{k_{1}}(n-i)\lambda_{i}q(d,n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}(n+j)\mu_{j}q(d-j,n+j,t),\ d\geq 0,\,n\geq 0,$		(3.7)

with initial conditions $q(0,1,0)=1$ and $q(d,n,0)=0$ for all $d\neq 0$ and $n\neq 1$ , where $q(d,n,t)=0$ for all $d<0$ or $n<0$ .

Proof.

In an infinitesimal interval of length $h$ such that $o(h)/h$ as $h\to 0$ , we have

\mathrm{Pr}\{D(t+h)=d+j,{N}(t+h)=n+k|D(t)=d,{N}(t)=n\}\\ =\begin{cases}1-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{i=1}^{k_{2}}\mu_{i}\bigg{)}nh+o(h),\ j=k=0,\vspace{0.1cm}\\ n\lambda_{k}h+o(h),\ j=0,\,k=1,2,\ldots,k_{1},\vspace{0.1cm}\\ n\mu_{j}h+o(h),\ k=-j,\,j=1,2,\ldots,k_{2},\vspace{0.1cm}\\ o(h),\ \mathrm{otherwise}.\end{cases}

So,

	$\displaystyle q(d,n,t+h)$	$\displaystyle=\left(1-\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\right)nh\right)q(d,n,t)+\sum_{i=1}^{k_{1}}(n-i)h\lambda_{i}p(d,n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}(n+j)h\mu_{j}q(d-j,n+j,t)+o(h).$

After rearranging the terms and on letting $h\to 0$ , we get the required result. ∎

The proof of the next result is along the similar line to that of Proposition 3.2, thus it is omitted.

Proposition 3.5.

The pgf $\mathcal{G}(u,v,t)=\mathbb{E}\left(u^{D(t)}v^{{N}(t)}\right)$ , $|u|\leq 1$ , $|v|\leq 1$ solves

\frac{\partial}{\partial t}\mathcal{G}(u,v,t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}v(v^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}v((uv^{-1})^{j}-1)\right)\frac{\partial}{\partial v}\mathcal{G}(u,v,t)

with $\mathcal{G}(u,v,0)=v.$

Thus, the cgf $\mathcal{K}(\theta_{u},\theta_{v},t)=\ln\mathcal{G}(e^{\theta_{u}},e^{\theta_{v}},t)$ , $\theta_{u}\leq 0$ , $\theta_{v}\leq 0$ solves

\frac{\partial}{\partial t}\mathcal{K}(\theta_{u},\theta_{v},t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}(e^{i\theta_{v}}-1)+\sum_{j=1}^{k_{2}}\mu_{j}(e^{j(\theta_{u}-\theta_{v})}-1)\right)\frac{\partial}{\partial\theta_{v}}\mathcal{K}(\theta_{u},\theta_{v},t)

(3.8)

with $\mathcal{K}(\theta_{u},\theta_{v},0)={\theta_{v}}.$

On using (3.5), we get

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\big{(}\theta_{u}\mathbb{E}(D(t))$	$\displaystyle+\theta_{v}\mathbb{E}({N}(t))+\frac{\theta_{u}^{2}}{2}\mathrm{Var}(D(t))+\frac{\theta_{v}^{2}}{2}\mathrm{Var}({N}(t))+\theta_{u}\theta_{v}\mathrm{Cov}(D(t),{N}(t))+\cdots\big{)}$
	$\displaystyle=\sum_{l=1}^{\infty}\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}\frac{(i\theta_{v})^{l}}{l!}+\sum_{j=1}^{k_{2}}\mu_{j}\frac{(j(\theta_{u}-\theta_{v}))^{l}}{l!}\bigg{)}\big{(}\mathbb{E}({N}(t))+\theta_{v}\mathrm{Var}({N}(t))$
	$\displaystyle\ \ +\theta_{u}\mathrm{Cov}(D(t),{N}(t))+\cdots\big{)}$	(3.9)

On comparing the coefficients on the both side of the (3), we get

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}(D(t))$	$\displaystyle=\begin{cases}\sum_{j=1}^{k_{2}}j\mu_{j}\mathbb{E}(N(t)),\ \eta\neq 0,\vspace{0.1cm}\\ \sum_{j=1}^{k_{2}}j\mu_{j},\ \eta=0,\end{cases}$
	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{Cov}(D(t),{N}(t))$	$\displaystyle=\begin{cases}\sum_{j=1}^{k_{2}}\left(-j^{2}\mu_{j}\mathbb{E}({N}(t))+j\mu_{j}\mathrm{Var}({N}(t))\right)+\eta\mathrm{Cov}(D(t),{N}(t)),\ \eta\neq 0,\vspace{0.1cm}\\ \sum_{j=1}^{k_{2}}\left(-j^{2}\mu_{j}+j\mu_{j}\left(\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}\right)t\right),\ \eta=0,\end{cases}$

and

\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{Var}(D(t))=\sum_{j=1}^{k_{2}}j^{2}\mu_{j}\mathbb{E}(N(t))+2\sum_{j=1}^{k_{2}}j\mu_{j}\mathrm{Cov}(D(t),N(t)).

On solving these differential equations with initial conditions $\mathbb{E}(D(0))=0$ , $\mathrm{Cov}(D(0),{N}(0))=0$ and $\mathrm{Var}(D(0))=0$ , respectively, we get the following result.

Proposition 3.6.

For $t\geq 0$ , we have

(i) $\mathbb{E}(D(t))=\begin{cases}\frac{\sum_{j=1}^{k_{2}}j\mu_{j}}{\eta}(e^{\eta t}-1),\ \eta\neq 0,\vspace{0.1cm}\\ \sum_{j=1}^{k_{2}}j\mu_{j}t,\ \eta=0,\end{cases}$
(ii) $\mathrm{Cov}(D(t),{N}(t))=\begin{cases}\frac{\zeta\sum_{j=1}^{k_{2}}j\mu_{j}}{\eta^{2}}\left(e^{2\eta t}-e^{\eta t}\right)-\frac{\xi}{\eta}te^{\eta t},\ \eta\neq 0,\vspace{0.1cm}\\ -\sum_{j=1}^{k_{2}}j^{2}\mu_{j}t+\zeta\sum_{j=1}^{k_{2}}j\mu_{j}\frac{t^{2}}{2},\ \eta=0,\end{cases}$
(iii) $\mathrm{Var}(D(t))=\begin{cases}\frac{\sum_{j=1}^{k_{2}}j^{2}\mu_{j}}{\eta}(e^{\eta t}-1)+2\sum_{j=1}^{k_{2}}j\mu_{j}\bigg{(}\frac{\zeta\sum_{j=1}^{k_{2}}j\mu_{j}}{2\eta^{3}}\left(e^{\eta t}-1\right)^{2}-\frac{\xi}{\eta^{3}}\left(\eta te^{\eta t}-e^{\eta t}+1\right)\bigg{)},\ \eta\neq 0,\vspace{0.1cm}\\ \sum_{j=1}^{k_{2}}j^{2}\mu_{j}t+\sum_{j=1}^{k_{2}}j\mu_{j}\left(-\sum_{j=1}^{k_{2}}j^{2}\mu_{j}{t^{2}}+\zeta\sum_{j=1}^{k_{2}}j\mu_{j}\frac{t^{3}}{3}\right),\ \eta=0.\end{cases}$

Here, again using the results of Proposition 3.2 and Proposition 3.6, we can calculate the correlation coefficient ${R_{D}}(t)=\mathrm{Corr}(D(t),N(t))$ between $D(t)$ and $N(t)$ . For $k_{1}=k_{2}=2$ , some numerical examples of cumulative deaths are given in Figure 3.

From Figure 4, we can observe that for $\eta>0$ , as $t$ increases the correlation coefficient $R_{D}(t)$ changes from negative to positive and it is strong for large values of $t$ . For $\eta<0$ there is no correlation between $N(t)$ and $D(t)$ .

Proposition 3.7.

The joint pmf $p(d,b,n,t)=\mathrm{Pr}\{D(t)=d,B(t)=b,{N}(t)=n\}$ is the solution of the following system of differential equations:

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p(d,b,n,t)$	$\displaystyle=-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\bigg{)}np(d,b,n,t)+\sum_{i=1}^{k_{1}}(n-i)\lambda_{i}p(d,b-i,n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}(n+j)\mu_{j}p(d-j,b,n+j,t),\ d\geq 0,\,b\geq 1,\,n\geq 0,$		(3.10)

with initial conditions $p(0,1,1,0)=1$ and $p(d,b,n,0)=0$ for all $d\neq 0$ , $b\neq 1$ and $n\neq 1$ , where $p(d,b,n,t)=0$ for $d<0$ or $b<0$ or $n<0$ .

Proof.

Let $h$ be the length of an infinitesimal interval such that $o(h)/h\to 0$ as $h\to 0$ . Then,

\mathrm{Pr}\{D(t+h)=d+j,B(t+h)=b+k,{N}(t+h)=n+m|D(t)=d,B(t)=b,{N}(t)=n\}\\ =\begin{cases}1-\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{i=1}^{k_{2}}\mu_{i}\right)nh+o(h),\ j=k=m=0,\vspace{0.1cm}\\ n\lambda_{k}h+o(h),\ j=0,\,k=m=1,2,\ldots,k_{1},\vspace{0.1cm}\\ n\mu_{j}h+o(h),\ k=0,\,m=-j,\,j=1,2,\ldots,k_{2},\vspace{0.1cm}\\ o(h),\ \mathrm{otherwise}.\end{cases}

So,

	$\displaystyle p(d,b,n,t+h)$	$\displaystyle=\bigg{(}1-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\bigg{)}nh\bigg{)}p(d,b,n,t)+\sum_{i=1}^{k_{1}}(n-i)\lambda_{i}hp(d,b-i,n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}(n+j)\mu_{j}hp(d-j,b,n+j,t)+o(h).$

After rearranging the terms and on letting $h\to 0$ , we get the required result. ∎

The pgf $U(u,v,w,t)=\mathbb{E}(u^{D(t)}v^{B(t)}w^{{N}(t)})$ , $|u|\leq 1$ , $|v|\leq 1$ , $|w|\leq 1$ of $(D(t),B(t),N(t))$ satisfies the following differential equation:

\frac{\partial}{\partial t}U(u,v,w,t)=\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}w((vw)^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}w((uw^{-1})^{j}-1)\bigg{)}\frac{\partial}{\partial w}U(u,v,w,t)

with initial condition $U(u,v,w,0)=vw.$

As $N(t)=B(t)-D(t)$ , $t\geq 0$ , we have $\mathrm{Var}(N(t))=\mathrm{Var}(D(t))+\mathrm{Var}(B(t))-2\mathrm{Cov}(D(t),B(t)).$ So, the covariance of $B(t)$ and $D(t)$ , $t\geq 0$ is

	$\displaystyle\mathrm{Cov}(D(t),B(t))$	$\displaystyle=\frac{1}{2}(\mathrm{Var}(D(t))+\mathrm{Var}(B(t))-\mathrm{Var}(N(t)))$
		$\displaystyle=\frac{1}{2}\Bigg{(}\frac{\zeta}{\eta^{3}}\Bigg{(}\bigg{(}\sum_{i=1}^{k_{1}}i\lambda_{i}\bigg{)}^{2}+\bigg{(}\sum_{j=1}^{k_{2}}j\mu_{j}\bigg{)}^{2}\Bigg{)}\left(e^{\eta t}-1\right)^{2}-\frac{\zeta}{\eta}\left(e^{\eta t}-1\right)^{2}$
		$\displaystyle\ \ -\frac{2\xi}{\eta^{3}}\Bigg{(}\sum_{i=1}^{k_{1}}i\lambda_{i}+\sum_{j=1}^{k_{2}}j\mu_{j}\Bigg{)}\left(\eta te^{\eta t}-e^{\eta t}+1\right)\Bigg{)}$
		$\displaystyle=\frac{\zeta\sum_{i=1}^{k_{1}}\sum_{j=1}^{k_{2}}\lambda_{i}\mu_{j}}{\eta^{3}}\left(e^{\eta t}-1\right)^{2}-\frac{\xi}{\eta^{3}}\Bigg{(}\sum_{i=1}^{k_{1}}i\lambda_{i}+\sum_{j=1}^{k_{2}}j\mu_{j}\Bigg{)}\left(\eta te^{\eta t}-e^{\eta t}+1\right),\ \eta\neq 0,$

which can further be used to calculate the correlation coefficient $R(t)=\mathrm{Corr}(D(t),B(t))$ .

Next, we do the similar study for $\{N^{*}(t)\}_{t\geq 0}$ , that is, the GBDP with constant birth and death rates.

3.1. GBDP with constant birth and death rates

Let $B^{*}(t)$ and $D^{*}(t)$ be the total number of births and deaths up to time $t$ in the process $\{N^{*}(t)\}_{t\geq 0}$ , respectively. Similar to the result of Proposition 3.7, the joint pmf $p^{*}(d,b,n,t)=\mathrm{Pr}\{D^{*}(t)=d,B^{*}(t)=b,N^{*}(t)=n\}$ solves the following system of differential equations:

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p^{*}(d,b,n,t)$	$\displaystyle=-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\bigg{)}p^{}(d,b,n,t)+\sum_{i=1}^{k_{1}}\lambda_{i}p^{}(d,b-i,n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}\mu_{j}p^{*}(d-j,b,n+j,t),\ d\geq 0,\,b\geq 1,\,n\geq 0$

with initial conditions $p^{*}(0,1,1,0)=1$ and $p^{*}(d,b,n,t)=0$ for all $d\neq 0$ , $b\neq 1$ and $n\neq 1$ .

So, the pgf $U^{*}(u,v,w,t)=\mathbb{E}(u^{D^{*}(t)}v^{B^{*}(t)}w^{N^{*}(t)})$ , $|u|\leq 1$ , $|v|\leq 1$ , $|w|\leq 1$ is the solution of the following differential equation:

\frac{\mathrm{d}}{\mathrm{d}t}U^{*}(u,v,w,t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}\left((vw)^{i}-1\right)+\sum_{j=1}^{k_{2}}\mu_{j}\left((uw^{-1})^{j}-1\right)\right)U^{*}(u,v,w,t)

with $U^{*}(u,v,w,0)=vw$ . It is given by

U^{*}(u,v,w,t)=vwe^{-\Lambda t}\exp{\left(\sum_{i=1}^{k_{1}}\lambda_{i}(vw)^{i}+\sum_{j=1}^{k_{2}}\mu_{j}(uw^{-1})^{j}\right)t}.

(3.11)

Theorem 3.1.

The joint pmf $p^{*}(d,b,n,t)=\mathrm{Pr}\{D^{*}(t)=d,B^{*}(t)=b,N^{*}(t)=n\}$ is given by

p^{*}(d,b,n,t)=\sum_{\Omega^{*}(k_{1},k_{2},d,b,n)}\prod_{i=1}^{k_{1}}\prod_{j=1}^{k_{2}}\frac{\lambda_{i}^{x_{i}}\mu_{j}^{y_{j}}t^{x_{i}+y_{j}}}{x_{i}!y_{j}!}e^{-\Lambda t},\ d\geq 0,\,b\geq 1,\,n\geq 0,

(3.12)

where $\Omega^{*}(k_{1},k_{2},d,b,n)=\bigg{\{}(x_{1},\ldots,x_{k_{1}},y_{1},\ldots,y_{k_{2}})\in\mathbb{N}_{0}^{k_{1}+k_{2}}:\sum_{j=1}^{k_{2}}jy_{j}=d,\,\sum_{i=1}^{k_{1}}ix_{i}=b-1,\,b-d=n\bigg{\}}.$

Proof.

From (3.11), we have

U^{*}(u,v,w,t)=vwe^{-\Lambda t}\sum_{l=0}^{\infty}\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{i}(vw)^{i}+\sum_{j=1}^{k_{2}}\mu_{j}(uw^{-1})^{j}\bigg{)}^{l}\frac{t^{l}}{l!}.

Equivalently,

U^{*}(u,v,w,t)=e^{-\Lambda t}\sum_{l=0}^{\infty}\sum_{S^{*}(k_{1},k_{2},l)}\prod_{i=1}^{k_{1}}\prod_{j=1}^{k_{2}}\frac{\lambda_{i}^{x_{i}}\mu_{j}^{y_{j}}t^{x_{i}+y_{j}}}{x_{i}!y_{j}!}u^{jy_{j}}v^{ix_{i}+1}w^{ix_{i}-jy_{j}+1},

where $S^{*}(k_{1},k_{2},l)=\big{\{}(x_{1},\ldots,x_{k_{1}},y_{1},\ldots,y_{k_{2}}):x_{i},y_{j}\in\{0,1,\ldots,l\},\,\sum_{i=1}^{k_{1}}x_{i}+\sum_{j=1}^{k_{2}}y_{j}=l\big{\}}.$ So, the coefficients of $u^{d}v^{b}w^{n}$ in $U^{*}(u,v,w,t)$ is the required result. ∎

Remark 3.2.

From (3.12), we get the following joint pmfs

	$\displaystyle p^{*}(b,n,t)$	$\displaystyle=\mathrm{Pr}\{B^{}(t)=b,N^{}(t)=n\}=\sum_{\Omega^{*}(k_{1},k_{2},b,n)}\prod_{i=1}^{k_{1}}\prod_{j=1}^{k_{2}}\frac{\lambda_{i}^{x_{i}}\mu_{j}^{y_{j}}t^{x_{i}+y_{j}}}{x_{i}!y_{j}!}e^{-\Lambda t},\ b\geq 1,\,n\geq 0,$
	$\displaystyle q^{*}(d,n,t)$	$\displaystyle=\mathrm{Pr}\{D^{}(t)=b,N^{}(t)=n\}=\sum_{\Omega^{**}(k_{1},k_{2},d,n)}\prod_{i=1}^{k_{1}}\prod_{j=1}^{k_{2}}\frac{\lambda_{i}^{x_{i}}\mu_{j}^{y_{j}}t^{x_{i}+y_{j}}}{x_{i}!y_{j}!}e^{-\Lambda t},\ d\geq 0,\,n\geq 0,$

and

\tilde{p}^{*}(d,b,t)=\mathrm{Pr}\{D^{*}(t)=d,B^{*}(t)=b\}=\sum_{\Omega^{***}(k_{1},k_{2},d,b)}\prod_{i=1}^{k_{1}}\prod_{j=1}^{k_{2}}\frac{\lambda_{i}^{x_{i}}\mu_{j}^{y_{j}}t^{x_{i}+y_{j}}}{x_{i}!y_{j}!}e^{-\Lambda t},\ d\geq 0,\,b\geq 1,

where

	$\displaystyle\Omega^{*}(k_{1},k_{2},b,n)$	$\displaystyle=\bigg{\{}(x_{1},\ldots,x_{k_{1}},y_{1},\ldots,y_{k_{2}})\in\mathbb{N}_{0}^{k_{1}+k_{2}}:\sum_{i=1}^{k_{1}}ix_{i}=b-1,\,\sum_{i=1}^{k_{1}}ix_{i}-\sum_{j=1}^{k_{2}}jy_{j}=n-1\bigg{\}},$
	$\displaystyle\Omega^{**}(k_{1},k_{2},d,n)$	$\displaystyle=\bigg{\{}(x_{1},\ldots,x_{k_{1}},y_{1},\ldots,y_{k_{2}})\in\mathbb{N}_{0}^{k_{1}+k_{2}}:\sum_{j=1}^{k_{2}}jy_{j}=d,\,\sum_{i=1}^{k_{1}}ix_{i}-\sum_{j=1}^{k_{2}}jy_{j}=n-1\bigg{\}},$

and

\Omega^{***}(k_{1},k_{2},d,b)=\bigg{\{}(x_{1},\ldots,x_{k_{1}},y_{1},\ldots,y_{k_{2}})\in\mathbb{N}_{0}^{k_{1}+k_{2}}:\sum_{j=1}^{k_{2}}jy_{j}=d,\,\sum_{i=1}^{k_{1}}ix_{i}=b-1\bigg{\}}.

Remark 3.3.

From (3.12), we obtain the following marginal pmfs:

\tilde{p}^{*}(d,t)=\mathrm{Pr}\{D^{*}(t)=d\}=\sum_{\Omega^{*}(k_{2},d)}\prod_{j=1}^{k_{2}}\frac{(\mu_{j}t)^{y_{j}}}{y_{j}!}e^{-\sum_{j=1}^{k_{2}}\mu_{j}t},\ d\geq 0,

and

{p}^{*}(b,t)=\mathrm{Pr}\{B^{*}(t)=b\}=\sum_{\Omega^{**}(k_{1},b)}\prod_{i=1}^{k_{1}}\frac{(\lambda_{i}t)^{x_{i}}}{x_{i}!}e^{-\sum_{i=1}^{k_{1}}\lambda_{i}t},\ b\geq 1,

where $\Omega^{*}(k_{2},d)=\bigg{\{}(y_{1},\ldots,y_{k_{2}})\in\mathbb{N}_{0}^{k_{2}}:\sum_{j=1}^{k_{2}}jy_{j}=d\bigg{\}}$ and $\Omega^{**}(k_{1},b)=\bigg{\{}(x_{1},\ldots,x_{k_{1}})\in\mathbb{N}_{0}^{k_{1}}:\sum_{i=1}^{k_{1}}ix_{i}=b-1\bigg{\}}.$

Remark 3.4.

Note that the joint pmf of $(D^{*}(t),B^{*}(t))$ is the product of their marginal pmfs. Hence, ${\{D^{*}(t)\}}_{t\geq 0}$ and ${\{B^{*}(t)\}}_{t\geq 0}$ are two independent random processes.

Proposition 3.8.

For $t\geq 0$ , we have
(i) $\mathbb{E}(B^{*}(t))=1+\sum_{i=1}^{k_{1}}i\lambda_{i}t$ ,
(ii) $\mathbb{E}(D^{*}(t))=\sum_{j=1}^{k_{2}}j\mu_{j}t$ ,
(iii) $\mathrm{Var}(B^{*}(t))=\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}t$ ,
(iv) $\mathrm{Var}(D^{*}(t))=\sum_{j=1}^{k_{2}}j^{2}\mu_{j}t$ .

Proof.

On taking derivatives of (3.11) with respect to $u$ and $v$ , we get

	$\displaystyle\frac{\partial}{\partial u}U^{*}(u,v,w,t)$	$\displaystyle=vwe^{-\Lambda t}e^{\left(\sum_{i=1}^{k_{1}}\lambda_{i}(vw)^{i}+\sum_{j=1}^{k_{2}}\mu_{j}(uw^{-1})^{j}\right)t}\sum_{j=1}^{k_{2}}j\mu_{j}u^{j-1}w^{-j}t,$
	$\displaystyle\frac{\partial^{2}}{\partial u^{2}}U^{*}(u,v,w,t)$	$\displaystyle=vwe^{-\Lambda t}e^{\left(\sum_{i=1}^{k_{1}}\lambda_{i}(vw)^{i}+\sum_{j=1}^{k_{2}}\mu_{j}(uw^{-1})^{j}\right)t}\Bigg{(}\Bigg{(}\sum_{j=1}^{k_{2}}j\mu_{j}u^{j-1}w^{-j}\Bigg{)}^{2}t^{2}$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}j(j-1)\mu_{j}u^{j-2}w^{-j}t\Bigg{)},$
	$\displaystyle\frac{\partial}{\partial v}U^{*}(u,v,w,t)$	$\displaystyle=we^{-\Lambda t}e^{\left(\sum_{i=1}^{k_{1}}\lambda_{i}(vw)^{i}+\sum_{j=1}^{k_{2}}\mu_{j}(uw^{-1})^{j}\right)t}\left(1+\sum_{i=1}^{k_{1}}i\lambda_{i}(vw)^{i}t\right),$

and

	$\displaystyle\frac{\partial}{\partial v^{2}}U^{*}(u,v,w,t)$	$\displaystyle=we^{-\Lambda t}e^{\left(\sum_{i=1}^{k_{1}}\lambda_{i}(vw)^{i}+\sum_{j=1}^{k_{2}}\mu_{j}(uw^{-1})^{j}\right)t}\Bigg{(}\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}v^{i-1}w^{i}t$
		$\displaystyle\ \ +\Bigg{(}1+\sum_{i=1}^{k_{1}}i\lambda_{i}(vw)^{i}t\Bigg{)}\sum_{i=1}^{k_{1}}i\lambda_{i}v^{i-1}w^{i}t\Bigg{)}.$

So,

	$\displaystyle\mathbb{E}(D^{*}(t))$	$\displaystyle=\frac{\partial}{\partial u}U^{*}(u,v,w,t)\bigg{\|}_{u=1,v=1,w=1}=\sum_{j=1}^{k_{2}}j\mu_{j}t,$
	$\displaystyle\mathbb{E}(B^{*}(t))$	$\displaystyle=\frac{\partial}{\partial v}U^{*}(u,v,w,t)\bigg{\|}_{u=1,v=1,w=1}=1+\sum_{i=1}^{k_{1}}i\lambda_{i}t,$
	$\displaystyle\mathrm{Var}(D^{*}(t))$	$\displaystyle=\frac{\partial^{2}}{\partial u^{2}}U^{}(u,v,w,t)\bigg{\|}_{u=1,v=1,w=1}+\mathbb{E}(D^{}(t))-(\mathbb{E}(D^{*}(t)))^{2}=\sum_{j=1}^{k_{2}}j^{2}\mu_{j}t,$

and

\mathrm{Var}(B^{*}(t))=\frac{\partial^{2}}{\partial v^{2}}U^{*}(u,v,w,t)\bigg{|}_{u=1,v=1,w=1}+\mathbb{E}(N^{*}(t))-(\mathbb{E}(B^{*}(t)))^{2}=\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}t.

This completes the proof. ∎

Proposition 3.9.

The variance-covariance matrix $\Sigma$ of $(D^{*}(t),B^{*}(t),N^{*}(t))$ is given by

\Sigma=\begin{pmatrix}\sum_{j=1}^{k_{2}}j^{2}\mu_{j}t&0&-\sum_{j=1}^{k_{2}}j^{2}\mu_{j}t\vspace{0.2cm}\\ 0&\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}t&\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}t\vspace{0.2cm}\\ -\sum_{j=1}^{k_{2}}j^{2}\mu_{j}t&\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}t&\left(\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}\right)t\end{pmatrix}

for all $t\geq 0$ .

Proof.

On taking derivatives of (3.11) with respect to $u,\,v$ and $w$ , we get

	$\displaystyle\frac{\partial^{2}}{\partial u\partial v}U^{*}(u,v,w,t)$	$\displaystyle=we^{-\Lambda t}e^{\left(\sum_{i=1}^{k_{1}}\lambda_{i}(vw)^{i}+\sum_{j=1}^{k_{2}}\mu_{j}(uw^{-1})^{j}\right)t}\left(1+\sum_{i=1}^{k_{1}}i\lambda_{i}(vw)^{i}t\right)\sum_{j=1}^{k_{2}}j\mu_{j}u^{j-1}w^{-j}t,$
	$\displaystyle\frac{\partial^{2}}{\partial u\partial w}U^{*}(u,v,w,t)$	$\displaystyle=ve^{-\Lambda t}e^{\left(\sum_{i=1}^{k_{1}}\lambda_{i}(vw)^{i}+\sum_{j=1}^{k_{2}}\mu_{j}(uw^{-1})^{j}\right)t}\bigg{(}\sum_{j=1}^{k_{2}}j\mu_{j}u^{j-1}w^{-j}t\bigg{(}\sum_{i=1}^{k_{1}}i\lambda_{i}(vw)^{i}$
		$\displaystyle\ \ -\sum_{j=1}^{k_{2}}j\mu_{j}(uw^{-1})^{j}\bigg{)}t-\sum_{j=1}^{k_{2}}j(j-1)\mu_{j}u^{j-1}w^{-j}t\bigg{)},$
	$\displaystyle\frac{\partial^{2}}{\partial v\partial w}U^{*}(u,v,w,t)$	$\displaystyle=e^{-\Lambda t}e^{\left(\sum_{i=1}^{k_{1}}\lambda_{i}(vw)^{i}+\sum_{j=1}^{k_{2}}\mu_{j}(uw^{-1})^{j}\right)t}\Bigg{(}\Bigg{(}1+\sum_{i=1}^{k_{1}}i(i+1)\lambda_{i}(vw)^{i}t\Bigg{)}$
		$\displaystyle\ \ +\Bigg{(}1+\sum_{i=1}^{k_{1}}i\lambda_{i}(vw)^{i}t\Bigg{)}\Bigg{(}\sum_{i=1}^{k_{1}}i\lambda_{i}(vw)^{i}-\sum_{j=1}^{k_{2}}j\mu_{j}(uw^{-1})^{j}\Bigg{)}t\Bigg{)}.$

So,

	$\displaystyle\mathrm{Cov}(D^{}(t),B^{}(t))$	$\displaystyle=\frac{\partial^{2}}{\partial u\partial v}U^{}(u,v,w,t)\bigg{\|}_{u=1,v=1,w=1}-\mathbb{E}(D^{}(t))\mathbb{E}(B^{*}(t))$
		$\displaystyle=\left(1+\sum_{i=1}^{k_{1}}i\lambda_{i}t\right)\sum_{j=1}^{k_{2}}j\mu_{j}t-\left(1+\sum_{i=1}^{k_{1}}i\lambda_{i}t\right)\sum_{j=1}^{k_{2}}j\mu_{j}t=0,$

	$\displaystyle\mathrm{Cov}(D^{}(t),N^{}(t))$	$\displaystyle=\frac{\partial^{2}}{\partial u\partial w}U^{}(u,v,w,t)\bigg{\|}_{u=1,v=1,w=1}-\mathbb{E}(D^{}(t))\mathbb{E}(N^{*}(t))$
		$\displaystyle=\sum_{j=1}^{k_{2}}j\mu_{j}t\left(\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}\right)t-\sum_{j=1}^{k_{2}}j(j-1)\mu_{j}t$
		$\displaystyle\ \ -\sum_{j=1}^{k_{2}}j\mu_{j}t\left(1+\left(\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}\right)t\right)$
		$\displaystyle=-\sum_{j=1}^{k_{2}}j^{2}\mu_{j}t,$

and

	$\displaystyle\mathrm{Cov}(B^{}(t),N^{}(t))$	$\displaystyle=\frac{\partial^{2}}{\partial v\partial w}U^{}(u,v,w,t)\bigg{\|}_{u=1,v=1,w=1}-\mathbb{E}(B^{}(t))\mathbb{E}(N^{*}(t))$
		$\displaystyle=1+\sum_{i=1}^{k_{1}}i(i+1)\lambda_{i}t+\left(1+\sum_{i=1}^{k_{1}}i\lambda_{i}t\right)\left(\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}\right)t$
		$\displaystyle\ \ -\left(1+\sum_{i=1}^{k_{1}}i\lambda_{i}t\right)\left(1+\left(\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}\right)t\right)$
		$\displaystyle=\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}t.$

This completes the proof. ∎

4. Stochastic integral of GBDP

Puri (1966) and Puri (1968) studied the joint distribution of linear birth-death process and its integral. McNeil (1970) obtained the distribution of the integral functional of BDP. Moreover, Gani and McNeil (1971) studied the joint distributions of random variables and their integrals for certain birth-death and diffusion processes. Doubleday (1973) has done similar study for the linear birth-death process with multiple births.

In GBDP, let $\mathcal{N}(0)=k$ , that is, there are $k$ individuals present in a population at time $t=0$ . The hitting time of $\{\mathcal{N}(t)\}_{t\geq}$ to state $0$ given that $\mathcal{N}(0)=k$ is defined as

Z_{k}=\mathrm{inf}\{t\geq 0:\mathcal{N}(t)=0\}.

(4.1)

Let $p_{k}(t)=\mathrm{Pr}\{\mathcal{N}(t)=0|\mathcal{N}(0)=k\}.$ Suppose that there is no possibility of immigration, that is, $\lambda_{(0)_{i}}=0$ for all $i=1,2,\ldots,k_{1}$ . Then, $\mathrm{Pr}\{Z_{k}\leq s\}=p_{k}(s),\ s\geq 0.$

Let us consider an infinitesimal time interval of length $h$ . Then, we have

	$\displaystyle p_{k}(t+h)$	$\displaystyle=\mathrm{Pr}\{\mathcal{N}(h)=k\|\mathcal{N}(0)=k\}p_{k}(t)+\sum_{i=1}^{k_{1}}\mathrm{Pr}\{\mathcal{N}(h)=k+i\|\mathcal{N}(0)=k\}p_{k+i}(t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}\mathrm{Pr}\{\mathcal{N}(h)=k-j\|\mathcal{N}(0)=k\}p_{k-j}(t).$

On using (2.1), we get

\displaystyle p_{k}(t+h)

\displaystyle=\bigg{(}1-\sum_{i=1}^{k_{1}}\lambda_{(k)_{i}}h+\sum_{j=1}^{k_{2}}\mu_{(k)_{j}}h\bigg{)}p_{k}(t)+\sum_{i=1}^{k_{1}}\lambda_{(k)_{i}}hp_{k+i}(t)+\sum_{j=1}^{k_{2}}\mu_{(k)_{j}}hp_{k-j}(t)+o(h),

which on rearranging the terms and letting $h\to 0$ reduces to

\frac{\mathrm{d}}{\mathrm{d}t}p_{k}(t)=-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{(k)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(k)_{j}}\bigg{)}p_{k}(t)+\sum_{i=1}^{k_{1}}\lambda_{(k)_{i}}p_{k+i}(t)+\sum_{j=1}^{k_{2}}\mu_{(k)_{j}}p_{k-j}(t).

(4.2)

Let

\mathcal{W}_{k}=\int_{0}^{Z_{k}}g(\mathcal{N}(t))\,\mathrm{d}t,

(4.3)

where $g(\cdot)$ is a non-negative real valued function on the state-space of GBDP. Observe that if $g(u)=1$ then $\mathcal{W}_{k}=Z_{k}$ and if $g(u)=u$ then $\mathcal{W}_{k}$ is the area under $\mathcal{N}(t)$ up to time of first extinction of the process given that it start from state $k$ .

Now, let us define $\mathcal{W}_{k}(\theta)=\mathbb{E}(e^{-\theta\mathcal{W}_{k}}),\ \theta\geq 0.$ Recall that $T_{s}$ denotes the total time spent by the process $\{\mathcal{N}(t)\}_{t\geq 0}$ in state $s$ before the first transition takes place given that it start from state $s$ , and let $B_{i}$ denotes the event that the first transition is $i$ many births and $D_{j}$ be the event that the first transition is $j$ many deaths given that we start from $k$ . Then, we have

	$\displaystyle\mathcal{W}_{k}(\theta)$	$\displaystyle=\sum_{i=1}^{k_{1}}\mathbb{E}(e^{-\theta\mathcal{W}_{k}}\|B_{i})\mathrm{Pr}(B_{i})+\sum_{j=1}^{k_{2}}\mathbb{E}(e^{-\theta\mathcal{W}_{k}}\|D_{j})\mathrm{Pr}(D_{j})$
		$\displaystyle=\sum_{i=1}^{k_{1}}\mathbb{E}(e^{-\theta(\mathcal{W}_{k+i}+T_{k}g(k))})\mathrm{Pr}(B_{i})+\sum_{j=1}^{k_{2}}\mathbb{E}(e^{-\theta(\mathcal{W}_{k-j}+T_{k}g(k))})\mathrm{Pr}(D_{j})$
		$\displaystyle=\sum_{i=1}^{k_{1}}\mathcal{W}_{k+i}(\theta)\mathbb{E}(e^{-\theta T_{k}g(k)})\mathrm{Pr}(B_{i})+\sum_{j=1}^{k_{2}}\mathcal{W}_{k-j}(\theta)\mathbb{E}(e^{-\theta T_{k}g(k)})\mathrm{Pr}(D_{j}),$

where we used strong Markov property to get the last two steps. From the Markovian property of GBDP, we have $\mathrm{Pr}(B_{i})={\lambda_{(k)_{i}}}/{\rho}$ and $\mathrm{Pr}(D_{j})={\mu_{(k)_{j}}}/{\rho},$ where $\rho=\sum_{i=1}^{k_{1}}\lambda_{(k)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(k)_{j}}$ . On using Remark 2.2, we get

	$\displaystyle\mathcal{W}_{k}(\theta)$	$\displaystyle=\bigg{(}\sum_{i=1}^{k_{1}}\mathcal{W}_{k+i}(\theta)\lambda_{(k)_{i}}+\sum_{j=1}^{k_{2}}\mathcal{W}_{k-j}(\theta)\mu_{(k)_{j}}\bigg{)}\int_{0}^{\infty}e^{-\theta sg(k)}e^{-\rho s}\,\mathrm{d}s$
		$\displaystyle=\bigg{(}\sum_{i=1}^{k_{1}}\mathcal{W}_{k+i}(\theta)\lambda_{(k)_{i}}+\sum_{j=1}^{k_{2}}\mathcal{W}_{k-j}(\theta)\mu_{(k)_{j}}\bigg{)}\frac{1}{\theta g(k)+\rho}.$

Hence,

\theta\mathcal{W}_{k}(\theta)=\sum_{i=1}^{k_{1}}\lambda_{ik}^{*}\mathcal{W}_{k+i}(\theta)+\sum_{j=1}^{k_{2}}\mu_{jk}^{*}\mathcal{W}_{k-j}(\theta)-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{ik}^{*}+\sum_{j=1}^{k_{2}}\mu_{jk}^{*}\bigg{)}\mathcal{W}_{k}(\theta)

(4.4)

with $\mathcal{W}_{0}=\int_{0}^{Z_{0}}g(\mathcal{N}(t))\,\mathrm{d}t=0$ and $\mathcal{W}_{0}(\theta)=1.$ Here, $\lambda_{ik}^{*}={\lambda_{(k)_{i}}}/{g(k)}$ and $\mu_{jk}^{*}={\mu_{(k)_{j}}}/{g(k)}.$

As McNeil (1970) pointed out that Eq. (4.4) correspond to finding the Laplace transform of (4.2) with parameters $\lambda_{ik}^{*}$ and $\mu_{jk}^{*}$ in place of $\lambda_{(k)_{i}}$ and $\mu_{(k)_{j}}$ , respectively, that is, on taking the Laplace transform of (4.2), we get

\theta\mathcal{L}_{k}(\theta)=\sum_{i=1}^{k_{1}}\lambda_{(k)_{i}}\mathcal{L}_{k+i}(\theta)+\sum_{j=1}^{k_{2}}\mu_{(k)_{j}}\mathcal{L}_{k-j}(\theta)-\bigg{(}\sum_{i=1}^{k_{1}}\lambda_{(k)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(k)_{j}}\bigg{)}\mathcal{L}_{k}(\theta),

(4.5)

where

\mathcal{L}_{k}(\theta)=\int_{0}^{\infty}e^{-\theta t}\bigg{(}\frac{\mathrm{d}}{\mathrm{d}t}p_{k}(t)\bigg{)}\,\mathrm{d}t=\mathbb{E}(e^{-\theta Z_{k}})

with $\mathcal{L}_{0}(\theta)=1$ . Hence, the Eq. (4.4) is nothing but the backward equation for GBDP with scaled parameters $\lambda_{ik}^{*}$ and $\mu_{jk}^{*}$ and the distribution of $\mathcal{W}_{k}$ is known whenever the distribution of $Z_{k}$ is known.

Remark 4.1.

For $k_{1}=k_{2}=1$ , (4.4) reduces to

\theta\mathcal{W}_{k}(\theta)=\lambda_{1k}^{*}\mathcal{W}_{k+1}(\theta)+\mu_{1k}^{*}\mathcal{W}_{k-1}(\theta)-(\lambda_{1k}^{*}+\mu_{1k}^{*})\mathcal{W}_{k}(\theta),

with $\mathcal{W}_{0}(\theta)=1$ , which coincide with result obtained in McNeil (1970).

Remark 4.2.

In the cases of GLBDP, that is, $\lambda_{(k)_{i}}=k\lambda_{i}$ and $\mu_{(k)_{j}}=k\mu_{j}$ , (4.4) reduces to

\theta\mathcal{W}_{k}(\theta)=\sum_{i=1}^{k_{1}}\tilde{\lambda}_{ik}^{*}\mathcal{W}_{k+i}(\theta)+\sum_{j=1}^{k_{2}}\tilde{\mu}_{jk}^{*}\mathcal{W}_{k-j}(\theta)-\bigg{(}\sum_{i=1}^{k_{1}}\tilde{\lambda}_{ik}^{*}+\sum_{j=1}^{k_{2}}\tilde{\mu}_{jk}^{*}\bigg{)}\mathcal{W}_{k}(\theta),

where $\tilde{\lambda}_{ik}^{*}={k\lambda_{i}}/{g(k)}$ and $\tilde{\mu}_{jk}^{*}={k\mu_{j}}/{g(k)}.$

4.1. Stochastic path integral of GBDP

Let us consider the stochastic path integral $\mathcal{X}(t)=\int_{0}^{t}g(\mathcal{N}(s))\,\mathrm{d}s$ of GBDP, where $g(\cdot)$ is a non-negative real valued function. Then, conditional on the event $\{\mathcal{N}(0)=m\}$ , the joint distribution function $q_{m}(n,x,t)=\mathrm{Pr}\{\mathcal{N}(t)=n,\mathcal{X}(t)\leq x|\mathcal{N}(0)=m\}$ , $t\geq 0$ of $\mathcal{N}(t)$ and $\mathcal{X}(t)$ solves

	$\displaystyle\frac{\partial}{\partial t}q_{m}(n,x,t)+g(n)\frac{\partial}{\partial x}q_{m}(n,x,t)$	$\displaystyle=-\left(\sum_{i=1}^{k_{1}}\lambda_{(n)_{i}}+\sum_{j=1}^{k_{2}}\mu_{(n)_{j}}\right)q_{m}(n,x,t)+\sum_{i=1}^{k_{1}}\lambda_{(n-i)_{i}}q_{m}(n-i,x,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}\mu_{(n+j)_{j}}q_{m}(n+j,x,t),\ n\geq 0,$		(4.6)

where $q_{m}(n-i,x,t)=0$ for all $i>n$ . This follows from the transitions probabilities of GBDP and its Markov property.

Remark 4.3.

Let $\mathcal{T}_{k}(t)$ be the total time in $(0,t)$ during which the population size is $k$ . Then, for $g(x)=x$ , the stochastic path integral of GBDP can be written as

X(t)=\int_{0}^{t}\mathcal{N}(s)\,\mathrm{d}s=\sum_{k=1}^{\infty}k\mathcal{T}_{k}(t).

Next, we consider the stochastic path integral of GLBDP.

4.1.1. A linear case of GBDP

Here, we obtain some results related to the stochastic path integral of GLBDP. Let $X(t)$ be the stochastic path integral of GLBDP defined as

X(t)=\int_{0}^{t}N(s)\,\mathrm{d}t,\ t\geq 0.

From (4.6), the system of differential equations that governs the joint distribution $p_{m}(n,x,t)=\mathrm{Pr}\{N(t)=n,X(t)<x|N(0)=m\}$ , $t\geq 0$ of $N(t)$ and $X(t)$ is given by

	$\displaystyle\frac{\partial}{\partial t}p_{m}(n,x,t)+n\frac{\partial}{\partial x}p_{m}(n,x,t)$	$\displaystyle=-n\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\right)p_{m}(n,x,t)+\sum_{i=1}^{k_{1}}(n-i)\lambda_{i}p_{m}(n-i,x,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}(n+j)\mu_{j}p_{m}(n+j,x,t),\ n\geq 0.$		(4.7)

Thus, their joint pgf $\tilde{G}(u,v,t)=\int_{0}^{\infty}\sum_{n=0}^{\infty}u^{n}v^{x}p_{m}(n,x,t)\,\mathrm{d}x$ , $|u|\leq 1$ , $|v|\leq 1$ solves

\displaystyle\frac{\partial}{\partial t}\tilde{G}(u,v,t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}(u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}(u^{-j}-1)+\ln v\right)u\frac{\partial}{\partial u}\tilde{G}(u,v,t)

with initial condition $\tilde{G}(u,v,0)=u^{m}$ . Hence, the cgf $\tilde{K}(\theta_{u},\theta_{v},t)=\ln\tilde{G}(e^{\theta_{u}},e^{\theta_{v}},t)$ solves

\frac{\partial}{\partial t}\tilde{K}(\theta_{u},\theta_{v},t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}(e^{i\theta_{u}}-1)+\sum_{j=1}^{k_{2}}\mu_{j}(e^{-j\theta_{u}}-1)+\theta_{v}\right)\frac{\partial}{\partial\theta_{u}}\tilde{K}(\theta_{u},\theta_{v},t)

(4.8)

with $\tilde{K}(\theta_{u},\theta_{v},0)=m\theta_{u}.$

Next, we obtain the mean and variance of the stochastic path integral of GLBDP and its covariance with $N(t)$ .

Proposition 4.1.

For $m=1$ and $t\geq 0$ , we have
(i) $\mathbb{E}(X(t))=\begin{cases}\frac{(e^{\eta t}-1)}{\eta},\ \eta\neq 0,\vspace{0.1cm}\\ t,\ \eta=0,\end{cases}$
(ii) $\mathrm{Cov}(N(t),X(t))=\begin{cases}\frac{\zeta}{\eta^{2}}(e^{2\eta t}-e^{\eta t}-\eta te^{\eta t}),\ \eta\neq 0,\vspace{0.1cm}\\ \zeta\frac{t^{2}}{2},\ \eta=0,\end{cases},$
(iii) $\mathrm{Var}(X(t))=\begin{cases}\frac{2\zeta}{\eta^{2}}\left(\frac{e^{2\eta t}-1}{2\eta}-te^{\eta t}\right),\ \eta\neq 0,\vspace{0.1cm}\\ \zeta\frac{t^{3}}{3},\ \eta=0.\end{cases}$

Proof.

On using (3.5) in (4.8), we have

$\displaystyle\frac{\partial}{\partial t}(\theta_{u}\mathbb{E}(N(t))$	$\displaystyle+\theta_{v}\mathbb{E}(X(t))+\frac{\theta_{u}^{2}}{2}\mathrm{Var}(N(t))+\frac{\theta_{v}^{2}}{2}\mathrm{Var}(X(t))+\theta_{u}\theta_{v}\mathrm{Cov}(N(t),X(t))+\dots)$
	$\displaystyle-\theta_{v}\sum_{l=0}^{\infty}\frac{(-\theta_{u})^{l}}{l!}(\mathbb{E}(N(t))+\theta_{u}\mathrm{Var}(N(t))+\theta_{v}\mathrm{Cov}(N(t),X(t))+\dots)$
	$\displaystyle=\sum_{l=1}^{\infty}\left(\sum_{i=1}^{k_{1}}\lambda_{i}\frac{(i\theta_{u})^{l}}{l!}+\sum_{j=1}^{k_{2}}\mu_{j}\frac{(-j\theta_{u})^{l}}{l!}\right).$	(4.9)

On comparing the coefficients on both sides of (4.9), we get

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}(X(t))$	$\displaystyle=\mathbb{E}(N(t)),$
	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{Cov}(N(t),X(t))$	$\displaystyle=\mathrm{Var}(N(t))+\eta\mathrm{Cov}(N(t),X(t))$

and

\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{Var}(X(t))=2\mathrm{Cov}(N(t),X(t)).

On solving the above differential equations with initial conditions $\mathbb{E}(X(0))=0$ , $\mathrm{Cov}(N(0),X(0))=0$ and $\mathrm{Var}(X(0))=0$ , respectively, we get the required results. ∎

Remark 4.4.

The limiting behavior of the mean of $X(t)$ is given as follows:

\lim_{t\to\infty}\mathbb{E}(X(t))=\begin{cases}-\frac{1}{\eta},\ \eta<0,\vspace{0.1cm}\\ \infty,\ \eta\geq 0.\end{cases}

The correlation coefficient $R_{X}(1)=\mathrm{Corr}(N(1),X(1))$ for different values of parameter are illustrated in Figure 5.

The variation of correlation coefficient between $N(t)$ and its stochastic path integral with time is shown in Figure 6. For $\eta>0$ , $N(t)$ and $X(t)$ have strong positive correlation and as $t$ increases it converges to $+1$ . For $\eta<0$ , the correlation is positive but it is getting weaker with time.

4.2. Stochastic path integral of GBDP in the case of constant birth and death rates

Let $X^{*}(t)=\int_{0}^{t}N^{*}(s)\,\mathrm{d}s$ be the stochastic path integral of the process $\{N^{*}(t)\}_{t\geq 0}$ defined in Section 2.2. From (4.6), the joint distribution $p^{*}_{1}(n,x,t)=\mathrm{Pr}\{N^{*}(t)=n,X^{*}(t)\leq x|N^{*}(0)=1\}$ solves

	$\displaystyle\frac{\partial}{\partial t}p_{1}^{}(n,x,t)+n\frac{\partial}{\partial x}p_{1}^{}(n,x,t)$	$\displaystyle=-\left(\sum_{i=1}^{k_{1}}\lambda_{i}+\sum_{j=1}^{k_{2}}\mu_{j}\right)p_{1}^{}(n,x,t)+\sum_{i=1}^{k_{1}}\lambda_{i}p_{1}^{}(n-i,x,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}\mu_{j}p_{1}^{*}(n+j,x,t),\ n\geq 0$		(4.10)

with $p_{1}^{*}(1,0,0)=1$ , where $p^{*}(n,x,t)=0$ for all $n<0$ .

The pgf $\tilde{G}^{*}(u,v,t)=\int_{0}^{\infty}\sum_{n=0}^{\infty}u^{n}v^{x}{p}_{1}^{*}(n,x,t)\,\mathrm{d}x$ , $0<u<1,\ 0<v<1$ of $(N^{*}(t),X^{*}(t))$ solves the following differential equation:

\frac{\partial}{\partial t}\tilde{G}^{*}(u,v,t)-\ln v\frac{\partial}{\partial u}\tilde{G}^{*}(u,v,t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}(u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}(u^{-j}-1)\right)\tilde{G}^{*}(u,v,t)

(4.11)

with initial condition $\tilde{G}^{*}(u,v,0)=u$ . And, its cgf $\tilde{K}^{*}(\theta_{u},\theta_{v},t)=\ln\tilde{G}^{*}(e^{\theta_{u}},e^{\theta_{v}},t)$ solves

\frac{\partial}{\partial t}\tilde{K}^{*}(\theta_{u},\theta_{v},t)-\theta_{v}e^{-\theta_{u}}\frac{\partial}{\partial\theta_{u}}\tilde{K}^{*}(\theta_{u},\theta_{v},t)=\sum_{i=1}^{k_{1}}\lambda_{i}(e^{i\theta_{u}}-1)+\sum_{j=1}^{k_{2}}\mu_{j}(e^{-j\theta_{u}}-1)

(4.12)

with $\tilde{K}^{*}(\theta_{u},\theta_{v},0)=\theta_{u}.$ On following the similar lines as in the proof of Proposition 4.1, we get

	$\displaystyle\mathbb{E}(X^{*}(t))$	$\displaystyle=t+\frac{\eta t^{2}}{2},$
	$\displaystyle\mathrm{Cov}(N^{}(t),X^{}(t))$	$\displaystyle=t+\left(\eta-\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}\right)\frac{t^{2}}{2}$

and

\mathrm{Var}(X^{*}(t))=t^{2}+\left(\eta-\sum_{i=1}^{k_{1}}i^{2}\lambda_{i}+\sum_{j=1}^{k_{2}}j^{2}\mu_{j}\right)\frac{t^{3}}{3}.

Next, we solve differential equation (4.11) to obtain the joint pgf of $N^{*}(t)$ and $X^{*}(t)$ .

The subsidiary equation corresponding to (4.11) is as follows:

\frac{\partial t}{1}=\frac{\partial u}{-\ln v}=\frac{1}{\sum_{i=1}^{k_{1}}\lambda_{i}(u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}(u^{-j}-1)}\frac{\partial\tilde{G}^{*}}{\tilde{G}^{*}(u,v,t)}.

(4.13)

So, we get

t\ln v+u=\text{constant}

and

u\left(\sum_{i=1}^{k_{1}}\lambda_{i}\left(\frac{u^{i}}{i+1}-1\right)+\sum_{j=1}^{k_{2}}\mu_{j}\left(\frac{u^{-j}}{1-j}-1\right)\right)+\ln v\ln\tilde{G}^{*}(u,v,t)=\text{constant}.

Thus, the general solution of (4.11) is given by

u\left(\sum_{i=1}^{k_{1}}\lambda_{i}\left(\frac{u^{i}}{i+1}-1\right)+\sum_{j=1}^{k_{2}}\mu_{j}\left(\frac{u^{-j}}{1-j}-1\right)\right)+\ln v\ln\tilde{G}^{*}(u,v,t)=\Phi(t\ln v+u),

where $\Phi$ is a real valued function. By using the initial condition $\tilde{G}^{*}(u,v,0)=u$ , we obtain

\Phi(u)=u\left(\sum_{i=1}^{k_{1}}\lambda_{i}\left(\frac{u^{i}}{i+1}-1\right)+\sum_{j=1}^{k_{2}}\mu_{j}\left(\frac{u^{-j}}{1-j}-1\right)\right)+\ln v\ln u.

Thus,

\tilde{G}^{*}(u,v,t)=\exp\left(\frac{\Phi(t\ln v+u)}{\ln v}-\frac{u}{\ln v}\left(\sum_{i=1}^{k_{1}}\lambda_{i}\left(\frac{u^{i}}{i+1}-1\right)+\sum_{j=1}^{k_{2}}\mu_{j}\left(\frac{u^{-j}}{1-j}-1\right)\right)\right).

5. Immigration effect in GLBDP

In GLBDP, the state $n=0$ is the absorbing state of the process, that is, once the population get extinct it can not be revive again. Here, we consider a possibility that the process can revive again once it reaches state $n=0$ . To capture the immigration effect in GLBDP, we introduce a linear version of GBDP, namely, the generalized linear birth-death process with immigration (GLBDPwI).

Let $\nu>0$ be the rate of population growth due to immigration. Here, we consider two different cases of immigration as follows:

5.1. Immigration effect at state $n=0$ only

In this case, the immigration happens only when no alive individual is present in the population. So, the immigration rate is $\lambda_{(0)_{i}}=\nu$ , and once the immigration happens, the birth and death rates are $\lambda_{(n)_{i}}=n\lambda_{i}$ , $i\in\{1,2,\dots,k_{1}\}$ and $\mu_{(n)_{j}}$ , $j\in\{1,2,\dots,k_{2}\}$ , respectively, for all $n\geq 1$ . From (2.2), the state probabilities $p(n,t)$ , $n\geq 0$ of GLBDPwI solve the following system of differential equations:

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p(0,t)$	$\displaystyle=-k_{1}\nu p(0,t)+\sum_{j=1}^{k_{2}}j\mu_{j}p(j,t),$
	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p(n,t)$	$\displaystyle=-\left(\sum_{i=1}^{k_{1}}n\lambda_{i}+\sum_{j=1}^{k_{2}}n\mu_{j}\right)p(n,t)+\nu p(0,t)+\sum_{i=1}^{k_{1}}(n-i)\lambda_{i}p(n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}(n+j)\mu_{j}p(n+j,t),\ 1\leq n\leq k_{1}$

and

\frac{\mathrm{d}}{\mathrm{d}t}p(n,t)=-\left(\sum_{i=1}^{k_{1}}n\lambda_{i}+\sum_{j=1}^{k_{2}}n\mu_{j}\right)p(n,t)+\sum_{i=1}^{k_{1}}(n-i)\lambda_{i}p(n-i,t)+\sum_{j=1}^{k_{2}}(n+j)\mu_{j}p(n+j,t),\ n>k_{1},

with $p(1,0)=1$ and $p(n,t)=0$ for all $n\neq 1$ . Here, $p(n,t)=0$ for all $n<0$ . Its pgf $\mathcal{H}(u,t)=\sum_{n=0}^{\infty}u^{n}p(n,t)$ solves

\frac{\partial}{\partial t}\mathcal{H}(u,t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}u(u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}u(u^{-j}-1)\right)\frac{\partial}{\partial u}\mathcal{H}(u,t)+\sum_{i=1}^{k_{1}}\nu(u^{i}-1)p(0,t)

(5.1)

with $\mathcal{H}(u,0)=u$ .

On taking the derivative of (5.1) with respect to $u$ , we get

	$\displaystyle\frac{\partial^{2}}{\partial u\partial t}\mathcal{H}(u,t)$	$\displaystyle=\left(\sum_{i=1}^{k_{1}}\lambda_{i}u(u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}u(u^{-j}-1)\right)\frac{\partial^{2}}{\partial u^{2}}\mathcal{H}(u,t)$
		$\displaystyle\ \ +\left(\sum_{i=1}^{k_{1}}\lambda_{i}u((i+1)u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}((1-j)u^{-j}-1)\right)\frac{\partial}{\partial u}\mathcal{H}(u,t)+\sum_{i=1}^{k_{1}}i\nu u^{i-1}p(0,t).$		(5.2)

On taking $u=1$ in (5.2), we get the following differential equation that governs the mean $\mathscr{M}(t)=\sum_{n=0}^{\infty}np(n,t)$ of GLBDPwI:

\frac{\mathrm{d}}{\mathrm{d}t}\mathscr{M}(t)=\eta\mathscr{M}(t)+\frac{k_{1}(k_{1}+1)}{2}\nu p(0,t)

with $\mathscr{M}(0)=1$ , where $\eta=\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}.$ Thus,

\mathscr{M}(t)=e^{\eta t}\left(1+\nu\frac{k_{1}(k_{1}+1)}{2}\int_{0}^{t}e^{-\eta s}p(0,s)\,\mathrm{d}s\right).

5.2. Immigration effect at any state

In this case, the immigration effect is always present in the population at a constant rate. So, the birth and death rates are $\lambda_{(n)_{i}}=n\lambda_{i}+\nu$ , $i\in\{1,2,\dots,k_{1}\}$ and $\mu_{(n)_{j}}$ , $j\in\{1,2,\dots,k_{2}\}$ , respectively, for all $n\geq 0$ . From (2.2), the state probabilities $p(n,t)$ of GLBDPwI solve

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p(n,t)$	$\displaystyle=-\left(k_{1}\nu+n\sum_{i=1}^{k_{1}}\lambda_{i}+n\sum_{j=1}^{k_{2}}\mu_{j}\right)p(n,t)+\sum_{i=1}^{k_{1}}(\nu+(n-i)\lambda_{i})p(n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{k_{2}}(n+j)\mu_{j}p(n+j,t),\ n\geq 0,$

where $p(1,0)=1$ and $p(n,t)=0$ for all $n\neq 1$ . So, the differential equation that governs the pgf $\mathcal{H}(u,t)=\sum_{n=0}^{\infty}u^{n}p(n,t)$ of GLBDPwI is given by

\frac{\partial}{\partial t}\mathcal{H}(u,t)=\left(\sum_{i=1}^{k_{1}}\lambda_{i}u(u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}u(u^{-j}-1)\right)\frac{\partial}{\partial u}\mathcal{H}(u,t)+\left(\sum_{i=1}^{k_{1}}\nu(u^{i}-1)\right)\mathcal{H}(u,t)

(5.3)

with $\mathcal{H}(u,0)=u$ .

On taking the derivative with respect to $u$ on both sides of (5.3), we have

$\displaystyle\frac{\partial}{\partial t}\frac{\partial}{\partial u}\mathcal{H}(u,t)$	$\displaystyle=\left(\sum_{i=1}^{k_{1}}\lambda_{i}u(u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}u(u^{-j}-1)\right)\frac{\partial^{2}}{\partial u^{2}}\mathcal{H}(u,t)$
	$\displaystyle\ \ +\left(\sum_{i=1}^{k_{1}}\lambda_{i}((i+1)u^{i}-1)+\sum_{j=1}^{k_{2}}\mu_{j}((-j+1)u^{-j}-1)\right)\frac{\partial}{\partial u}\mathcal{H}(u,t)+\sum_{i=1}^{k_{1}}i\nu u^{i-1}\mathcal{H}(u,t)$
	$\displaystyle\ \ +\sum_{i=1}^{k_{1}}\nu(u^{i}-1)\frac{\partial}{\partial u}\mathcal{H}(u,t).$	(5.4)

On taking $u=1$ in (5.4), we get the mean $\mathscr{M}(t)$ of GLBDPwI as a solution of the following differential equation:

\frac{\partial}{\partial t}\mathscr{M}(t)=\eta\mathscr{M}(t)+\frac{k_{1}(k_{1}+1)\nu}{2}

with initial condition $\mathscr{M}(0)=1$ . It is given by

\mathscr{M}(t)=\begin{cases}e^{\eta t}+\frac{k_{1}(k_{1}+1)\nu}{2\eta}(e^{\eta t}-1),\ \eta\neq 0,\vspace{0.1cm}\\ 1+\frac{k_{1}(k_{1}+1)\nu t}{2},\ \eta=0,\end{cases}

where $\eta=\sum_{i=1}^{k_{1}}i\lambda_{i}-\sum_{j=1}^{k_{2}}j\mu_{j}$ .

Remark 5.1.

For $\nu=0$ , the mean of GLBDPwI reduces to the mean of GLBDP given in (2.15).

For different values of $\eta$ the variation of expected values of GLBDP and GLBDPwI with time is shown in Figure 7.

In GLBDPwI, for $\eta>0$ , the population size grow exponentially fast for the large value of $t$ and for $\eta<0$ , the limiting mean value for large $t$ is $-k_{1}(k_{1}+1)\nu/2\eta$ , that is,

\lim_{t\to\infty}\mathscr{M}(t)=\begin{cases}\infty,\ \eta\geq 0,\vspace{0.1cm}\\ -\frac{k_{1}(k_{1}+1)\nu}{2\eta},\ \eta<0.\end{cases}

Thus, as expected in GLBDPwI, the population will never go extinct.

Remark 5.2.

On comparing (2.6) and (5.3), we note that they differs in the term $\sum_{i=1}^{k_{1}}\nu(u^{i}-1)\mathcal{H}(u,t)$ . It is a generalized counting process component with constant transition rates. Suppose that the birth rates in GLBDPwI tend to zero, that is, $\lambda_{i}\to 0$ for all $i\in\{1,2,\dots,k_{1}\}$ then it can be used to study the fluctuation in the number of particles enclosed inside very small volume. Here, the immigration of particle from surrounding happens as GCP with constant rates and the emigration of finitely many particles from closed volume is considered as death whose rates are linearly proportional to the number of particle present.

6. Application to vehicles parking management system

Suppose a vehicles parking lot has maximum capacity of $K$ parking spots. If we assume that the gate of the parking area is wide enough such that multiple vehicles can arrive or depart together and the chance of simultaneous arrival and departure of the vehicles is negligible. Let $N(t)$ denotes the total number of vehicles parked at time $t$ . Then, we can use a linear version of GBDP to model this problem, where $\{N(t)\}_{t\geq 0}$ is a GBDP with parameters

\lambda_{(n)_{i}}=(K-n)\lambda_{i},\ i\in\{1,2,\ldots,K_{1}\}\ \text{and}\ \mu_{(n)_{j}}=n\mu_{j},\ j\in\{1,2,\ldots,K_{2}\}

for all $0\leq n\leq K$ , $K_{1}<K$ and $K_{2}<K$ . Thus, the rates of arrival and departure of vehicles are proportional to the number of vacant spots and the total number of vehicles parked in the parking lot, respectively. Let $A(t)$ and $D(t)$ denote the number of vehicles arrived and departed by time $t$ . Here, we are interested in knowing the quantities $\mathbb{E}(N(t))$ , $\mathbb{E}(A(t))$ and $\mathbb{E}(D(t))$ . Let us assume that there are no vehicles parked at time $t=0$ . As obtained in the case of GLBDP, it can be shown that the pmf $p(n,t)=\mathrm{Pr}\{N(t)=n\}$ solves the following system of differential equations:

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p(n,t)$	$\displaystyle=-\bigg{(}\sum_{i=1}^{K_{1}}(K-n)\lambda_{i}+\sum_{j=1}^{K_{2}}n\mu_{j}\bigg{)}p(n,t)+\sum_{i=1}^{K_{1}}(K-n+i)\lambda_{i}p(n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{K_{2}}(n+j)\mu_{j}p(n+j,t),\ 0\leq n\leq K,$

with initial condition $p(0,0)=1$ . Here, $p(n,t)=0$ for all $n<0$ and $n>K$ . Hence, the governing system of differential equation for its mean is given by

\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}(N(t))=-\left(\sum_{i=1}^{K_{1}}i\lambda_{i}+\sum_{j=1}^{K_{2}}j\mu_{j}\right)\mathbb{E}(N(t))+K\sum_{i=1}^{K_{1}}i\lambda_{i}

(6.1)

with initial condition $\mathbb{E}(N(0))=0$ . On solving (6.1), we get

\mathbb{E}(N(t))=\frac{K\sum_{i=1}^{K_{1}}i\lambda_{i}}{\beta}\left(1-e^{-\beta t}\right),

where $\beta=\sum_{i=1}^{K_{1}}i\lambda_{i}+\sum_{j=1}^{K_{2}}j\mu_{j}$ . In the long run, the expected number of vehicles parked is given by

\lim_{t\to\infty}\mathbb{E}(N(t))=\frac{K\sum_{i=1}^{K_{1}}i\lambda_{i}}{\beta}<K.

The joint pmf $p(a,n,t)=\mathrm{Pr}\{A(t)=a,N(t)=n\}$ solves

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p(a,n,t)$	$\displaystyle=-\left(\sum_{i=1}^{K_{1}}(K-n)\lambda_{i}+\sum_{j=1}^{K_{2}}n\mu_{j}\right)p(a,n,t)+\sum_{i=1}^{K_{1}}(K-n+i)\lambda_{i}p(a-i,n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{K_{2}}(n+j)\mu_{j}p(a,n+j,t),\ a,n\in\{0,1,\dots,K\}$

with $p(0,0,0)=1$ . Its proof follows similar lines to that of Proposition 3.1.

The joint cgf $K(\theta_{u},\theta_{v},t)=\ln\mathbb{E}(e^{\theta_{u}A(t)+\theta_{v}N(t)})$ satisfies the following differential equation:

\frac{\partial}{\partial t}K(\theta_{u},\theta_{v},t)+\left(\sum_{i=1}^{K_{1}}\lambda_{i}(e^{i(\theta_{u}+\theta_{v})}-1)-\sum_{j=1}^{K_{2}}\mu_{j}(e^{-j\theta_{v}}-1)\right)\frac{\partial}{\partial\theta_{v}}K(\theta_{u},\theta_{v},t)=K\sum_{i=1}^{K_{1}}\lambda_{i}(e^{i(\theta_{u}+\theta_{v})}-1)

with initial condition $K(\theta_{u},\theta_{v},0)=0$ . On using (3.5), we get

\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}(A(t))=\sum_{i=1}^{K_{1}}i\lambda_{i}(K-\mathbb{E}(N(t)))

with $\mathbb{E}(A(0))=0$ . Thus,

\mathbb{E}(A(t))=\sum_{i=1}^{K_{1}}i\lambda_{i}\left(Kt-\frac{K\sum_{i=1}^{K_{1}}i\lambda_{i}}{\beta}\bigg{(}t+\frac{e^{-t\beta}-1}{\beta}\bigg{)}\right).

Similarly, the joint pmf $q(d,n,t)=\mathrm{Pr}\{D(t)=d,N(t)=n\}$ is the solution of following system of differential equations:

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}q(d,n,t)$	$\displaystyle=-\left(\sum_{i=1}^{K_{1}}(K-n)\lambda_{i}+\sum_{j=1}^{K_{2}}n\mu_{j}\right)q(d,n,t)+\sum_{i=1}^{K_{1}}(K-n+i)\lambda_{i}q(d,n-i,t)$
		$\displaystyle\ \ +\sum_{j=1}^{K_{2}}(n+j)\mu_{j}q(d-j,n+j,t),\ d,n\in\{0,1,\dots,K\}.$

Hence, the cgf $\mathcal{K}(\theta_{u},\theta_{v},t)=\ln\mathbb{E}(e^{\theta_{u}D(t)+\theta_{v}N(t)})$ solves

\frac{\partial}{\partial t}\mathcal{K}(\theta_{u},\theta_{v},t)+\left(\sum_{i=1}^{K_{1}}\lambda_{i}(e^{i\theta_{v}}-1)+\sum_{j=1}^{K_{2}}\mu_{j}(e^{j(\theta_{u}-\theta_{v})}-1)\right)\frac{\partial}{\partial\theta_{v}}\mathcal{K}(\theta_{u},\theta_{v},t)=K\sum_{i=1}^{K_{1}}\lambda_{i}(e^{i\theta_{v}}-1)

with $\mathcal{K}(\theta_{u},\theta_{v},0)=0$ . Thus, we get

\mathbb{E}(D(t))=K\sum_{i=1}^{K_{1}}i\lambda_{i}\sum_{j=1}^{K_{2}}j\mu_{j}\left(\frac{t}{\beta}+\frac{e^{-\beta t}-1}{\beta^{2}}\right).

Let us consider the stochastic path integral $X(t)=\int_{0}^{t}N(s)\,\mathrm{d}s$ . Then, we can find the average occupancy $O(t)$ of the parking lot at any given time $t\geq 0$ which we defined as $O(t)=t^{-1}\mathbb{E}(X(t))$ .

From (4.6), the system of differential equations that governs the joint distribution $p_{0}(n,x,t)$ is given by

	$\displaystyle\frac{\partial}{\partial t}p_{0}(n,x,t)+n\frac{\partial}{\partial x}p_{0}(n,x,t)$	$\displaystyle=-\left(\sum_{i=1}^{K_{1}}(K-n)\lambda_{i}+\sum_{j=1}^{K_{2}}n\mu_{j}\right)p_{0}(n,x,t)$
		$\displaystyle\ \ +\sum_{i=1}^{K_{1}}(K-n+i)\lambda_{i}p_{0}(n-i,x,t)$
		$\displaystyle\ \ +\sum_{j=1}^{K_{2}}(n+j)\mu_{j}p_{0}(n+j,x,t),\ 0\leq n\leq K,$

where $p_{0}(n,x,t)=0$ for all $n<0$ and $n>K$ .

So, the joint cgf $\tilde{K}(\theta_{u},\theta_{v},t)=\ln\mathbb{E}(e^{\theta_{u}N(t)+\theta_{v}X(t)})$ solves

\frac{\partial}{\partial t}\tilde{K}(\theta_{u},\theta_{v},t)+\left(\sum_{i=1}^{K_{1}}\lambda_{i}(e^{i\theta_{u}}-1)-\sum_{j=1}^{K_{2}}\mu_{j}(e^{-j\theta_{u}}-1)-\theta_{v}\right)\frac{\partial}{\partial\theta_{u}}\tilde{K}(\theta_{u},\theta_{v},t)=K\sum_{i=1}^{K_{1}}\lambda_{i}(e^{i\theta_{u}}-1)

(6.2)

with $\tilde{K}(\theta_{u},\theta_{v},0)=0$ .

On using (3.5), we get the following differential equation:

\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}(X(t))=\mathbb{E}(N(t))

with $\mathbb{E}(X(0))=0$ . Hence,

\mathbb{E}(X(t))=\frac{K\sum_{i=1}^{K_{1}}i\lambda_{i}}{\beta}\left(t+\frac{e^{-\beta t}-1}{\beta}\right).

Thus, average occupancy is

O(t)=\frac{K\sum_{i=1}^{K_{1}}i\lambda_{i}}{\beta}\left(1+\frac{e^{-\beta t}-1}{\beta t}\right).

References

[1] Bailey, N.T.J., 1964. The Elements of Stochastic Processes with Applications to the Natural Sciences. Wiley, New York.
[2] Doubleday, W.G., 1973. On linear birth-death processes with multiple births. Math Biosci. 17, 43-56.
[3] Di Crescenzo, A., Martinucci, B. and Meoli, A. 2016. A fractional counting process and its connection with the Poisson process. ALEA Lat. Am. J. Probab. Math. Stat. 13(1), 291-307.
[4] Feller, W., 1968. An Introduction to Probability Theory and Its Applications, Vol. 1, 3rd ed. Wiley, New York.
[5] Gani, J., McNeil, D.R., 1971. Joint distributions of random variables and their integrals for certain birth-death and diffusion processes. Advances Appl. Prob. 3, 339-352.
[6] Kendall, D.G., 1948. On the Generalized ”Birth-and-Death” Process. Ann. Math. Statist. 19(1), 1-15.
[7] Karlin, S., Taylor, H.M., 1975. A first course in stochastic processes. Academic Press.
[8] Moran, P.A.P, (1951). Estimation methods for evolutive processes. J. R. Statist. Soc. (B)13, 141-146.
[9] Moran, P.A.P, (1953). Estimation of parameters of a birth and death process. J. R. Statist. Soc. (B)15, 241-245.
[10] McNeil, D.R., 1970. Integral functions of birth and death processes and related limiting distributions. Ann. Math. Stat. 41, 480-485.
[11] Puri, P.S., 1966. On the homogeneous birth and death process and its integral. Biometrika. 53, 61-67.
[12] Puri, P.S., 1968. Some further results on the birth-and-death process and its integral. Proc. Camb. Phil. Sot. 67, 141-154.

On the Generalized Birth-Death Process and its Linear Versions

Abstract.

Key words and phrases:

2010 Mathematics Subject Classification:

1. Introduction

2. Generalized birth-death process

Remark 2.1.

Proposition 2.1.

Proof.

Remark 2.2.

2.1. A linear case of GBDP

Remark 2.3.

Remark 2.4.

Proposition 2.2.

Proof.

Remark 2.5.

Remark 2.6.

Remark 2.7.

Remark 2.8.

Remark 2.9.

Theorem 2.1.

Proof.

2.1.1. Parameter estimation in GLBDP

2.2. GBDP with constant birth and death rates

Remark 2.10.

Remark 2.11.

Theorem 2.2.

Proof.

Theorem 2.3.

Proof.

3. Cumulative births and deaths in GLBDP

Proposition 3.1.

Proof.

Proposition 3.2.

Proof.

Proposition 3.3.

Proof.

Remark 3.1.

Proposition 3.4.

Proof.

Proposition 3.5.

Proposition 3.6.

Proposition 3.7.

Proof.

3.1. GBDP with constant birth and death rates

Theorem 3.1.

Proof.

Remark 3.2.

Remark 3.3.

Remark 3.4.

Proposition 3.8.

Proof.

Proposition 3.9.

Proof.

4. Stochastic integral of GBDP

Remark 4.1.

Remark 4.2.

4.1. Stochastic path integral of GBDP

Remark 4.3.

4.1.1. A linear case of GBDP

Proposition 4.1.

Proof.

Remark 4.4.

4.2. Stochastic path integral of GBDP in the case of constant birth and death rates

5. Immigration effect in GLBDP

5.1. Immigration effect at state n=0n=0 only

5.2. Immigration effect at any state

Remark 5.1.

Remark 5.2.

6. Application to vehicles parking management system

References

5.1. Immigration effect at state $n=0$ only