The morphing of fluid queues into Markov-modulated Brownian motion

Ramaswami showed recently that standard Brownian motion arises as the limit of a family of Markov-modulated linear fluid processes. We pursue this analysis with a fluid approximation for Markov-modulated Brownian motion. Furthermore, we prove that the stationary distribution of a Markov-modulated Brownian motion reflected at zero is the limit from the well-analyzed stationary distribution of approximating linear fluid processes. Key matrices in the limiting stationary distribution are shown to be solutions of a new quadratic equation, and we describe how this equation can be efficiently solved. Our results open the way to the analysis of more complex Markov-modulated processes.


Introduction
Our purpose is to construct and analyse a family of fluid queues converging to Markov-modulated Brownian motion (MMBM) with the intention of adapting, to the analysis of MMBM, tools and methods which have been developed in the context of fluid queues.
Fluid queues are two-dimensional processes {X(t), ϕ(t) : t ≥ 0}, where {ϕ(t) : * The authors had interesting discussions with V. Ramaswami over the questions examined in this paper. They thank three anonymous referees for their careful analysis and insightful suggestions, which helped streamline some of the proofs and improve the overall presentation of the paper. They also thank the Ministère de la Communauté française de Belgique for funding this research through the ARC grant AUWB-08/13-ULB 5, and acknowledge the financial support of the Australian Research Council through the Discovery Grant DP110101663.
t ≥ 0} is a continous-time Markov chain on a finite state space S, and c i for i ∈ M are arbitrary real numbers. These are also known as Markovmodulated linear fluid processes, with X referred to as the fluid level and ϕ as the phase: during intervals of time where the phase ϕ remains in some state i ∈ M, the fluid level varies linearly at the rate c i . The associated process reflected at zero is denoted by { X(t), ϕ(t) : t ≥ 0}, where assuming that X(0) = 0. Also referred to as first-order fluid processes, these stochastic models are useful when the relevant rates of change can be welldescribed by their first moments.
Markov-modulated Brownian motion (MMBM), or the family of second-order fluid processes, takes into account first and second moments of the rates of change. In particular, the fluid level Y (t) of a Markov-modulated Brownian motion {Y (t), κ(t) : t ≥ 0} is a Brownian motion with drift µ i and variance σ 2 i during time intervals where κ(t) = i; one sometimes writes that where {W (t)} is a standard Brownian motion.
Three papers appeared in close succession on the stationary distribution of MMBM reflected at zero: Rogers [25], Asmussen [4], and Karandikar and Kulkarni [17]. The focus in the third paper is on solving partial differential equations, and it is not of further concern to us in the present paper. In [4,25], on the other hand, the authors obtain the stationary distribution in a form which is suitable for calculations with linear algebraic procedures. These results crucially depend on the technique of reversing time, a method already used in Loynes [20] whereby the stationary distribution of the process is obtained from the distribution of the maximum of a random walk with negative drift. More recent work for obtaining the stationary distribution of Markov-modulated Lévy processes with reflecting boundaries, as in Asmussen and Kella [5], Ivanovs [16], D'Auria et al. [11], and D'Auria and Kella [12], also uses the reverse-time approach.
For fluid processes without a Brownian component, another line of investigation was open in Ramaswami [23], based on renewal-type arguments similar to the ones used in the analysis of quasi-birth-and-death processes (Neuts [21,Chapter 3], Latouche and Ramaswami [19,Chapter 6]). This eventually led, in addition to interesting algorithmic procedures, to theoretical developments for fluid processes in finite time, and for systems with more complex interactions between phase and level. There, the flow of time is not reversed and this creates a significant difference in the methods of analysis.
Our intention is to establish a link between the results for fluid queues and those for MMBM. In Section 2, we extend the argument from Ramaswami [24] and define a parameterised family of linear fluid processes that converge weakly, as the parameter tends to infinity, to a Markov-modulated Brownian motion. Ahn and Ramaswami [1] are independently using matrix-analytic methods to analyze Markov-modulated Brownian motions, with different approaches to ours. We determine in Sections 3 the limiting structure of key matrices and quadratic matrix equations, and we establish the connection between the stationary distribution so obtained, and the one which follows from the time-reversed approach. We present in Section 4 a computational procedure for solving the quadratic equation efficiently.

Markov-modulated Brownian motion
We show here that a family of linear fluid processes converges weakly to a Markov-modulated Brownian motion {Y (t), κ(t) : t ≥ 0}, where the phase process κ is a Markov chain with state space M = {1, . . . , m}, and Y is a Brownian motion with drift µ i and variance σ 2 i whenever κ(t) = i ∈ M. We denote by D the drift matrix diag(µ 1 , . . . , µ m ), by V the variance matrix diag(σ 2 1 , . . . , σ 2 m ), and by Q the generator of κ, and we assume that Q is irreducible.
Formally, for t ≥ 0 the process Y (t) can be defined recursively as where Y (0) = 0, the random variable T is the last jump epoch of κ before t (T = 0 if there has yet to be a jump), the process Y i is a Brownian motion with mean µ i and variance σ 2 i , and 1 {·} denotes the indicator function. The processes Y i for al i ∈ M and ϕ are assumed to be mutually independent.
Intuitively speaking, we duplicate the state space M in the Markov-modulated Brownian motion {Y (t), κ(t)} and the auxiliary process β λ (t) keeps track of which copy is in use. Note that for λ sufficiently large, the phases in the copy with β λ (t) = 1 have all positive rates while the phases in the other copy have all negative rates. With this construction, we show that the conditional moment generating function of {L λ (t), ϕ λ (t)} converges to that of {Y (t), κ(t)}. Remark 2.3. Based on the recursive definition (1) of Y , an alternative interpretation for our fluid-based approximation is that for each phase i ∈ M we approximate the process Y i by a two-phase fluid process {L i λ (t), β i λ (t)} the phase β i λ is a Markov chain with state space S i = {1, 2} and generator and the fluid rate matrix C i λ is given by the approximating processes being independent.
Denote by ∆ Y (s) the m × m Laplace matrix exponent of {Y (t), κ(t)}, bỹ ∆ λ (s) the 2m × 2m Laplace matrix exponent of {L λ (t), β λ (t), ϕ λ (t)}, and by ∆ λ (s) the m × m Laplace matrix exponent of {L λ (t), ϕ λ (t)}. These matrices are such that By Asmussen and Kella [5], the Laplace matrix exponent of a Markov-modulated Lévy process {Z(t), ξ(t)} with jumps is given by where φ i (s) is the Laplace exponent of an unmodulated Lévy process with parameters defined for phase i ∈ {1, . . . , p}, Q is the phase-transition matrix of ξ(t), R is the matrix with components [R(s)] ij = E[e sWij ], which are the Laplace transforms of the jumps W ij for Z(t) when ξ(t) moves from i to j, and • indicates the Hadamard product.
As the Laplace exponent of an unmodulated Brownian motion with drift µ and variance σ 2 is given by µs + σ 2 s 2 /2, one can verify that where 1 is an appropriately-sized column vector of ones. The next lemma gives a technical property of the matrix exponential, it is used in the proof of Theorem 2.5.

The matrix H(t) is the solution of
Proof. We decompose S as the sum S = S A + S E , with By Higham [15,Equations (10:13) and (10:40)], we obtain and e SAt = e S11t t 0 e S11(t−u) S 12 e S22u du 0 e S22t Thus, which proves (4).
Proof. We proceed in three steps. First, we observe that where e = (1, −1) T and with The proof of (6) is by induction: that equation trivially holds for k = 0, with A 0 = I and B 0 = 0 and, if it also holds for a given value of k, then we easily verify that∆ k+1 Equation (7) readily follows.
To simplify the notation, we define H λ (t) = e ∆ λ (s)t for the remainder of the proof. By (3), we have Integrating by parts the inner integral, we find This completes the proof.
The Laplace matrix exponent uniquely characterizes the finite-dimensional distributions of the process and therefore Theorem 2.5 implies the following result.
Proof. In this proof, we use alternative interpretation of our fluid-based approximation as outlined in Remark 2.3.
The proof of [24,Theorem 5] shows that for each i ∈ M the family of marginal processes {L i λ } is tight on any compact interval of [0, ∞). By [26, Corollary 7], we can extend this result to show that the family of marginal processes As the family of {L i λ } converges weakly to the Brownian motion {Y i } with µ i and σ i , so do the processes {sup s<t |L i λ (t)|} to the corresponding supremum process {sup s<t |Y i (t)|}, by the continuous mapping theorem. Thus, the family As L λ is stochastically dominated by L λ , the tightness property of the family of {L λ } follows from the tightness property of the family {sup The next theorem follows from Corollary 2.6 and Theorem 2.7.

Stationary distribution
We consider again the Markov-modulated Brownian motion {Y (t), κ(t) : t ≥ 0} described in Section 2, but with a reflection at level zero. The reflected process For notational convenience, we define ε = 1/ √ λ. With this, the reflected fluid processes is written as { L ε (t), β ε (t), ϕ ε (t) : t ≥ 0} and our purpose is to show that the stationary distribution of { Y (t), κ(t)} is the limit, as ε → 0, of the stationary distribution of the reflected fluid process { L ε (t), ϕ ε (t)}. We emphasize that the processes { L λ (t), ϕ λ (t)} and { L ε (t), ϕ ε (t)} are the same. The change in subscripts only reflects the notational change in our perturbation parameter.
Assumption 3.1. The mean drift αD1 is negative, so that all reflected processes are positive recurrent.
This assumption ensures the existence of Θ −1 , which we need later on.
The following result is a direct corollary of Theorem 2.8, by the Skorokhod mapping theorem.
and the associated stationary density vector by π ε (x), with components , 2} and i ∈ M, and we partition the generator and the fluid rate matrices as We assume that ε is sufficiently small that the diagonal elements of C + ε are all positive, and those of C − ε are all negative, and we write |C − ε | for the matrix of absolute values of the elements of C − ε . Let ξ + ε (x) = inf{t < ∞ : L ε (t) > x} and ξ − ε (x) = inf{t < ∞ : L ε (t) < x} be the first passage times to the level x, respectively from below and from above, of the unbounded process L ε . A key component of the stationary distribution of { L ε (t), β ε (t), ϕ ε (t)} is the matrix Ψ ε of first passage probability from above, that is, where and ζ ε is the unique solution of Probabilistically, ζ ε is up to a multiplicative constant the stationary distribution of the process censored at level 0, and e Kεx is the matrix of expected number of crossings of level x in the various phases (1, i), starting from level 0, before the first return to level 0. In view of (11), we need to analyze Ψ ε , K ε and ζ ε as ε → 0, and it is obvious from (12) and (13,14) that we should focus on the matrix Ψ ε first. The next lemma is the key to our analysis. One might expect (15,16) to have a simple proof but the one we have is lengthy and tedious. We place it in Appendix to preserve the flow of the paper. Lemma 3.6 and Theorem 3.7 easily follow.
where Θ −1 Ψ 1 and −Θ −1 Ψ * 1 are both solutions to and are irreducible. Furthermore, the roots θ 1 , θ 2 , . . . , θ 2m of the polynomial associated to (17), numbered in increasing order of their real parts, satisfy the inequalities Finally, Θ −1 Ψ 1 has one eigenvalue equal to zero and m − 1 eigenvalues with strictly negative real parts, and it is the unique such solution; −Θ −1 Ψ * 1 has m eigenvalues with strictly positive real parts, and is the unique such solution.
By Lemma 3.4, the matrices Ψ 1 and Ψ * 1 are uniquely identified through (17). We now turn to the matrix K ε , and, to complete the picture, to the matrix K * ε defined as Lemma 3.6. For ε ≥ 0, The matrices −K 0 and K * 0 are solutions of the equation and are irreducible. Furthermore, K 0 has m eigenvalues with strictly negative real parts, and it is the unique such solution; K * 0 has one eigenvalue equal to zero and m − 1 eigenvalues with strictly negative real parts, and is the unique such solution.
Proof. We write (12) as by (40), (41) and (15), which proves (20); equation (21) follows in a similar manner. It is easy to verify that K 0 and −K * 0 both satisfy (22), of which the associated polynomial is by (52), after some simple manipulations. This, together with Lemma 3.4, shows that the eigenvalues of −K 0 are the roots θ m+1 , . . . , θ 2m of γ(z) with strictly positive real parts. Finally, using (53) we write and conclude that the eigenvalues of K * 0 are the roots θ 1 to θ m . Finally, as Ψ 1 and Ψ * 1 are irreducible, and Θ, V , and D are diagonal matrices, we conclude that K 0 and K * 0 are irreducible, and this completes the proof. The next theorem states that the limit, as λ → ∞, of the stationary distributions of the approximating fluid processes is indeed the stationary distribution of the limiting process {Y (t), κ(t)}. We prove this result by showing that the limiting distribution (23)  Theorem 3.7. The limiting distribution of { L λ (t), ϕ λ (t)} converges, as λ goes to infinity, to the stationary distribution of {Y (t), κ(t)}, and is given by where Proof. The solution of (13) is of the form ζ ε = ζ 0 +εζ 1 +o(ε) ([18, Theorem 5.4]) and (14) becomes We equate the coefficients of 1/ε on both sides of (27), using (15,20) and we find that (10) implies that ζ ε ≥ 0 and, by continuity, ζ 0 ≥ 0. Furthermore, Θ −1 is a diagonal matrix with strictly positive diagonal. Finally, K 0 is irreducible and has eigenvalues with strictly negative real part by Lemma 3.6, so that ∞ 0 e K0u du converges to −K −1 0 and is strictly positive. This implies that ζ 0 = 0, which proves (24). Using ζ 0 = 0, and equating coefficients of ε on both sides of (13) gives (25), while equating the coefficients of ε 0 on both sides of (27) leads to (26).
To verify that the limiting distribution in Theorem 3.7 is the stationary density vector g(x) of the Markov-modulated Brownian motion, we use Asmussen [4]. By [4, Theorem 2.1 and Corollary 4.1], where α is the stationary distribution vector of Q, ∆ α = diag(α), and Λ is a defective generator matrix satisfying Define Z = ∆ 1/α Λ T ∆ α and rewrite (29) as By (30), we find that (1/2)Z 2 V − ZD + Q = 0 and that Θ −1 ZΘ is a solution of (22). It is similar to Z and so to Λ T , therefore, the eigenvalues of Θ −1 ZΘ all have strictly negative real parts and we have Substituting (32) into (23) gives lim ε→0 π ε (1 ⊗ I) = 2ζ 1 Θ −1 e Zx . Finally, it is straightforward to verify that ζ 1 = −α(Θ −1 D + (1/2)ΘΨ 1 Θ −1 ), and consequently Remark 3.8. An alternative way to show that the stationary distribution of the approximating fluid process { L λ (t), ϕ λ (t)} converges, as λ → ∞, to the stationary distribution of the limiting Markov-modulated Brownian motion {Y (t), κ(t)} is via the maximum representation of the relevant processes. Asmussen [4] derives the stationary distribution, both for fluid queues and for the MMBM in this manner, linking these to the distribution of the maximum of the time-reversed process. Following the arguments in Enikeeva et al. [13] and in Stenflo [22], one might show that there is continuity of the maximum distributions of the backward processes, as λ → ∞, and consequently obtain the continuity of stationary distributions. This would lead to a time reversal-based proof of convergence. As stated in the introduction, we aim at following the forward-time approach and so obtain a different representation of the stationary distribution. In addition, we obtain limiting properties for key matrices, and these results will be proved useful in future work.

Computational procedure
Theorem 3.7 indicates that the matrix Ψ 1 is the central ingredient in evaluating the stationary distribution of the Markov-modulated Brownian motion {Y (t), ϕ(t)}. We describe here how to use the splitting property (18) in numerically solving for Ψ 1 and Ψ * 1 . Bini and Gemignani [6] consider quadratic matrix equations C +AX +BX 2 = 0 where the roots of the associated polynomial det(C + zA + z 2 B) are split by a circle in C, half being inside the disk and half outside. The problem in [6] is to find the minimal solution, that is, the solution matrix with all eigenvalues inside the disk.
In our case, the roots are split between the negative and the positive halfplanes and we need to apply some transformation, such as the one described in Bini et al. [8] and based on the inverse Möbius mapping [3, Chapter 2.1] This inverse mapping applies the open unit disk |z| < 1 onto the negative halfplane C − , the unit circle |z| = 1 minus the point z = −1 onto the imaginary axis C 0 , the outside |z| > 1 of the closed unit disk onto the positive half-plane C + , and the imaginary axis C 0 onto the unit circle |w| = 1 minus the point w = 1.
The roots of det(H(z)) are given by ω i = w −1 (τ i ) = (1 + τ i )/(1 − τ i ) for i = 1, . . . , 2m, and satisfy the splitting property We note that Bini and Geminiani [6] requires |ω i | > 0 for all i, but as seen in Several iterative algorithms to compute Z 0 are discussed in [6]. Some have superlinear convergence, such as Cyclic reduction [7], Logarithmic reduction [19,Chapter 8], subspace iteration [2] or Graeffe iteration [7]. This means that the approximation error at the ith iteration is O(σ 2 i ) with σ = 1/|ω m+1 | < 1. These algorithms are globally convergent but they are not self-correcting. Since the coefficient matrices in P (X) are of mixed signs, one might prefer the algorithm developed in [6]: it is self-correcting and the approximation error is O(σ i2 k ) for arbitrary k, which makes it arbitrarily fast.
As for Ψ * 1 , if we define Z 1 as the root of H(Z) such that all of its eigenvalues are outside the closed unit disk, then W (Z 1 ) has all its eigenvalues in the halfplane C + , and Θ −1 Ψ * 1 = −W (Z 1 ). The algorithms in [6], however, do not seem to be well adapted to the computation of Z 1 , and we suggest to use a different transformation, in order to bring the eigenvalues of −Θ −1 Ψ * 1 inside the unit disk. This transformation is based on Möbius' mapping [3, Chapter 2.1] Appendix A: Proof of Lemma 3.4 The proof goes in four main steps. The matrix Ψ ε is the stochastic (or substochastic) solution of the Riccati equation [25]), equation that we write as For ε = 0, we see that F 0 (I) = 0. It is tempting to invoke the Implicit Function Theorem and claim that Ψ ε is an analytic function of ε in a neighborhood of ε = 0. Unfortunately, the operator ∂/∂XF ε (X) is singular at the point (ε = 0, X = I), the Implicit Function Theorem does not apply, and we follow a longer, more tortuous path.
Now, we observe that Ψ ε belongs to the compact set {M : M ≥ 0, M 1 ≤ 1} of (sub)stochastic matrices; therefore, for every sequence {Ψ ε } ε→0 there exist subsequences that converge. LetΨ be the limit of one such convergent subsequence, and {ε i } i=1,2,... be a sequence such that ε i → 0 and Ψ εi →Ψ as i → ∞. Since G(ε i ) remains bounded as i → ∞, necessarily = 0, and thusΨ = I. This follows from the facts that Θ −1 (I−Ψ) is a nilpotent matrix, that the trace of every nilpotent matrix is zero, and thatΨ is a (sub)stochastic matrix while Θ −1 is a strictly positive diagonal matrix.
All convergent subsequences having the same limit, the conclusion is that Ψ ε converges to I as ε → 0, and (38) follows. The proof of (39) is by analogous arguments.
Case 2: (1/ε)Φ ε is unbounded in any neighborhood of ε = 0. Then, there exists a sequence {ε k } k=1,2... such that ε k → 0 and max ij |Φ ε k | ij /ε k → ∞. In this case, we may write (1/ε)Φ ε = u ε B ε , where u ε is a scalar function such that lim k→∞ u ε k = ∞ while B ε k remains bounded and does not converge to zero: since Φ ε is an irreducible generator, its maximum element is on the diagonal and we take u ε = max j |Φ ε | jj /ε, then B ε is an irreducible generator with at least one diagonal element equal to −1, and with |B ij | ≤ 1 for all i and j.
Next, for ε in the sequence {ε k }, we replace (1/ε)Φ ε in (46) by u ε B ε and divide both sides of the equation by u 2 ε to obtain This implies that lim Now, take any converging subsequence of B ε and denote its limit as B. By construction, the trace of Θ −1 B ε is at most equal to min j (−σ −1 j ) < 0, independently of ε. Thus, the trace of Θ −1 B is strictly negative, the matrix Θ −1 B is not nilpotent, and Θ −1 BΘ −1 B = 0, which contradicts (49).
The polynomial det Γ(z) has at most 2m roots, which concludes the proof. Now we are ready to conclude. By the properties of the roots of det Γ(z) given in Proposition A.4, (17) has one unique solution suitable for the role ofΨ 1 and another unique solution suitable for the role ofΨ * 1 . Consequently, all convergent subsequences give the same limit Ψ 1 for (1/ε k )Φ ε k , and Ψ * 1 for (1/ε k )Φ * ε k .