TEXT SIZE

search for



CrossRef (0)
Application of GTH-like algorithm to Markov modulated Brownian motion with jumps
Communications for Statistical Applications and Methods 2021;28:477-491
Published online September 30, 2021
© 2021 Korean Statistical Society.

Sung-Chul Honga, Soohan Ahn1,a

aDepartment of Statistics, University of Seoul, Korea
Correspondence to: 1 Department of Statistics, University of Seoul, 163 Seoulsiripdaero, Dongdaemun-gu, Seoul 02504, South Korea. E-mail: sahn@uos.ac.kr
This work was supported by the 2018 Research Fund of the University of Seoul.
Received February 15, 2021; Revised April 12, 2021; Accepted April 13, 2021.
 Abstract
The Markov modulated Brownian motion is a substantial generalization of the classical Brownian Motion. On the other hand, the Markovian arrival process (MAP) is a point process whose family is dense for any stochastic point process and is used to approximate complex stochastic counting processes. In this paper, we consider a superposition of the Markov modulated Brownian motion (MMBM) and the Markovian arrival process of jumps which are distributed as the bilateral ph-type distribution, the class of which is also dense in the space of distribution functions defined on the whole real line. In the model, we assume that the inter-arrival times of the MAP depend on the underlying Markov process of the MMBM. One of the subjects of this paper is introducing how to obtain the first passage probabilities of the superposed process using a stochastic doubling algorithm designed for getting the minimal solution of a nonsymmetric algebraic Riccatti equation. The other is to provide eigenvalue and eigenvector results on the superposed process to make it possible to apply the GTH-like algorithm, which improves the accuracy of the doubling algorithm.
Keywords : Markov modulated Brownian motion, Markovian arrival process, first passage time, doubling algorithm, nonsymmetric algebraic Riccati equation, GTH-like algorithm
1. Introduction

This paper is concerned with the superposed process of the Markov modulated Brownian motion (MMBM) and Markovian arrival process (MAP), and a technique to enhance the performance of an algorithm for computation of the first passage probabilities of the superposed process.

The MMBM, a generalization of the Brownian motion, is a bivariate Markov process of which the drift and diffusion parameters vary according to the states of a underlying Markov process. On the other hand, the MAP is a point process whose jump epochs are also governed by an underlying Markov process. Since the family of MAP is dense for any stochastic point process, it is used to approximate complex stochastic counting processes. Concerning the sizes of jumps of the MAP, we assume that they are mutually independent and distributed as the bilateral ph-type (BPH) distributions with their parameters depending on the states of the underlying Markov processes of both the MMBM and MAP. We note that the bilateral ph-type distribution, introduced by Ahn and Ramaswami (2005), can be represented as a mixture of two independent ph-type distributions, and the class of bilateral ph-type distributions is also dense in the space of distribution functions defined on the whole real line.

Then, the superposition of the MMBM and MAP becomes a generalized process that can be used for modelling complex systems in numerous academic fields including queues, finance, and insurance risk theories. For example, of applications, see Belhaj (2010), Jiang and Pistorius (2008), Asmussen (2003), and the references therein. We call this superposed process the generalized Markov modulated Brownian motion and use G-MMBM for abbreviation. In relation to G-MMBM, we deal with two main subjects in this paper. The first one is to provide Markovian representation of the G-MMBM so that we can obtain the first passage probabilities of the G-MMBM using the approaches of Ahn (2016), in which they show that the first passage probabilities can be achieved by using the minimal non-negative solution of a nonsymmetric algebraic Riccati equation (NARE) of the form

AZ+ZB+ZCZ+D=0.

The other is to show eigenvalue and eigenvector results so that we can use the GTH-like algorithm introduced by Alfa et al. (2001), which enhances the accuracy of the algorithms for computing the minimal nonnegative solution of the NARE in relation to the first passage probabilities of the G-MMBM.

The remainder of this paper is organized as follows. In Section 2, we introduce the G-MMBM and also the relation between the first passage probabilities of the process and the minimal nonnegative solution of a nonsymmetric algebraic Riccati equation. In Section 3, we provide parameter representation of the G-MMBM which can be obtained using embedding technique. The doubling algorithm and eigenvalue and eigenvector results on the G-MMBM to be used for GTH-like algorithm are introduced in Section 3. Finally, we illustrate the effects of the GTH-like algorithm by using numerical examples.

2. Description of the G-MMBM and related NARE

The G-MMBM is the superposed process of an MMBM and an MAP. We note that the MMBM consists of a modulating irreducible Markov process JB with a state space S(B)=Sb(B)Su(B)Sd(B), a drift row-vector μ(B)=(μb(B)μu(B)μd(B)), and a diffusion row-vector σ(B)=(σb(B)00); the MAP is modulated by an irreducible Markov process JM with a state space S(M) = {1, …, m} and its transition and jump rates depend on JB. The resulting G-MMBM is a two-dimensional Markov process and use (X, J) to denote the process, where we call X and J the level and the phase process of the G-MMBM, respectively. Here, the phase process J is the superposed process of JB and JM of which the state space S is represented as S = S bS uS d with Sb=Sb(B)S(M),Su=Su(B)S(M), and Sd=Sd(B)S(M). Then the level process X can be represented as

X(t)=a+0tμJ(u)du+0tσJ(u)dB(u)+0tdM(u),

where (i) when J(t) = (i, j) ∈ S, μ(i, j) = [μ(B)]i and σ(i, j) = [σ(B)]i; (ii) M is the jump process with the jump epochs being either the transition epochs of J, or the arrival epochs of the compound Poisson processes whose intensity depends on the state of J. Moreover, the sizes of the incurred jumps are mutually independent and distributed as the BPH distribution of the form

G(dx)=β+eB+x(-B+1)dxχ(x>0)+β-e-B-x(-B-1)dxχ(x<0),

in which nonnegative probability vectors β+ and β satisfy β+1 + β+1 = 1, and B+ and B are sub-stochastic matrices which have negative diagonal elements, nonnegative off-diagonal elements, and non-positive row-sums; (iii) and B = {B(t), t ≥ 0} is a standard Brownian motion that is independent of J and M.

In connection with the state spaces S, we denote the drift and diffusion row vectors by μ = (μi, iS ) and σ = (σi, iS ), respectively, and we define the sub-vectors σk and μk as σk = (σi, iS k) and μk = (μi, iS k) for k = b, u, d. In addition to the state space S, it is necessary to adopt jump-state spaces S+ and S modulating bilateral ph-type jumps, details of which are to be given in Section 3. For these state spaces, we define the drift and diffusion row vectors as μ+=1+,μ-=1-,σj+=0+, σj+ = 0+, and σ = 0 with the prime denoting the transpose operator. Note that, for k ∈ {b, u, d, +, −}, we use 0k, 1k, Ik, and 0k,l to denote the sk(the dimension of S k)-dimensional row vector of 0s, column vector of 1s, sk-dimensional identity matrix, and sk × sl-dimensional zero matrix, respectively. We also use I and 0 to denote the identity matrix and the zero matrix of appropriate dimension which can be identified from the context. Throughout the paper, for a given vector v, we use the notation Δv to denote the diagonal matrix with the elements of v on its diagonal.

2.1. Representation of the G-MMBM through the embedding technique and completed graph

Let (, ) denote the embedded process of (X, J) to be obtained by using an embedding technique, which replaces the positive and negative jumps in (X, J) with linear stretches having slopes +1 and −1 respectively in , and lets track the states (phases) in BPH jumps in addition to the states in J. In Figure 1, we illustrate how a sample path of (X, J) is transformed to its corresponding sample path of (, ). Then, the embedded process (, ) becomes an MMBM without a jump, in which its state space Se, diffusion vector σe, and drift vector μe are respectively given as Se = S bS uS+S dS, σe = (σb0u0+0d0), and μe = (μbμuμ+μdμ). We denote by T the infinitesimal generator of and represent it in a partitioned form as

T=(Ti,j,i,j=b,u,+,d,-)with T+,-=0and T-,+=0,

where Ti, j is a sub-matrix of T containing [T]lm, the (l,m)th element of T with lS i and mS j. We note that the restriction of T on S should be Q, and that (−[T]ii)−1, iS, represents the mean dwelling time of J in i until a jump or a transition to other states in S occurs, whereas (−[Q]ii)−1 is the mean dwelling time of J in i until a transition of J to other states in S. We let πe = (πbπuπ+πdπ) be the stationary probability vector of T satisfying πeT = 0 and πe1 = 1. Defining the average drift μ̃ of the embedded process as μ̃ = ∑iSe [πe]i[μe]i, we say that the G-MMBM is positively recurrent, transient, and null recurrent when μ̃ < 0, > 0, and = 0, respectively. We describe the detailed constructions of Se and T in Section 3.

All the information for analyzing the G-MMBM is contained in the parameter triple (T, μe,σe). So, when we call the parameter of the G-MMBM, it indicates the parameter triple. With regard to the level process , it is represented as

X˜(t)=a+0t[μe]J˜(u)du+0t[σe]J˜(u)dB(u).

For description of G-MMBM, it is useful to introduce a terminology called the completed graph of the process and its phase process to be defined for . We note that the completed graph is a subset of the plane ℝ2 containing (X(t), t) for all t and the line segment joining (X(t), t) and (X(t−), t) for all discontinuity points t. The completed graph, introduced in Whitt (1980), is useful for representing quantities in relation to a process with discontinuous sample paths. For more details, refer to Whitt (1980) and Ahn (2016, 2017). We illustrate a G-MMBM, its completed graph, and the phase process defined on the completed graph in Figure 2.

2.2. First passage time of the G-MMBM and a nonsymmetric algebraic Riccati equation

We define the first passage times of X to 0 as τ = inf{t > 0 : X(t) < 0}. In relation to τ, we define an probability matrix (a) such that its (i, j)th element, i, jSe, satisfies

[f^(a)]ij=P[J(0,τ)=jJ(a,0)=i],         a0.

Then, the (sb + sd + s)-dimensional squared sub-matrix bd−,bd(a) = ([ (a)]i j; i, jS bS dS) has the matrix exponential form eHa, where H is a sub-stochastic matrix.

Concerning H, Ahn (2016) showed that it can be obtained by using the minimal solution of a non-symmetric algebraic Riccati equation of the form

AZ+ZB+ZCZ+D=0,

where, with definition of the following relevant matrices

Λ=diag{-[T]ii,iSb},Σ=Δσb,D1=Σ-2Δμb,and D2=Σ-1(2Λ+Σ-2Δμb2)12,

the coefficient matrices satisfy

A=(D1-D22Σ-1Tb,u2Σ-1Tb,+0u,bΔμu-1Tu,uΔμu-1Tu,+0+,bT+,uT+,+),         C=(Σ-10b,u0b,+0d,bΔ-μd-1Td,uΔ-μd-1Td,+0-,bT-,u0-,+),B=(-D1-D20b,d0b,-Δ-μd-1Td,bΔ-μd-1Td,dΔ-μd-1Td,-T-,bT-,dT-,-),         D=(2Σ-1(Tb,b+Λ)2Σ-1Tb,d2Σ-1Tb,-Δμu-1Tu,bΔμu-1Tu,dΔμu-1Tu,-T+,bT+,d0+,-).

Furthermore, from the results of Section 2.1.2 of Bini et al. (2012), Liu and Xue (2012), and Theorem 3.1.1 and 3.1.2 of Ahn (2016), we can derive the following result:

Theorem 1

Let σ(·) denote the set of eigenvalues of a given squared matrix ·.

(a) The NARE (2.4) has a minimal nonnegative solution Z. Moreover, Z is the unique solution such that σ(B + CZ) ⊂ ℂ ∪ {0} and σ(A + ZC) ⊂ ℂ ∪ {0}.

(b) Let Z be the minimal non-negative solution of the NARE (2.4). ThenH satisfies the relation

H=B+CZ.
3. Representation of the superposed process as a G-MMBM

We introduce how to represent the G-MMBM which is the superposed process of an MMBM and Markovian arrival processes (MAP) with BPH-distributed jumps. Since the transition rates and jump rates of the MAP are assumed to be dependent on the process JB modulating the MMBM, we denote the MAP parameters by (Γ0JB(t),Γ1JB(t)). We also assume that MAP jumps are mutually independent and distributed as the BPH distribution of which the parameters depend on the state of the superposed Markov process J. We note that the representation form of the G-MMBM depend on dependence structures of the jumps on the Markov process J.

3.1. BPH jumps depending on only the Markov states before their corresponding jump epochs

In this section, we assume that the jump sizes are mutually independent and distributed as the BPH distribution, of which the parameters depend on the state of J just before their jump epochs. We begin with simple examples for illustration.

Example 1

We consider the superposition of a two-state MMBM and two-state MAP, where

  • (i) S(B) = {1, 2}, Q(B) = (qi j), μ(B) = (μ1μ2) and σ(B) = (σ1σ2),

  • (ii) two-state MAP(Γ0JB(t),Γ1JB(t)) with S(M) = {1, 2}, Γ0i=(γ0,kli;k,lS(M)) and Γ1i=(γ1,kli;k,lS(M)) for iS(B),

  • (iii) an MAP jump occurring at the transition epoch of JM from k to l while JB stays in iS(B) is assumed to be distributed as ph-type distribution of the form

    βkieBkix(-Bki1)χ(x>0).

Then, the infinitesimal generator of the embedded process of the superposition can be represented as

(q11+γ0;111γ0;121q120(γ1;111+γ1;121)β11000γ0;211q11+γ0;2210q120(γ1;211+γ1;221)β2100q210q22+γ0;112γ0;12200(γ1;112+γ1;122)β1200q21γ0;212q22+γ0;222000(γ1;212+γ1;222)β22(-B111)γ¯1;111(-B111)γ¯1;12100B11000(-B211)γ¯1;211(-B211)γ¯1;221000B210000(-B121)γ¯1;112(-B121)γ¯1;12200B12000(-B221)γ¯1;212(-B221)γ¯1;222000B22)

where γ¯1;kli=γ1;kli/(γ1;k1i+γ1;k2i) for iS(B) and k, lS(M). The drift and diffusion vectors are given as (μ1μ1μ2μ21111) and (σ1σ1σ2σ20 0 0 0), respectively.

Example 2

We consider the special case of Example 3.1. In addition to the assumption (i) and (ii) of Example 3.1, we assume that (ii)′ γ1:ki/S(M)γ1:ki=γ¯1:i for all kS(M); (iii) an MAP jump occurring at the transition epoch of JM from k to l while JB stays in iS(B) is assumed to be distributed as BPH distribution of the form

βki-eBi-x(-Bi-1)χ(x<0)+βki+eBi+x(-Bi+1)χ(x>0).

Then, the drift and diffusion vectors are given as (μ1μ1μ2μ211′−1′−1′) and (σ1σ1σ2σ20 0 0 0), respectively, and the infinitesimal generator of the embedded process can be represented as

(q11+γ0;111γ0;121q120(γ1;111+γ1;121)β11+0(γ1;211+γ1;221)β11-0γ0;211q11+γ0;2210q12(γ1;211+γ1;221)β21+0(γ1;211+γ1;221)β21-0q210q22+γ0;112γ0;1220(γ1;112+γ1;122)β12+0(γ1;112+γ1;122)β12-0q21γ0;212q22+γ0;2220(γ1;112+γ1;122)β22+0(γ1;212+γ1;222)β22-(-B1+1)γ¯1;11(-B1+1)γ¯1;2100B1+00000(-B2+1)γ¯1;12(-B2+1)γ¯1;220B2+00(-B1-1)γ¯1;11(-B1-1)γ¯1;210000B1-000(-B2-1)γ¯1;12(-B2-1)γ¯1;22000B2-)
Example 3

We consider the superposition of a two-state MMBM and inhomogeneous Poisson process, where,

  • (i) S(B) = {1, 2}, Q(B) = (qi j), μ(B) = (μ1μ2) and σ(B) = (σ1σ2),

  • (ii) Poisson process with intensity γ(JB(t)),

  • (iii) the size of a jump which occurs while JB stays in iS(B) is assumed to be distributed as BPH distribution of the form

    βi+eBi+x(-Bi+1)χ(x>0)+βi-eBi-x(-Bi-1)χ(x<0).

Then, the infinitesimal generator of the embedded process of the superposition can be represented as

(q11-γ(1)q12γ(1)β1+0γ(1)β1-0q21q22-γ(2)0γ(2)β2+0γ(2)β2-(-B1+1)0B1+0000(-B2+1)0B2+00(-B1-1)000B1-00(-B2-1)000B2-).

The drift and diffusion vectors of the embedded process are given as (μ1μ21′− 1′) and (σ1σ20 0), respectively.

Now, we consider the general case of the superposition. For the MAP jump occurring when J = (i, ), we let (βi+,Bi+,βi-,Bi-) denote the parameters of the BPH distribution. If we denote by Si+ and Si- the state spaces corresponding to Bi+ and Bi-, respectively, then S+=iS(B),S(M)Si+ and S-=iS(B),S(M)Si-. The drift and diffusion vectors μe=(μbμuμd1j+1j-) and σe = (σb0 0 0 0) satisfy μb=μb(B)1m,μu=μu(B)1m,μd=μd(B)1m and σb=σb(B)1m. Moreover, the infinitesimal generator T of the embedded process given in the following partitioned form

T=(TbbTbuTbdTb+Tb-TubTuuTudTu+Tu-TdbTduTddTd+Td-T+bT+uT+dT++0T-bT-uT-d0T--),

in which, with χ(·) denoting the indicator function and matrices aijk for k, ∈ {b, u, d, +, −} and iS k, jS being defined as

aijkk=qijIm+χ(i=j)Γ0iand aijk=qijImfor k,

the submatrices satisfy that

Tbb=(aijbb,i,jSb),         Tbu=(aijbu,iSb,jSu),         Tbd=(aijbd,iSb,jSd),Tuu=(aijuu,i,jSu),         Tub=(aijub,iSu,jSb),         Tud=(aijuu,iSu,jSd),Tdd=(aijdd,i,jSd),         Tdb=(aijdb,iSd,jSb),         Tdu=(aijdu,iSd,jSu),

with γ1;ki=lSMγ1;kli for iS and matrices bik+ and bik- for k ∈ {b, u, d}, iS k being defined as

bik+=diag {γ1;jiβji+,jS(M)}   and bik-=diag {γ1;jiβji-,jS(M)},Tb+=(diag {bib+,iSb}00),Tu+=(0diag {biu+,iSu}0),Td+=(00diag {bid+,iSd}),Tb-=(diag {bib-,iSb}00),Tu-=(0diag {biu-,iSu}0),Td-=(00diag {bid-,iSd}),

with matrices ci+ and ci- for ∈ {b, u, d}, iS being defined as

ci+=diag {(-Bji+1)(γ1;ji)-1,jS(M)}Γ1iand ci-=diag {(-Bji-1)(γ1;ji)-1,jS(M)}Γ1i,T+b=(diag {ci+b,iSb}00),T+u=(0diag {ci+u,iSu}0),T+d=(00diag {ci+d,iSd}),T-b=(diag {ci+b,iSb}00),T-u=(0diag {ci-u,iSu}0),T-d=(00diag {ci-d,iSd}),

and with di++=diag{Bki+,kS(M)} and di--=diag{Bki-,kS(M)} for iS,

T++=(diag {di++,iSb}000diag {di++,iSu}000diag {di++,iSd}),T--=(diag {di--,iSb}000diag {di--,iSu}000diag {di--,iSd}).

3.2. BPH jumps depending on both the Markov states before and after their corresponding jump epochs

In this section, we consider the superposition in which the MAP jumps depend on the states of J both before and after the jump epoch of the process. We also take a simple example for illustration of the structure of the resulting embedded process.

Example 4

Consider the superposition of a two-state MMBM and two-state MAP, where,

  • (i) S(B) = {1, 2}, Q(B) = (qi j), μ(B) = (μ1μ2) and σ(B) = (σ1σ2),

  • (ii) two-state MAP(Γ0JB(t),Γ1JB(t)) with S(M) = {1, 2}, Γ0i=(γ0,kli;k,lS(M)) and Γ0i=(γ0,kli;k,lS(M)) for iS(B),

  • (iii) an MAP jump occurring at the transition epoch of J from (i, k) to (i, ) is assumed to be distributed as ph-type distribution of the form

    βkieBkix(-Bki1)χ(x>0).

Then, the infinitesimal generator of the embedded process of the superposition can be represented as

(q11+γ0;111γ0;121q120γ1;111β111γ1;121β121000000γ0;211q11+γ0;2210q1200γ1;211β211γ1;221β2210000q210q22+γ0;112γ0;1220000γ1;112β112γ1;122β122000q21γ0;212q22+γ0;222000000γ1;212β212γ1;222β222-B1111000B11100000000-B1211000B121000000-B121100000B211000000-B22110000B221000000-B112100000B112000000-B122100000B1220000-B21210000000B2120000-B22210000000B222)

The drift and diffusion vectors are given as (μ1μ1μ2μ211111111′) and (σ1σ1σ2σ20 0 0 0 0 0 0 0), respectively.

Now, we consider the general case of the superposition. For the MAP jump which occurs at the transition epoch of J from (i, j) to (i, k), we let (βjki+,Bjki+,βjki-,Bjki-) be the parameters of the BPH distribution. If we denote by Sjki+ and Sjki- the state space corresponding to Bjki+ and Bjki-, respectively, then S+=iS(B)jS(M)kS(M)Sjki+, and S-=iS(B)jS(M)kS(M)Sjki-. The drift and diffusion vectors μe=(μbμuμd1+1-) and σe = (σb0 0 0 0) satisfy μb=μb(B)1m,μu=μu(B)1m,μd=μd(B)1m and σb=σb(B)1m. Moreover, the infinitesimal generator T of the embedded process given in the following partitioned form

(TbbTbuTbdTb+Tb-TubTuuTudTu+Tu-TdbTduTddTd+Td-T+bT+uT+dT++0T-bT-uT-d0T--),

in which, with matrices aijk for k, ∈ {b, u, d, +, −} and iSk, jS being defined as

aijkk=qijIm+χ(i=j)Γ0iand aijk=qijImfor k,

the submatrices satisfy that,

Tbb=(aijbb,i,jSb),         Tbu=(aijbu,iSb,jSu),         Tbd=(aijbd,iSb,jSd),Tuu=(aijuu,i,jSu),         Tub=(aijub,iSu,jSb),         Tud=(aijud,iSu,jSd),Tdd=(aijdd,i,jSd),         Tdb=(aijdb,iSd,jSb),         Tdu=(aijdu,iSd,jSu),

with matrices bi+ and bi- for ∈ {b, u, d}, iS being defined as

bi+=bi+=diag {(γ1;k1iβk1i+,,γ1;kmiβkmi+),kSM}   and bi-=diag {(γ1;k1iβk1i-,,γ1;kmiβkmi-),kSM},Tb+=(diag {bib+,iSb}00),Tu+=(0diag {biu+,iSu}0),Td+=(00diag {bid+,iSd}),Tb-=(diag {bib-,iSb}00),Tu-=(0diag {biu-,iSu}0),Td-=(00diag {bid-,iSd}),

with matrices ci+ and ci+ for ∈ {b, u, d}, iS being defined as

ci+={diag {-B1,ki+1,kS(M)}diag {-Bm,ki+1,kS(M)}}and ci-={diag {-B1,ki-1,kS(M)}diag -{-Bm,ki-1,kS(M)}},T+b=(diag {ci+b,iSb}00),T+u=(0diag {ci+u,iSu}0),T+d=(00diag {ci+d,iSd}),T-b=(diag {ci-b,iSb}00),T-u=(0diag {ci-u,iSu}0),T-d=(00diag {ci-d,iSd}),

and with di++=diag{diag{Bk,li+,lS(M)},kS(M)} and di--=diag{diag{Bk,li-,lS(M)},kS(M)} for iSbSuSd,

T++=(diag {di++,iSb}000diag {di++,iSu}000diag {di++,iSd}),T--=(diag {di--,iSb}000diag {di--,iSu}000diag {di--,iSd}).
4. The structure-preserving doubling algorithm and GTH-like algorithm

We first introduce some relevant definitions and notations. For any matrices A, B ∈ ℝm×n, we write AB(A > B) if [A]i j ≥ [A]i j([A]i j > [B]i j for all i, j, where [A]i j denotes the (i, j)th element of A. A real square matrix A is called a Z-matrix if all its off-diagonal elements are non-positive. Any Z-matrix A can be written as sIB with B ≥ 0. A Z-matrix A is called an M-matrix if sρ(B), where ρ(·) is the spectral radius; it is called a singular M-matrix if s = ρ(B) and a non-singular M-matrix if s > ρ(B).

4.1. Structure-preserving doubling algorithm

We recall the parameter triple (T, μe,σe) of the G-MMBM given in the previous section and the stationary probability vector πe of the infinitesimal generator. Then the first passage probability of the G-MMBM can be obtained using the minimal solution of the NARE AZ + ZB + ZCZ + D = 0 given in (2.4) and (2.5) as mentioned in Section 2.2. In connection with the NARE (2.4), we define an M-matrix M as

M=(-B-C-D-A),

which is irreducible by the irreducibility of JB and JM. Moreover, M is a singular matrix.

To get the nonnegative minimal solution of the NARE, we can use the structure-preserving doubling algorithm (SDA) presented in Guo et al. (2006), which is as follow,

Structure-preserving doubling (SDA) algorithm for an NARE AZ + ZB + ZCZ + D = 0
  • Step 1. Choose γ ≥ max{−[A]ii, −[B]ii, iS } and set (E0G0H0F0)=(γI-B-C-DγI-A)-1(γI+BCDγI+A)

  • Step 2. Ek+1 = Ek(IGkHk)−1Ek; Fk+1 = Fk(IHkGk)−1Fk,

    Gk+1 = Gk + Ek(IGkHk)−1GkFk; Hk+1 = Hk + Fk(IHkGk)−1HkEk,

  • Step 3. Z = H,

The convergence results of the SDA when M is an irreducible M-matrix is given in the following theorem, which can be found in Theorem 3.1 and 3.2 of Wang et al. (2012) and Theorem 3.1 of Xue and Li (2017). For the theorem, we define the Cayley transform of R = −BCX and S = −ADY such that

Rγ=(R+γI)-1(R-γI)and Sγ=(S+γI)-1(S-γI).
Theorem 2

Note that {Ek}, {Gk}, {Hk}, and {Fk} denote the matrices in the kth iteration of the SDA.

(a) For all k ≥ 0, the matrices Ek, Gk, Hk, and Fk are non-negative. Moreover, for all k ≥ 0, IGkHk and IHkGk are non-singular M-matrices.

(b) In the positive recurrent case, then ρ(Rγ) = 1 and ρ(Sγ) < 1. Furthermore, {Hk} converges to X quadratically with,

lim supkHk-X2kρ(Sγ),

{Fk} converges to0quadratically withlim supkFk2kρ(Sγ), and {Ek} is bounded. The notation ||A|| denotes the maximum of the absolute values of the elements in a Matrix A.

(c) In the transient case, then ρ(Rγ) < 1 and ρ(S γ) = 1. Furthermore, {Hk} converges to X quadratically with

lim supkHk-X2kρ(Rγ),

{Ek} converges to0quadratically withlim supkEk2kρ(Rγ), and {Fk} is bounded.

(d) In the null recurrent case, then ρ(Rγ) = 1 and ρ(S γ) = 1. In this case, {Hk} converges to X and {Ek}, {Fk} are bounded.

As we can see in SDA algorithm, the initialization phase and the iteration kernel of the SDA involve inverting non-singular M-matrices. In particular, in the null recurrent case, the M-matrices, although non-singular, move towards singular M-matrices at convergence. These inversions are causes of concerns on entry-wise relative accuracy of the computed minimal nonnegative solution. Fortunately, a non-singular M-matrices can be inverted by the GTH-like algorithm of Alfa et al. (2001) to almost full entry-wise relative accuracy, which enhances the entry-wise relative accuracy of the computed minimal nonnegative solution.

4.2. GTH-like algorithm

The GTH-like algorithm, proposed by Alfa et al. (2001), is used for solving linear systems Ax = b when A is an M-matrix satisfying Av = c with v > 0 and c ≥ 0. It is a subtraction-free algorithm and compute the inverse of A with relative errors in the magnitude of the machine precision. For the details on the algorithm, refer to (2001). When it holds that Av = c with v > 0 and c ≥ 0, we say that the triplet representation of A is available and represent as A = (NA, v, c) with NA = diag(A) − A.

Proposition 1

Recall thatΛ = diag{−[Tb,b]i,i, iSb}. Then, the left and right eigenvectorsuandvcorresponding to the eigenvalue 0 of M are given as

u:=(u1u2)=([(D1+D2)-1ΛπbΔ-μdπdπ-][0.5ΣπbΔμuπuπ+])andv:=(v1v2)=([1b1d1-][(D1+D2)Σ1b1u1+]),

which are unique up to a scalar multiple.

Proof

We can prove by showing that uM = 0 and Mv = 0. Recall that

A=(D1-D22Σ-1Tb,u2Σ-1Tb,+0u,bΔμu-1Tu,uΔμu-1Tu,+0+,bT+,uT+,+),         C=(Σ-10b,u0b,+0d,bΔ-μd-1Td,uΔ-μd-1Td,+0-,bT-,u0-,+),B=(-D1-D20b,d0b,-Δ-μd-1Td,bΔ-μd-1Td,dΔ-μd-1Td,-T-,bT-,dT-,-),         D=(2Σ-1(Tb,b+Λ)2Σ-1Tb,d2Σ-1Tb,-Δμu-1Tu,bΔμu-1Tu,dΔμu-1Tu,-T+,bT+,d0+,-).

Note that uM=(-u1B-u2D-u1C-u2A). If we temporarily -u1B-u2D=(zbzdz-), then

zb=πbΛ(D1+D2)-1(D1+D2)-πdTd,b-π-T-,b-πb(0.5Σ)(2Σ-1)(Tb,b+Λ)-πuTu,b-π+T+,b=-πbTb,b-πuTu,b-π+T+,b-πdTd,b-π-T-,b=0.

Similarly, we can show that

zd=-πbTb,d-πuTu,d-π+T+,d-πdTd,d-π-T-,d=0,z-=-πbTb,--πuTu,--π+T+,--πdTd,--π-T-,-=0.

If we also temporarily let -u1C-u2A=(zbzuz+), then

zb=-πbΛ(D1+D2)-1Σ-1-πb(0.5Σ)(D1-D2)=-πb(D1+D2)-1[ΛΣ-1+0.5Σ(D12-D22)]=-πb(D1+D2)-1[ΛΣ-1-Σ-1Λ]=0.

Note that the matrices D1, D2, ∑, and Λ are all diagonal. We can also show that

zu=-πbTb,u-πuTu,u-π+T+,u-πdTd,u-π-T-,u=0,z+=-πbTb,+-πuTu,+-π+T+,+-πdTd,+-π-T-,+=0.

We left the proof of Mv = 0 to the reader, which is arithmetic.

By examining the iteration kernel of the SDA algorithm, we see the inverse matrices IHkGk and IGkHk, which we can apply the GTH-like algorithm provided that their triplet representation are available. In general, the triplet representation are not available. However, the spectral result in Proposition 1 make it possible, which is shown in the following proposition.

Proposition 2

Recall the matrices Ek,Gk, Hk, Fk in the kth iteration of SDA algorithm. Then, it satisfies that, for all k ≥ 0,

(I-HkGk)v2=Fkv2+HkEkv1and(I-GkHk)v1=Ekv1+GkFkv2,

where v=(v1v2)is the right eigenvector corresponding to 0 of the matrix M of which the form is given in (4.3).

Proof

Since Mv=(-B-C-D-A)v=0 from Proposition 1 and M is an M-matrix, the matrix (γI-B-C-DγI-A) is non-singular M-matrix. Moreover, since

(γI+BCDγI+A)v=γv+(BCDA)v=γv-(BCDA)v=(γI-B-C-DγI-A)v,

it holds that

(E0G0H0F0)v=(γI-B-C-DγI-A)-1(γI+BCDγI+A)v=v.

We can also show that

(EkGkHkFk)v=v

for all k > 0 by induction using the same way in Theorem 3.2 of Xue and Li (2017), which yields that for all k ≥ 0

(I-Gk-HkI)v=(Ek00Fk)v.

Multiplying (I Gk) to the right of both sides, we can obtain the result of (IGkHk). Similarly, we can get the result of (IHkGk) by using the equation

(FkHkGkEk)(v2v1)=(v2v1).

Since IHkGk and IGkHk are non-singular M-matrices for each k ≥ 0, Proposition 2 makes it possible to compute their inverses using the GTH-like algorithm.

5. Numerical study

For the numerical study, we used MATLAB (R2019a) running on a Windows 10 64 bit in DELL Server with dual processor Intel (R) Xeon (R) Gold 6126 CPU @ 2.60GHz and 256GB of main memory.

To see the effect of GTH-like algorithm, we compare the original SDA algorithm(A1) with the subtraction-free SDA(A2) which uses the GTH-like algorithm. For the study, we take one of the examples in Example 3.1, of which S(B) = {1}, S(M) = {1, …, n}, Γ01=-(15/8)In,Γ11=-(15/8n)Jn,βk1-=βk1-=(1/2n)1n, for all kS(M), Bk1-=Bk1+=-2In for all kS(M), μb=(3/2)1n, and σb=7/41n. Then, the infinitesimal generator of the embedded process of the superposition is represented as

T=(Tb,bTb,+Tb,-T+,bT+,+T+,-T-,bT-,+T-,-)=(-158In1516nJn1516Jn2nJn-2In0n×n2nJn0n×n-2In),

where Jn denotes the n-dimensional square matrix of 1’s.

We note that the average drift value of (F, J), average value of aggregated amount of fluid per unit time, is (3/2) for all n. To investigate the effects of the average drift value, we consider another parameter m = (2/3) × μ̄m, which makes the average drift value of (dm × F, J) be μ̄m. For the MMBM (dm × F, J) with τm and Hm being its first passage time to 0 and the H-matrix, it holds that for all m and x

P[τ<J(0)=1,F(0)=x]=P[τm<J(0)=1,dmF(0)=dmx]=e1e(dmx)Hm1=34e-x+14e-3x,

where e1 is a unit vector with appropriate dimension, in which its first element is 1 and the others are all 0. For the exact function of x, we refer to Example 6.1 of Asmussen (1995).

We choose the values of dm so that the average drift value μ̄m of the MMBM’s (dm × F, J) is to be 0.00001, 0.0001, 5, and then compute P[τm < ∞|J(0) = 1, dmF(0) = 3dm] to compare with the exact value in (5.1). Note that the diffusion and drift vector of (dm × F, J) are given as (dmσb00) and (dmμb1-1). To check our computations even more, we consider the MMBM (F, J) for which S = {1, …, 3n} with n = 10, 100, 500. For the case of μ̄ = 0, we used the MMBM (F0, J), of which all the parameters are same as (F, J) except μb being given as 0n×1. We compute P[τ0 < ∞|J(0) = 1, F0(0) = 3] = e1e3H01, which is equal to 1. Here, τ0 and H0 denotes the first passage time to 0 and the H-matrix of (F0, J). The numerical results are presented in Tables 1. From the table, we can observe that the GTH-like algorithm makes the doubling algorithm more stable and accurate, especially when the G-MMBM is critical, that is, the average value of drift is 0.

6. Concluding remarks

In this paper, we introduce how to parameterize the superposition of the MMBM and MAP with bilateral ph-type distributed jumps and also the non-symmetric algebraic Riccati equation to be used to compute the first passage probability of the superposed process. We note that the coefficient matrices of the equation can be obtained from the parametrization. Furthermore, we provide the spectral results which makes it possible to apply the GTH-like algorithm to a structure-preserving doubling algorithm for solving the equation. Through a numerical study, we have shown that using the GTH-like algorithm enhances the performance of the structure-preserving algorithm and makes it remarkably more stable and accurate especially when the superposed process is in the critical case.

Figures
Fig. 1. Illustration of how a sample path of (, ) is transformed to a sample path of its embedded process (, ). Note that , 늿 , 늿 , 늿 , , 늿 , and 늿 .
Fig. 2. (a) Level process and phases of . (b) Completed graph and phases of .
TABLES

Table 1

Comparison of error values

(μ̄, n)(0, 10)(0,100)(0, 500)(10−5, 10)(10−5,100)(10−5, 500)
A1NaN4.98E–08NaN2.84E-161.70E–155.07E–16
A21.55E-151.71E–141.78E–141.94E–161.43E–159.71E–17
(μ̄, n)(10−4, 10)(10−4,100)(10−4, 500)(5, 10)(5,100)(5, 500)
A12.84E–161.69E-155.27E–162.08E–161.71E–155.69E–16
A22.08E–161.51E-154.84E–171.32E–161.54E–156.94E–17

References
  1. Ahn S and Ramaswami V (2005). Bilateral phase type distributions. Stochastic Models, 21, 239-259.
    CrossRef
  2. Ahn S (2016). Total shift during the first passages of Markov modulated Brownian motion with bilateral ph-type jumps: Formulas driven by the minimal solution matrix of a Riccati equation. Stochastic Models, 32, 433-459.
    CrossRef
  3. Ahn S (2017). Time-dependent and stationary analyses of the two-sided reflected Markov modulated Brownian motion with bilateral ph-type jumps. Journal of the Korean Statistical Society, 46, 45-69.
    CrossRef
  4. Alfa AS, Xue J, and Ye Q (2001). Accurate computation of the smallest eigenvalue of a diagonally dominant M-matrix. Math Comp, 71, 217-236.
    CrossRef
  5. Asmussen S (1995). Stationary distributions for fluid flow models with or without Brownian noise. Stochastic Models, 11, 1-20.
  6. Asmussen S (2003). Applied Probability and Queues (2nd ed), Springer Verlag.
  7. Belhaj M (2010). Optimal dividend payments when risk reserves follow a jump diffusion process. Mathematical Finance, 20, 313-325.
    CrossRef
  8. Bini DA, Iannazzo B, and Meini B (2012). Numerical Solution of Algebraic Riccati Equations, Philadelphia, Siam.
  9. Guo XX, Lin WW, and Xu SF (2006). A structure-preserving doubling algorithm for nonsymmetric algebraic Riccati equation. Numer Math, 103, 393-412.
    CrossRef
  10. Jiang Z and Pistorius MR (2008). On perpetual American put valuation and first passage in a regime-switching model with jumps. Finance Stoch, 12, 331-355.
    CrossRef
  11. Liu C and Xue J (2012). Complex nonsymmetric algebraic Riccati equations arising in Markov modulated fluid flows. SIAM J Matrix Anal Appl, 33, 569-596.
    CrossRef
  12. Wang W, Wang W, and Li R (2012). Alternating-directional doubling algorithm for M-matrix algebraic Riccati equations. SIAM J Matrix Anal & Appl, 33, 170-194.
    CrossRef
  13. Whitt W (1980). Some useful functions for functional limit theorems. Mathematics of Operations Research, 5, 67-85.
    CrossRef
  14. Xue J and Li R (2017). Highly accurate doubling algorithm for M-matrix algebraic Riccati equations. Numer Math, 135, 733-767.
    CrossRef