Phase-type Approximation of the Gerber-Shiu Function

TL;DR

This paper introduces a phase-type approximation method for the Gerber-Shiu function, enabling more tractable risk measurement in insurance models by fitting surplus processes with phase-type Levy processes.

Contribution

It presents a novel closed-form approximation of the Gerber-Shiu function using phase-type Levy processes, advancing computational methods in risk theory.

Findings

Provides a closed-form approximation for the Gerber-Shiu function.

Demonstrates the effectiveness of the approximation through numerical results.

Enhances computational efficiency in risk assessment models.

Abstract

The Gerber-Shiu function provides a way of measuring the risk of an insurance company. It is given by the expected value of a function that depends on the ruin time, the deficit at ruin, and the surplus prior to ruin. Its computation requires the evaluation of the overshoot/undershoot distributions of the surplus process at ruin. In this paper, we use the recent developments of the fluctuation theory and approximate it in a closed form by fitting the underlying process by phase-type Levy processes. A sequence of numerical results are given.

Tables

Table 1. Table 1: Parameters of the hyperexponential distribution fitted to Weibull( 0.6 0.6 0.6 , 0.665 0.665 0.665 ) and Pareto(1.2,5) (taken from Tables 3 of [ 14 ] ).

$i$ ${\alpha }_i$ ${\eta }_i$ 1 0.029931 676.178 2 0.093283 38.7090 3 0.332195 4.27400 4 0.476233 0.76100 5 0.068340 0.24800 6 0.000018 0.09700 – – –	$i$ ${\alpha }_i$ ${\eta }_i$ $i$ ${\alpha }_i$ ${\eta }_i$ 1 8.37E-11 8.3E-09 8 0.000147 0.0020 2 7.18E-10 6.8E-08 9 0.001122 0.0100 3 5.56E-09 3.9E-07 10 0.008462 0.0570 4 4.27E-08 2.2E-06 11 0.059768 0.3060 5 3.27E-07 1.2E-05 12 0.307218 1.5460 6 2.50E-06 6.5E-05 13 0.533823 6.5160 7 1.92E-05 3.5E-04 14 0.089437 23.304
(i) Weibull$(0.6,0.665)$	(ii) Pareto$(1.2,5)$

Equations

$$\psi (s):=\log \mathbb{E}\left[e^{s{X}_1}\right]=cs+\frac{1}{2}{\sigma }^2{s}^2+{\int }_{(-\infty ,0)}(e^{sz}-1-sz{1}_{\{z>-1\}})\Pi (\mathrm{d}z),s\ge 0,$$

$$\sigma =0\text{and}{\int }_{(-\infty ,0)}(1\wedge |z|)\Pi (\mathrm{d}z)<\infty ;$$

$$\psi (s)=\mu s+{\int }_{(-\infty ,0)}(e^{sz}-1)\Pi (\mathrm{d}z),$$

$$\mu :=c-{\int }_{(-1,0)}z\Pi (\mathrm{d}z).$$

$${\tau }_0^{-}:=\inf \left\{t\ge 0:X_t<0\right\}.$$

$$\displaystyle GS_{f}(x,q):=\mathbb{E}^{x}\left[e^{-q\tau_{0}^{-}}f\big{(}X_{\tau_{0}^{-}},X_{\tau_{0}^{-}-}\big{)};\tau_{0}^{-}<\infty\right].$$

$$K^{(q)}(x,\mathrm{d}y,\mathrm{d}z):={\mathbb{E}}^x\left[e^{-q{\tau }_0^{-}};X_{{\tau }_0^{-}}\in \mathrm{d}y,X_{{\tau }_0^{-}-}\in \mathrm{d}z,{\tau }_0^{-}<\infty \right],x,z>0,y<0,$$

$$G{S}_f(x,q)={\int }_{(0,\infty )}{\int }_{(-\infty ,0)}f(y,z)K^{(q)}(x,\mathrm{d}y,\mathrm{d}z);$$

$$R^{(q)}(x,\mathrm{d}z):={\int }_0^{\infty }{e}^{-qt}{\mathbb{P}}^x\{X_t\in \mathrm{d}z,{\tau }_0^{-}>t\}\mathrm{d}t,x,z>0,$$

$$R^{(q)}(x,\mathrm{d}z)=r^{(q)}(x,z)\mathrm{d}z,x,z>0;$$

$$K^{(q)}(x,\mathrm{d}y,\mathrm{d}z)=\Pi (\mathrm{d}y-z)r^{(q)}(x,z)\mathrm{d}z.$$

$${\int }_0^{\infty }{e}^{-sx}{W}^{(q)}(x)\mathrm{d}x=\frac{1}{\psi (s)-q},s>\Phi (q)$$

$$\Phi (q):=\sup \{s\ge 0:\psi (s)=q\},q\ge 0.$$

$$\displaystyle\begin{split}W^{(q)}(0)&=\left\{\begin{array}[]{ll}0,&\textrm{if $X$ is of unbounded variation}\\ \frac{1}{\mu},&\textrm{if $X$ is of bounded variation}\end{array}\right\},\\ W^{(q)^{\prime}}(0+)&=\left\{\begin{array}[]{ll}\frac{2}{\sigma^{2}},&\textrm{if }\sigma>0\\ \infty,&\textrm{if }\sigma=0\;\textrm{and}\;\Pi(-\infty,0)=\infty\\ \frac{q+\Pi(-\infty,0)}{\mu^{2}},&\textrm{if $X$ is compound Poisson}\end{array}\right\}.\end{split}$$

$${\tau }_b^{-}:=\inf \left\{t>0:X_t<b\right\}\text{and}{\tau }_b^{+}:=\inf \left\{t>0:X_t>b\right\},b\in \mathbb{R}.$$

$$\begin{array}{rl}{\mathbb{E}}^x\left[e^{-q{\tau }_b^{+}}{1}_{\left\{{\tau }_b^{+}<{\tau }_0^{-}\right\}}\right]& =\frac{W^{(q)}(x)}{W^{(q)}(b)},\\ {\mathbb{E}}^x\left[e^{-q{\tau }_0^{-}}{1}_{\left\{{\tau }_b^{+}>{\tau }_0^{-}\right\}}\right]& =Z^{(q)}(x)-Z^{(q)}(b)\frac{W^{(q)}(x)}{W^{(q)}(b)},\\ {\mathbb{E}}^x\left[e^{-q{\tau }_0^{-}}\right]& =Z^{(q)}(x)-\frac{q}{\Phi (q)}{W}^{(q)}(x),\end{array}$$

$$Z^{(q)}(x)$$

$$r^{(q)}(x,z)=e^{-\Phi (q)z}{W}^{(q)}(z)-W^{(q)}(x-z),x,z>0.$$

$$Y=\{Y_t;t\ge 0\}$$

$$\mathbf{Q}:=\left[\begin{array}{cc}\mathbf{T}& \mathbf{t}\\ 0& 0\end{array}\right].$$

$$F(z;\mathbf{\alpha },\mathbf{T})=1-\mathbf{\alpha }{e}^{\mathbf{T}z}1\text{and}f(z;\mathbf{\alpha },\mathbf{T})=\mathbf{\alpha }{e}^{\mathbf{T}z}\mathbf{t},z>0.$$

$$X_t-X_0=\mu t+\sigma B_t-\sum _{n=1}^{N_t}{Z}_n,0\le t<\infty ,$$

$$\psi (s)=\mu s+\frac{1}{2}{\sigma }^2{s}^2+\lambda \left(\mathbf{\alpha }(s\mathbf{I}-\mathbf{T}{)}^{-1}\mathbf{t}-1\right),$$

$${\mathcal{I}}_q$$

$$W^{(q)}(x)=\frac{e^{\Phi (q)x}}{{\psi }^{\prime }(\Phi (q))}-\sum _{i=1}^n\sum _{k=1}^{m_i}{C}_{i,q}^{(k)}\frac{x^{k-1}}{(k-1)!}{e}^{-{\xi }_{i,q}x},x\ge 0,$$

$$C_{i,q}^{(k)}$$

$$W^{(q)}(x)=\frac{e^{\Phi (q)x}}{{\psi }^{\prime }(\Phi (q))}-\sum _{i=1}^n{C}_{i,q}{e}^{-{\xi }_{i,q}x},x\ge 0,$$

$$C_{i,q}$$

$$f(z)=\sum _{j=1}^m{\alpha }_j{\eta }_j{e}^{-{\eta }_{j}z},z>0,$$

$$\psi (s)=\mu s+\frac{1}{2}{\sigma }^2{s}^2-\lambda \sum _{j=1}^m{\alpha }_j\frac{s}{{\eta }_j+s}.$$

Full text

Keywords: Risk management, applied probability

1 Introduction

The fundamental objective of the actuarial ruin theory is to measure the vulnerability to insolvency. Typically, the surplus of an insurance company is modeled by a stochastic process, and ruin occurs at the first time it goes below a certain threshold. A most classical and important quantity of interest is the ruin probability, and the Gerber-Shiu function is its generalization; it is given as an expected discounted value of the cost function that is dependent on the ruin time, the deficit at ruin, and the surplus prior to ruin. The evaluation of the Gerber-Shiu function involves that of the overshoot and undershoot distributions at the first down-crossing time that do not admit explicit expressions. Hence its computation is in general a challenging task.

In ruin theory, the surplus process is commonly modeled by a stochastic process with downward jumps. Due to the premiums received from the insured persons, the surplus tends to increase constantly. On the other hand, it experiences sudden downward jumps due to the insurance payments. The classical Cramér-Lundberg model uses a compound Poisson process with downward jumps. Its generalization called the Sparre-Andersen model modifies it by allowing the arrival of the claims to follow a general renewal process.

In the last decade, significant progress has been made in insurance mathematics and ruin theory, thanks to the development of the theory of Lévy processes [9, 19]. In particular, many results on the Cramér-Lundberg model have been generalized for a general spectrally negative Lévy process, or the Lévy process with only downward jumps; see, e.g., [4, 6, 7, 24]. This generalization enables one to construct more realistic models; one can, for example, introduce noise by including Brownian motion and/or infinitesimal jumps of infinite activity/variation. Using the so-called scale function, one can express concisely many quantities of interest for a general spectrally negative Lévy process.

The objective of this paper is to give an approximation to the Gerber-Shiu function using the theory of scale functions. By the compensation formula of the Lévy process, the Gerber-Shiu function admits an expression as a (double) integral with respect to the resolvent measure and the Lévy measure. Because the resolvent can be written using the scale function, at least in theory the computation of the Gerber-Shiu function boils down to that of the scale function.

However, a major obstacle still remains in putting these in practice because scale functions are in general known only up to their Laplace transforms, and only a few cases admit explicit expressions. The most straightforward approach of computing the scale function is to apply numerical Laplace inversion as in [18, 33]. However, this approach is not suitable for the computation of the Gerber-Shiu function because the numerically approximated scale functions need to be further integrated with respect to the Lévy measure. In particular, the undershoot density that is essential in the computation of the Gerber-Shiu function tends to have a very peculiar form with a possible spike. The approximation hence requires a high precision in computing the scale function and minimal discretization errors in numerical integration.

In this paper, we adopt a phase-type fitting of [13] by using the scale function for the class of spectrally negative phase-type Lévy processes, or Lévy processes with negative phase-type-distributed jumps. Consider a continuous-time Markov chain with some initial distribution and state space consisting of a single absorbing state and a finite number of transient states. A phase-type distribution is that of the first entry time to the absorbing state. As has been discussed in [13, 18], the scale function of this process becomes the sum of (possibly complex) exponentials; it can be integrated with respect to the Lévy measure analytically to obtain an explicit form of the Gerber-Shiu function. More importantly, the class of phase-type distributions is dense in the class of all positive-valued distributions. Consequently, the Gerber-Shiu function of any given spectrally negative Lévy process can be approximated in a closed form by that of an approximating spectrally negative phase-type Lévy process. Our aim is to evaluate numerically the practicability of this approach.

In our numerical results, we focus on the case where the Lévy measure is finite and has a completely monotone density. In this case, the jump size distribution can be approximated by a special class of phase-type distributions, called the hyperexponential distributions. While fitting a phase-type distribution for a general distribution is often difficult, fitting a hyperexponential distribution for the ones with completely monotone densities can be efficiently done, for example by the algorithm by [14] that is guaranteed to converge to the desired distribution. The class of Lévy processes with completely monotone Lévy densities includes, for example, a subset of compound Poisson processes, variance gamma [26, 27], CGMY [10], generalized hyperbolic [12] and normal inverse Gaussian [5] processes.

In order to evaluate our approach, we obtain, for the hyperexponential case, the closed-form expressions of the (discounted) overshoot/undershoot distributions at the first down-crossing time, and use it to approximate those for the processes with Weibull/Pareto-type jumps. The obtained results are then compared with those obtained by Monte Carlo simulation.

To our best knowledge, this is the first paper on the numerical evaluation of the Gerber-Shiu function via phase-type fitting. As the Gerber-Shiu measure is sensitive to approximation errors, it is important to evaluate its numerical performance. Recently, the resolvents of related extensions of the Lévy process have been developed, and it is reasonable to conjecture that the Gerber-Shiu function of these can be approximated precisely in the same way.

The rest of the paper is organized as follows. Section 2 reviews the spectrally negative Lévy process, the Gerber-Shiu function, and the scale function. Section 3 gives a summary of [13] on the spectrally negative phase-type Lévy process and its scale function. Section 4 computes the Gerber-Shiu measure of the hyperexponential Lévy process as an approximation for the case the jump size distribution admits a completely monotone density. We evaluate the performance using numerical results in Section 5. Section 6 concludes with remarks on the cases of other variants of spectrally negative Lévy processes.

2 Gerber-Shiu Functions for Spectrally Negative Lévy Processes

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space hosting a spectrally negative Lévy process $X=\left\{X_{t};t\geq 0\right\}$ that models the surplus of a company. Let $\mathbb{P}^{x}$ be the conditional probability under which $X_{0}=x$ (and also $\mathbb{P}\equiv\mathbb{P}^{0}$), and $\mathbb{F}:=\left\{\mathcal{F}_{t}:t\geq 0\right\}$ the filtration generated by $X$. The process $X$ is uniquely characterized by its Laplace exponent

[TABLE]

where $\sigma\geq 0$ is the diffusion (Brownian motion) coefficient and $\Pi$ is a Lévy measure with the support $(-\infty,0)$ that satisfies the integrability condition $\int_{(-\infty,0)}(1\wedge|z|^{2})\Pi({\rm d}z)<\infty$. It has paths of bounded variation if and only if

[TABLE]

see, for example, Lemma 2.12 of [19]. In this case, we can rewrite the Laplace exponent (2.1) by

[TABLE]

with

[TABLE]

We disregard the case when $X$ is the negative of a subordinator (or decreasing a.s.).

2.1 Gerber-Shiu functions

Define the ruin time as the first time the surplus goes below zero:

[TABLE]

Here and throughout the paper, we use the convention that $\inf\emptyset=\infty$. On the event $\{\tau_{0}^{-}<\infty\}$, the random variables $X_{\tau_{0}^{-}}$ and $X_{\tau_{0}^{-}-}$ model, respectively, the deficit at ruin and surplus immediately before ruin.

Fix $f:(-\infty,0]\times[0,\infty)\rightarrow[0,\infty)$ bounded and measurable. We define the Gerber-Shiu function

[TABLE]

With the Gerber-Shiu measure

[TABLE]

we can write

[TABLE]

see pages 4 and 5 of [20].

Using the compensation formula (see Theorem 4.4 of [19]), this can be written in terms of the $q$-resolvent measure of $X$ killed on exiting $[0,\infty)$:

[TABLE]

which is known to admit a density $r^{(q)}$ for the case of spectrally negative Lévy processes such that

[TABLE]

see (2.10) below for the form of $r^{(q)}$. As in Section 1.3 of [20], we can write

[TABLE]

Hence the computation of the Gerber-Shiu measure boils down to that of the resolvent measure.

2.2 Scale functions

In order to compute the resolvent measure (2.2), we shall introduce the scale function.

Fix $q\geq 0$. The scale function $W^{(q)}:\mathbb{R}\rightarrow[0,\infty)$ of $X$ is a function whose Laplace transform is given by

[TABLE]

where

[TABLE]

On the negative half line, it is assumed that $W^{(q)}(x)=0$.

Regarding the smoothness of the scale function, if the Lévy measure has no atoms or $X$ is of unbounded variation, then $W^{(q)}\in C^{1}(0,\infty)$; if it has a Gaussian component ($\sigma>0$), then $W^{(q)}\in C^{2}(0,\infty)$. See [11] for other known results on the smoothness.

The behavior in the neighborhood of zero is given as follows. As in Lemmas 4.3 and 4.4 of [23], for every $q\geq 0$,

[TABLE]

We plot in Figure 1 sample plots of the scale function and its derivative for the bounded variation case with $\sigma=0$ and the unbounded variation case with $\sigma>0$.

The most well-known application of the scale function can be found in the two-sided exit identities. Let us define the first down- and up-crossing times, respectively, of $X$ by

[TABLE]

Then, for any $b>0$ and $x\leq b$,

[TABLE]

where

[TABLE]

For a comprehensive account on the scale function, see among others [18, 19].

2.3 Resolvent via the scale function

Recall, as in our discussion in Section 2.1, that the Gerber-Shiu function can be written in terms of the resolvent measure. This can be concisely written in terms of the scale function as follows: by Corollary 8.8 of [19], the resolvent density (2.3) can be written

[TABLE]

Now in view of the identity (2.4), the computation of the Gerber-Shiu measure boils down to that of the scale function.

3 Scale Functions for Spectrally Negative Phase-type Lévy Processes

As discussed in the previous section, the scale function is defined by its Laplace transform and for its computation the Laplace transform (2.5) must be inverted. Here, we review the results on [13] for a special class of Lévy processes where it can be inverted analytically.

3.1 Phase-type distribution

Consider a continuous-time Markov chain

[TABLE]

with a finite state space $\{1,\ldots,m\}\cup\{\Delta\}$ where $1,\ldots,m$ are transient and $\Delta$ is absorbing. Its initial distribution is given by a simplex ${\bm{\alpha}}=[\alpha_{1},\ldots,\alpha_{m}]$ such that $\alpha_{i}=\mathbb{P}\left\{Y_{0}=i\right\}$ for every $i=1,\ldots,m$. The intensity matrix ${\bm{Q}}$ is partitioned into the $m$ transient states and the absorbing state $\Delta$, and is given by

[TABLE]

Here ${\bm{T}}$ is an $m\times m$-matrix called the phase-type-generator, and ${\bm{t}}=-{\bm{T}}{\bm{1}}$ where ${\bm{1}}=[1,\ldots,1]^{\prime}$. A distribution is called phase-type with representation $(m,{\bm{\alpha}},{\bm{T}})$ if it is the distribution of the absorption time to $\Delta$ in the Markov chain described above. It is known that ${\bm{T}}$ is non-singular and thus invertible; see [2]. Its distribution and density functions are given, respectively, by

[TABLE]

3.2 Phase-type Lévy processes

Let $X=\left\{X_{t};t\geq 0\right\}$ be a spectrally negative Lévy process of the form

[TABLE]

for some $\mu\in\mathbb{R}$ and $\sigma\geq 0$ (with $\mu>0$ when $\sigma=0$ so that it is not the negative of a subordinator). Here $B=\{B_{t};t\geq 0\}$ is a standard Brownian motion, $N=\{N_{t};t\geq 0\}$ is a Poisson process with arrival rate $\lambda$, and $Z=\left\{Z_{n};n=1,2,\ldots\right\}$ is an i.i.d. sequence of phase-type-distributed random variables with representation $(m,{\bm{\alpha}},{\bm{T}})$. These processes are assumed mutually independent. Its Laplace exponent is then

[TABLE]

which is analytic for every $s\in\mathbb{C}$ except at the eigenvalues of ${\bm{T}}$.

Let $\mathcal{I}_{q}$ be the set of (the sign-changed) negative roots:

[TABLE]

where $\mathcal{R}(z)$ is the real part of $z\in\mathbb{C}$. Let $n$ denote the number of different roots in $\mathcal{I}_{q}$ and $m_{i}$ the multiplicity of a root $\xi_{i,q}$ for $i=1,\ldots,n$.

Because the Laplace exponent (3.2) has a rational form, it can be inverted analytically by partial fraction decomposition. Hence, in view of (2.5), the scale function can be obtained.

Proposition 3.1** (Section 5.4 of [18] and Proposition 2.1 of [13]).**

Suppose $q\geq 0$ and $\psi^{\prime}(0+)<0$ if $q=0$. Then the scale function is written

[TABLE]

where

[TABLE]

In particular, if all the roots in $\mathcal{I}_{q}$ are distinct, then

[TABLE]

where

[TABLE]

3.3 Approximation results

It is known that the class of phase-type distributions is dense in the class of all positive-valued distributions. Using this, Proposition 1 of [3] shows that, there exists, for any spectrally negative Lévy process $X$, a sequence of spectrally negative phase-type Lévy processes $X^{(n)}$ converging to $X$ in $D[0,\infty)$. In other words, $X_{1}^{(n)}\rightarrow X_{1}$ in distribution by Corollary VII 3.6 of [15]; see also [32].

Using these results, Egami and Yamazaki [13] study the convergence of the corresponding scale function, and show that the approximation is in most cases very accurate. On the other hand, as they also point out, no existing algorithm is guaranteed to construct a converging sequences and fitting phase-type distributions can get difficult; see, e.g., the case of fitting for the uniform distributed jump size in [13].

On the other hand, as we shall discuss next, it is guaranteed to work for the case the jump size admits a completely monotone density.

3.4 Hyperexponential case

As an important example where all the roots in $\mathcal{I}_{q}$ are distinct and real, we consider the case where $Z$ has a hyperexponential distribution with a density function

[TABLE]

for some $0<\eta_{1}<\cdots<\eta_{m}<\infty$ and $\alpha_{j}>0$ for $1\leq j\leq m$ such that $\alpha_{1}+\cdots+\alpha_{m}=1$; this is the phase-type distribution with its Markov chain such that the $m$ transient states are connected only to $\Delta$. Its Laplace exponent (2.1) is then

[TABLE]

Notice in this case that $-\eta_{1}$, …, $-\eta_{m}$ are the poles of the Laplace exponent. Furthermore, all the roots in $\mathcal{I}_{q}$ are distinct and real and satisfy the following interlacing condition for every $q>0$:

for $\sigma>0$, there are $m+1$ roots $-\xi_{1,q},\ldots,-\xi_{m+1,q}$ such that

[TABLE] 2. 2.

for $\sigma=0$, there are $m$ roots $-\xi_{1,q},\ldots,-\xi_{m,q}$ such that

[TABLE]

Because all roots are real and distinct, the scale function can be written as (3.5).

Recall that a density function $f$ of a positive valued random variable is called completely monotone if all the derivatives exist and, for every $n\geq 1$,

[TABLE]

where $f^{(n)}$ denotes the $n^{th}$ derivative of $f$.

Feldmann and Whitt [14] showed that if a density function is completely monotone, then it can be approximated by those of hyperexponential distributions. As shown by [8], every completely monotone density function is a mixture of exponential density functions, and this implies that, for any distribution with a completely monotone density, there exists a sequence of hyperexponential distributions converging to it. The class of distributions with completely monotone densities contains a number of distributions such as the (subset of) Pareto distribution, the Weibull distribution, and the gamma distribution. Feldmann and Whitt [14] took advantage of this fact and proposed a recursive algorithm for fitting hyperexponential distributions to these distributions. We refer the reader to [1, 17] for other approximation methods.

4 Computation of the Gerber-Shiu Function

We shall now consider the approximation of the Gerber-Shiu function using the fitted scale functions. Here, we consider the hyperexponential case as discussed in Section 3.4; in this case, because the scale function and Lévy measure are both written as mixtures of exponential functions, we attain closed-form expressions.

4.1 Overshoot and undershoot distributions

Recall from (2.4) regarding the equivalence of the Gerber-Shiu measure and the product of the Lévy measure and the resolvent measure. If we define for all $A\in\mathcal{B}(-\infty,0)$ and $B\in\mathcal{B}(0,\infty)$,

[TABLE]

Combining (2.4) and (2.10), we can write

[TABLE]

where $\overline{\Pi}$ is the Lévy measure of the dual process $-X$.

When $X$ is a phase-type Lévy process, as we have studied in the previous section $W^{(q)}(x)$ can be written as a sum of (possibly complex) exponentials. In particular, if it is hyperexpontial, we can write $\overline{\Pi}({\rm d}u)=\lambda\sum_{j=1}^{m}\alpha_{j}\eta_{j}e^{-\eta_{j}u}{\rm d}u$ for all $u\in(0,\infty)$, and hence $h_{q}$ can be obtained analytically.

Here we assume that $X$ is a hyperexponential Lévy process and $q>0$. The following results are immediate by straightforward integration.

Proposition 4.1**.**

Suppose $B=(\underline{b},\overline{b})$ and $A=(-\overline{a},-\underline{a})$ for some $0\leq\underline{a}\leq\overline{a}$ and $0\leq\underline{b}\leq\overline{b}$. Then

[TABLE]

where, for each $1\leq j\leq m$,

[TABLE] 2. 2.

We have

[TABLE]

In particular, by setting $B=(0,\infty)$ ($A=(-\infty,0)$),

[TABLE]

where

[TABLE]

5 Numerical Results

Using the identities obtained in Proposition 4.1, we shall evaluate the efficiency of the phase-type fitting approach of the Gerber-Shiu function.

Here we consider the spectrally negative Lévy processes $X^{(\textrm{weibull})}$ and $X^{(\textrm{pareto})}$ in the form (3.1) where $Z$ is

1. (i)

Weibull$(0.6,0.665)$, and 2. (ii)

Pareto$(1.2,5)$

respectively. Recall that the Weibull distribution with parameters $c$ and $a$ (denoted Weibull$(c,a)$) is give by $F(t)=1-e^{-(t/a)^{c}}$, $t\geq 0$, and the Pareto distribution with positive parameters $a$ and $b$ (denoted Pareto$(a,b)$) is given by $F(t)=1-(1+bt)^{-a}$, $t\geq 0$. See [16] for more details about these distributions. With the choice of our parameters, the corresponding Lévy densities are thus completely monotone.

As has been noted in Section 3.4, any spectrally negative Lévy process with a completely monotone Lévy density can be approximated arbitrarily closely by fitting hyperexponential distributions. Here, we use the fitted data computed by [14] to approximate the scale function for $X^{(\textrm{weibull})}$ and $X^{(\textrm{pareto})}$ (with or without a Brownian motion component). Tables 3 and 9, respectively, of [14] show the parameters of the hyperexponential distributions obtained by [14] fitted to (i) with $m=6$ and to (ii) with $m=14$. We use these parameters to construct hyperexponential Lévy processes $\widetilde{X}^{(\textrm{weibull})}$ and $\widetilde{X}^{(\textrm{pareto})}$ (see Section 3.4) that will be used to approximate $X^{(\textrm{weibull})}$ and $X^{(\textrm{pareto})}$, respectively.

In order to evaluate the errors associated with the phase-type-fitting of the Gerber-Shiu function, we shall consider approximating the overshoot/undershoot density

[TABLE]

for $X^{(\textrm{weibull})}$ and $X^{(\textrm{pareto})}$. The phase-type-fitting approach approximates them by computing those for the approximating hyperexponential Lévy processes $\widetilde{X}^{(\textrm{weibull})}$ and $\widetilde{X}^{(\textrm{pareto})}$. These can be done analytically by the identity (4.2).

We evaluate these results by comparing with the simulated results. For $X^{(\textrm{weibull})}$ and $X^{(\textrm{pareto})}$, we simulate

[TABLE]

with $\Delta a=\Delta b=0.1$ by Monte Carlo simulation with $500,000$ samples.

In order to confirm the accuracy of the simulated results, we also compute the results using these two methods for $X^{(\textrm{exp})}$ corresponding to the process where the jump size is exponential with parameter $1$; this is a special case of the phase-type (and hyperexponential) Lévy process and hence obtained results are exact.

Figures 2 and 3 show the results for the cases $\sigma=1$ and $\sigma=0$ with common parameters $x=5$, $\mu=1$, $\lambda=10$ and $q=0.05$. From the results on the exponential case, we can confirm that the Monte Carlo simulated results are accurate. In Figure 3, the density has a jump at the initial position $x=5$ for the case $\sigma=0$ (while it is continuous for the case $\sigma=1$) due to the fact that $x=5$ is irregular for $(-\infty,5)$ (see Definition 6.4 of [19]). As can be seen from these figures, the approximation accurately captures the overshoot/undershoot densities for $X^{(\textrm{weibull})}$ and $X^{(\textrm{pareto})}$. The spike at the initial position in Figure 3 is precisely realized thanks to the closed-form expression (4.2); this would be difficult to realize if the scale function is approximated via numerical Laplace inversion.

6 Concluding Remarks

In this paper, we studied and evaluated the performance of the phase-type fitting approach (motivated by [13]) for computing the Gerber-Shiu function for the spectrally negative Lévy process. The method is overall accurate and is powerful in that it can obtain an approximation in a closed form.

While this paper focused on the case of a spectrally negative Lévy process, the proposed method can easily be generalized. By the compensation formula, the decomposition (2.4) holds for a more general class of stochastic processes. Therefore, by simply replacing the resolvent measure $r^{(q)}$, the Gerber-Shiu function can be computed in the same way.

Thanks to the recent developments of the fluctuation theories, the resolvents of various extensions of the spectrally negative Lévy process are now available in terms of the scale function. Here, we list several known examples of the resolvent.

The fluctuation theories of reflected spectrally negative Lévy processes are well-developed. In the optimal dividend problem, where one wants to maximize the expected net present value of dividends until ruin, it is in many cases shown to be optimal to reflect the surplus process at a suitable boundary. It is therefore of interest to investigate the Gerber-Shiu function of the reflected process to evaluate the risk of an dividend paying company. As given in [31], the resolvent requires the derivative or the integral of the scale function, depending on whether the reflection barrier is upper or lower. For the doubly reflected case with both upper and lower barriers, see [30]. 2. 2.

As a variant of the reflected process, the refracted spectrally negative Lévy process of [21] changes its drift by $\delta>0$ whenever it is above a threshold $b$ – it is the unique strong solution $U$ to the stochastic differential equation

[TABLE]

In insurance, this can be used to model the surplus of a dividend paying company when the dividend rate must be bounded from above by $\delta$ (see [22]). The resolvent is given in [21]. The resolvents for the cases with additional classical reflection have recently been obtained by [28, 29]. 3. 3.

Given two levels $s$ and $S$, the classical $(s,S)$-policy controls the process by pushing the process to $S$ immediately when it goes above or below $s$. It is, under a suitable condition, an optimal strategy in the optimal dividend problem in the presence of a fixed cost (see [7, 25]). The resolvent is obtained in [35]. For its two-sided cases with four-parameter $(d,D,U,u)$ policy, see [34].

In these examples, the resolvent $r^{(q)}$ is written in terms of the scale function, and hence as in the cases considered in this paper, those for the phase-type case can be written as (a linear combination of) exponential forms. Consequently, the Gerber-Shiu function can be analytically computed in the same way. In view of the results obtained in this paper, the phase-type error is minimal and the same procedure is expected to give an accurate approximation of the Gerber-Shiu function.