Numerical Ultimate Ruin Probabilities under Interest Force

This work addresses the issue of ruin of an insurer whose portfolio is exposed to insurance risk arising from the classical surplus process. Availability of a positive interest rate in the financial world forces the insurer to invest into a risk free asset. We derive a linear Volterra integral equation of the second kind and apply an order four Block-by-block method in conjuction with the Simpson rule to solve the Volterra equation for ultimate ruin. This probability is arrived at by taking a linear combination of some two solutions to the Volterra integral equation. The several numerical examples given show that our results are excellent and reliable.


INTRODUCTION
The problem of finding the probability of ultimate ruin was first considered by Lundberg [1] . Since then, the problem has received much attention up to present day. In his thesis, Lundberg considered a surplus model of the type: Surplus = Initial reserve + Income -outflow. Among the earlier authors who gave a rigorous mathematical basis of Lundberg's work was Cram r [2,3] . His contributions were presented in his monograph 'Collective Risk Theory'. Lundberg's model, expounded by Cram r, is termed the Cram r-Lundberg model or the surplus model.
In this model, at time t, the surplus t Y of an insurance company is given by where, 0 0 ≥ = Y y is the initial reserve, p > 0 is the premium rate, { } + is a Poisson process with intensity λ , modelling the number of claims in (0, t] is an independent and identically distributed sequence of positive random variables (with distribution F ) independent of N, modelling the claim sizes. The distribution F has finite expectation µ and finite variance 2 σ . In the literature, the process Y in equation (1) is commonly known as the classical risk model. A critical look at the process in (1) raises a couple of questions. One question that has received much attention is 'what is the probability that Y ever becomes negative?' The first time when this happens is termed the time of ruin and the associated probability is the probability of ruin. Ruin is considered as a technical term. It does not mean that the company is bankrupt. However, if ruin occurs, it is interpreted as meaning that the company has to take action in order to make the business profitable.
The Cram r-Lundberg model serves as a skeleton for more realistic models that have been studied in the insurance literature. This standard model for nonlife insurance is simple enough to calculate probabilities of interest, but too simple to be realistic. For example, it does not include interest earned on the invested surplus. There are several papers treating this model in many directions and forms, all with a view of finding the probability of ruin. By far the majority of these papers are concentrated on the analytical aspects of the problem but there is also a quite considerable number that deal with numerical methods to calculate this probability. More on the history of this problem can be traced from Segerdahl [4,5] , Andersen [6] , Davidson [7] , Thorin [8] , Wikstad [9] , Gerber [10] , Harrison [11] , De Vylder [12] . For a general background to ruin theory, we refer to B hlmann [13] . In this study, we shall be concerned with ruin under interest force and our emphasis will be on numerical methods.

THE MODEL AND THEORETICAL RESULTS
All processes and random variables are assumed to be defined on the stochastic basis , t t f f satisfying the usual conditions, i.e. t f is right continuous and P-complete. Here, Ω is an abstract sample space whose elements are denoted as ω ; f is aσ -algebra on Ω ; Ρ is a probability measure and is a filtration. A filtration means an increasing and right continuous family of sub σ - The risk model for the surplus t Y of an insurance company at time t is given by equation (1). A natural generalisation of this model is to allow reserves to earn interests. If there is a positive constant force of interest r, the surplus now become Τ is the time of ruin, equal to infinity if ruin never occurs. Then ) ( y ψ is the probability of ultimate ruin when initial capital is y and it is well known [14] that satisfies the integrodifferential equation Obviously ( ) 1 = y ψ for y < 0 and it is well known [15] Theorem 3.1, that where the kernel and ( ) ( ) If we knew ( )  (14) in Sundt and Teugels [16] , but this is rather complicated to calculate, also numerically. Sundt and Teugels give easy to calculate inequalities for ( ) 0 ψ , but they are less useful if one is interested in exact values for ( ) y ψ . Asmussen and Nielsen [17] provide an efficient simulation procedure for ( ) Simpler simulation procedures are suggested in Paulsen and Rasmussen [18] , but they are less precise than those of Asmussen and Nielsen. They also do not work well when p ≥

λµ
. The method we present in this work works well p ≥ λµ .
We will now show how equation (2) where, ( ) x y K , is as in (3). Note that in particular . It follows by Linz [19] , Theorem 3.2 p.32 and dominated convergence that (4) has a unique continuous solution for all y and all a. The following result will be proved at the end.
. From this it is clear that any solution of (4) can be written in the form (5) for a suitable choice of α . In particular, the choice ( ) hence, We thus have Corollary: Let 2 1 a a ≠ and α be given by (6). Then The procedure to numerically find ( ) By Moiseiwitsch [20] , (9) has a unique integrable solution on any interval (0, T]. Letting y = 0 in (9) but may be infinite. Again by Lemma 1, to prove that the limit is finite, it is sufficient to prove that for some constant C and a function ε which satisfies , hence by Lemma 1, for all y. By the version of a Tauberian theorem given in Linz [19] , Theorem 6.8, (10) follows provided we can prove that Where, 1 c , 2 c and 3 c are constants. Replacing t by a complex z in (12), it is clear that the right hand side of (12) is analytical for Re(z) > 0. Hence, by Widder [21] , exists for all t > 0. Furthermore

NUMERICAL METHODS
Here, we will discuss numerical solutions of (2) (and its associated equation (4)).
To fix ideas, we write a general form of both equations thus Where, K(y, x) is as in (3) where, i g is the numerical approximation to The i ω are the integration weights. Here we shall use the block-by-block method in conjuction with the Simpson integration rule, thus obtaining solutions in blocks of 2. The block-by-block is recommended [22] as the best of the higher order methods for numerically solving (13). To briefly explain the method, Simpson's rule gives Quadratic interpolation gives that 2 8 and inserting this into (16) yields ( ) Equations (15) and (17)  It follows from results by Linz [19] that for fixed y so that y nh = , the solution satisfies ( ) ( ), provided g is four times continuously differentiable as is the case here by Theorem 2.1 in Paulsen et al. [23] . On the other hand, Dickson and Waters [24] presented several numerical methods for this case, but without any effort to analyse the error rate. Of all the methods they discuss, the one that is closest to ours is a method by Sundt and Teugels [16] . They showed that it performs very well, although it is clear from the algorithm of their method that it is of order no higher than 2. The results i show that the block-by-block method performs extremely well.

NUMERICAL RESULTS
We now report some numerical results obtained using the method described earlier.
For a given The implementation of the method we described earlier was done using FORTRAN and taking advantage of the Double Precision feature to get satisfactory accuracy. Of course slower programs like Splus, R, Gauss, Matlab, Maple or Mathematica could have been used, but at the expense of considerably longer computing time.
The process parameters 05 . 0 , 2 , 3 = = = r p λ are used and S is assumed to be exponentially distributed with expectation 2. The true ruin probability was found way back by Segerdahl [5] and later reappearing in many papers. The exact ruin probability is calculated to a high degree of accuracy from Segerdahl's formula using the Splus program integrate. Tables 1, 2 and 3 show that we can achieve the same results irrespective of what a1 we choose. The value of (0) is very accurate! As expected, the smaller the step size h the better the results. This is the case as shown in Table 1-3  compared to the results in Table 4. Note that the stepsize in Table 4 is 10 times that in Table 1-3.

CONCLUSION
We have been able to use a numerical method which is simple and straight forward to calculate ultimate ruin probabilities under interest force. The results indicate the quality of our numerical method. An extension of this problem is to allow a diffusion in (1) thus With another diffusion, the return on investments process will no longer be rt t R = but For a detailed description of such a model, see e.g Paulsen et al. [23] . The same method is likely to work once the corresponding Volterra equation is obtained. The task will be to identify the right kernel and the forcing function to use in the FORTRAN program. Also, apart from the exponential distribution, other distributions can be used. In the event that the distribution is difficult to achieve analytically, Simpson's rule can be employed. Such an approach was used in Paulsen et al. [23] .