On First-Passage Times and Sojourn Times in Finite QBD Processes and Their Applications in Epidemics

: In this paper, we revisit level-dependent quasi-birth-death processes with ﬁnitely many possible values of the level and phase variables by complementing the work of Gaver, Jacobs, and Latouche (Adv. Appl. Probab. 1984), where the emphasis is upon obtaining numerical methods for evaluating stationary probabilities and moments of ﬁrst-passage times to higher and lower levels. We provide a matrix-analytic scheme for numerically computing hitting probabilities, the number of upcrossings, sojourn time analysis, and the random area under the level trajectory. Our algorithmic solution is inspired from Gaussian elimination, which is applicable in all our descriptors since the underlying rate matrices have a block-structured form. Using the results obtained, numerical examples are given in the context of varicella-zoster virus infections.


Introduction
In this paper, we shall be concerned with time-homogeneous continuous-time Markov chains X = {(I(t), J(t)) : t ∈ [0, ∞)} that take values in a bivariate state space S = {(i, j) : i ∈ {0, . . . , N}, j ∈ {0, . . . , M i }}, for integers M i ∈ N 0 with i ∈ {0, . . . , N} and N ∈ N, and have an irreducible infinitesimal generator Q with block-tridiagonal representation where the sub-matrices Q i,i are of dimension (1 + M i ) × (1 + M i ), for i ∈ {max{0, i − 1}, i, min{i + 1, N}}, i.e., from any state (i, j), process X can access only states (i , k) with |i − i| ≤ 1 and k ∈ M i , whence I(t) can change its value by {−1, 0, 1} only. These Markov chains are referred to as finite quasi-birth-death (QBD) processes, being its first variable I(t) called the level, and the second one J(t) the phase; see e.g. Gaver et al. [1]. They are capable of modeling, in a unified and algorithmically tractable manner, relevant systems, such as queueing models (Baumann and Sandmann [2]; Gun and Makowski [3]; Perel and Yechiali [4]; Ye and Li [5]), communication networks (Artalejo and Gómez-Corral [6]; Artalejo and Lopez-Herrero [7]), reliability and manufacturing systems (Chakravarthy and Gómez-Corral [8]; Moghaddass et al. [9]), and epidemic models (Amador and Gómez-Corral [10]; Artalejo et al. [11]; Economou et al. [12]; Gamboa and Lopez-Herrero [13]), among others. The properties of finite QBD processes have been studied extensively in continuous and discrete time. For a comprehensive discussion, see, e.g., the monograph of Latouche and Ramaswami (Chapter 10 in Reference [14]) and references therein. Most of the theoretical work so far carried out in this setting is concerned with the stationary probability vector. More concretely, Hajek [15] observes that, under certain conditions for the underlying transition matrices, the stationary vector of a finite level-independent QBD process has a mixed matrix-geometric form; for a related numerical method that exploits invariant subspace computations, see the paper by Akar et al. [16]. Gaver et al. (Section 2 in Reference [1]) use linear level reduction, which is a refined probabilistic justification for block-Gaussian elimination (Latouche and Ramaswami (Remark 10.1.6 in Reference [14]); Stewart [17]). Ye and Li [5] (see also Latouche and Ramaswami (Section 10.2 in Reference [14])) propose an efficient computational approach, termed folding algorithm, which is seen to be a form of odd-even reduction, and Naoumov [18] presents an algebraically matrix-multiplicative approach for a finite QBD process. LU-type and UL-type RG-factorizations are used by Li and Cao (Section 3.1 in Reference [19]) in an attempt to clarify the relations among those results in Gaver et al. [1], Hajek [15], and Naoumov [18]. Elhafsi and Molle [20] present a matrix-geometric expression for the stationary vector of a finite level-independent QBD process in terms of a set of rate matrices and compare their solution with several existing solution methods, including the approaches of de Nitto Persone and Grassi [21], Gun and Makowski [3], and Hajek [15], and a more general solution of Le Boudec [22] for upper block Hessenberg matrices. Other approaches are related to the group inverse (da Silva Soares and Latouche [23]) and the matrix continued fraction algorithm (Baumann and Sandmann [24]).
Despite their simplicity and their success in stochastic modeling, there have been few contributions dealing with random descriptors different from the stationary probability vector in the theory of finite QBD processes.  in Reference [1]) use Laplace-Stieltjes transforms of first-passage times from level l(i − 1) to level l(i), where the i th level is defined as l(i ) = {(i , j) : j ∈ {0, . . . , M i }}, to derive a product-form expression for the Laplace-Stieltjes transform of the first-passage time from level l(i − 1) to higher levels l(i ), with i ∈ {i, . . . , N}, and similar results for the first-passage times to lower levels. An up-integral functional is defined in the paper by Li and Cao [19] (Section 4) and explicit expressions for the Laplace transforms of its conditional distributions and moments are written in terms of the underlying Rand G-measures. Li and Sheng [25] propose and implement a generalized folding algorithm for transient analysis of a finite QBD process, which entails the solution of xQ = a, where a is a predetermined boundary condition vector, and present comprehensive numerical examples on transient queueing behavior with correlated arrivals. Based on the Rand G-measures, Shin [26] derives an algorithm for computing the fundamental matrix of a finite QBD process with absorbing states and shows how the resulting algorithm can be applied to absorption probabilities, sojourn time distributions, the busy period and the stationary vector in the multiprogramming model discussed by Neuts [27] (pp. 245-253) and the M/M/c retrial queue with finite orbit capacity (Artalejo and Gómez-Corral (Chapter 4 in Reference [28])). Recently, Gómez-Corral and López-García [29] use matrix calculus and exploit the specific structure of the infinitesimal generator Q in (1) to provide the sensitivities and the elasticities of finite QBD processes with respect to model parameters, with special emphasis on first-passage times to level l(0) and hitting probabilities, the maximum level visited by the process before reaching states in level l(0) and the stationary vector, and illustrate the approach by means of applying multi-type versions of the susceptible-infective (SI) and susceptible-infective-susceptible (SIS) epidemic models to the spread of antibiotic-sensitive and antibiotic-resistant bacterial strains in a hospital ward.
The objective of this paper was to complement the work of Gaver et al. [1] by showing here that, in the spirit of algorithmic probability, the fundamental ideas of Gaver et al. [1] extend beyond the stationary distribution and first-passage times to higher and lower levels, and apply to more complex descriptors arising from epidemic models. The methods we shall follow in this paper are suggested by some related work of Amador et al. [30], Baumann and Sandmann [31], Lefèvre and Simon [32], and Neuts [33], where, for instance, the random length of an outbreak and the final size of the epidemic in SIR-type and SIS-type epidemic models correspond to the absorption time and the hitting probabilities in a suitably defined absorbing QBD process, and related block-structured Markov processes. In order to simplify the presentation throughout the paper, we shall assume that (1, j) is the initial state of process X , for phases j ∈ {0, . . . , M 1 }, which is related to an invasion time in the spread of a disease.
The paper is organized as follows. In Section 2, we derive recursive expressions for the hitting probabilities of restricted versions of process X . In Section 3, we characterize the probability laws of the number of upcrossings from level l(K − 1) to states in l(K), for a predetermined threshold K ∈ {2, . . . , N}, and the total time spent at states in higher levels before reaching states in level l(0). In Section 4, we study the random area under the trajectory of the integer-valued random variable I(t) before process X visits states in l(0) for the first time. Numerical examples are given in Section 5 in the setting of varicella-zoster virus infections, and concluding remarks follow in Section 6.
In order to derive expressions for the restricted hitting probabilities {P (1,j),K (n) : n ∈ {0, . . . , M K }}, it is useful to consider an auxiliary absorbing process where 0 and K are assumed to be absorbing states and K−1 i=1 l(i) is a class of transient states. The infinitesimal generator of X (K) has the block-structured form where the column vector 0 a is the null vector of order a, T denotes transposition, the sub-matrix T(K) is obtained from Q by removing rows and columns associated with states in l(0) and N i=K l(i), the column vectors t 0 (K) and t K (K) have the form and the column vector 1 a is the unit vector of order a. It is clear that, starting from states in K−1 i=1 l(i), the process X (K) defined by (2) is the restriction of process X to the states in K−1 i=1 l(i) before the first entrance into states in l(0) ∪ N i=K l(i). In particular, it is found that where e a (b) is a row vector of order a, in which entries are all equal to 0, except the bth entry which equals 1. The proof of (3) is based on the fact that the transition function P(X(t) = K|X(0) = (1, j)) of process X (K) is given by which yields P(η (1,j),K <η (1,j),0 ) by taking the limit as t → ∞, since the sub-matrix T(K) is stable and the real parts of its eigenvalues are all strictly negative. By appealing to statement (b) of Lemma 1 in the paper by Gaver et al. [1], it is noted that the matrix −T −1 (K) records expected total times that process X spends at states in K−1 i=1 l(i) before X visits states in either l(0) or l(K) for the first time, provided that X starts from any initial state in K−1 i=1 l(i). Furthermore, a slight modification of statement (c) of Lemma 1 in Reference [1] permits to establish that, by denoting by T K (K) the matrix recording infinitesimal rates associated with transitions to states in l(K) from K−1 i=1 l(i), the matrix −T −1 (K)T K (K) consists of first passage probabilities to level l(K), given that X starts from any state in K−1 i=1 l(i). This means that the (1 + j)th row of −T −1 (K)T K (K) consists of the restricted hitting probabilities P (1,j),K (n), for phases n ∈ {0, . . . , M K }, and consequently, recursive relations for these probabilities can be derived from Theorem 4.2.4 of the book by Hunter [34] by observing that T(K) can be written in terms of the sub-matrix T(K − 1). Theorem 1. For a predetermined integer K ∈ {2, . . . , N} and an initial state (1, j) with j ∈ {0, . . . , M 1 }, the restricted hitting probability P (1,j),K (n) is specified by for phases n ∈ {0, . . . , Moreover, the sub-matrices V 1,2 (K), for integers K ∈ {3, . . . , N}, are given by and satisfy This has the following immediate consequence.

Remark 2.
For the sake of completeness, we note that the restricted hitting probability P (i,j),K (n), for any initial state (i, j) with j ∈ {0, . . . , M i } and i ∈ {2, . . . , K − 1}, can be similarly derived from the probabilistic interpretation of the matrix −T −1 (K)T K (K), in such a way that for phases n ∈ {0, . . . , M K }, where we denote by (A) l,l the (l, l )-th entry of matrix A. Therefore, a general version of Algorithm A1 (Appendix A) can be written by exploiting the block-structured form of −T −1 (K) in (5), which is left to the interested reader.

Number of Upcrossings and Sojourn Times
We shall examine here the number µ (1,j),K of upcrossings from states in level l(K − 1) to level l(K) (i.e., the number of times the level variable I(t) increases its value from K − 1 to K) and the total time γ (1,j),K that process X spends at states in N i=K l(i) before the first entrance of X into states in l(0), provided that (1, j) is the initial state.
We note that these descriptors can have different interpretations depending on the particular QBD process of interest. For example, in an epidemic model (see Section 5), if the level l(K) consists of states in which exactly K infected hosts are present in the population, the number µ (1,j),K of upcrossings is equivalent to the number of times the threshold number K of infected hosts is reached during a given outbreak. This can be relevant when, in particular, K represents a threshold level of infectives at which some mitigation strategies are implemented, or when K represents a maximum number of infectives for whom there is actual capacity in the healthcare system to deal with. In these scenarios, it is also relevant to quantify the total amount of time γ (1,j),K that the epidemic stays at or above level l(K), which would quantify the time during which the mitigation strategies are in place, as a measure of cost, or the time period during which the healthcare capacity is exceeded.
To begin with, we establish the connection between µ (1,j),K and the numbers µ (i,l),K of upcrossings from level l(K − 1) to level l(K) before a visit to l(0), provided that (i, l) is the initial state, for phases l ∈ {0, . . . , M i } and integers i ∈ {K − 1, K}.

Proof.
Let us consider what must happen in order for the event {µ (1,j),K = m}, for integers m ∈ N, to occur. Conditioning on the first entrance of X into l(K) before a visit to l(0), the process X is forced to return to level l(K), starting from states in l(K), by means of m − 1 upcrossings from l(K − 1) to l(K) before a visit to l(0) and a number of downcrossings from l(K + 1) to l(K), which are not relevant for the event {µ (1,j),K = m} to occur. Then, we may obtain (7) by conditioning on the state (K, n) visited by X in its first entrance into l(K), starting from (1, j), avoiding l(0). Note that, in the special case m = 0, process X must enter level l(0), starting from (1, j), before moving up to l(K), which occurs with probability Equation (8) is similarly obtained by conditioning on the state (K − 1, l) visited by process X in its first entrance into level l(K − 1), starting from state (K, n), whereas (9) is derived from the above argument by replacing the initial state (1, j) by (K − 1, l). This completes the proof.
In evaluating the sequence {Q (K,n),K−1 (l) : l ∈ {0, . . . , M K−1 }} of first-passage probabilities, we may observe that they correspond to the absorption probabilities into states in level l(K − 1) in an auxiliary absorbing process X (K) that, starting from (K, n), takes values in the state space S(K) = l(K − 1) ∪ N i=K l(i) and has the following block-structured infinitesimal generator: , and the sub-matrix T(K) is obtained from Q by removing rows and columns associated with states in K−1 i=0 l(i). Specifically, it is readily seen that the matrix − T −1 (K) T K−1 (K) consists of first-passage probabilities to states in level l(K − 1), provided that process X starts from states in N i=K l(i). This yields for phases n ∈ {0, . . . , is the Kronecker's delta, and the matrices U 1,2 (K) and U 2,1 (K) are given by . Algorithm A2 (Appendix B) uses the above recursive relations and iteratively computes probabilities ,K can be derived from (7)-(9) and is given by for |z| ≤ 1, where the generating functions ϕ (K,n),K (z), for phases n ∈ {0, . . . , M K }, satisfy for |z| ≤ 1 and phases l ∈ {0, . . . , M K−1 }. By differentiating these equalities successively with respect to z and setting z = 1, we may obtain the column vectors m where the (1 + j, 1 + n)-th entry of matrix P 1,K is given by P (1,j),K (n), for phases j ∈ {0, . . . , M 1 } and n ∈ {0, . . . , M K }, and the column vectors m (k) K are obtained iteratively from the equality In this equality, the (1 + n, 1 + l)-th entry of P K,K−1 and the (1 + l, 1 + n)-th entry of P K−1,K are specified by the probabilities Q (K,n),K−1 (l) in (10) and P (K−1,l),K (n) in (6), respectively, for phases l ∈ {0, . . . , M K−1 } and n ∈ {0, . . . , M K }.
Let us now turn to the probability law of the random variable γ (1,j),K . The probability law of γ (1,j),K has the discrete contribution , and a continuous contribution, which is linked the restricted Laplace-Stieltjes transform By conditioning on the first state visited in l(K) before a visit to l(0), it is seen that for Re(s) ≥ 0, where γ (K−1,l),K is the total time that process X spends at states in N i=K l(i) before a visit to l(0), provided that (K − 1, l) is the initial state, and θ (K,n),K−1 (s; l) is the restricted Laplace-Stieltjes transform of the first-passage time ζ (K,n),K−1 (l) to level l(K − 1), provided that (K, n) is the initial state and the first state visited in level Theorem 3 tells us how to compute the column vectors Θ K,K−1 (s; l), for l ∈ {0, . . . , M K−1 }, with (1 + n)th entry given by θ (K,n),K−1 (s; l), for n ∈ {0, . . . , M K }, in terms of the Laplace-Stieltjes transforms θ (i,n),K−1 (s; l) of the first-passage time ζ (i,n),K−1 (s; l) to level l(K − 1), provided that (i, n) is the initial state and the first state visited by process X in level l(K − 1) is (K − 1, l), for phases n ∈ {0, . . . , M K } and l ∈ {0, . . . , M K−1 }, and integers i ∈ {K + 1, . . . , N}. From now on we denote by q (i,j),(i ,j ) the (1 + j, 1 + j )-th entry of the sub-matrix Q i,i and by −q (i,j) the (1 + j, 1 + j)-th entry of sub-matrix Q i,i , for phases j ∈ {0, . . . , M i } and j ∈ {0, . . . , M i }, and integers i ∈ {max{0, i − 1}, i, min{i + 1, N}}. Theorem 3. For phases l ∈ {0, . . . , M K−1 }, the column vector Θ K,K−1 (s; l) is determined from the equality for , for i ∈ {K + 1, . . . , N − 1}, and the matrices G i (s) and the column vectors g i (s; l) satisfy for integers i ∈ {K, . . . , N}.
From (14) and (15), it is found that the following relations hold for the column vectors Γ (k) , so that the kth moment of the random variable γ (1,j),K is given by for phases j ∈ {0, . . . , M 1 }. We only briefly indicate here that the matrix Θ

Area under the Level Trajectory
Consider the trajectory of the level variable I(t), for which (I(0), J(0)) = (i, j) with (i, j) ∈ S \ l(0), and define the area under the level trajectory as the stochastic integral where, in a similar manner toη (1,j),0 , the random variableη (i,j),0 denotes the first-passage time to states in level l(0), provided that process X starts from state (i, j).

Remark 3.
In the setting of epidemic models, α (1,j),0 can be thought of as the total amount of individual time units of infection during the course of the epidemic. Its probability law, which is inherently linked to the duration η (1,j),0 of an outbreak, is directly related to the cost of the epidemic, as defined by Downton [35] and Gani and Jerwood [36]; for a related work, also see Ball [37].

An Application to Varicella-Zoster Virus Infections
We illustrate our analysis by focusing here on a mathematical model for varicella-zoster virus in a nursing home. Varicella-zoster virus (VZV) causes chickenpox (varicella) and shingles (herpes zoster). In childhood, chickenpox disease produces itchy blisters but rarely causes serious problems. However, in adults who have not suffered the disease as children, chickenpox can lead to serious complications. After primary infection, VZV can remain dormant within dorsal root ganglia for life, and it can reactivate depending on a number of factors including the host immune system. Herpes zoster is the reactivation of VZV. Factors associated with recurrent disease include, among others, aging, immunosupression, intrauterine exposure to VZV, or having had varicella at a young age. Shingles is characterized by a rash of blisters, and it can be very painful but is not life-threatening. Varicella is highly contagious because the transmission occurs by direct contact, by air (coughing and sneezing) and by areosolization of virus from skin lesions. A person with active herpes zoster can transmit the VZV, through direct contact, to a person who never had chickenpox; we refer the reader to Reference [38] for a flow diagram showing the different stages for an individual during an outbreak.
Our interest is in illustrating the applicability of our methodology by focusing on a Markovian model for the spread of VZV infection in a closed population, such as a nursing home for elderly people, subject to repeated outbreaks. The mathematical model assumes an homogeneous and finite population with constant size N, where residents who die or abandon the nursing home are immediately replaced by the admission of a new resident, reflecting high demand for nursing resources. Newly admitted residents are assumed to be either susceptible or asymptomatic, but never showing signs or symptoms of varicella or herpes zoster. We define the continuous-time Markov chain Y = {(I(t), Z(t), A(t)) : t ≥ 0}, where I(t), Z(t) and A(t) represent the number of varicella infectious (i.e., residents infected and infectious with varicella), zoster infectious (i.e., residents showing dermal rash of herpes zoster) and asymptomatic (i.e., residents with dormant VZV) residents at time t. Since we assume constant population size, it is clear that the number of susceptible residents is given by We assume that recovery times, as well as transmission times and the inter-event times for all other events in this process are exponentially distributed and mutually independent, leading to the Markovian hypothesis. We make the following further assumptions: • Varicella infectious residents recover at rate γ. • Zoster infectious residents recover at rate .

•
Residents are removed (either by dying or abandoning the nursing facility) at rate µ.

•
Newly admitted residents are (asymptomatically) infected with probability q.

•
Transmission between susceptible and varicella infectious residents occurs at rate β, while transmission between susceptible and zoster infectious residents occurs at rate η.

•
Reactivation of zoster for asymptomatic residents occurs at rate δ.
These hypotheses lead to a number of state transitions for process Y with infinitesimal rates if (i , z , a ) = (i, z, a + 1) and a + i + z < N, z, a) and a + i + z < N, aδ, if (i , z , a ) = (i, z + 1, a − 1) and a > 0, for any states (i, z, a) and (i , z , a ) with (i , z , a ) = (i, z, a). We point out that some events are omitted here since they are unobserved by process Y, meaning that they do not cause a change of state (e.g., if an asymptomatic resident dies and a newly admitted asymptomatic resident arrives, which occurs with rate aqµ but does not affect the state of the process). The continuous-time Markov chain Y is defined on the state space Even though process Y is three-dimensional (while process X in Section 1 is two-dimensional), one can obtain a QBD formulation for process Y by organizing S Y in levels as This organization of states in terms of levels and sub-levels, as well as the transitions described above, leads to the following block-tridiagonal representation for the infinitesimal generator of Y: We note that this tridiagonal structure means that process Y is a finite QBD process, so that our arguments in Sections 2-5 directly apply. Variable I(t) in Y is the level, representing the number of varicella infectious individuals at any given time. The phase J(t) in our original process X in Section 1 has a bivariate representation here J(t) = (Z(t), A(t)). These phases are organized in terms of sub-levels as described above. In particular, one orders sub-levels inside l(i) as while states inside each sub-level l(i, z) are ordered as Thus, the first phase (denoted by j = 0 in Section 1) in level l(i) is equivalent to the pair (z, a) = (0, 0), while the last phase is the j = ( (N+1−i)(N+2−i) 2 − 1) th one which corresponds to (z, a) = (N − i, 0), since the number of phases in each level l(i) (denoted by M i + 1 in Section 1) is given by The organization of phases for each level l(i) in terms of sub-levels l(i, z) means that sub-matrices A i,i are given as follows: • For i ∈ {1, . . . N − 1}, sub-matrix A i,i−1 is associated with jumps of process Y from states in level l(i) to states in level l(i − 1). In particular, it is given by Matrices B i,i−1 (z, z) consist of transition rates from sub-level l(i, z) to sub-level l(i − 1, z) and are given in Appendix E. We can deal separately with the case i = N, for which • For i ∈ {0, . . . N − 1}, sub-matrix A i,i is related to transitions from level l(i) to level l(i) and has the structured form where matrices B i,i (z, z ) consist of transition rates from states in sub-level l(i, z) to states in sub-level l(i, z ), for z ∈ {z − 1, z, z + 1}, and are detailed in Appendix E. It is worth noting the particular case i = N, which leads to the single-element matrix A N,N = −N(µ + γ).

•
For i ∈ {1, . . . , N − 1}, sub-matrix A i,i+1 corresponds to transitions from states in level l(i) to states in level l(i + 1) and is given by where matrices B i,i+1 (z, z) contain transition rates from states in sub-level l(i, z) to states in sub-level l(i + 1, z), and are given in Appendix E. It is worth noting the particular case i = 0, which leads to matrix For process Y, the parameter values (γ, , µ, q, β, δ) in our numerical experiments have been chosen according to the following assumptions: We consider a long-term nursing facility where the average length of stay of residents is 800 days (i.e., 2.202 years) [39,40] and set µ = 1/2.202 = 0.454 years −1 . (ii) Residents recover from varicella after an average of 7 days [38]; thus, we set γ = 52.142 years −1 . (iii) The average recovery time for herpes zoster is 20 days [38], so we set = 18.250 years −1 . (iv) Based on existing data for the reactivation rate for zoster [38,41,42], we set δ = 0.065 years −1 .
(v) Since it is estimated that around 90% of adults in the USA carry VZV and are at risk of developing herpes zoster [41], and a similar percentage is estimated for other countries in temperate climates [43], we explore values q ∈ (0.8, 0.9) in our numerical results. (vi) Varicella is highly infectious, but the particular transmissibility depends on a number of factors [38,[41][42][43] (e.g., amount of contact in the nursing facility, control strategies in place, such as cohorting of staff and residents, etc). Thus, we will explore two parameter values β ∈ {38.0, 190.0}. In particular, we set β = 38.0 to broadly represent a situation where the nursing home has a range of control strategies in place when the outbreak starts, while β = 190.0 represents the situation where no control strategies are in place. Note that the value β = 190.0 is selected by assuming a basic reproduction number R 0 = β/γ = 3.644, in accordance with estimated results for VZV transmission displayed in Reference [44].
Finally, herpes zoster is transmitted significantly less than varicella itself [41], but its particular transmission rate has not been precisely quantified to the best of our knowledge, so we explore values η ∈ (0, β) in our numerical results.
For purely illustrative purposes, we consider a nursing home with N = 50 residents and focus on four different scenarios which vary in potential initial conditions and which allow us to carry out a sensitivity analysis of the summary statistics described in Sections 2-4 on the parameter values discussed above. In particular, these scenarios are specified as follows: • Scenario 1. Outbreak initiated by a varicella infectious resident and prevalence of asymptomatics given by q = 0.9, i.e., initial condition (i, z, a, s) = (1, 0, 45, 4). • Scenario 2. Outbreak initiated by a zoster infectious resident and prevalence of asymptomatics given by q = 0.9, i.e., initial condition (i, z, a, s) = (0, 1, 45, 4). • Scenario 3. Outbreak initiated by a varicella infectious resident and prevalence of asymptomatics given by q = 0.8, i.e., initial condition (i, z, a, s) = (1, 0, 40,9). • Scenario 4. Outbreak initiated by a zoster infectious resident and prevalence of asymptomatics given by q = 0.8, i.e., initial condition (i, z, a, s) = (0, 1,40,9).
For each scenario, our interest is in computing, by means of the analysis in Sections 2-4, the following summary statistics:

•
The probability P(η (i,z,a),K <η (i,z,a),0 ) of reaching level l(K) (that is, K varicella infectious residents in the nursing home) before reaching level l(0) (end of the varicella outbreak), for different values of K.

•
The mean number E[µ (i,z,a),K ] of upcrossings to level l(K).

•
The proportion of expected time E[γ (i,z,a),K ]/E[η (i,z,a),0 ] that the process Y stays at or above level l(K) (so, with K or more varicella infectious residents).

•
The ratio E[α (i,z,a),0 ]/(N · E[η (i,z,a),0 ]), which represents the area under the curve of varicella infectives normalized by the total amount of time units that N residents stay in the nursing home during the varicella outbreak.
In Figure 1, we plot the probability P(η (i,z,a),K <η (i,z,a),0 ) of reaching a particular threshold value K of varicella infectives (which might represent, for example, a particular number of infectives which would trigger the implementation of particular control strategies in the nursing home), for K ∈ {3, 5}, β ∈ {38.0, 190.0} and versus different values of η ∈ (0, β), for Scenarios 1-4. It is clear that this probability, for those scenarios related to a varicella infective starting the outbreak (Scenarios 1 and 3), does not significantly depend on the reactivation rate η of zoster infection as one would expect (since the initially infectious resident is a varicella infective). For Scenarios 2 and 4, in which a zoster infective starts the outbreak, the probability of reaching K varicella infective residents sharply increases for initially increasing values of η (better allowing for the initially infective zoster to transmit before recovering), and then slowly stabilizes. On the other hand, it is worth noting that, due to the highly infectious nature of VZV, the maximum number of varicella infective residents reached before the end of the outbreak -not explicitly computed here, but which clearly affects the probability of reaching the particular threshold levels l(5) and l(10)directly depends on the number of susceptible residents in the nursing home when the outbreak starts. This number of susceptibles directly depends on parameter q (q = 0.9 in Scenarios 1 and 2, leading to 4 susceptible residents; q = 0.8 in Scenarios 3 and 4, leading to 9 susceptible residents) and directly explains the asymptotic behavior of curves observed in Figure 1 for K = 3 and 5. In Figure 2, we plot the mean number of upcrossings to level l(K), for K = 3, as a function of η. We note that we always obtain, regardless of the parameters varied here, a relatively small mean number of upcrossings (e.g., less than 1.5 in all scenarios). We point out that a mean number of upcrossings significantly below 1 means that the threshold number K = 3 of varicella infective residents is not actually reached in most of the stochastic realizations of process Y. This happens, for example, when a zoster infective resident starts the outbreak and there are control measures in place, so that the transmissibility of herpes zoster is low (i.e., relatively small values of η in Figure 2). For large enough values of η, or for those scenarios in which a varicella infective resident starts the outbreak, a mean number of upcrossings relatively close to 1 illustrates how level l(K) with K = 3 is in fact reached during the early stages of the outbreak (i.e., while the number of infective residents is increasing over early times, moving towards its peak beyond K = 3). Values of E[µ (i,z,a),K ] slightly above 1 illustrate the stochasticity of the epidemic process under analysis and might be explained by stochastic contributions from some particular realizations of the process (for example, if process Y reaches K = 3 varicella infective residents, and one of these residents recovers before another infection immediately occurs, moving the process to level l(2) before visiting level l(3) again, and causing then two upcrossings). In Figure 3, we plot the proportion of expected time E[γ (i,z,a),K ]/E[η (i,z,a),0 ] that the process Y stays at or above level l(K). This might be of interest since, for example, if K varicella infective residents represent an alert or warning level from which the nursing home needs to incur in some particular costs per unit time (e.g., particular enhanced control measures), the quantity E[γ (i,z,a),K ]/E[η (i,z,a),0 ] would become an implicit quantification of this cost. We plot E[γ (i,z,a),K ]/E[η (i,z,a),0 ] for a nursing home with some control strategies initially in place (so β = 38.0), and for K = 3 for illustrative purposes. The proportion of expected time spent above level l(K) directly depends on who starts the outbreak (from either a varicella infective resident or a zoster infective resident) in the first place, the number of susceptible residents in the nursing home (q = 0.9 versus q = 0.8) and the zoster infectivity rate η. Proportions near to zero for Scenarios 2 and 4 and small values of η are explained by the fact that level l(3) is likely not reached in these situations. On the other hand, for increasing values of the infectivity rate η, the proportion of expected time spent above level l(3) directly depends on how many susceptible residents stay in the nursing home (which allows for larger epidemic peaks beyond K = 3 to be reached, likely leading to larger intervals of time until the epidemic curve decreases-by recoveries of infective residents-and goes again below level l(3)). Similar comments explain the behavior of Scenarios 1 and 3, but, in these, the impact of parameter η is negligible since a varicella infective resident starts the outbreak. In Figure 4, the interest is in the ratio E[α (i,z,a),0 ]/(N · E[η (i,z,a),0 ]), which represents the area under the curve of varicella infectious residents normalized by the total amount of time units that N residents stay at the nursing home during the varicella outbreak. We note here that this descriptor can be seen as a measure of cost, if, for example, the nursing home incurs in some costs per infective resident and per unit of time they remain infected. The behavior of this descriptor is similar to that observed in previous plots, whence we do not repeat preceding observations. However, we should emphasize here that the normalized area under the number of infectives is always less than the 7% of the total cost during the outbreak -measured in terms of N · E[η (i,z,a),0 ]-, in such a way that the smallest and the highest values of the ratio E[α (i,z,a),0 ]/(N · E[η (i,z,a),0 ]) correspond to Scenarios 2 and 3, respectively, in our numerical experiments.

Conclusions
In this paper, we extended the theory of finite QBD processes, given by Gaver et al. [1] by studying a number of random descriptors of interest which are defined inspired in their application to epidemic models, namely first-passage times and hitting probabilities to higher levels, number of upcrossings, and sojourn times, as well as the area under the level trajectory. We noted that these descriptors can have specific meanings when focusing on epidemic models: first-passage times and hitting probabilities can, for example, relate to the duration of an outbreak and the size of the epidemic, respectively; the number of upcrossings to a particular high level can represent the number of instances a particular threshold number of infectives is reached (which might lead to the introduction of particular control measures), while the area under the level trajectory can be translated into the area under the curve of infectives, a descriptor that has been analyzed before in the theory of mathematical epidemiology due to its potential interpretation in terms of the cost of a given outbreak. Although our results have been derived under the assumption that the finite QBD process is irreducible, they also hold under the less restrictive condition that l(0) is accessible from the class ∪ N i=1 l(i) of communicating states. In order to analyze these stochastic descriptors or summary statistics, we considered auxiliary bi-dimensional absorbing processes in which infinitesimal generators are described in terms of levels and phases, and we presented a block-structured form of these generators. We were able to exploit this structure to determine recursive relations involving sub-matrices of these infinitesimal generators, in order to compute the quantities of interest in an iterative and efficient way. Our algorithmic solution involves recursive procedures and block-Gaussian elimination, whence the computational complexity of Algorithms A1-A6 (Appendices A-D) can be written in a similar manner to the computational complexity of the linear level reduction algorithm of Gaver et al. [1], and Algorithms 1.A-3.A of Gómez-Corral and López-García [29]. This makes the memory requirements of Algorithms A1-A6, especially for storing auxiliary matrices in Algorithm A1, very demanding even for moderate values of N depending of the number of phases per level.
To illustrate the applicability of both our analytical and computational results, we have presented a numerical study of an epidemic model for the transmission of varicella-zoster virus within a nursing home. There is clearly future work to be done on the joint distribution of the random vector (I max , τ max ), where I max is the maximum number of residents who are simultaneously varicella infectious (i.e., the maximum of the level variable) and τ max is the time taken to reach this maximum number during an outbreak, and on how to relate herd immunity to the phase variable in the study of vaccination strategies in the context of varicella-zoster virus infections. Among other interesting problems to be addressed in epidemic models are also the extension to a population with variable size or the assumption of non-exponential events. These would lead to QBD processes with infinitely many possible values of the level variable and/or the phase variable, and piecewise-deterministic Markov processes, respectively. In these frameworks, the study of the stochastic descriptors analyzed in Sections 2-4 should require an analytical treatment and related numerical procedures that will be substantially different from those used in this paper; in particular, the analysis of first-passage times and sojourn times could benefit from the approach described by Dolgopyat and Goldsheid [45] , who obtain necessary and sufficient conditions for the existence of the density of the invariant measure for random walks when the environment is ergodic in both the transient and recurrent regimes. Another aspect that deserves further exploration is how the finite-dimensional linear algebra ideas in Sections 2 and 3 could be replaced by some suitable extension to the QBD case of the Karlin-McGregor orthogonal polynomial/spectral representation by translating the arguments in Reference [46][47][48] for discrete time processes to continuous time.