Ruin Probabilities and Complex Analysis

This paper considers the solution of the equations for ruin probabilities in inﬁnite continuous time. Using the Fourier Transform and certain results from the theory of complex functions, these solutions are obtained as complex integrals in a form which may be evaluated numerically by means of the inverse Fourier Transform. In addition the relationship between the results obtained for the continuous time cases, and those in the literature, are compared. Closed form ruin probabilities for the heavy tailed distributions: mixed exponential; Gamma (including Erlang); Lognormal; Weibull; and Pareto, are derived as a result (or computed to any degree of accuracy, and without the use of simulations).


Introduction
Throughout the history of the collective theory of risk the concept of ruin probability has played a pre-eminent role in quantifying the risk of an insurance business in terms of premiums and reserves. This is evidenced by a considerable body of literature devoted to the calculation of ruin probabilities. Surveys of results may be found [1], [2], [3]. The last mentioned includes an extensive bibliography, which we do not reproduce in this paper. Other papers [4], [5], [6] have discussed numerical or simulation methods for calculating ruin probabilities.
This paper is concerned with application of the theory of complex functions to deriving exact results for ruin probabilities (or as precise as may be desired). We are specifically concerned with solutions to the ruin equations in terms of complex integrals, or with numerical methods of evaluating those integrals to any required degree of accuracy. It is stated that [3, 1.4b] exact solutions are known in only a few cases, mainly for exponential or mixed exponential claim densities. Even then this paper shows how and why the results are incomplete, or even erroneous. It is not the purpose of this paper to review the expansive literature on this subject -see [3] for voluminous references. [7] puts it more eloquently: "Although for most of the general claim amount distributions, e.g., heavy-tailed, the Laplace transform technique does not work, explicit expressions under other assumptions, such as Pareto distributions, have been obtained but they are too complicated and require large computation to calculate the values of the ultimate ruin probability. For example, Garcia (2005) derived complicated exact solutions under series representation and Seal (1980) and Wei and Yang (2004) under integral representations. Grandell and Segerdahl (1971) showed that for the gamma claim amount distribution Risks 2019, 7, 68; /risks7020068 www.mdpi.com/journal/risks Risks 2019, 7, 68 2 of 16 under some restrictions on the parameters, the exact value of ruin probability can be computed via a formula which involves a complicated integral. In Ramsay (2003), an expression based on numerical integration was derived for the probability of ultimate ruin under the classical compound Poisson risk model, given an initial reserve of u in the case of Pareto individual claim amount distributions. Furthermore, Albrecher et al. (2011) have obtained closed-form expressions for ruin probability functions under some kind of dependence assumption also using the mixing representation. In this regard, as Asmussen and Albrecher (2010) pointed out, the ideal situation is to come up with closed-form solutions for the ruin probabilities; however, these are limited. More recently, Tamturk  It is apparent that exact solutions to the ruin equation are multifarious and rely almost always on approximations (via orthogonal functions) or asymptotics. The practical consequence is that claim density relies on the claim shape.
We present a different and novel approach using complex analysis, that allows the accuracy of the results to be specified. As we shall see, however, the theory of complex functions also provides means for exact evaluation of the solutions derived in certain natural cases. Furthermore there exist new methods of integration, arising in physical and other contexts, which may be of use in numerical evaluation for other cases. In deriving our results we demonstrate the contention that the theory of complex functions provides a natural framework for the investigation of ruin probabilities. Where appropriate, we demonstrate how these results relate to those which are already known.
The application of these results may be evident. Insurance is the primary and motivating factor. However as [8]suggests, water storage, solar electricity batteries, data retention are relevant. Ordnance and academic research itself are also candidates.
Complex function analysis is introduced by use of the Fourier Transform (FT). This has many advantages over the Laplace Transform, though the two are functionally equivalent. First, as noted by [9], the FT of a probability density always exists, unlike its Laplace counterpart. Second, analysis of probability distributions is aided by classical results in analysis in the complex plane. Third, the FT is numerically computable using modern software 1 The plan of this paper is as follows: • Section 2 introduces some preliminary concepts regarding ruin and derives the ruin equations from first principles; • Section 3 presents some preliminary mathematical results which are applied in later sections; • Section 4 discusses the solution of the ruin equation for various claim distributions; • Section 5 charts the ruin probabilities in a comparable way.
The general plan of this paper is to highlight results via a sequence of Propositions. Proofs are given in detail only where the references are insufficient for the particular point in question.

Basic Equations
Consider a risk business involving the following parameters: • P is the rate of premium received per unit time; • ζ is the stochastic variable measuring the amount of claim (given that a claim has occurred) with probability density function p (ζ); • u is the reserve held at any time t • ψ(u) is the probability of ruin of the business at any time after time t, where the initial reserve is u at time t; • ϕ(u) is the corresponding probability of survival, with ψ+ϕ = 1.
It will be seen that, since claim amount ζ is part of the change in reserve u, the symbols ζ and u will be used interchangeably.
In this paper we are concerned only with claim processes which are compound Poisson distributed, that is, where the probability of a single claim follows a Poisson distribution, and the amount of each claim ζ is identically and independently distributed according to the density p (ζ). Without loss of generality we may assume that time may be scaled so that the Poisson parameter is 1.
We define L 1 (R) as the space of Lebesgue integrable functions with finite norm Throughout this paper all integrals are taken to be defined in the sense of Lebesgue unless otherwise specified. We also consider later (in Appendix A, in connection with the inverse Fourier Transform) the space of square integrable functions L 2 (R) with norm If f ∈ L 1 (R) and f is bounded, then it is clear that f ∈ L 2 (R).
In general we require that the claim amount density p (ζ) satisfy the conditions p ≥ 0 for ζ ≥ 0, and p = 0 for ζ < 0. In addition we require that the claim amount density generally satisfy These conditions are to ensure that the probability density of claims is sensible, and that it has a finite mean and variance. Additional restrictions on p (ζ) will be imposed as required.
Without loss of generality we may scale the claim amount ζ so that the mean claim is 1 and ζ f 1 = 1. In practice this means that we take always the gross premium rate P > 1.

Ruin in Continuous Time,
Let ϕ(t, u) be the probability of survival. For a small time interval dt, we have: since claim frequency is Poisson distributed with parameter 1. In addition we have a finite probability e −dt that there will be no claims. Approximating e −dt = 1 − dt for small dt, and ignoring terms of order O dt 2 or higher, we get from the following equality: This is justified if there are no claims the reserve increases with the premium Pdt.
If, on the other hand there are claims, with probability dt, then the reserve u > v will fall to v by an amount of claim u − v with probability p(u − v).Passing to the limit dt→ ∞, rearranging and dividing through by dt, we get the integrodifferential equation This is the equation governing ruin in finite time with continuous testing.
If we move to the limit t → ∞ then we should expect ϕ t → 0 for ϕ to be stationary at ∞. Hence the equation governing ruin in infinite time with continuous testing is: This is the same equation as in [10, 8.8] and is the equivalent to the integral equation often given for ruin [2, ex 2.5.11], [11, eqn 2.4].
Hence in terms of the ruin probability ψ (u) = 1 − ϕ(u) is where g (u) = v>u p(v)dv and 1 − g (u) is the cumulative claim density.

Proposition 1
If the FTp has a finite mean, the ruin for zero reserve is ψ (0) = 1 P .
The ruin probability at u = 0 turns out to be critical. First integrate 3 for u over (0, ∞) , (using the convention that only the area of integration needs to be specified, as Fubini's theorem makes their order irrelevant). The double integral becomes The above shows that for any p ∈ L 1 (R).
We denote ψ (0) as ψ 0 for convenience. The FT of g (u) can be found from the general identityf Thus ψ (0) = 1 P for all p.
It remains to show thatψ ∈ L 2 (R) for application of the inverse FT for real u = 0, by Proposition 2. We first show that is bounded at z = 0. This follows from 5.
In fact the value ofψ at z = 0 is lim involves the second moment of p.
Appendix B provides some simple implications of these results, which are well known in the literature (albeit with more complicated proofs).

Exact Evaluation of ψ(u)
We now present a complex analysis approach to respond to the comments in [7] (see Introduction.) To solve equation 3 for the ruin probability ψ we need to have suitable boundary conditions. It is plausible that survival must be certain as the initial reserve u → ∞, so that a suitable condition would be ψ (∞) = 0, which is well known . This condition is designed to ensure that the resulting solution is physically meaningful.
Taking the FT of 1 and using the relationsψ In addition This may be formalized as follows.

Proposition 3
The solution ψ ∈ L 1 (R) to the ruin equation 5 for infinite time is given by the inverse FT for ψ, wherê The FTψ may be computed in several ways.

The Hermitian Property
The Hermitian property (HP) for a complex function f (z) satisfies There are many examples given later in this paper. Most them are the FT transform of a real function or one of them (e.g.p (z) , η(z), g (z)) Most important is the FT of the ruin probabilityψ (z) .
We note thatψ (z) has the Hermitian property provide and its the following this simplification in numerical and other contexts.
so that ψ is always real.

Direct Integration
The inverse FT integrals given above do not exist in the L 1 (R) sense but need to be interpreted as principal values or improper integrals [13], 1966], for example: where the integral is taken as improper Thus ψ may be numerically computed as

Cauchy residues
Asĝ = 1 has a simple zero at z = 0, the poles of 1 η(z) depend on the zeros of η(z), unless z = 0 is a double zero (which does not arise in the cases considered in this paper). It remains to show thatψ ∈ L 2 (R) for application of the inverse FT for real u = 0, by Proposition 2. We first show that is bounded at z = 0.
Appendix B provides some simple implications of these results, which are well known in the literature (albeit with more complicated proofs).
In many cases the inverse FT for ψ may be evaluated exactly by means of the Cauchy residue theorem [14, ch 7].
Thus ψ(u) is achieved by writing the inverse FT from 7 as: We may then integrate the function e −iuz η(z) (z) along a contour consisting of the semicircle Γ with center z = 0 and radius R in the lower half plane C − , together with that part of the real axis bounded by Γ.
The functionψ, is well behaved at z = 0, and moreover ψ → 0 uniformly along Γ as R→ ∞, sinceψ is analytic at ∞ and tends to 0 there. Hence by Jordan's lemma [14, 7.9] the contribution of the integral along Γ vanishes as r→ ∞.
Note that ω = 0 may be regarded as a removable singularity as the term η(z) may be expanded around z = 0 as: Integrating along a small semi-circle of radius r around z = 0 and using z = re iθ 1 2π which tends to 0 as the radius r of γ tends to 0. This produces the following simplification: The Cauchy residue theorem implies that ψ is given in terms of the residues of the integrand in 10 as follows: where the sum is taken over the singularities ω ∈ C − , ω = 0 of e −iuz η(z) with corresponding residues ρ (ω) .
If η(z) is rational, the singularities ω can be only poles, which occur precisely at the zeroes of η(z). If z = ω is a simple zero of η(z), then the residue at this point is defined as The formula for ψ in equation 11 must give real results, even if complex residues are evident. If z = ω is a root of η(z), then so too is z = −ω with residue ρ (−ω) = −ρ (ω). Hence ψ may be written alternatively as: where the first sum is taken over all purely imaginary roots iξ and the second sum is taken over all roots ω ∈ C − with (ω) > 0.
It is important to note the standard techniques applied above, involving Jordan's lemma and Cauchy's theorem, for they will be applied without further ado in the rest of this paper.

Partial Factions
In some cases, finding all the roots of η(z) is prone to error; the Matlab software, or its Chebfun extension, is not always reliable. In this case we express η(z) as a polynomial, or the truncated portion of a fast converging series. The software is very reliable with finding roots of polynomials, but less so with non-linear functions. The less the truncation, the more accurate the results, but the greater the number of roots to contend with.
An alternative technique may be used 3 using partial fractions for η(z) if it can be rationalized: .where A and B are polynomials. This may be expressed as partial fractions: using, for example the residue function in Matlab. Then the inverse FT of each term is, using Jordan's lemma:  The direct integration technique may always be used, though it is computationally cumbersome. The Cauchy and partial fraction techniques are similar in that both require finding the roots of η(z). However the partial fraction technique requires η(z) to be expressed as polynomials (or approximations thereof). The Cauchy technique requires only that the roots of η(z) be found, for which many methods, and many roots, are possible.

Various Claim Densities
In this section we consider claim densities with varying tails. Numerical results are charted in section 5 for comparison. Where possible. the parameters of the density are chosen to give unit mean and variance.
In addition, the roots of η and their parameters are shown only for those in C − with (z) > 0. The Hermitian principle provides the other roots.

Exponential Claims
As the first application of equation 12 consider the case of exponentially distributed claims, so that p = e −u andp(z) = i i+z . In this case η(z) has zeros at 0 and another at iξ : It is also clear that η(z) has a simple zero at z = 0 and at z = −i 1 − 1 P where i + z = i P andp(z) = P. Thus

Mixed Exponentials
As a slight extension of 13 consider the case of mixtures of exponentials, that is for some range of k. This corresponds to: Thus for a probability distribution with unit mean and unit variancê and for unit mean We take β k = [1 2 3] and find α k accordingly. We thus have to find the roots of This equation can be rationalized by multiplying by R (z) = ∏ k (z + β k i) , which is a well defined polynomial, so that the equation becomes which again is a well defined polynomial. Standard Matlab functions may then be applied. The solution for ψ is then a sum of the form i(P − 1) e −iuω η (ω) , the coefficients of e −iuω being evaluated in terms of the zeros ω [6].
The zeros of η then correspond to the α non-zero roots of (1 + iPz)(1 − iz/α) α = 1. In this case ψ may also be expressed as a finite sum of exponentials, as in §13.
We now turn to more complicated examples, with heavier tails.

The Gamma Case
[15] provides a complicated and approximate approach using Mittag-Leffler functions to the Gamma case. However another simple application of equation 12 may be found for the gamma distribution, namely from which we obtainp = (1 − iz/β) −α , which again is a rational function. For integer α > 1, it is known as the Erlang distribution.
For unit mean and variance we need to have α = β = 1. However this is the same as the exponential. However this can be approached generally. The roots of If α = m n can be rationalized, we need the roots of (1 + iPz) m (1 − iz/β) n = 1 As a simple and particular example, take α = β = 2. We need the roots of Disregarding the root z = 0, we find the non-zero roots and residues: This provides an interesting contrast to the result in 13.

Remark 6
The roots of η correspond are those of 1 1−iz − 1 − iPz in this example. However the Chebfun package in Matlab produces several more roots for this function, some of which are spurious. This suggests that rootfinding using software needs to be approached with caution. The safest approach is to employ rationalization, as Matlab is reliable in finding all roots of polynomials, including those with multiplicity.
Remark 7 [15] also illustrates this case, but with different parameters. While we have used (Basis A) Under our method, this is analyzed as Basis B. For completeness' sake, the roots are as follows.
A comparison of results is as follows

The Weibull Case
The Weibull density for claims is For unit mean and variance, we have λ = 0.8767 and k = 0.80647. The partial fractions are shown below.

The Lévy Case
The Lévy distribution has infinite mean and variance, with density In this case direct integration by 12 is possible and may be computed below. As it has infinite mean, we take µ = 0 and ψ (0) = 1 P , leading to comparability with the other densities in this section.

The Lognormal case
For a variance parameter σ 2 and mean µ the pdf of the lognormal distribution may be written asp For unit mean and variance. we must µ = 1 2 and σ = 1. Its FT is [16] which may be expressed aŝ Note that In this casep (z) and thus η (z) satisfy the Hermitian property. However it has a branch point at z = 0; In additionp Hence the ruin probability is which is difficult to evaluate numerically as it involves the integral of a non-linear function of an integral. Hence we resort to contour integration over R + , indented by a small circle γ in C − − at z = 0 of radius r = 10 −9 , as shown below.
Remark 8 It is tempting to evaluate ψ (u) by finding the zeros of η (z) . However it has possibly thousands of zeros, depending on the tolerance given to 0. [17] also employ complex analysis to derive approximations for the ruin probability. However they lack a simple expression for the FT.

The Pareto Case
The Pareto density for claims 2 may be written simply as where α > 1 and q = α−1 α , so as give a unit mean. The variance is given by The parameter α > 1 defines its shape. For unit variance we must have Its FT is known [9] and is given bŷ where Γ (−α, −iqz) is the upper incomplete gamma function.
Γ (α, z) admits the expansion for α not a non-negative integer.
Unfortunately it is well known that p does not have finite moments of all orders, so thatp cannot be analytically continued to z = 0. In fact, it has an essential singularity there, which makes evaluation of residues problematic where α is not an integer.
Fortunately Γ (−α, −iqz) has a series expansion 3 for α not a positive integer: where M = αΓ (−α) is a constant in a similar manner to the contour integration of the lognormal. To simplify the notation, we use the transformation z = iv, or −iz = v in 15 to givê There are various ways to calculate 15. Direct integration of the FT of 14 is one way. Another is to employ the Matlab function for the upper incomplete gamma, known as igamma(z). However this is not employed in rootfindng, as it requires an inordinate amount of computer memory. The third way to truncate the summation in 15,either version, with or without expansion of the term e −iqz . Fortunately. all these methods lead to the same computational results, as shown in Figure 2.
We take α = 1 + √ 2, ensuring that p has unit mean and variance. As the factorial appearing in the expressions above allows very rapid convergence of the series expansions, we may truncate the summations for k, κ = 1..K, and assume K = 5, say. Its derivative is Thus using truncation of the infinite series; dx.
The roots of the denominator as follows:   is similar to that in [18], which uses exponential integrals. These underpin the upper incomplete gamma function. However their version of the Pareto density has support on x > 0 and has cdf for positive integer m, which is inconsistent with 14.

Results
The ruin probabilities for the various claim densities may be compared in the following chart, with parameters chosen for each density to maximize comparability.
There is evidently uncertainty at u = 0, which is an essential singularity for some claim densities.The results are not altogether surprising; the two extreme case are those of the Pareto and Lévy densities; the tail of the former is greatest, and the variance of the latter greatest.

Conclusions
In this paper we demonstrate how complex function theory enables a coherent approach to the solution of ruin probability problems. This has involved a heavy application of the Cauchy theorem for analytic functions. It provides exact solutions in many cases, or to any desired degree of accuracy for more complicated distributions (such as the Gamma and the Pareto). Further, extensions (or rather contractions) of these results using complex analysis may be made for finite or discrete time ruin probabilities in future papers.
(f) for p ∈ L 2 (R) and z ∈ C + ,p (z) is a holomorphic function. Moreover, by Plancherel's theorem, one has |p(x + iy)| 2 dx ≤ Property (b) is the Riemann-Lebesgue lemma for FTs. Property (e) follows directly from the fact that p is real. The differentiability of p in property (a) follows from taking the derivative as a limit and applying the Lebesgue dominated convergence theorem [19, §l.l(g)]. Property (f) is the Paley Wiener theorem.
To apply the FT to the ruin equations it will be necessary to recover a function from its FT. For the purposes of this paper the forms of the Fourier inversion theorem in Appendix A will suffice: Proposition 11 If either p ∈ L 1 (R) or p ∈ L 2 (R) then the inverse FT ofp exists and, moreover, the above inversion formula holds in the appropriate norm. If p ∈ L 1 (R) has a discontinuity at z, then so that at a point of continuity 1 2π e −iuzp (z) = p (u) . . As a consequence the Proposition implies that the inversion formula holds pointwise if p is continuous. We remark that p ∈ L 2 (R) if and only ifp ∈ L 2 (R), as stated in Proposition l(f). However the corresponding property for L 1 (R) does not generally hold.
In later sections of this paper we consider the function defined as η(z) =p − 1 − iPz for a given p ∈ L 1 (R). This is clearly an analytic function whereverp is defined and analytic. Other relevant properties of this function are summarized below.
(b) Ifp can be analytically continued to a neighborhood of z = 0, then there exists a root of η(z) with z ∈ C − and (z) = 0. In addition this root has the smallest modulus of all roots in C − .
(a) We first demonstrate the proposition for θ = 0. The function η clearly has a root at z = 0. To show that it is unique in C + define the function We have g 1 = 1 andĝ =p −1 iz so that applying Proposition l(e) to the function g:  If θ ∈ C + then the circle γ lies completely within C + , whereas if θ ∈ R then it touches the real axis at x = re θ P .In either case, |p| ≤ 1 for z ∈ C + ∪ R by Proposition l(e), so that η(z) cannot have a zero outside γ.
In the case of θ ∈ C + it is clear that a closed curve Γ ⊆ C + may be constructed surrounding γ, on which holds the inequality: Hence by Rouché's theorem [14, §8,2], the function 1 + iPz − iθ has precisely the same number of zeros within Γ as 1 + iPz − iθ−p. But it is easily shown that the former function has precisely one such zero, from which the result follows for θ ∈ C + . (Note that this also gives the proof where θ ∈ R, but only ifp θ P < 1.) In the case of θ ∈ R we use a continuity argument to establish the existence of a root of η(z) = iθ. Let {iθ n ∈ C + , n = 1, 2, 3 . . .} be a sequence such that θ n → θ.
Then from the previous case, there exist {z n } such that η(z n ) = iθ n . Now the sequence {z n } is bounded and hence must have a limit point z with z ∈ C + ∪ R. If necessary we can construct a convergent subsequence so that z n → z say. Since the function η is continuous, we have η(z) = iθ, which proves existence of a root. To show uniqueness, let ω be another zero, so that we have: .Using the same argument as for the proof of part (a), we consider in place of p(u) the function e iuω p(u) ∈ L 1 (R) and the related function g(u) = ∞ u e ivω p(v)dv ∈ L 1 (R).

B Some Immediate Implications
The above formula forψ in Proposition 12 is by no means new. It has been attributed variously to Khintchine, Lundberg, and other authors [20, eqn 4.44]. Unlike previous versions, the formula above is expressed in terms of the FT rather than the Laplace Transform. This enables some easy implications to be drawn from complex function theory.
The well known relation ψ(0) = 1 P for zero initial reserve has been derived quite simply in §4.2.

Proposition 13
The convolution formula for ψ [4] may be shown from Proposition 12 as follows.
Since P > 1 we have ĝ P < 1 for z ∈ C + ∪ R, so that we may expand 1 −ĝ P as an infinite geometric series, giving: It may be seen, from using the dominated convergence theorem to interchange an integral and an infinite sum, that this is just the FT of the convolution formula where G is the function ∞ u g (v) dv and G n * = G (n−1) * * G, for n > 1.
For u = 0, this becomes Proposition 14 Under the conditions of proposition 12, the asymptotic formula of Lundberg is evident.
We know that η has at least one zero iξ ∈ C − , as in 12. Since iξ is the zero of smallest modulus in C − ,This may be proved by expressing ψ as a (possibly infinite) sum of exponential terms. In the literature it is also known as the adjustment coefficient.