Computation of the probability of ultimate ruin and some other actuarial quantities under the classical risk model via Fast Fourier Transform

Abstract

The Probability of ultimate ruin under the classical risk model is obtained as a solution of an integro -differential equation involving convolutions and we have used Fast Fourier Transform (FFT) to obtain the approximate values of the probability of ultimate ruin from this integro -differential equation under the situation when the claim severity is modelled by the Mixture of 3 Exponentials and the Weibull distribution. Another application of FFT in ruin theory is shown by means of applying it to obtain the quantiles of the aggregate claim distribution under these claim severity distributions. Extension of the application of FFT is shown by using it to obtain the first moment of the time to ruin under the classical risk model for these distributions. The distributions which have been used are such that one is light tailed and the another is heavy tailed so that a comparison can be made between them on the precision of the actuarial quantities obtained through FFT. FFT has been found to be efficient in obtaining these actuarial quantities when used in conjunction with certain modifications like exponential tilting to control the aliasing error.

Keywords

Fast Fourier Transform aggregate claim modelling time to ruin

1. Introduction

The use of Fast Fourier Transform (FFT) for the evaluation of the probability of ruin is discussed in Embrechts, Grubel and Pitts (1993). Methods of reducing the error that creeps into the use of FFT have been discussed in Grubel and Hermesmeir (1999, 2000). Lima et al. (2002) also embodies a detailed discussion on the use of Fourier transform for the evaluation of the probability of ultimate ruin. The naive application of FFT in computing the probability of ultimate ruin, in the evaluation of the quantiles of the aggregate loss distribution and in computing the first moment of the time to ruin is illustrated in this paper and the methodologies associated with the technical background of the work are cited sequentially as given below.

2. The classical risk model

Let $\left\{{U\left(t\right)}\right\}_{t\geqslant 0}$ denote the surplus process of an insurer as

$\displaystyle U\left(t\right)=u+ct-S\left(t\right)$ (1)

where $u\geqslant 0$ is the initial surplus, $c$ is the rate of premium income per unit time and $\left\{{S\left(t\right)}\right\}_{t\geqslant 0}$ is the aggregate claim process and we have $S\left(t\right)=\sum_{i=1}^{M\left(t\right)}X_{i},$ where $\left\{{M\left(t\right)}\right\}_{t\geqslant 0}$ is a homogeneous Poisson process with parameter $\lambda$ , $X_{i}$ denotes the amount of the ith claim and $\left\{{X_{i}}\right\}_{i=1}^{\infty}$ is a sequence of iid random variables with distribution function $F$ such that $F\left(0\right)=$ 0 and probability density function $f .$ We denote $E\left({X_{1}^{k}}\right)$ by $p_{k}$ . Also, we have $c=\left({1+\theta_{1}}\right)\lambda p_{1}$ , where $\theta_{1}$ is the security loading factor.

Let $T_{u}$ denote the time to ruin from initial surplus $u$ so that $T_{u}=\text{inf}\{t:U\left(t\right)<0\}$ and define $\psi\left(u\right)=$ pr $\{T_{u}<\infty\}=1-\delta\left(u\right)$ and $\psi\left({u,t}\right)=$ Pr ( $T_{u}\leqslant t).\psi\left(u\right)$ is known as the ultimate ruin probability whereas $\psi\left({u,t}\right)$ is the finite time ruin probability. For a detailed discussion on the classical risk model and the probability of ruin, one can refer to Grandell (1991), Panjer and Willmot (1992), Klugman et al. (1998) and Asmussen (2000).

3. The fast fourier transform in ruin theory

The probability of ultimate ruin satisfies the following integro-differential equation (Klugman et al., 1998)

$\displaystyle\psi^{/}\left(u\right)=\frac{\lambda}{c}\psi\left(u\right)-\frac{% \lambda}{c}\int_{0}^{u}\psi\left({u-x}\right)dF\left(x\right)-\frac{\lambda}{c% }\left\{{1-F\left(u\right)}\right\},u\geqslant 0$ (2)

The following is a description on the use of FFT to compute probability of ultimate ruin following the ideas as given in Pitts (2006).

The above equation can also be put in the form of a defective renewal equation (see Lin and Willmot, 1999) as given by

$\displaystyle\psi\left(u\right)=\frac{1}{1+\theta_{1}}\int_{0}^{u}\psi\left({u% -x}\right)dF_{X,e}\left(x\right)+\frac{1}{1+\theta_{1}}\bar{F}_{X,e}(u)$ (3)

and $F_{x,e}\left(u\right)=1-\bar{F}_{x,e}(u)$ is the claim size equilibrium distribution given by

$\displaystyle F_{X,e}\left(x\right)=\frac{1}{p_{1}}\int_{0}^{x}\left\{{1-F% \left(t\right)}\right\}dt$ (4)

It has been shown in Lin and Willmot, (1999), that a solution to Eq. (3) is given by

$\displaystyle\psi\left(u\right)=1-G\left(u\right),G=\sum_{k=0}^{\infty}\left({% 1-\phi}\right)\phi^{k}F_{X,e}^{\ast\left(k\right)}$ (5)

where $\phi=\frac{1}{\left({1+\theta_{1}}\right)}$ .

If $T$ is the time to ruin, then

$\displaystyle E\left\{{T|I(T<\infty}\right\}=\frac{\psi_{1}\left(u\right)}{% \psi\left(u\right)}$ (6)

where $\psi_{1}\left(u\right)$ ,as derived in Eq. (6.22), under corollary (6.1) of (Lin & Willmot, 2000) satisfies the following defective renewal equation.

$\displaystyle\psi_{1}\left(u\right)=\frac{1}{1+\theta_{1}}\int_{0}^{u}\psi_{1}% \left({u-x}\right)dF_{X,e}\left(x\right)+\frac{1}{c}　\int_{u}^{\infty}\psi% \left(x\right)dx$ (7)

If we look into Eq. (3) and Eq. (7), we see that both involve convolutions, viz in Eq. (3), it is of the form $\int_{0}^{u}\psi\left({u-x}\right)dF_{X,e}\left(x\right)$ and in Eq. (7), it is of the form $\int_{0}^{u}\psi_{1}\left({u-x}\right)dF_{X,e}\left(x\right)$ .

This renders us valid scope for the use of FFT for the calculation of the required quantities.

4. Fast fourier transform in applied probability

For a sequence $f_{0},f_{1},\ldots.f_{M-1}$ , the discrete Fourier transformation is defined as

$\displaystyle\chi_{k}=\sum_{m=0}^{M-1}f_{m}\exp\left({\frac{2\pi m}{M}ik}% \right),k=0,1,\ldots M-1$ (8)

and the original sequence $f_{k}$ can be recovered from $\chi_{k}$ by the inverse transformation

$\displaystyle f_{k}=\frac{1}{M}\sum_{m=0}^{M-1}\chi_{m}\exp\left({-\frac{2\pi m% }{M}ik}\right),k=0,1,2,\ldots sM-1$ (9)

Here $M$ is some truncation point.

5. Discretization of continuous distributions

Typically, severity distributions are continuous and since FFT allows only discrete severities as input, discretization of the continuous severity is necessary. To concentrate severity, whose continuous distribution is $F\left(x\right)$ on $\left\{{0,\delta,2\delta,\ldots}\right\},\delta>0$ , the central difference approximation might be used. The central difference approximation is given by

$\displaystyle f_{0}=F\left({\frac{\delta}{2}}\right)$ (10) $\displaystyle f_{n}=F\left({n\delta+\left({\frac{\delta}{2}}\right)}\right)-F% \left({n\delta-\left({\frac{\delta}{2}}\right)}\right)$ (11)

where $\delta$ is a small positive discretization parameter. Likewise, the other modes of discretization are the forward difference discretization given by

$\displaystyle f_{n}^{u}=F\left({n\delta+\delta}\right)-F\left({n\delta}\right)$ (12)

and the backward difference discretization given by

$\displaystyle f_{n}^{L}=F\left({n\delta}\right)-F\left({n\delta-\delta}\right)% ,n=1,2,3\ldots$ (13)

The implication of the discretized distribution is that it assigns to the non – negative integer $`n^{\prime}$ , a probability mass equal to that assigned by $F$ to the interval $\left(\left(n-\left(\frac{1}{2}\right)\right)\delta,\left({n+\left({\frac{1}{2% }}\right)}\right)\delta\right]$ .

Proper choice of the truncation parameter $M$ is to be made so that $\left({f_{0},f_{1},\ldots f_{M-1}}\right)$ can be used as the input array of the FFT. For the effective use of the FFT, we choose $M$ to be a power of 2. To reduce both the discretization and the aliasing error, we need to choose $M$ and $\delta$ such that $\bar{F}((M-0.5)\delta)$ is negligible. (Grubel and Hermesmeir (1999, 2000)).

We now require taking Fourier transform on both sides of $G=\sum_{k=0}^{\infty}\left({1-\phi}\right)\phi^{k}F_{X,e}^{\ast\left(k\right)}$ .

Note that in accordance with Eq. (1.25) of Das (2017), G is the distribution function of the maximal aggregate loss random variable given by $L=Y_{1}+Y_{2}+\ldots.+Y_{K}$ , where $K$ has a geometric distribution with parameter $\phi=\frac{1}{1+\theta_{1}}$ and each of $Y_{i}(i=1,2,\ldots K)$ has the distribution whose cdf is given by $F_{X,e}$ .

6. Computing the characteristic function of the maximal aggregate loss random variable

To compute the FFT of $L$ which is the discrete counterpart of the characteristic function, we proceed as follows. The characteristic function of $L$ is given by

$\displaystyle\chi_{L}\left(t\right)=E\left({e^{itL}}\right)=E\left\{{E\left({e% ^{itL}|K}\right)}\right\}.\qquad\text{(Where $K$ has a compound geometric % distribution with parameter $\phi=\frac{1}{1+\theta_{1}}.)$}$

Therefore, we have

$\displaystyle\chi_{L}\left(t\right)=E\left\{{\left({\chi_{Y}\left(t\right)}% \right)^{K}}\right\}.$

Furthermore, since the characteristic function of the sum of $K$ independent random variables is equal to the product of the characteristic function of each of the random variables. we have

$\displaystyle\chi_{L}\left(t\right)=E\left({e^{-i^{2}K\text{log}\left({\chi_{Y% }\left(t\right)}\right)}}\right)=E\left({e^{iK\left({-i\text{log}(\chi_{Y}% \left(t\right)}\right)}}\right)=\chi_{K}\left({-i\text{log}\left({\chi_{Y}% \left(t\right)}\right)}\right)$ (14)

But $K$ has a geometric distribution with parameter $\phi=\frac{1}{1+\theta_{1}}$ and consequently,

$\displaystyle\chi_{K}\left(t\right)=\frac{1-\phi}{1-\phi e^{it}}$ (15)

Therefore,

$\displaystyle\chi_{L}\left(t\right)=\chi_{K}\left({-i\text{log}\left({\chi_{Y}% \left(t\right)}\right)}\right)=\frac{1-\phi}{1-\phi\chi_{Y}\left(t\right)}$ (16)

Hence, to compute the FFT of $G$ , we first need to discretize the equilibrium distribution as

$\displaystyle f_{e,n}=F_{X,e}\left({n\delta+\left({\frac{\delta}{2}}\right)}% \right)-F_{X,e}\left({n\delta-\left({\frac{\delta}{2}}\right)}\right)$ (17)

To calculate an approximation to the compound Geometric distribution $G$ , we first apply the FFT to the array $\left\{{f_{e,0},f_{e,1},f_{e,2},\ldots,f_{e,M-1}}\right\}$ to obtain the array $\textit{fft}\left({f_{e}}\right)$ .

Then we calculate the array

$\displaystyle g_{\textit{fft}}=\frac{1-\phi}{1-\phi\textit{fft}\left({f_{e}}% \right)}$ (18)

The operations are to be carried point wise. The inverse FFT is then applied to give an array $g=\left\{{g_{0},g_{1},g_{2},\ldots g_{M-1}}\right\}$ , where $g_{n}$ is an approximation to the mass assigned by $G$ to the interval $\left(\left({n-\left({\frac{1}{2}}\right)}\right)\delta\right.$ , $\left.\left({n+\left({\frac{1}{2}}\right)}\right)\delta\right]$ .

Hence an approximation to $\psi\left({\left({j+0.5}\right)h}\right)$ is given by

$\displaystyle\psi_{j}=1-\sum_{i=0}^{j}g_{i}$ (19)

7. Computing the equilibrium distribution for Weibull and mixture of 3 exponentials

For Weibull distribution, the equilibrium distribution is computed as follows.

We have

$\displaystyle F_{X,e}\left(x\right)=\frac{1}{p_{1}}\int_{0}^{x}\left\{{1-F% \left(y\right)}\right\}dy$ (20)

Here $F\left(x\right)=1-e^{-\left({\frac{x}{\theta}}\right)^{\beta}},x>0,\theta>0,% \beta>0\ p_{1}=\theta{\Gamma}1+\frac{1}{\beta}$ .

Therefore

$\displaystyle\int_{0}^{x}\left\{{1-F\left(y\right)}\right\}dy=\int_{0}^{x}e^{-% \left({\frac{y}{\theta}}\right)^{\beta}}dy=\frac{\theta}{\beta}\int_{0}^{x_{1}% }e^{-z}z^{\frac{1}{\beta}-1}dz,x_{1}=\left({\frac{x}{\theta}}\right)^{\beta}.$ (21) $\displaystyle\quad=\frac{\theta}{\beta}{\Gamma}\frac{1}{\beta}\textit{pgamma}% \left({x_{1},\frac{1}{\beta},1}\right),x_{1}=\left({\frac{x}{\theta}}\right)^{\beta}$

where $\textit{pgamma}\left({x,a,s}\right)$ in R is given by

$\displaystyle\textit{pgamma}\left({x,a,s}\right)=\frac{1}{s^{a}{\Gamma}a}\int_% {0}^{x}y^{a-1}e^{-\frac{1}{s}y}dy$ (22)

For mixture of 3 Exponentials, the equilibrium distribution is computed as follows.

Here

$\displaystyle F\left(x\right)=1-\left({w_{1}e^{-\lambda_{1}x}+w_{2}e^{-\lambda% _{2}x}+w_{3}e^{-\lambda_{3}x}}\right),x>0,w_{i}>0,\lambda_{i}>0\left({i=1,2,3}% \right)\sum_{i=1}^{3}w_{i}=1.$ $\displaystyle p_{1}=\frac{w_{1}}{\lambda_{1}}+\frac{w_{2}}{\lambda_{2}}+\frac{% w_{3}}{\lambda_{3}}.$

Therefore,

$\displaystyle\int_{0}^{x}\left\{{1-F\left(y\right)}\right\}dy=\frac{w_{1}}{% \lambda_{1}}\left({1-e^{-\lambda_{1}x}}\right)+\frac{w_{2}}{\lambda_{2}}\left(% {1-e^{-\lambda_{2}x}}\right)+\frac{w_{3}}{\lambda_{3}}\left({1-e^{-\lambda_{3}% x}}\right)$ (23)

8. Introduction to the aggregate loss model

One of the classical problems in risk theory is the computation of the aggregate loss distribution. For the distributions typically used in loss modeling, the closed form expressions for the aggregate claims are not available. Initially, to evaluate the aggregate loss distribution, a number of approximation methods were in use but with the availability of the modern computer processing power, numerical methods like the recursions or the numerical inversion of the Fourier transform are becoming more important and provide excellent results.

The model for the aggregate loss is given by

$\displaystyle z=X_{1}+X_{2}+\ldots+X_{N}$ (24)

where the frequency $N$ is a discrete random variable and $X_{i}^{\prime}s$ are continuous random variables representing the claim severities. Chernoboi et al. (2007), Shevchenko (2010, a) and Shevchenko (2010, b) are some good reviews on this topic. Among the various methods in use for determining the aggregate claim models, one can mention Monte Carlo methods, Panjer recursion and Fourier inversion techniques as most widely used methods. Monte Carlo methods are rather slow whereas the Panjer recursions and Fast Fourier Transforms are better alternatives, although there exists evidence for the existence of many pitfalls in them when they are applied to compute high quantiles of the compound distribution, specially, when the underlying distribution is a heavy tailed distribution (Panjer & Wang (1993), Panjer and Willmot (1986) and Robertson (1992).

The notations and the assumptions of the classical risk model are used in defining the technical background required for introducing the algorithms for the evaluation of the aggregate loss model. The distribution function and the density function of the aggregate loss random variable $Z$ are denoted as ${\cal H}\left(z\right)$ and $h\left(z\right)$ respectively.

The algorithm for implementing FFT to compute the aggregate loss distribution is stated as given below.

(1)

Discretize severity to obtain $f_{0},f_{1},\ldots f_{M-1}$ where $M=2^{r}$ with integer $r$ and $M$ is some truncation point.

(2)

Using FFT, calculate the characteristic function of the severity $\chi_{0},\chi_{1},\ldots\chi_{M-1}$ .

(3)

Using $\xi\left(t\right)=\sum_{k=0}^{\infty}(\chi\left(t\right))^{k}b_{k}=\gamma\left% ({\chi\left(t\right)}\right)$ , calculate the FFT of the compound distribution given by $\xi_{m}=\gamma(\chi_{m})$ ; $b_{k}=P\left({N=k}\right)$ .

(4)

Now perform inverse FFT to $\xi_{0},\xi_{1},\xi_{2}\ldots\xi_{M-1}$ for obtaining the compound distribution $h_{0},h_{1},h_{2},\ldots h_{m-1}$ .

9. Aliasing error and tilting associated with fast fourier transform

If $\sum_{m=1}^{M}f_{m}=1$ , then FFT calculates the compound distribution on $m=0,1,2,\ldots,M-1$ , that is, the mass of the compound distribution beyond $M$ is wrapped and is thus ignored. However, in case of heavy tailed distributions, this practice can lead to large error. This is called the aliasing error. One way to reduce this error is to apply some transformation to increase the tail decay (the so-called tilting).

The exponential tilting for reducing the aliasing error in context of calculating compound distribution was first investigated by Grubel and Hermesmeier (1999), following which and many subsequent research, the following tilting transformation was suggested

$\displaystyle\tilde{f}_{j}=\exp(-j\omega)f_{j},j=0,1,2,\ldots,M-1\ \text{Where% }\ \omega>0$ (25)

One advantage of using this transformation is that this transformation commutes with convolution in a sense that the convolution of two functions $f\left(x\right)$ and $g\left(x\right)$ equals the convolution of their transformation functions $\tilde{f}(x)=f(x)\exp(-\omega x)$ and $\tilde{g}(x)=g(x)\exp(-\omega x)$ multiplied by $e^{\omega x}$ , i.e., $\left({f\ast g}\right)\left(x\right)=e^{\omega x}\left(\tilde{f}*\tilde{g}% \right)(x)$ .

The following is the algorithm for computing the compound distribution via FFT with tilting.

(1)

Define $f_{0},f_{1},\ldots f_{M-1}$ for some large $M$ .

(2)

Perform tilting i.e., calculate the transformation function $\tilde{f}_{j}=\exp(-j\omega)f_{j}$ , $j=$ 0, 1, 2, $\ldots$ $M-1$ where $\omega>0$ .

(3)

Apply FFT to the set $\tilde{f}_{0}$ , $\tilde{f}_{1}$ , $\tilde{f}_{2}$ , $\ldots$ , $\tilde{f}_{M-1}$ to obtain $\tilde{\chi}_{0}$ , $\tilde{\chi}_{1}$ , $\tilde{\chi}_{2}$ , $\ldots$ , $\tilde{\chi}_{M-1}.$

(4)

Calculate $\tilde{\xi}_{m}=\gamma(\tilde{\chi}_{m})$ , $m=$ 0, 1, 2, $\ldots$ , $M-1$ .

(5)

Apply the inverse FFT to the set $\tilde{\xi}_{0}$ , $\tilde{\xi}_{1}$ , $\tilde{\xi}_{2}$ , $\ldots$ , $\tilde{\xi}_{M-1}$ , to obtain $\tilde{h}_{0}$ , $\tilde{h}_{1}$ , $\tilde{h}_{2}$ , $\ldots$ , $\tilde{h}_{M-1}$ ,

(6)

Untilt to calculate the final compound distribution as $h_{j}=\tilde{h}_{j}\exp(\omega j)$ , $j=$ 0, 1, 2, $\ldots$ , $M-1$ .

The tilting procedure is found to be very effective in reducing the aliasing error. The value of $\omega$ should be large but care should also be taken that it does not lead to underflow or overflow. Embrechts and Frei (2009) stated that the choice $M\omega\approx 20$ works well for standard double precision (eight bytes) calculation. Furthermore, the evaluation of the p.g.f. $\gamma\left(.\right)$ of the frequency distribution may lead to the problem of underflow in case of large frequencies which can be resolved using the methods described in Shevchenko, (2010).

10. FFT for calculating the first moment of the time to ruin

The defective renewal equation satisfied by the time to ruin is given by

$\displaystyle\psi_{1}\left(u\right)=\frac{1}{1+\theta_{1}}\int_{0}^{u}\psi_{1}% \left({u-x}\right)dF_{X,e}\left(x\right)+\frac{1}{c}\int_{u}^{\infty}\psi\left% (x\right)dx$ (26)

If we take Fourier transform on both sides of the equation, it leads to the following methodology adopted to compute the first moment of the time to ruin. The array $g=\left\{{g_{0},g_{1},g_{2},\ldots g_{M-1}}\right\}$ that was used to compute $\psi\left(u\right)$ will be used to compute $\psi_{1}\left(u\right)$ as described in the following algorithm.

Step A: Compute $w_{j}=\sum_{i=1}^{j}g_{i}$ . Step B: Compute $V_{j}=\frac{1}{C}\sum_{i=1}^{j}w_{j}$ .

Since, $\psi_{1}=\frac{1}{1-\varphi}\left({V\ast w}\right)$ , Compute

Step C: $d_{1}=\textit{fft}\left(V\right),d_{2}=\textit{fft}\left(w\right),d_{3}=\frac{% 1}{1-\phi}\left({d_{1}\ast d_{2}}\right)$ . Step D: Invert $d_{3}$ to get $\psi_{1}^{\textit{Computed}}=\frac{1}{1-\phi}\left({V\ast w}\right)$ .

An approximation to $\psi_{1}\left({jh}\right)$ is given by $\psi_{1,j}=\frac{1}{h}\psi_{1}^{\textit{Computed}}\left(j\right)$ .

However, it is observed that the values of $\psi_{1,j}$ thus obtained apparently lack accuracy which can be ascertained from the fact they are significantly different from the values of $E\left(T\right)$ as found in Das and Nath (2019) and Das (2017). Hence, we have devised an alternative method to compute $\psi_{1}\left(u\right)$ using a combination of numerical integration and FFT and the method is described as below.

11. Computing the first moment of the time to ruin using numerical integration along with probability of ultimate ruin computed through FFT

Lin and Willmot (2000) showed that the $k^{th}$ moment of the distribution of the time to ruin is given by

$\displaystyle E\left({T^{k}}\right)=\frac{\psi_{k}\left(u\right)}{\psi\left(u% \right)},k=1,2,3\ldots$ (27)

where

$\displaystyle\psi_{k}\left(u\right)=\frac{k}{\lambda p_{1}\theta_{1}}\left[{% \int_{0}^{u}\psi\left({u-x}\right)\psi_{k-1}\left(x\right)dx+\int_{u}^{\infty}% \psi_{k-1}\left(x\right)dx-\psi\left(u\right)\int_{0}^{\infty}\psi_{k-1}\left(% x\right)dx}\right]$

Let $L :$ the maximum of the aggregate loss process so that $\psi(u)=pr(L>u)$ (see, Bowers et al. (1998), formula (13.6.2)).

In Gerber (1979), it has been shown that

$\displaystyle E\left(L\right)=\int_{0}^{\infty}\psi\left(x\right)dx=\frac{p_{2% }}{2\theta_{1}p_{1}}$ (28)

Formula (6.2.1) of Lin and Willmot (2000) has been simplified in Dickson and Waters (2002) as

$\displaystyle\psi_{1}\left(u\right)=\frac{1}{\lambda p_{1}\theta_{1}}\left\{E% \left(L\right)\delta\left(u\right)-\int_{0}^{u}\psi\left(x\right)\delta\left({% u-x}\right)dx\right\}$ (29)

where $\psi\left(u\right)$ can be evaluated using numerical integration in case, no closed form expression is available for $\psi\left(u\right)$ . The quantity $E\left(L\right)$ is obtained from Eq. (28) and $\delta\left(x\right)=1-\psi\left(x\right)$ .

We have used Eq. (29) to obtain $\psi_{1}\left(u\right)$ . The quantities $\psi\left(x\right)$ and $\delta\left(x\right)$ appearing in the integrand are obtained through FFT as discussed in Section 2. Then Simpson’s 1/3 rd rule for numerical integration is used to solve the integral $\int_{0}^{u}\psi\left(x\right)\delta\left({u-x}\right)dx$ to obtain $\psi_{1}\left(u\right)$ .

12. Results and discussions

Two distributions have been used for modelling the claim severity, namely the Weibull distribution which is a heavy tailed distribution and the Mixture of 3 Exponentials, which is a light tailed distribution. The parameters involved in these distributions have been estimated using a dataset which have been discussed in Das (2017) and Nath and Das (2017).

In computing the FFT, we have used the function fft of the R software. The computational procedure discussed in Section 2 and some parts in Section 4 can be implemented easily via the fft command of R software. In computing the probability of ultimate ruin via FFT, the main difficulty lies in obtaining the equilibrium distribution of the corresponding claim severity distribution and then to discretize it with the objective of inserting it as an input to FFT. In our case, as it is evident, a closed form expression can be found for the equilibrium distribution of Mixture of 3 Exponentials but no closed form expression can be found for the equilibrium distribution of Weibull and it has to be computed numerically.

For discretizing the equilibrium distribution, we have used the central difference discretization. We have taken an illustrative value of the security loading factor as $\theta_{1}=$ 0.3. As used in Pitts, (2006), we have truncated the discretization at $M=$ 65536 and the interval of discretisation is taken as $\delta=$ 0.001 for both Weibull and Mixture of 3 Exponentials distributions. The reference Pitts (2006) has used the interval of discretization as $\delta=$ 4/6 for all the claim severity distributions. As can be realized, the lower the value of $\delta$ , better is the estimate of the probability of ultimate ruin obtained through FFT.

It may be recalled from the earlier discussion, that if we need to compute the probability of ultimate ruin for a value of the initial surplus say, $u=u^{\ast}$ , then the array $\left({g_{0},g_{1},g_{2},\ldots.g_{M-1}}\right)$ computed through FFT is summed upto $j,i.e.,\sum_{i-0}^{j}g_{i}=k\left(\textit{say}\right)$ is obtained where $j=\left({\frac{u^{\ast}}{\delta}}\right)-0.5$ and the probability of ruin for this value of the initial surplus is given by $\psi\left({u^{\ast}}\right)=1-k$ .

Table 1
Ultimate ruin probabilities for the Weibull distribution obtained through FFT and product integration

Value of the initial surplus $u\left(\textit{in Rs}\right)$	$\psi\left(u\right)$ (Through FFT)	$\psi\left(u\right)$ (Through Product Integration)
10	0.7691	0.7688
20	0.7690	0.7684
30	0.7689	0.7679
40	0.7688	0.7675
50	0.7687	0.7671
60	0.7686	0.7667
70	0.7685	0.7662
80	0.7684	0.7658
90	0.7683	0.7654
100	0.7680	0.7650
200	0.7672	0.7607
500	0.7642	0.7480
1000	0.7575	0.7274

Table 2

Ultimate ruin probabilities for the Mixture of 3 Exponentials distribution obtained through FFT and the method of product integration

$u$ (in Rs)	$\psi\left(u\right)$ (Through FFT)	$\psi\left(u\right)$ (Through product integration)
10	0.7691	0.7691
20	0.7690	0.7690
30	0.7690	0.7690
40	0.7689	0.7689
50	0.7688	0.7688
60	0.7687	0.7687
70	0.7686	0.7686
80	0.7685	0.7685
90	0.7684	0.7684
100	0.7683	0.7683
200	0.7674	0.7675
500	0.7647	0.7648
1000	0.7592	0.7605

Table 1 gives the probability of ultimate ruin for the Weibull distribution with parameters $\hat{\theta}=$ 18058.838 and $\hat{\beta}=$ 1.0197 through FFT. As mentioned earlier, these parameters are estimated using the dataset given in Das (2017) and Nath and Das (2017). Probability of ultimate ruin as expected, is found to be decreasing with an increase in initial surplus. This is intuitively logical for with the increase in the initial surplus of the insurance company, the company has more reserves to accommodate claims, thereby reducing its chance of ruin. Similarly, the Table 2, shows the probability of ultimate ruin for Mixture of 3 Exponentials (The estimated parameters of this distribution being $\hat{\lambda}_{1}=$ 1.0667e-05, $\hat{\lambda}_{2}=$ 7.979e-05, $\hat{\lambda}_{3}=$ 1.0056e-05, $\hat{w}_{1}=$ 0.098545, $\hat{w}_{2}=$ 6.7670e-01 and $\hat{w}_{3}=$ 2.245e-01) computed through FFT and even, in this case, the same trend as seen in Weibull is observed in the probability of ultimate ruin, as they decrease with an increase in the initial surplus. These results are consistent with the probability of ultimate ruin for these two distributions obtained through the method of product integration found in Nath and Das (2017) and Das and Nath (2019).

Table 3

FFT calculating the quantile of the compound distribution Weibull $\sim$ Poisson using central difference discretization with the discretization step $\delta=$ 1

$r$	Order of the quantile	Value of the quantile	Order of the quantile	Value of the quantile
	Without tilting		with tilting
14	2.021e-06	16384	2.204e-11	16384
15	5.677e-03	32768	7.264e-10	32768
16	0.455	65536	6.644e-08	65536
17	0.983	131072	1.538e-05	131072
18	0.999	262144	4.803e-03	262144
19	1.000	524288	3.656e-01	524288

Table 4

FFT calculating the quantile of the compound distribution for the Mixture of 3 Exponentials distribution $\sim$ Poisson using central difference discretization with the discretization step $\delta=$ 1

$r$	Order of the quantile	Value of the quantile	Order of the quantile	Value of the quantile
	Without tilting		With tilting
14	0.00004	16384	1.106e-10	16384
15	0.01610	32768	5.360e-09	32768
16	0.17983	65536	6.459e-07	65536
17	0.45389	131072	0.00014	131072
18	0.82290	262144	0.01860	262144
19	0.98819	524288	0.36049	524288

Table 3 shows the quantiles of the compound distribution, Weibull compounded with Poisson (32.427), obtained through FFT whereas the Table 4 shows the quantiles of the compound distribution, Mixture of 3 Exponentials compounded with Poisson (32.427). One very distinct advantage of FFT over Panjer recursion is that it is comparatively much faster than the Panjer Recursion. For example, it executes (2^19 $=$ 524288) computations in an instance of 5 seconds whereas the Panjer Recursion is much slower as can be ascertained from the fact that it takes more than 10 minutes, even to execute 65536 computations.

In computing the quantiles of the aggregate claim amount distribution, it needs to be mentioned that the equilibrium distributions are discretized at $M$ ,with $M$ ranging among $M=2^{n}$ , $n=$ 14, 15, 16, 17, 18 and 19. However, since, areas up to these points given by $\sum_{i=1}^{M}f_{i}$ for Weibull distribution as well as for Mixture of 3 exponentials is less than 1, there is scope for the occurrence of some amount of truncation error in both the cases. It is evident that the value of M should have been much larger but at the same time it needs to be pointed out that executing 65536 iterations is itself very time consuming and requires much memory space and hence, increasing M beyond it, is extremely inconvenient. The interval of discretization is taken to be $\delta=1$ for both of these aggregate claim amount distributions.

Also, it is evident from Table 3, in case of the compound distribution generated by Weibull with these set of parameters, without tilting, 45.65 percent of the distribution lies below 65536( $2^{16}$ ) and with tilting, only 6.643e-06 percent of the distribution lies below this point. Whereas, for this same distribution, without tilting, 100 percent of the distribution lies below 524288 ( $2^{19}$ ) and with tilting, a mere 36.56 percent of this compound distribution lies below this point. Similarly, from Table 4 for the compound distribution formed from the Mixture of 3 Exponentials, i.e., the aggregate claim distribution obtained in case of Mixture of 3 Exponentials, without tilting, 17.98 percent of the distribution lies below 65536 whereas, 98.81 percent of the distribution lies below 524288. With tilting, for this compound distribution, 6.46e-05 percent of the distribution lies below 65536 and 36.049 percent of the distribution lies below 524288. For the full coverage of the compound distributions, a much larger truncation point is actually required, however as already mentioned in earlier discussion, implementing this is not practicable, considering the constraints of time and memory space. Nevertheless, the quantiles, even, if they are of the low order, are casting some light on the corresponding compound distributions obtained from the Weibull and the Mixture of 3 Exponentials.

The results obtained through FFT are comparable to those obtained through Panjer recursion only when it is subjected to tilting. The results obtained though Panjer Recursion are considered to be the actual values and the FFT is compatible to it only under tilting. This justifies the efficiency of the tilting procedure. As discussed in Shevchenko, (2010 a, b), quantiles computed via FFT are comparable with the results obtained through Panjer recursion scheme, only, when they are subjected to tilting because the tilting has the impact of reducing the decay rate of the tail of the claim severity distribution.(see Table 9.14 and Table 9.10 of Das (2017) for the quantiles obtained through Panjer Recursion for these distributions). FFT, as a technique to compute the quantiles for all of the compound distributions is found to be remarkably fast but as evident, it leads to a significant loss of accuracy unless it is subjected to tilting. The contribution of tilting in improving the efficiency of the naive FFT is found by comparing the results obtained through FFT with those obtained through the Panjer Recursive algorithm.

Table 5

First moment of the time to ruin for the Weibull distribution

Initial surplus (u in Rs)	First moment of the time to ruin $E\left({T{\|}T<\infty}\right)$ , (Mean in years) (By FFT)	First moment of the time to ruin, $E(T\|T<\infty)$ (Mean in years) (By combination of FFT & numerical integration)	First moment of the time to ruin, $E(T\|T<\infty)$ (Mean in years) (As extracted from Das & Nath (2019))
10	30.005	0.10088	0.10088
20	30.009	0.10093	0.10092
30	30.013	0.10097	0.10097
40	30.017	0.10102	0.10101
50	30.021	0.10106	0.10105
60	30.025	0.10111	0.10109
70	30.029	0.10115	0.10114
80	30.033	0.10119	0.10118
90	30.037	0.10124	0.10122
100	30.042	0.10128	0.10127
200	30.082	0.10173	0.10170
500	30.2073	…	0.10299

Table 6

First moments of the time to ruin for the Mixture of 3 Exponentials

Initial surplus (u in Rs)

First moment of the time to ruin

E\left(T|T<\infty\right)

, (Mean in years)

(By FFT)

First moment of the time to ruin,

E(T|T<\infty)

(Mean in years)

(By combination of FFT &

numerical integration)

First moment of the time to ruin,

E(T|T<\infty)

(Mean in years)

(As extracted from

Das (2017))

30.004

0.25684

0.25681

30.008

0.25696

0.25690

30.012

0.25708

0.25698

30.015

0.25720

0.25707

30.019

0.25732

0.25716

30.023

0.25744

0.25724

30.026

0.25756

0.25733

30.030

0.25768

0.25742

30.034

0.25780

0.25751

100

30.037

0.25792

0.25760

200

30.074

0.25911

0.25847

500

30.182

…

0.26107

Table 5 gives the first moment of the time to ruin for the Weibull distribution with parameters $\hat{\theta}=$ 18058.838 and $\hat{\beta}=$ 1.0197 and $\lambda=$ 32.427. The first column gives the first moment of the time to ruin computed by using FFT for solving the renewal Eq. (7) as described in Section 4.1. The trend observed in these computed moments are consistent with what can be expected in practice. The mean of the time to ruin is found to be increasing with an increase in the initial surplus. This is intuitively logical for the induction of larger surpluses should prolong the time to ruin, if ruin at all, occurs. However, it needs to be mentioned that these are just the approximations to the first moment of the time to ruin for Weibull since the renewal Eq. (7) can’t be solved analytically, in case, the claim severity is Weibull. The numerical computation of these moments is done in Das and Nath (2019) and the results are cited in the fourth column of Table 5. Since, the values computed through FFT and the values obtained earlier in literature seem to be significantly different, we sort ways to increase the efficiency of FFT in computing first moment of the time to ruin by combining it with a numerical integration procedure as discussed in Section 4.2. These values are cited in the third column of Table 5 and they are remarkably similar to the values computed earlier in literature and this leads us to conclude that FFT as a technique, to compute the first moment of the time to ruin works more efficiently when combined with a prescribed numerical integration procedure. Similar results are shown for the Mixture of 3 Exponentials in Table 6. It may be noted that for the Mixture of 3 Exponentials, exact values for the first moment of the time to ruin are available in Das (2017) and they are tabulated in fourth column of the Table 6 for comparison.

An interpretation of a typical value in Table 5 is that for the claim severity distributed as the Weibull, starting with an initial surplus of Rs 100 and when computed through the combination of FFT and numerical integration, the first moment suggests that it would on the average, take 0.10128 years for the surplus process to be less than or equal to zero for the first time, thereby leading to ruin in the sense of its definition. However, with the use of FFT alone, this value comes out to be almost 30 years.

The occurrence of numerical error in this procedure of finding the first moment through FFT is evident since, there is scope for evolution of error from multiple sources. For example, in using just FFT for computing the first moment, there is error in computing the $g_{i}$ s which get transferred to the series on $w_{i}$ and $v_{i}$ . In addition to this, the inherent error in computing the FFT of $w_{i}$ and $v_{i}$ is inevitable and it does come into existence. Again, inverting the FFT of the product of the FFT of $w_{i}$ and the FFT of $v_{i}$ are done to get the final output. This entire process raises the possibility of the accumulation of error from various sources.

Likewise, when FFT is combined with numerical integration to get the first moment of the time to ruin, we may observe that there are two distinguishable sources of error. Observing the numerical integral in Eq. (29), we see that apart from inheriting the error of integration in using the Simpson’s 1/3 rd rule, the other source for error is the fact that the integrand itself is a function of $\psi\left(x\right)$ and $\delta\left(x\right)$ , which are themselves computed numerically and then, used as inputs for the subsequent numerical integration. Since, two numerical procedures are merged, one into the another, this computation also takes significantly large time as the integrand has to be computed at a large no. of points and then combined to get the final output. It is also evident that as $u$ increases, the computational time increases exponentially thereby rendering us almost impossible to compute the first moment of the time to ruin beyond $u$ equal to 200 in a 64 bit personal computer. All the computations are done numerically using the R software (R core team, 2013).

13. Concluding remarks

In this work, we have tried to illustrate the applications of FFT in obtaining some of the important quantities in the actuarial domain under the classical risk model. The probability of ultimate ruin is an important quantity for an insurance company in predicting its prospect for insolvency and there are many methods cited in literature to compute it. Basically, the fact that an integral equation involving convolution can be solved using Fourier transform has been the base, on which the discrete counterpart of Fourier transform i.e., FFT has been used to compute the probability of ultimate ruin which is obtained as a solution of an integro-differential equation. Similarly, an insight into the aggregate claim amount distribution is essential for an insurance company to articulate proper decisions on the amount of reserves, it needs to hold to pay the claims made by its customers. Here, the aggregate claim amount distribution is basically a compound distribution with the claim severity distribution and the counting distribution as its ingredients. The exact form of the compound distribution is not known analytically for most of the claim severity distributions but its Fourier transform or its characteristics function is either known explicitly or can be found numerically. Hence, the compound distribution can be identified by inverting its Fourier transform and this form the base on which FFT is exploited to compute the quantiles of the aggregate claim distribution. Similarly, the first moment of the time to ruin gives some insight into the probable time, the ruin (insolvency) might occur, which in most cases, is infinite for a stable company. In this case also, it may be noted that the first moment of the time to ruin can be obtained as a solution of an integro-differential equation and hence, FFT can be utilized to solve it numerically to get an approximate value for this moment.

References

Asmussen

(2000). Ruin Probabilities, World Scientific, Singapore.

Bowers

N.L.

Gerber

H.U.

Hickman

J.C.

Jones

D.A.

, & Nesbitt

C.J.

(1998). Actuarial Mathematics, Society of Actuaries, Itasca, Illinois.

Chernobai

A.S.

Rachev

S.T.

Fabozzi

F.J.

(2007). Operational Risk: A Guide to Basel II Capital Requirements, Models and Analysis, John Wiley and Sons, New Jersey.

Das

(2017). Actuarial Risk Process and the Probability of Ruin for Non Life Insurance [Unpublished PHD Thesis], Gauhati University, available online at https://shodhganga.inflibnet.ac.in.

Das

, & Nath

D.C.

(2019). Weibull Distribution as an Actuarial Risk Model: Computation of its Probability of ultimate Ruin and the moments of the time to ruin, deficit at ruin and Surplus prior to Ruin. Journal Of Data Science, 17(1), 161-194.

Dickson

D.C.M.

, & Waters

H.R.

(2002). The distribution of the time to ruin in the classical risk model. ASTIN Bulletin, 32, 299-313.

Embrechts

Grubel

, & Pitts

S.M.

(1993). Some applications of the Fast Fourier Transform algorithm in Insurance Mathematics. Statistica Neerlandica, 47, 59-75.

Embrechts

, & Frei

(2009). Panjer recursion Versus FFT for Compound distributions. Mathematical Methods of Operations Research, 69(3), 497-508.

Gerber

H.U.

(1979). An Introduction to Mathematical Risk Theory, Homewood, Richard, D. Irwin Inc.

10.

Grandell

(1991). Aspects of Risk Theory, Springer, New York.

11.

Grubel

, & Hermesmeir

(1999). Computation of Compound distributions 1: aliasing error and Exponential tilting. ASTIN Bulletin, 29(2), 197-214.

12.

Grubel

, & Hermesmeir

(2000). Computation of compound distributions II: discretization errors and Richardson extrapolation. ASTIN Bulletin, 30, 309-331.

13.

Klugman

S.A.

Panjer

H.H.

, & Willmot

G.E.

(1998). Loss Models: From Data to Decisions, Wiley, New York.

14.

Lima

F.D.P.

Garcia

J.M.A.

, & Dos Reis

A.D.E.

(2002). Fourier/Laplace Transforms and Ruin Probabilities. ASTIN Bulletin, 32, 91-105.

15.

Lin

X.S.

, & Willmot

G.E.

(1999). Analysis of a Defective Renewal Equation Arising in Ruin Theory. Insurance Mathematics and Economics, 25, 63-84.

16.

Lin

X.S.

, & 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.

17.

Nath

D.C.

, & Das

(2017). Modelling of Claim Severity through the mixture of Exponential Distribution and Computation of its Probability of Ultimate Ruin. Thailand Statistician, 15(2), 128-148.

18.

Panjer

H.H.

, & Willmot

G.E.

(1986). Computational Aspects Of Recursive Evaluation of Compound Distributions. Insurance: Mathematics and Economics, 5, 113-116.

19.

Panjer

H.H.

, & Willmot

G.E.

(1992). Insurance Risk Models, Society of Actuaries, Schaumburg.

20.

Panjer

H.H.

, & Wang

(1993). On the Stability of Recursive formulas. ASTIN Bulletin, 23(2), 227-258.

21.

Pitts

S.M.

(2006). The Fast Fourier Transform Algorithm in Ruin theory for the Classical Risk Model. HERMIS: An International Journal of Computer Mathematics and its Applications, 7, 80-94.

22.

R Core Team (2013). R: A language and environment for Statistical computing, R Foundation for Statistical computing, Vienna, Austria, ISBN-3-900051-07-0, URL http://www.R-project.org/.

23.

Robertson

(1992). The Computation of aggregate Loss distributions, Proceedings of the Casualty Actuarial Society, pp. 57-133.

24.

Shevchenko

P.V.

(2008). Estimation of Operational risk capital charge under parameter uncertainty. The Journal of Operational Risk, 3(1)), 51-53.

25.

Shevchenko

P.V.

(2010a). Implementing Loss Distribution Approach for Operational Risk. Applied Stochastic Models in Business and Industry, 26(3), 277-307.

26.

Shevchenko

P.V.

(2010b). Calculation Of Aggregate Loss Distributions. The Journal of Operational Risk, 5(2), 3-40.

Computation of the probability of ultimate ruin and some other actuarial quantities under the classical risk model via Fast Fourier Transform

Abstract

Keywords

1. Introduction

2. The classical risk model

Table 1 Ultimate ruin probabilities for the Weibull distribution obtained through FFT and product integration

References

Table 1
Ultimate ruin probabilities for the Weibull distribution obtained through FFT and product integration