Ruin probabilities as functions of the roots of a polynomial

A new formula for the ultimate ruin probability in the Cramér–Lundberg risk process is provided when the claims are assumed to follow a ﬁnite mixture of m Erlang distributions. Using the theory of recurrence sequences, the method proposed here shifts the problem of ﬁnding the ruin probability to the study of an associated characteristic polynomial and its roots. The found formula is given by a ﬁnite sum of terms, one for each root of the polynomial, and allows for yet another approximation of the ruin probability. No constraints are assumed on the multiplicity of the roots and that is illustrated via a couple of numerical examples


Introduction
We here study the classical Cramér-Lundberg risk process {R(t) : t ≥ 0}, see [2] or [14], given by The constant u ≥ 0 is the initial capital, c > 0 is the premium per unit of time which is also assumed to be constant, N(t) is a Poisson process with intensity λ > 0, and X 1 , X 2 , . . . are i.i.d. random variables representing the claims' amounts. It is usually assumed that the process of claims' amounts is independent of the Poisson process. The model (1) is one of the simplest representations of the time evolution of the capital of an insurance company. A historical exposition of it can be found in [11]. It was introduced by F. Lundberg in 1903, see [28], and was later revised in 1926, see [29]. The model is now known in terms of the theory of stochastic processes and that was due to the work of 1930 by H. Cramér, see [10]. The time of ruin τ is defined as the first instant of time when a negative value is taken by the process (1), that is, τ = inf {t > 0 : R(t) < 0}. In the case when this set is empty we define τ = ∞, i.e. the ruin never occurs.
The ultimate ruin probability, denoted by ψ(u), is the probability of the occurrence of the event (R(t) < 0) for some finite t > 0. One of the major problems in the mathematical theory of ruin is to find ψ(u) for a given distribution of the claims, not only for the Cramér-Lundberg model (1) but also for the numerous generalizations that have been proposed and studied since the work of Cramér, e.g., [1,8,9,13,15,17,22,24,26,27,40,41]. Unfortunately, closed-form expressions for ψ(u) are hard to find in general, see [2], and much effort has been directed to approximate, bound, numerically compute or study the limiting behavior of ψ(u) as u → ∞, see, e.g., [3,6,7,12,16,18,19,31,37,46]. In particular, see the recent results of A. Grigutis in [20] for a closed-form expression of ψ(u) for the discrete version of the model (1). On the other hand, the Erlang mixture distribution has been frequently considered as a claims' distribution because of its flexibility, see [25,30,32,[43][44][45].
Our aim here is to give an alternative approach to finding ψ(u) in the Cramér-Lundberg model (1) under the hypothesis that claims take values according to a finite mixture of Erlang distributions. We will see that this assumption makes the connection with the theory of recurrence sequences possible. Theorem 1 stated below is our main result. It is a generalization of our previous findings in [35] and states that ψ(u) admits a representation as a linear combination of exponential functions. For this to be possible, one needs to find the zeroes of a polynomial and solve a system of linear equations. The case when the polynomial has only simple roots is studied in [35], whereas in this work we allow for any multiplicity of the roots. Furthermore, in [36] we have successfully applied the same procedure for a discrete-time risk process with finite negative binomial mixture claims and found an expression for ψ(u) in both cases: simple and general roots.
We first define the Erlang mixture distributions of our interest and then use the known formula (6) below to prove our main formula (21). Two numerical examples are given to empirically show our results. The safety loading for the Cramér-Lundberg model (1), also known as the net profit condition, is given by c > λμ. This critical condition can be written as c = (1 + θ)λμ for some θ > 0, and it is also known as the premium loading factor. It is well known, see [2], that when the net profit condition is not satisfied, the ultimate ruin is certain, which is not desirable. It is also known, see [14], that for null initial capital, i.e. u = 0, we have ψ(0) = λμ/c = 1/(1 + θ) for θ > 0. As it is customary, the survival probability P(X > x) is denoted by F (x), and the letters d.f. refer to the distribution function of a random variable.

Erlang mixture distributions
The density of the Erlang distribution is denoted by erlang(k, β)(x), that is, where k ≥ 1 is an integer and β > 0. We use Erlang(k, β)(x) for its d.f. In the case k = 1, the above reduces to the exp(β) distribution. We denote by π the probability function (f N (1), f N (2), . . .) of a r.v. N with values 1, 2, . . . On the other hand, a r.v. S is distributed according to an Erlang mixture with parameters π and β > 0 when its density is We denote this by writing S ∼ ErM(π , β). The identity (2) defines an infinite mixture of Erlang distributions with the vector π being the discrete mixing distribution. The variable N is referred to as the counting random variable of S and it is known that is a continuous d.f. on [0, ∞) and for each natural number n, one defines with f n (k) defined as the difference between F (k/n) and F ((k − 1)/n), for k ≥ 1, then a theorem by Schassberger, see [38], states that F n (x) → F (x) pointwise as n → ∞. Furthermore, when F (x) has bounded support, there is uniformity in the convergence. A proof of these claims can be found in Tijms [42] and in Lee and Lin [25]. On the other hand, some methods used to calculate ψ(u) require knowing the equilibrium distribution of the random variables used in a given risk model. More specifically, given a nonnegative r.v. X with d.f. F (x) and mean μ < ∞, its equilibrium distribution is defined as We use the letters e.d. for the term equilibrium distribution. It is straightforward to show that f e (x) is a genuine probability or density function. For example, it can be shown, see [37], that the e.d. of ErM(π , β) is ErM(π e , β), where π e is the distribution of N e given by the sequence where F N (k − 1) = P(N > k − 1).

Recurrence sequences applied to ruin probabilities
As in our previous work [35], the starting point in our analysis is the ruin probability formula for the Cramér-Lundberg model (1), where Z βu ∼ Poisson(βu) and with claims having a distribution as defined in (2). The formula (6) can be found in Klugman et al. [23] and a detailed proof is given in [35]. It is clear that the sequence C 0 , C 1 , . . . plays a crucial role in the determination of the ruin probability. This sequence of positive numbers is defined recursively as shown in (8) and clearly depends on N through the distribution of N e . From the recursive relation (8), it can be shown that C 0 > C 1 > · · · and that C n → 0 as n → ∞.
For the method proposed in this work to successfully operate it is essential to assume that N has finite support {1, . . . , m} for some fixed integer m ≥ 1, that is, π = (f N (1), . . . , f N (m), 0, . . . ), where f N (1) + · · · + f N (m) = 1. In particular, we assume that f N (m) > 0. For distributions π of this type we call (2) a finite Erlang mixture distribution. That is denoted by finiteErM(π , β) and its density is given by This will be the distribution assumed for the claims in the Cramér-Lundberg model (1). The distribution of the associated r.v. N e is given by the numbers f N e (1) ≥ f N e (2) ≥ · · · ≥ f N e (m), which are the first m terms of (5). for n ≥ 0. Substituting in (6) and some further simplifications yield the well-known formula, see [2] or [14], of the ruin probability for Exp(β) claims, Our goal is to solve the recurrence relation (8) for C n in the finiteErM(π , β) case, where the first m terms C 0 , . . . , C m−1 are taken as initial data. Once this is done we can substitute the values of C n in (6) to find ψ(u). The following procedure is the same as that used in [35]. Define the decreasing and strictly positive sequence where, clearly, 1 > α 1 ≥ α 2 ≥ · · · ≥ α m−1 ≥ α m > 0 and α 1 + · · · + α m = 1/(1 + θ) < 1. One can then prove, see [35], that the recursive relation (8) can be written as the following m-order recurrence sequence where C 0 , . . . , C m−1 are the initial values. The characteristic polynomial associated with the recurrence sequence (12), see [5], being It is known that, see [5] or [39], the solution of (12) can be written as a function of the zeroes of (13). Using standard basic results from the theory of polynomials, see [35], it can be shown that the polynomial (13) has a unique root z 1 which is positive and is such that 0 < |z| < z 1 < 1 for any other root z. Also, the roots of (13) occur in conjugate pairs since (13) has real coefficients. In the following section we will solve (12) and write ψ(u) in terms of the zeroes of (13).

A new formula for ψ(u)
Suppose the solutions of p(y) = 0 are not necessarily simple. Let z 1 , . . . , z be the roots and assume their multiplicities are n 1 , . . . , n , respectively, where 1 ≤ ≤ m and n 1 + · · · + n = m. Suppose z 1 > 0 is the unique positive root which we know has multiplicity n 1 = 1. Thus, (13) can be written as It can be shown, see [39, p. 55], that the solution of (12) is where the m constants b k,j are chosen so that (15) meets the initial data C 0 , . . . , C m−1 . This yields the system of linear equations and the matrix Z has the form where the vertical lines help to visually separate the column vectors for k = 2, . . . , and j = 1, . . . , n k . Observe that the multiplicity of the roots determines the blocks of column vectors. For a better understanding of the system (16), see Examples 2, 3 and 4 below for some particular cases of Z. The m initial values C 0 > · · · > C m−1 are all real, positive and are computed using (8). On the other hand, the numbers b k,1 , . . . , b k,n k are associated with the root z k , for k = 1, . . . , , and can be real or complex depending on the nature of z k . This is explained in the following statement. (13), the following holds.

Proposition 1. For the characteristic polynomial
1. If two roots z k and z j are complex conjugates and have multiplicity n, then the associated coefficients b k,1 , . . . , b k,n and b j,1 , . . . , b j,n , are also complex conjugates, respectively.
2. If a root z k is real with multiplicity n k , then all the associated coefficients b k,1 , . . . , b k,n k are also real.
Proof. Taking complex conjugate of (16) yields the system Z·b = C, where the righthand side is the same real vector C, here the upper bar does not denote conjugate. The matrix Z is the same as Z but now every complex root z k is interchanged with its conjugate z k . Denote by (b k,1 , . . . , b k,n k ) and (b k+1,1 , . . . , b k+1,n k ) the subvectors of b associated with the roots z k and z k , respectively. Rearranging the columns of the matrix Z so that the original matrix Z is again reconstructed, the positions of the subvectors (b k,1 , . . . , b k,n k ) and (b k+1,1 , . . . , b k+1,n k ) are interchanged. Assuming that the solution of (16) is unique, we have the equality This means that the entries of b associated with a complex root z k and those of its conjugate z k are also conjugates. For the second statement, when a root z k is real, the same argument of taking conjugate of (16) and the assumption of a unique solution for the linear system can be used to obtain that the associated values b k,1 , . . . , b k,n k are such that b k,1 = b k,1 , . . . , b k,n k = b k,n k , that is, they are all real.
In Examples 2, 3 and 4 below we show particular instances of the system (16) where the two statements of the last proposition can be verified. In particular, for the unique positive real root z 1 we know that it is a simple root and its associated term b 1,1 , which can be written simply as b 1 , is real. We suspect b 1 is a positive number, however, a proof of this claim does not seem to be immediate and would require further technical analysis of the linear system (16).

Conjecture 1.
The coefficient b 1 associated with the unique positive root z 1 of the characteristic polynomial (13) and determined by the linear system (16) This claim will be numerically verified for some examples in Section 6 and will be used to propose an approximation for the ruin probability in Section 5. (17) below shows a section of the system (16) related with two conjugate roots z k and z k+1 = z k with common multiplicity n = 3. The entries of the vector b m×1 related with those roots are written as (b k,1 , b k,2 , b k, 3 ) and

Example 2. The equation
As explained before, taking conjugate of this system and rearranging the columns of the matrix so that the original matrix Z is again obtained, the position of the values b k,1 , b k,2 , b k,3 and b k+1,1 , b k+1,2 , b k+1,3 are interchanged. Then the assumption of a unique solution to the system implies that the identity

Example 3.
Here we consider the case m = 7 with = 3 roots z 1 , z 2 and z 3 , with multiplicities n 1 = 1, n 2 = 4 and n 3 = 2, respectively. Then the sysem (16) reads The vertical and horizontal lines on the left-hand side of (18) separate the terms associated with the 3 different roots. Observe that the roots z 2 and z 3 are not conjugates, otherwise, they would have the same multiplicity. Moreover, z 2 and z 3 are not complex since their conjugates do not appear in the linear system.

Example 4.
Consider the case when all the roots have multiplicity 1, i.e. = m and n 1 = · · · = n = 1. The relation (15) reduces to C n = m k=1 b k z n k , n ≥ 0, where b k is b k,1 and the linear system (16) is shown below in the equation (19). Observe the matrix of this system is the transpose of a complex Vandermonde matrix, see [21], ⎛ This simple case is studied in [35], where some theoretical and numerical examples are given. The formula derived for the ruin probability, in this case, has the rather manageable form This is a ruin probability formula for the risk process (1) with claims following a finite Erlang mixture distribution defined as in (9) and under the assumption that the underlying characteristic polynomial (13) only has simple roots.
We are now ready to state and prove our main result. This is a general version of (20) where now the roots of the characteristic polynomial have no constraints on their multiplicity. (1), let the claim sizes have a distribution as in (9). Let z 1 , . . . , z , 1 ≤ ≤ m, be the zeroes of (13) with multiplicities  n 1 , . . . , n , respectively, and let b 1 , b 2,1 , . . . , b 2,n 2 , . . . , b ,1 , . . . , b ,n be the solution of (16). Then

Theorem 1. For the Cramér-Lundberg model
where H j (z k , β, u) = 1 for j = 0, and and

Proof.
We make use of the following extension of Philipson's formula, see [33], which provides an expression for the moments of a Poisson distribution. Let λ be a real or complex number, for any integer j ≥ 1, where (23) can be found in the Appendix. Substituting the expression (15) for C n in (6) yields Observe that we have separated the cases k = 1 and k ≥ 2. In the first case, since n 1 = 1, the sum over j reduces to the constant b 1,1 . Using equation (23), and F k (x, β, u) = x v=0 (−z k βu) v /v! for x = 0, 1, . . . . Formula (21) and its simple case (20) show that the ruin probability can be expressed as a finite sum of exponentials of functions of the roots z k , multiplied by a nontrivial polynomial in z k . It is very interesting to observe that those exponentials are the probability generating functions of a Poisson(βu) distribution evaluated at z k given by H j (z k , β, u), and that there are as many exponentials as there are different roots in the characteristic polynomial (13). Each root z k (real or complex) contributes to the ruin probability ψ(u) and it is straightforward to show that the right-hand side of (21) is equal to its conjugate so that the whole formula yields a real number. In particular, the contribution of the unique positive root z 1 is given by the real number b 1 e −βu(1−z 1 ) , which is a positive number provided the conjecture b 1 > 0 holds. This will turn out to be the leading term in (21) and will be proposed as an approximation for ψ(u) in Section 5 below. It is also reassuring to notice that in the case when all the roots z k are simple, we have = m and n 1 = · · · = n m = 1 and the formula (21) reduces to the known solution (20), where the constants b k,1 are the solution b k to the system of equations (19). Observe also that, even with the assumption that the counting r.v. N has bounded support {1, . . . , m}, the number of summands in the ruin probability formula (6) is always infinite, whereas the new formulas (20) and (21) involve only m terms.
As an example, it is shown in [35] that a finiteErM(π, β) distribution can be seen as a phase-type distribution and our solution (20), in the case of simple roots, is the same as the well-known formula for ψ(u) for phase-type claims, see [2]. In the general case, when not all roots are simple, it is still an open problem to show that the standard formula for ψ(u) for phase-type claims is the same as our formula (21). The difficulty here lies in finding the exponential of a matrix that is not necessarily similar to a diagonal matrix.

An approximation
Formula (21) looks cumbersome as the whole expression contains four nested sums. However, it is an explicit expression of the ruin probability and its computer implementation is rather straightforward. In Section 6 below we give some numerical results to show particular values of the roots of the characteristic polynomial, the solution to the linear system, and the exact probability of ruin. The empirical results obtained from those computer experiments suggest the following ideas: consider an approximate solution to the linear system (16) of the form (b 1 , 0, . . . , 0), see (18) for a particular example of the matrix Z, then the first two equations of the system yields More explicitly, using (8) and the identity f N e (1) = 1/E(N ), the approximate value of z 1 is Thus, when knowing the exact values of z 1 and b 1 we have the following approximation.

Corollary 1. For the Cramér-Lundberg model (1) with finiteErM(π N , β) claims,
This is a Lundberg-type approximation and, as mentioned before, is the contribution of the unique positive root z 1 to ψ(u) in the formula (21). On the other hand, when z 1 and b 1 are unknown, the approximations (24) and (25) can be used to obtain the following estimate.

Corollary 2. For the Cramér-Lundber model (1) with finiteErM(π N , β) claims,
It is interesting to observe that in the case of Exp(β) claims (N = 1), the approximation (27) is the same as the exact ruin probability ψ(u) given before in (10).

Numerical examples
In this section, we show some numerical results of the calculation of the exact values of ψ(u) using our main formula (21) and the two approximationsψ 1 (u) andψ 2 (u) given in (26) and (27). We use R commands, see [4] and [34], to find the zeroes of (13) and to solve (16). Due to lack of space, we present only two examples. All of them show the case when some of the roots are complex and have a multiplicity larger than 1. In [35]  The distribution of the counting r.v. N is = (0.99858, 1.4073 · 10 −3 , 1.2457 · 10 −5 , 9.9605 · 10 −9 , 3.1827 · 10 −11 ), (29) where α 6 := 0. Also, C 0 = 5 j =1 α j = 0.01294389, then θ = C  5) with multiplicities, respectively, n 1 = 1, n 2 = 2 and n 3 = 2. The system (16) and the vector C 5×1 and the solution b 5×1 are given by Observe that the first term b 1 is positive and is much more significant than any of the other entries of b. As z 2 and z 3 are complex conjugates, their associated values b 2,1 , b 2,2 and b 3,1 , b 3,2 are also complex conjugates, respectively. By Theorem 1, where the function H is obtained using (22). Table 1 below shows the exact values of ψ(u) for different values of u, the values taken by the function H and the approximationsψ Table 1. (Example 5) Ruin probabilities ψ(u), the values of the function H and the approximationsψ 1 (u) andψ 2 (u) in the case of finiteErM(π, 1/10) claims and π given by (29) u The distribution of the counting r.v. N is The characteristic polynomial is with five distinct roots z 1 = 2/3, 3) and z 5 = −1/6 with multiplicities n 1 = 1, n 2 = 1, n 3 = 2, n 4 = 2 and n 5 = 1, respectively. The system of linear equations (16) is  (13) in the complex plane and the ruin probabilities for α 1 , . . . , α 5 given in (28) where and the vector C 7×1 and the solution b 7×1 are given by Observe again that the first term b 1 is positive and is much more significant than any of the other entries of b. As z 3 and z 4 are complex conjugates, their associated values b 3,1 , b 3,2 and b 4,1 , b 4,2 are also complex conjugates, respectively. By Theorem 1, the ruin probability reads where the function H is obtained using (22). Table 2 below shows the exact values of ψ(u) for different values of u, the values taken by the function H and the approximationsψ These approximations are again fairly accurate as the values of b k,1 , b k,2 multiplied by e −βu(1−z k ) , for k = 2, . . . , 5 are rather small. Table 2. (Example 6) Ruin probabilities ψ(u), the values of the function H and the approximationsψ 1 (u) andψ 2 (u) in the case of finiteErM(π, 1/4) claims and π given by (31) u

Remark 1.
As a subjective comparison of the two approximationsψ 1 (u) andψ 2 (u) and the numerical examples shown above, we observe that the approximationψ 2 (u) seems to be better thanψ 1 (u) for small values of the initial capital u, the opposite occurs for large values of u. Also, there seems to be no apparent difference in the quality of the two approximations when the values of the parameters θ and m are modified. As expected, the ruin probability is larger in the case when E(N ) is large and it seems there is no significant dependence on the value of m. It is also numerically observed that the ruin probability decreases as the surcharge factor θ increases.

Conclusions
We have given an alternative approach to calculating ψ(u) in the classical Cramér-Lundberg risk process (1) when claims are distributed according to a finite mixture of Erlang distributions. The method is based on the theory of recurrence sequences and requires to solve such a sequence, which in turn called for finding the roots of its characteristic polynomial and the solution of a linear system. One advantage of this procedure is that ψ(u) can be seen as a finite sum of as many exponential terms as the degree of the characteristic polynomial, see formulas (20) and (21), where there is one term for each root. A mathematical computer program such as R, see [34], can be used to solve the system (16) and to find the zeroes of (13), although the problem might be computationally challenging when the parameter m is large. Indeed, this is the case when an arbitrary claim distribution is approximated by a finiteErM(π, β) distribution. This is an approximation problem we have left aside for future work. On the other hand, further studies are also needed to quantify the accuracy of the approximations (26) and (27). It must be emphasized that the present work is a sequel of [35]. In that article, we presented the case of simple roots, whereas in this work we have considered the general case. On the other hand, using the same procedure shown here, it is possible to find a closed-form formula for the ruin probability in a discrete-time risk process with claims following a finite mixture of negative binomial distributions, see [36]. It is interesting to observe that all the ruin probabilities formulas found in those works are similar in the sense that they are expressed as linear combinations of probability generating functions. For the continuous risk model, the probability-generating function is that of a Poisson(βu) distribution evaluated at each of the roots z k , whereas for the discrete-time risk process, the probability-generating function is that of a negative binomial distribution.
Since finding the roots of polynomials of the form (13) is crucial for the ruin probability formulas (20) and (21), it is only natural to further the study of this type of polynomial and the solution to the linear system (16). Any knowledge about those roots might lead to an approximation or bound for ψ(u). Finally, it must be recalled that some quantities of interest in the mathematical theory of risk, such as the ruin probability, are difficult to calculate but they sometimes satisfy a recursive relation. Hence, the theory of recurrence sequences might as well be applied in many more instances within the theory of risk following the procedure shown in this work.