1. Introduction
The risks of an epidemic are difficult to assess and to guarantee. Nowadays, some insurance companies offer financial protection against such exposures. However, the use of classical epidemic models in insurance is still in its infancy.
There is a vast literature on mathematical models for the spread of infectious diseases. Much of the theory is taken up in the books of Diekmann and Hesterbeek [Reference Diekmann and Heesterbeek8] and Diekmann et al. [Reference Diekmann, Heesterbeek and Britton9] for the deterministic approach and Daley and Gani [Reference Daley and Gani7] and Andersson and Britton [Reference Andersson and Britton1] for the stochastic approach.
Susceptible–infective–removed (SIR) models form an important class of epidemic models. In the simplest formulation of a SIR-type model, a closed population is subdivided into three groups of individuals: the susceptibles (S), who are in good health but are susceptible to contracting the disease; the infectives (I), who have contracted the disease and are able to transmit it; and the removed cases (R), who have permanently left the infected status through recovery or death. The only possible transitions between groups are S 
$\to$ I (infection of suceptibles) and I 
$\to$ R (removal of infectives). Let S(t), I(t), and R(t) be the numbers of susceptibles, infectives, and removals at time 
$t\geq 0$, with the respective initial values 
$S(0)=n$, 
$I(0)=m$, 
$R(0)=0$. The population being closed, 
$R(t)=n+m-S(t)-I(t)$, so that only the process 
$\{S(t), I(t)\}$ may be followed. The epidemic terminates at the first moment 
$T_e$ when there are no more infectives present in the population, i.e. such that 
$I(T_e)=0$, which occurs almost surely in finite time. We will focus specifically on epidemic models of this type.
For insurance coverage of epidemics, the pioneering work is by Feng and Garrido [Reference Feng and Garrido10], who deal with the classical SIR model, called a general epidemic, in its deterministic version. Lefèvre et al. [Reference Lefèvre, Picard and Simon18] pursued that study by examining the Markovian version of the general epidemic model. The central idea is that the insurance receives a premium from each susceptible at rate 
$a>0$ and reimburses medical expenses to each infective at rate 
$b>0$. The objective there is to set the premium level a by applying the well-established equivalence principle in life insurance.
In this paper, we consider a similar but much broader framework, and this time, we are interested in the risk of ruin for the insurance company. For example, for the general epidemic model mentioned before, we write the associated reserves process 
$\{F(t)\}$ as
\begin{equation}F(t) = u + P(t) - E(t), \quad t\geq 0,\end{equation}
where u is the initial amount of reserves, and P(t) and E(t) respectively denote the income received and expenses incurred by the insurance up to time t. These can be represented by the area under functionals of the epidemic trajectory, i.e.
\begin{equation}P(t) = \int_0^t a(S(y),I(y)) \, dy, \qquad E(t) = \int_0^t b(S(y),I(y)) \, dy,\end{equation}
where 
$a(.)$ and 
$b(.)$ are appropriate positive functions. A simple case for P(t) and E(t) is when
\begin{equation*}P(t) = a \int_0^t S(y)\, dy, \qquad E(t) = b\int_0^t I(y)\,dy,\end{equation*}
with 
$a,b>0$, which is the assumption made in Lefèvre et al. [Reference Lefèvre, Picard and Simon18].
The company is ruined at the first moment 
$T_r$ when the reserves become zero, i.e. such that 
$F (T_r) = 0$, insofar as this occurs, which requires that 
$T_r < T_e$. Our objective is to discuss various issues related to this risk of ruin. Such a subject does not seem to have been studied so far. It is worth comparing (1.1) with the traditional reserves process in insurance (see e.g. the book of Asmussen and Albrecher [Reference Asmussen and Albrecher3]). In both cases, the reserves vary according to the arrival of premiums and claims. The crucial difference is that these quantities are not modeled here exogenously, but come from the very evolution of the epidemic.
The SIR epidemic class considered covers and generalizes a variety of standard models, such as general and fatal epidemics. It was introduced in Lefèvre and Simon [Reference Lefèvre and Simon19] as a block-structured Markov process. The insurance premiums and costs are not necessarily constant per individual, but are provided by arbitrary functions. We will mainly obtain the distribution of the time of ruin, the probability of ruin, and the final amount of reserves. To this end, the joint evolution of the epidemic and reserves processes will be described as a Markov-modulated fluid flow model, which then allows us to work using matrix-analytic methods (see e.g. the books of Latouche and Ramaswami [Reference Latouche and Ramaswami17] and He [Reference He13]).
Fluid flow models have been studied by many researchers, including Asmussen [Reference Asmussen2], Ramaswami [Reference Ramaswami, Smith, Hey and Elsevier21], Bean et al. [Reference Bean, O’Reilly and Taylor5], and Latouche and Nguyen [Reference Latouche and Nguyen16], amongst others. They were first introduced in telecommunications and operations research, but have been used later in various other fields of application. This is the case in the theory of risk and ruin in the actuarial sciences. Much can be found in the review of Badescu and Landriault [Reference Badescu and Landriault4] and the references therein. The use of fluid flows in the context of epidemic insurance is an original approach to our knowledge.
The paper is organized as follows. In Section 2, we present the generalized SIR epidemic model considered and the fluid model describing the insurance risk process. In Section 3, we derive the joint distribution of the time of ruin and the associated epidemic state, with the probability of ruin as a particular case. In Section 4, we determine the distribution of the amount of reserves when the epidemic dies out, first regardless of the occurrence of a possible ruin, then when ruin occurs (or not) before the end of the epidemic. The expected amounts of the corresponding final reserves are also deduced. In Section 5, we add a Brownian noise to the fluid model for the reserves and show that the distributions of the ruin time and the final reserves can be obtained in rather similar ways. In Section 6, we illustrate and comment on some of our results by examining two particular standard SIR models: the general epidemic and the fatal epidemic.
2. Epidemics and insurance
We begin by defining the extended SIR epidemic model constructed by Lefèvre and Simon [Reference Lefèvre and Simon19] (Section 2.1). Next, we introduce a fluid model to describe the evolution of the reserves of the insurance company (Section 2.2). Finally, we present the matrix framework adapted to the proposed approach (Section 2.3).
2.1. A class of SIR-type epidemics
We are concerned with SIR epidemic models whose propagation is described by a bivariate Markov process 
$\{Y(t)\}=\{S(t), \varphi(t)\}$.
The variable S(t) counts the susceptibles present in the population at time t. Initially, 
$S(0)=n$. If ever infected, a susceptible becomes an infective before being finally removed. Thus, S(t) decreases with time and the susceptible state space is 
$\{n,n-1, \ldots, 0\}$.
The variable 
$\varphi(t)$ represents the state of the infection process at time t. In classical models such as the general epidemic, 
$\varphi(t)$ is simply the number I(t) of infectives in the population. However, 
$\varphi(t)$ may be more general. It can be multivariate, of the form 
$\varphi(t)=(I_1(t), \ldots, I_v(t))$, to represent the number of infectives in several groups which differ in their behaviors or infectivity levels. Individual transitions from one group to another are then possible. For example, an individual could enter the 
$I_1$ group when he has been exposed to the infection but is not yet contagious, and then move into the 
$I_2$ group when he becomes contagious (as in a so-called SEIR model). The infection process can also be of the form 
$\varphi(t)=(I(t), E(t))$ where 
$\{E(t)\}$ is an external Markov environment, climatic or economic, which affects the model parameters.
The states of 
$\varphi(t)$ which can be reached at a given time depend, of course, on the number of susceptibles present. When there are s susceptibles, the infection state space is a finite set 
${\mathcal{E}}_s$ of transient states and an absorbing state denoted by 
$\{0\}$. The epidemic process ends in this state 
$\{0\}$ when all sources of infection are extinguished. This will occur at the (almost surely finite) time 
$T_e$; i.e.
\begin{equation}T_e=\inf\{t>0 \,|\, \varphi(t)=0\}.\end{equation}
The initial infection state 
$\varphi(0)$ is often difficult to evaluate. It is therefore assumed to be random, with probability mass function given by the row vector 
$\boldsymbol{\pi}$ of components
\begin{equation}\pi_i = {\mathbb{P}\left( {\varphi(0) = i} \right)}, \quad i \in {\mathcal{E}}_n.\end{equation}
Returning to the bivariate Markov process 
$\{Y(t)\}$, its state space is 
$\{0\} \cup V_n \cup V_{n-1} \cup \cdots \cup V_0$, where 
$\{0\}$ is still the absorbing state, and each set 
$V_s = \{(s,i), \, i \in {\mathcal{E}}_s\}$ contains the possible transient states when there are s susceptibles. The associated generator is defined by
\begin{equation}\begin{bmatrix} 0 &\quad \vline &\quad &\quad \textbf{0}^{\tau} &\quad &\quad \\\hline& \vline & & & & \\\, \textbi{v} \!\!\! & \vline & & Q & & \\& \vline & & & &\end{bmatrix},\end{equation}
where 
$\textbf{0}^{\tau}$, the transpose of 
$\textbf{0}$, is a row vector with all entries equal to 0, and
\begin{equation}\textbi{v} =\begin{bmatrix}\textbi{v}(n) \\[3pt] \textbi{v}(n\!-\!1)\\[3pt] \textbi{v}(n\!-\!2) \\[3pt] \textbi{v}(n\!-\!3)\\[3pt] \vdots \\[3pt] \textbi{v}(1) \\[3pt] \textbi{v}(0)\end{bmatrix},~ Q =\begin{bmatrix} A(n) &\quad B(n) &\quad 0 &\quad 0 &\quad \cdots &\quad 0 &\quad 0 \\[3pt]0&\quad A(n\!-\!1) &\quad B(n\!-\!1) &\quad 0 &\quad \cdots &\quad 0 &\quad 0 \\[3pt]0&\quad0&\quad A(n\!-\!2) &\quad B(n\!-\!2) &\quad \cdots &\quad 0 &\quad 0 \\[3pt]0&\quad0&\quad0&\quad A(n\!-\!3) &\quad \cdots &\quad 0 &\quad 0 \\[3pt] \vdots &\quad\vdots &\quad\vdots &\quad \vdots &\quad \ddots &\quad \vdots &\quad \vdots \\[3pt]0&\quad0&\quad0&\quad 0 &\quad \cdots &\quad A(1) &\quad B(1) \\[3pt]0&\quad0&\quad0&\quad 0 &\quad \cdots &\quad 0 &\quad A(0)\end{bmatrix},\end{equation}
in which, for each 
$s=0, \ldots, n$,
•
$\textbi{v}(s)$ is a column vector of dimension $|{\mathcal{E}}_s|$ such that $v_i(s)$, $i \in {\mathcal{E}}_s$, is the rate of absorption in $\{0\}$ from the state (s,i);• B(s) is a rectangular matrix of dimensions
$|{\mathcal{E}}_s|\times |{\mathcal{E}}_{s-1}|$ such that $B_{i,j}(s)$ is the rate at which one of the s suceptibles is infected with a change of infectious state from $i \in {\mathcal{E}}_s$ to $j\in {\mathcal{E}}_{s-1}$;• A(s) is a square matrix of dimension
$|{\mathcal{E}}_s|$ such that $A_{i,j}(s)$ for $i\neq j$ is the rate at which the infection jumps from $i \in {\mathcal{E}}_s$ to $j \in {\mathcal{E}}_s$ without any change in the number s of susceptibles, and $A_{i,i}(s)$, $i \in {\mathcal{E}}_s$, is defined to satisfy the relations
(2.5)where
\begin{equation}\textbi{v}(0) + A(0) \textbf{1} = \textbf{0} \,\,\, \mbox{ and } \,\,\, \textbi{v}(s) + A(s) \textbf{1}+ B(s) \textbf{1} = \textbf{0}, \,\,\,\, 1\le s \le n,\end{equation}
$\textbf{1}$ is a column vector with all entries equal to 1.
We recall that in the traditional general epidemic model, 
$\varphi(t)$ denotes the number of infectives, each remaining infectious during an exponential period and able to infect any susceptible according to a Poisson process. The current epidemic model generalizes the general epidemic insofar as (i) 
$\varphi(t)$ counts the infectives but may also incorporate other infection sources, (ii) all the transitions between phases are possible and the sojourn distributions are of phase-type, and (iii) the infection process is general and an infection can arbitrarily affect the current phase. The flexibility of the modeling considered here was underlined by Lefèvre and Simon [Reference Lefèvre and Simon19] by examining four extensions of the general epidemic, namely with arbitrary contagion rates, a detection process of infectious agents, the influence of a random environment, and interference between disease strains.
2.2. Fluid model for the risk process
We now propose a simple way to insert an insurance risk process as part of an epidemic process that is spreading within a company.
The insurance plan. The goal of the insurance is to collect premiums from people in good health and provide care coverage to those infected. On the one hand, the insurance receives premiums at a continuous rate that depends on the state occupied by 
$\{Y(t)\}$. Specifically, when 
$Y(t)=(s,i)$, the amount of the premium collected is assumed to be a(s, i) per time unit. For example, 
$a(s,i) = as$ means that the premium is charged at the rate a for the susceptibles; this is the simple case presented in the introduction. If the removed individuals also contribute to the premium at a continuous rate c, then 
$a(s,i) = as + c(n+m-s-i)$. The function a(s, i) can be nonlinear to take into account, for example, group effects in premium collection.
On the other hand, the insurance reimburses medical expenses at a continuous rate which also depends on the state occupied by 
$\{Y(t)\}$. Specifically, when 
$Y(t)=(s,i)$, the amount reimbursed is assumed to be b(s, i) per time unit. For instance, 
$b(s,i) = bi$ means that coverage is provided at the rate b for the infectives, as considered in the introduction. If the company also pays a logistic cost at rate l when the number of infectives exceeds a fixed threshold 
$i^{\star}$, then 
$b(s,i) = bi + l1_{(i>i^{\star})}$. The function b(s,i) can depend on s in some special cases; for example, the level of coverage could become limited if there are not enough healthy people.
Therefore, the difference between the two rates,
\begin{equation}c(s,i) = a(s,i) - b(s,i),\end{equation}
represents the net instantaneous income rate of the insurance company when 
$Y(t)=(s,i)$. It may be positive or negative, depending on the values of s and i.
The reserves process. To the epidemic model 
$\{Y(t)\}$, here called the phase process, we associate a new process 
$\{F(t)\}$, called fluid flow, where the variable F(t) represents the reserves of the insurance at time t. Initially 
$F(0)=u\geq 0$, the starting capital of the company. The process 
$\{F (t)\}$ varies according to 
$\{Y(t)\} $ through the function of income rate c(Y(t)). More precisely,
\begin{equation}\frac{d}{dt}F(t) = \left\{\begin{array}{l@{\quad}l}c(s,i) &\textrm{if } t
< T_e\, \mbox{ and } \,Y(t)=(s,i),\\[3pt] 0 &\textrm{if } t \ge T_e.\end{array}\right.\end{equation}
Note that F(t) takes values on the real line and can increase or decrease over time. Its trajectories are continuous and piecewise linear.
We can now look at the joint process 
$\{F(t), Y(t)\}$. It corresponds to a Markov-modulated fluid flow as studied by Bean et al. [Reference Bean, O’Reilly and Taylor5] and others mentioned in the introduction. The process is characterized by the previous sub-generator Q, defined in (2.3), (2.4), and (2.5), and the diagonal matrix C of the income rates (2.6), given by
\begin{equation}C = \begin{bmatrix} C(n) &\quad 0 &\quad \cdots &\quad 0 \\[3pt]0&\quad C(n\!-\!1) &\quad \cdots &\quad 0 \\[3pt]\vdots &\quad\vdots &\quad \ddots &\quad \vdots \\[3pt]0&\quad0 &\quad \cdots &\quad C(0)\end{bmatrix},\end{equation}
where C(s), 
$s=0,\ldots,n$, is the diagonal matrix of dimension 
$|{\mathcal{E}}_s|$ such that 
$C_{i,i}(s) = c(s,i)$, 
$i \in {\mathcal{E}}_s$.
It is clear that the reserves F(t) can become negative before the ending time 
$T_e$ of the epidemic (see (2.1)), thus leading to the ruin of the insurance company. We denote by 
$T_r$ the time of the ruin, if it actually occurs, i.e.
\begin{equation}T_r=\inf\{t>0 \,|\, F(t)<0\}.\end{equation}
A main objective of our work is to discuss the two situations in which the insurance company is ruined or not during the epidemic.
2.3. Preliminaries
It is obvious that ruin can only occur when the infection process 
$\{\varphi (t)\}$ is in a state (s, i) associated with a negative rate c (s, i). Consequently, it is appropriate in what follows to separate the states (s, i) in which the reserves F (t) increase and those in which they decrease. For each s, we thus define 
${\mathcal{E}}_s = {\mathcal{E}}^+_s \cup {\mathcal{E}}^-_s$, where
\begin{equation*}{\mathcal{E}}^+_s = \{i \in {\mathcal{E}}_s \,|\, c(s,i)>0\} \,\,\,\mbox{ and } \,\,\, {\mathcal{E}}^-_s = \{i \in {\mathcal{E}}_s \,|\, c(s,i)<0\}.\end{equation*}
To simplify, we assume that c(s,i) is never zero, but this assumption can be lifted without difficulty. By adopting this subdivision, we can rewrite the row vector 
$\boldsymbol{\pi}$ in (2.2) as 
$[\boldsymbol{\pi}_+,\boldsymbol{\pi}_-]$ and partition the matrices A(s), B(s), C(s) in (2.4), (2.5) into the form
\begin{equation*}A(s) = \begin{bmatrix}A_{++}(s) &\quad A_{+-}(s) \\[3pt]A_{-+}(s) &\quad A_{--}(s)\end{bmatrix}, ~~B(s) = \begin{bmatrix}B_{++}(s) &\quad B_{+-}(s) \\[3pt]B_{-+}(s) &\quad B_{--}(s)\end{bmatrix}, ~~C(s) = \begin{bmatrix}C_{+}(s) &\quad 0 &\quad \\[3pt]0 &\quad C_{-}(s)\end{bmatrix}.\end{equation*}
Next, defining
\begin{equation*}V^+_s = \{(s,i) \,|\, i \in {\mathcal{E}}^+_s\} \,\,\,\mbox{ and } \,\,\, V^-_s = \{(s,i) \,|\, i \in {\mathcal{E}}^-_s\},\end{equation*}
we reorganize the matrices Q, C according to the order of states 
$(V^+_n \cup V^+_{n-1} \cup \cdots \cup V^+_0) \cup (V^-_n \cup V^-_{n-1} \cup \cdots \cup V^-_0)$. This leads to the partition
\begin{equation}Q = \begin{bmatrix}Q_{++} &\quad Q_{+-} \\[3pt]Q_{-+} &\quad Q_{--}\end{bmatrix}, ~~~~C = \begin{bmatrix}C_{+} &\quad 0 \\[3pt]0 &\quad C_{-}\end{bmatrix},\end{equation}
where
\begin{align}Q_{++} = \begin{bmatrix}A_{++}(n) &\quad B_{++}(n) &\quad 0 &\quad 0 &\quad \cdots &\quad 0 &\quad 0 \\[3pt]0&\quad A_{++}(n\!-\!1) &\quad B_{++}(n\!-\!1) &\quad 0 &\quad \cdots &\quad 0 &\quad 0 \\[3pt]0&\quad0&\quad A_{++}(n\!-\!2) &\quad B_{++}(n\!-\!2) &\quad \cdots &\quad 0 &\quad 0 \\[3pt]0&\quad0&\quad0&\quad A_{++}(n\!-\!3) &\quad \cdots &\quad 0 &\quad 0 \\[3pt]\vdots &\quad\vdots &\quad\vdots &\quad \vdots &\quad \ddots &\quad \vdots &\quad \vdots \\[3pt]0&\quad0&\quad0&\quad 0 &\quad \cdots &\quad A_{++}(1) &\quad B_{++}(1) \\[3pt]0&\quad0&\quad0&\quad 0 &\quad \cdots &\quad 0 &\quad A_{++}(0)\end{bmatrix};\end{align}
$Q_{+-}, Q_{-+}, Q_{--}$ are defined in the same way but with the appropriate indices; and
\begin{equation}C_{+} = \begin{bmatrix} C_{+}(n) &\quad 0&\quad \cdots &\quad 0 \\[3pt]0&\quad C_{+}(n\!-\!1) &\quad \cdots &\quad 0 \\[3pt]\vdots &\quad\vdots &\quad \ddots &\quad \vdots \\[3pt]0&\quad0 &\quad \cdots &\quad C_{+}(0)\end{bmatrix},\end{equation}
with 
$C_{-}$ again built similarly. It is with this matrix structure (2.10), (2.11), (2.12) that we will work.
3. Ruin time and probability of ruin
This section is concerned with the case where 
$T_r<T_e$. Our goal is to obtain the joint distribution of the time of ruin 
$T_r$ and the corresponding epidemic state 
$Y(T_r)$. As usual, 
${\mathds{1}_\mathcal{A}}$ denotes the indicator variable of a random event 
$\mathcal{A}$. We will determine the following Laplace transform:
\begin{equation*}{\mathbb{E}\left[ {e^{-\theta T_r} \,{\mathds{1}_{T_r
< T_e,\, Y(T_r)=(s,j)}}} \right]},\end{equation*}
where 
$\theta \geq 0$, 
$0\leq s \leq n$, and 
$j\in {\mathcal{E}}^-_s$ (by the very construction of the model).
To begin, we construct two first-passage matrices 
$\Psi_{\theta}$ and 
$\Phi_{\theta}(x)$, 
$x \ge 0$, which will play a central role in the analysis. For 
$0 \le s \le v \le n$, 
$i\in {\mathcal{E}}^+_v$, 
$j\in {\mathcal{E}}^-_s$, we define
\begin{equation*}\left(\Psi_{\theta}\right)_{(v,i),(s,j)} \equiv \left(\Psi_{\theta}(v,s)\right)_{i,j} = {\mathbb{E}\left[ {e^{-\theta T_r} \,{\mathds{1}_{T_r
< T_e, Y(T_r)=(s,j)}} \,|\,F(0)=0, Y(0)=(v,i)} \right]},\end{equation*}
and for 
$0 \le s \le v \le n$, 
$i\in {\mathcal{E}}^-_v$, 
$j\in {\mathcal{E}}^-_s$, we define
\begin{align*}\left(\Phi_{\theta}\right)_{(v,i),(s,j)}(x) &\equiv \left(\Phi_{\theta}(x;v,s)\right)_{i,j}\\[3pt]&= {\mathbb{E}\left[ {e^{-\theta T_r} \,{\mathds{1}_{T_r
< T_e, Y(T_r)=(s,j)}} \,|\,F(0)=x, Y(0)=(v,i)} \right]}.\end{align*}
These expectations give the Laplace transforms of the ruin time 
$T_r$ with the state 
$Y(T_r)$ when ruin arises before the end time 
$T_e$, under the condition that the reserves process starts for 
$\Psi_{\theta}$ at the level 0 in an ascending phase (
$i\in {\mathcal{E}}^+_v$), and for 
$\Phi_{\theta}(x)$ at the level x in a descending phase (
$i\in {\mathcal{E}}^-_v$).
As the number of susceptibles can only decrease, the matrices 
$\Psi_{\theta}$ and 
$\Phi_{\theta}(x)$ have an upper triangular block structure of the form
\begin{equation}\Psi_{\theta} = \begin{bmatrix}\Psi_{\theta}(n,n) &\quad \Psi_{\theta}(n,n\!-\!1) &\quad \Psi_{\theta}(n,n\!-\!2) &\quad \cdots &\quad \Psi_{\theta}(n,0) \\[3pt]0&\quad \Psi_{\theta}(n\!-\!1,n\!-\!1) &\quad \Psi_{\theta}(n\!-\!1,n\!-\!2) &\quad \cdots &\quad \Psi_{\theta}(n\!-\!1,0)\\[3pt]0&\quad0&\quad \Psi_{\theta}(n\!-\!2,n\!-\!2) &\quad \cdots &\quad \Psi_{\theta}(n\!-\!2,0)\\[3pt]\vdots &\quad\vdots &\quad \vdots &\quad \ddots &\quad \vdots \\[3pt]0&\quad0&\quad0&\quad \cdots &\quad \Psi_{\theta}(0,0)\end{bmatrix},\end{equation}
\begin{equation}\Phi_{\theta}(x) = \begin{bmatrix}\Phi_{\theta}(x;n,n) &\quad \Phi_{\theta}(x;n,n\!-\!1) &\quad \Phi_{\theta}(x;n,n\!-\!2) &\quad \cdots &\quad \Phi_{\theta}(x;n,0) \\[3pt]0&\quad \Phi_{\theta}(x;n\!-\!1,n\!-\!1) &\quad \Phi_{\theta}(x;n\!-\!1,n\!-\!2) &\quad \cdots &\quad \Phi_{\theta}(x;n\!-\!1,0)\\[3pt]0&\quad0&\quad \Phi_{\theta}(x;n\!-\!2,n\!-\!2) &\quad \cdots &\quad \Phi_{\theta}(x;n\!-\!2,0)\\[3pt]\vdots &\quad\vdots &\quad \vdots &\quad \ddots &\quad \vdots \\[3pt]0&\quad0&\quad0&\quad \cdots &\quad \Phi_{\theta}(x;0,0)\end{bmatrix}.\end{equation}
Note that a subspace 
${\mathcal{E}}^+_v$ may be empty, so that the corresponding matrix 
$\Psi_{\theta}(v,s)$ is not defined. Similarly, if 
${\mathcal{E}}^-_v$ is empty, the matrix 
$\Phi_{\theta}(x;v,s)$ does not exist. This also occurs with the matrix 
$U_\theta(v,s)$ in (3.5) below and all the matrices defined in Section 4. In these cases, the parts of the equations that contain them will be ignored.
The Laplace transform mentioned above can be expressed from the block matrices of 
$\Phi_{\theta}(u)$ and of the first row of 
$\Psi_{\theta}$. The notation 
$\{\ldots\}_j$ in (3.3) below means the jth component of the row vector 
$\{\ldots\}$.
Proposition 3.1. For 
$0\leq r \leq n$, 
$j \in {\mathcal{E}}_{n-r}^-$,
\begin{eqnarray}&&{\mathbb{E}\left[ {e^{-\theta T_r} \,{\mathds{1}_{T_r
< T_e,\, Y(T_r)=(n-r,j)}}} \right]}\nonumber\\&&\qquad\quad = \left\{\boldsymbol{\pi}_- \Phi_{\theta}(u;n,n-r) + \boldsymbol{\pi}_+ \sum_{l=0}^{r} \Psi_{\theta}(n,n-l) \Phi_{\theta}(u;n-l,n-r)\right\}_j.\end{eqnarray}
Proof. The result is easily derived from the definition of 
$\Psi_{\theta}$ and 
$\Phi_{\theta}(x)$. With the probabilities in 
$\boldsymbol{\pi}_-$, the fluid process F(t) starts from the level u in a descending phase, and the matrix 
$\Phi_{\theta}(u;n,n-r)$ gives the Laplace transform of 
$T_r$ with 
$S(T_r)=n-r$ when 
$T_r < T_e$.
With the probabilities in 
$\boldsymbol{\pi}_+$, the process F(t) starts from the level u in an ascending phase. Thus, it must first return to the level u before 
$T_e$, and the susceptibles present at that time will number, say, 
$n-l$. The process F(t) is then at the level u in a descending phase and the previous argument is again applicable.
To apply Proposition 3.1, we must first determine the blocks in the matrices of 
$\Psi_{\theta}$ and 
$\Phi_{\theta}(x)$. From the definition of 
$\Phi_{\theta}(x)$, we notice (see e.g. Ramaswami [Reference Ramaswami, Smith, Hey and Elsevier21]) that 
$\Phi_{\theta}(x)$ can be expressed in the exponential form
\begin{equation}\Phi_{\theta}(x) = e^{U_{\theta} x},\end{equation}
where 
$U_{\theta}$ is a sub-generator with the same block structure as 
$\Psi_{\theta}$, i.e.
\begin{equation}U_{\theta} = \begin{bmatrix}U_{\theta}(n,n) &\quad U_{\theta}(n,n\!-\!1) &\quad U_{\theta}(n,n\!-\!2) &\quad \cdots &\quad U_{\theta}(n,0) \\[3pt]0&\quad U_{\theta}(n\!-\!1,n\!-\!1) &\quad U_{\theta}(n\!-\!1,n\!-\!2) &\quad \cdots &\quad U_{\theta}(n\!-\!1,0)\\[3pt]0&\quad0&\quad U_{\theta}(n\!-\!2,n\!-\!2) &\quad \cdots &\quad U_{\theta}(n\!-\!2,0)\\[3pt]\vdots &\quad\vdots &\quad \vdots &\quad \ddots &\quad \vdots \\[3pt]0&\quad0&\quad0&\quad \cdots &\quad U_{\theta}(0,0)\end{bmatrix}.\end{equation}
We now show that the blocks of the matrices 
$\Psi_{\theta}$ and 
$U_{\theta}$ can be obtained recursively. Let 
$1_\mathcal{B}$ be the indicator function of a non-random set 
$\mathcal{B}$, and denote by 
$|C_-(s)|$ the entrywise absolute value of 
$C_-(s)$.
Proposition 3.2. For 
$0 \le k \le s \le n$, the matrices 
$\Psi(s,s-k)$ and 
$U(s,s-k)$ satisfy the identities
\begin{align}& (C_+(s))^{-1}[A_{++}(s)-\theta I] \Psi_{\theta}(s,s-k) + (C_+(s))^{-1}B_{++}(s) \Psi_{\theta}(s-1,s-k)\, 1_{s,k>0} \nonumber \\[3pt]& \qquad \qquad + \sum_{l=0}^{k} \Psi_{\theta}(s,s-l) U_{\theta}(s-l,s-k) \nonumber\\[3pt]& \qquad \qquad = -(C_+(s))^{-1}A_{+-}(s)\, 1_{k=0} - (C_+(s))^{-1}B_{+-}(s)\, 1_{s>0,k=1}\end{align}
and
\begin{align}U_{\theta}(s,s-k) = &|C_-(s)|^{-1} [A_{--}(s)-\theta I] \, 1_{k=0} + |C_-(s)|^{-1} B_{--}(s) \, 1_{s>0,k=1} \nonumber \\[3pt]&+ |C_-(s)|^{-1} A_{-+}(s)\Psi_{\theta}(s,s-k) \nonumber\\[3pt]&+ |C_-(s)|^{-1} B_{-+}(s)\Psi_{\theta}(s-1,s-k) \, 1_{s,k>0}.\end{align}
Proof. To derive (3.6), we start by expressing 
$\Psi_{\theta}(s, s-k)$ with respect to the level of the reserves at time 
$\tau=\min(\tau_1, \tau_2)$, where 
$\tau_1 $ is the end of the first ascent and 
$\tau_2$ the moment of the first contamination. This gives
\begin{align}&\Psi_{\theta}(s, s-k) = \int_0^\infty e^{(C_+(s))^{-1}[A_{++}(s)-\theta I]y} (C_+(s))^{-1}A_{+-}(s) \Phi_{\theta}(y;s,s-k) dy \nonumber\\[3pt]&\,\,\, + 1_{s>0}\int_0^\infty e^{(C_+(s))^{-1}[A_{++}(s)-\theta I]y} (C_+(s))^{-1}B_{+-}(s) \Phi_{\theta}(y;s-1,s-k) dy \nonumber \\[3pt]&\,\,\, + 1_{s>0}\int_0^\infty e^{(C_+(s))^{-1}[A_{++}(s)-\theta I]y} (C_+(s))^{-1}B_{++}(s)\nonumber \\[3pt]& \qquad \qquad \cdot \sum_{l=1}^{k} \Psi (s-1, s-l) \Phi_{\theta}(y;s-l,s-k) dy.\end{align}
Let us multiply both sides of (3.8) by 
$(C_+(s))^{-1}[A_{++}(s)-\theta I]$ and integrate by parts the three integrals in the right-hand side. Note from (3.4) that 
$\Phi_{\theta}'(y) = \Phi_{\theta}(y) U_{\theta}$, so that using the block structures (3.2) and (3.5),
\begin{equation}\Phi_{\theta}'(y;s-l,s-k) = \sum_{v=l}^k \Phi_{\theta}(y;s-l,s-v) U_{\theta}(s-v,s-k).\end{equation}
From (3.8) and using (3.9), we then obtain
\begin{align}& (C_+(s))^{-1}[A_{++}(s)-\theta I]\Psi_{\theta}(s, s-k) = -(C_+(s))^{-1}A_{+-}(s) \Phi_{\theta}(0;s,s-k)\nonumber\\[3pt]& \qquad\qquad\qquad - 1_{s>0}\, (C_+(s))^{-1}B_{+-}(s) \Phi_{\theta}(0;s-1,s-k)\nonumber\\[3pt]& \qquad\qquad\qquad - 1_{s>0}\, (C_+(s))^{-1}B_{++}(s) \sum_{l=1}^{k} \Psi_{\theta}(s-1, s-l) \Phi_{\theta}(0;s-l,s-k)\nonumber\\[3pt]& \qquad\qquad\qquad - \sum_{v=0}^{k} \Psi_{\theta}(s, s-v) U_{\theta}(s-v,s-k).\end{align}
Since 
$\Phi_{\theta}(0;s-l,s-k)=I \, 1_{k=l}$, (3.10) gives the desired relations (3.6).
To obtain (3.7), we express 
$\Phi_{\theta}(x;s,s-k)$ with respect to the level of the reserves at time 
$\tau^*=\min(\tau_1^*, \tau_2)$, where 
$\tau_1^*$ is now the end of the first descent and 
$\tau_2$ is still the moment of the first contamination. This yields
\begin{equation}\Phi_{\theta}(x;s,s-k) = e^{|C_-(s)|^{-1}[A_{--}(s)-\theta I]x} \, [I\, 1_{k=0} + F(x)],\end{equation}
where
\begin{align}F(x) =& 1_{s>0} \int_0^x e^{-|C_-(s)|^{-1}[A_{--}(s)-\theta I]y} |C_-(s)|^{-1} B_{--}(s) \Phi_{\theta}(y;s-1,s-k) dy \nonumber \\[3pt]&+ \int_0^x e^{-|C_-(s)|^{-1}[A_{--}(s)-\theta I]y} |C_-(s)|^{-1} A_{-+}(s)\nonumber \\[3pt]&\qquad \cdot \sum_{l=0}^{k} \Psi_{\theta}(s, s-l) \Phi_{\theta}(y;s-l,s-k) dy \nonumber\\[3pt]&+ 1_{s>0} \int_0^x e^{-|C_-(s)|^{-1}[A_{--}(s)-\theta I]y} |C_-(s)|^{-1} B_{-+}(s)\nonumber\\[3pt]&\qquad \qquad \cdot \sum_{l=1}^{k} \Psi_{\theta}(s-1, s-l)\Phi_{\theta}(y;s-l,s-k) dy.\end{align}
Now, (3.4) implies that 
$\Phi_{\theta}'(0) = U_{\theta}$; hence 
$U_{\theta}(s,s-k)=\Phi_{\theta}'(0;s,s-k)$. Thus, differentiating (3.12) with 
$x=0$ gives
\begin{eqnarray}U_{\theta}(s,s-k) = |C_-(s)|^{-1}[A_{--}(s)-\theta I] [I \, 1_{k=0} + F(0)] + F'(0).\end{eqnarray}
By definition, 
$F(0)=0$ and we obtain from (3.11) that
\begin{align}F'(0) & = 1_{s>0} \, |C_-(s)|^{-1} B_{--}(s) \Phi_{\theta}(0;s-1,s-k) \nonumber\\[3pt]&\quad + |C_-(s)|^{-1} A_{-+}(s) \sum_{l=0}^{k}\Psi_{\theta}(s,s-l) \Phi_{\theta}(0;s-l,s-k) \nonumber\\[3pt]&\quad + 1_{s>0} \, |C_-(s)|^{-1} B_{-+}(s) \sum_{l=1}^{k}\Psi_{\theta}(s-1,s-l) \Phi_{\theta}(0;s-l,s-k).\end{align}
From (3.13) and (3.14), we deduce (3.7) using again that 
$\Phi_{\theta}(0;s-l,s-k)=I \, 1_{k=l}$.
The blocks inside 
$\Psi_{\theta}$ and 
$U_{\theta}$ can be computed recursively from the relations (3.6) and (3.7) as follows.
Induced algorithm.
(i) First, for each
$0\leq s \leq n$, we take $k=0$, so that (3.7) gives
(3.15)and after insertion in (3.6),
\begin{equation}U_{\theta}(s,s) = |C_-(s)|^{-1} [A_{--}(s)-\theta I] + |C_-(s)|^{-1} A_{-+}(s)\Psi_{\theta}(s,s),\end{equation}
$Z\equiv \Psi_{\theta}(s,s)$ satisfies the Riccati equation
(3.16)It can be proved that
\begin{eqnarray}&& (C_+(s))^{-1}[A_{++}(s)-\theta I]Z + Z|C_-(s)|^{-1} [A_{--}(s)-\theta I] + Z|C_-(s)|^{-1} A_{-+}(s)Z \nonumber \\&& \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad = -(C_+(s))^{-1}A_{+-}(s).\end{eqnarray}
$\Psi_{\theta}(s,s)$ is the minimal nonnegative solution of the equation (3.16) (see [Reference Bean, O’Reilly and Taylor5]). This equation can be solved numerically in an efficient way using the algorithms proposed e.g. by Asmussen [Reference Asmussen2], Bini et al. [Reference Bini, Iannazzo, Latouche and Meini6], or Guo [Reference Guo12]. Once $\Psi_{\theta}(s,s)$ is known, $U_{\theta}(s,s)$ is provided by (3.15).
-
(ii) Next, for each
$1\leq s \leq n$, we take $k=1$, so that from (3.7),
(3.17)and by injecting this into (3.6), we find that
\begin{equation}{U_\theta }(s,s - 1) = |C\_(s){|^{ - 1}}[{B_{ - - }}(s) + {A_{ - + }}(s){\Psi _\theta }(s,s - 1) + {B_{ - + }}(s){\Psi _\theta }(s - 1,s - 1)],\end{equation}
$Z\equiv \Psi_{\theta}(s,s-1)$ is the unique solution of the Sylvester equation
(3.18)where G(s) is defined as
\begin{eqnarray}&& G(s)Z + ZU_{\theta}(s-1,s-1) = - (C_+(s))^{-1} [B_{+-}(s) + B_{++}(s) \Psi_{\theta}(s-1,s-1)] \nonumber\\[3pt]&& \qquad - \Psi_{\theta}(s,s)|C_-(s)|^{-1} [B_{--}(s) + B_{-+}(s)\Psi_{\theta}(s-1,s-1)],\end{eqnarray}
(3.19)Remember that
\begin{equation}G(s) = (C_+(s))^{-1} [A_{++}(s)-\theta I] + \Psi_{\theta}(s,s) |C_-(s)|^{-1} A_{-+}(s).\end{equation}
$\Psi_{\theta}(s,s)$ and $U_{\theta}(s-1,s-1)$ were computed in Step (i). The equation (3.18) can be efficiently solved numerically (see e.g. Gardiner et al. [Reference Gardiner, Laub, Amato and Moler11]). When $\Psi_{\theta}(s,s-1)$ is obtained, $U_{\theta}(s,s-1)$ follows from (3.17).
-
(iii) The same procedure applies successively to each
$k=2, \ldots, s$. From (3.7), $U_{\theta}(s,s-k)$ is given by
(3.20)From (3.6) and (3.20),
\begin{equation}U_{\theta}(s,s-k) = |C_-(s)|^{-1} \left[A_{-+}(s)\Psi_{\theta}(s,s-k) + B_{-+}(s)\Psi_{\theta}(s-1,s-k)\right].\end{equation}
$Z\equiv \Psi_{\theta}(s,s-k)$ is the unique solution of the Sylvester equation
(3.21)where G(s) is defined in (3.19). After solving (3.21), we get
\begin{eqnarray}&& G(s)Z + ZU_{\theta}(s-k,s-k) = - \sum_{l=1}^{k-1} \Psi_{\theta}(s,s-l) U_{\theta}(s-l,s-k) \nonumber\\&& \qquad\quad - [(C_+(s))^{-1}B_{++}(s) + \Psi_{\theta}(s,s)|C_-(s)|^{-1}B_{-+}(s)] \Psi_{\theta}(s-1,s-k),\end{eqnarray}
$U_{\theta}(s,s-k)$ from (3.20).
Once the blocks of 
$U_{\theta}$ are known, the exponential 
$\Phi_{\theta}(x)$ is efficiently evaluated from (3.4) using the block triangular structure of 
$U_{\theta}$ (see e.g. Kressner et al. [Reference Kressner, Luce and Statti14]).
Ruin probability. Obviously, (3.3) with 
$\theta=0$ gives an expression for
\begin{equation*}{\mathbb{P}\left( {T_r
< T_e,\, Y(T_r)=(n-r,j)} \right)},\end{equation*}
the probability that ruin will occur and that there are 
$n-r$ susceptibles and j infectives present at that time in the population.
Differentiating (3.3) with respect to 
$\theta$ and taking 
$\theta = 0$ yields the moments of the ruin time when this arises. In particular,
\begin{align}&{\mathbb{E}\left[ {T_r\,{\mathds{1}_{T_r
< T_e,\, Y(\tau)=(n-r,j)}}} \right]} = \{\boldsymbol{\pi}_- \Phi'_0(u;n,n-r) \nonumber\\&\quad +\boldsymbol{\pi}_+ \sum_{l=0}^{r} [\Psi'_0(n,n-l) \Phi_0(u;n-l,n-r) + \Psi_0(n,n-l) \Phi'_0(u;n-l,n-r)]\}_j.\end{align}
Since from (3.4)
\begin{equation*}\Phi_{\theta}'(x) = x \Phi_{\theta}(x) U_{\theta}',\end{equation*}
the application of (3.22) only requires the calculation of the blocks of 
$\Psi_{\theta}'$ and 
$U_{\theta}'$. This can be done by differentiating the relations (3.6) and (3.7).
4. Final amount of reserves
Our objective in this section is twofold. First, we want to determine the distribution of 
$F(T_e)$, the amount of reserves at the end of the epidemic, regardless of whether ruin has occurred before time 
$T_e$. Second, we will determine the distribution of 
$F(T_e)$ when ruin occurs (or not) before 
$T_e$. As a corollary, we will obtain the corresponding expectations of the final amount of reserves.
We start by introducing two matrices 
$\Psi^{\ast}$ and 
$\Phi^{\ast}(x)$, 
$x\geq 0$, which involve the stopping times (4.1) below for the reserves:
\begin{equation}T_x = \inf\{t>0 \,|\, F(t)=x\}, \quad x\geq 0.\end{equation}
Note that 
$T_0\equiv T_r$. For 
$0 \le s \le v \le n$, 
$i\in {\mathcal{E}}^-_v$, 
$j\in {\mathcal{E}}^+_s$, we define
\begin{equation*}\left(\Psi^{\ast}\right)_{(v,i),(s,j)} \equiv \left(\Psi^{\ast}(v,s)\right)_{i,j} = {\mathbb{P}\left( {T_0
< T_e, Y(T_0)=(s,j) \,|\,F(0)=0, Y(0)=(v,i)} \right)},\end{equation*}
and for 
$0 \le s \le v \le n$, 
$i\in {\mathcal{E}}^+_v$, 
$j\in {\mathcal{E}}^+_s$, we define
\begin{align*}\left(\Phi^{\ast}\right)_{(v,i),(s,j)}(x) &\equiv \left(\Phi^{\ast}(x;v,s)\right)_{i,j} \\ &= {\mathbb{P}\left( {T_{x}
< T_e, Y(T_{x})=(s,j) \,|\,F(0)=0, Y(0)=(v,i)} \right)}.\end{align*}
The former gives the probability of first return to the reserves level 0 with the state 
$Y(T_0)$ before the epidemic end 
$T_e$, under the condition that the reserves process starts at the level 0 in a descending phase (
$i\in {\mathcal{E}}^-_v$). The latter gives the probability of first passage to the higher reserves level x before time 
$T_e$, under the condition that the reserves process starts at the level 0 in an ascending phase (
$i\in {\mathcal{E}}^+_v$).
Like (3.4) for 
$\Phi_{\theta}(x)$, we have
\begin{equation}\Phi^{\ast}(x) = e^{U^{\ast}x}\end{equation}
for a sub-generator 
$U^{\ast}$. The matrices 
$\Psi^{\ast}, \Phi^{\ast}(x), U^{\ast}$ have the same block structures (3.1), (3.2), (3.5) as 
$\Psi_0, \Phi_0(x), U_0$ in Section 3, but are of different dimensions: the block 
$\Psi^{\ast}(v,s)$ is of dimensions 
$|{\mathcal{E}}^-_v|\times |{\mathcal{E}}^+_s|$, and the blocks 
$\Phi^{\ast}(x;v,s), U^{\ast}(v,s)$ are of dimensions 
$|{\mathcal{E}}^+_v|\times |{\mathcal{E}}^+_s|$.
Note that the fluid process is homogeneous with respect to the level, which means that only the distance between the starting level and the target level is relevant. In this sense, the matrices 
$\Psi^{\ast}, \Phi^{\ast}(x), U^{\ast}$ can be considered as complementary to 
$\Psi_0, \Phi_0(x), U_0$. As a result, the blocks of 
$\Psi^{\ast}, U^{\ast}$ can be computed by applying Proposition 3 to the level-reversed fluid model 
$\{\widehat{F}(t), Y(t)\}$, which is defined as 
$\{F(t), Y(t)\}$ but with the rates c(s,i) this time having reversed signs.
4.1. Distribution of $F(T_e)$
Initially, the level of reserves is 
$u\geq 0$, the population contains n susceptibles, and the infection phase is of distribution (
$\sim$) 
$\boldsymbol{\pi}$. At the end time 
$T_e$, the amount of reserves 
$F(T_e)$ may be positive or not, and we aim to calculate 
${\mathbb{P}\left( {F(T_e)>x} \right)}$ where x is any real number.
Since the fluid is homogeneous with respect to the level, we can write
\begin{equation*}{\mathbb{P}\left( {F(T_e)>x \,|\, F(0)=u, S(0)=n, \varphi(0)\sim \boldsymbol{\pi}} \right)} = \boldsymbol{\pi} \textbi{M}(x-u;n), \quad x\in \mathbb{R},\end{equation*}
where each 
$\textbi{M}(x;s)$, for 
$x\in \mathbb{R}$, 
$0 \le s \le n$, denotes a column vector of dimension 
$|{\mathcal{E}}_s|$ with components
\begin{equation*}M_i(x;s) = {\mathbb{P}\left( {F(T_e)>x \,|\, F(0)=0, Y(0)=(s,i)} \right)}, \quad i \in {\mathcal{E}}_s.\end{equation*}
According to the subdivision 
${\mathcal{E}}_s = {\mathcal{E}}_s^+ \cup {\mathcal{E}}_s^-$, the vector 
$\textbi{M}(x;s)$ is partitioned as
\begin{equation*}\textbi{M}(x;s) = \begin{bmatrix} \textbi{M}_+(x;s) \\[3pt] \textbi{M}_-(x;s) \end{bmatrix}.\end{equation*}
We show first that 
$\textbi{M}(x;s)$ can be expressed in terms of only the vectors 
$\textbi{M}(0;l)$ for 
$0\leq l \leq s$.
Proposition 4.1. For 
$x \ge 0$,
\begin{eqnarray}\textbi{M}_+(x;s) &=& \sum_{l=0}^{s} \Phi^{\ast}(x;s,l) \textbi{M}_+(0;l), \end{eqnarray}
\begin{eqnarray}\textbi{M}_-(x;s) &=& \sum_{l=0}^{s} \Psi^{\ast}(s,l) \textbi{M}_+(x;l), \end{eqnarray}
while for 
$x < 0$,
\begin{eqnarray}\textbi{M}_-(x;s) &=& \textbf{1}- \sum_{l=0}^{s} \Phi_0(-x;s,l) \textbf{1}+ \sum_{l=0}^{s} \Phi_0(-x;s,l) \textbi{M}_-(0;l), \end{eqnarray}
\begin{eqnarray}\textbi{M}_+(x;s) &=& \textbf{1}- \sum_{l=0}^{s} \Psi_0(s,l) \textbf{1}+ \sum_{l=0}^{s} \Psi_0(s,l) \textbi{M}_-(x;l). \end{eqnarray}
Proof. Consider (4.3). Starting from 
$F(0)=0$ with s susceptibles and phase in 
${\mathcal{E}}_s^+$, the fluid process F(t) must reach the level x before time 
$T_e$, and at that time there will remain l susceptibles, with 
$0 \le l\le s$. This yields the factor 
$\Phi^{\ast}(x;s,l)$. Next, starting from the reserves x with l susceptibles, the fluid process F(t) must be above the level x at time 
$T_e$, hence the factor 
$\textbi{M}_+(0;l)$.
For (4.4), the argument is similar. Starting from 
$F(0)=0$ with s susceptibles and phase in 
${\mathcal{E}}_s^-$, the reserves process F(t) must first return to zero before time 
$T_e$, say with l susceptibles, hence 
$\Psi^{\ast}(s,l)$. Then, it must be above the level x at time 
$T_e$, hence 
$\textbi{M}_+(x;l)$.
For (4.5), the difference of the first two terms in the right-hand side gives the probabilities that starting from 
$F(0)=0$ with s susceptibles and phase in 
${\mathcal{E}}_s^-$, the reserves process does not reach the level x before time 
$T_e$. The third term contains the probabilities 
$\Phi_0(-x;s,l)$ of reaching the level x with l susceptibles, multiplied by the probabilities 
$\textbi{M}_-(0;l)$ of then being above x at 
$T_e$.
The reasoning to derive (4.6) is analogous.
In fact, we observe from Proposition 4.1 that the computation of all the 
$\textbi{M}(x;s)$ only requires the determination of the vectors 
$\textbi{M}_+(0;l)$. Let us show now that these can be obtained by recursion.
Proposition 4.2. For 
$s=0, \ldots,n$,
\begin{equation}\textbi{M}_+(0;s) = \left[I-\Psi(s,s)\Psi^{\ast}(s,s)\right]^{-1} \textbi{g}(s),\end{equation}
where the vectors 
$ \textbi{g}(s)$ are given by
\begin{eqnarray}\textbi{g}(s) && = \textbf{1}- \sum_{l=0}^{s} \Psi_0(s,l) \textbf{1}+ 1_{s>0} \, \Psi_0(s,s)\sum_{l=0}^{s-1} \Psi^{\ast}(s,l) \textbi{M}_+(0;l)\nonumber\\[3pt]&& + 1_{s>0} \,\sum_{l=0}^{s-1} \sum_{v=0}^{l} \Psi_0(s,l) \Psi^{\ast}(l,v) \textbi{M}_+(0;v).\end{eqnarray}
Proof. Arguing as for (4.4), we first express 
$\textbi{M}_+(0;s)$ as
\begin{equation}\textbi{M}_+(0;s) = \textbf{1}- \sum_{l=0}^{s} \Psi_0(s,l) \textbf{1}+ \sum_{l=0}^{s} \sum_{v=0}^{l} \Psi_0(s,l) \Psi^{\ast}(l,v) \textbi{M}_+(0;v).\end{equation}
Indeed, the difference of the first two terms in the right-hand side gives the probabilities that starting from 
$F(0)=0$ with s susceptibles and phase in 
${\mathcal{E}}_s^+$, the reserves process does not return to the level zero before time 
$T_e$. In the third term, 
$\Psi_0(s,l)$ represents the probabilities of returning to the level zero with l susceptibles. The process then is at level zero in a descending phase, and 
$\Psi^{\ast}(l,v)$ gives the probabilities of reaching zero again (from below) before 
$T_e$ with v susceptibles. Finally, it remains to compute the probabilities 
$\textbi{M}_+(0;v)$.
By rearranging (4.9), we then obtain
\begin{equation*}\left[I-\Psi_0(s,s)\Psi^{\ast}(s,s)\right] \textbi{M}_+(0;s)= \textbi{g}(s),\end{equation*}
with 
$\textbi{g}(s)$ defined by (4.8). Since the matrix 
$\Psi_0(s,s)\Psi^{\ast}(s,s)$ is strictly sub-stochastic, the matrix 
$I-\Psi_0(s,s)\Psi_0^{\ast}(s,s)$ is invertible, hence the formula (4.7) for 
$\textbi{M}_+(0;s)$.
4.2. Distribution of $F(T_e)$ when $T_e>(<) \,T_r$
Under the same initial conditions, we now determine the distribution of the amount of reserves 
$F(T_e)$ in case the company is ruined or not before the final time 
$T_e$. Note that 
$F(T_e)$ is positive if the epidemic ends before ruin (i.e. when 
$T_e<T_r$).
We start by writing
\begin{align*}&{\mathbb{P}\left( {F(T_e)>x, T_r<T_e \,|\, F(0)=u, S(0)=n, \varphi(0)\sim \boldsymbol{\pi}} \right)} = \boldsymbol{\pi} \textbi{V}(x;u,n), \quad x\in \mathbb{R},\\[3pt]&{\mathbb{P}\left( {F(T_e)>x, T_e<T_r \,|\, F(0)=u, S(0)=n, \varphi(0)\sim \boldsymbol{\pi}} \right)} = \boldsymbol{\pi} \textbi{W}(x;u,n),\quad x\geq 0,\end{align*}
where 
$\textbi{V}(x;y,s)$ and 
$\textbi{W}(x;y,s)$, for 
$y\geq 0$ and 
$0 \le s \le n$, denote column vectors of dimension 
$|{\mathcal{E}}_s|$ with components
\begin{align*}& V_i(x;y,s) = {\mathbb{P}\left( {F(T_e)>x, T_r<T_e \,|\, F(0)=y, Y(0)=(s,i)} \right)}, \quad i \in {\mathcal{E}}_s, \\[3pt]& W_i(x;y,s) = {\mathbb{P}\left( {F(T_e)>x, T_e<T_r \,|\, F(0)=y, Y(0)=(s,i)} \right)}, \quad i \in {\mathcal{E}}_s.\end{align*}
Like 
$\textbi{M}(x;s)$ above, 
$\textbi{V}(x;y,s)$ and 
$\textbi{W}(x;y,s)$ can be partitioned according to the subdivision 
${\mathcal{E}}_s = {\mathcal{E}}_s^+ \cup {\mathcal{E}}_s^-$.
Let us determine these vectors. By definition, we directly see that for 
$x\geq 0$,
\begin{equation}\textbi{W}(x;y,s) = \textbi{M}(x-y;s) - \textbi{V}(x;y,s).\end{equation}
On the other hand, 
$\textbi{V}(x;y,s)$, 
$x\in \mathbb{R}$, is provided by the formulas (4.11) and (4.12) below once the vectors 
$\textbi{M}(x;s)$ have been obtained from Propositions 4.1 and 4.2.
Proposition 4.3. For 
$x\in \mathbb{R}$,
\begin{eqnarray}\textbi{V}_-(x;y,s) &=& \sum_{l=0}^{s} \Phi_0(y;s,l) \textbi{M}_-(x;l), \end{eqnarray}
\begin{eqnarray}\textbi{V}_+(x;y,s) &=& \sum_{l=0}^{s} \sum_{r=0}^{l} \Psi_0(s,l) \Phi_0(y;l,r) \textbi{M}_-(x;r). \end{eqnarray}
Proof. From 
$\textbi{V}(x;y,s)$, we know that F(t) starts from y with s susceptibles and phase in 
${\mathcal{E}}_s^-$ or 
${\mathcal{E}}_s^+$, then must go to the level 0 before 
$T_e$, and finally must be above the level x at time 
$T_e$. In probabilistic terms, the first return to 0 is given by the matrices 
$\Phi_0(y;s,l)$ if the initial phase is in 
${\mathcal{E}}_s^-$, and by the matrices 
$\Psi_0(s,l) \Phi_0(y,l,r) $ if the initial phase is in 
${\mathcal{E}}_s^+$. We thus deduce the desired expansions (4.11) and (4.12), respectively.
4.3. Expectation of $F(T_e)$ in both cases
As a corollary, we will determine the expected amount of reserves at the end time 
$T_e$, either taking into account the possible ruin before, or not doing so.
For that, we define by 
$\textbi{h}(s)$, 
$0 \le s \le n$, the column vector of dimension 
$|{\mathcal{E}}_s|$ with components
\begin{equation*}h_i(s) = {\mathbb{E}\left[ {F(T_e) \,|\, F(0)=0, Y(0)=(s,i)} \right]}, \quad i \in {\mathcal{E}}_s.\end{equation*}
Here too, the vector is partitioned according to the subdivision 
${\mathcal{E}}_s = {\mathcal{E}}_s^+ \cup {\mathcal{E}}_s^-$, which gives 
$\textbi{h}(s) \equiv [\textbi{h}_+(s), \textbi{h}_-(s)]^{\tau}$.
Proposition 4.4. Given 
$F(0)=u\geq 0$, 
$S(0)=n$, 
$\varphi(0)\sim \boldsymbol{\pi}$, we have
\begin{align}{\mathbb{E}\left[ {F(T_e)} \right]} &= u + \boldsymbol{\pi}\textbi{h}(n), \end{align}
\begin{align}{\mathbb{E}\left[ {F(T_e) {\mathds{1}_{T_r>T_e}}} \right]} &= {\mathbb{E}\left[ {F(T_e)} \right]} - {\mathbb{E}\left[ {F(T_e) {\mathds{1}_{T_r<T_e}}} \right]}, \end{align}
\begin{align}{\mathbb{E}\left[ {F(T_e) {\mathds{1}_{T_r<T_e}}} \right]} &= \boldsymbol{\pi}_-\sum_{l=0}^{n} \Phi_0(u;n,l) \textbi{h}_-(l) + \boldsymbol{\pi}_+\sum_{l=0}^{n} \sum_{v=0}^{l} \Psi_0(s,l) \Phi_0(u;l,v) \textbi{h}_+(v), \end{align}
where the vectors 
$\textbi{h}_+(s)$, 
$\textbi{h}_-(s)$ are given by
\begin{eqnarray*}\textbi{h}_+(s) &=& \sum_{l=0}^{s} \,[(-U^{\ast})^{-1}]_{(s,l)} \, \textbi{M}_+(0;l) - \sum_{l=0}^{s} \sum_{v=0}^{l} \Psi_0(s,l) \,[(-U_0)^{-1}]_{(l,v)} \,[\textbf{1}-\textbi{M}_-(0;v)], \\[3pt]\textbi{h}_-(s) &=& \sum_{l=0}^{s} \sum_{v=0}^{l} \Psi^{\ast}(s,l) \,[(-U^{\ast})^{-1}]_{(l,v)}\, \textbi{M}_+(0;v) - \sum_{l=0}^{s} \, [(-U_0)^{-1}]_{(s,l)} \,[\textbf{1}-\textbi{M}_-(0;l)].\end{eqnarray*}
Proof. The equality (4.13) is immediate from the homogeneity of the fluid flow with respect to the level. The formula (4.14) is obvious, as was the case for (4.10). The formula (4.15) follows directly from Proposition 4.3. Thus, it remains to determine 
$\textbi{h}(s)$.
First, consider 
$\textbi{h}_+(s)$. The expectation of a random variable X can be represented as
\begin{equation*}{\mathbb{E}\left[ {X} \right]} = \int_0^{\infty}{\mathbb{P}\left( {X>x} \right)} dx - \int_0^{\infty}{\mathbb{P}\left( {X<-x} \right)} dx.\end{equation*}
We apply this identity to 
$X=F(T_e)$ when 
$F(0)=0$ and 
$Y(0)=(s,i)$, 
$i \in {\mathcal{E}}_s^+$. Using the previous notation, this easily leads to the formula
\begin{align*}\textbi{h}_+(s) = \int_0^{\infty}\textbi{M}_+(x;s)\, dx - \int_0^{\infty} \sum_{l=0}^{s} \sum_{v=0}^{l} \Psi_0(s,l) \Phi_0(x;l,v) \,[\textbf{1}-\textbi{M}_-(0;v)] dx.\end{align*}
From Proposition 4.1, we then obtain
\begin{eqnarray}\textbi{h}_+(s) &=& \sum_{l=0}^{s} \left(\int_0^{\infty}\Phi^{\ast}(x;s,l) dx\right) \textbi{M}_+(0;l) \nonumber\\[3pt] \qquad &-& \sum_{l=0}^{s} \sum_{v=0}^{l} \Psi_0(s,l)\left(\int_0^{\infty} \Phi_0(x;l,v) dx \right) [\textbf{1}-\textbi{M}_-(0;v)].\end{eqnarray}
By (3.4) and (4.2), 
$\Phi_0(x) = \exp(U_0x)$ and 
$\Phi^{\ast}(x) = \exp(U^{\ast}x)$, so that
\begin{equation*}\int_0^{\infty}\Phi^{\ast}(x;s,l) \, dx = [(-U^{\ast})^{-1}]_{(s,l)} \quad \mbox{and} \quad \int_0^{\infty}\Phi_0(x;l,v) \, dx = [(-U_0)^{-1}]_{(l,v)}.\end{equation*}
Substituting this in (4.16), we deduce the stated formula. We can proceed in the same way for 
$\textbi{h}_-(s)$.
The blocks 
$[(-U_0)^{-1}]_{(l,v)}$ of 
$(-U_0)^{-1}$, as well as the blocks 
$[(-U^{\ast})^{-1}]_{(l,v)}$ of 
$(-U^{\ast})^{-1}$, can be efficiently calculated using the identities 
$U_0 U_0^{-1}=I$ and 
$U^{\ast}(U^{\ast})^{-1}=I$ together with the block structure (3.5).
5. Fluid flow with Brownian noise
Consider again a SIR epidemic model whose propagation is described by the Markov process 
$\{Y(t)\}$ with generator given by (2.3)–(2.4). For the fluid flow of the reserves, we suppose now that the process 
$\{F(t)\}$ is no longer piecewise linear but is perturbed by a Brownian noise. More precisely, if 
$Y(t)=(s,i)$ at time 
$t<T_e$, we modify (2.7) by defining
\begin{equation}dF(t) = c(s,i) \,dt + \nu(s,i) \, dW(t),\end{equation}
where the c(s,i) are the rates given by (2.6), the 
$\nu(s,i)$ are strictly positive coefficients, and 
$\{W(t)\}$ is a standard Brownian motion. Thus, the reserves process now evolves as a Brownian motion with drift c(s,i) and standard deviation 
$\nu(s,i)$. The addition of a Brownian noise component has the advantage of incorporating some variability in the premiums collected and the medical expenses reimbursed by the insurance.
In the Brownian setting, it is no longer necessary to separate the states according to the sign of the drift. We will work with the generator Q and the matrix of rates C structured according to the partition of states 
$V_n \cup V_{n-1} \cup \cdots \cup V_0$, i.e. of the previous forms (2.4) and (2.8). We also define the matrix 
$\Sigma$ of volatility parameters,
\begin{equation}\Sigma = \begin{bmatrix} \Sigma(n) &\quad 0 &\quad \cdots &\quad 0 \\[3pt]0&\quad \Sigma(n\!-\!1) &\quad \cdots &\quad 0 \\[3pt]\vdots &\quad\vdots &\quad \ddots &\quad \vdots \\[3pt]0&\quad0 &\quad \cdots &\quad \Sigma(0)\end{bmatrix},\end{equation}
where 
$\Sigma(s)$, 
$s=0, \ldots,n$, is the diagonal matrix of dimension 
$|{\mathcal{E}}_s|$ such that 
$\Sigma_{i,i}(s) = \nu(s,i)$, 
$i\in {\mathcal{E}}_s$.
5.1. Ruin time distribution
Let 
$T_x$, 
$x\geq 0$, be the stopping times defined by (4.1), with 
$T_0\equiv T_r$ the time of ruin. We introduce two matrices 
$\Phi_{\theta}(x)$ and 
$\Phi_{\theta}^{\ast}(x)$, 
$\theta \geq 0$, defined by
\begin{align*}(\Phi_{\theta})_{(v,i),(s,j)}(x) &\equiv \left(\Phi_{\theta}(x;v,s)\right)_{i,j}\\[3pt]&= {\mathbb{E}\left[ {e^{-\theta T_0} \,{\mathds{1}_{T_0 < T_e, Y(T_0)=(s,j)}} \,|\,F(0)=x, Y(0)=(v,i)} \right]}, \\[3pt](\Phi_{\theta}^{\ast})_{(v,i),(s,j)}(x) &\equiv \left(\Phi_{\theta}^{\ast}(x;v,s)\right)_{i,j}\\[3pt]&= {\mathbb{E}\left[ {e^{-\theta T_0} \,{\mathds{1}_{T_0
< T_e, Y(T_0)=(s,j)}} \,|\,F(0)=-x, Y(0)=(v,i)} \right]},\end{align*}
where 
$0 \le s \leq v \le n$, 
$i \in {\mathcal{E}}_v$, 
$j \in {\mathcal{E}}_s$. As with (3.4) and (4.2), it is easy to see that
\begin{equation}\Phi_{\theta}(x) = e^{U_{\theta} x} \,\,\,\mbox{ and } \,\,\, \Phi_{\theta}^{\ast}(x) = e^{U_{\theta}^{\ast} x}\end{equation}
for some (sub-)generators 
$U_{\theta}$ and 
$U_{\theta}^{\ast}$. The matrices 
$\Phi_{\theta}(x)$, 
$\Phi_{\theta}^{\ast}(x)$ and 
$U_{\theta}$, 
$U_{\theta}^{\ast}$ can still be written in the block-structured triangular forms (3.2) and (3.5), but here they are of different dimensions 
$|{\mathcal{E}}_v| \times |{\mathcal{E}}_s|$.
Consider 
$T_r< T_e$. The Laplace transform of the time of ruin 
$T_r$ jointly with the epidemic state 
$Y(T_r)$ is given by (5.4) below in terms of the blocks of the first row of 
$\Phi_{\theta}(u)$. Using (5.3), these matrices can be calculated efficiently once the blocks of 
$U_{\theta}$ are obtained from (5.5).
Proposition 5.1. For 
$0\leq r \leq n$, 
$j \in {\mathcal{E}}_{n-r}$, we have
\begin{equation}{\mathbb{E}\left[ {e^{-\theta T_r} \,{\mathds{1}_{T_r
< T_e,\, Y(T_r)=(n-r,j)}}} \right]} = \left\{\boldsymbol{\pi} \Phi_{\theta}(u;n,n-r)\right\}_j,\end{equation}
where for 
$0\leq k\leq s \leq n$, the matrices 
$U_{\theta}(s,s-k)$ satisfy the identities
\begin{equation}\begin{aligned}& 1_{k=0}\,\, 2 (\Sigma(s))^{-2} (A(s)-\theta I) + 1_{k=1}\,\, 2 (\Sigma(s))^{-2}B(s) + 2(\Sigma(s))^{-2}C(s) U_{\theta}(s,s-k) \\& \qquad\qquad\qquad\qquad + \sum_{l=0}^{k} U_{\theta}(s,s-l) U_{\theta}(s-l, s-k) = 0.\end{aligned}\end{equation}
Proof. The formula (5.4) is obtained by the same reasoning as (3.3) in Proposition 3. To derive (5.5), we will use a well-known Wiener–Hopf factorization for the Brownian motion, which is stated below (see e.g. Sato [22, Corollary 45.8]).
Let 
$\{B(t)\}$ be a Brownian motion with drift d and standard deviation 
$\sigma$, and let 
$\tau$ be an independent exponential variable of parameter 
$\mu$. Define the two random variables
\begin{equation}Z_1 = - \inf_{0 \le t \le \tau} B(t) \,\,\,\mbox{ and } \,\,\, Z_2 = B(\tau) - \inf_{0 \le t \le \tau} B(t).\end{equation}
Then 
$Z_1$ and 
$Z_2$ are independent exponential variables with respective parameters 
$\omega$ and 
$\eta$ given by
\begin{equation}\omega = \frac{d}{\sigma^2} + \sqrt{\frac{d^2}{\sigma^4} + \frac{2 \mu}{\sigma^2}} \,\,\,\mbox{ and } \,\,\, \eta= \frac{-d}{\sigma^2} + \sqrt{\frac{d^2}{\sigma^4} + \frac{2 \mu}{\sigma^2}}.\end{equation}
Let us fix s and k. To begin, we note that 
$\Phi_{\theta}(x;s,s-k)$ can be expressed as the matrix of first passage probabilities to the level 0, starting from the level 
$x>0$, under the constraint that this first passage occurs before the ringing of an exponential clock of parameter 
$\theta$. In other words,
\begin{equation*}\left(\Phi_{\theta}(x;s,s-k)\right)_{i,j} = {\mathbb{P}\left( {T_r
< \min(T_e,\kappa_\theta), Y(T_r)=(s-k,j) \,|\,F(0)=x, Y(0)=(s,i)} \right)},\end{equation*}
where 
$\kappa_\theta$ is an independent exponential variable of parameter 
$\theta$. Thanks to this interpretation, we will determine 
$(\Phi_{\theta}(x;s,s-k))_{i,j}$ using the Wiener–Hopf factorization above to analyze the possible paths before the first moment 
$\tau$ where either there is a phase transition or the clock rings. A similar argument was applied by Asmussen [Reference Asmussen2] in a different context.
For each 
$i \in {\mathcal{E}}_s$, we define the parameters 
$\omega_i$ and 
$\eta_i$ provided by (5.7) for a Brownian motion of drift 
$d = c(s,i)$ and standard deviation 
$\sigma = \nu(s,i)$, and an exponential variable 
$\tau$ of parameter 
$\mu = \theta - A_{i,i}(s)$. Starting at the level x in phase i, the reserves process may reach the level 0 before time 
$\tau$, or not. According to the result related to (5.6), the first case occurs with a probability equal to 
${\mathbb{P}\left( {Z_1>x} \right)} = \exp(-\omega_i x)$. In the second case, we use this result again, so that by conditioning on 
$Z_1$, 
$Z_2$, and the new phase after time 
$\tau$, we get
\begin{align}& \left(\Phi_{\theta}(x;s,s-k)\right)_{i,j} =1_{i=j} \, 1_{k=0} \, e^{-\omega_i x} \nonumber \\[3pt]& \qquad + \sum_{r \in {\mathcal{E}}_s} \int_{0}^{x} \int_{0}^{\infty} \omega_i e^{-\omega_i u} \eta_i e^{-\eta_i v} P_{i,r} \, \left(\Phi_{\theta}(x-u+v;s,s-k)\right)_{r,j} \, dv \, du \\[3pt] & \qquad + 1_{s>0} \, \sum_{r \in {\mathcal{E}}_{s-1}} \int_{0}^{x} \int_{0}^{\infty} \omega_i e^{-\omega_i u} \eta_i e^{-\eta_i v} \hat{P}_{i,r} \, \left(\Phi_{\theta}(x-u+v;s-1,s-k)\right)_{r,j} \, dv \, du, \nonumber \end{align}
where P and 
$\hat{P}$ are the transition probability matrices of the jump chain from 
${\mathcal{E}}_{s}$ to 
${\mathcal{E}}_{s}$ and 
${\mathcal{E}}_{s-1}$, i.e.
\begin{align*}&P_{i,r} = {\mathbb{P}\left( {\zeta < \kappa_\theta, Y(\zeta)=(s,r) \,|\, Y(0)=(s,i)} \right)}, \quad i,r \in {\mathcal{E}}_{s}, \\[3pt]&\hat{P}_{i,r} = {\mathbb{P}\left( {\zeta
< \kappa_\theta, Y(\zeta)=(s-1,r) \,|\, Y(0)=(s,i)} \right)}, \quad i\in {\mathcal{E}}_{s}, r\in {\mathcal{E}}_{s-1},\end{align*}
where 
$\zeta$ is the first time that there is a transition in the process 
$\{Y(t)\}$. Therefore, we have from (2.3) and (2.4) that
\begin{equation}P = I + \Lambda^{-1} (A(s)-\theta I) \,\,\,\mbox{ and } \,\,\, \hat{P} = \Lambda^{-1} B(s),\end{equation}
where 
$\Lambda$ is the diagonal matrix of terms 
$\theta-A_{i,i}(s)$. Let us differentiate each side of (5.8) with respect to x and take 
$x=0$, also using (5.3). We then obtain, in matrix notation,
\begin{equation}U_{\theta}(s,s-k) = -\Delta_{\omega} \, 1_{k=0} + \Delta_{\omega} L,\end{equation}
where
\begin{equation*}L = \int_{0}^{\infty} \Delta_{\eta} e^{-\Delta_{\eta}v} \left(P \Phi_{\theta}(v;s,s - k) + \hat{P} \Phi_{\theta}(x;s-1,s-k) 1_{s>0} \right) dv,\end{equation*}
and 
$\Delta_{\omega}$ and 
$\Delta_{\eta}$ are the diagonal matrices of the terms 
$\omega_i$ and 
$\eta_i$, respectively. Now, integrating L by parts, with (3.9) and 
$\Phi_{\theta}(0;s-l,s-k)=I \, 1_{k=l}$, we can rewrite (5.10) as
\begin{align*}& U_{\theta}(s,s-k) = 1_{k=0} \, \Delta_{\omega} (P-I) + 1_{k=1} \, \Delta_{\omega} \hat{P} + \Delta_{\omega} \int_{0}^{\infty} e^{-\Delta_{\eta}v} P \Phi_{\theta}(v;s,s) dv \, U_{\theta}(s,s-k)\\[3pt]& \quad + \sum_{l=1}^{k} \Delta_{\omega} \int_{0}^{\infty} e^{-\Delta_{\eta}v} \left(P \Phi_{\theta}(v;s,s - l) + 1_{s>0}\, \hat{P} \Phi_{\theta}(x;s-1,s-l) \right) dv \, U_{\theta}(s-l,s-k).\end{align*}
Using (5.10) with l substituted for k, this relation can be expressed as
\begin{eqnarray}&& 1_{k=0}\, \Delta_{\eta} \Delta_{\omega} (P-I) + 1_{k=1}\, \Delta_{\eta}\Delta_{\omega} +\hat{P}(\Delta_{\omega}-\Delta_{\eta}) U_{\theta}(s,s-k) \nonumber \\[3pt]&& \qquad\qquad + \sum_{l=0}^{k} U_{\theta}(s,s-l) \, U_{\theta}(s-l,s-k) = 0.\end{eqnarray}
Finally, inserting in (5.11) the identities (5.9),
\begin{equation*}\Delta_{\omega} - \Delta_{\eta} = 2(\Sigma(s))^{-2} C(s), \,\,\,\mbox{ and } \,\,\, \Delta_{\eta} \Delta_{\omega} = 2(\Sigma(s))^{-2} \Lambda,\end{equation*}
we deduce the desired relations (5.5).
The blocks of 
$U_{\theta}$ can be computed recursively from (5.5) as follows.
Induced algorithm.
(i) First, for
$0\leq s \leq n$, taking $k=0$ implies that $Z\equiv U_{\theta}(s,s)$ is the unique sub-generator satisfying the quadratic equation
Various numerical procedures have been developed in the literature to solve such an equation numerically (see e.g. Asmussen [Reference Asmussen2] or Latouche and Nguyen [Reference Latouche and Nguyen15]).
\begin{equation*}Z^2 + 2 (\Sigma(s))^{-2}C(s)Z + 2 (\Sigma)^{-2}(A(s)-\theta I)= 0.\end{equation*}
-
(ii) Next, for
$1 \le s \le n$, taking $k=1$ shows that $Z\equiv U_{\theta}(s,s-1)$ is the unique solution of the Sylvester equation
\begin{equation*}[U_{\theta}(s,s) + 2(\Sigma(s))^{-2}C(s)] Z + ZU_{\theta}(s-1,s-1) = -2 (\Sigma(s))^{-2}B(s).\end{equation*}
-
(iii) Finally, for
$k=2, \ldots, s$, we find that $Z\equiv U_{\theta}(s,s-k)$ is the unique solution of the Sylvester equation
\begin{equation*}[U_{\theta}(s,s) + 2(\Sigma(s))^{-2}C(s)] Z + ZU_{\theta}(s-k,s-k) = -\sum_{l=1}^{k-1}U_{\theta}(s,s-l)U_{\theta}(s-l,s-k).\end{equation*}
5.2. Distribution of $F(T_e)$
Let us turn to the final amount of reserves. Similarly to the fluid model, we have 
${\mathbb{P}\left( {F(T_e)>x} \right)} = \boldsymbol{\pi} \textbi{M}(x-u;n)$ for 
$x\in \mathbb{R}$, where 
$\textbi{M}(x;s)$ is a column vector with components
\begin{equation*}M_i(x;s) = {\mathbb{P}\left( {F(T_e)>x \,|\, F(0)=0, \, Y(0)=(s,i)} \right)}, \qquad i \in {\mathcal{E}}_s.\end{equation*}
Here again, we no longer separate 
$\textbi{M}(x;s)$ according to the sign of the drift.
We will express 
$\textbi{M}(x;s)$ in terms of only the vectors 
$\textbi{M}(0;l)$, 
$0\leq l \leq s$, which can then be determined by recursion. Before that, we notice that, as in the fluid case, 
$\Phi^{\ast}_0(x)$ is defined as 
$\Phi_0(x)$ but for the level-reversed model 
$\{{\hat F}(t), Y(t)\}$ built with rates of reversed signs. Thus, the blocks of 
$U_0^{\ast}$ are computed from (5.5) after replacing C(s) by 
$-C(s)$.
Proposition 5.2. For 
$x \ge 0$,
\begin{equation}\textbi{M}(x;s) = \sum_{l=0}^{s} \Phi_0^{\ast}(x;s,l) \textbi{M}(0;l),\end{equation}
while for 
$x \le 0$,
\begin{equation}\textbi{M}(x;s) = \textbf{1}- \sum_{l=0}^{s} \Phi_0(-x;s,l) \textbf{1}+ \sum_{l=0}^{s} \Phi_0(-x;s,l) \textbi{M}(0;l).\end{equation}
The vectors 
$\textbi{M}(0;s)$, 
$0\leq s \leq n$, are provided recursively from
\begin{align}\textbi{M}(0;s) = \left[U_0(s,s)+U_0^{\ast}(s,s)\right]^{-1}\, \{\sum_{l=0}^{s} U_0(s,l) \textbf{1} - \sum_{l=0}^{s-1} [U_0(s,l)+U_0^{\ast}(s,l)] \textbi{M}(0;l)\}.\end{align}
Proof. The formulas (5.12) and (5.13) are obtained using the same argument as for Proposition 4.1. To get (5.14), we will proceed by approximation.
Fix 
$\varepsilon >0$ and define the vectors 
$\textbi{M}_{\varepsilon}(0;s)$ as 
$\textbi{M}(0;s)$ but with the following additional constraint: the reserves process, which starts at level 0, each time it reaches 0 before 
$T_e$, must also reach the level 
$\varepsilon$ before 
$T_e$. In other words,
\begin{equation*}\textbi{M}_{\varepsilon, i}(0;s) = {\mathbb{P}\left( {F(T_e)>x, \, \kappa_{\varepsilon} \ge \kappa_0 \,|\, F(0)=0, Y(0)=(s,i)} \right)},\end{equation*}
where 
$\kappa_y$ denotes the time of the last visit at level y before 
$T_e$. Let us first show that 
$\textbi{M}_{\varepsilon}(0;s)$ can be expressed as
\begin{equation} \textbi{M}_{\varepsilon}(0;s) = \sum_{l=0}^{s} \Phi_0^{\ast}(\varepsilon;s,l) [\textbf{1} - \sum_{r=0}^{l} \Phi_0(\varepsilon;l,r) \textbf{1}] + \sum_{l=0}^{s} \sum_{r=0}^{l} \Phi_0^{\ast}(\varepsilon;s,l) \Phi_0(\varepsilon,l,r) \textbi{M}_{\varepsilon}(0;r).\end{equation}
Indeed, in the right-hand side, the first term gives the probability that the reserves process reaches 
$\varepsilon$ and does not return to 0 before 
$T_e$, while the second term gives the probability that the process reaches 
$\varepsilon$ and return to 0 before 
$T_e$. By grouping the 
$\textbi{M}_{\varepsilon}(0;s)$ on both sides of (5.15), we get
\begin{align}\textbi{M}_{\varepsilon}(0;s)&= [I-\Phi_0^{\ast}(\varepsilon;s,s) \Phi_0(\varepsilon;s,s)]^{-1} \{ \sum_{l=0}^{s} \Phi_0^{\ast}(\varepsilon;s,l) [\textbf{1} - \sum_{r=0}^{l} \Phi_0(\varepsilon;l,r) \textbf{1}] \nonumber\\& \qquad\qquad + \sum_{l=0}^{s} \sum_{r=0}^{l}\Phi_0^{\ast}(\varepsilon;s,l) \Phi_0(\varepsilon;l,r) \textbi{M}_{\varepsilon}(0;r) (1-\delta_{s,l,r})\},\end{align}
where 
$\delta_{s,l,r}$ equals 1 when 
$s=l=r$, and 0 otherwise.
Now, from (5.3) with 
$\theta=0$, for 
$\epsilon$ small we have
\begin{equation}\Phi_0(\varepsilon;s,l) = I \delta_{s,l} + \varepsilon U_0(s,l) + o(\varepsilon) \,\,\,\mbox{ and } \,\,\, \Phi_0^{\ast}(\varepsilon;s,l) = I \delta_{s,l} + \varepsilon U_0^{\ast}(s,l) + o(\varepsilon).\end{equation}
By inserting (5.17), we obtain in (5.16) the following approximations:
\begin{align*}I-\Phi_0^{\ast}(\varepsilon;s,s) \Phi_0(\varepsilon;s,s) &= -\varepsilon [U_0(s,s)+U_0^{\ast}(s,s)] + o(\varepsilon),\\\sum_{l=0}^{s} \Phi_0^{\ast}(\varepsilon;s,l) [\textbf{1} - \sum_{r=0}^{l} \Phi_0(\varepsilon;l,r) \textbf{1}] &= -\varepsilon \sum_{l=0}^{s} U_0(s,l) \textbf{1} + o(\varepsilon),\\\sum_{l=0}^{s} \sum_{r=0}^{l} \Phi_0^{\ast}(\varepsilon;s,l) \Phi_0(\varepsilon;l,r) \textbi{M}_{\varepsilon}(0;r) (1-\delta_{s,l,r}) &= \varepsilon \sum_{l=0}^{s-1} [U_0(s,l)+U_0^{\ast}(s,l)]\textbi{M}(0;l) + o(\varepsilon).\end{align*}
Consequently, when 
$\varepsilon \to 0^+$, the identity (5.14) is reduced to the formula (5.14).
Furthermore, we can determine the distribution of the amount 
$F(T_e)$ and its expectation, when ruin does or does not occur before 
$T_e$, in terms of the vectors 
$\textbi{M}(x;s)$. The results are similar to those in Propositions 4.3 and 4.4; the proofs are omitted.
Proposition 5.3. Given 
$F(0)=u\geq 0$, 
$S(0)=n$, 
$\varphi(0)\sim \boldsymbol{\pi}$,
\begin{align*}&{\mathbb{P}\left( {F(T_e)>x, T_r<T_e} \right)} = \boldsymbol{\pi} \sum_{l=0}^{n} \Phi_0(u;n,l) \textbi{M}(x;l), \\&{\mathbb{P}\left( {F(T_e)>x, T_e<T_r} \right)} = \boldsymbol{\pi} \textbi{M}(x-u;n) - \boldsymbol{\pi} \sum_{l=0}^{n} \Phi_0(u;n,l) \textbi{M}(x;l).\end{align*}
For the means, 
${\mathbb{E}\left[ {F(T_e)} \right]}$ and 
${\mathbb{E}\left[ {F(T_e) {\mathds{1}_{T_r>T_e}}} \right]}$ are provided by (4.13)–(4.14), and
\begin{equation*}{\mathbb{E}\left[ {F(T_e) {\mathds{1}_{T_r<T_e}}} \right]} = \boldsymbol{\pi}\sum_{l=0}^{n} \Phi_0(u;n,l) \textbi{h}(l),\end{equation*}
where the vectors 
$\textbi{h}(s)$ are given by
\begin{align*}\textbi{h}(s) &= \sum_{l=0}^{s} \left[(-U_0^{\ast})^{-1}\right]_{(s,l)} \textbi{M}(0;l) - \sum_{l=0}^{s} \left[(-U_0)^{-1}\right]_{(s,l)} [\textbf{1}-\textbi{M}(0;l)].\end{align*}
6. Numerical illustrations
Let us consider an extension of the general epidemic 
$\{S(t),I(t)\}$ where the infinitesimal transition probabilities are
\begin{equation}\begin{aligned}& {\mathbb{P}\left( {S(t+dt)=s, I(t+dt)=i-1 \,|\, S(t)=s, I(t)=i} \right)}= \mu(s,i) \, dt + o(dt), \\& {\mathbb{P}\left( {S(t+dt)=s-1, I(t+dt)=i+1 \,|\, S(t)=s, I(t)=i} \right)} =\beta(s,i)\, dt + o(dt),\end{aligned}\end{equation}
starting from 
$S(0)=n$, 
$I(0)=m$. Thus, the removal rates 
$\mu(s,i)$ and the infection rates 
$\beta(s,i)$ are arbitrary functions of the current state (s,i) with, of course, 
$\mu(s,0)=\beta(s,0)=\beta(0,i)=0$. For the standard general epidemic,
\begin{equation}\mu(s,i)=\mu i \,\,\,\mbox{ and } \,\,\, \beta(s,i)=\beta si/(n+m),\end{equation}
where the infection rate has a factor 
$1/(n+m)$ because permanent immunity occurs after infection. For the so-called fatal epidemic,
\begin{equation}\mu(s,i)=\mu i \,\,\,\mbox{ and } \,\,\, \beta(s,i)=\beta si/(s+i),\end{equation}
where the infection rate has now a factor 
$1/(s+i)$ because death inevitably follows infection (see e.g. Picard and Lefèvre [Reference Picard and Lefèvre20]). These two epidemic models are used below to illustrate some of the theoretical results obtained.
The epidemic process (6.1) is a particular case of the model built in Section 2.1. The variable 
$\varphi(t)$ here represents the number I(t) of infectives; therefore 
${\mathcal{E}}_s = \{1, \ldots, n+m-s\}$. The generator of the Markovian process 
$\{S(t),I(t)\}$ is then defined by (2.3) and (2.4) in which 
$\textbi{v}(s) = [\mu(s,1), 0, \ldots, 0]^{\tau}$ is a column vector of dimension 
$n+m-s$, B(s) is the 
$(n+m-s)\times (n+m-s+1)$ matrix
\begin{equation*}B(s) = \begin{bmatrix}[0.8]0 &\quad \beta(s,1) &\quad 0 &\quad 0 &\quad 0 &\quad \cdots &\quad 0 \\[2pt]0 &\quad 0 &\quad \beta(s,2) &\quad 0 &\quad 0 &\quad\cdots &\quad 0 \\[2pt]0 &\quad 0 &\quad 0 &\quad \beta(s,3) &\quad 0 &\quad\cdots &\quad 0 \\[2pt]0 &\quad 0 &\quad 0 &\quad 0 &\quad\beta(s,4) &\quad\cdots &\quad 0 \\[2pt]\vdots &\quad \vdots &\quad \vdots &\quad \vdots &\quad \vdots &\quad \ddots &\quad \vdots &\quad \\[2pt]0 &\quad 0 &\quad 0 &\quad 0 &\quad 0 &\quad \cdots &\quad \beta(s,n\! +\! m \! - \! s)\end{bmatrix},\end{equation*}
and A(s) is the 
$(n+m-s)\times (n+m-s)$ matrix
\begin{equation*}A(s) = \begin{bmatrix}[0.8]-\gamma(s,1) &\quad 0 &\quad 0 &\quad 0 &\quad \cdots &\quad 0 \\[2pt]\mu(s,2) &\quad -\gamma(s,2) &\quad 0 &\quad 0 &\quad\cdots &\quad 0 \\[2pt]0 &\quad \mu(s,3) &\quad -\gamma(s,3) &\quad 0 &\quad\cdots &\quad 0 \\[2pt]0 &\quad 0 &\quad \mu(s,4) &\quad -\gamma(s,4) &\quad\cdots &\quad 0 \\[2pt]\vdots &\quad \vdots &\quad \vdots &\quad \vdots &\quad \ddots &\quad \vdots &\quad \\[2pt]0 &\quad 0 &\quad 0 &\quad 0 &\quad \cdots &\quad -\gamma(s,n\! +\! m \! - \! s)\end{bmatrix},\end{equation*}
where 
$\gamma(s,i) = \beta(s,i) + \mu(s,i)$ for all 
$i\geq 1$. From (2.2), 
$\boldsymbol{\pi}$ is the row vector of dimension m given by 
$[0, \ldots, 0, 1]$.
The associated risk process is a particular case of the fluid flow model of Section 2.2. The reserves process 
$\{F(t)\}$ evolves as a function of the epidemic 
$\{S(t),I(t)\}$ via the (conditional) differential equation (2.7) with 
$c(s,i)=a(s,i)-b(s,i)$ as net income rate (see (2.6)). For the examples below, each susceptible is assumed to pay premiums at a fixed rate a, so that 
$a(s,i)=as$. On the other side, we suppose that each infective is reimbursed for their healthcare costs either at a constant rate b, which yields 
$b(s,i)=bi$, or at a variable rate bi, which yields 
$b(s,i)=bi^2$. Indeed, the number of infected cases could influence the level of care provided by medical units, and therefore the associated cost. With a rate bi, the cost of care per individual increases linearly with the number of infectives present. A reverse influence could also be possible, depending on the situation.
The premium level a per susceptible can be set according to various criteria. For example, the usual equivalence principle in life insurance requires that 
${\mathbb{E}\left[ {\mbox{benefit outgo}} \right]} = {\mathbb{E}\left[ {\mbox{premium income}} \right]}$, which gives the net premium rate
\begin{equation}a^* = \frac{{\mathbb{E}\left[ {\int_0^{T_e} b(I(u))\, du} \right]}}{{\mathbb{E}\left[ {\int_0^{T_e} S(u)\, du} \right]}},\end{equation}
where b(I(u)) is equal to bI(u) or 
$b(I(u))^2$. This case is examined below, as it was in Lefèvre et al. [Reference Lefèvre, Picard and Simon18] and Lefèvre and Simon [Reference Lefèvre and Simon19]. In addition, we also show the influence of the premium level a by considering it as a free parameter with values around 
$a^*$.
Consider a population initially containing 
$n=47$ susceptibles and 
$m=3$ infectives, and suppose that the epidemic is described by a general or fatal model with removal and infection rates (6.2) or (6.3) with 
$\mu = 1$ and 
$\beta = 1.5$. Figure 1 shows the graphs of the distribution of 
$S(T_e)$, the number of suceptibles who ultimately escaped infection. This distribution was determined by applying the results obtained in Section 2.2 of Lefèvre and Simon [Reference Lefèvre and Simon19]. As is well known, 
$S(T_e)$ typically has a bimodal distribution, for both models. The fatal epidemic is expected to be more serious than the general epidemic because the infection rates are larger (since 
$1/(s + i)\geq 1/(n + m)$). This explains the fact that 
$S(T_e)$ presents a mode close to 0 and of high probability.
Figure 1. Distribution of 
$S(T_e)$, the final number of susceptibles, for the general epidemic (left) and the fatal epidemic (right) when 
$n = 47,\,\,m = 3$ and 
$m = 1,\,\,\beta = 1.5$.
In the following, we focus on the reserves process assuming different initial values u, the premium rate 
$a(s,i)=as$, and the reimbursement rate for care 
$b(s;i)=bi$ or 
$b(s;i)=bi^2$. We choose 
$b=0.5$ and a, when it varies, in the interval 
$[0, 2a^{\ast}]$ where 
$a^{\ast}$ is the reference premium level (6.4) provided in Table 1.
Table 1. Values of 
$a^{\ast}$ for the general and fatal epidemics when 
$n=47$, 
$m=3$ and 
$\mu=1$, 
$\beta=1.5$, for 
$b(s,i)=bi$ or 
$b(s,i)=bi^2$ with 
$b=0.5$.
Figure 2. Ruin probability, 
${\mathbb{P}\left( {T_r<T_e} \right)}$, in the general epidemic (
$n=47$, 
$m=3$, 
$\mu=1$, 
$\beta=1.5$), as a function of the premium rate a when 
$b(s,i)=bi$ (left) and 
$b(s,i)=bi^2$ (right) with 
$b=0.5$, for 
$u=5, 10, 20$.
Figure 3. Expected reserves in case of non-ruin, 
${\mathbb{E}\left[ {F(T_e) {\mathds{1}_{T_r>T_e}}} \right]}$, in the general epidemic (
$n=47$, 
$m=3$, 
$\mu=1$, 
$\beta=1.5$), as a function of the premium rate a when 
$b(s,i)=bi$ (left) or 
$b(s,i)=bi^2$ (right) with 
$b=0.5$, for 
$u=5, 10, 20$.
Figures 2, 3, and 4 relate to the general epidemic. Figure 2 gives the graphs of 
${\mathbb{P}\left( {T_r<T_e} \right)}$, the probability that ruin occurs, as a function of the premium rate a for the two previous reimbursement rates with 
$b=0.5$ and when 
$u=5, 10, 20$. The calculations are performed using the formula (3.3) derived in Section 3 (with 
$\theta=0$). We observe that this probability decreases with a, and is higher when u is low or when the reimbursement is large (the case with 
$bi^2$), which is in agreement with intuition. Figure 3 gives the graphs of 
${\mathbb{E}\left[ {F(T_e) {\mathds{1}_{T_r>T_e}}} \right]}$, the expected reserves when ruin does not occur, under the same parametric conditions. The expectations are calculated using Proposition 4.4 in Section 4. As guessed, 
${\mathbb{E}\left[ {F(T_e) {\mathds{1}_{T_r>T_e}}} \right]}$ increases with a, and is higher when u is high or when the reimbursement is small (the case with bi). Figure 4 gives the graphs of 
${\mathbb{P}\left( {T_r<T_e} \right)}$ as in Figure 2, but this time as a function of the infection rate 
$\beta$ and when 
$a = a^{\ast}(\beta)$ is provided by (6.4). We see that the probability function is bell-shaped, i.e. it first increases and then decreases with 
$\beta$, and the decrease is weaker when the reimbursement rate is large. A likely reason is that the premium rate 
$a^{\ast}(\beta)$ increases with the infection rate 
$\beta$, and after a certain level of 
$\beta$, 
$a^{\ast}(\beta)$ is high enough to decrease the risk of ruin.
Figure 4. Ruin probability, 
${\mathbb{P}\left( {T_r<T_e} \right)}$, in the general epidemic (
$n=47$, 
$m=3$, 
$\mu=1$, 
$\beta=1.5$), as a function of the infection rate 
$\beta$ when 
$b(s,i)=bi$ (left) and 
$b(s,i)=bi^2$ (right) with 
$b=0.5$, for 
$a = a^{\ast}(\beta)$ and 
$u=5, 10, 20$.
Figure 5. Ruin probability, 
${\mathbb{P}\left( {T_r<T_e} \right)}$, in the fatal epidemic (
$n=47$, 
$m=3$, 
$\mu=1$, 
$\beta=1.5$), as a function of the premium rate a when 
$b(s,i)=bi$ (left) and 
$b(s,i)=bi^2$ (right) with 
$b=0.5$, for 
$u=5, 10, 20$.
Figure 6. Expected reserves in case of non-ruin, 
${\mathbb{E}\left[ {F(T_e) {\mathds{1}_{T_r>T_e}}} \right]}$, in the fatal epidemic (
$n=47$, 
$m=3$, 
$\mu=1$, 
$\beta=1.5$) as a function of the premium rate a when 
$b(s,i)=bi$ (left) and 
$b(s,i)=bi^2$ (right) with 
$b=0.5$, for 
$u=5, 10, 20$.
Figure 7. 
${\mathbb{P}\left( {T_r<T_e} \right)}$ and 
${\mathbb{E}\left[ {F(T_e) {\mathds{1}_{T_r>T_e}}} \right]}$ in the fatal epidemic (
$n=47$, 
$m=3$, 
$\mu=1$, 
$\beta=1.5$) with Brownian noise, as a function of the initial reserves u when 
$b(s,i)=bi$ with 
$b = 0.5$, for 
$a=a^*=0.07$ and 
$\nu(s,i)= \nu i$ with 
$\nu = 0, 0.5, 1$.
Figures 5, 6, and 7 relate to the fatal epidemic. Figures 5 and 6 show the graphs of 
${\mathbb{P}\left( {T_r<T_e} \right)}$ and 
${\mathbb{E}\left[ {F(T_e) {\mathds{1}_{T_r>T_e}}} \right]}$ for the fatal epidemic, as did Figures 2 and 3 for the general epidemic. Numerical observations give rise to similar comments. Note that the infection is significantly more severe in the fatal model, so the probability of ruin is much higher and the expected reserves are much lower. Figure 7 gives the graph of these two statistics as functions of the initial reserves u for the model of Section 5 with a Brownian noise of standard deviation 
$\nu(s,i)= \nu i$, say, with 
$\nu = 0, 0.2, 0.5$. An interesting point is that a greater noise dispersion implies, when u is large, a greater probability of ruin and larger expected reserves. Note also that the amounts of expected reserves are very close to the initial level of reserves, which is a consequence of our choice 
$a=a^{\ast}$ as the premium rate.
Acknowledgements
We thank the editor and the referees for their valuable comments and suggestions. C. Lefèvre received partial support from the DIALog Excellence Chair, sponsored by CNP Assurances. M. Simon was supported by the Australian Research Council Center of Excellence for Mathematical and Statistical Frontiers (ACEMS).
