We here study the classical Cramér–Lundberg risk process { R ( t ) : t ≥ 0 }, see [2] or [14], given by (1) R ( t ) = u + c t − ∑ j = 1 N ( t ) X j , t ≥ 0 . The constant u ≥ 0 is the initial capital, c > 00$]]> is the premium per unit of time which is also assumed to be constant, N ( t ) is a Poisson process with intensity λ > 00$]]>, 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 }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 > 00$]]>. 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–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 > λ μ\lambda \mu $]]>. This critical condition can be written as c = ( 1 + θ ) λ μ for some θ > 00$]]>, 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 θ > 00$]]>. As it is customary, the survival probability P ( X > x )x)$]]> is denoted by F ‾ ( x ), and the letters d.f. refer to the distribution function of a random variable.

The density of the Erlang distribution is denoted by erlang ( k , β ) ( x ), that is, erlang ( k , β ) ( x ) = ( β x ) k − 1 ( k − 1 ) ! β e − β x , x > 0 . 0.\]]]> where k ≥ 1 is an integer and β > 00$]]>. 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 β > 00$]]> when its density is (2) f ( x ) = ∑ k = 1 ∞ f N ( k ) · erlang ( k , β ) ( x ) , x > 0 . 0.\]]]> 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 E ( S ) = E ( N ) / β. If F ( x ) is a continuous d.f. on [ 0 , ∞ ) and for each natural number n, one defines (3) F n ( x ) = ∑ k = 1 ∞ f n ( k ) Erlang ( k , n ) ( x ) , x > 0 . 0.\]]]> 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 (4) f e ( x ) = 1 μ [ 1 − F X ( x ) ] , x ≥ 0 .

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 (5) f N e ( k ) = F ‾ N ( k − 1 ) E ( N ) , k = 1 , 2 , … where F ‾ N ( k − 1 ) = P ( N > k − 1 )k-1)$]]>.

As in our previous work [35], the starting point in our analysis is the ruin probability formula (6) ψ ( u ) = ∑ n = 0 ∞ C ‾ n P ( Z β u = n ) , u > 0 , 0,\]]]> for the Cramér–Lundberg model (1), where Z β u ∼ Poisson ( β u ) and (7) C ‾ 0 = Ψ ( 0 ) = 1 / ( 1 + θ ) , (8) C ‾ n = C ‾ 0 ∑ j = 1 n f N e ( j ) C ‾ n − j + F ‾ N e ( n ) , n ≥ 1 , 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 > ⋯ {\overline{C}_{1}}>\cdots \hspace{0.1667em}$]]> 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 ) > 00$]]>. For distributions π of this type we call (2) a finite Erlang mixture distribution. That is denoted by finiteErM ( π , β ) and its density is given by (9) f ( x ) = ∑ k = 1 m f N ( k ) · erlang ( k , β ) ( x ) , x > 0 . 0.\]]]> 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).

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 (11) α j = C ‾ 0 f N e ( j ) = F ‾ N ( j − 1 ) ( 1 + θ ) E ( N ) , j = 1 , … , m , where, clearly, 1 > α 1 ≥ α 2 ≥ ⋯ ≥ α m − 1 ≥ α m > 0{\alpha _{1}}\ge {\alpha _{2}}\ge \cdots \ge {\alpha _{m-1}}\ge {\alpha _{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 (12) C ‾ m + n − ∑ j = 0 m − 1 α m − j C ‾ j + n = 0 , n = 0 , 1 , … where C ‾ 0 , … , C ‾ m − 1 are the initial values. The characteristic polynomial associated with the recurrence sequence (12), see [5], being (13) p ( y ) = y m − ∑ j = 0 m − 1 α m − j y j . 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).

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 > 00$]]> is the unique positive root which we know has multiplicity n 1 = 1. Thus, (13) can be written as (14) p ( y ) = ( y − z 1 ) 1 ( y − z 2 ) n 2 ⋯ ( y − z ℓ ) n ℓ . It can be shown, see [39, p. 55], that the solution of (12) is (15) C ‾ n = ∑ k = 1 ℓ ∑ j = 1 n k b k , j n j − 1 z k n , n = 0 , 1 , … , 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 (16) Z m × m · b m × 1 = C ‾ m × 1 , where b m × 1 = b 1 , 1 b 2 , 1 ⋮ b 2 , n 2 ⋮ b ℓ , 1 ⋮ b ℓ , n ℓ , C ‾ m × 1 = C ‾ 0 C ‾ 1 ⋮ C ‾ m − 1 , and the matrix Z has the form Z m × m = Z _ 1 Z _ 2 , 1 Z _ 2 , 2 ⋯ Z _ 2 , n 2 ⋯ Z _ ℓ , 1 Z _ ℓ , 2 ⋯ Z _ ℓ , n ℓ , where the vertical lines help to visually separate the column vectors Z _ 1 = 1 z 1 z 1 2 ⋮ z 1 m − 1 , Z _ k , j = 1 ( j = 1 ) z k 2 j − 1 z k 2 ⋮ ( m − 1 ) j − 1 z k m − 1 , 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 \cdots >{\overline{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.

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).

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.

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.

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 > 00$]]> 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.

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 (24) b 1 ≈ C ‾ 0 = 1 / ( 1 + θ ) , (25) z 1 ≈ C ‾ 1 / C ‾ 0 . More explicitly, using (8) and the identity f N e ( 1 ) = 1 / E ( N ), the approximate value of z 1 is C ‾ 1 / C ‾ 0 = f N e ( 1 ) C ‾ 0 + F ‾ N e ( 1 ) = 1 ( 1 + θ ) E ( N ) + 1 − 1 E ( N ) = 1 − θ ( 1 + θ ) E ( N ) . Thus, when knowing the exact values of z 1 and b 1 we have the following approximation.

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.

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).

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] we provided numerical examples in the case when all roots are simple. In all the examples shown below the conjecture b 1 > 00$]]> is fulfilled and the relations stated in Proposition 1 regarding the roots z k and their associated numbers b k , 1 , … , b k , n k are verified.