# A computational approach to extreme values and related hitting probabilities
in level-dependent quasi-birth-death processes

A. Di Crescenzo<EMAIL_ADDRESS>[ A. Gómez-Corral<EMAIL_ADDRESS>[ D.
Taipe<EMAIL_ADDRESS>Dipartimento di Matematica, Università degli Studi di
Salerno, Via Giovanni Paolo II, 132, Fisciano, I-84084, Italy  Departamento de
Estadística e Investigación Operativa, Facultad de Ciencias Matemáticas,
Universidad Complutense de Madrid, Plaza de Ciencias 3, 28040 Madrid, Spain
Departamento de Estadística y Ciencia de Datos, Facultad de Estudios
Estadísticos, Universidad Complutense de Madrid, Avda. Puerta de Hierro s/n,
28040 Madrid, Spain

###### Abstract

This paper analyzes the dynamics of a level-dependent quasi-birth-death
process ${\cal X}=\\{(I(t),J(t)):t\geq 0\\}$, i.e., a bi-variate Markov chain
defined on the countable state space $\cup_{i=0}^{\infty}l(i)$ with
$l(i)=\\{(i,j):j\in\\{0,...,M_{i}\\}\\}$, for integers
$M_{i}\in\mathbb{N}_{0}$ and $i\in\mathbb{N}_{0}$, which has the special
property that its $q$-matrix has a block-tridiagonal form. Under the
assumption that the first passage to the subset $l(0)$ occurs in a finite time
with certainty, we characterize the probability law of
$(\tau_{\max},I_{\max},J(\tau_{\max}))$, where $I_{\max}$ is the running
maximum level attained by process ${\cal X}$ before its first visit to states
in $l(0)$, $\tau_{\max}$ is the first time that the level process
$\\{I(t):t\geq 0\\}$ reaches the running maximum $I_{\max}$, and
$J(\tau_{\max})$ is the phase at time $\tau_{\max}$. Our methods rely on the
use of restricted Laplace-Stieltjes transforms of $\tau_{\max}$ on the set of
sample paths $\\{I_{\max}=i,J(\tau_{\max})=j\\}$, and related processes under
taboo of certain subsets of states. The utility of the resulting computational
algorithms is demonstrated in two epidemic models: the SIS model for
horizontally and vertically transmitted diseases; and the SIR model with
constant population size.

###### keywords:

epidemic model , first-passage time , hitting probability , quasi-birth-death
process

###### MSC:

[2010] 60J28 , 92B05

††journal: Mathematics and Computers in Simulation

url]https://docenti.unisa.it/005305/home

url]blogs.mat.ucm.es/agomez-corral

##  1 Introduction

In this work, we are interested in continuous-time Markov chains ${\cal
X}=\\{(I(t),J(t)):t\geq 0\\}$ defined on a countable state space ${\cal S}$,
which we partition as $\cup_{i=0}^{\infty}l(i)$ with subsets
$l(i)=\\{(i,j):j\in\\{0,...,M_{i}\\}\\}$, for integers
$M_{i}\in\mathbb{N}_{0}$ and $i\in\mathbb{N}_{0}$. The $q$-matrix of ${\cal
X}$ is assumed to be conservative and irreducible, and exhibits the structured
form

$\displaystyle Q$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}Q_{0,0}&Q_{0,1}&&&&\\\
Q_{1,0}&Q_{1,1}&Q_{1,2}&&&\\\ &\ddots&\ddots&\ddots&&\\\
&&Q_{i,i-1}&Q_{i,i}&Q_{i,i+1}&\\\ &&&\ddots&\ddots&\ddots\end{array}\right),$
(6)

where the sub-matrices $Q_{i,i^{\prime}}$ have entries
$q_{(i,j),(i^{\prime},j^{\prime})}$, for pairs $(i,j)$ with
$i\in\mathbb{N}_{0}$ and $j\in\\{0,...,M_{i}\\}$, and
$(i^{\prime},j^{\prime})$ with $i^{\prime}\in\\{\max\\{0,i-1\\},i,i+1\\}$ and
$j^{\prime}\in\\{0,...,M_{i^{\prime}}\\}$. This means that process ${\cal X}$
can move in one step from state $(i,j)$, for $i\in\mathbb{N}_{0}$ and
$j\in\\{0,...,M_{i}\\}$, to state $(i^{\prime},j^{\prime})$ only if
$|i^{\prime}-i|\leq 1$ and $j^{\prime}\in\\{0,...,M_{i^{\prime}}\\}$, and,
consequently, the values of $I(t)$ can only change by $-1$, $0$ and $1$ units
at jump times of ${\cal X}$. Process ${\cal X}$ is commonly referred to as a
level-dependent quasi-birth-death (LD-QBD) process, with the first variable
$I(t)$ termed the _level_ , and the second one $J(t)$ the _phase_ ; for
convenience, we also refer to the subset $l(i)$ as the $i$th _level_. Markov
chains with LD-QBD structure can be thought of as the natural generalization
of the uni-variate birth-death process (see, e.g., Anderson [3, Chapter 8]),
competition processes (Iglehart [29]; Reuter [46]) and related models,
including prey-predator processes (Hitchcock [27]; Ridle-Rowe [48]), two-
species competition processes (Gómez-Corral and López García [20]; Ridler-Rowe
[47]), and birth/birth-death processes with biological applications (Ho et al.
[28]).

LD-QBD processes have been extensively investigated by several researchers,
among whom are Baumann and Sandmann [5, 6], Bright and Taylor [10], Phung-Duc
et al. [42], and Takine [53]. Readers are also referred to Chapter 12 of
Latouche and Ramaswami [36], the overviews of Bright and Taylor [11], and
Kharoufeh [33], and references therein. The focus in these references is
mainly on the stationary vector of a regular process ${\cal X}$ and its
expression in terms of _rate_ matrices $\\{R_{i}:i\in\mathbb{N}_{0}\\}$, where
$R_{i}$ consists of expected rates of visit to state $(i+1,j^{\prime})$ in
level $l(i+1)$ per unit of the local time of $l(i)$ spent in $(i,j)$. To be
concrete, the solution of Bright and Taylor [10] is developed as an extension
of previous matrix-geometric methods of Neuts [40] and the logarithmic
reduction solution in Ref. [35] —which are both related to a level-independent
QBD process— for the level-dependent case. From the probabilistic
interpretation of $R_{i}$ and a suitable continued fraction representation,
Phung-Duc et al. [42] provide a novel algorithmic procedure to evaluate the
underlying rate matrices, which is found to be more efficient than that in
Ref. [10]; for a related work, see Baumann and Sandmann [5]. In Ref. [6],
Baumann and Sandmann present a matrix-analytic solution for numerically
computing stationary expectations for long-run averages of additive
functionals, such as long-run average costs and reward rates, without at first
computing the stationary vector of ${\cal X}$. The approach of Baumann [8]
(see also Baumann and Sandmann [6]) relies on truncating the state space
${\cal S}$ at a suitably selected level $l(N)$, which results in inevitable
errors for which lower and upper bounds on stationary expectations are
efficiently computed; see Ref. [8, Section 4]. QBD processes on a finite state
space have been extensively studied in the literature, specifically in Refs.
[1, 17, 19, 22, 23], among others. In the setting of LD-QBD processes with an
_explosive_ state space (i.e., as the cardinality $1+M_{i}$ of level $l(i)$
increases exponentially with $i$), Takine [53] presents approximation
techniques for computing performance measures. For the analysis of quasi-
stationary regime, see the work of Bean et al. [9], and Ramaswami and Taylor
[45].

In this paper, the focus is on the extreme values of the LD-QBD process ${\cal
X}$ before the first entrance to states in $l(0)$, and hitting probabilities
of the underlying phase process $\\{J(t):t\geq 0\\}$ as the maximum level of
${\cal X}$ is visited for the first time. A first motivation to study extreme
values and related hitting probabilities in a LD-QBD process comes from
applications. For instance, the structured form of the $q$-matrix in (6) is
often obtained directly in the case of retrial queues (Artalejo and Gómez-
Corral [4]), in such a way that the length of a busy period is equivalent to
the time to reach level $l(0)$ in a suitably defined LD-QBD process; for a
variety of retrial queues, see the papers by Dayar and Orhan [16], Jeganathan
et al. [31], Liu and Zhao [38], Phung-Duc and Kawanishi [43], Phung-Duc [44],
Sanga and Charan [49], and Saravanan et al. [51], among others. We shall
demonstrate here that our approach (Section 2) does not need to exploit
further properties of the retrial queue under consideration, except for the
specific expressions of sub-matrices $Q_{i,i^{\prime}}$, for
$i\in\mathbb{N}_{0}$ and $i^{\prime}\in\\{\max\\{0,i-1\\},i,i+1\\}$, in Eq.
(6) defining the model. This is also applicable to queueing models with
impatient customers and flexible matching mechanism (Chai et al. [13]). In
epidemics (Baumann and Sandmann [7]; Lefèvre and Simon [37]), the first visit
of ${\cal X}$ to level $l(0)$ usually amounts to the duration of an outbreak,
and the maximum level of ${\cal X}$ and the state of $\\{J(t):t\geq 0\\}$ at
the first access to $l(0)$ are related to the maximum number of simultaneously
infected hosts during the outbreak and the state of the population at the peak
of infection, respectively. As in the case of retrial queues, our approach
(Section 2) should be seen to be a unified procedure in the analysis of
epidemic models at the first time of maximum incidence of the pathogen,
including multi-strain models (Chalub et al. [14]), stochastic chemical
kinetics (Dayar et al. [15]), gene family evolution (Diao et al. [18]),
resource-limited environments (Gómez-Corral et al. [26]), and recurrence in
tumors (Santana et al. [50]).

As a second motivation, let us recall that extreme values of the LD-QBD
process ${\cal X}$ in (6) are analyzed by Mandjes and Taylor [39] (see also
Javier and Fralix [30]) in terms of the _running maximum level_
$\bar{I}(T)=\sup\\{I(t):t\in[0,T]\\}$ attained by ${\cal X}$ at an independent
exponentially distributed time $T$ and, via Erlangization, at a predetermined
time $T$; for a related work in retrial queues, see Gómez-Corral and López
García [21]. We shall complement here the work of Mandjes and Taylor [39] by
characterizing the running maximum level of a LD-QBD process at time $T$ that
amounts to the first time process ${\cal X}$ makes a transition to states
contained in $l(0)$, and a detailed description of the phase process
$\\{J(t):t\geq 0\\}$ as the running maximum is first attained. Interested
readers in seeing a simpler description in epidemic models are referred to
Amador et al. [2].

This paper proceeds as follows. In Section 2, we first introduce the basic
notation for first-passage times of process ${\cal X}$, and give the formal
definition of extreme value and phase as the maximum level is first attained.
The distributional properties of the first time to reach the maximum level are
analyzed in terms of taboo Laplace-Stieltjes transforms of certain suitably
defined first-passage times (Section 2.2), and the probability law of the
maximum level visited by process ${\cal X}$ (Section 2.3). In Section 3, the
computational approach is exemplified with two epidemic models: the first one
is used to illustrate the solution in the case of a LD-QBD process with
infinitely many states (Section 3.1), and the second one to do so in the case
of a finite QBD process (Section 3.2). Section 4 concludes the paper with a
brief discussion.

##  2 Extreme values of a LD-QBD process and related hitting probabilities

Throughout this section, the LD-QBD process ${\cal X}$ is assumed to be
regular and we study the dynamics of ${\cal X}$ before the first visit to
level $l(0)$, under the assumption that the first passage from any initial
state $(i_{0},j_{0})\in{\cal S}\setminus l(0)$ to level $l(0)$ occurs in a
finite time with certainty. This assumption amounts to the fact that the null
column vector is the largest non-negative solution bounded above by 1 to
$Q^{*}v=0$, where $Q^{*}$ is the matrix obtained by removing from $Q$ in (6)
the rows and columns corresponding to states in $l(0)$; see Kulkarni [34,
Theorem 6.13].

For ease of presentation, we adopt the notation and terminology used in our
recent work [22, 23, 24] and related papers [10, 30, 39]. More particularly,
we define the first-passage time to level $l(i)$ as $\tau_{l(i)}=\inf\\{t\geq
0:I(t)=i\\}$, for integers $i\in\mathbb{N}_{0}$, and we let $I_{\max}$ denote
the running maximum level attained by ${\cal X}$ before the first entrance
into any state in level $l(0)$, and $\tau_{\max}$ be the first time to reach
the value $I_{\max}$. As a result, $J(\tau_{\max})$ is the phase variable at
time $\tau_{\max}$, and $I_{\max}=I(\tau_{\max})$.

Remark 1. Readers should note that Javier and Fralix [30] analyze the joint
probability law of both the state $(I(t),J(t))$ of process ${\cal X}$ and its
associated running maximum level $\bar{I}(t)$ at any arbitrary time $t$,
whereas the state $(I_{\max},J(\tau_{\max}))$ is inherently linked to the
random time $\tau_{\max}$, provided that the entrance of ${\cal X}$ to any
state in level $l(0)$ is certain.

For later use, we introduce some notation as follows: $0_{a}$ and $1_{a}$ are
the column vectors of order $a$ such that all their entries are equal to $0$
and $1$, respectively; $e_{a}(b)$ denotes the column vector of order $a$ with
its $b$th entry equals $1$, and $0$ elsewhere; $0_{a\times b}$ and $I_{a}$ are
the null matrix of dimension $a\times b$ and the identity matrix of order $a$,
respectively; T denotes transposition; and $\delta_{a,b}$ denotes the
Kronecker’s delta.

###  2.1 Statement of the problem

Our objective is to determine the probability law of $\tau_{\max}$ on those
sample paths of process ${\cal X}$ satisfying
$\\{I_{\max}=i,J(\tau_{\max})=j\\}$, for any state
$(i,j)\in\cup_{k=i_{0}}^{\infty}l(k)$ and a predetermined initial state
$(i_{0},j_{0})\in{\cal S}\setminus l(0)$. To start with, we notice that

$\displaystyle
P\left(\left.\tau_{\max}=0,I_{\max}=i,J(\tau_{\max})=j\right|(I(0),J(0))=(i_{0},j_{0})\right)$
$\displaystyle=\ \delta_{(i,j),(i_{0},j_{0})}P_{(i_{0},j_{0})}(i_{0}),$

where
$P_{(i_{0},j_{0})}(i)=P\left(\left.I_{\max}=i\right|(I(0),J(0))=(i_{0},j_{0})\right)$,
and then derive iterative schemes for evaluating the _restricted_ Laplace-
Stieltjes transform

$\displaystyle\psi_{(i_{0},j_{0})}(\theta;i,j)$ $\displaystyle=$
$\displaystyle
E\left[\left.e^{-\theta\tau_{\max}}1\\{I_{\max}=i,J(\tau_{\max})=j\\}\right|(I(0),J(0))=(i_{0},j_{0})\right],$

for $\Re(\theta)\geq 0$ and states $(i,j)\in\cup_{k=i_{0}+1}^{\infty}l(k)$.

In terms of the _taboo_ Laplace-Stieltjes transform

$\displaystyle\varphi_{(i_{0},j_{0})}(\theta;i,j)$ $\displaystyle=$
$\displaystyle
E\left[\left.e^{-\theta\tau_{l(i)}}1\\{\tau_{l(i)}<\tau_{l(0)},J(\tau_{l(i)})=j\\}\right|(I(0),J(0))=(i_{0},j_{0})\right],$

for $\Re(\theta)\geq 0$, the restricted Laplace-Stieltjes transform of
$\tau_{\max}$ on $\\{\tau_{\max}>0,I_{\max}=i,J(\tau_{\max})=j\\}$ is
determined by

$\displaystyle\psi_{(i_{0},j_{0})}(\theta;i,j)$ $\displaystyle=$
$\displaystyle\varphi_{(i_{0},j_{0})}(\theta;i,j)P_{(i,j)}(i),$ (7)

for $\Re(\theta)\geq 0$, states $(i,j)\in\cup_{k=i_{0}+1}^{\infty}l(k)$ and
the initial state $(i_{0},j_{0})\in{\cal S}\setminus l(0)$. To prove Eq. (7),
we observe that $\psi_{(i_{0},j_{0})}(\theta;i,j)$ is linked to the time
interval $[0,\tau_{\max})$ and the behavior of ${\cal X}$ during the residual
interval $[\tau_{\max},\tau_{l(0)})$ during which the subset of states
$\cup_{k=i+1}^{\infty}l(k)$ is avoided. More concretely, the random time
$\tau_{\max}$ on $\\{I_{\max}=i,J(\tau_{\max})=j\\}$ is identically
distributed to the first-passage time $\tau_{l(i)}$ on
$\\{J(\tau_{l(i)})=j\\}$ under the taboo of states in $l(0)$, which yields the
first term $\varphi_{(i_{0},j_{0})}(\theta;i,j)$ in the right-hand side of
(7). The second term $P_{(i,j)}(i)$ in the right-hand side of (7) guarantees
that, provided that $(I(\tau_{\max}),J(\tau_{\max}))=(i,j)$, the maximum level
that ${\cal X}$ visits during the residual interval
$[\tau_{\max},\tau_{l(0)})$ is not greater than the level variable at time
$\tau_{\max}$. Then, the proof of (7) follows by applying the strong Markov
property to the stopping time $\tau_{(i,j)}$ —defined as the time of visit to
state $(i,j)$— and the natural filtration.

###  2.2 Taboo Laplace-Stieltjes transforms of $\tau_{l(i)}$

In deriving expressions for $\varphi_{(i_{0},j_{0})}(\theta;i,j)$, for states
$(i,j)\in\cup_{k=i_{0}+1}^{\infty}l(k)$, we introduce taboo Laplace-Stieltjes
transforms $\varphi_{(i^{\prime},j^{\prime})}(\theta;i,j)$ in a more general
setting by replacing the initial state $(i_{0},j_{0})\in{\cal S}\setminus
l(0)$ by the _current_ state $(i^{\prime},j^{\prime})$ of process ${\cal X}$
at any arbitrary time, for $(i^{\prime},j^{\prime})\in\cup_{k=1}^{i-1}l(k)$.

Theorem 1. For any fixed state $(i,j)\in\cup_{k=i_{0}+1}^{\infty}l(k)$, the
column vectors
$\Phi_{i^{\prime}}(\theta;i,j)=(\varphi_{(i^{\prime},j^{\prime})}(\theta;i,j):j^{\prime}\in\\{0,...,M_{i^{\prime}}\\})$,
for integers $i^{\prime}\in\\{1,...,i-1\\}$, are uniquely characterized as the
solution of the system of linear equations

$\displaystyle\left(\theta
I_{1+M_{i^{\prime}}}-Q_{i^{\prime},i^{\prime}}\right)\Phi_{i^{\prime}}(\theta;i,j)$
$\displaystyle=$
$\displaystyle(1-\delta_{i^{\prime},1})Q_{i^{\prime},i^{\prime}-1}\Phi_{i^{\prime}-1}(\theta;i,j)$
(8)
$\displaystyle+Q_{i^{\prime},i^{\prime}+1}\Phi_{i^{\prime}+1}(\theta;i,j),\quad
i^{\prime}\in\\{1,...,i-2\\},~{}~{}~{}$

with

$\displaystyle\left(\theta
I_{1+M_{i-1}}-Q_{i-1,i-1}\right)\Phi_{i-1}(\theta;i,j)$ $\displaystyle=$
$\displaystyle Q_{i-1,i-2}\Phi_{i-2}(\theta;i,j)$ (9)
$\displaystyle+Q_{i-1,i}e_{1+M_{i}}(1+j).$

Moreover, the column vectors $m_{i^{\prime}}^{(n)}(i,j)$ consisting of the nth
moments of $\tau_{l(i)}$ on $\\{\tau_{l(i)}<\tau_{l(0)},J(\tau_{l(i)})=j\\}$,
for $n\in\mathbb{N}$ and integers $i^{\prime}\in\\{1,...,i-1\\}$, satisfy

$\displaystyle-Q_{i^{\prime},i^{\prime}}m_{i^{\prime}}^{(n)}(i,j)$
$\displaystyle=$ $\displaystyle
nm_{i^{\prime}}^{(n-1)}(i,j)+(1-\delta_{i,1})Q_{i^{\prime},i^{\prime}-1}m_{i^{\prime}-1}^{(n)}(i,j)$
(10)
$\displaystyle+Q_{i^{\prime},i^{\prime}+1}m_{i^{\prime}+1}^{(n)}(i,j),\quad
i^{\prime}\in\\{1,...,i-2\\},$ $\displaystyle-Q_{i-1,i-1}m_{i-1}^{(n)}(i,j)$
$\displaystyle=$ $\displaystyle
nm_{i-1}^{(n-1)}(i,j)+Q_{i-1,i-2}m_{i-2}^{(n)}(i,j),$ (11)

with
$m_{i^{\prime}}^{(0)}(i,j)=\left.\Phi_{i^{\prime}}(\theta;i,j)\right|_{\theta=0}$.

Proof. Let us fix the state $(i,j)\in\cup_{k=i_{0}+1}^{\infty}l(k)$. By
conditioning on the first jump of process ${\cal X}$ when
$(i^{\prime},j^{\prime})$ is its current state, it is found that

(i)

For states $(i^{\prime},j^{\prime})\in\cup_{k=1}^{i-2}l(k)$,

$\displaystyle\varphi_{(i^{\prime},j^{\prime})}(\theta;i,j)$ $\displaystyle=$
$\displaystyle(1-\delta_{i^{\prime},1})\sum_{j^{\prime\prime}=0}^{M_{i^{\prime}-1}}\frac{q_{(i^{\prime},j^{\prime}),(i^{\prime}-1,j^{\prime\prime})}}{\theta+q_{(i^{\prime},j^{\prime})}}\varphi_{(i^{\prime}-1,j^{\prime\prime})}(\theta;i,j)$
(12) $\displaystyle+\sum_{j^{\prime\prime}=0,j^{\prime\prime}\neq
j^{\prime}}^{M_{i^{\prime}}}\frac{q_{(i^{\prime},j^{\prime}),(i^{\prime},j^{\prime\prime})}}{\theta+q_{(i^{\prime},j^{\prime})}}\varphi_{(i^{\prime},j^{\prime\prime})}(\theta;i,j)$
$\displaystyle+\sum_{j^{\prime\prime}=0}^{M_{i^{\prime}+1}}\frac{q_{(i^{\prime},j^{\prime}),(i^{\prime}+1,j^{\prime\prime})}}{\theta+q_{(i^{\prime},j^{\prime})}}\varphi_{(i^{\prime}+1,j^{\prime\prime})}(\theta;i,j),$

where
$q_{(i^{\prime},j^{\prime})}=-q_{(i^{\prime},j^{\prime}),(i^{\prime},j^{\prime})}$.

(ii)

For integers $j^{\prime}\in\\{0,...,M_{i-1}\\}$,

$\displaystyle\varphi_{(i-1,j^{\prime})}(\theta;i,j)$ $\displaystyle=$
$\displaystyle\sum_{j^{\prime\prime}=0}^{M_{i-2}}\frac{q_{(i-1,j^{\prime}),(i-2,j^{\prime\prime})}}{\theta+q_{(i-1,j^{\prime})}}\varphi_{(i-2,j^{\prime\prime})}(\theta;i,j)$
(13) $\displaystyle+\sum_{j^{\prime\prime}=0,j^{\prime\prime}\neq
j^{\prime}}^{M_{i-1}}\frac{q_{(i-1,j^{\prime}),(i-1,j^{\prime\prime})}}{\theta+q_{(i-1,j^{\prime})}}\varphi_{(i-1,j^{\prime\prime})}(\theta;i,j)$
$\displaystyle+\frac{q_{(i-1,j^{\prime}),(i,j)}}{\theta+q_{(i-1,j^{\prime})}}.$

Straightforward algebra allows us to verify that, by multiplying (12)-(13) by
$\theta+q_{(i^{\prime},j^{\prime})}$ and $\theta+q_{(i-1,j^{\prime})}$,
respectively, the matrix versions of the resulting equations are (8)-(9).
Then, Eqs. (10)-(11) are derived by taking derivatives in (8)-(9) at
$\theta=0$, and observing that the column vector
$m_{i^{\prime}}^{(n)}(i,j)=(-1)^{n}\left.d^{n}\Phi_{i^{\prime}}(\theta;i,j)/d\theta^{n}\right|_{\theta=0}$
contains the taboo expectations

$\displaystyle
E\left[\left.\tau_{l(i)}^{n}1\\{\tau_{l(i)}<\tau_{l(0)},J(\tau_{l(i)})=j\\}\right|(I(0),J(0))=(i^{\prime},j^{\prime})\right],$

for any current state $(i^{\prime},j^{\prime})\in\cup_{k=1}^{i-1}l(k)$. This
completes the proof. $\square$

An important observation is that Eqs. (8)-(9) (respectively, (10)-(11)) become
a structured system of linear equations in the unknowns
$\Phi_{i^{\prime}}(\theta;i,j)$ (respectively, $m_{i^{\prime}}^{(n)}(i,j)$)
with $i^{\prime}\in\\{1,...,i-1\\}$, whence its solution can be evaluated by
applying any general-purpose procedure, such as block Gaussian elimination;
see, e.g., Stewart [52]. This is summarized in Algorithms A-B, where the
column vectors $m_{i^{\prime}}^{(0)}(i,j)$ with $i^{\prime}\in\\{1,...,i-1\\}$
record the values
$\left.\varphi_{(i^{\prime},j^{\prime})}(\theta;i,j)\right|_{\theta=0}$, for
integers $j^{\prime}\in\\{0,...,M_{i^{\prime}}\\}$, and matrices
$H_{i^{\prime}}(i,j)$ with $i^{\prime}\in\\{1,...,i-1\\}$ are defined by
$\left.H_{i^{\prime}}(\theta;i,j)\right|_{\theta=0}$; as a result,
$m_{i^{\prime}}^{(0)}(i,j)$ and $H_{i^{\prime}}(i,j)$ in Algorithm B are
derived from Algorithm A by selecting $\theta=0$.

Algorithm A: Computation of the column vectors $\Phi_{i^{\prime}}(\theta;i,j)$
with $i^{\prime}\in\\{1,...,i-1\\}$, for $\Re(\theta)\geq 0$ and a fixed state
$(i,j)\in\cup_{k=i_{0}+1}^{\infty}l(k)$.

Step 1: $i^{\prime}:=1$;

$H_{i^{\prime}}(\theta;i,j):=(\theta
I_{1+M_{i^{\prime}}}-Q_{i^{\prime},i^{\prime}})^{-1}Q_{i^{\prime},i^{\prime}+1}$;

while $i^{\prime}<i-2$, repeat

$i^{\prime}:=i^{\prime}+1$;

$H_{i^{\prime}}(\theta;i,j):=(\theta
I_{1+M_{i^{\prime}}}-Q_{i^{\prime},i^{\prime}}-Q_{i^{\prime},i^{\prime}-1}H_{i^{\prime}-1}(\theta;i,j))^{-1}Q_{i^{\prime},i^{\prime}+1}$;

$H_{i^{\prime}}(\theta;i,j):=H_{i^{\prime}}(\theta;i,j)e_{1+M_{i^{\prime}+1}}(1+j)$.

Step 2: $\Phi_{i^{\prime}}(\theta;i,j):=H_{i^{\prime}}(\theta;i,j)$;

while $i^{\prime}>1$, repeat

$i^{\prime}:=i^{\prime}-1$;

$\Phi_{i^{\prime}}(\theta;i,j):=H_{i^{\prime}}(\theta;i,j)\Phi_{i^{\prime}+1}(\theta;i,j)$.

Algorithm B: Computation of the column vectors $m_{i^{\prime}}^{(n)}(i,j)$
with $i^{\prime}\in\\{1,...,i-1\\}$, for $n\in\mathbb{N}$ and a fixed state
$(i,j)\in\cup_{k=i_{0}+1}^{\infty}l(k)$, in terms of
$m_{i^{\prime}}^{(n-1)}(i,j)$ with $i^{\prime}\in\\{1,...,i-1\\}$.

Step 1: $i^{\prime}:=1$;

$h_{i^{\prime}}(i,j):=-nQ_{i^{\prime},i^{\prime}}^{-1}m_{i^{\prime}}^{(n-1)}(i,j)$;

while $i^{\prime}<i-1$, repeat

$i^{\prime}:=i^{\prime}+1$;

$h_{i^{\prime}}(i,j):=(-Q_{i^{\prime},i^{\prime}}-Q_{i^{\prime},i^{\prime}-1}H_{i^{\prime}-1}(i,j))^{-1}$;

$h_{i^{\prime}}(i,j):=h_{i^{\prime}}(i,j)(Q_{i^{\prime},i^{\prime}-1}h_{i^{\prime}-1}(i,j)+nm_{i^{\prime}}^{(n-1)}(i,j))$.

Step 2: $m_{i^{\prime}}^{(n)}(i,j):=h_{i^{\prime}}(i,j)$;

while $i^{\prime}>1$, repeat

$i^{\prime}:=i^{\prime}-1$;

$m_{i^{\prime}}^{(n)}(i,j):=H_{i^{\prime}}(i,j)m_{i^{\prime}+1}^{(n)}(i,j)+h_{i^{\prime}}(i,j)$.

From Algorithms A-B, the taboo Laplace-Stieltjes transform of $\tau_{l(i)}$ on
$\\{\tau_{l(i)}<\tau_{l(0)},J(\tau_{l(i)})=j\\}$ is easily obtained as

$\displaystyle\varphi_{(i_{0},j_{0})}(\theta;i,j)$ $\displaystyle=$
$\displaystyle e^{T}_{1+M_{i_{0}}}(1+j_{0})\Phi_{i_{0}}(\theta;i,j),$

and the expectation
$E\left[\left.(\tau_{\max})^{n}1\\{I_{\max}=i,J(\tau_{\max})=j\\}\right|(I(0),J(0))=(i_{0},j_{0})\right]$
is evaluated as

$\displaystyle e^{T}_{1+M_{i_{0}}}(1+j_{0})m_{i_{0}}^{(n)}(i,j)P_{(i,j)}(i),$

for any state $(i,j)\in\cup_{k=i_{0}+1}^{\infty}l(k)$.

###  2.3 Marginal distribution of $I_{\max}$

Using processes under a taboo (see, e.g., Section 2.1 in Gómez-Corral and
López García [20]), the mass function
$\\{P_{(i_{0},j_{0})}(i):i\in\mathbb{N}_{0}\\}$ of $I_{\max}$ is found to be
$P_{(i_{0},j_{0})}(i)=0$ if $i\in\\{0,...,i_{0}-1\\}$; otherwise,
$P_{(i_{0},j_{0})}(i)=F_{\max}(i;i_{0}j_{0})-(1-\delta_{i,i_{0}})F_{\max}(i-1;i_{0},j_{0})$,
where the probability distribution function
$F_{\max}(i;i_{0},j_{0})=P\left(\left.I_{\max}\leq
i\right|(I(0),J(0))=(i_{0},j_{0})\right)$ has the form

$\displaystyle F_{\max}(i;i_{0},j_{0})$ $\displaystyle=$ $\displaystyle
e^{T}_{c(i)}((1-\delta_{i_{0},1})c(i_{0}-1)+1+j_{0})\left(-T^{-1}(i)t_{0}(i)\right),$
(14)

for $i\in\\{i_{0},i_{0}+1,...\\}$. In this expression,
$c(i)=\sum_{i^{\prime}=1}^{i}(1+M_{i^{\prime}})$ is the cardinality of the
subset of states $\cup_{k=1}^{i}l(k)$, the matrix $T(i)$ consists of
infinitesimal rates $q_{(i,j),(i^{\prime},j^{\prime})}$, for states
$(i,j),(i^{\prime},j^{\prime})\in\cup_{k=1}^{i}l(k)$, and $t_{0}(i)$ is the
column vector

$\displaystyle t_{0}(i)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}Q_{1,0}1_{1+M_{0}}\\\
0_{c(i)-(1+M_{1})}\end{array}\right).$

It is worth noting that $-T^{-1}(i)$ in (14) is a matrix of expected total
times spent in states of $\cup_{k=1}^{i}l(k)$, whence $-T^{-1}(i)$ is a non-
negative matrix and $-T^{-1}(i)t_{0}(i)$ contains probabilities that the first
visit of ${\cal X}$ to $l(0)$ occurs in a finite time, provided that process
${\cal X}$ does not visit states in the subset $\cup_{k=i+1}^{\infty}l(k)$,
from any initial state in $\cup_{k=1}^{i}l(k)$.

Algorithm C derives the mass function of $I_{\max}$ in a recursive manner. Its
proof mostly repeats arguments of Ref. [20, Section 2.1], whence we omit the
details.

Algorithm C: Computation of the non-null probabilities $P_{(i_{0},j_{0})}(i)$,
for $i\in\\{i_{0},i_{0}+1,...\\}$ and a predetermined initial state
$(i_{0},j_{0})\in{\cal S}\setminus l(0)$.

Step 1: $i:=1$;

$-T^{-1}(i):=-Q_{i,i}^{-1}$;

if $i=i_{0}$, then

$F_{\max}(i;i_{o},j_{0}):=e^{T}_{c(i)}(1+j_{0})(-T^{-1}(i)t_{0}(i))$;

$P_{(i_{0},j_{0})}(i):=F_{\max}(i;i_{0},j_{0})$.

Step 2: While $F_{\max}(i;i_{0},j_{0})<1-\varepsilon$, repeat

$i:=i+1$;

$A_{1,2}:=\left(\begin{array}[]{c}0_{c(i-2)\times(1+M_{i})}\\\
Q_{i-1,i}\end{array}\right)$;

$A_{2,1}:=\left(0_{(1+M_{i})\times c(i-2)},Q_{i,i-1}\right)$;

$B_{2,2}:=\left(-Q_{i,i}-A_{2,1}(-T^{-1}(i-1))A_{1,2}\right)^{-1}$;

$B_{2,1}:=B_{2,2}A_{2,1}(-T^{-1}(i-1))$;

$B_{1,2}:=-T^{-1}(i-1)A_{1,2}B_{2,2}$;

$B_{1,1}:=-T^{-1}(i-1)\left(I_{c(i-1)}+A_{1,2}B_{2,1}\right)$;

$-T^{-1}(i):=\left(\begin{array}[]{cc}B_{1,1}&B_{1,2}\\\
B_{2,1}&B_{2,2}\end{array}\right)$;

if $i\geq i_{0}$, then

$F_{\max}(i;i_{0},j_{0}):=e^{T}_{c(i)}(c(i_{0}-1)+1+j_{0})\left(-T^{-1}(i)t_{0}(i)\right)$;

$P_{(i_{0},j_{0})}(i):=F_{\max}(i;i_{0},j_{0})-F_{\max}(i-1;i_{0},j_{0})$.

Note that, from an appropriate use of Algorithms A and C, it is possible to
evaluate the following expressions for the hitting probabilities at time
$\tau_{\max}$ and the $n$th moment of $J(\tau_{\max})$ on $\\{I_{\max}=i\\}$:

(i)

For states $(i,j)\in\cup_{k=i_{0}+1}^{\infty}l(k)$,

$\displaystyle
P\left(\left.I_{\max}=i,J(\tau_{\max})=j\right|(I(0),J(0))=(i_{0},j_{0})\right)$
$\displaystyle=\
\left.\varphi_{(i_{0},j_{0})}(\theta;i,j)\right|_{\theta=0}P_{(i,j)}(i).$

(ii)

For $n\in\mathbb{N}_{0}$ and integers $i\in\\{i_{0}+1,i_{0}+2,...\\}$,

$\displaystyle
E\left[\left.(J(\tau_{\max}))^{n}1\\{I_{\max}=i\\}\right|(I(0),J(0))=(i_{0},j_{0})\right]$
$\displaystyle=\
\sum_{j=0}^{M_{i}}j^{n}\left.\varphi_{(i_{0},j_{0})}(\theta;i,j)\right|_{\theta=0}P_{(i,j)}(i).$

To conclude the section, we briefly present the solution for a finite LD-QBD
process.

Remark 2. For a process ${\cal X}$ on the finite space state ${\cal
S}=\cup_{i=0}^{N}l(i)$, Algorithms A-B should be executed progressively for
any state $(i,j)\in\cup_{k=i_{0}+1}^{N}l(k)$ and a fixed initial state
$(i_{0},j_{0})\in\cup_{k=1}^{N-1}l(k)$. In Algorithm C, the probabilities
$P_{(i_{0},j_{0})}(i)$ can be evaluated for integers
$i\in\\{i_{0},i_{0}+1,...,N\\}$ by modifying Step 2 with the condition $i<N-1$
and noting that $F_{\max}(N;i_{0},i_{0})=1$.

##  3 Epidemic models

In this section, the analytical solution in Section 2 is exemplified with two
epidemic models: the SIS model for a pathogen transmitted both vertically and
horizontally; and the SIR model with constant population size. For
convenience, the number $I(t)$ is related to infectious individuals in
Subsections 3.1-3.2, and the level of process ${\cal X}$ is suitably defined
according to the extinction of the population (Subsection 3.1) and the
duration of an outbreak (Subsection 3.2).

###  3.1 The SIS epidemic model for vertically and horizontally transmitted
pathogens

In the terminology of Reuter [46], the SIS model with vertical and horizontal
transmission is defined in Ref. [25] as a competition process ${\cal
Y}=\\{(I(t),S(t)):t\geq 0\\}$, where $I(t)$ and $S(t)$ are related to the
number of infectious and susceptible individuals at time $t$, respectively.
For states
$(i,s),(i^{\prime},s^{\prime})\in\mathbb{N}_{0}\times\mathbb{N}_{0}$ with
$(i^{\prime},s^{\prime})\neq(i,s)$, the dynamics of ${\cal Y}$ are specified
by the infinitesimal rates

$\displaystyle q_{(i,s),(i^{\prime},s^{\prime})}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\beta is,&\hbox{if $i^{\prime}=i+1$
and $s^{\prime}=s-1$,}\\\ (1-p)\beta_{I}i,&\hbox{if $i^{\prime}=i+1$ and
$s^{\prime}=s$,}\\\ \delta_{S}s,&\hbox{if $i^{\prime}=i$ and
$s^{\prime}=s-1$,}\\\ p\beta_{I}i+\beta_{S}s,&\hbox{if $i^{\prime}=i$ and
$s^{\prime}=s+1$,}\\\ \delta_{I}i,&\hbox{if $i^{\prime}=i-1$ and
$s^{\prime}=s$,}\\\ \gamma i,&\hbox{if $i^{\prime}=i-1$ and
$s^{\prime}=s+1$,}\end{array}\right.$ (22)

and $q_{(i,s),(i,s)}=-(\beta
is+(\delta_{I}+\beta_{I}+\gamma)i+(\delta_{S}+\beta_{S})s)$, for
$(i,j)\in\mathbb{N}_{0}\times\mathbb{N}_{0}$, where $\beta$, $\beta_{S}$,
$\beta_{I}$, $\delta_{S}$, $\delta_{I}$ and $\gamma$ are strictly positive,
and $p\in(0,1)$. In particular, the value $\beta$ in (22) is the contact rate;
parameters $\beta_{S}$ and $\delta_{S}$ (respectively, $\beta_{I}$ and
$\delta_{I}$) are the birth and death rates of susceptible (respectively,
infectious) individuals; $\gamma$ is the recovery rate of an infectious
individual; and $p$ represents the proportion of the offspring of infectious
parents that is born as susceptible.

In Figures 1-3, we focus on the life cycle of a population initially
consisting of a single infectious individual and twenty susceptible ones,
i.e., $I(0)=1$ and $S(0)=20$. Thus, the analytical results in Section 2 have
to be applied here to process ${\cal X}=\\{(N(t),I(t)):t\geq 0\\}$ with
$N(t)=I(t)+S(t)$, instead of ${\cal Y}$. Process ${\cal X}$ can be routinely
formulated as a LD-QBD process111Specifically, we let $N(t)$ and $I(t)$ be the
level and the phase variables of ${\cal X}$, respectively. Hence,
$N_{\max}=N(\tau_{\max})$ is the maximum population size before extinction.
defined on ${\cal S}=\cup_{n=0}^{\infty}l(n)$ with levels
$l(n)=\\{(n,i):i\in\\{0,...,n\\}\\}$, for $n\in\mathbb{N}_{0}$, in such a way
that the time taken to access the absorbing state $(0,0)$ (i.e., level $l(0)$)
corresponds to the extinction time. This allows quantifying the incidence of
the disease at the time of maximum population size through the number
$I(\tau_{\max})$ of infectious individuals, which acts as the phase of ${\cal
X}$ at time $\tau_{\max}$.

Figure 1: Conditional probabilities
$P\left(\left.N_{\max}=n,I(\tau_{\max})\leq
i\right|(N(0),I(0))=(21,1)\right)$, for integers $n\in\\{21,...,n_{0.95}\\}$
and $i\in\\{\delta_{n,21},...,n\\}$, in scenarios with rates
$\beta\in\\{1.0,2.5,5.0\\}$ and $\beta_{I}\in\\{0.625,1.25,2.5\\}$.

Similarly to Ref. [25, Section 4], birth and death rates in our experiments
are assumed to be related to each other through the equalities
$\delta_{S}=\beta_{S}$, $\beta_{S}=2\beta_{I}$ and $\delta_{I}=2\delta_{S}$.
Hence, the compartment of susceptible individuals is expected to be
demographically stable, offspring from infectious parents are expected to be
born less frequently than those from susceptible parents, and susceptible
individuals have a longer life expectancy than infectious ones. With the
selections $p=0.2$ and $\gamma^{-1}=1.0$, it is supposed that $20\%$ of
newborn offspring of infectious parents will be susceptible and that an
infectious individual will be recovered, on average, in $1$ unit of time.
Contact rates $\beta\in\\{1.0,2.5,5.0\\}$ and birth rates
$\beta_{I}\in\\{0.625,1.25,2.5\\}$ yield nine scenarios in Figures 1-3, for
which the extinction of the population is certain and occurs in a finite mean
time222Since process ${\cal X}$ is regular, the inequality
$\min\\{\beta_{S},\beta_{I}\\}<\delta_{S}+\delta_{I}$ can be seen as a
sufficient condition for ${\cal X}$ to be absorbed almost surely, provided
that $\min\\{\delta_{S},\delta_{I}\\}=\min\\{\beta_{S},\beta_{I}\\}$; see Ref.
[24, Theorems 1 and 3]. Indeed, the absorption of ${\cal X}$ occurs in a
finite mean time from any initial non-absorbing state.. For convenience,
intervals on the axes _ox_ and _oy_ in Figure 1, and on the _ox_ axis in
Figures 2-3 are chosen to concentrate at least $0.95$ of the probability,
whence the integer $n_{0.95}$ in these figures satisfies
$F_{\max}(n_{0.95}-1;21,1)<0.95\leq F_{\max}(n_{0.95};21,1)$.

Figure 2: The mass function
$\\{P\left(\left.N_{\max}=n\right|(N(0),I(0))=(21,1)\right):n\in\\{21,...,n_{0.95}\\}\\}$
in scenarios with rates $\beta\in\\{1.0,2.5,5.0\\}$ and
$\beta_{I}\in\\{0.625,1.25,2.5\\}$. Figure 3: Conditional probabilities
$\\{P\left(\left.I(\tau_{\max})=i\right|(N(0),I(0))=(21,1),N_{\max}=n\right):i\in\\{0,...,n\\}\\}$,
for selected numbers $n\in\\{22,24,26\\}$, in scenarios with rates
$\beta\in\\{1.0,2.5,5.0\\}$ and $\beta_{I}\in\\{0.625,1.25,2.5\\}$.

It is seen that, in our experiments, the transmission dynamics of the pathogen
are limited and the duration of the outbreak is expected to be short. This is
mainly related to the values of $P\left(\left.N_{\max}=n,I(\tau_{\max})\leq
i\right|(N(0),I(0))=(21,1)\right)$ in Figure 1, especially because the most
significant ones correspond to the integers $n=21$ and $i=1$ (i.e., the
initial population size and the initially infectious individual), regardless
of the choice of $(\beta,\beta_{I})$. In Figure 2, it is seen that the mass
function of $N_{\max}$ is unimodal and is mainly concentrated on the initial
population size and the closest values to it. Note that, for a fixed value of
$\beta_{I}$, shorter life cycles are expected to occur when isolation measures
between compartments are relaxed (i.e., increasing values of $\beta$), which
is explained by seeing that infectious individuals die at a higher rate than
susceptible ones. Although a higher birth rate from infectious parents leads
to an increased number of infectious newborns, it will also generate a greater
number of susceptible newborns, who are born in a lower proportion than
infectious ones but with longer life spans and a higher birth rate. Therefore,
heavier-tailed distributions of $N_{\max}$ —which would predict longer
durations of a life cycle— are observed when isolation measures between
compartments become more significant and, in particular, the birth rate from
infectious parents increases.

A detailed description of the compartment of infectious individuals at the
time of maximum population size is displayed in Figure 3 in terms of the
conditional probabilities
$P\left(\left.I(\tau_{\max})=i\right|(N(0),I(0))=(21,1),N_{\max}=n\right)$,
for integers $i\in\\{0,...,n\\}$ and selections $n\in\\{22,24,26\\}$. The
unimodal conditional distribution of $I(\tau_{\max})$, provided that
$N_{\max}=n$, reflects how a more accurate estimation of the maximum
population size $N_{\max}$ will allow a better prediction of the incidence of
the disease in terms of $I(\tau_{\max})$, as observed, for example, in the
case $n=22$ and any of our scenarios. For some values of $(\beta,\beta_{I})$,
the selection $n=26$ leads to a bimodal conditional distribution of
$I(\tau_{\max})$, with a first peak at a small integer $i$ and a second one at
$i=n$. This behavior is linked to the inherent stochasticity of the Markov
chain model, for which sample paths of process ${\cal X}$ on
$\\{N_{\max}=n\\}$ in the case $n=26$ and some pairs $(\beta,\beta_{I})$ may
make it significantly plausible that the outbreak has a short duration in time
and a minimal incidence in terms of the number $I(\tau_{\max})$ of infectious
individuals, and at the same time a major outbreak with values of
$I(\tau_{\max})$ comparable to the maximum population size can also be
observed.

###  3.2 The SIR epidemic model

In the SIR model (Kermack and McKendrick [32]; Neuts and Li [41]), the process
${\cal X}=\\{(I(t),S(t)):t\geq 0\\}$ records the number $I(t)$ of infectious
individuals and the number $S(t)$ of susceptible individuals, whereas the
number $R(t)$ of removed or recovered individuals is given by
$R(t)=N-I(t)-S(t)$, where $N$ is the (constant) population size. The process
${\cal X}$ is then defined on the state space

$\displaystyle{\cal S}$ $\displaystyle=$
$\displaystyle\\{(i,s):i\in\\{0,...,N_{I}+N_{S}-s\\},s\in\\{0,...,N_{S}\\}\\},$

for initial numbers of $N_{I}$ infectious individuals and of $N_{S}$
susceptible individuals —consequently, $N=N_{I}+N_{S}$ with $R(0)=0$—, and its
$q$-matrix $Q$ has entries

$\displaystyle q_{(i,s),(i^{\prime},s^{\prime})}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}N^{-1}\beta is,&\hbox{if
$i^{\prime}=i+1$ and $s^{\prime}=s-1$,}\\\ -(N^{-1}\beta is+\gamma
i),&\hbox{if $i^{\prime}=i$ and $s^{\prime}=s$,}\\\ \gamma i,&\hbox{if
$i^{\prime}=i-1$ and $s^{\prime}=s$,}\\\
0,&\hbox{otherwise,}\end{array}\right.$

for states $(i,s),(i^{\prime},s^{\prime})\in{\cal S}$, where the parameters
$\beta$ and $\gamma$ are the contact and recovery rates, respectively, and are
assumed to be strictly positive.

Figure 4:  The mass function
$\\{P\left(\left.I_{\max}=i\right|(I(0),S(0))=(1,24)\right):i\in\\{1,...,25\\}\\}$
for populations with value ${\cal R}_{0}=1.5$ (left) and $3.8$ (right).
Figure 5: Conditional probabilities
$\\{P\left(\left.R(\tau_{\max})=r\right|(I(0),S(0))=(1,24),I_{\max}=i\right):r\in\\{0,...,25-i\\}\\}$,
for selected numbers $i\in\\{3,8,13,18,23\\}$, and populations with value
${\cal R}_{0}=1.5$ (left) and $3.8$ (right).

In Figures 4-5 and Table 1, we consider a small community consisting initially
of one infectious and twenty-four susceptible individuals, i.e., $I(0)=1$ and
$S(0)=24$. An infectious individual is assumed to recover, on average, in
$\gamma^{-1}=1$ unit of time, and the contact rate $\beta$ is determined from
the basic reproductive number ${\cal R}_{0}=\gamma^{-1}\beta$. The choice of
${\cal R}_{0}$ values corresponds to the estimated mean values in Ref. [54,
Table 1] for the 1918 Spanish influenza. Figure 4 illustrates the transmission
dynamics in terms of the mass function333Since $I(t)$ and $S(t)$ specify the
level and the phase variables of ${\cal X}$, respectively,
$I_{\max}=I(\tau_{\max})$ represents the maximum size of the compartment of
infectious individuals during an outbreak.
$\\{P\left(\left.I_{\max}=i\right|(I(0),S(0))=(1,24)\right):i\in\\{1,...,25\\}\\}$,
which is plotted for values ${\cal R}_{0}=1.5$ (spring wave) and $3.8$ (fall
wave), and Table 1 gives the expected values
$E\left[\left.\tau_{\max}\right|(I(0),S(0))=(1,24),I_{\max}=i\right]$, for
selected numbers $i\in\\{3,8,13,18,23\\}$. We notice that, not surprisingly,
the less interaction between individuals, the less spread of the pathogen will
be observed. In particular, the maximum number $I_{\max}$ of infectious
individuals in the spring wave has an unimodal distribution with a clear peak
at point $i=1$, meaning that the event $\\{I_{\max}=1\\}$ occurs more
frequently than the others. As a result, the disease is more likely to be
eradicated as soon as the initially infectious individual recovers than to
spread to the subpopulation of susceptible individuals. However, the latter
may also occur with significant probability. A more relevant interaction
between individuals in the fall wave yields a bimodal distribution of
$I_{\max}$, with a first peak at point $i=1$ and a second one at $i=13$. This
illustrates how increasing values of ${\cal R}_{0}$ may cause major outbreaks
—i.e., outbreaks with a significant number of infectious cases— at the same
time that the disease could also disappear at an early stage of the epidemic.
This fact highlights the need to analyze the maximum number $I_{\max}$ of
infectious individuals as a random variable (i.e., in terms of its mass
function) rather than its expected value or its deterministic counterpart.

Table 1: Expected values $E\left[\left.\tau_{\max}\right|(I(0),S(0))=(1,24),I_{\max}=i\right]$ for selected numbers $i\in\\{3,8,13,18,23\\}$, and populations with ${\cal R}_{0}=1.5$ and $3.8$. ${\cal R}_{0}$ | $i$ | $E\left[\left.\tau_{\max}\right|(I(0),S(0))=(1,24),I_{\max}=i\right]$
---|---|---
1.5 | 3 | 0.84726
| 8 | 2.44037
| 13 | 2.27654
| 18 | 2.07524
| 23 | 1.91487
3.8 | 3 | 0.41006
| 8 | 1.73395
| 13 | 1.45888
| 18 | 1.26555
| 23 | 1.13595

To further investigate how the values of $I_{\max}$ affect the subpopulation
of recovered individuals at the time $\tau_{\max}$ of greatest incidence of
the pathogen, a detailed description of $R(\tau_{\max})$ is displayed in
Figure 5, in terms of the mass function
$\\{P\left(\left.R(\tau_{\max})=r\right|(I(0),S(0))=(1,24),I_{\max}=i\right):r\in\\{0,...,25-i\\}\\}$,
for integers $i\in\\{3,8,13,18,23\\}$. It is important to note that the most
likely event $\\{I_{\max}=i\\}$ in Figure 4 is linked to the choice $i=1$,
irrespectively of ${\cal R}_{0}$, which leads to $R(\tau_{\max})=0$ almost
surely (since $\tau_{\max}=0$), provided that $(I(0),S(0))=(1,24)$ and
$I_{\max}=1$. The particular selection $i^{*}=13$ is related to the second
peak of the bimodal distribution of $I_{\max}$ in the fall wave. Therefore,
the conditional probabilities
$\\{P\left(\left.R(\tau_{\max})=r\right|(I(0),S(0))=(1,24),I_{\max}=i^{*}\right):r\in\\{0,...,25-i^{*}\\}\\}$
must be seen as an accurate description of the variability of $R(\tau_{\max})$
for a major outbreak in the case ${\cal R}_{0}=3.8$, whereas they should not
be usually appropriate when ${\cal R}_{0}=1.5$ as $i^{*}=13$ leads to an event
$\\{I_{\max}=i^{*}\\}$ that it will hardly ever occur. Our experiments always
found that the conditional distribution of $R(\tau_{\max})$ is unimodal,
regardless of the selected integer $i$.

##  4 Conclusions

For a regular LD-QBD process ${\cal X}=\\{(I(t),J(t)):t\geq 0\\}$ with
countably many states, the probability law of the random vector
$(\tau_{\max},I_{\max},J(\tau_{\max}))$ has been characterized under the
assumption that the first passage from any initial state
$(i_{0},j_{0})\in{\cal S}\setminus l(0)$ to level $l(0)$ is certain. For this
purpose, efficient algorithms for computing the marginal distribution of
$I_{\max}$, and restricted Laplace-Stieltjes transforms of $\tau_{\max}$ on
the sample paths of ${\cal X}$ such that $\\{I_{\max}=i,J(\tau_{\max})=j\\}$,
have been developed as an extension of the well-known block-Gaussian
elimination technique, whence the algorithmic complexity of the solution in
Subsections 2.2 and 2.3 is similar to that in Refs. [19, Algorithm A] and [22,
Algorithms 1A-3A] in the finite case.

The computational approach can be thought of as an alternative one to that of
Mandjes and Taylor [39], and Javier and Fralix [30] for the running maximum
$\bar{I}(T)$ attained by process ${\cal X}$ during a time interval $[0,T]$
when $T$ is either independent and exponentially distributed time, or a fixed
time. Here, time $T$ is described as the first time that process ${\cal X}$
visits states in level $l(0)$, so $T$ is closely linked to the dynamics of
${\cal X}$ and, for practical use, amounts to the duration of a life cycle in
a population, of an outbreak in epidemics, and of a busy period in a queueing
model. In the context of epidemics, an aspect that deserves further
exploration is how phase $J(\tau_{\max})$ and its probability law —when level
$I(t)$ in process ${\cal X}$ is defined as the number of infectious
individuals— could be used to characterize a stochastic version of the herd
immunity threshold (see, e.g., Britton et al. [12]), for which it is crucial
to notice that no more endemic transmissions are expected to occur after time
$\tau_{\max}$.

## Acknowledgements

This research was supported by Ministerio de Ciencia e Innovación (Project
PID2021-125871NB-B-I00), as well as by GNCS-INdAM and “European Union – Next
Generation EU” PRIN 2022 PNRR (Project P2022XSF5H). D. Taipe also acknowledges
the support of Banco Santander-Universidad Complutense de Madrid (Pre-doctoral
Researcher Contract CT63/19-CT64/19).

## Data availability

No data was used for the research described in the article. Details on
numerical experiments in Section 3 will be made available from the
corresponding author upon reasonable request.

## Authors’ contributions

All the authors equally contributed to this work.

## References

  * [1] Akar N, Oğuz NC, Sohraby K. A novel computational method for solving finite QBD processes. Commun Stat Stoch Model. 2000;16:273–311. doi.org/10.1080/15326340008807588
  * [2] Amador J, Armesto D, Gómez-Corral A. Extreme values in SIR epidemic models with two strains and cross-immunity. Math Biosci Eng. 2019;16:1992–2022. doi.org/10.3934/mbe.2019098
  * [3] Anderson WJ. Continuous-Time Markov Chains. An Applications-Oriented Approach. New York: Springer-Verlag; 1991.
  * [4] Artalejo JR, Gómez-Corral A. Retrial Queueing Systems. A Computational Approach. Berlin Heidelberg: Springer-Verlag; 2008.
  * [5] Baumann H, Sandmann W. Numerical solution of level dependent quasi-birth-and-death processes. Procedia Comput Sci. 2010;1(1):1561–1569. doi.org/10.1016/j.procs.2010.04.175
  * [6] Baumann H, Sandmann W. Computing stationary expectations in level-dependent QBD processes. J Appl Probab. 2013;50(1):151–165. doi.org/10.1239/jap/1363784430
  * [7] Baumann H, Sandmann W. Structured modeling and analysis of stochastic epidemics with immigration and demographic effects. PLoS ONE 2016;11(3):e0152144. doi.org/10.1371/journal.pone.0152144
  * [8] Baumann H. Finite-state-space truncations for infinite quasi-birth-death processes. J Appl Math. 2020;2020:2678374. doi.org/10.1155/2020/2678374
  * [9] Bean NG, Pollett PK, Taylor PG. Quasistationary distributions for level-dependent quasi-birth-and-death processes. Commun Stat Stoch Models. 2000;16(5):511–541. doi.org/10.1080/15326340008807602
  * [10] Bright L, Taylor PG. Calculating the equilibrium distribution in level dependent quasi-birth-and-death processes. Commun Stat Stoch Models. 1995;11:497–525. doi.org/10.1080/15326349508807357
  * [11] Bright LW, Taylor PG. Equilibrium distributions for level-dependent quasi-birth-and-death processes. In: Chakravarthy SR, Alfa AS, eds. Matrix Analytic Methods in Stochastic Models: Proceedings of the 1st International Conference. New Jersey: Marcel Dekker; 1997:359–375.
  * [12] Britton T, Ball F, Trapman P. A mathematical model reveals the influence of population heterogeneity on herd immunity to SARS-CoV-2. Science 2020;369:846–849. doi.org/10.1126/science.abc6810
  * [13] Chai X, Jiang T, Li L, Xu W, Lin L. On a many-to-many matched queueing system with flexible matching mechanism and impatient customers. J Comput Appl Math. 2022;416:114573. doi.org/10.1016/j.cam.2022.114573
  * [14] Chalub FACC, Gómez-Corral A, López-García M, Palacios-Rodríguez F. A Markov chain model to investigate the spread of antibiotic-resistant bacteria in hospitals. Stud Appl Math. 2023;151:1498–1524. doi.org/10.1111/sapm.12637
  * [15] Dayar T, Sandmann W, Spieler D, Wolf V. Infinite level-dependent QBD processes and matrix-analytic solutions for stochastic chemical kinetics. Adv Appl Probab. 2011;43:1005-1026. doi.org/10.1239/aap/1324045696
  * [16] Dayar T, Orhan MC. Steady-state analysis of a multiclass MAP/PH/c queue with acyclic PH retrials. J Appl Probab. 2016;53(4):1098–1110. doi.org/10.1017/jpr.2016.67
  * [17] De Nitto Personè V, Grassi V. Solution of finite QBD processes. J Appl Probab. 1996;33:1003–1010. doi.org/10.2307/3214981
  * [18] Diao J, Stark TL, Liberles DA, O’Reilly MM, Holland BR. Level-dependent QBD models for the evolution of a family of gene duplicates. Stoch Models. 2020;36(2):285-311. doi.org/10.1080/15326349.2019.1680296
  * [19] Gaver DP, Jacobs PA, Latouche G. Finite birth-and-death models in randomly changing environments. Adv Appl Probab. 1984;16:715–731. doi.org/10.2307/1427338
  * [20] Gómez-Corral A, López García M. Extinction times and size of the surviving species in a two-species competition process. J Math Biol. 2012;64:255–289. doi.org/10.1007/s00285-011-0414-8
  * [21] Gómez-Corral A, López García M. Maximum queue lengths during a fixed time interval in the M/M/c retrial queue. Appl Math Comput. 2014;235:124–136. doi.org/10.1016/j.amc.2014.02.074
  * [22] Gómez-Corral A, López-García M. Perturbation analysis in finite LD-QBD processes and applications to epidemic models. Numer Linear Algebra Appl. 2018;25:2160. doi.org/10.1002/nla.2160
  * [23] Gómez-Corral A, López-García M, Lopez-Herrero MJ, Taipe D. On first-passage times and sojourn times in finite QBD processes and their applications in epidemics. Mathematics 2020;8:1718. doi.org/10.3390/math8101718
  * [24] Gómez-Corral A, Langwade J, López-García M, Molina-París C. Sufficient conditions for regularity, positive recurrence and absorption in level-dependent QBD processes and related block-structured Markov chains. Math Meth Appl Sci. 2023;46(6):6756-6766. doi.org/10.1002/mma.8938
  * [25] Gómez-Corral A, Palacios-Rodríguez F, Rodríguez-Bernal MT. On the exact reproduction number in SIS epidemic models with vertical transmission. Comp Appl Math. 2023;42:291. doi.org/10.1007/s40314-023-02424-5
  * [26] Gómez-Corral A, Lopez-Herrero MJ, Taipe D. A Markovian epidemic model in a resource-limited environment. Appl Math Comput. 2023;458:128252. doi.org/10.1016/j.amc.2023.128252
  * [27] Hitchcock SE. Extinction probabilities in predator-prey models. J Appl Probab. 1986;23:1-13. doi.org/10.2307/3214112
  * [28] Ho LST, Xu J, Crawford FW, Minin VN, Suchard MA. Birth/birth-death processes and their computable transition probabilities with biological applications. J Math Biol. 2018;76:911-944. doi.org/10.1007/s00285-017-1160-3
  * [29] Iglehart DL. Multivariate competition processes. Ann Math Stat. 1964;35:350-361. doi.org/10.1214/aoms/1177703758
  * [30] Javier K, Fralix B. On the study of the running maximum and minimum level of level-dependent quasi-birth-death processes and related models. J Appl Probab. 2023;60(1):14–29. doi.org/10.1017/jpr.2022.22
  * [31] Jeganathan K, Abdul Reiyas M, Prasanna Lakshmi K, Saravanan S. Two server Markovian inventory systems with server interruptions: Heterogeneous vs. homogeneous servers. Math Comput Simul. 2019;155:177–200. doi.org/10.1016/j.matcom.2018.03.001
  * [32] Kermack WO, McKendrick AG. A contribution to the mathematical theory of epidemics. Proc R Soc Lond. A 1927;115:700–721. doi.org/10.1098/rspa.1927.0118
  * [33] Kharoufeh JP. Level-dependent quasi-birth-and-death processes. In: Cochran JJ, Cox Jr LA, Keskinocak P, Kharoufeh JP, Smith JC, eds. Wiley Encyclopedia of Operations Research and Management Science. John Wiley and Sons; 2011:1–9. doi.org/10.1002/9780470400531.eorms0460
  * [34] Kulkarni VG. Modeling and Analysis of Stochastic Systems, 3rd ed. Boca Raton: Chapman & Hall/CRC Press; 2017.
  * [35] Latouche G, Ramaswami V. A logarithmic reduction algorithm for quasi-birth-and-death processes. J Appl Probab. 1993;30:650–674. doi.org/10.2307/3214773
  * [36] Latouche G, Ramaswami V. Introduction to Matrix Analytic Methods in Stochastic Modeling. Philadelphia: ASA-SIAM; 1999.
  * [37] Lefèvre C, Simon M. SIR-type epidemic models as block-structured Markov processes. Methodol Comput Appl Probab. 2020;22:433–453. doi.org/10.1007/s11009-019-09710-y
  * [38] Liu B, Zhao YQ. Analyzing retrial queues by censoring. Queueing Syst. 2010;64:203–225. doi.org/10.1007/s11134-009-9163-4
  * [39] Mandjes M, Taylor P. The running maximum of a level-dependent quasi-birth-death process. Probab Eng Inf Sci. 2016;30(2):212–223. doi.org/10.1017/S026996481500039X
  * [40] Neuts MF. Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach. Baltimore: The Johns Hopkins University Press; 1981.
  * [41] Neuts MF, Li JM. An algorithmic study of S-I-R stochastic epidemic models. In: Heyde CC, Prohorov YV, Pyke R, Rachev ST, eds. Athens Conference on Applied Probability and Time Series Analysis, Lecture Notes in Statistics, vol. 114. New York: Springer; 1996:295–306. doi.org/10.1007/978-1-4612-0749-8_21
  * [42] Phung-Duc T, Masuyama H, Kasahara S, Takahashi Y. A simple algorithm for the rate matrices of level-dependent QBD processes. In: Wang J, Yue W, eds. Proceedings of the 5th International Conference on Queueing Theory and Network Applications. New York: Association for Computing Machinery; 2010:46–52. doi.org/10.1145/1837856.1837864
  * [43] Phung-Duc T, Kawanishi K. Performance analysis of call centers with abandonment, retrial and after-call work. Perform Eval. 2014;80:43–62. doi.org/10.1016/j.peva.2014.03.001
  * [44] Phung-Duc T. Asymptotic analysis for Markovian queues with two types of nonpersistent retrial customers. Appl Math Comput. 2015;265:768–784. doi.org/10.1016/j.amc.2015.05.133
  * [45] Ramaswami V, Taylor PG. Some properties of the rate operators in level dependent quasi-birth-and-death processes with countable number of phases. Commun Stat Stoch Model. 1996;12:143–164. doi.org/10.1080/15326349608807377
  * [46] Reuter GEH. Competition processes. In: Neyman J, ed. Proceedings of the 4th Berkeley Symposium on Mathematical Statistics and Probability, vol. II: Contributions to probability theory. Berkeley: University of California Press; 1961:421–430.
  * [47] Ridler-Rowe CJ. On competition between two species. J Appl Probab. 1978;15:457-465. doi.org/10.2307/3213109
  * [48] Ridler-Rowe CJ. Extinction times for certain predator-prey processes. J Appl Probab. 1988;25:612-616. doi.org/10.2307/3213988
  * [49] Sanga SS, Charan GS. Fuzzy modeling and cost optimization for machine repair problem with retrial under admission control F-policy and feedback. Math Comput Simul. 2023;211:214–240. doi.org/10.1016/j.matcom.2023.03.036
  * [50] Santana LM, Ganesan S, Bhanot G. A quasi birth-and-death model for tumor recurrence. J Theor Biol. 2019;480:175–191. doi.org/10.1016/j.jtbi.2019.07.017
  * [51] Saravanan V, Poongothai V, Godhandaraman P. Performance analysis of a multi server retrial queueing system with unreliable server, discouragement and vacation model. Math Comput Simul. 2023;214:204–226. doi.org/10.1016/j.matcom.2023.07.008
  * [52] Stewart WJ. Introduction to the Numerical Solution of Markov Chains. Princeton: Princeton University Press; 1994.
  * [53] Takine T. On level-dependent QBD processes with explosive state space. Queueing Syst. 2022;100:353–355. doi.org/10.1007/s11134-022-09796-1
  * [54] van den Driessche P. Reproduction numbers of infectious disease models. Infect Dis Model. 2017;2(3):288–303. doi.org/10.1016/j.idm.2017.06.002