Ruin Probability Approximations in Sparre Andersen Models with Completely Monotone Claims

We consider the Sparre Andersen risk process with interclaim times that belong to the class of distributions with rational Laplace transform. We construct error bounds for the ruin probability based on the Pollaczek–Khintchine formula, and develop an efficient algorithm to approximate the ruin probability for completely monotone claim size distributions. Our algorithm improves earlier results and can be tailored towards achieving a predetermined accuracy of the approximation.


Introduction
The Sparre Andersen model is a classical object of study in insurance risk theory, see e.g. [17,19,21,22,27] and [5] for an overview. In this model, claims occur according to a renewal process, which generalizes the Cramér-Lundberg model, where claims arrive according to a Poisson process. Ruin probabilities in such a general setting are typically expressed as solutions of defective renewal equations, differential equations, the so-called Wiener-Hopf factorisation, etc, but the latter are typically inadequate to be used for numerical computations. However, if either the inter-claim times or the claim sizes belong to the class of phase-type distributions, then ruin-related quantities can be found in an explicit form; see e.g. [1,9,19] and [18], respectively. However, in many relevant situations in practice, the behaviour of the claim sizes is better captured by heavy-tailed distributions [13], but in that case explicit expressions are hard or impossible to evaluate even in terms of Laplace transforms. Under a heavy-tailed setting, a standard approach is hence to seek for asymptotic approximations [2,11,25], for initial capital levels being very large. At the same time, this capital level typically has to be very large in order to be reasonably accurate, when actual magnitudes matter. One mathematically appealing solution is then to look for higher-order approximations (see e.g. [4]), but also then an actual error bound for fixed values can not be given.
Another alternative is to approximate the actual heavy-tailed claim distribution by a tractable lighttailed one and control the introduced error in some way. Spectral approximations in this spirit were recently developed in [24] for the classical Cramér-Lundberg model.
The present paper proposes an extension of techniques in [24] to the more general Sparre Andersen model, and at the same time improves the bound derived there and the efficiency of the algorithm to establish it. Using the geometric compound tail representation of the ruin probability, we derive our error bound in terms of the ladder height distribution, which is explicitly available when the distribution of the inter-claim times has a rational Laplace transform. We focus on heavy-tailed claim sizes, where numerical evaluations of ruin probabilities are typically challenging and we develop an algorithm for the class of completely monotone distributions. Concretely, we approximate the ladder height distribution by a hyper-exponential distribution, and we are able to prescribe the number of required phases for a desired resulting accuracy for the ruin probability.
The rest of the paper is organised as follows. In Section 2, we introduce the model and provide the exact formula for the ladder height distribution. As a next step, we derive in Section 3 the error bound for the ruin probability and we construct our approximation algorithm. Section 4 compares our approximations with existing asymptotic approximations. In Section 5 we then perform an extensive numerical analysis to check the tightness of the bound and the quality of the derived approximations. Finally, we conclude in Section 6.

Model description
Consider the Sparre Andersen risk model for an insurance surplus process defined as where u ≥ 0 is the initial capital, c > 0 is the constant premium rate, and the i.i.d. positive random variables {X i } i≥1 with distribution function F X represent the claim sizes. The counting process {N (t), t ≥ 0} denotes the number of claims within [0, t] and is defined as N (t) = max{n ∈ N : where the inter-claim times W i are assumed to be i.i.d. with common distribution function K, independent of the claim sizes; see e.g. [5]. We also assume cEW > EX, providing a positive safety loading condition. Now, let T = inf{t ≥ 0 : R(t) < 0} be the time of ultimate ruin. Then, the ruin probability is defined as The ruin probability satisfies the defective renewal equation where φ = ψ(0), H(u) is the distribution of the ascending ladder height associated with the surplus process S(t) := u − R(t), andH(u) = 1 − H(u), for u ≥ 0; see e.g. [29]. The solution to Equation (3) is i.e. ψ(u) is a geometric compound tail with geometric parameter φ; see Section 1.2.3 in [28] for details.
Although Equation (4) provides a closed-form formula for the ruin probability, it is impractical because the ladder height distribution H(u) is not available in most cases of interest. However, when the distribution K of the inter-claim times has a rational Laplace transform, H(u) has an explicit form [19], which we recall in the next subsection. In the sequel we will then use this as a starting point for developing highly accurate approximations for ψ(u), which is of particular interest for heavy-tailed claim sizes.

The ladder height distribution with inter-claim times of rational Laplace transform
We assume now that the Laplace transform of the inter-claim times is a rational function of the form where µ n > 0, ∀n = 1, . . . , N , µ * = N n=1 µ n , and β(s) is a polynomial of degree N −2 or less. Obviously, of the claim sizes, it was shown in [19] that the generalised Lundberg equation has exactly N roots ρ 1 , ρ 2 ,. . . , ρ N , with ρ N = 0 and (ρ n ) > 0, n = 1, 2, . . . , N − 1. These roots play an important role in the evaluation of the ladder height distribution and the geometric parameter φ.
Denote withF X (x) the complementary cumulative distribution function (ccdf) of the claim sizes and consider the Dickson-Hipp operator for a function f (x) (see [10]). Moreover, let ρ * = N −1 n=1 ρ n . Then, as shown in [19], the ccdf of the ascending ladder heights is calculated via the formulā where Although the ladder height distribution in this model has an explicit formula, it is hard to evaluate ψ(u) either via Equation (4) or by taking Laplace transforms (an equivalent formula to the Pollaczek-Khinchine in the Cramér-Lundberg model). This is in particular the case when the claim sizes follow a heavy-tailed distribution, as already mentioned in Section 1. As a result, in such cases, opting for approximations seems a natural solution.
In the next section, we will study error bounds for ψ(u) when the ladder height distribution is approximated by a phase-type distribution. In particular, we will provide an efficient algorithm to construct approximations for ψ(u) when approximating H(u) by the subclass of hyper-exponential distributions.

Spectral approximation for the ruin probability
Starting point for the approximation of ψ(u) is its geometric compound tail representation in Equation (4). Note that this representation is similar to the Pollaczek-Khintchine formula for ψ(u) in the Cramér-Lundberg model where φ is replaced by the average amount of claim per unit time ρ < 1 and the ladder height distribution is equal to the stationary excess claim size distribution. Therefore, following the reasoning in [24], we will approximate the ladder height distribution by a hyper-exponential distribution (which has a rational Laplace transform), to construct approximations for the ruin probability.

Error bound for the ruin probability
LetĤ(u) be an approximation of the ladder height distribution H(u) andψ(u) be the exact result we obtain from (4) when we useĤ(u). From Equation (4) and the triangle inequality, the error between the ruin probability and its approximation then is If we define the sup norm distance between two distribution functions F 1 and F 2 as x ≥ 0, the following result holds: A bound for the approximation error of the ruin probability is Proof. The result is a direct application of Theorem 4.1 of [20] by (i) choosing the functionsF 1 ,F 2 to be H andĤ, (ii) taking ρ = φ, and (iii) recognising that sup y<u H(y) −Ĥ(y) ≤ D(H,Ĥ).
To sum up, when the ladder height distribution is approximated with some desired accuracy, then a bound for the ruin probability is guaranteed by Theorem 3.1. While this result holds for any approximationĤ of H, we will in the sequel focus on hyper-exponential approximations since these lead to very tractable expressions and at the same time are sufficiently accurate for the purpose.
Consequently, our next goal is to construct an algorithm to approximate the ladder height distribution by a hyper-exponential distribution.

Completely monotone claim sizes
We are mostly interested in evaluating ruin probabilities when the claim sizes follow a heavy-tailed distribution, such as Pareto or Weibull. These two distributions belong to the class of completely monotone distributions. Completely monotone distributions can be approximated arbitrarily closely by hyper-exponentials; see e.g. [14]. Here, we provide a method to approximate a completely monotone ladder height distribution with a hyper-exponential one to achieve any desired accuracy for the ruin probability. The following result is standard; see [15].
We call S the spectral cdf.
Remark 3.5. With a slight abuse of terminology, we will say that a function S is the spectral cdf of a distribution if it is the spectral cdf of its ccdf.
Note that Theorem 3.4 extends also to the case where S(y) is not a distribution but simply a finite measure on the positive half line, i.e. a function f is completely monotone if and only if it can be expressed as the Laplace-Stieltjes integral of such a finite measure S(y). We will show that under the assumption that the claim size distribution is c.m., the ladder height distribution is c.m. too. We first need the following intermediate result. Proof. Assume that the claim sizes are completely monotone, i.e.F X (u) = ∞ 0 e −uy dS(y), for some spectral cdf S(y). In this case, it holds that where dS Tρ n (y) = dS(y) y+ρn , n = 1, . . . , N , is a finite measure on the positive half line with S Tρ n (+∞) = (1 −f X (ρ n ))/ρ n , n = 1, . . . , N − 1, and S T 0 (+∞) = EX.
We can now state the following result.
Proof. It was proven in [8] that the ascending ladder height distribution in the Sparre Andersen model is c.m. if the claim size distribution is c.m. This means thatH(u) can be represented as the Laplace- Stieltjes transform of some spectal cdf S H (y). Due to the uniqueness of Laplace transforms, it hence suffices to find the formula of the spectral cdf S H (y) by applying Lemma 3.6 to (8).
We show in the next section how to utilise the above results to construct approximations for the ruin probability ψ(u) that have a guaranteed error bound given by Theorem 3.1.

Approximation algorithm
Following the proof of Lemma 2 in [24], we can directly deduce the following result. The above lemma states that if we want to approximate a c.m. ladder height distribution with a hyperexponential one with some fixed accuracy , it suffices to approximate its spectral cdf with a step function with the same accuracy. As pointed out in Remark 1 of [24], we could approximate S H with a step function having k jumps that occur at the quantiles λ i , such that S H (λ i ) = i/(k + 1), i = 1, . . . , k, and are all of size 1/k to achieve D(H,Ĥ) ≤ = 1/(k + 1). Another possibility is to use the step function in Step 4d of our Algorithm; see also Figure 1 for a graphical representation of the approximate step function and its corresponding hyper-exponential distribution. Clearly, this new step function leads to D(H,Ĥ) ≤ = 1/2(k − 1).
The error bound for the approximate ruin probabilityψ(u) can be calculated afterwards through Theorem 3.1. An interesting question in this context is how many phases k for the approximate ladder height distribution suffice to guarantee an error bound ψ(u) −ψ(u) ≤ δ for some pre-determined δ > 0. We answer this question in the next lemma. Lemma 3.9. To achieve ψ(u) −ψ(u) ≤ δ for some pre-determined δ > 0, the ladder height distribution H(u) must be approximated by a hyper-exponential one with at least k phases such that where x is the integer that is greater than or equal to x but smaller than x + 1.
Proof. Observe that the error bound in Theorem 3.1 depends on the approximate hyper-exponential distributionĤ(u), which means that one should first determineĤ(u) and then calculate the error bound. However, when D(H,Ĥ) ≤ , this translates to H(u) − ≤Ĥ(u) ≤ H(u) + . Therefore, the worst case scenario for the bound is whenĤ(u) = H(u) + and consequently D(H,Ĥ) = . As a result, if we want to achieve ψ(u) −ψ(u) ≤ δ for all possible scenarios ofĤ(u), we should solve the inequality with respect to . By substituting = 1/2(k − 1), we calculate In addition, the bound is asymptotically equal to φ/(1 − φ) according to Remark 3.2. Consequently, it must also hold that Finally, as the number of phases k must be an integer, the smallest possible integer that satisfies at least one of the inequalities is the one described in Equation (12).
After this, we present our algorithm under the setting that we fix the desired accuracy δ for the approximation of the ruin probabilityψ(u).
2. Find the spectral cdf S(y) ofF X (x).
3. Use Proposition 3.7 to calculate the spectral cdf S H (y) ofH(u).
4. ApproximateH(u) by a hyper-exponential distribution with k phases.
(a) Choose the accuracy of the ruin probability δ for a fixed u > 0.
(b) Calculate k required to achieve this accuracy using Lemma 3.9 and set = 1 2(k − 1) .

The accuracy forψ(u) is then
Remark 3.10. The decomposition of L ψ (u) (s) at Step 6 is guaranteed by [6] who showed that the ruin probability in the Sparre Andersen model has a phase-type representation when the claim sizes are phase-type. Moreover, the particular hyper-exponential representation ofψ(u) at Step 7 is due to the fact that the poles of L ψ (u) (s) are exactly the roots of the polynomial function , where δ ij is the Kronecker delta. It is immediate from perturbation theory that P φ (s) has exactly k simple roots analytic in φ; see [7] for details.
Remark 3.11. The above algorithm is an extension of the one developed for the Cramér-Lundberg model in [24], to which we refer for further details on technical implementation.

Asymptotic approximation
In many cases, it is of importance to investigate the asymptotic behaviour of the ruin probability when the initial risk reserve tends to infinity. This question is particularly interesting in the case of heavytailed claim sizes. Towards this direction, when the claim sizes belong to the class of subexponential distributions S [23], e.g. Pareto, Weibull, Lognormal, etc, the following asymptotic approximation is classical (see e.g. [12]): Theorem 4.1. Suppose in the general Sparre Andersen model that the claim sizes and inter-claim times have both finite means, EX and EW, respectively, such that cEW > EX.
Note that the heavy tail approximation ψ S (u) holds for any inter-claim time distribution. However, further modifications have been attained in [26] when the Laplace transform of the inter-claim times is a rational function of the form (5) with β(s) = β and F X belongs to the subclass of regularly varying distributions, i.e.F X (u) ∼ L(u)u −α−1 e −γu , u → +∞, where L(u) a slowly varying function and α > 0, γ ≥ 0. For example, the Pareto(a, b) distribution (see Section 5.2.1) belongs to the class of regularly varying distributions with L(u) = b + 1/u −a , α = a − 1, and γ = 0, and its modified asymptotic approximation is then given by: , which is smaller than ψ S (u) by a factor bu bu + 1 that converges to 1 as u → +∞; see [26] for details.
Clearly, the heavy tail approximation admits a simple formula whenever the expectations of the interclaim times and claim sizes are finite. Its drawback though is that when cEW ≈ EX the approximation is useful only for extremely large values of u.
We compare in the next section the accuracy of the spectral approximation to the accuracy of the heavy tail one, i.e. ψ S (u). An interesting observation is that the spectral approximation converges faster to zero than any heavy-tailed distribution due to the exponential decay rate of the former. Thus, the heavy tail approximation is expected to outperform the spectral approximation in the far tail, but for medium values this new approximation can be very competitive.

Numerical analysis
The goal of this section is to implement our algorithm in order to check the accuracy of the spectral approximation and the tightness of its accompanying bound, which is given in Theorem 3.1. To perform the numerical examples, we need to make a selection for the distribution K of the inter-claim times as well as the claim size distribution F X .

Claim sizes
For the claim sizes, we consider the Pareto(a, b) distribution with shape parameter a > 0 and scale parameter b > 0 and the Weibull(c, a) distribution with c and a positive shape and scale parameters, respectively.

Pareto
This distribution is c.m. since its ccdfF X (x) = (1 + bu) −a can be written as the LST of the Gamma distribution with shape and scale parameters a and b respectively, i.e.
The nth moment of the Pareto distribution exists if and only if the shape parameter is greater than n. Since we are interested in comparing the spectral approximation to the asymptotic approximation of Section 4, it is necessary to have a finite first moment for the claim sizes. Therefore, the shape parameter a must be chosen to be greater than 1.
Using Proposition 3.7, we can easily verify that:

Weibull
It can be verified that the ccdfF X (x) = e −(u/a) c with fixed shape parameter c = 1/2 arises as a c.m.
distribution [16], where the mixing measure (measure of the spectral function) S is given by 2 aπy 3 dy.

Numerical results
The goal of this section is to implement our algorithm to check the accuracy of the spectral approximation and the tightness of its accompanying bound, which is given in Theorem 3.1. Impact of phases. It is intuitively true that the spectral approximation becomes more accurate as the number of phases increases. To test this hypothesis, we compare three different spectral approximations with number of phases 10, 30, and 100, respectively, with the exact value of the ruin probability (which we obtain through simulation). We display our results in Table 1 only for Pareto claim sizes. The conclusion is that indeed a more accurate spectral approximation is achieved as the number of phases increases for every fixed initial capital u, which is in line with expectations. Quality of the bound. A compelling question regarding the bound is if it is strict or pessimistic, i.e. how far it is from the true error of the spectral approximation. To answer this question, we first need to determine the accuracy δ we would like to guarantee for the ruin probability. Using Lemma 3.9, we present in Figure 2  As we can observe in Figure  Finally, notice that the predicted error bound is almost 4 times smaller than δ = 0.02 in the Pareto case.
This happens because D(H,Ĥ) could be a lot smaller than ; see also Figure 1 where D(H,Ĥ) < 0.1.
However, most importantly, the true error is close to the predicted bound and thus, we can say that Lemma 3.9 provides a good proxy for the necessary number of phases k to achieve it.
Comparison between spectral and heavy tail approximations. As we pointed out in Section 4, the spectral approximation is expected to underestimate both the exact ruin probability and the asymptotic approximation ψ S (u) in Theorem 4.1 for large u due to its exponential decay rate. It is of interest to see the magnitude of u for which the asymptotic approximation outperforms the spectral approximation.
We select the spectral approximations with k = 67 phases for Pareto(2, 3) claim sizes and k = 11 phases for Weibull(0.5, 3) claim sizes as in the previous experiment and present the distributions in a graph. The pink shadow in Figure 4 enfolding the spectral approximation represents its bound. We observe that for small values of u, the spectral approximation is more accurate than the heavy tail approximation, where the second provides a rough estimate of the ruin probability. On the other hand, the heavy tail approximation is slightly more accurate than the spectral approximation in the tail, i.e.
for u > 25, under Pareto claim sizes. However, for the Weibull case, we observe that even for values of u around 300, the spectral approximation still outperforms the heavy tail approximation.

Conclusions
In this paper, we considered the ruin probability of the Sparre Andersen model with heavy-tailed claim sizes and inter-claim times with rational Laplace transform. Using the geometric random sum representation, we developed an explicit bound and also constructed a spectral approximation by approximating the c.m. ladder height distribution with a hyper-exponential one. Our spectral approximation algorithm advances on the algorithm established in [24] in various aspects. We provide below a summary of our conclusions both for the spectral approximation and the bound: • When comparing with the technique proposed in [24], the strategic selection of the quantiles in Step 4d reduces the number of phases to almost a half in order to guarantee a certain accuracy for the ladder height distribution.
• Since the bound depends on the initial capital, this allows us to focus on one area and optimise the required number of phases to achieve a desired accuracy, e.g. we would need 110 phases for u = 5 and 132 phases for u = 30 to guarantee accuracy of at most δ = 0.01 in our example.
• The step function is constructed to guarantee D(H,Ĥ) ≤ , but in most applications D(H,Ĥ) is a lot smaller than . Thus, the use of D(H,Ĥ) in the bound makes it tighter.
To sum up, the spectral approximation is highly accurate for all values of u as opposed to the heavy tail approximation, which fails to provide a good fit for small values. Moreover, it is accompanied by a rather tight bound.
Finally, note that the results of this paper are also valid for the risk model with two-sided jumps, i.e.
where u, c, and X i are defined as before, while N + (t) and N − (t) are independent Poisson processes with intensities λ + and λ − , respectively; see e.g. [3]. In addition, the sequence {Y j } i≥1 of i.i.d. r.v.'s, independent of {X i } j≥1 , N + (t), and N − (t), and having the common d.f. G Y that belongs to the class of distributions with rational Laplace transform, are the sizes of premium payments. The positive security loading condition in this model becomes: c + λ + EY > λ − EX.