# Computing Hitting Probabilities of Markov Chains: Structural Results with regard to the Solution Space of the Corresponding System of Equations

Authors: Hendrik Baumann; Thomas Hanschke
Journal: Journal of Applied Mathematics
 (2020)

---

## Abstract

In a previous paper, we have shown that forward use of the steady-state difference equations arising from homogeneous discrete-state space Markov chains may be subject to inherent numerical instability. More precisely, we have proven that, under some appropriate assumptions on the transition probability matrixP, the solution space S of the difference equation may be partitioned into two subspaces S=S1⊕S2, where the stationary measure of P is an element of S1, and all solutions in S1 are asymptotically dominated by the solutions corresponding to S2. In this paper, we discuss the analogous problem of computing hitting probabilities of Markov chains, which is affected by the same numerical phenomenon. In addition, we have to fulfill a somewhat complicated side condition which essentially differs from those conditions one is usually confronted with when solving initial and boundary value problems. To extract the desired solution, an efficient and numerically stable generalized-continued-fraction-based algorithm is developed.


---

## Body

## 1. Introduction

This paper is dedicated to homogeneous discrete-time Markov chainsX=Xkk≥0 with a countably infinite set of states (labelled as 0,1,2,…) and a one-step transition matrix P=Pij0≤i,j<∞ of the form(1)P=1000……P10P11P120…P20P21P22P230…⋮⋮⋱⋱⋱⋱Pn0Pn1…Pn,n−1PnnPn,n+10…0Pn+1,1Pn+1,2…Pn+1,nPn+1,n+1Pn+1,n+20…⋱⋱⋱⋱⋱⋱,with some n≥1 and Pk,k+1≠0 for all k∈ℕ=1,2,3,… and Pk,k−n≠0 for all k∈n,n+1,n+2,…. P0ℓ=δ0ℓ implies that 0 is an absorbing state, and the assumptions on Pk,k+1 and Pk,k−n guarantee that state 0 is accessible from any state k>0; that is, for all k>0, there is some m∈ℕ with pk0m>0, where pijm=Pm.We shall calculate the sequenceαkk=0∞, where αk denotes the probability that starting from k>0 the Markov chain ultimately enters the absorbing state 0 (note that we only have assumed αk>0 for all k∈ℕ. Although there is only one absorbing state, due to the infinite number of states, there is no guarantee that αk=1.); e.g., in a population model, this corresponds to eventual extinction of the population. Note that the choice P0j=δ0j is only due to simplicity. For arbitrary P0j, state 0 is not absorbing anymore, but these values do not have any impact on the probabilities αk. Since they still reflect the probabilities that the state 0 is eventually reached, they are referred to as hitting probabilities.The sequenceαkk=0∞ is the uniquely determined minimal nonnegative solution to(2)xk=∑ℓ=0∞Pkℓxℓ,k=1,2,3,…,x0=1,for example, see Section 8.4 in [1]. Due to the structure of P given in (1), we conclude that αkk=0∞ is the minimal nonnegative solution to the n+1th-order homogeneous linear recurrence relation(3)Pk+n,k+n+1xk+n+1+Pk+n,k+n−1xk+n+⋯+Pk+n,kxk=0,k=0,1,2,…,subject to the side conditions(4)x0=1,x1=P10x0+P11x1+P12x2,x2=P20x0+P21x1+P22x2+P23x3,⋮xn−1=Pn−1,0x0+⋯+Pn−1,n−1xn−1+Pn−1,nxn.Since we assumePk+n+1,k+n≠0 for all k, system (4) has rank n, such that the desired solution αkk=0∞ is uniquely determined by α1 ((4) and (3)).The solution of the difference scheme involves two problems:(i)
To start the sequenceαkk=0∞, we need α1. Put αk=xk1α1+xk2 for k=0,1,2,…, where x01=0, x11=1, x02=1, and x12=0 (we will refer to this representation of αkk=0∞ as (FW1)). Due to the linear structure of (3) and (4), xk1k=0∞ and xk2k=0∞ satisfy (3) and (4) up to the first two equations. Another representation of αkk=0∞ is αk=x˜k1c+x˜k21−c with xk1=1 for k=0,1,2,…, and another solution xk2k=0∞ to (3) and (4) being linearly independent of xk1k=0∞. This representation (which we will refer to as (FW2)) results from the fact that P is stochastic.If the difference equation can be solved analytically,α1 and c, respectively, can be determined by the requirement that αkk=0∞ becomes nonnegative and minimal. To prove this assertion, knowledge of all solutions to (3) and (4) (including their asymptotic behaviour) is necessary. But if the determination of αkk=0∞ is based on numerical calculations such that only a finite number of components of αkk=0∞ can be checked, these proofs fail. In this situation, the determination of α1 and c, respectively, is an open problem.(ii)
From the theory of linear difference equations (see Miller [2]), it is known that there exists n+1 linearly independent functions x0,x1,…,xn defined on ℕ0 which assume prescribed values at k=0,1,…,n and satisfy the homogeneous recurrence relation (3) for all k≥0. In the sequel, it is shown that the solution space S of the recurrence relation (3) is the direct sum(5)S=S1⊕S2,of two subspacesS1 and S2 with the property that every solution y∈S2 dominates over every solution x∈S1, i.e.,(6)limk⟶∞xkyk=0.In addition, we havedimS1=n,dimS2=1, and α=αkk=0∞∈S1. From the literature [3–7], it is known that computing solutions x∈S1 by means of forward computation is not a meaningful procedure since forward computation of a subdominant solution becomes numerically unstable. The reason is as follows: if we replace xk0,…,xk0+n−1 by approximate values zk0,…,zk0+n−1 (due to rounding errors, for example), we would have zk=αxk+βyk for all k≥k0 with some β>0 and y∈S2. Even with infinite precision, this would imply zk−xk/xk⟶∞ as k⟶∞; that is, the relative error becomes arbitrarily large (see [4] for more details in the case n=2). Since rounding errors are inevitable, forward computation of x∈S1 will always fail.In [8, 9], it was shown that the problem of numerical instability also arises in the context of computing invariant measures of time-homogeneous and discrete-state space Markov chains with band-structured transition probability matrices. The idea was to construct a decoupled recursion for the solutions x∈S1 in which the dominating solutions y∈S2 do not appear. To carry out the reduction of order for the linear difference equation (3), we apply a technique which relies on generalized-continued fractions. Since the resulting linear difference equation is of order n, instead of n+1, the problem of finding the initial value α1 (or the value c) becomes obsolete.


## 2. Generalized-Continued Fractions and Linear Difference Equations

For the following, we need some results on generalized-continued fractions (GCFs). To point out the analogy to our results in [8], we use the same notations and terminologies.Definition 1.
A generalized-continued fraction (GCF) of dimensionn is an n+1-tuple a1,…,an,b of real-valued sequences and a convergence structure as follows. Let A1,…,An,B be the principal solution of the corresponding homogeneous linear difference equation:(7)xk+n+1=bkxk+n+aknxk+n−1+⋯+ak1xk,k=0,1,2,…,satisfying(8)Aji=δi,j+1=1fori=j+10fori≠j+1,i=1,…,n;j=0,…,n,Bj=0,forj=0,…,n−1,Bn=1.
The GCF is said to converge iff all the limits(9)ξi=limk⟶∞AkiBk,i=1,…,n,do exist.
Convergence of a GCF is indicated by the notation(10)ξ1⋮ξn=ka01a11a21……………a0na1na2n…b0b1b2….
Terminating a GCF after theN-th column, we get the so-called N-th approximants:(11)AN+n+11BN+n+1⋮AN+n+1nBN+n+1=a01a11…aN1…………a0na1n…aNnb0b1…bN.
The term “generalized-continued fraction” becomes obvious if forward computation of theN-th approximant of a GCF is replaced by the equivalent backward algorithm. It is known that the recurrence relation (7) can be converted into a first-order vector recursion of the form(12)uk+1=Wkuk,k=0,1,2,…,by putting(13)uk=xkxk+1⋮xk+n−1xk+n,Wk=010…00001…00………………00……01ak1ak2……aknbk,k=0,1,2,….
Comparing (7) and (12), it is seen that the numerators AN+n+11,…,AN+n+1n and the denominators BN+n+1 of the GCF appear in the last row of the matrix:(14)WNWN−1⋯W0.
Consider now the backward recurrence scheme(15)Nvk=Nvk+1Wk,k=N,N−1,…,0,with initial vector NvN+1=0,…,0,1. Then,(16)Nv0=NvN+1WNWN−1⋯W0=AN+n+11,…,AN+n+1n,BN+n+1.
Writing (15) in the expanded form and inserting the structure of Wk, we get(17)Nvk1=Nvk+1n+1ak1,k=N,N−1,…,0,Nvki=Nvk+1i−1+Nvk+1n+1aki,i=2,…,n,k=N,N−1,…,0,Nvkn+1=Nvk+1n+Nvk+1n+1bk,k=N,N−1,…,0,or equivalently(18)Nrk1≔Nvk1Nvkn+1=ak1bk+Nrk+1n,k=N,N−1,…,0,Nrki≔NvkiNvkn+1=aki+Nrk+1i−1bk+Nrk+1n,i=2,…,n;k=N,N−1,…,0,with initial values NrN+1i=0 for i=1,…,n. Hence, (19)AN+n+1iBN+n+1=Nr0i,i=1,…,n,ξi=limN⟶∞Nr0i,i=1,…,n.
Alternative procedures for computing GCFs are described in [10].
The relations between GCFs and linear difference equations have been first recognized by Perron [11]. Perron’s results were generalized by Van der Cruyssen [10], Hanschke [12, 13], and Levrie and Bultheel [14]. The following theorem is due to Van der Cruyssen [10].Theorem 1.
A GCFa1,…,an,b converges iff there are n+1 linearly independent solutions x1,…,xn,y of the recurrence relation (7) satisfying(20)limk⟶∞xkiyk=0,i=1,…,n,x01…xn−11⋮x0n…xn−1n≠0.
As pointed out by Gautschi [4, 5] and Van der Cruyssen [10], forward computation of a solution x∈spanx1,…,xn is numerically unstable. Van der Cruyssen [10] establishes that if the limits(21)ξl1⋮ξln=al1al+11al+21……………alnal+1nal+2n…blbl+1bl+2…,exist for all l≥0, then x∈spanx1,…,xn iff(22)xl+n=−ξlnxl+n−1−⋯−ξl1xl,l=0,1,2,….
Combining (18), (19), and (22), one obtains an efficient algorithm (which we will refer to as Van der Cruyssen’s algorithm) for approximating the first L+1 components x0,…,xL of an element x∈spanx1,…,xn with prescribed values x0,…,xn−1:Step 1.
SelectN>L−n and define Nrli for i=1,…,n;l=0,1,…,N+1 by(23)NrN+11=⋯=NrN+1n=0,Nrl1=al1bl+Nrl+1n,l=N,N−1,…,0,Nrli=ali+Nrl+1i−1bl+Nrl+1n,i=2,…,n;l=N,N−1,…,0.Step 2.
Set(24)x0N=x0,…,xn−1N=xn−1,xl+nN=−Nrlnxl+n−1N−⋯−Nrl1xlN,for l=0,1,…,L−n.Step 3.
IncreaseN until the accuracy of the iterates is within prescribed limits.
For anyl, the vector Nrl1,…,NrlnT is an approximant of a GCF. Hence, convergence of the algorithm is related to convergence of GCFs. For the latter, we cite a result which can be interpreted as a generalization of Pringsheim’s convergence criterion ([15], p. 58) for ordinary continued fractions.Theorem 2 (Levrie [16] and Perron [17]).
The GCF(25)a01a11a21……………a0na1na2n…b0b1b2…,converges if it satisfies the so-called dominance condition, i.e.,(26)ak1+ak2+⋯+akn+1≤bk,for allk≥0.


## 3. Main Results

To make (3) congruent with (7), we put(27)aki=−Pk+n,k+i−1Pk+n,k+n+1,i=1,…,n,k=0,1,2,…,(28)bk=−Pk+n,k+n−1Pk+n,k+n+1,k=0,1,2,….Theorem 3.
Letaki,bk be defined as in (27) and (28), and let Pk,k+1≠0 for all k∈ℕ and Pk,k−n≠0 for all k≥n. Then, the limits(29)ξl1⋮ξln=al1al+11al+21……………alnal+1nal+2n…blbl+1bl+2…,exist for all l≥0. The vector αkk=0∞ of the absorption probabilities with respect to state 0 is a dominated solution of (7) and satisfies the decoupled recursion(30)xl+n=−ξlnxl+n−1−⋯−ξl1xl,l=0,1,2,….Proof.
As in the corresponding statement Theorem 3 in [8], the basic idea for the proof is that αkk=0∞ is the limit of the solutions of appropriately truncated systems of equations. We give some details on the probabilistic interpretation of these systems although these considerations are well-known from standard literature on Markov chains (e.g., [18]).
ByNP=Pik0≤i,k≤N, denote the northwest corner truncation of P, and by XN=xkNk≥0, denote a discrete-time Markov chain with state space 0,1,…,N+1 and block-transition probability matrix(31)NPq01,where q=qii=0N is a column vector with qi=∑k=N+1∞Pik. For XN, both 0 and N+1 are absorbing states. By αkN, denote the probability that XN will be absorbed in state 0 when starting in state k. Obviously, we have αN+1N=0, and hence, the values αkN define a solution to the truncated system:(32)Nx0=1,Nxk=∑ℓ=0NPkℓ⋅Nxℓ,k=1,…,N.
Since for the original Markov chainX, state 0 was assumed to be accessible from all k>0, either state 0 or N+1 is accessible from all states k∈1,…,N for XN. Hence, every state ∈1,…,N is transient, implying that NI−Pkℓk,ℓ=1N is invertible, where IN is the N×N-identity matrix. Thus, αkNk=0N is the unique solution to (32).
Under the given assumptions on the structure ofP, (32) coincides with the homogeneous linear difference equation (3) (or (7) in terms of bk and aki), coupled to the original side conditions (4) (assuming that N>n) and the additional side condition:(33)NxN+1=0.
For the GCFξl1⋮ξln defined in (29), let us denote the numerators by Al,k1,…,Al,kn and the denominators by Bl,k. Then, Al,kik and Bl,kk satisfy(34)xk+n+1=bk+lxk+n+ak+lnxk+n−1+⋯+ak+l1xk,k=0,1,2,…,and build up a fundamental system of solutions to equation (7) for k≥l. Hence, any solution of the recurrence relation (7) can be expressed in terms of these functions. In view of their initial values(35)Al,ji=δi,j+1=1,i=j+10,i≠j+1,i=1,…,n;j=0,…,n,Bl,j=0,forj=0,…,n−1,Bl,l+n=1,we get(36)Nxk=NxlAl,k−l1+⋯+Nxl+n−1Al,k−ln+Nxl+nBl,k−l,k=l,l+1,….
Choosingk=N+1 and utilizing (33), equation (36) reduces to(37)Nxl+nBl,N−l+1=−NxlAl,N−l+11−⋯−Nxl+n−1Al,N−l+1n.
Dividing both sides of (37) by Bl,N−l+1 and Nx0=1, we formally arrive at(38)Nxl+n=−NxlAl,N−l+11Bl,N−l+1−⋯−Nxl+n−1Al,N−l+1nBl,N−l+1.
Passing to the limitN⟶∞, we get the decoupled recursion (30), provided that the limits limN⟶∞Al,N−l+1i/Bl,N−l+1 do exist for i=0,1,…,n, and that limN⟶∞xkN=αk for all k∈ℕ0.(i)
Using the definitions ofaki and bk and the fact that P is stochastic, we immediately find(39)ak1+ak2+⋯+akn+1=∑i=1nPk+n,k+i−1+Pk+n,k+n+1Pk+n,k+n+1=1−Pk+n,k+nPk+n,k+n+1=bk.(ii)
Hence, Theorem2 guarantees that the limits limN⟶∞Al,N−l+1i/Bl,N−l+1 do exist.(iii)
As pointed out above, the unique solution to (3) and (4) with the additional side condition (33) is given by Nxk=αkN, where for the truncated Markov chain XN, αkN is the probability to enter state 0 before N+1 when starting in k. For the original Markov chain X, let τ0 and τN+1 be the first hitting times of 0 and N+1, respectively, that is(40)τi=infm∈ℕ:Xm=i.(iv)
Then, we can writeαkN=ℙτ0<τN+1X0=k. For N⟶∞, these probabilities increase monotonically up to the limit ℙτ0<∞X0=k=αk. Hence,(41)limN⟶∞Nxk=limN⟶∞αkN=αk.Remark 1.
Notice that the statement of Theorem3 consists of three parts:(i)
The sequence of absorption probabilities is dominated by other solutions of the difference equation (3).(ii)
By means of GCFs, we are able to construct a decoupled recursion which is not affected by numerical instability.(iii)
Since the decoupled recursion is a linear difference equation of ordern (instead of order n+1), the problem of finding the initial value α1 (or the value c) is eliminated since α1,…,αn can be obtained by exploiting (30) for l=0 and the side conditions (4). In other words, the calculation of α1 or c is incorporated in Van der Cruyssen’s algorithm.
The fact that the solution of the truncated system converge to the desired minimal solution is far from new. But to the best of our knowledge, although the numerical solution of the system is associated with massive problems, both the derivation of the structural resultS1⊕S2 and the relationship to generalized-continued fractions have not been addressed in the literature before (up to the corresponding considerations in [8] concerning the computation of invariant measures of Markov chains).


## 4. Continuous-Time Markov Chains

The results of the preceding section can easily be extended to continuous-time Markov chainsY=Ytt≥0 generated by a conservative, irreducible, and regular Q-matrix of the form(42)Q=0000……Q10Q11Q120…Q20Q21Q22Q230…⋮⋮⋱⋱⋱⋱Qn0Qn1…Qn,n−1QnnQn,n+10…0Qn+1,1Qn+1,2…Qn+1,nQn+1,n+1Qn+1,n+20…⋱⋱⋱⋱⋱⋱,where n≥1, Qkk≠0, and Qk,k+1≠0 for all k∈1,2,3,… and Qk,k−n≠0 for all k∈n,n+1,n+2,….For the continuous-time Markov chainY, the probability αk of eventually reaching state 0 when starting in state k is equal to the probability of eventually reaching state 0 for the embedded jump chain X=Xmm=0∞ with transition probability matrix P=Pkℓk,ℓ=0∞ where(43)Pkℓ=δkℓ,Qkℓ=0⟺k=0,Qkℓ−Qkk+δkℓ,Qkℓ≠0.Hence,αk solves the homogeneous system (3) with the side conditions (4), where Pkℓ is defined as above; that is, we obtain the modified homogeneous linear difference equation:(44)Qk+n,k+n+1xk+n+1+Qk+n,k+nxk+n+⋯+Qk+n,kxk=0,k=0,1,2,…,subject to the side conditions(45)x0=1,0=Q10x0+Q11x1+Q12x2,0=Q20x0+Q21x1+Q22x2+Q23x3,⋮0=Qn−1,0x0+⋯+Qn−1,n−1xn−1+Qn−1,nxn.In view of the difference equation (3) on which the definition of GCFs relies, we set(46)aki=−Qk+n,k+i−1Qk+n,k+n+1,i=1,…,n,bk=−Qk+n,k+nQk+n,k+n+1.Using (43), we obtain that this corresponds exactly to the choice of aki and bk given by (27) and (28), respectively. Therefore, Theorem 3 can be directly transferred to continuous-time Markov chains.Theorem 4.
Letaki,bk be defined as in (46) and let Qk,k+1≠0 for all k∈ℕ and Qk,k−n≠0 for all k≥n. Then, the limits(47)ξl1⋮ξln=al1al+11al+21……………alnal+1nal+2n…blbl+1bl+2…,exist for all l≥0. The vector αkk=0∞ of absorption probabilities with respect to state 0 is a dominated solution of (7) and satisfies the decoupled recursion:(48)xl+n=−ξlnxl+n−1−⋯−ξl1xl,l=0,1,2,….


## 5. Numerical Examples

As a simple example, we consider a modification of the famous gambler’s ruin problem: A gambler gambles against a bank (with infinite capital) and wins or loses a coin in each round with probabilityp and 1−p, respectively. Let Xk be the number of coins the gambler owns after n rounds. Then, Xkk≥0 is a discrete-time Markov chain with the above transition structure where n=1, Pi,i−1=1−p, Pii=0, and Pi,i+1=p for all i≥1. Due to n=1, (4) simplifies to α0=1, and (3) reduces to(49)pxk+2−xk+1+1−pxk=0,k=0,1,2,….We briefly demonstrate how to determine the exact solution: As pointed out above (representation (FW2)), a solution is given by the constantx˜k1=1 for k=0,1,2,.... For p=1/2, another linearly independent solution is given by x˜k2=k, and hence, αk=c+1−ck for k=0,1,2,…. Due to nonnegativity and minimality, we have c=1, and hence, αk=1 for all k=0,1,2,…. For p<1/2, another solution is given by x˜k2=1−p/pk, and hence, αk=c+1−c1−p/pk for k=0,1,2,…. Due to 1−p/p>1, minimality and nonnegativity imply c=1, and hence, still αk=1 for k=0,1,2,…. For p>1/2, we still have αk=c+1−c1−p/pk, but here, due to 1−p/p<1, we have to choose c=0 and obtain αk=1−p/pk for k=0,1,2,…. Regarding the analytic derivation of the exact solution, we could use the representation (FW1) in a similar way.In most applications, it is not possible to find the explicit representation ofαk, and therefore, we have to use numerical calculations. For the sake of simplicity and in order to be able to compare our results with the exact solution, we apply numerical computation first to the gamblers’s-ruin example. Using data type double in C++, for p=1/3, for the values xk1k=0∞ and xk2k=0∞ in the representation (FW1), we obtain the results depicted in Table 1.Table 1
Forward computation applied to the gambler’s ruin problem.
k
x
k
1
x
k
2
α
k
=
x
k
1
+
x
k
2
0
0
1
1
1
1
0
1
2
3
−
2
1
3
7
−
6
1
4
15
−
14
1
10
1023
−
1022
1
20
1048575
−
1048574
1
30
1073741823
−
1073741822
1
40
1099511627774.999878
−
1099511627774.000244
0.999634
50
1125899906842622.875
−
1125899906842622.250
0.625
75
·
·
−
25165824
100
·
·
−
844424930131968Note that the knowledge thatαk=xk1+xk2 relies on α1=1. So even if we knew the exact value for α1 (which we usually do not know if we use numerical computations), we can get extremely wrong results from applying forward computation of the representation (FW1). In the representation (FW2), we can use x˜k2=xk1 or x˜k2=xk2. For both choices, the numerically computed values suggest that x˜k2 increases exponentially, and hence, for obtaining the minimal nonnegative solution, we should choose c=1 and c=0. So here, using numerically computed values for FW2 works. However, it is obvious that this argument is susceptible to misinterpretations.The GCF-based approach requires a truncation levelN. For N=120, we obtain αk=1 for k≤100 up to errors <10−6.We demonstrate the use of the GCF-based approach in more detail by applying it in a situation where we do not know the exact solution. We will see that the method provides plausible results. Letn=2 and(50)Pk,k+1=12+gk,k=1,2,…,Pkk=14,k=1,2,…,P10=14−gk,Pk,k−1=18−gk,k=2,3,…,Pk,k−2=18,k=2,3,….Forgk=0, gk=1/2k+2, and gk=1/8k, some values can be found in Table 2. For N=120 and N=150, the first six digits of the values α0,…,α100 coincide which shows that our results are trustworthy. Furthermore, as pointed out above, we can check the results with respect to plausibility: For gk=0, we can compute the exact solution by means of the theory of linear difference equations with constant coefficients and obtain(51)αk=2+5251+54k−2−5251−54k,k=0,1,2,…,see first column of Table 2. In general, the term gk may be interpreted as a kind of disturbance. For gk>0, the effect of the Markov chain “drifting away” from state 0 is enforced. In particular, since we have limk⟶∞gk=0 for our choices of g, this effect is strong for small states k. Hence, for small k, αk decreases quite slowly for gk=0 and more quickly for the other choices. Since 1/8k>1/2k+2, we observe that αk decreases faster for gk=1/8k.Table 2
Numerical computation of the absorption probabilities by means of the GCF-based approach.
k
α
k for gk=0 (exact)
α
k for gk=0
α
k for gk=1/2k+2
α
k for gk=1/8k
0
1.000000
1.000000
1.000000
1.000000
1
0.750000
0.750000
0.579198
0.545488
2
0.625000
0.625000
0.495038
0.454585
3
0.500000
0.500000
0.373473
0.323282
4
0.406250
0.406250
0.303615
0.251803
5
0.328125
0.328125
0.242391
0.191475
6
0.265625
0.265625
0.195995
0.148602
8
0.173828
0.173828
0.127894
0.089963
10
0.113770
0.113770
0.083650
0.055233
15
0.039429
0.039429
0.028984
0.016904
20
0.013665
0.013665
0.010045
0.005333
50
0.000024
0.000024
0.000017
0.000007


## 6. Conclusion and Further Research

This paper deals with the evaluation of higher-order linear difference equations associated with the computation of hitting probabilities of discrete-state space Markov chains. To extinguish effects of inherent instability, a generalized-continued-fraction algorithm is established. Although aspects of inherent instability have been widely discussed in the literature (e.g., [3–7]), their relevance to Markov chains has not been considered yet. It is to be expected that much more problems in stochastic modeling are subject to this phenomenon.Similar to the computation of absorption or hitting probabilities, mean hitting times can be calculated: letαi=1 and X0=i. Then, one may inquire for the mean time mi which elapses until state 0 is assumed. With m0=0, we have(52)mi=1+∑j=0∞pijmj,or0=1+∑j=0∞qijmj.Obviously, system (52) for mii=0∞ is an inhomogeneous version of the homogeneous system for αii=0∞. If the transition structure is given by (1) or (42), respectively, we obtain an inhomogeneous version of the difference equation (3), which is subject to numerical instability. Hence, the same is true for the inhomogeneous system so that forward computation is not a meaningful procedure. To overcome this problem, further research is required.Further research is also required for extending our results to block band transition matrices because there is no counterpart to Pringsheim’s convergence criterion (Theorem2) in the matrix-driven case.



---

*Source: 9874072-2020-01-13.xml*