Received 12 January 2016; accepted 26 February 2016; published 29 February 2016
An actuarial risk model is concerned with the study of the mathematical aspects observed in the behavior of a collection of risks generated by an insurance portfolio. In general insurance risk modeling, the two quantities of paramount importance are the number of claims arriving in a fixed time period and size of each claim. Modeling of the former aspect is done in terms of a discrete distribution, more specifically a counting distribution whereas the claim severity is modeled through what is known as a loss distribution.
There are compelling reasons to use mathematical models to describe insurance loss amounts (claim severities). As specified, in  in the most general sense, all of actuarial science is about loss distributions. For any insurance company, the sound statistical analysis of the underlying claim scenario, reflected in terms of loss modeling is of uttermost importance because it is the basis on which lies the subsequent determination of various actuarial quantities of interest like probability of ruin, premium loading, pure premium, expected profits, reserved to be maintained and the impact of reinsurance and deductibles.
A good introduction to the subject of fitting distribution to losses is given in  . Most data in general insurance are skewed to the right and therefore distributions with high degree of positive skewness such as Lognormal, Pareto, Gamma, Weibull and Burr had been used by actuaries to fit claim sizes  . However as stated in  , in different classes of insurance business, it is not clear which distributions are suitable for modeling claims arising in different portfolios.
The three-parameter Burr XII distribution was originally used in the analysis of lifetime data and is becoming increasingly useful in the context of actuarial science  . The data used in this paper are on motor insurance where one typical characteristics is the occurrence of large but infrequent claims and hence there is a need to fit and use a statistical distribution which is heavy tailed and highly skewed towards the right and this justifies why heavy tailed distributions such as Pareto, Weibull and Burr are suitable candidates for loss modeling in motor insurance (see   ). Although there can be other suitable models for loss modeling in general insurance, we are primarily concerned with the Burr distribution as a loss model for our claim data and have concentrated on the computation of various actuarial quantities like the probability of ruin and the moments of the time to ruin when the loss model or claim severity model is Burr XII. Literature (   ) reveals that no closed form expression is available for the determination of these quantities in case of Burr XII distributed claim amounts and hence, we resort to numerical techniques to compute them. Briefly our objectives for this paper are:
1) To fit the Burr XII distribution to a set of insurance data through an algorithm mentioned in  and to assess the goodness of fit through some statistics based on the empirical distribution functions (EDF statistics).
2) To compute the probability of ultimate ruin when the claim severity is distributed as our fitted Burr XII distribution as well as a Burr XII distribution with a set of illustrative values for its parameters, using two numerical algorithms namely a stable recursive algorithm and the method of Product Integration.
3) To compute the first two moments of the time to ruin in case of Burr XII distributed claim severity.
The first part of the paper deals with the Watkins  algorithm for obtaining the MLE of the parameters of the Burr XII distribution, followed by testing the goodness of fit through some statistics based on the empirical distribution function (EDF). This is followed by the computation of the probability of ultimate ruin for various values of the initial surplus through the adaptation of the two above mentioned numerical methods. The concluding section deals with the computation of the moments of the time to ruin in case of Burr XII distributed claim severity.
The illustrative Burr XII distribution that is being used is the one that is being fitted to the Property Claim Services (PCS) dataset covering losses resulting from natural catastrophic events in USA that occurred between 1990 and 1999  .
2.1. Fitting of the Burr Distribution
The pdf of the three parameter Burr XII distribution is given by
The algorithm for finding the maximum likelihood estimators (MLE) for the parameters of the Burr XII distribution is taken from  and this algorithm exploits the link between the three parameter Burr XII distribution and the two parameter Weibull distribution with the latter emerging as the limiting case of the former.
The cumulative distribution function for the two parameter Weibull distribution is given by
in which the positive parameters are respectively the shape and the scale parameters.
The basic two parameter Burr XII distribution with shape parameters α and τ has the cumulative distribution function
An scale parameter is introduced into (2.1.3) by substituting thereby giving the cdf of y as
Letting with remaining finite, it is seen that the Burr XII distribution emerges as the limiting
distribution for the Weibull distribution with shape parameter and scale parameter
If we consider a sample of “m” items from the Weibull distribution whose cdf is given by (2.1.2), the log likelihood function is given by
And the log likelihood function of the Burr XII distribution is given by
The main steps of the algorithm are:
Step 1: First, we find the maximum likelihood of the parameters appearing in (2.1.2) using the Multi parameter Newton Raphson Iterative method yielding the two values and.
Step 2: Then, we rescale the original data by so that in implementing the Newton Raphson for determin-
ing the MLEs of the parameters of the Burr XII distribution, the utilized values are the rescaled values.
The argument in  leads us to conclude that rescaling the data introduces a large amount of stability into the algorithm. After the parameter estimates have been obtained, the MLE for the original observations are obtained by undoing the effect of scaling on the estimated values. In Appendix, we have given a very brief introduction to the Multi Parameter Newton Raphson method and have obtained the gradient and hessian matrices for the Weibull and the Burr XII distribution which are required for obtaining the maximum likelihood estimators for the parameters of the Weibull and the Burr XII distributions respectively.
2.2. Classical Risk Model
Let denote the surplus process of an insurer as
where is the initial surplus, c is the rate of premium income per unit time and is the aggregate claim process and we have where is a homogeneous Poisson process with parameter, denotes the amount of the ith claim and is a sequence of iid random variables with
distribution function F such that and probability density function f. We denote by. Also we have, where is the security loading factor.
Let denote the time to ruin from initial surplus u so that and define
and. is known as the ultimate ruin probability
whereas is the finite time ruin probability. For a detailed discussion on the Classical Risk model and the probability of ruin refer to    - .
However it needs to be noted that the Classical Risk Model involves many simplication criterions which might make it deviate from real life situations. For example, the classical risk model assumes that the intensity parameter λ is independent of time, independence between claim severity distribution and claim number distribution, premium is received continuously in time, surplus earns no interest and is neither subjected to tax, no effect of inflation etc. Despite such assumptions, Classical Risk model still constitute the basis of many models in insurance mathematics.
Probability of Ruin is a very important component of the operational risk theory. It reflects the volatility in the business and can serve as a useful tool in long range planning for the use of insurer’s funds. Ruin, in some sense, corresponds to the insolvency of the insurance company although; solvency/insolvency of an insurance company involves many other complicated considerations.
2.3. A Stable Recursive Algorithm for the Evaluation of the Ultimate Ruin Probabilities
Probability of ruin can be obtained as the solution of an integro differential equation  . Stable Recursive Algorithm as the name suggests, is a recursive algorithm which is used to solve the convolution part of the integro differential equation for the probability of ultimate ruin and this in turn leads to the bounds of the ultimate ruin probability within a prescribed tolerance level.
According to  , for calculating the infinite time ruin probability numerically, one has to solve the following integral equation
Let and (2.3.2)
The usual approach is to apply a discretization technique to approximate the integral in (2.3.1) (see  - ). A discretization technique is said to be effective if it does not lead to the propagation of errors thereby producing stable numerical results.
In Reference  , an efficient discretization technique for the convolution integral have been introduced and it is given in the form of a theorem stated therein (Section 2, Equation (5)).
Reference  have used another version of this theorem to obtain analytically upper and lower bounds of the infinite time ruin probabilities in case of constraints in the claim size distributions whereas  have used it to obtain a stable recursive algorithm for deducing the numerical bounds on the infinite time ruin probabilities and as mentioned there, a practical procedure for implementing the stable recursive algorithm is indicated as given below.
1) First carry out the sub division of the interval as given where
“n” the number of intervals is chosen to be sufficiently large.
2) For every let and then for every and for every since is a decreasing function of y.
3) Then the upper bound and the lower bound to the probability of ruin is given by.
Upper bound is
and the lower bound is
with of course,
can be approximated by
An upper bound for the error in the estimation of is given by
The stability of this numerical procedure is justified from the fact that there is no cumulative effect of propa-
gation of error as it can be shown that if are calculated with an error then the error on is bounded by
For the derivation of the bounds and the full justification on the stability of this method, see  .
Computing the Function h(x) for the Burr XII Distribution
We have from (2.3.2),
For the Burr XII distribution given in (2.2.1), we have
Now, it can be shown that,
Using this in Equation (2.3.5), we have,
And this can be computed using the pbeta function of R Software.
Also, we have,
2.4. Product Integration
Product integration which was traditionally used to numerically solve Volterra Integral equation of the second type (see  ) can also be used to compute the ultimate ruin probabilities, especially while dealing with heavy tailed claim severity distributions  .
As stated in Section (2.3), for calculating the infinite time ruin probability numerically, one has to solve the integral Equation (2.3.1) which can also be put in the form (of a Volterra integral equation of the second kind) as shown below
Since, the early 1980’s, the numerical methods for the evaluation of were based on the discretization of the Risk process and then computing recursively using some initial conditions. However, the recursive schemes usually suffer from the drawback of being slow and less accurate because the quadrature rule employed in the recursive schemes are usually of low order. The method of Product integration used to evaluate numerically the Probability of Ultimate Ruin as given in  and as justified there are fast and accurate and are more suitable for dealing with heavy tailed distributions such as Burr XII distribution. As such, in this paper, as a second method, we have used product integration to compute the probability of ultimate ruin as prescribed in  . For a detailed description of this method refer to  , although we have highlighted the execution procedure of this method.
The Volterra integral equation of the second kind is given by
where is the kernel (and is known) and is the unknown function to be determined. If or one of it’s low order derivative is badly behaved in one of its arguments, Newton Cotes integration formulae may produce inaccurate results or converge slowly. To deal with such situations, when is badly behaved,  and  recommend the use of Product Integration.
We first factorize as,
where is smooth and well behaved and can be accurately approximated by a suitable Langrange’s Interpolation Polynomial and is badly behaved.
The interval is divided into n subintervals where
A quadrature rule of the form
where for is used to approximate the integral appearing in (2.4.1). The weights are determined by ensuring that the rule of Equation (2.4.6) is exact when is a polynomial in t of degree.
It is assumed that moments exist for this is necessary for the application of the method of Product Integration and for each i, the moments can be calculated as
Assuming, is linear in t i.e., it can be shown that
And finally, the estimate of is given by, where are obtained recursively by
with and the weights are given by
For accelerating the convergence, as mentioned in  we have used the Richardson’s Extrapolation technique (Also see   ).
2.4.1. Product Integration for the Computation of the Ultimate Probability of Ruin for Burr Distributed Claims
We have used product integration to compute the Ultimate Probability of Ruin for Burr XII distributed claims taking an illustrative value of as
We have considered the Burr distribution which has been fitted to our data as well as a Burr distribution with a set of illustrative values for its parameters.
As derived in (2.3.7),
From (2.4.3), we have
As for this distribution, all the moments exist for any finite s, Product Integration can be used.
2.4.2. Computation of the Weights When the Claim Severity Distribution Is Burr XII
And (w is any upper limit of the integral, w > 0,) can be computed using the
pbeta function of the R software.
Similarly, from (2.4.13), we have
2.5. The Moment of the Time to Ruin
The distribution of the time to ruin is an interesting quantity related to the probability of ruin and plays a vital role in warning the management for possible adverse situations. There is no closed form expression developed for the distribution of the time to ruin except for the Exponential distribution and the Erlang group of distributions  and hence the computation of the moments except for the mentioned distributions has to be done numerically.
Initial ideas on this aspect can be found in  and working on those ideas,  presented methods from which explicit solutions for the moments of the time to ruin can be found recursively provided that an explicit solution exist for the ultimate ruin probability. Reference  simplified the results of  to make them mathematically tractable for numerical computation and have used them to calculate the approximate values for the moments of the time to ruin when explicit solutions for the probability of ultimate ruin do not exist. In their numerical computations, values of have been calculated from the stable algorithms described in  . (Also see   .)
Reference  showed that the moment of the distribution of the time to ruin T is given by
Let L: the maximum of the aggregate loss process so that see,  , formula (12.6.2).
In  , it had been shown that
Formula (6.2.1) of  has been simplified in  as
Hence, can be evaluated using numerical integration.
Similarly, appearing in  has been simplified in  as
3. Results and Discussion
Data: Our data is a set of 160,000 claim amounts spread over a period of 6 months i.e. from April, 2013 to September, 2013 obtained from Bajaj Allianz General Insurance company, India from its motor insurance portfolio covering all its branches in India. No adjustment was made for inflation for the time horizon is narrow. It needs to be mentioned that the data is utilized more for the illustration of the various methodologies rather than for the extraction of any concrete meaningful conclusion. Since the inter arrival time of claim was difficult to track, the intensity parameter was estimated on the basis of the number of claims arriving per day during the period.
Summary statistics of the data as shown in Table 1 reveal the existence of high coefficient of skewness which suggests that a highly skewed right tailed distribution such as the Burr XII can be a probable candidate for modeling this data. The histogram of the data plotted in Figure 1 and the empirical probability density function plotted in Figure 2 show the same trend. These figures too indicate a high degree of skewness towards the right which in a way justifies the use of Burr XII distribution for modeling our data.
For finding the maximum likelihood estimators for the parameters of the Burr XII distribution, the use of the algorithm mentioned in  has been made. The log likelihood got maximized at the 30th iteration thereby giving the estimated values of the parameters as shown in Table 2. Initial assessment of the fit was done through some
Table 1. Summary statistics for the Insurance claim data.
Table 2. Parameter estimates for the Burr XII distribution obtained through the Watkin algorithm and the value of the EDF statistics along with their p-values indicated in parentheses.
Figure 1. Histogram of the observed claim data on motor insurance.
Figure 2. Estimate of the probability density function for the claim data on motor insurance.
graphical displays. Figure 3 show the histogram for a set of data simulated from the Burr XII distribution with the values of the parameters as estimated using the algorithm. This histogram has some resemblance with the histogram for the observed data as shown in Figure 1. The QQ plot displayed in Figure 4 indicate moderate deviation from the straight line passing through the origin which leads us to conclude that the fit is moderately good. The assessment of the fit was also done through the computation of two statistics namely the Anderson Darling statistics and the Cramer von statistics which are based on the empirical distribution function  . Table 2 shows the values of the Anderson Darling and Cramer Von statistics for testing the goodness of fit along with their p-values which were obtained through the Monte-Carlo simulation based on 100 iterations  .
Hence, we have little evidence to believe that the Burr XII distribution adequately describes the claim data. However, in the subsequent sections, we have used this fitted Burr distribution along with another Burr XII distribution with a set of illustrative values for its parameters mainly with the objective of depicting the computational methodologies associated with the Burr XII distribution in obtaining some of the important Actuarial Quantities.
Table 3 and Table 4 display the upper and lower bounds to the probability of ultimate ruin in case of Burr XII distributed claim severity using the stable Recursive algorithm. In Table 3, the claim severity is the fitted Burr XII distribution whereas the Table 4 is constructed with an illustrative Burr XII distribution. The number of intervals for the stable recursive algorithms has been taken to be n = 160 and an illustrative value for the security loading factor has been taken as θ = 0.3. The upper bounds to the error of estimation have also been indicated in both the tables.
In both the tables, it has been observed that the probability of ultimate ruin is decreasing with an increase in the initial capital which is as expected. In case of our fitted Burr XII distribution, the difference between the upper
Figure 3. Histogram for a data set simulated from for Burr XII distribution with, and.
Figure 4. QQ Plot between the empirical quantiles estimated from the motor insurance data and the theoretical quantiles for Burr XII distribution with, and.
Table 3. Upper and lower bounds to the ultimate ruin probabilities for Burr distribution with, and computed through the stable recursive algorithm.
Table 4. Upper and lower bounds to the ultimate ruin probabilities for Burr distribution with, and computed through the stable recursive algorithm.
and the lower bounds to the probability of ultimate ruin seems to be decreasing in the absolute sense which has lead to the decline in the upper bound to the error of estimation. In case of the Burr XII distribution with a set of illustrative values for its parameters (Table 4), which we would choose to call the illustrative Burr XII in the subsequent sections, it is observed that there is no visible difference between the upper and lower bounds to the probability of ultimate ruin and hence the algorithm seems to be giving more accurate results in this case. The probable reasons for the equality of both the bounds can be that for this set of illustrative values of the parameters for the Burr XII distribution, the function h(x) is getting stabilized (approaching p1) more rapidly.
Table 5 shows the probability of ultimate ruin for the fitted Burr XII distribution obtained through the method of product integration whereas the Table 6 displays the corresponding values for the illustrative Burr XII distribution. For accelerating the convergence, we have used the Richardson’s extrapolation technique as mentioned in  with and thereby giving. Both the tables reveal that the values of the ultimate ruin probabilities obtained through the method of product integration are fairly consistent with those obtained through the stable recursive algorithm. One notable trend observed in the computation of the probability
Table 5. Ultimate ruin probabilities for Burr XII distribution with, and obtained through product integration.
Table 6. Ultimate ruin probabilities for Burr XII distribution with, and obtained through product integration.
of ultimate ruin through both the algorithms is that in case of our fitted Burr XII distribution with an increase in the value of the initial capital, the values for the probability of ultimate ruin are decreasing at a significantly high rate whereas this declined at a moderate rate in case of the illustrative Burr XII distribution.
In computing the moments of the time to ruin, in contrast to  , where has been obtained through the stable algorithms described in  , in our case, has been computed through the stable recursive algorithm mentioned above. For the illustrative Burr distribution, we have also computed the second moment for a few values of the initial capital since the computing time for the evaluation of the second moment is significantly high. The high execution time is attributed to the fact that the computation of the second moment has taken
the first moment as an input. For numerical integration, we have used Simpson’s rule for numerical in-
tegration. All of the computations have been done using the R software  .
From Table 7 showing the first moment (mean) of the time to ruin for the fitted Burr XII distribution as well as for the illustrative Burr XII, we make the following observations;
1) Mean (in years) of the time to ruin for the illustrative Burr XII i.e. Burr XII with, and as indicated in column (3) for the various values of u is consistent with practical rationalism as the value goes on increasing with an increase in the initial capital for it is expected that induction of more capital should prolong the occurrence of ruin (if any).
2) Mean (in years) of the time to ruin for the fitted Burr XII distribution i.e. Burr XII with, and is not consistent in the above sense (column (2)). Probable reasons for this inconsistency can be the occurrence of error in the evaluation of and also in numerical integration. As seen in Equation (2.5.4), numerator of the expression for the mean of the time to ruin is the difference of two quantities and it is the latter part which is evaluated numerically and although not shown explicitly, both the parts were of very low order (of the order 1e−04) and were highly affected by the occurrence of numerical error leading to inconsistent results. One important point to be noted is that compared to the illustrative Burr XII distribution, the fitted Burr XII distribution has very low second order moment () and the second order moment of the underlying claim severity distribution is an important factor in determining the mean of the time to ruin, in fact the latter exists if and only if the second order moment of the claim severity distribution exists  . This low value of the second order moment for the fitted Burr XII distribution might have been the cause for this inconsistency in the pattern of the mean of the time to ruin.
Table 8 displays the second order moments of the time to Ruin for our illustrative Burr XII distribution for a few values of the initial capital.
In obtaining the mean of the time to ruin, an illustrative value of λ has been taken as λ = 32.78 and the same value is retained in obtaining the second moment of the time to ruin.
Considering the fact that heavy tailed right skewed distribution like Burr XII arises frequently in case of risk modeling in general insurance, our work may be useful for insurance practitioners and experts from the financial industry. Our main objective was to compute the probability of ultimate ruin in case the claim severity is distributed as Burr XII distribution and this was implemented through the application of two numerical algorithms.
Table 7. First moment (mean) of the time to ruin in case of Burr XII distributed claim severity distribution.
Table 8. Second moment of the time to ruin in case of Burr distributed claim severity distribution.
However, no effort has been made to judge which method gives the better estimate to the probability of ultimate ruin as there is no exact expression for it in case of Burr XII claim severity distribution. Although, there are instances in literature  , where the probability of ultimate ruin obtained through the Pollachez Khinchin formula has been used as a baseline method and the accuracies of other methods judged in terms of it, yet we have not gone to the extent of making this comparison between the two methods. The work reveals fairly good amount of consistencies in the approximate values of the probability of ultimate ruin obtained by the two numerical algorithms under consideration.
Also, in obtaining the moments of the time to ruin, it was found that in case of the illustrative Burr, the first moment (mean) of the time to ruin is increasing with an increase in the value of the initial surplus. This is to be expected in practice, because the induction of larger surpluses tends to prolong the time to ruin (if it ever happens). However in case of our fitted Burr XII distribution, there was deviation from this intuitive logic, implying that the mean of the time to ruin was found to be decreasing with an increase in initial surplus. The numerical error accumulated via the two numerical algorithms namely the numerical computation of the value of ultimate ruin probability through the stable recursive algorithm and then inserting it as an input into another numerical algorithm to compute the mean of the time to ruin might be the cause of this deviation. The executing time for computing the second moment was too high thereby limiting us just to the computation of this moment for a very few values of the initial surplus.
Extension of this work can be directed towards the computation of other actuarial quantities like aggregate claim models, number of claims until ruin etc in case of Burr XII claim severity. Further analysis is required to give more explicit error bounds to the solutions generated via the two numerical algorithms and the control of error in the numerical computation of the moments of the time to ruin.
The authors would like to thank the anonymous referees of the Journal for imparting valuable suggestions that helped us to improve the paper.
A.1. The Newton Raphson Method: The Multiparameter Situation
One of the most used methods for optimization in the Multi Parameter situation in Statistics is the Newton- Raphson method which is described briefly as given below:
Assume is a vector of p (say) unknown parameters and the log likelihood of the distri-
bution involving is given by. Then the MLE for are obtained by solving the equations
Let us now define what is known as the gradient matrix and the Hessian matrix given by.
The gradient matrix is given by
And the Hessian matrix is given by
Then the iterative relationship for the multi parameter Newton Raphson method is given by
where is the estimated value of at the iteration. The iteration is carried out until there is no significant difference between and
A.2. Multi Parameter Newton Raphson for Weibull Distribution
The log likelihood of the Weibull distribution is given by (2.1.5).
The Gradient matrix for Weibull is given by
and the Hessian matrix is given by
A.3. Multiparameter Newton Raphson for Burr XII Distribution
The Log likelihood of Burr XII distribution is given by (2.1.6).
Its Gradient matrix is given by
The hessian matrix is given by
Tests Based on Empirical Distribution Function
A statistics measuring the difference between the Empirical and the fitted distribution function is called an EDF (empirical distribution function) Statistics and is based on the vertical distances between the distributions.
A class of measures of discrepancy given by the Cramer-Von Mises Family is
where is a suitable function which gives weights to the squared differences When
, we obtain the statistic of Cramer Von Mises.
When, we have the statistic of Anderson and Darling. Here n is the sample
size  .
 Hogg, R. and Klugman, S.A. (1984) Loss Distributions. Wiley, New York. http://dx.doi.org/10.1002/9780470316634
 Klugman, S.A. (1986) Loss Distributions. Proceedings of Symposia in Applied Mathematics, Actuarial Mathematics, 35, 31-55. http://dx.doi.org/10.1090/psapm/035/9873
 Embrechts, P., Klüppelberg, C. and Mikosh, T. (1997) Modelling Extremal Events for Insurance and Finance. Springer, Berlin. http://dx.doi.org/10.1007/978-3-642-33483-2
 Grandell, J. (1991) Aspects of Risk Theory. Springer, New York. http://dx.doi.org/10.1007/978-1-4613-9058-9
 Watkins, J.A. (1999) An Algorithm for Maximum Likelihood Estimation in the Three Parameter Burr XII Distribution. Computational Statistics and Data Analysis, 32, 19-27. http://dx.doi.org/10.1016/S0167-9473(99)00024-9
 Burnecki, K., Misiorek, A. and Weron, R. (2005) Ruin Probabilities in Finite and Infinite Times, Statistical Tools for Finance and Insurance. Springer, Berlin Heidelberg. http://dx.doi.org/10.1007/3-540-27395-6_15
 Gerber, H.U. (1982) On the Numerical Evaluation of the Distribution of Aggregate Claims and Its Stop Loss Premiums. Insurance: Mathematics and Economics, 1, 13-18. http://dx.doi.org/10.1016/0167-6687(82)90016-6
 Panjer, H.H. (1986) Direct Calculations of Ruin Probabilities. Journal of Risk and Insurance, 53, 521-529.http://dx.doi.org/10.2307/252402
 Goovaerts, M. and De Vylder, F. (1983) Upper and Lower Bounds to Infinite Time Ruin Probabilities in Case of Constraints on Claim Size Distributions. Journal of Econometrics, 23, 77-90.http://dx.doi.org/10.1016/0304-4076(83)90076-3
 Ramsay, C.M. and Usabel, M.A. (1997) Calculating Ruin Probabilities via Product Integration. ASTIN Bulletin, 27, 263-271. http://dx.doi.org/10.2143/AST.27.2.542051
 Linz, P. (1985) Analytical and Numerical Methods for Volterra Equations. SIAM Studies in Applied Mathematics, SIAM, Philadelphia. http://dx.doi.org/10.1137/1.9781611970852
 Ramsay, C.M. (1992) Improving Goovaert’s and De Vylder’s Stable Recursive Algorithm. ASTIN Bulletin, 22, 51-59.http://dx.doi.org/10.2143/AST.22.1.2005126
 Lin, X.S. and Willmot, G.E. (1999) Analysis of a Defective Renewal Equation Arising in Ruin Theory. Insurance: Mathematics and Economics, 25, 63-84. http://dx.doi.org/10.1016/S0167-6687(99)00026-8
 Gerber, H.U. and Shiu, E.S.W. (1998) On the Time Value of Ruin. North American Actuarial Journal, 2, 48-78.http://dx.doi.org/10.1080/10920277.1998.10595671
 Dickson, D.C.M. and Waters, H.R. (2002) The Distribution of the Time to Ruin in the Classical Risk Model. ASTIN Bulletin, 32, 299-313. http://dx.doi.org/10.2143/AST.32.2.103
 Lin, X.S. and Willmot, G.E. (2000) The Moments of the Time to Ruin, the Surplus before Ruin and the Deficit at Ruin. Insurance: Mathematics and Economics, 27, 19-44. http://dx.doi.org/10.1016/S0167-6687(00)00038-X
 Dickson, D.C.M., Egideo dos Reis, A.D. and Waters, H.R. (1995) Some Stable Algorithms in Ruin Theory and Their Applications. ASTIN Bulletin, 25, 153-175. http://dx.doi.org/10.2143/AST.25.2.563245
 Dickson, D.C.M. (1989) Recursive Calculations of the Probability and Severity of Ruin. Insurance: Mathematics and Economics, 8, 145-148. http://dx.doi.org/10.1016/0167-6687(89)90028-0
 Dickson, D.C.M. and Waters, H.R. (1991) Recursive Calculation of Survival Probabilities. ASTIN Bulletin, 21, 199-221. http://dx.doi.org/10.2143/AST.21.2.2005364
 R Core Team (2013) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna. http://www.R-project.org/