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 matrix P, the solution space S of the difference equation may be partitioned into two subspaces , where the stationary measure of P is an element of , and all solutions in are asymptotically dominated by the solutions corresponding to . 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.
1. Introduction
This paper is dedicated to homogeneous discrete-time Markov chains with a countably infinite set of states (labelled as ) and a one-step transition matrix of the formwith some and for all and for all . implies that 0 is an absorbing state, and the assumptions on and guarantee that state 0 is accessible from any state ; that is, for all , there is some with , where .
We shall calculate the sequence , where denotes the probability that starting from the Markov chain ultimately enters the absorbing state 0 (note that we only have assumed for all . Although there is only one absorbing state, due to the infinite number of states, there is no guarantee that .); e.g., in a population model, this corresponds to eventual extinction of the population. Note that the choice is only due to simplicity. For arbitrary , state 0 is not absorbing anymore, but these values do not have any impact on the probabilities . Since they still reflect the probabilities that the state 0 is eventually reached, they are referred to as hitting probabilities.
The sequence is the uniquely determined minimal nonnegative solution tofor example, see Section 8.4 in [1]. Due to the structure of P given in (1), we conclude that is the minimal nonnegative solution to the th-order homogeneous linear recurrence relationsubject to the side conditions
Since we assume for all k, system (4) has rank n, such that the desired solution is uniquely determined by ((4) and (3)).
The solution of the difference scheme involves two problems:(i)To start the sequence , we need . Put for , where , , , and (we will refer to this representation of as (FW1)). Due to the linear structure of (3) and (4), and satisfy (3) and (4) up to the first two equations. Another representation of is with for , and another solution to (3) and (4) being linearly independent of . 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, and c, respectively, can be determined by the requirement that 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 is based on numerical calculations such that only a finite number of components of can be checked, these proofs fail. In this situation, the determination of and c, respectively, is an open problem.(ii)From the theory of linear difference equations (see Miller [2]), it is known that there exists linearly independent functions defined on which assume prescribed values at and satisfy the homogeneous recurrence relation (3) for all . In the sequel, it is shown that the solution space S of the recurrence relation (3) is the direct sum of two subspaces and with the property that every solution dominates over every solution , i.e.,
In addition, we have , and . From the literature [3–7], it is known that computing solutions 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 by approximate values (due to rounding errors, for example), we would have for all with some and . Even with infinite precision, this would imply as ; that is, the relative error becomes arbitrarily large (see [4] for more details in the case ). Since rounding errors are inevitable, forward computation of 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 in which the dominating solutions 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 , the problem of finding the initial value (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 dimension n is an -tuple of real-valued sequences and a convergence structure as follows. Let be the principal solution of the corresponding homogeneous linear difference equation:satisfyingThe GCF is said to converge iff all the limitsdo exist.
Convergence of a GCF is indicated by the notationTerminating a GCF after the N-th column, we get the so-called N-th approximants:The term “generalized-continued fraction” becomes obvious if forward computation of the N-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 formby puttingComparing (7) and (12), it is seen that the numerators and the denominators of the GCF appear in the last row of the matrix:Consider now the backward recurrence schemewith initial vector . Then,Writing (15) in the expanded form and inserting the structure of , we getor equivalentlywith initial values for . Hence, 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 GCF converges iff there are linearly independent solutions of the recurrence relation (7) satisfyingAs pointed out by Gautschi [4, 5] and Van der Cruyssen [10], forward computation of a solution is numerically unstable. Van der Cruyssen [10] establishes that if the limitsexist for all , then iffCombining (18), (19), and (22), one obtains an efficient algorithm (which we will refer to as Van der Cruyssen’s algorithm) for approximating the first components of an element with prescribed values :
Step 1. Select and define for by
Step 2. Setfor
Step 3. Increase N until the accuracy of the iterates is within prescribed limits.
For any l, the vector 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 GCFconverges if it satisfies the so-called dominance condition, i.e.,
3. Main Results
To make (3) congruent with (7), we put
Theorem 3. Let be defined as in (27) and (28), and let for all and for all . Then, the limitsexist for all . The vector of the absorption probabilities with respect to state 0 is a dominated solution of (7) and satisfies the decoupled recursion
Proof. As in the corresponding statement Theorem 3 in [8], the basic idea for the proof is that 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]).
By , denote the northwest corner truncation of P, and by , denote a discrete-time Markov chain with state space and block-transition probability matrixwhere is a column vector with . For , both 0 and are absorbing states. By , denote the probability that will be absorbed in state 0 when starting in state k. Obviously, we have , and hence, the values define a solution to the truncated system:Since for the original Markov chain X, state 0 was assumed to be accessible from all , either state 0 or is accessible from all states for . Hence, every state is transient, implying that is invertible, where is the -identity matrix. Thus, is the unique solution to (32).
Under the given assumptions on the structure of P, (32) coincides with the homogeneous linear difference equation (3) (or (7) in terms of and ), coupled to the original side conditions (4) (assuming that ) and the additional side condition:For the GCF defined in (29), let us denote the numerators by and the denominators by . Then, and satisfyand build up a fundamental system of solutions to equation (7) for . Hence, any solution of the recurrence relation (7) can be expressed in terms of these functions. In view of their initial valueswe getChoosing and utilizing (33), equation (36) reduces toDividing both sides of (37) by and , we formally arrive atPassing to the limit , we get the decoupled recursion (30), provided that the limits do exist for , and that for all .(i)Using the definitions of and and the fact that P is stochastic, we immediately find(ii)Hence, Theorem 2 guarantees that the limits do exist.(iii)As pointed out above, the unique solution to (3) and (4) with the additional side condition (33) is given by , where for the truncated Markov chain , is the probability to enter state 0 before when starting in k. For the original Markov chain X, let and be the first hitting times of 0 and , respectively, that is(iv)Then, we can write . For , these probabilities increase monotonically up to the limit . Hence,
Remark 1. Notice that the statement of Theorem 3 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 order n (instead of order ), the problem of finding the initial value (or the value c) is eliminated since can be obtained by exploiting (30) for and the side conditions (4). In other words, the calculation of 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 result 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 chains generated by a conservative, irreducible, and regular Q-matrix of the formwhere , , and for all and for all .
For the continuous-time Markov chain Y, the probability 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 with transition probability matrix where
Hence, solves the homogeneous system (3) with the side conditions (4), where is defined as above; that is, we obtain the modified homogeneous linear difference equation:subject to the side conditions
In view of the difference equation (3) on which the definition of GCFs relies, we set
Using (43), we obtain that this corresponds exactly to the choice of and given by (27) and (28), respectively. Therefore, Theorem 3 can be directly transferred to continuous-time Markov chains.
Theorem 4. Let be defined as in (46) and let for all and for all . Then, the limitsexist for all . The vector of absorption probabilities with respect to state 0 is a dominated solution of (7) and satisfies the decoupled recursion:
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 probability p and , respectively. Let be the number of coins the gambler owns after n rounds. Then, is a discrete-time Markov chain with the above transition structure where , , , and for all . Due to , (4) simplifies to , and (3) reduces to
We briefly demonstrate how to determine the exact solution: As pointed out above (representation (FW2)), a solution is given by the constant for . For , another linearly independent solution is given by , and hence, for . Due to nonnegativity and minimality, we have , and hence, for all . For , another solution is given by , and hence, for . Due to , minimality and nonnegativity imply , and hence, still for . For , we still have , but here, due to , we have to choose and obtain for . 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 , 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 , for the values and in the representation (FW1), we obtain the results depicted in Table 1.
Note that the knowledge that relies on . So even if we knew the exact value for (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 or . For both choices, the numerically computed values suggest that increases exponentially, and hence, for obtaining the minimal nonnegative solution, we should choose and . 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 level N. For , we obtain for up to errors .
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. Let and
For , , and , some values can be found in Table 2. For and , the first six digits of the values coincide which shows that our results are trustworthy. Furthermore, as pointed out above, we can check the results with respect to plausibility: For , we can compute the exact solution by means of the theory of linear difference equations with constant coefficients and obtainsee first column of Table 2. In general, the term may be interpreted as a kind of disturbance. For , the effect of the Markov chain “drifting away” from state 0 is enforced. In particular, since we have for our choices of , this effect is strong for small states k. Hence, for small k, decreases quite slowly for and more quickly for the other choices. Since , we observe that decreases faster for .
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 and . Then, one may inquire for the mean time which elapses until state 0 is assumed. With , we have
Obviously, system (52) for is an inhomogeneous version of the homogeneous system for . 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 (Theorem 2) in the matrix-driven case.
Data Availability
No additional data were used to support this study.
Additional Points
Highlights. (i) Inherent instability of forward computation of absorption probabilities for certain transition structures is rigorously proven by establishing a structural result. (ii) Generalized-continued fractions are used to derive a decoupled difference equation. (iii) Based on this decoupled equation, a numerically stable algorithm is established.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge support by the Open Access Publishing Fund of Clausthal University of Technology.