# Analysis of fluid flow models

Guy Latouche Université libre de Bruxelles, Département d’informatique, CP
212, Boulevard du Triomphe, 1050 Bruxelles, Belgium<EMAIL_ADDRESS>Giang
T. Nguyen111The University of Adelaide, School of Mathematical Sciences, SA
5005, Australia<EMAIL_ADDRESS>
###### Abstract

Markov-modulated fluids have a long history. They form a simple class of
Markov additive processes, and were initially developed in the 1950s as models
for dams and reservoirs, before gaining much popularity in the 1980s as models
for buffers in telecommunication systems, when they became known as fluid
queues. More recent applications are in risk theory and in environmental
studies.

In telecommunication systems modelling, the attention focuses on determining
the stationary distribution of the buffer content. Early ODE resolution
techniques have progressively given way to approaches grounded in the analysis
of the physical evolution of the system, and one only needs now to solve a
Riccati equation in order to obtain several quantities of interest. To the
early algorithms proposed in the Applied Probability literature, numerical
analysts have added new algorithms, improved in terms of convergence speed,
numerical accuracy, and domain of applicability.

We give here a high-level presentation of the matrix-analytic approach to the
analysis of fluid queues, briefly address computational issues, and conclude
by indicating how this has been extended to more general processes.

##  1 Introduction

Markov-modulated processes are popular because they can be used to describe
the evolution of simple systems under conditions that vary in time. The fluid
flow processes presented here have long found applications as models for dams
and reservoirs (Loynes [42]), and for buffers in telecommunication systems,
Anick et al. [2] being a famous early paper. It is in the latter context that
the term fluid queue was coined. Later, the domain of applicability of fluid
flows has been extended to risk theory (Avram and Usábel [5], Badescu et al.
[6]), operations management (Bean et al. [12]) and others.

Fluid flows are two-dimensional processes
$\\{X(t),\varphi(t):t\in\mathbb{R}^{+}\\}$, where $\\{\varphi(t)\\}$ is a
continuous-time, irreducible Markov chain on some finite state space
$\mathcal{S}=\\{1,\ldots,m\\}$, and $X(t)$ takes values in $\mathbb{R}$ under
the control of $\varphi$. In the simplest form,

$X(t)=X(0)+\int_{0}^{t}c_{\varphi(s)}\,\mathrm{d}s,$ (1)

with $\boldsymbol{c}=[c_{i}:i\in\mathcal{S}]$ being a vector of arbitrary real
constants. The component $X$ is called the level, $\varphi$ is called the
phase, and the level is a piecewise-linear function, with constant slope over
intervals when the phase is constant. An example with four phases is shown on
Figure 1.

It will be useful in the sequel to partition $\mathcal{S}$ into three subsets
according to the sign of the rates $c_{i}$:

$\mathcal{S}_{\mbox{\tiny$+$}}=\\{i\in\mathcal{S}:c_{i}>0\\},\quad\mathcal{S}_{\mbox{\tiny$-$}}=\\{i\in\mathcal{S}:c_{i}<0\\},\quad\mathcal{S}_{0}=\\{i\in\mathcal{S}:c_{i}=0\\}.$
(2)

$t$$t$$X(t)$$\varphi(t)$$4$$3$$2$$1$

Figure 1: Sample trajectory of a fluid flow (top) and of its controlling
Markov chain (bottom). The fluid rates vector is
$\boldsymbol{c}=[-0.8,-1.4,2,1]$ so that different states correspond to
different slopes.

In many applications, $X$ is a model for a physical quantity (water in a
reservoir, packets in a buffer, energy level of a battery, etc.) and may not
take negative values. In such cases, one might assume that whenever the level
becomes equal to 0 and the rate is negative, the level remains equal to 0
until there is a change to a phase with positive rate. This is illustrated in
Figure 2: at time $\theta_{1}$ the fluid hits level 0, the phase is equal to 1
with $c_{1}=-0.8$. The fluid remains equal to 0 until time $\delta_{1}$ where
the phase process switches to state 3, with $c_{3}=2$.

$t$$\delta_{0}$$\theta_{1}$$\delta_{1}$$\theta_{2}$$Y(t)$

Figure 2:  Sample trajectory of a regulated fluid flow

Such a mechanism is justified by the fact that the rates $c_{i}$ often result
from the superposition of different effects, some which remove fluid from the
buffer and some which add fluid. If at some time the buffer is empty and the
output rate remains greater than the input rate, fluid does not accumulate and
the buffer remains empty.

Formally, we define the regulator $R(t)=-\min(0,\min_{0\leq s\leq t}X(s))$ and
the fluid queue (the regulated fluid flow) is $\\{Y(t),\varphi(t)\\}$, with

$Y(t)=X(t)+R(t).$ (3)

In the first part of this paper, we focus on characterising the stationary
distribution of the regulated fluid queue, when it exists. Although the
details of our presentation are very much inspired by Ramaswami [44] and da
Silva Soares and Latouche [22], we follow a slightly different path and give
explicit reference to semi-regenerative processes (Çinlar [20, Chapter 10]);
this allows us to interpret in a unified manner the ad-hoc analysis of several
published variants of our basic model.

Let us assume that $Y(0)=0$ and define the sequences $\\{\delta_{n}:n\geq
0\\}$ and $\\{\theta_{n}:n\geq 1\\}$ as follows:

$\displaystyle\delta_{0}$ $\displaystyle=\inf\\{t>0:Y(t)>0\\},$
$\displaystyle\theta_{n}$
$\displaystyle=\inf\\{t>\delta_{n-1}:Y(t)=0\\},\qquad\delta_{n}=\inf\\{t>\theta_{n}:Y(t)>0\\}$

(see an illustration in Figure 2). It is easily seen that $\\{\theta_{n}:n\geq
1\\}$ is a set of regenerative epochs for the process $\\{Y(t),\varphi(t)\\}$:

  * •

$\\{Y(\theta_{n}),\varphi(\theta_{n})\\}$ is a Markov chain on the state space
$\\{0\\}\times\mathcal{S}_{\mbox{\tiny$-$}}$ as the fluid is reaching down to
level zero at these epochs, in a phase of $\mathcal{S}_{\mbox{\tiny$-$}}$,

  * •

the process over the interval $[\theta_{n},\infty)$ is independent of the
process over the interval $[0,\theta_{n})$, given $\varphi(\theta_{n})$, for
all $n$, and

  * •

the distribution of the process over the interval $[\theta_{n},\infty)$, given
$\varphi(\theta_{n})=s$, is the same as the distribution of the process over
the interval $[\theta_{1},\infty)$, given $\varphi(\theta_{1})=s$ for all $n$
and $s\in\mathcal{S}_{\mbox{\tiny$-$}}$.

In consequence, we may immediately write that the stationary distribution
$\boldsymbol{G}(x)=\begin{bmatrix}G_{i}(x):i\in\mathcal{S}\end{bmatrix}$,
defined as222We use boldface letters to represent vectors, and capital letters
for matrices.

$G_{j}(x)=\lim_{t\rightarrow\infty}{\rm P}[Y(t)\leq x,\varphi(t)=j],\qquad
j\in\mathcal{S},$

is given by

$\boldsymbol{G}^{\textnormal{\tiny
T}}(x)=(\boldsymbol{\rho}^{\textnormal{\tiny
T}}\boldsymbol{m})^{-1}\boldsymbol{\rho}^{\textnormal{\tiny T}}M(x)$ (4)

where

  * •

$\boldsymbol{\rho}$ is the stationary probability vector of
$\\{(Y(\theta_{n}),\varphi(\theta_{n}))\\}$,

  * •

$M(x)$ is an $\mathcal{S}_{\mbox{\tiny$-$}}$ by $\mathcal{S}$ matrix, with
components $M_{ij}(x)$ equal to the conditional expected sojourn time of
$(Y(t),\varphi(t))$ in $[0,x]\times\\{j\\}$ during a regeneration interval
$[\theta_{n},\theta_{n+1})$, given that $\varphi(\theta_{n})=i$, and

  * •

$\boldsymbol{m}=M(\infty)\boldsymbol{1}$, where $\boldsymbol{1}$ represents a
vector of ones; the components of $\boldsymbol{m}$ are the conditional
expected lengths of intervals between regeneration points, given the phase at
the beginning of the interval.

We give in Section 2 a few basic characteristics of the process at level 0.
The vector $\boldsymbol{\rho}$ is determined in Section 3, where we define and
analyse two important first passage probability matrices, and $M(x)$ is
determined in Section 4 through the number of crossings of a given level
during regenerative intervals. The results in these two sections are brought
together in Section 5 to give the stationary distribution of the fluid queue.

The key matrices defined in Sections 3 and 4 have very distinct physical
significance but they are algebraically closely related. We show this in
Section 6, using results originally proved in Rogers [46].

In Section 7, we characterise the distribution of first passage times to a
given level, and we analyse in Section 8 the first exit from an interval. The
two sections come as complements to Section 3.

In many applications, in particular in telecommunication modeling, the buffer
level is not allowed to grow without bounds. Often, the evolution of the
process changes as the upper or lower boundary is reached. These, and other
modifications of the basic fluid flow model, are briefly discussed in Section
10 where we show how the regenerative approach may be readily adapted to more
complex assumptions.

One of the nice features of the matrix-analytic approach is that computational
algorithms are easily constructed, following the development of the
theoretical results. As an illustration, we give in Section 9 two of the
simplest, and yet very efficient, algorithms for the numerical computation of
the key matrix $\Psi$ identified in Section 3.

Most of the results presented here have appeared earlier. For that reason, we
give explanatory justifications mostly, and we refer to published sources for
formal justifications. In a few cases, however, we give formal proofs: in
Section 8 we give a new treatment of escape probabilities for null recurrent
processes (Lemma 8.3 and Theorem 8.5), and we offer with Theorem 9.1 a novel
justification for a nice computational procedure for $\Psi$.

##  2 Preliminaries

We partition the generator $Q$ of the Markov process
$\\{\varphi(t):t\in\mathbb{R}^{\mbox{\tiny$+$}}\\}$ in a manner conformant to
the partition (2) of $\mathcal{S}$ and write, possibly after a permutation of
rows and columns,333We write $\boldsymbol{v}_{\mbox{\tiny$+$}}$ for the
subvector $[v_{i}:i\in\mathcal{S}_{\mbox{\tiny$+$}}]$ of any vector
$\boldsymbol{v}$, and $M_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}$ for the
submatrix of any matrix $M$ at the intersection of the rows and columns in
$\mathcal{S}_{\mbox{\tiny$+$}}$. Other sub-vectors and sub-blocks are
similarly defined.

$Q=\begin{bmatrix}Q_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}&Q_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}&Q_{{\mbox{\tiny$+$}}0}\\\
Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}&Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}&Q_{{\mbox{\tiny$-$}}0}\\\
Q_{0{\mbox{\tiny$+$}}}&Q_{0{\mbox{\tiny$-$}}}&Q_{00}\end{bmatrix}.$ (5)

We assume that $Q$ is irreducible. We also define the diagonal matrix $C$ of
fluid rates, $C=\mathrm{diag}(\boldsymbol{c})$, and we partition it as

$C=\begin{bmatrix}C_{\mbox{\tiny$+$}}\\\ &C_{\mbox{\tiny$-$}}\\\
&&0\end{bmatrix}.$

The fluid process $\\{X(t)\\}$ moves up and down in a random manner but its
general direction is determined by the stationary drift
$\mu=\boldsymbol{\alpha}^{\textnormal{\tiny T}}\boldsymbol{c}$, where
$\boldsymbol{\alpha}$ is the stationary probability vector of $Q$:
$\boldsymbol{\alpha}^{\textnormal{\tiny T}}Q=\boldsymbol{0}$,
$\boldsymbol{\alpha}^{\textnormal{\tiny T}}\boldsymbol{1}=1$. If $\mu>0$, the
process eventually drift to $+\infty$, that is,
$\lim_{t\rightarrow\infty}X(t)=\infty$, if $\mu<0$, then
$\lim_{t\rightarrow\infty}X(t)=-\infty$; in both cases the process is
transient. If $\mu=0$, the process is null-recurrent and $\lim\sup
X(t)=\infty$ while $\lim\inf X(t)=-\infty$ (Asmussen [4, Page 314, Proposition
2.10]).

Things are slightly different for the regulated process $\\{Y(t)\\}$. If
$\mu<0$, the process repeatedly alternates between intervals of time where
$Y(t)>0$ and intervals where $Y(t)=0$. As we find in Section 6, the length of
each cycle has finite expectation, the regulated process is positive
recurrent, and we may determine its stationary distribution. If $\mu>0$, then
$Y(t)$ might not return to level 0 and so the process is transient. Not
surprisingly, it is null-recurrent if $\mu=0$.

We define $\bar{\varphi}_{n}=\varphi(\theta_{n})$, $n\geq 1$. The process
$\\{\bar{\varphi}_{n}\\}$ of the phases visited at epochs of regeneration is a
Markov chain on the state space $\mathcal{S}_{\mbox{\tiny$-$}}$ by definition,
and its transition matrix is $H$, with

$H_{ij}={\rm
P}[\theta_{n+1}-\theta_{n}<\infty,\bar{\varphi}_{n+1}=j|\bar{\varphi}_{n}=i],\qquad\mbox{$i$,
$j$ in $\mathcal{S}_{\mbox{\tiny$-$}}$},.$

To determine $H$, we split the interval $(\theta_{n},\theta_{n+1}]$ in two and
we condition on the phase occupied at time $\delta_{n}$. Thus, $H$ is the
product

$H=\Phi\,\Psi,$ (6)

where

$\displaystyle\Phi_{ik}$ $\displaystyle={\rm
P}[\delta_{n}-\theta_{n}<\infty,\varphi(\delta_{n})=k|\bar{\varphi}_{n}=i],$
$\displaystyle
i\in\mathcal{S}_{\mbox{\tiny$-$}},k\in\mathcal{S}_{\mbox{\tiny$+$}},$
$\displaystyle\Psi_{kj}$ $\displaystyle={\rm
P}[\theta_{n+1}-\delta_{n}<\infty,\bar{\varphi}_{n+1}=j|\varphi(\delta_{n})=k],$
$\displaystyle
k\in\mathcal{S}_{\mbox{\tiny$+$}},j\in\mathcal{S}_{\mbox{\tiny$-$}}.$ (7)

The matrix $\Phi$ is easily determined: $\varphi(t)$ remains in
$\mathcal{S}_{\mbox{\tiny$-$}}\cup\mathcal{S}_{0}$ during the interval
$(\theta_{n},\delta_{n})$, and thus

$\Phi=\begin{bmatrix}I&0\end{bmatrix}\left(-\begin{bmatrix}Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}&Q_{{\mbox{\tiny$-$}}0}\\\
Q_{0{\mbox{\tiny$-$}}}&Q_{00}\end{bmatrix}^{-1}\right)\begin{bmatrix}Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}\\\
Q_{0{\mbox{\tiny$+$}}}\end{bmatrix}$ (8)

(Latouche and Ramaswami [41, Section 5.5]). To determine the matrix $\Psi$
requires more effort, and we devote Section 3 to the determination of its
characteristic equation. Once $\Psi$ is known, the vector $\boldsymbol{\rho}$
in (4) is determined by the system $\boldsymbol{\rho}^{\textnormal{\tiny
T}}H=\boldsymbol{\rho}^{\textnormal{\tiny T}}$,
$\boldsymbol{\rho}^{\textnormal{\tiny T}}\boldsymbol{1}=1$.

In a similar manner, we decompose $M(x)$ as $M(x)=M(0)+M(0,x]$, where $M(0)$
is the expected time spent at level 0 between the epochs $\theta_{n}$ and
$\delta_{n}$, and $M(0,x]$ is the expected time spent in the semi-open
interval $(0,x]$. We decompose $M(0)$ as

$M(0)=\begin{bmatrix}M(0)_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}&M(0)_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}&M(0)_{{\mbox{\tiny$-$}}0}\end{bmatrix},$

and immediately note that $M(0)_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}=0$ as the
fluid queue does not spend any time at level 0 in a phase of
$\mathcal{S}_{\mbox{\tiny$+$}}$. Furthermore,

$\begin{bmatrix}M(0)_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}&M(0)_{{\mbox{\tiny$-$}}0}\end{bmatrix}=\begin{bmatrix}I&0\end{bmatrix}\left(-\begin{bmatrix}Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}&Q_{{\mbox{\tiny$-$}}0}\\\
Q_{0{\mbox{\tiny$-$}}}&Q_{00}\end{bmatrix}^{-1}\right),$ (9)

see [41, Section 5.5]. The second term in $M(x)$ is equal to

$M(0,x]=\Phi\,\widetilde{M}(x),$ (10)

where $\widetilde{M}(x)$ is the matrix of expected times spent in $(0,x]$
during the interval of time $(\delta_{n},\theta_{n+1}]$. It is determined in
Section 4.

##  3 First passage probabilities

We deal in this section with the fluid flow model $\\{X(t),\varphi(t)\\}$
without boundary. Its transition structure is independent of the level and so
we shall not always pay close attention to the exact value of $X$, but be more
interested in differences of level. For instance, the matrix $\Psi$ defined in
(7) might have been defined as

$\Psi_{ij}={\rm P}[\Theta<\infty,\varphi(\Theta)=j|\varphi(0)=i],\qquad
i\in\mathcal{S}_{\mbox{\tiny$+$}},j\in\mathcal{S}_{\mbox{\tiny$-$}}.$

independently of $X(0)$, where $\Theta=\inf\\{t>0:X(t)=X(0)\\}$ is the first
return time to the initial level, starting in a phase of
$\mathcal{S}_{\mbox{\tiny$+$}}$.

Furthermore, let

$\tau^{\mbox{\tiny$-$}}_{x}=\inf\\{t:X(t)<X(0)-x\\}$ (11)

be the first passage time to level $X(0)-x$, for $x>0$, and denote by
$\gamma^{\mbox{\tiny$-$}}(x)=\varphi(\tau^{\mbox{\tiny$-$}}_{x})$ the value of
the phase when $X(0)-x$ is reached for the first time.444As the trajectory of
$X(t)$ is continuous, we might have defined
$\tau^{\mbox{\tiny$-$}}_{x}=\inf\\{t:X(t)=X(0)-x\\}$, but the strict
inequality in (11) will be useful in Section 8.  For $\varphi(0)$ in
$\mathcal{S}_{\mbox{\tiny$-$}}$, the process
$\\{\gamma^{\mbox{\tiny$-$}}(x):x\geq 0\\}$ is a continuous-parameter Markov
process on the state space $\mathcal{S}_{\mbox{\tiny$-$}}$, and there exists a
matrix $U$ such that

$(e^{Ux})_{ij}={\rm
P}[\tau_{x}^{\mbox{\tiny$-$}}<\infty,\gamma^{\mbox{\tiny$-$}}(x)=j|\gamma^{\mbox{\tiny$-$}}(0)=i],\qquad
i,j\in\mathcal{S}_{\mbox{\tiny$-$}}.$

As $Q$ is irreducible, $\Psi$ is strictly positive, meaning that $\Psi_{ij}>0$
for all $i$ in $\mathcal{S}_{\mbox{\tiny$+$}}$, $j$ in
$\mathcal{S}_{\mbox{\tiny$-$}}$. To see this, imagine a trajectory of positive
probability such that, starting from $(0,i)$, the process returns at 0 for the
first time in phase $j$, in finite time. A formal proof is in Govorun et al.
[30, Lemma 4.3], or Guo [33, Theorem 5]. In consequence, the off-diagonal
elements of the generator $U$ are all strictly positive. We discuss at greater
length the algebraic properties of $U$ and $\Psi$ in Section 6, but mention
here the most important ones, in relation to the stationary drift $\mu$:

  * •

if $\mu\leq 0$, then $\Psi\boldsymbol{1}=\boldsymbol{1}$ and
$U\boldsymbol{1}=\boldsymbol{0}$, that is, $e^{Ux}$ is stochastic,

  * •

if $\mu>0$, then $\Psi\boldsymbol{1}<\boldsymbol{1}$ and
$U\boldsymbol{1}<\boldsymbol{0}$, that is, $e^{Ux}$ is substochastic.555With
$\boldsymbol{u}$ and $\boldsymbol{v}$ two vectors on the same set of indices,
we write $\boldsymbol{u}<\boldsymbol{v}$ if $u_{i}<v_{i}$ for all index $i$.

Finally, we define the matrix

$T=\begin{bmatrix}Q_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}&Q_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}\\\
Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}&Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}\end{bmatrix}+\begin{bmatrix}Q_{{\mbox{\tiny$+$}}0}\\\
Q_{{\mbox{\tiny$-$}}0}\end{bmatrix}(-Q_{00})^{-1}\begin{bmatrix}Q_{0{\mbox{\tiny$+$}}}&Q_{0{\mbox{\tiny$-$}}}\end{bmatrix}$
(12)

indexed by the states in
$\mathcal{S}_{\mbox{\tiny$+$}}\cup\mathcal{S}_{\mbox{\tiny$-$}}$ and we
partition it as

$T=\begin{bmatrix}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}&T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}\\\
T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}&T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}\end{bmatrix}.$

This is the generator of the censored process $\\{\varphi(t)\\}$ observed only
during the intervals of time spent in
$\mathcal{S}_{\mbox{\tiny$+$}}\cup\mathcal{S}_{\mbox{\tiny$-$}}$.

###### Theorem 3.1

The matrix $U$ is given by

$U=|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}+|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}\Psi,$
(13)

and $\Psi$, as defined in (7), is the minimal nonnegative solution of the
quadratic Riccati equation

$C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}+C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}\Psi+\Psi|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}+\Psi|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}\Psi=0,$
(14)

where $|C_{\mbox{\tiny$-$}}|$ is the matrix of absolute values of the elements
of $C_{\mbox{\tiny$-$}}$. $\square$

We give here a high-level justification, a formal proof is in Ahn and
Ramaswami [1] and da Silva Soares and Latouche [24]. To determine $\Psi$, it
is simpler to use the censored process on
$\mathcal{S}_{\mbox{\tiny$+$}}\cup\mathcal{S}_{\mbox{\tiny$-$}}$ since there
is no change in the level while $\varphi$ is in $\mathcal{S}_{0}$; that is why
we use in (14) the generator $T$ of (12) instead of the generator $Q$ of (5).

Furthermore, instead of tracking the evolution of the phase process in time as
one might be tempted to do, we track its evolution over changes of the level.
The parameters $|c_{i}|$, for
$i\in\mathcal{S}_{\mbox{\tiny$+$}}\cup\mathcal{S}_{\mbox{\tiny$-$}}$, are
conversion rates of time to fluid and their reciprocal $|c_{i}|^{-1}$ are
conversion rates of fluid to time, so that the matrix $|C|^{-1}T$ indicates
how the phase evolves as the fluid level is increasing or decreasing.

With this in mind, we write

$\Psi=\int_{0}^{\infty}e^{C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}y}\,C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}\,e^{Uy}\,\mathrm{d}y.$
(15)

The justification goes as follows (see the illustration in Figure 3). The
process starts in a phase of $\mathcal{S}_{\mbox{\tiny$+$}}$, and we assume
without loss of generality that $X(0)=0$. The factor
$e^{C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}y}$ in the
right-hand side of (15) is the probability that the phase remains in
$\mathcal{S}_{\mbox{\tiny$+$}}$ until the fluid has increased up to level $y$.
The factor
$C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}\mathrm{d}y$ is
the probability that between the levels $y$ and $y+\mathrm{d}y$, the phase
moves to $\mathcal{S}_{\mbox{\tiny$-$}}$ and starts to decrease. The factor
$e^{Uy}$ is the probability that the fluid eventually goes down to level 0.
Level $y$ is reached at some unspecified moment $t(y)$, and the process moves
without constraint between $t(y)$ and $\Theta$.

$t$$\Theta$$t(y)$$X(t)$$y$$0$

Figure 3: The fluid increases and reaches $y$ at some unspecified time $t(y)$.
Between $t(y)$ and $\Theta$ the fluid moves without constraint.

We pre-multiply both sides of (15) by
$C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}$ integrate by
part, and find

$C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}+C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}\Psi+\Psi
U=0.$ (16)

This is a nonsingular Sylvester equation (Lancaster and Tismenetsky [37]) and
we may characterise $\Psi$ as the unique solution of (16) if $U$ is known.

To prove (13), we write

$\displaystyle U_{ij}h$ $\displaystyle={\rm
P}[\gamma^{\mbox{\tiny$-$}}(u+h)=j|\gamma^{\mbox{\tiny$-$}}(u)=i]+o(h),$
$\displaystyle=(|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}})_{ij}h+(|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}\Psi)_{ij}h+o(h)$
(17)

for $i$, $j$ in $\mathcal{S}_{\mbox{\tiny$-$}}$. The first term is the
probability that $\varphi$ changes from $i$ to $j$, the second term is the
probability that $\varphi$ changes from $i$ to a phase in
$\mathcal{S}_{\mbox{\tiny$+$}}$ at some level $u$ and is in phase $j$ when the
process later returns to level $u$. We illustrate this in Figure 4, where we
plot the trajectory of the process from level $y$ to level 0, and draw with
solid lines the part that corresponds to $\gamma^{\mbox{\tiny$-$}}(x)$. At
level $y_{a}$ there is a simple change from phase 2 to phase 1; at level
$y_{b}$ there is a change from phase 2 to phase 4 in
$\mathcal{S}_{\mbox{\tiny$+$}}$ and upon return to level $y_{b}$ the process
is in phase 1.

$t$$\tau^{\mbox{\tiny$-$}}_{y}$$X(t)$f $y$$y_{b}$$y_{a}$

Figure 4: The sample path of $X(t)$ from $y$ to 0 is reproduced from Figure 3,
the path of $\gamma^{\mbox{\tiny$-$}}(x)$ corresponds to the solid downward
segments, $y_{a}$ and $y_{b}$ are two of the levels where the value of
$\gamma^{\mbox{\tiny$-$}}(x)$ changes.

The Riccati equation (14) is obtained by replacing in (16) $U$ by its
expression in (13). Over the years, very efficient algorithms have been
developed to solve (14), we describe two of these in Section 9.

Two other matrices, $\widehat{U}$ and $\widehat{\Psi}$, may be defined at this
stage. Let $\tau^{\mbox{\tiny$+$}}_{x}=\inf\\{t:X(t)>X(0)+x\\}$ be the first
passage time to level $X(0)+x$, for $x>0$, and denote by
$\gamma^{\mbox{\tiny$+$}}(x)=\varphi(\tau^{\mbox{\tiny$+$}}_{x})$ the value of
the phase when level $X(0)+x$ is reached for the first time. The matrix
$\widehat{U}$ is the generator of $\\{\gamma^{\mbox{\tiny$+$}}(x)\\}$ and
$\widehat{\Psi}$ is the matrix of first passage probability back to the
initial level, given that the phase at time 0 is in
$\mathcal{S}_{\mbox{\tiny$-$}}$:

$\widehat{\Psi}_{ij}={\rm
P}[\Theta<\infty,\varphi(\Theta)=j|\varphi(0)=i],\qquad
i\in\mathcal{S}_{\mbox{\tiny$-$}},j\in\mathcal{S}_{\mbox{\tiny$+$}}.$

We easily adapt the argument in Theorem 3.1 to prove the following.

###### Corollary 3.2

The matrix $\widehat{U}$ is given by

$\widehat{U}=C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}+C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}\widehat{\Psi}$
(18)

and $\widehat{\Psi}$ is the minimal nonnegative solution of the equation

$|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}+|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}\widehat{\Psi}+\widehat{\Psi}C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}+\widehat{\Psi}C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}\widehat{\Psi}=0.$
(19)

Furthermore, $\widehat{\Psi}\boldsymbol{1}<\boldsymbol{1}$ and
$\widehat{U}\boldsymbol{1}<\boldsymbol{0}$ if $\mu<0$, while
$\widehat{\Psi}\boldsymbol{1}=\boldsymbol{1}$ and
$\widehat{U}\boldsymbol{1}=\boldsymbol{0}$ if $\mu\geq 0$. $\square$

##  4 Number of crossings

As a preliminary step to determining the expected sojourn time
$\widetilde{M}(x)$ in $(0,x]$ during an interval $(\delta_{n},\theta_{n+1}]$,
we analyse the number of times the fluid crosses a given level during a
regenerative interval. We define $N_{ij}(x)$ to be the expected number of
times $(X,\varphi)=(X(0)+x,j)$ during the interval $(0,\Theta]$, given that
$\varphi(0)=i$, for $i$ and $j$ in $\mathcal{S}_{\mbox{\tiny$+$}}$.

###### Theorem 4.1

The matrix $N(x)$ is given by $N(x)=e^{Kx}$, where

$K=C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}+\Psi|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}$
(20)

is a matrix indexed by $\mathcal{S}_{\mbox{\tiny$+$}}$. $\square$

The formal proof is given in Ramaswami [44] and proceeds as follows.666The
argument is similar to the one used in Latouche and Ramaswami [41, Theorem
6.2.7] in the context of Quasi-Birth-and-Death processes. Assume without loss
of generality that $X(0)=0$, take $x$ and $y>0$, and count the expected number
of visits to $(x+y,j)$, starting from $(0,i)$, before the first return to
level 0. We group the visits to $(x+y,j)$ into subintervals between successive
up-crossings of level $x$, and write that

$N_{ij}(x+y)=\sum_{k\in\mathcal{S}_{\mbox{\tiny$+$}}}N_{ik}(x)N_{kj}(x,x+y),\qquad
i,j\in\mathcal{S}_{\mbox{\tiny$+$}},$

where $N_{kj}(x,x+y)$ is the expected number of visits to $(x+y,j)$ between
two successive visits to level $x$, starting from $(x,k)$. Remember that it is
the distance $y$ between the target $x+y$ and the starting level $x$ matters,
not the specific location of the latter. Thus, $N(x,x+y)=N(0,y)$ and so the
equation above is also written as $N(x+y)=N(x)N(y)$. From the semi-group
property, we conclude that there exists a matrix $K$ such that $N(x)=e^{Kx}$.

Next, we approximate $(e^{Kh})_{ij}$ for $h$ small and $i\not=j$ as

$\displaystyle(e^{Kh})_{ij}$ $\displaystyle=K_{ij}h+o(h)$
$\displaystyle=c_{i}^{-1}Q_{ij}h+\sum_{k\in\mathcal{S}_{0}}c_{i}^{-1}Q_{ik}((-Q_{00})^{-1}Q_{0{\mbox{\tiny$+$}}})_{kj}h$
$\displaystyle\quad+(1+c_{i}^{-1}Q_{ii}h)\sum_{\ell\in\mathcal{S}_{\mbox{\tiny$-$}}}\Psi_{i\ell}|c_{\ell}|^{-1}Q_{\ell
j}h$
$\displaystyle\quad+(1+c_{i}^{-1}Q_{ii}h)\sum_{\ell\in\mathcal{S}_{\mbox{\tiny$-$}}}\Psi_{i\ell}|c_{\ell}|^{-1}\sum_{k\in\mathcal{S}_{0}}|c_{\ell}|^{-1}Q_{\ell
k}((-Q_{00})^{-1}Q_{0{\mbox{\tiny$+$}}})_{kj}h+o(h).$

The first two terms are about the process being in phase $j$ at the first
crossing of level $h$;

  * •

the first is the probability that during the interval of time $c_{i}^{-1}h$
the phase process changes from $i$ to $j$ and is still in $j$ when the fluid
crosses level $h$;

  * •

the second term is the probability that during that interval of time, the
phase switches to $k$ in $\mathcal{S}_{0}$ at some unspecified level below
$h$, remains at that level until there is a jump to $j$, and is still in $j$
when the fluid crosses level $h$.

The next two terms cover the circumstances where the process does not leave
phase $i$ before crossing level $h$, returns to level $h$ in some phase $\ell$
in $\mathcal{S}_{\mbox{\tiny$-$}}$, switches to phase $j$ during the interval
of time $|c_{\ell}|^{-1}h$ and eventually crosses level $h$ in phase $j$. The
$o(h)$ term captures the negligible probabilities of crossing level $h$ more
than once. Simple manipulations give us

$\displaystyle(e^{Kh})_{ij}$
$\displaystyle=(c_{i}^{-1}T_{ij}+\sum_{\ell\in\mathcal{S}_{\mbox{\tiny$-$}}}\Psi_{i\ell}|c_{\ell}|^{-1}T_{\ell
j})h+o(h),\qquad i\not=j\in\mathcal{S}_{\mbox{\tiny$+$}}.$

A similar argument holds for $i=j$ and so is (20) justified.

We partition the matrix $\widetilde{M}(x)$ in three blocks, corresponding to
the three subsets $\mathcal{S}_{\mbox{\tiny$+$}}$,
$\mathcal{S}_{\mbox{\tiny$-$}}$ and $\mathcal{S}_{0}$ of phases, and we deal
with each one separately.

###### Theorem 4.2

The matrix of expected sojourn times in $(0,x]$ during a regenerative interval
is

$\widetilde{M}(x)=\begin{bmatrix}\widetilde{M}_{\mbox{\tiny$+$}}(x)&\widetilde{M}_{\mbox{\tiny$-$}}(x)&\widetilde{M}_{0}(x)\end{bmatrix}$

with

$\displaystyle\widetilde{M}_{\mbox{\tiny$+$}}(x)$
$\displaystyle=\mathcal{F}(K,x)C_{\mbox{\tiny$+$}}^{-1},$ (21)
$\displaystyle\widetilde{M}_{\mbox{\tiny$-$}}(x)$
$\displaystyle=\mathcal{F}(K,x)\Psi|C_{\mbox{\tiny$-$}}|^{-1},$ (22)
$\displaystyle\widetilde{M}_{0}(x)$
$\displaystyle=\mathcal{F}(K,x)\begin{bmatrix}C_{\mbox{\tiny$+$}}^{-1}&\Psi|C_{\mbox{\tiny$-$}}|^{-1}\end{bmatrix}\begin{bmatrix}Q_{{\mbox{\tiny$+$}}0}\\\
Q_{{\mbox{\tiny$-$}}0}\end{bmatrix}(-Q_{00})^{-1}.$ (23)

where $\mathcal{F}(K,x)=\int_{0}^{x}e^{Ku}\,\mathrm{d}u$.

Proof  A formal proof of (21, 22) is given in Latouche and Nguyen [38], we
give here a heuristic argument, and we give a justification for (23).

For $i$, $j$ in $\mathcal{S}_{\mbox{\tiny$+$}}$, we argue that

$\widetilde{M}_{ij}(x)=\int_{0}^{x}N_{ij}(u)c_{j}^{-1}\,\mathrm{d}u.$

To see this, we interpret $c_{j}^{-1}\mathrm{d}u$ as the expected time spent
by $X$ in $(u,u+\mathrm{d}u)$ each time there is a visit to $(u,j)$, and
multiply this by the expected number $N_{ij}(u)$ of such visits, so that
$N_{ij}(u)c_{j}^{-1}\,\mathrm{d}u$ is the total time spent by the process.
There only remains for us to integrate over the interval $(0,x)$ and use
Theorem 4.1.

Next, we define the matrix $N^{\prime}(x)$ of expected number of down-
crossings of level $x$ in a phase of $\mathcal{S}_{\mbox{\tiny$-$}}$. As
$X(0)=0<x$, of necessity each down-crossing of level $x$ is preceded by an up-
crossing of that same level. So, if we condition on the phase at the up-
crossing, we find that

$N^{\prime}_{ij}(x)=\sum_{k\in\mathcal{S}_{\mbox{\tiny$+$}}}N_{ik}(x)\Psi_{kj},\qquad
i\in\mathcal{S}_{\mbox{\tiny$+$}},j\in\mathcal{S}_{\mbox{\tiny$-$}}$

and (22) is proved by repeating the argument of (21).

Finally, let us use the expression excursion to $\mathcal{S}_{0}$ for
intervals of time spent in $(0,x)\times\mathcal{S}_{0}$ and separated by
visits to
$(0,x)\times\mathcal{S}_{\mbox{\tiny$+$}}\cup\mathcal{S}_{\mbox{\tiny$-$}}$.
Equation (23) expresses $\widetilde{M}_{0}(x)$ as the product of the expected
number of excursions by the time spent in individual phases during excursion.
Indeed, $Q_{{\mbox{\tiny$+$}}0}$ and $Q_{{\mbox{\tiny$-$}}0}$ are the matrices
of transition rates from $\mathcal{S}_{\mbox{\tiny$+$}}$ and
$\mathcal{S}_{\mbox{\tiny$-$}}$ to $\mathcal{S}_{0}$, and
$\mathcal{F}(K,x)\begin{bmatrix}C_{\mbox{\tiny$+$}}^{-1}&\Psi|C_{\mbox{\tiny$-$}}|^{-1}\end{bmatrix}$
is the matrix of expected time spent in the states of
$(0,x)\times\mathcal{S}_{\mbox{\tiny$+$}}\cup\mathcal{S}_{\mbox{\tiny$-$}}$.
Thus, the component $N^{\prime\prime}_{ik}(x)$ of the product

$N^{\prime\prime}(x)=\mathcal{F}(K,x)\begin{bmatrix}C_{\mbox{\tiny$+$}}^{-1}&\Psi|C_{\mbox{\tiny$-$}}|^{-1}\end{bmatrix}\begin{bmatrix}Q_{{\mbox{\tiny$+$}}0}\\\
Q_{{\mbox{\tiny$-$}}0}\end{bmatrix}$

is the expected number of excursions that start in phase $k$. As
$(-Q_{00}^{-1})_{kj}$ is the expected time spent in $j$ during an excursion
which starts in $k$, this proves (23). $\square$

As we shall see in Section 6, all the eigenvalues of $K$ are in
$\mathbb{C}_{\leq 0}$, that is, they have a negative real part. If $\mu<0$,
the real parts are all strictly negative, otherwise $K$ has one eigenvalue
equal to 0, the others being in $\mathbb{C}_{<0}$. With this, we may express
the integral $\int_{0}^{x}e^{Ku}\,\mathrm{d}u$ as follows.

  * •

If $\mu<0$, then

$\mathcal{F}(K,x)=(-K^{-1})(I-e^{Kx}),$ (24)

  * •

if $\mu\geq 0$, then

$\mathcal{F}(K,x)=(-K^{\\#})(I-e^{Kx})+x\boldsymbol{v}\boldsymbol{u}^{\textnormal{\tiny
T}},$

where $\boldsymbol{v}$ and $\boldsymbol{u}$ are the normalised right- and
left-eigenvectors of $K$ for the eigenvalue 0, and $K^{\\#}$ is the group
inverse777The group inverse $K^{\\#}$ of $K$ is the unique matrix such that
$K^{\\#}K=I-\boldsymbol{v}\boldsymbol{u}^{\textnormal{\tiny T}}$,
$K^{\\#}\boldsymbol{v}=\boldsymbol{0}$. See Campbell and Meyer [19, Chapter
7]. of $K$.

In the same manner as we defined $\widehat{U}$ and $\widehat{\Psi}$, we may
define $N^{*}(x)$ as the number of crossings of level $X(0)-x$ in a phase of
$\mathcal{S}_{\mbox{\tiny$-$}}$, starting from a phase in
$\mathcal{S}_{\mbox{\tiny$-$}}$, before the first return to level $X(0)$, and
we have the following corollary to Theorem 4.1:

###### Corollary 4.3

The matrix $N^{*}(x)$ is given by $N^{*}(x)=e^{\widehat{K}x}$, where

$\widehat{K}=|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}+\widehat{\Psi}C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}$
(25)

is indexed by $\mathcal{S}_{\mbox{\tiny$-$}}$. Furthermore, the eigenvalues of
$\widehat{K}$ are in $\mathbb{C}_{<0}$, with the exception of one eigenvalue
equal to 0 if $\mu\leq 0$. $\square$

##  5 Stationary distribution

We collect the results obtained in the preceding sections and we express as
follows the stationary distribution of the fluid queue, when it exists.

###### Theorem 5.1

If $\mu<0$, the regulated process $(Y(t),\varphi(t))$ has a stationary
distribution, given by

$\boldsymbol{G}^{\textnormal{\tiny
T}}(x)=c(\begin{bmatrix}0&\boldsymbol{p}^{\textnormal{\tiny
T}}_{\mbox{\tiny$-$}}&\boldsymbol{p}^{\textnormal{\tiny
T}}_{0}\end{bmatrix}+(\boldsymbol{p}^{\textnormal{\tiny
T}}_{\mbox{\tiny$-$}}Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}+\boldsymbol{p}^{\textnormal{\tiny
T}}_{0}Q_{0{\mbox{\tiny$+$}}})\widetilde{M}(x)),$

where

$\displaystyle\begin{bmatrix}\boldsymbol{p}^{\textnormal{\tiny
T}}_{\mbox{\tiny$-$}}&\boldsymbol{p}^{\textnormal{\tiny T}}_{0}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\boldsymbol{\rho}^{\textnormal{\tiny
T}}&0\end{bmatrix}\left(-\begin{bmatrix}Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}&Q_{{\mbox{\tiny$-$}}0}\\\
Q_{0{\mbox{\tiny$-$}}}&Q_{00}\end{bmatrix}^{-1}\right),$
$\displaystyle\widetilde{M}(x)$
$\displaystyle=(-K)^{-1}(I-e^{Kx})\begin{bmatrix}C_{\mbox{\tiny$+$}}^{-1}&\Psi|C_{\mbox{\tiny$-$}}|^{-1}&\Gamma\end{bmatrix},$
(26) $\displaystyle\Gamma$
$\displaystyle=(C_{\mbox{\tiny$+$}}^{-1}Q_{{\mbox{\tiny$+$}}0}+\Psi|C_{\mbox{\tiny$-$}}|^{-1}Q_{{\mbox{\tiny$-$}}0})(-Q_{00})^{-1},$

and $c$ is the normalising constant, with

$c^{-1}=\begin{bmatrix}\boldsymbol{p}^{\textnormal{\tiny
T}}_{\mbox{\tiny$-$}}&\boldsymbol{p}^{\textnormal{\tiny
T}}_{0}\end{bmatrix}(\boldsymbol{1}+\begin{bmatrix}Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}\\\
Q_{0{\mbox{\tiny$+$}}}\end{bmatrix}(-K)^{-1}(C_{\mbox{\tiny$+$}}^{-1}\boldsymbol{1}+\Psi|C_{\mbox{\tiny$-$}}|^{-1}\boldsymbol{1}+\Gamma\boldsymbol{1})).$

Proof The expression (26) for $\widetilde{M}(x)$ results from Theorem 4.2 and
Equation (24). Furthermore, the vector $\boldsymbol{m}$ of conditional
expected duration of of a regeneration interval, given the initial phase, is
given by

$\displaystyle\boldsymbol{m}$
$\displaystyle=\lim_{x\rightarrow\infty}M(x)\boldsymbol{1}$
$\displaystyle=M(0)\boldsymbol{1}+\lim_{x\rightarrow\infty}\Phi\widetilde{M}(x)\boldsymbol{1}\qquad\mbox{by
(\ref{e:mo}, \ref{e:phim}),}$
$\displaystyle=M(0)\boldsymbol{1}+\Phi(-K)^{-1}(C_{\mbox{\tiny$+$}}^{-1}\boldsymbol{1}+\Psi|C_{\mbox{\tiny$-$}}|^{-1}\boldsymbol{1}+\Gamma\boldsymbol{1})$
(27)

by (26). The normalising constant $c$ equals
$(\boldsymbol{\rho}^{\textnormal{\tiny T}}\boldsymbol{m})^{-1}$ by (4); simple
calculations complete the proof. $\square$

The stationary distribution appears under various forms in the literature
(Asmussen [3], Govorun et al. [30], Ramaswami [44], Rogers [46]); the matrix
$\widetilde{M}(x)$ is a common feature in papers that rely on matrix-analytic
methods, the vector of probability mass at zero and the left-factor of
$\widetilde{M}(x)$ are given very different expressions, depending on the
specific approach followed by the authors.

##  6 Wiener-Hopf factorisation

The matrices $U$, $\Psi$ and $K$ defined in Sections 3 and 4 are related in
many ways. Obviously, $U$ and $K$ are determined by $\Psi$ through (13) and
(20), respectively, but $\Psi$ may be seen as a function of $U$ by (16) and we
might combine (14) and (20) to write

$C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}+K\Psi+\Psi|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}=0,$

from which we conclude that $\Psi$ is a function of $K$. There exist also
numerous relations between the 3-tuples $(U,\Psi,K)$ and
$(\widehat{U},\widehat{\Psi},\widehat{K})$, as we briefly discuss below.

The starting point is the Wiener-Hopf factorisation

$C^{-1}T\begin{bmatrix}I&\Psi\\\
\widehat{\Psi}&I\end{bmatrix}=\begin{bmatrix}I&\Psi\\\
\widehat{\Psi}&I\end{bmatrix}\begin{bmatrix}\widehat{U}&0\\\
0&-U\end{bmatrix}$ (28)

proved in Rogers [46], with

$C=\begin{bmatrix}C_{\mbox{\tiny$+$}}\\\ &C_{\mbox{\tiny$-$}}\end{bmatrix}.$

Equation (28) may be proved by direct verification, starting from the Riccati
equations (14) and (19), and the expressions (13) and (18) for $U$ and
$\widehat{U}$. The equation below may also be proved by direct verification:

$\begin{bmatrix}I&-\Psi\\\
-\widehat{\Psi}&I\end{bmatrix}C^{-1}T.=\begin{bmatrix}K&0\\\
0&-\widehat{K}\end{bmatrix}\begin{bmatrix}I&-\Psi\\\
-\widehat{\Psi}&I\end{bmatrix},$ (29)

An immediate consequence is that $U$ and $\widehat{U}$ on the one hand, $K$
and $\widehat{K}$ on the other hand, share the eigenvalues of $C^{-1}T$. We
denote by $m_{\mbox{\tiny$+$}}$ and $m_{\mbox{\tiny$-$}}$ respectively the
number of phases in $\mathcal{S}_{\mbox{\tiny$+$}}$ and
$\mathcal{S}_{\mbox{\tiny$-$}}$, and we denote by $\\{\lambda_{i}:1\leq i\leq
m_{\mbox{\tiny$+$}}+m_{\mbox{\tiny$-$}}\\}$ the eigenvalues of $C^{-1}T$,
labeled in increasing value of their real part.

###### Theorem 6.1

The eigenvalues $\lambda_{{m_{\mbox{\tiny$+$}}}}$ and
$\lambda_{{m_{\mbox{\tiny$+$}}}+1}$ are real, and are distinct from the
others:

$\mathrm{Re}(\lambda_{1})\leq\cdots\leq\mathrm{Re}(\lambda_{{m_{\mbox{\tiny$+$}}}-1})<\lambda_{m_{\mbox{\tiny$+$}}}\leq\lambda_{{m_{\mbox{\tiny$+$}}}+1}<\mathrm{Re}(\lambda_{{m_{\mbox{\tiny$+$}}}+2}).\leq\cdots\leq\mathrm{Re}(\lambda_{m}).$
(30)

Furthermore, $\lambda_{i}$, $1\leq i\leq m_{\mbox{\tiny$+$}}$ are the
eigenvalues of $K$ and also the eigenvalues of $\widehat{U}$, while
$\lambda_{{m_{\mbox{\tiny$+$}}}+1}$,
…$\lambda_{m_{\mbox{\tiny$+$}}+m_{\mbox{\tiny$-$}}}$ are those of $-U$ and of
$-\widehat{K}$. Finally,

  * •

if $\mu<0$, then
$\lambda_{m_{\mbox{\tiny$+$}}}<0=\lambda_{{m_{\mbox{\tiny$+$}}}+1}$,

  * •

if $\mu=0$, then
$\lambda_{{m_{\mbox{\tiny$+$}}}}=0=\lambda_{{m_{\mbox{\tiny$+$}}}+1}$,

  * •

if $\mu>0$, then
$\lambda_{{m_{\mbox{\tiny$+$}}}}=0<\lambda_{{m_{\mbox{\tiny$+$}}}+1}$.

Details are given in [30]. The main properties are summarised in Table 1.

drift $\mu$ | $\,<0$ |  0 | $>0$
---|---|---|---
$\Psi\boldsymbol{1}$ | $\boldsymbol{1}$ | $\boldsymbol{1}$ | $<\boldsymbol{1}$
$\widehat{\Psi}\boldsymbol{1}$ | $<\boldsymbol{1}$ | $\boldsymbol{1}$ | $\boldsymbol{1}$
$U$, $\widehat{K}$ | 0 | 0 | $<0$
$\widehat{U}$, $K$ | $<0$ | 0 | 0
Table 1: Crossing probabilities, and maximal eigenvalues of $U$, $K$,
$\widehat{U}$, and $\widehat{K}$ in relation to the asymptotic drift $\mu$.

Actually, $K$ and $\widehat{U}$ are similar matrices, and $\widehat{K}$ is
similar to $U$. This is easy to show if the stationary drift is different from
0 for, we find from (28, 29) after some simple algebraic manipulation that

$K(I-\Psi\widehat{\Psi})=(I-\Psi\widehat{\Psi})\widehat{U}\quad\mbox{and}\quad\widehat{K}(I-\widehat{\Psi}\Psi)=(I-\widehat{\Psi}\Psi)U.$
(31)

If $\mu\not=0$, then both $\Psi\widehat{\Psi}$ and $\widehat{\Psi}\Psi$ are
strictly sub-stochastic matrices, $I-\Psi\widehat{\Psi}$ and
$I-\widehat{\Psi}\Psi$ are non-singular, and we may write

$\displaystyle\widehat{U}=(I-\Psi\widehat{\Psi})^{-1}K(I-\Psi\widehat{\Psi})$
$\displaystyle\widehat{K}=(I-\widehat{\Psi}\Psi)U(I-\widehat{\Psi}\Psi)^{-1}$

which shows that $K$ is similar to $\widehat{U}$ and $\widehat{K}$ is similar
to $U$.

If $\mu=0$, however, $\Psi\widehat{\Psi}\boldsymbol{1}=\boldsymbol{1}$ and
$\widehat{\Psi}\Psi\boldsymbol{1}=\boldsymbol{1}$, $I-\Psi\widehat{\Psi}$ and
$I-\widehat{\Psi}\Psi$ are both singular, and the argument above fails.
Instead, one must develop the Jordan chain argument from the proof of [30,
Lemma 4.6].

As a matter of fact, we often find that $\mu=0$ is a case that presents
additional difficulties. This will be seen in Theorem 8.5 about escape
probabilities — that is but one example. Moreover, the convergence of
computational algorithms is much slower (Guo [34]).

##  7 First passage times

The matrix $H$ analysed in Section 2 gives us the joint probability that a
regeneration interval is finite, and the phase at the end of the interval.
Here, we are interested in the distribution of the length of the regeneration
interval, that is, the distribution functions

$H_{ij}(t)={\rm P}[\theta_{n+1}-\theta_{n}\leq
t,\bar{\varphi}_{n+1}=j|\bar{\varphi}_{n}=i],$

for $i$, $j$ in $\mathcal{S}_{\mbox{\tiny$-$}}$. Using the same partition that
gave us (6), we write $H(t)$ as the convolution product $H(t)=\Phi*\Psi(t)$,
with

$\displaystyle\Phi_{ik}(t)$ $\displaystyle={\rm P}[\delta_{n}-\theta_{n}\leq
t,\varphi(\delta_{n})=k|\bar{\varphi}_{n}=i],$ $\displaystyle
i\in\mathcal{S}_{\mbox{\tiny$-$}},k\in\mathcal{S}_{\mbox{\tiny$+$}},$
$\displaystyle\Psi_{kj}(t)$ $\displaystyle={\rm P}[\theta_{n+1}-\delta_{n}\leq
t,\bar{\varphi}_{n+1}=j|\varphi(\delta_{n})=k],$ $\displaystyle
k\in\mathcal{S}_{\mbox{\tiny$+$}},j\in\mathcal{S}_{\mbox{\tiny$-$}}.$

As usual, it is easier to characterise the regenerative intervals through
their LS transforms

$\widehat{H}_{ij}(s)=\int_{0}^{\infty}e^{-st}\,\mathrm{d}H_{ij}(t),$ (32)

and to write $\widehat{H}(s)=\widehat{\Phi}(s)\widehat{\Psi}(s)$ where
$\widehat{\Phi}(s)$ and $\widehat{\Psi}(s)$ are the matrices of LS transforms
of $\Phi(t)$ and $\Psi(t)$, respectively. It is well-known that LS transforms
may be interpreted in probabilistic terms through the introduction of an
exponential random variable $V$ with parameter $s$, independent of the fluid
flow process: we rewrite (32) as

$\widehat{H}_{ij}(s)={\rm P}[\theta_{n+1}-\theta_{n}\leq
V,\bar{\varphi}_{n+1}=j|\bar{\varphi}_{n}=i],$

and similarly

$\displaystyle\widehat{\Phi}_{ik}(s)$ $\displaystyle={\rm
P}[\delta_{n}-\theta_{n}\leq V,\varphi(\delta_{n})=k|\bar{\varphi}_{n}=i],$
$\displaystyle\widehat{\Psi}_{kj}(s)$ $\displaystyle={\rm
P}[\theta_{n+1}-\delta_{n}\leq
V,\bar{\varphi}_{n+1}=j|\varphi(\delta_{n})=k].$

At this point, it is easy to verify that $\widehat{\Phi}(s)$ and
$\widehat{\Psi}(s)$ are given by slight modifications of (8) and (14):

$\widehat{\Phi}(s)=\begin{bmatrix}I&0\end{bmatrix}\left(-\begin{bmatrix}Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}-sI&Q_{{\mbox{\tiny$-$}}0}\\\
Q_{0{\mbox{\tiny$-$}}}&Q_{00}-sI\end{bmatrix}^{-1}\right)\begin{bmatrix}Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}\\\
Q_{0{\mbox{\tiny$+$}}}\end{bmatrix},$ (33)

and $\widehat{\Psi}(s)$ is the minimal nonnegative solution of

$C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}(s)+C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}(s)\widehat{\Psi}(s)+\widehat{\Psi}(s)|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}(s)+\widehat{\Psi}(s)|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}(s)\widehat{\Psi}(s)=0,$
(34)

where

$T(s)=\begin{bmatrix}Q_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}-sI&Q_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}\\\
Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}&Q_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}-sI\end{bmatrix}+\begin{bmatrix}Q_{{\mbox{\tiny$+$}}0}\\\
Q_{{\mbox{\tiny$-$}}0}\end{bmatrix}(sI-
Q_{00})^{-1}\begin{bmatrix}Q_{0{\mbox{\tiny$+$}}}&Q_{0{\mbox{\tiny$-$}}}\end{bmatrix}.$

The Riccati equation may be solved for any given $s$ by the same algorithms as
discussed in Section 9, and this makes it feasible to compute the
distributions themselves by numerical inversion procedures.

Moments of first passage times are obtained by taking the derivatives of
$\widehat{H}(s)=\widehat{\Phi}(s)\widehat{\Psi}(s)$ and evaluating it for
$s=0$. Derivatives of $\widehat{\Phi}(s)$ are easily obtained from (33) but
those of $\widehat{\Psi}(s)$ are more involved, as shown in Bean et al. [14],
and are expressed as solutions of nonsingular Sylvester equations. The first
moment may be obtained in a more straightforward manner, as we show in (27).

##  8 Escape from an interval

Section 3 is about the distribution of the phase upon the first passage of
$X(t)$ to a given level. Here, we deal with the first passage to the boundary
of a finite interval: assuming that $X(0)=0$, we look for the distribution of
the phase when the level escapes for the first time from the interval
$(-a,b)$, with $a$ and $b\geq 0$.

We define as follows the matrices $A^{(a,b)}$ and $B^{(a,b)}$ indexed by
$\mathcal{S}\times\mathcal{S}_{\mbox{\tiny$-$}}$ and
$\mathcal{S}\times\mathcal{S}_{\mbox{\tiny$+$}}$, respectively:

$\displaystyle A_{ij}^{(a,b)}$ $\displaystyle={\rm
P}[\tau_{a}^{\mbox{\tiny$-$}}<\tau_{b}^{\mbox{\tiny$+$}},\varphi(\tau_{a}^{\mbox{\tiny$-$}})=j|X(0)=0,\varphi(0)=i],\qquad\mbox{$i\in\mathcal{S}$,
$j\in\mathcal{S}_{\mbox{\tiny$-$}}$,}$ $\displaystyle B_{ij}^{(a,b)}$
$\displaystyle={\rm
P}[\tau_{b}^{\mbox{\tiny$+$}}<\tau_{a}^{\mbox{\tiny$-$}},\varphi(\tau_{b}^{\mbox{\tiny$+$}})=j|X(0)=0,\varphi(0)=i],\qquad\mbox{$i\in\mathcal{S}$,
$j\in\mathcal{S}_{\mbox{\tiny$+$}}$.}$

The matrices are partitioned into the usual subblocks:

$\begin{bmatrix}B^{(a,b)}&A^{(a,b)}\end{bmatrix}=\begin{bmatrix}B_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}^{(a,b)}\\\
B_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}^{(a,b)}\\\
B_{0{\mbox{\tiny$+$}}}^{(a,b)}&A_{0{\mbox{\tiny$-$}}}^{(a,b)}\end{bmatrix}.$

If $\varphi(0)$ is in $\mathcal{S}_{0}$, the process remains at level 0 for a
while, before jumping to a phase in either $\mathcal{S}_{\mbox{\tiny$+$}}$ or
$\mathcal{S}_{\mbox{\tiny$-$}}$. We condition on the first phase visited
either in $\mathcal{S}_{\mbox{\tiny$+$}}$ or $\mathcal{S}_{\mbox{\tiny$-$}}$
and we find that

$\begin{bmatrix}B_{0{\mbox{\tiny$+$}}}^{(a,b)}&A_{0{\mbox{\tiny$-$}}}^{(a,b)}\end{bmatrix}=(-Q_{00})^{-1}\begin{bmatrix}Q_{0{\mbox{\tiny$+$}}}&Q_{0{\mbox{\tiny$-$}}}\end{bmatrix}\begin{bmatrix}B_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}^{(a,b)}\\\
B_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}^{(a,b)}\end{bmatrix}.$
(35)

###### Remark 8.1

We see here why we defined $\tau_{x}^{\mbox{\tiny$-$}}$ in Section 3 as
$\tau^{\mbox{\tiny$-$}}_{x}=\inf\\{t:X(t)<X(0)-x\\}$. If we had defined it as
$\tau^{\mbox{\tiny$-$}}_{x}=\inf\\{t:X(t)=X(0)-x\\}$, then (35) would not have
held for $a=0$. Similarly, it would not have held for $b=0$ if we had defined
$\tau^{\mbox{\tiny$+$}}_{x}=\inf\\{t:X(t)\geq X(0)+x\\}$.

We need the following lemma to determine the remaining entries of $B^{(a,b)}$
and $A^{(a,b)}$.

###### Lemma 8.2

For $\varphi(0)$ in $\mathcal{S}_{\mbox{\tiny$+$}}$ or
$\mathcal{S}_{\mbox{\tiny$-$}}$, the escape probability matrices are solutions
of the linear system

$\begin{bmatrix}B_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}^{(a,b)}\\\
B_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}^{(a,b)}\end{bmatrix}(I+\mathcal{P})=\mathcal{U},$
(36)

where

$\displaystyle\mathcal{P}$ $\displaystyle=\begin{bmatrix}0&\Psi e^{U(a+b)}\\\
\widehat{\Psi}e^{\widehat{U}(a+b)}&0\end{bmatrix},$ $\displaystyle\mathcal{U}$
$\displaystyle=\begin{bmatrix}e^{\widehat{U}b}&\Psi e^{Ua}\\\
\widehat{\Psi}e^{\widehat{U}b}&e^{Ua}\end{bmatrix}.$

Proof  Assume $\varphi(0)$ is in $\mathcal{S}_{\mbox{\tiny$+$}}$. We have

$e^{\widehat{U}b}=B_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}^{(a,b)}+A_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}^{(a,b)}\widehat{\Psi}e^{\widehat{U}(a+b)}.$
(37)

Indeed, the left-hand side gives the distribution of the phase when the
process has moved up from level 0 to level $b$, it is decomposed in the right-
hand side as the sum of the probability that the process reaches $b$ without
going down to $-a$ and the probability that it goes down to $-a$ first, then
returns to level $-a$ from below, and eventually goes up by $a+b$ units, from
$-a$ to $b$. Similarly,

$\Psi
e^{Ua}=A_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}^{(a,b)}+B_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}^{(a,b)}\Psi
e^{U(a+b)}.$ (38)

Equations (37, 38) form the first row of the system (36). The argument for the
second row is similar. $\square$

If $\mu\not=0$, then (36) is nonsingular and it has a unique solution. The
reason is that either $\Psi e^{U(a+b)}$ or
$\widehat{\Psi}e^{\widehat{U}(a+b)}$ is substochastic, and so the series
$\sum_{\nu\geq 0}(-1)^{\nu}\mathcal{P}^{\nu}$ is converging to
$(I+\mathcal{P})^{-1}$ (see da Silva Soares and Latouche [24] for details). If
$\mu=0$, then both $\Psi e^{U(a+b)}$ and $\widehat{\Psi}e^{\widehat{U}(a+b)}$
are stochastic matrices, $I+\mathcal{P}$ is singular and we need to add one
equation.

Equation (39) is one such choice, as we prove in Theorem 8.5. This equation is
identical to the one given for Markov-modulated Brownian motion in Ivanovs
[36] in a comment after Theorem 3.1, referring to a result obtained in D’Auria
et al. [27, Section 7] by a spectral decomposition argument. The proof given
here is based on the analysis of the stochastic process itself, it is new and
for that reason we give all technical details.

###### Lemma 8.3

If $\mu=0$, then the escape probability matrix is such that

$\begin{bmatrix}B_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}^{(a,b)}\\\
B_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}^{(a,b)}\end{bmatrix}\boldsymbol{\beta}=\begin{bmatrix}\boldsymbol{h}_{\mbox{\tiny$+$}}\\\
\boldsymbol{h}_{\mbox{\tiny$-$}}\end{bmatrix},$ (39)

where $\boldsymbol{h}=-Q^{\\#}\boldsymbol{c}$ and

$\boldsymbol{\beta}=\begin{bmatrix}\ \
b\boldsymbol{1}_{\mbox{\tiny$+$}}+\boldsymbol{h}_{\mbox{\tiny$+$}}\\\
-a\boldsymbol{1}_{\mbox{\tiny$-$}}+\boldsymbol{h}_{\mbox{\tiny$-$}}\end{bmatrix}.$
(40)

Proof  We know from Coolen-Schrijner and van Doorn [21, Section 3] that the
deviation matrix of $Q$ is equal to the group inverse of $-Q$, that is,
$\int_{0}^{\infty}(e^{Qu}-\boldsymbol{1}\cdot\boldsymbol{\alpha}^{\textnormal{\tiny
T}})\,\mathrm{d}u=-Q^{\\#}$. Therefore,

$\displaystyle\boldsymbol{h}$
$\displaystyle=\int_{0}^{\infty}(e^{Qu}-\boldsymbol{1}\cdot\boldsymbol{\alpha}^{\textnormal{\tiny
T}})\,\mathrm{d}u\ \boldsymbol{c}$
$\displaystyle=\lim_{t\rightarrow\infty}\int_{0}^{t}(e^{Qu}-\boldsymbol{1}\cdot\boldsymbol{\alpha}^{\textnormal{\tiny
T}})\,\mathrm{d}u\ \boldsymbol{c}$
$\displaystyle=\lim_{t\rightarrow\infty}\int_{0}^{t}(e^{Qu}\boldsymbol{c}-\boldsymbol{1}\cdot\boldsymbol{\alpha}^{\textnormal{\tiny
T}}\boldsymbol{c})\,\mathrm{d}u$
$\displaystyle=\lim_{t\rightarrow\infty}\int_{0}^{t}e^{Qu}\boldsymbol{c}\,\mathrm{d}u$
(41)

as $\boldsymbol{\alpha}^{\textnormal{\tiny T}}\boldsymbol{c}=\mu$ and $\mu=0$
by assumption. Furthermore, $(e^{Qu})_{ij}$ is the probability ${\rm
P}[\varphi(u)=j|\varphi(0)=i]$ and so (41) may be interpreted by (1) as

$h_{i}=\lim_{t\rightarrow\infty}{\rm E}[X_{t}|X(0)=0,\varphi(0)=i],$

and we write $\boldsymbol{h}={\rm E}_{0}[X_{\infty}|\varphi(0)]$ for short. By
conditioning on the first time the process escapes from $(-a,b)$, it is
straightforward to verify that

$\displaystyle\boldsymbol{h}$ $\displaystyle={\rm
E}_{0}[X_{\infty}\mathds{1}\\{\tau_{b}^{\mbox{\tiny$+$}}<\tau_{a}^{\mbox{\tiny$-$}}\\}|\varphi(0)]+{\rm
E}_{0}[X_{\infty}\mathds{1}\\{\tau_{a}^{\mbox{\tiny$-$}}<\tau_{b}^{\mbox{\tiny$+$}}\\}|\varphi(0)]$
$\displaystyle=B^{(a,b)}(b\boldsymbol{1}_{\mbox{\tiny$+$}}+\boldsymbol{h}_{\mbox{\tiny$+$}})+A^{(a,b)}(-a\boldsymbol{1}_{\mbox{\tiny$-$}}+\boldsymbol{h}_{\mbox{\tiny$-$}})$

and this concludes the proof. $\square$

###### Remark 8.4

The vector $\boldsymbol{h}$ has a number of interesting properties that we
need later. We observe that $\boldsymbol{\alpha}^{\textnormal{\tiny
T}}\boldsymbol{h}=0$ since $\boldsymbol{\alpha}^{\textnormal{\tiny
T}}Q^{\\#}=\boldsymbol{0}$. Furthermore,

$\displaystyle\boldsymbol{h}_{\mbox{\tiny$+$}}$
$\displaystyle=\Psi\boldsymbol{h}_{\mbox{\tiny$-$}}\quad\mbox{if\ }\mu\leq 0,$
(42) $\displaystyle\boldsymbol{h}_{\mbox{\tiny$-$}}$
$\displaystyle=\Psi\boldsymbol{h}_{\mbox{\tiny$+$}}\quad\mbox{if\ }\mu\geq 0.$
(43)

To justify the first equation, we use the interpretation given to
$\boldsymbol{h}$, and we use the fact that if $\mu\leq 0$, starting from a
phase in $\mathcal{S}_{\mbox{\tiny$+$}}$, the process returns to 0 in finite
time with probability 1. The justification of (43) is similar.

In summary, the distribution of the phase upon escaping from the interval
$(-a,b)$ is given in the next theorem.

###### Theorem 8.5

The distribution of the phase at first escape from $(-a,b)$ is given by

$\begin{bmatrix}B_{0{\mbox{\tiny$+$}}}^{(a,b)}&A_{0{\mbox{\tiny$-$}}}^{(a,b)}\end{bmatrix}=(-Q_{00})^{-1}\begin{bmatrix}Q_{0{\mbox{\tiny$+$}}}&Q_{0{\mbox{\tiny$-$}}}\end{bmatrix}\begin{bmatrix}B_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}^{(a,b)}\\\
B_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}^{(a,b)}\end{bmatrix}$

if $\varphi(0)\in\mathcal{S}_{0}$. If $\mu\not=0$ and $\varphi(0)$ is in
$\mathcal{S}_{\mbox{\tiny$+$}}$ or $\mathcal{S}_{\mbox{\tiny$-$}}$, then

$\begin{bmatrix}B_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}^{(a,b)}\\\
B_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}^{(a,b)}\end{bmatrix}=\mathcal{U}(I+\mathcal{P})^{-1},$
(44)

where $\mathcal{U}$ and $\mathcal{P}$ are defined in Lemma 8.2.

If $\mu=0$ and $\varphi(0)$ is in $\mathcal{S}_{\mbox{\tiny$+$}}$ or
$\mathcal{S}_{\mbox{\tiny$-$}}$, then

$\begin{bmatrix}B_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}^{(a,b)}\\\
B_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}^{(a,b)}\end{bmatrix}=\mathcal{U}(I+\mathcal{P})^{\\#}+\boldsymbol{w}\cdot\boldsymbol{\eta}^{\textnormal{\tiny
T}},$ (45)

where $\boldsymbol{\eta}$ is the left eigenvector of $I+\mathcal{P}$ for the
eigenvalue 0,

$\boldsymbol{w}=(\boldsymbol{\eta}^{\textnormal{\tiny
T}}\boldsymbol{\beta})^{-1}(\boldsymbol{h}-\mathcal{U}(I+\mathcal{P})^{\\#}\boldsymbol{\beta}),$

and $\boldsymbol{\beta}$ is defined in (40).

Proof  The first two statements have been justified before, we include them
for completeness. If $\mu=0$, the matrix $\mathcal{P}$ is stochastic and
irreducible. As

$I-\mathcal{P}=\begin{bmatrix}I&-\Psi e^{U(a+b)}\\\
-\widehat{\Psi}e^{\widehat{U}(a+b)}&I\end{bmatrix}=\begin{bmatrix}I&0\\\
0&-I\end{bmatrix}(I+\mathcal{P})\begin{bmatrix}I&0\\\ 0&-I\end{bmatrix},$

the matrices $I-\mathcal{P}$ and $I+\mathcal{P}$ are similar and
$I+\mathcal{P}$ has a unique eigenvalue equal to 0. The corresponding left
eigenvector $\boldsymbol{\eta}$ is such that

$\displaystyle\boldsymbol{\eta}^{\textnormal{\tiny T}}_{\mbox{\tiny$-$}}$
$\displaystyle=-\boldsymbol{\eta}^{\textnormal{\tiny T}}_{\mbox{\tiny$+$}}\Psi
e^{U(a+b)},$ (46) $\displaystyle\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}$ $\displaystyle=\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}\Psi e^{U(a+b)}\widehat{\Psi}e^{\widehat{U}(a+b)},$

and we may choose $\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}>\boldsymbol{0}$. Thus, the system (36) has the solution

$\begin{bmatrix}B_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}^{(a,b)}\\\
B_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}^{(a,b)}&A_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}^{(a,b)}\end{bmatrix}=\mathcal{U}(I+\mathcal{P})^{\\#}+\boldsymbol{w}\cdot\boldsymbol{\eta}^{\textnormal{\tiny
T}}$ (47)

for some vector $\boldsymbol{w}$. We post-multiply (47) by
$\boldsymbol{\beta}$, the left-hand side is equal to $\boldsymbol{h}$ by Lemma
8.3 and we obtain

$\boldsymbol{w}=(\boldsymbol{\eta}^{\textnormal{\tiny
T}}\boldsymbol{\beta})^{-1}(\boldsymbol{h}-\mathcal{U}(I+\mathcal{P})^{\\#}\boldsymbol{\beta}),$

provided that $\boldsymbol{\eta}^{\textnormal{\tiny
T}}\boldsymbol{\beta}\not=0$. This is equivalent to showing that
$\boldsymbol{\beta}$ is, indeed, linearly independent of the columns of
$I+\mathcal{P}$, which in turn implies that (36, 39) is a non-singular system
when $\mu=0$.

Now,

$\displaystyle\boldsymbol{\eta}^{\textnormal{\tiny T}}\boldsymbol{\beta}$
$\displaystyle=b\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}\boldsymbol{1}_{\mbox{\tiny$+$}}+\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}\boldsymbol{h}_{\mbox{\tiny$+$}}-a\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$-$}}\boldsymbol{1}_{\mbox{\tiny$-$}}+\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$-$}}\boldsymbol{h}_{\mbox{\tiny$-$}}$
$\displaystyle=b\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}\boldsymbol{1}_{\mbox{\tiny$+$}}+a\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}\boldsymbol{1}_{\mbox{\tiny$+$}}+\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}(\boldsymbol{h}_{\mbox{\tiny$+$}}-\Psi
e^{U(a+b)}\boldsymbol{h}_{\mbox{\tiny$-$}})$ (48)

by (46). The vector $\boldsymbol{h}_{\mbox{\tiny$-$}}$ is indexed by phases in
$\mathcal{S}_{\mbox{\tiny$-$}}$ and we write, for short,

$\displaystyle\boldsymbol{h}_{\mbox{\tiny$-$}}$
$\displaystyle=\lim_{t\rightarrow\infty}{\rm
E}_{0}[X(t)|\varphi(0)\in\mathcal{S}_{\mbox{\tiny$-$}}]$
$\displaystyle=\lim_{t\rightarrow\infty}({\rm
E}_{0}[X(t)\mathbb{1}\\{\tau_{0}^{\mbox{\tiny$+$}}<t\\}|\varphi(0)\in\mathcal{S}_{\mbox{\tiny$-$}}]+{\rm
E}_{0}[X(t)\mathbb{1}\\{\tau_{0}^{\mbox{\tiny$+$}}\geq
t\\}|\varphi(0)\in\mathcal{S}_{\mbox{\tiny$-$}}])$
$\displaystyle=\lim_{t\rightarrow\infty}({\rm
P}[\tau_{0}^{\mbox{\tiny$+$}}<t,\varphi(\tau_{0}^{\mbox{\tiny$+$}})|\varphi(0)\in\mathcal{S}_{\mbox{\tiny$-$}}]\,{\rm
E}_{0}[X(t)|\tau_{0}^{\mbox{\tiny$+$}}<t,X(\tau_{0}^{\mbox{\tiny$+$}}),\varphi(\tau_{0}^{\mbox{\tiny$+$}})]$
$\displaystyle\quad+{\rm P}[\tau_{0}^{\mbox{\tiny$+$}}\geq
t|\varphi(0)\in\mathcal{S}_{\mbox{\tiny$-$}}]\,{\rm
E}_{0}[X(t)|\tau_{0}^{\mbox{\tiny$+$}}\geq
t,\varphi(0)\in\mathcal{S}_{\mbox{\tiny$-$}}])$
$\displaystyle=\widehat{\Psi}\boldsymbol{h}_{\mbox{\tiny$+$}}+\lim_{t\rightarrow\infty}{\rm
P}[\tau_{0}^{\mbox{\tiny$+$}}\geq
t|\varphi(0)\in\mathcal{S}_{\mbox{\tiny$-$}}]\,{\rm
E}_{0}[X(t)|\tau_{0}^{\mbox{\tiny$+$}}\geq
t,\varphi(0)\in\mathcal{S}_{\mbox{\tiny$-$}}])$ as
$X(\tau_{0}^{\mbox{\tiny$+$}})=0$ and
$\varphi(\tau_{0}^{\mbox{\tiny$+$}})\in\mathcal{S}_{\mbox{\tiny$+$}}$,
$\displaystyle\leq\widehat{\Psi}\boldsymbol{h}_{\mbox{\tiny$+$}}$

since $X(t)<0$ for $t<\tau_{0}^{\mbox{\tiny$+$}}$. The inequality is strict
for at least one component of $\boldsymbol{h}_{\mbox{\tiny$-$}}$. Otherwise,
by Remark 8.4, we would get
$\boldsymbol{h}_{\mbox{\tiny$-$}}=\widehat{\Psi}\boldsymbol{h}_{\mbox{\tiny$+$}}=\widehat{\Psi}\Psi\boldsymbol{h}_{\mbox{\tiny$-$}}$,
from which we would successively conclude that
$\boldsymbol{h}_{\mbox{\tiny$-$}}=c\boldsymbol{1}$ for some scalar $c$, that
$\boldsymbol{h}_{\mbox{\tiny$+$}}=c\boldsymbol{1}$, and that
$\boldsymbol{h}=c\boldsymbol{1}$, which would be in contradiction with
$\boldsymbol{\alpha}^{\textnormal{\tiny T}}\boldsymbol{h}=0$, by Remark 8.4.
Consequently, (48) becomes

$\boldsymbol{\eta}^{\textnormal{\tiny
T}}\boldsymbol{\beta}>(a+b)\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}\boldsymbol{1}_{\mbox{\tiny$+$}}+\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}(I-\Psi
e^{U(a+b)}\widehat{\Psi})\boldsymbol{h}_{\mbox{\tiny$+$}}.$ (49)

By an argument similar to the one held above, we find that

$\displaystyle\boldsymbol{h}_{\mbox{\tiny$+$}}$
$\displaystyle=\lim_{t\rightarrow\infty}{\rm
E}_{0}[X(t)-(a+b)|\varphi(0)\in\mathcal{S}_{\mbox{\tiny$+$}}]+(a+b)\boldsymbol{1}$
$\displaystyle=\lim_{t\rightarrow\infty}({\rm
P}_{0}[\tau_{a+b}^{\mbox{\tiny$+$}}<t,\varphi(\tau_{a+b}^{\mbox{\tiny$+$}})|\varphi(0)\in\mathcal{S}_{\mbox{\tiny$+$}}]{\rm
E}_{0}[X(t)-(a+b)|\tau_{a+b}^{\mbox{\tiny$+$}}<t,\varphi(\tau_{a+b}^{\mbox{\tiny$+$}})]$
$\displaystyle\quad+(a+b)\boldsymbol{1}$ $\displaystyle\quad+{\rm
P}_{0}[\tau_{a+b}^{\mbox{\tiny$+$}}\geq
t|\varphi(0)\in\mathcal{S}_{\mbox{\tiny$+$}}]{\rm
E}_{0}[(X(t)-(a+b)|\tau_{a+b}^{\mbox{\tiny$+$}}\geq
t,\varphi(0)\in\mathcal{S}_{\mbox{\tiny$-$}}])$ $\displaystyle\leq
e^{\widehat{U}(a+b)}\boldsymbol{h}_{\mbox{\tiny$+$}}+(a+b)\boldsymbol{1}$

as $X(\tau_{a+b}^{\mbox{\tiny$+$}})-(a+b)=0$,
$\varphi(\tau_{a+b}^{\mbox{\tiny$+$}})\in\mathcal{S}_{\mbox{\tiny$+$}}$, and
$X(t)-(a+b)<0$ for $t<\tau_{a+b}^{\mbox{\tiny$+$}}$. With this, (49) becomes

$\displaystyle\boldsymbol{\eta}^{\textnormal{\tiny T}}\boldsymbol{\beta}$
$\displaystyle>(a+b)\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}\boldsymbol{1}_{\mbox{\tiny$+$}}+\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}(I-\Psi
e^{U(a+b)}\widehat{\Psi}e^{\widehat{U}(a+b)})\boldsymbol{h}_{\mbox{\tiny$+$}}-(a+b)\boldsymbol{\eta}^{\textnormal{\tiny
T}}_{\mbox{\tiny$+$}}\boldsymbol{1}_{\mbox{\tiny$+$}}$ $\displaystyle=0.$

This completes the proof. $\square$

In the same manner as in Section 7, we may determine the LS transform of the
random variable $\min(\tau_{a}^{\mbox{\tiny$-$}},\tau_{b}^{\mbox{\tiny$+$}})$,
details are in Bean et al. [13].

##  9 Numerical procedures

It should be clear by this point that the numerical evaluation of many
quantities of interest is dependent on being able to compute the matrices
$\Psi$ and $\widehat{\Psi}$. If we replace $U$ in (15) by the right-hand side
of (13) and write

$\Psi=\int_{0}^{\infty}e^{C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}y}\,C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}\,e^{(|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}+|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}\Psi)y}\,\mathrm{d}y,$

then an obvious approach to compute $\Psi$ is to proceed by successive
substitution: define iteratively

$\Psi_{n+1}=\int_{0}^{\infty}e^{C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}y}\,C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}\,e^{(|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}+|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}\Psi_{n})y}\,\mathrm{d}y.$
(50)

for $n\geq 0$, starting from $\Psi_{0}=0$, The resulting sequence is
monotonically convergent to $\Psi$ as we show in Theorem 9.1. The proof is
new, and we give it in detail. It is based on the evolution of a stack
$\sigma$ associated to the fluid queue.

At the epochs when the phase process enters $\mathcal{S}_{\mbox{\tiny$+$}}$
after a sojourn in $\mathcal{S}_{\mbox{\tiny$-$}}$, we put the level at the
top of the stack; the value recorded on top of the stack is removed when the
fluid decreases to that level. Formally, we define the sequence
$\\{s_{k}:0\leq k\leq L\\}$ of epochs where the stack increases during the
interval $[0,\Theta]$:

$\displaystyle s_{0}$ $\displaystyle=0,$ $\displaystyle s_{k}$
$\displaystyle=\inf\\{t>f_{k}:\varphi(t)\in\mathcal{S}_{\mbox{\tiny$+$}}\\}\qquad\mbox{with
$f_{k}=\inf\\{t>s_{k-1}:\varphi(t)\in\mathcal{S}_{\mbox{\tiny$-$}}\\}$,}$

and $L=\sup\\{k:s_{k}<\Theta\\}$. On the sample path of Figure 5, $L=6$ and we
have marked $s_{0}$ to $s_{6}$.

$t$$s^{\prime}_{0}$$s_{0}$$s_{1}$$s_{2}$$s_{3}$$s_{4}$$s_{5}$$s_{6}$$X(t)$$y$$0$

Figure 5: The fluid level $X(s_{i})$ is pushed on the stack at time $s_{i}$,
$i=0,1\ldots$. The epochs $s^{\prime}_{1}$ to $s^{\prime}_{6}$ of removal of
the top are marked but not labeled, they are ordered as
$s^{\prime}_{3}<s^{\prime}_{2}<s^{\prime}_{1}<s^{\prime}_{4}<s^{\prime}_{6}<s^{\prime}_{5}<s^{\prime}_{0}$.
The maximum size of the stack is 4, reached at time $s_{3}$.

Next, we define $\\{s^{\prime}_{k}:0\leq k\leq L\\}$:

$s^{\prime}_{k}=\inf\\{t>s_{k}:X(t)=X(s_{k})\\}.$

On Figure 5, the epochs $s^{\prime}_{1}$ to $s^{\prime}_{6}$ are marked, but
not labeled so as not to clutter the graph.

At time $s_{k}$, the size $|\sigma|$ increases by one and we record $X(s_{k})$
on top of $\sigma$, at time $s^{\prime}_{k}$ we remove the top of $\sigma$ and
$|\sigma|$ decreases by one. Note that $s^{\prime}_{0}=\Theta$, and that the
stack becomes empty for the first time.

###### Theorem 9.1

The sequences $\Psi_{n}$, $n\geq 0$ defined by

$(\Psi_{n})_{ij}={\rm P}[\Theta<\infty,\varphi(\Theta)=j,\max_{0\leq
t\leq\Theta}|\sigma(t)|\leq n|\varphi(0)=i],$ (51)

for $i$ in $\mathcal{S}_{\mbox{\tiny$+$}}$, $j$ in
$\mathcal{S}_{\mbox{\tiny$-$}}$, satisfies (50). Furthermore, $\Psi_{n}$ is
the unique solution of the linear equation

$C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}+C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}\Psi_{n}+\Psi_{n}U_{n-1}=0,$
(52)

where

$U_{n-1}=|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}+|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}\Psi_{n-1}$
(53)

for $n\geq 1$.

The sequence $\\{\Psi_{n}\\}$ converges monotonically to $\Psi$ and
$\\{U_{n}\\}$ converges to $U$, as $n\rightarrow\infty$.

Proof  It is obvious that the sequence defined by (51) is monotone and
converges to $\Psi$ as there are fewer constraints on the trajectories for
increasing $n$ until there is none in the limit.

Next, we show that the return probabilities defined in (51) are solutions of
(50). For $n=1$, we have

$\Psi_{1}=\int_{0}^{\infty}e^{C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}y}\,C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}\,e^{|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}y}\,\mathrm{d}y.$

This means that the fluid queue spends some time in
$\mathcal{S}_{\mbox{\tiny$+$}}$ and grows up to some level $y$, then switches
to $\mathcal{S}_{\mbox{\tiny$-$}}$ and never returns to
$\mathcal{S}_{\mbox{\tiny$+$}}$ until hitting level 0. We push the value 0 on
$\sigma$ at time 0, remove it at time $\Theta$ and so $|\sigma|=1$ over the
whole interval. This shows that $\Psi_{1}$ is a solution of (50) with
$\Psi_{0}=0$.

For the general case, we illustrate on Figure 6 the physical meaning of the
right-hand side of (50): the fluid process grows up to some level $y$, then
goes down to 0 with occasional episodes of growth; during those episodes, the
trajectories followed by the process are constrained by the definition of
$\Psi_{n}$.

$t$$s^{\prime}_{0}$$s_{0}$$a_{1}$$b_{1}$$\Psi_{n}$$a_{2}$$b_{2}$$\Psi_{n}$$a_{3}$$b_{3}$$\Psi_{n}$$X(t)$$y$$0$

Figure 6: Approximation $\Psi_{n+1}$ by functional iteration: the size of the
stack increases at most by $n$ over the intervals $[a_{1},b_{1}]$,
$[a_{2},b_{2}]$ and $[a_{3},b_{3}]$ and is at most $n+1$ over the whole
interval $[s_{0},s^{\prime}_{0}]$.

Such episodes, if any, occur during intervals $(a_{i},b_{i})$, $1\leq i\leq
L^{\prime}$, with

$\displaystyle a_{1}$ $\displaystyle=s_{1},$ $\displaystyle b_{i}$
$\displaystyle=\inf\\{t>a_{i}:X(t)=X(a_{i})\\}$ $\displaystyle a_{i+1}$
$\displaystyle=\inf\\{t>b_{i}:\varphi(t)\in\mathcal{S}_{\mbox{\tiny$+$}}\\},$

$L^{\prime}=\sup\\{i:a_{i}<\Theta\\}$. During the intervals
$(s_{0},a_{1})\cup(\cup_{1\leq i\leq
L^{\prime}-1}(b_{i},a_{i+1})\cup(b_{L^{\prime}},s^{\prime}_{0})$, $\sigma$
contains only the value 0 which was pushed at time 0, and $|\sigma|=1$. During
the intervals $(a_{i},b_{i})$, the stack may increase at most by $n$ units.
Thus, $|\sigma|\leq n+1$ over the whole interval $[0,\Theta]$ and the right-
hand side of (50) is equal to $\Psi_{n+1}$.

With $U_{n}$ defined in (53), (50) becomes

$\Psi_{n+1}=\int_{0}^{\infty}e^{C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}y}\,C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}\,e^{U_{n}y}\,\mathrm{d}y.$
(54)

and (52) follows in the same manner as we proved (16). The coefficients
$C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}$ and $U_{n}$
are both defective generators, all of their eigenvalues are in
$\mathbb{C}_{<0}$, and so the system (52) has a unique solution $\Psi_{n}$.
Therefore, the sequence $\\{U_{n},\Psi_{n}\\}$ is well-defined and the theorem
follows. $\square$

The algorithm defined by (52, 53) is easily implemented and is the most
efficient among several linearly convergent algorithms, as shown in Bean et
al. [8]. Several other procedures have been proposed in Guo [32] and Bini et
al. [16, 17]. A special mention should be made of the Newton method as it is
easily implemented and much faster than functional iteration: the sequence
$\\{\Psi^{(N)}_{n}\\}$ defined by

$\displaystyle C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$-$}}}$
$\displaystyle+(C_{\mbox{\tiny$+$}}^{-1}T_{{\mbox{\tiny$+$}}{\mbox{\tiny$+$}}}+\Psi^{(N)}_{n}|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}})\Psi^{(N)}_{n+1}$
$\displaystyle+\Psi^{(N)}_{n+1}(|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$-$}}}+|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}\Psi^{(N)}_{n})=\Psi^{(N)}_{n}|C_{\mbox{\tiny$-$}}|^{-1}T_{{\mbox{\tiny$-$}}{\mbox{\tiny$+$}}}\Psi^{(N)}_{n}$

for $n\geq 0$, with $\Psi^{(N)}_{0}=0$, is well-defined and converges
quadratically if $\mu\not=0$.

The most efficient algorithms today form the family of doubling algorithms,
which solve simultaneously for $\Psi$ and for $\widehat{\Psi}$; they include
the structure-preserving doubling algorithm (SDA, Guo et al. [35]), SDA
shrink-and-shift (Bini et al. [15]), componentwise-accurate doubling
algorithms (Nguyen and Poloni [43]), and alternating-directional doubling
algorithm (Wang et al. [48]). These algorithms are quadratically convergent if
$\mu\neq 0$ and, furthermore, each iteration is faster due to fewer
computations.

##  10 Extensions

After the publication of Ramaswami’s seminal paper [44], the basic fluid queue
model defined through Equations (1) and (3) has been extended in many ways. We
cite some of these, without getting into details, the list is far from
exhaustive.

#### Finite buffers

In some applications, $Y(t)$ represents the content of some finite buffer, so
that the level may only take values in some interval $0\leq Y(t)\leq B<\infty$
(see da Silva Soares and Latouche [23, 24] and Bean et al. [13]). In that
case, it is natural to choose as regeneration points the epochs of return to
the upper boundary in addition to the returns to level 0, and the transition
matrix $H$ between regeneration points makes use of first passage
probabilities from boundary to boundary. These may be obtained from Lemma 8.2
with $a=0$ and $b=B$, if the process starts from the lower boundary, or $a=B$
and $b=0$, if the process starts from the upper boundary.

The final expression for the stationary density is a mixture of two matrix
exponentials, $e^{Kx}$ and $e^{\widehat{K}(B-x)}$. Details may be found in
[25] and in [38].

#### Level feedback

In our presentation so far, the level is driven by the phase, subject to the
boundary constrains that $Y(t)\geq 0$, or $0\leq Y(t)\leq B$, while the phase
evolves independently of the level. In many cases, the level has a direct
influence on the evolution of the phase; for instance, one might reduce the
flow into the buffer as the level gets nears the upper boundary, so as to
avoid spillage.

In Bean et al. [9] and da Silva Soares [25] (and other references cited
there), one defines a number of threshold values $0\leq
b_{1}<b_{2}<\cdots<b_{N}\leq\infty$ such that the parameters $C$ and $Q$ of
the phase process change upon the crossing of a threshold. Here, regenerations
occur when the fluid reaches any of the threshold, and the analysis of such
systems may be based on a systematic extension of the results for the system
with a finite buffer.

#### Fluid with jumps

In Remiche [45] $Y(t)$ represents the supply of tokens in a leaky bucket
system. It increases linearly in time and drops by a positive amount each time
a file is transmitted. In Bean et al. [12], $Y(t)$ represents the amount of
wear of a power generator and it may jump from $B$ (indicating that the
generator is unusable) to 0 (indicating that it has been replaced by a new
equipment). As discussed in Badescu et al. [7] and Stanford et al. [47], risk
processes may be analysed as fluid queues with jumps.

If the jumps have a phase-type distribution, then the analysis of the process
requires very little adaptation from the material presented in the present
paper.

#### Markov modulated Brownian motion

The definition of these processes is very similar to that of fluid flows. The
difference is that the fluid evolves like a Brownian motion with parameters
(drift and variance) which depend on $\varphi(t)$. Recent references are
d’Auria et al. [26, 28], Ivanovs [36], Gribaudo et al. [31] where the authors
focused on obtaining time-dependent distributions and first hitting times
using different approaches: stochastic ODE resolution, spectral decomposition
and martingale theory. Breuer [18] determined the occupation time of the
process in an interval before a one- or two-sided exit. Latouche and Nguyen
[38, 39] are two recent papers that follow a regenerative approach similar to
the one developed here.

#### Two-dimensional fluid

A few authors have considered systems where the component $X$ is two-
dimensional: Bean and O’Reilly [10, 11], Foss and Miyazawa [29] and Latouche
et al. [40] among others. The area of two- or higher-dimensional fluids is
still wide open, with many exciting unanswered questions.

### Acknowledgements

The authors would like to acknowledge the support of ACEMS (Australian
Research Council Centre of Excellence for Mathematical and Statistical
Frontiers).

## References

  * [1] S. Ahn and V. Ramaswami.  Fluid flow models and queues: A connection by stochastic coupling.  Stochastic Models, 19:325–348, 2003.
  * [2] D. Anick, D. Mitra, and M. M. Sondhi.  Stochastic theory of a data-handling system with multiple sources.  The Bell System Technical Journal, 61:1871–1894, 1982.
  * [3] S. Asmussen.  Stationary distributions for fluid flow models with or without Brownian noise.  Comm. Statist. Stochastic Models, 11:21–49, 1995.
  * [4] S. Asmussen.  Applied Probability and Queues.  Springer, New York, second edition, 2003.
  * [5] F. Avram and M. Usábel.  Ruin probabilities and deficit for the renewal risk model with phase-type interarrival times.  Astin Bulletin, 34:315–332, 2004.
  * [6] A. Badescu, L. Breuer, A. da Silva Soares, G. Latouche, M.-A. Remiche, and D. Stanford.  Risk processes analyzed as fluid queues.  Scandinavian Actuarial Journal, 105(2):127–141, 2005.
  * [7] A. Badescu, L. Breuer, S. Drekic, G. Latouche, and D. Stanford.  The surplus prior to ruin and the deficit at ruin for a correlated risk process.  Scandinavian Actuarial Journal, pages 433–445, 2005.
  * [8] N. Bean, M. O’Reilly, and P. Taylor.  Algorithms for return probabilities for stochastic fluid flows.  Stochastic Models, 21:149–184, 2005.
  * [9] N. G. Bean and M. M. O’Reilly.  Performance measures of a multi-layer Markovian fluid model.  Ann. Oper. Res., 160:99–120, 2008.
  * [10] N. G. Bean and M. M. O’Reilly.  A stochastic two-dimensional fluid model.  Stochastic Models, 29:31–63, 2013.
  * [11] N. G. Bean and M. M. O’Reilly.  The stochastic fluid-fluid model: A stochastic fluid model driven by an uncountable-state process, which is a stochastic fluid model itself.  Stoch. Proc. Appl., 124:1741–1772, 2014.
  * [12] N. G. Bean, M. M. O’Reilly, and J. E. Sargison.  A stochastic flow model of the operation and maintenance of power generation systems.  IEEE Transactions on Power Systems, 25:1361–1374, 2010.
  * [13] N. G. Bean, M. M. O’Reilly, and P. Taylor.  Hitting probabilities and hitting times for stochastic fluid flows. the bounded model.  Probability in the Engineering and Information Sciences, 23:121–147, 2009.  doi:10.1017/S0269964809000102.
  * [14] N. G. Bean, M. M. O’Reilly, and P. G. Taylor.  Hitting probabilities and hitting times for stochastic fluid flows.  Stoch. Proc. Appl., 115:1530–1556, 2005.
  * [15] D. Bini, B. Meini, and F. Poloni.  Transforming algebraic Riccati equations into unilateral quadratic matrix equations.  Numer. Math., 2010.
  * [16] D. A. Bini, B. Iannazzo, G. Latouche, and B. Meini.  On the solution of Riccati equations arising in fluid queues.  Linear Algebra Appl., 413:474–494, 2006.
  * [17] D. A. Bini, B. Iannazzo, and B. Meini.  Numerical Solution of Algebraic Riccati Equations.  Number 9 in Fundamentals of Algorithms. SIAM, Philadelphia PA, 2012.
  * [18] L. Breuer.  Occupation times for Markov-modulated Brownian motion.  J. Appl. Probab., 49:549–565, 2012.
  * [19] S. L. Campbell and C. D. Meyer.  Generalized Inverses of Linear Transformations.  Dover Publications, New York, 1991.  Republication.
  * [20] E. Çinlar.  Introduction to Stochastic Processes.  Prentice–Hall, Englewood Cliffs, NJ, 1975.
  * [21] P. Coolen-Schrijner and E. A. van Doorn.  The deviation matrix of a continuous-time Markov chain.  Probab. Engrg. Informational Sci., 16:351–366, 2002.
  * [22] A. da Silva Soares and G. Latouche.  Further results on the similarity between fluid queues and QBDs.  In G. Latouche and P. Taylor, editors, Proceedings of the 4th International Conference on Matrix-Analytic Methods, pages 89–106. World Scientific, River Edge, NJ, 2002.
  * [23] A. da Silva Soares and G. Latouche.  A matrix-analytic approach to fluid queues with feedback control.  I. J. of Simulation, 6:4–12, 2005.
  * [24] A. da Silva Soares and G. Latouche.  Matrix-analytic methods for fluid queues with finite buffers.  Performance Evaluation, 63:295–314, 2006.
  * [25] A. da Silva Soares and G. Latouche.  Fluid queues with level dependent evolution.  European J. Oper. Res., 196:1041–1048, 2009.
  * [26] B. D’Auria, J. Ivanovs, O. Kella, and M. Mandjes.  First passage of a Markov additive process and generalized Jordan chains.  J. Appl. Probab., 47:1048–1057, 2010.  doi:10.1239/jap/1294170518.
  * [27] B. D’Auria, J. Ivanovs, O. Kella, and M. Mandjes.  First passage process of a Markov additive process with applications to reflection problems.  ArXiv e-prints, arXiv:1006.2965v1 [math.PR], 2010.
  * [28] B. D’Auria, J. Ivanovs, O. Kella, and M. Mandjes.  Two-sided reflection of Markov-modulated Brownian motion.  Stochastic Models, 28:316–332, 2012.  doi:10.1080/15326349.2012.672285.
  * [29] S. Foss and M. Miyazawa.  Two-node fluid network with a heavy-tailed random input: the strong stability case.  J. Appl. Probab., 51A:249–265, 2014.
  * [30] M. Govorun, G. Latouche, and M.-A. Remiche.  Stability for fluid queues: characteristic inequalities.  Stochastic Models, 29:64–88, 2013.
doi: 10.1080/15326349.2013.750533.

  * [31] M. Gribaudo, D. Manini, B. Sericola, and M. Telek.  Second order fluid models with general boundary behaviour.  Annals of Operations Research, 160:69–82, 2008.
  * [32] C.-H. Guo.  Nonsymmetric algebraic Riccati equations and Wiener-Hopf factorization for M-matrices.  SIAM J. Matrix Anal. Appl., 23:225–242, 2001.
  * [33] C.-H. Guo.  A note on the minimal nonnegative solution of a nonsymmetric algebraic Riccati equation.  Linear Algebra Appl., 357:299–302, 2002.
  * [34] C.-H. Guo.  Efficient methods for solving a nonsymmetric algebraic Riccati equation arising in stochastic fluid models.  J. of Comput. Appl. Math., 192:353–373, 2006.
  * [35] X.-X. Guo, W.-W. Lin, and S.-F. Xu.  A structure-preserving doubling algorithm for nonsymmetric algebraic Riccati equation.  Numer. Math., 103:393–412, 2006.
  * [36] J. Ivanovs.  Markov-modulated Brownian motion with two reflecting barriers.  J. Appl. Probab., 47:1034–1047, 2010.
  * [37] P. Lancaster and M. Tismenetsky.  The Theory of Matrices: With Applications.  Academic Press, New York, 1985.
  * [38] G. Latouche and G. T. Nguyen.  Feedback control: Markov-modulated Brownian motion with instantaneous change of phase.  Performance Evaluation, 106:30–49, 2016.
doi:http://dx.doi.org/10.1016/j.peva.2016.09.004; arXiv: 1603.01945.

  * [39] G. Latouche and G. T. Nguyen.  Slowing time: Markov-modulated Brownian motion with a sticky boundary.  Stochastic Models, 33:297–321, 2017.
doi:10.1080/15326349.2017.1284000, arXiv: 1508.00922.

  * [40] G. Latouche, G. T. Nguyen, and Z. Palmowski.  Two-dimensional fluid queues with temporary assistance.  In G. Latouche, V. Ramaswami, J. Sethuraman, K. Sigman, M. S. Squillante, and D. D. Yao, editors, Matrix-Analytic Methods in Stochastic Models, volume 27 of Springer Proceedings in Mathematics & Statistics, chapter 9, pages 187–207. Springer, New York, NY, 2013.
doi : 10.1007/978-1-4614-4909-6_9; arXiv: 1404.3996.

  * [41] G. Latouche and V. Ramaswami.  Introduction to Matrix Analytic Methods in Stochastic Modeling.  ASA-SIAM Series on Statistics and Applied Probability. SIAM, Philadelphia PA, 1999.  Second printing 2011.
  * [42] R. M. Loynes.  A continuous-time treatment of certain queues and infinite dams.  J. Austral. Math. Soc., 2:484–498, 1962.
  * [43] G. T. Nguyen and F. Poloni.  Componentwise accurate fluid queue computations using doubling algorithms.  Numerische Mathematik, 130:763–792, 2015.  arXiv:1406.7301.
  * [44] V. Ramaswami.  Matrix analytic methods for stochastic fluid flows.  In D. Smith and P. Hey, editors, Teletraffic Engineering in a Competitive World (Proceedings of the 16th International Teletraffic Congress), pages 1019–1030. Elsevier Science B.V., Edinburgh, UK, 1999.
  * [45] M.-A. Remiche.  Compliance of the token bucket model for Markovian type traffic.  Stochastic Models, 21:615–630, 2005.
  * [46] L. C. G. Rogers.  Fluid models in queueing theory and Wiener-Hopf factorization of Markov chains.  Ann. Appl. Probab., 4:390–413, 1994.
  * [47] D. Stanford, F. Avram, A. Badescu, L. Breuer, A. da Silva Soares, and G. Latouche.  Phase-type approximations to finite-time ruin probabilities in the Sparre Andersen and stationary renewal risk models.  ASTIN Bulletin, 35:131–144, 2005.
  * [48] W.-G. Wang, W.-C. Wang, and R.-C. Li.  Alternating-directional doubling algorithm for M-matrix algebraic Riccati equations.  SIAM J. Matrix Anal. Appl., 2012.