**1. Introduction**

Let us first recall the Cramér–Lundberg risk model extended with Brownian perturbation Albrecher and Asmussen (2010); Dufresne and Gerber (1991b)

$$X\_t = x + ct + \sigma B(t) - S\_{t\prime} S\_t = \sum\_{i=1}^{N\_h(t)} \mathbb{C}\_i. \tag{1}$$

Here *x* ≥ 0 is the initial surplus, *c* ≥ 0 is the linear premium rate. The *Ci*'s, *i* = 1, 2, ... are independent identically distributed (i.i.d) random variables with distribution *<sup>F</sup>*(*z*) = *FC*(*z*) representing nonnegative jumps arriving after independent exponentially distributed times with mean 1/*λ*, and *<sup>N</sup>λ*(*t*) denotes the associated Poisson process counting the arrivals of claims on the interval [0, *t*]. Finally, *<sup>σ</sup><sup>B</sup>*(*t*), *σ* > 0 is an independent Brownian perturbation. Ruin happens when, for the first time, a jump takes *Xt* below 0.

Risk theory revolved initially around evaluating and minimizing the probability of ruin. Insurance companies are also interested in maximizing company value. This lead to the study of optimal dividend policies. As suggested by de Finetti in the 1950s de Finetti (1957)—see also Miller and Modigliani (1961)—an interesting objective is that of maximizing the expected value of the sum of discounted future dividend payments until the time of ruin.

The most important class of dividend policies is that of a constant barrier at *b*, which modifies the surplus only when *Xt* > *b*, by a lump paymen<sup>t</sup> bringing the surplus to *b*, and then keeps it there by Skorokhod reflection until the next negative jump. In financial terms, in the absence of a Brownian component, this amounts to paying out all the income while at *b*. In the case of Brownian perturbation, Skorokhod reflection means keeping the process above the barrier by minimal capital injections (whenever necessary), or below a barrier, by taking out dividends (if necessary) Skorokhod (1962).

In the presence of the barrier at *b*, the de Finetti objective (the expected value of the sum of discounted future dividend payments until ruin) has a simple expression Avram et al. (2007) in terms of the so-called "scale function" *W*:

$$V^{b^{\left[b\right]}}(\mathbf{x}) = \mathbb{E}\_{\mathbf{x}}^{b^{\left[b\right]}} \left[ \int\_{\left[0, T\_0^{b\right]} \right]} e^{-qt} d\mathcal{U}\_t^{b} \bigg|\_{} = \begin{cases} \frac{W\_q(\mathbf{x})}{W\_q(\mathbf{b})}, & \mathbf{x} \le \mathbf{b} \\ \mathbf{x} - \mathbf{b} + \frac{W\_q(\mathbf{b})}{W\_q(\mathbf{b})}, & \mathbf{x} > \mathbf{b}' \end{cases} \tag{2}$$

where *Tb*] 0 is the time of ruin, *q* denotes the discount rate, *Ubt* the total local time at *b* before time *t*, and E*b*] the law of the process reflected from above at *b* and absorbed at 0 and below.

The scale function Bertoin (1998); Kyprianou (2014); Landriault and Willmot (2019); Suprun (1976) *Wq*(*x*) : **R** → [0, <sup>∞</sup>), *q* ≥ 0 is defined on the positive half-line by the Laplace transform

$$
\hat{\mathcal{W}}\_{\emptyset}(\mathbf{s}) := \int\_0^\infty \mathbf{e}^{-s\mathbf{x}} \mathcal{W}\_{\emptyset}(\mathbf{x}) d\mathbf{x} = \frac{1}{\kappa(s) - q'} \quad \forall s > \Phi\_{\emptyset} \tag{3}
$$

where the "symbol" *κ*(*s*) (also-called the cumulant generating function) is defined in Equation (8) in Section 2 where we provide the necessary background information, and <sup>Φ</sup>*q* is the unique nonnegative root of the Cramér–Lundberg equation

$$\Phi\_{\emptyset} := \sup \{ s \ge 0 : \kappa(s) - q = 0 \}, \quad q \ge 0. \tag{4}$$

The scale function *Wq*(*x*) is continuous and increasing on [0, ∞) Bingham (1976), (Bertoin 1998, Thm. VII.8), (Kyprianou 2014, Thm. 8.1). It may have, however, many inflection points (such an example is depicted in Figure 1), and these play an important role in the optimization of dividends Avram et al. (2007, 2015); Schmidli (2007). For convenience, *Wq*(*x*) is extended to be 0 on **R**−. An important fact that will be exploited is that the Laplace transform of our function has a unique non-negative pole <sup>Φ</sup>*q*, see Equations (3) and (4).

This paper aims at computing/approximating the scale function *Wq*(*x*), using its moments. The techniques being used are classic: Padé approximation and Laguerre expansions. The order (*<sup>m</sup>*, *n*) Padé approximation of a function *g*(*x*) is a rational function in the form

$$R(\mathbf{x}) = \frac{a\_0 + a\_1 \mathbf{x} + a\_2 \mathbf{x}^2 + \dots + a\_m \mathbf{x}^m}{1 + b\_1 \mathbf{x} + b\_2 \mathbf{x}^2 + \dots + b\_n \mathbf{x}^n}$$

for which *R*(0) = *g*(0), *R*(0) = *g*(0), *R*(0) = *g*(0), ... , *R*(*m*+*n*)(0) = *<sup>g</sup>*(*m*+*n*)(0). In the context of probability distributions, given a density function *f*(*x*) and its Laplace transform *f* (*s*), the inverse Laplace transform of the order (*<sup>m</sup>*, *m*) Padé approximant of *f* (*s*) provides a matrix exponential approximation of *f*(*x*) that matches the first 2*m* moments of *f*(*x*) (including *m*0). In Avram et al. (2011) this approach was used to approximate ruin probabilities. In this paper we develop the same approach to approximate the scale function *Wq*(*x*) (Section 3). An extension of the above idea is the so-called two-point Padé approximation, which allows to match not only the moments of *Wq*(*x*) but also the behavior of the function at 0, i.e., to match *Wq*(0), *<sup>W</sup>q*(0), ... (Section 4). For more details on this extension see Avram et al. (2018) where ruin probabilities are approximated.

Let us draw attention now to several numeric challenges which were absent in the ruin probability problem.

1. Optimizing dividends starts by optimizing the so-called "barrier function"

$$H\_D(b) := \frac{1}{\mathcal{W}'\_q(b)'} \ b \ge 0,\tag{5}$$

and the optimal dividend policy is often simply a barrier strategy at its maximum. This is the case in particular when the barrier function *HD*(*b*) is differentiable with

$$\mathcal{H}\_D'(0) > 0 \Leftrightarrow \mathcal{W}\_q''(0) < 0\tag{6}$$

and has a unique local maximum *b*∗ > 0 =⇒ *W q* (*b*∗) = 0; then this *b*∗ yields the optimal dividend policy, and the optimal barrier function,

$$V(\mathbf{x}) := \sup p\_{b \ge 0} V^{b \restriction}(\mathbf{x}) = V^{b^\* \restriction}(\mathbf{x}),\tag{7}$$

turns out to be the largest concave minorant of *Wq*(*x*).<sup>1</sup>

2. The challenge of multiple inflection points. In the presence of several inflection points, however, the optimal policy is multiband Azcue and Muler (2005); Schmidli (2007); Avram et al. (2015); Loeffen (2008). The first numerical examples of multiband policies were produced in Azcue and Muler (2005); Loeffen (2008), with Erlang claims *Erl*2,1. However, it was shown in Loeffen (2008) that multibands cannot occur when *<sup>W</sup>q*(*x*) is increasing after its last global minimum *b*∗ (i.e., when no local minima are allowed after the global minimum).

Loeffen (2008) further made the interesting observation that for Erlang claims *ER*2,1 (which are non-monotone), multiband policies may occur for volatility parameters *σ* smaller than a threshold value, but barrier policies (with a non-concave value function) will occur when *σ* is large enough.

Figure 1 displays the first derivative *<sup>W</sup>q*(*x*), for *σ*2/2 ∈ { 12 , 1, 32 , <sup>2</sup>}. The last two values yield barrier policies with a non-concave value function, due to the presence of an inflection point in the interior of the interval [0, *b*∗].

1 10, *σ*2/2 ∈ {1/2, 1, 3/2, <sup>2</sup>}.

**Figure 1.**

1

Even when barrier strategies do not achieve the optimum, and multi-band policies must be used instead, constructing the solution must start by determining the global maximum of the barrier function Avram et al. (2015); Azcue and Muler (2005); Schmidli (2007).

Below we will investigate whether our approximations are precise enough to yield reasonable approximations for *W q* (0) and the root(s) of *W q* (·).

Special features. While our methods consist essentially of Padé approximation and Laguerre–Tricomi–Weeks Laplace transform inversion, we found that exploiting the special features of our problem is useful. These are:


Contents. We briefly review classical ruin theory in Section 2. Padé approximations are provided in Section 3, where we also spell out the simplest algorithm for the computation of the scale function. In Section 4, we derive low-order Padé and two-point Padé approximations of *Wq*(*x*), reminiscent of the de Vylder approximation of the ruin probability. Some of these approximations appeared already in Gerber et al. (2008), where however the Padé method and the fact that they can be easily extended to higher orders is not mentioned. Section 5 offers our personal strategy for inverting Laplace transforms of interest in probability, in the presence of uncertainty. Section 5.1 implements the Laguerre–Tricomi–Weeks Laplace transform inversion with a certain judicious choice of the exponential decay parameter (40), which is believed to be new. Section 6 presents numeric experiments with mixed exponential claims. Section 7 presents experiments with Pareto claims; since these have a heavy tail and, consequently, finitely many moments, we apply a "shifted" Padé approximation of the claim distribution. Section 8 includes a computer program required to obtain test cases with exact rational answers, using the Wiener–Hopf factorization; of course, this is quite convenient for the initial testing of the precision of our algorithms. Finally, Section 9 reviews a more general version of the Laguerre–Tricomi–Weeks Laplace inversion method, which may be of interest for further experiments.

### **2. A Short Review of Classical Ruin Theory**

The process defined in Equation (1) is a particular example of spectrally negative Lévy processes, with finite mean, which are defined by assuming instead of Equation (1) that *St* is a subordinator with *σ*-finite Lévy measure *ν*(*dx*) that integrates *x*, but having possibly infinite activity near the origin *ν*(0, ∞) = ∞ Bertoin (1998) (for Equation (1), the Lévy measure is given by *ν*(*dx*) = *λF*(*dx*) =⇒ *ν*(0, ∞) = *λ*). A spectrally negative Lévy process is characterized by its Lévy–Khintchine/Laplace exponent/cumulant generating function/symbol *κ*(*s*) defined by *<sup>E</sup>*0*<sup>e</sup>sX*(*t*) = *<sup>e</sup>tκ*(*s*), with *κ*(*s*) of the particular form

$$\kappa(\mathbf{s}) = c\mathbf{s} + \int\_0^\infty (e^{-s\mathbf{x}} - 1)\nu(d\mathbf{x}) + \frac{\sigma^2 \mathbf{s}^2}{2} \coloneqq s\left(c - \widehat{\overline{\nu}}(\mathbf{s}) + \frac{\sigma^2}{2}\mathbf{s}\right). \tag{8}$$

Some concepts of interest in classical risk theory are:

• First passage times below and above a level *a*

$$T\_{a,-(+)} := \inf \{ t > 0 : X(t) < (>)a \}.$$

• The first first passage quantity to be studied historically was the eventual ruin probability:

$$\Psi(\mathbf{x}) := P[T\_{0,-} < \infty | X(0) = \mathbf{x}].\tag{9}$$

In order that the eventual ruin probability not be identically 1, the parameter

$$p := c - \lambda m\_1 = \kappa'(0), \text{ where } m\_1 = \int\_0^\infty z F(dz),$$

which is called drift or profit rate, must be assumed positive.

The Laplace transform of the ruin probability is explicit, given by the so-called Pollaczek–Khinchine formula, which states that the Laplace transform of <sup>Ψ</sup>(*x*) = 1 − <sup>Ψ</sup>(*x*) is:

$$
\hat{\Psi}(s) = \int\_0^\infty e^{-sx} \overline{\Psi}(x) dx = \frac{c - \lambda m\_1}{\kappa(s)} = \frac{\kappa'(0)}{\kappa(s)}.\tag{10}
$$

The roots with negative real part of the Cramér–Lundberg equation

$$
\kappa(s) = q, q \ge 0 \tag{11}
$$

are important, when such roots exist. They will be denoted by −*γ*1, −*γ*2, ··· , <sup>−</sup>*γN*, *N* ≥ 0, and ordered by their absolute values |*<sup>γ</sup>*1|≤|*<sup>γ</sup>*2| ≤ ... , ≤ |*<sup>γ</sup>N*|. *γ*1 > 0 is called the adjustment coefficient, and furnishes the Cramér–Lundberg asymptotic approximation

$$\Psi(u) \sim \frac{\kappa'(0)}{-\kappa'(-\gamma\_1)} e^{-\gamma\_1 u}.$$

**Laplace transform inversion**. As explained here, the first passage theory for Lévy processes with one-sided jumps reduces essentially to inverting the Laplace transform Equation (3). This applies not only to the well-known ruin probabilities, but also to intricate optimization problems involving dividends, capital gains, liquidation of subsidiary companies, etc.

**Remark 1.** *If the claims have a phase-type or matrix exponential distribution, then the Cramér–Lundberg equation has a finite number of roots, and Laplace transform inversion reduces to finding roots of denominators and to partial fractions, operations which are essentially exact with the current computing resources. For example, the ruin probability is provided by an exact formula (which extends the Cramér–Lundberg asymptotic approximation). With distinct roots, this is*

$$\Psi(\mu) = \sum\_{i=1}^{N} \frac{\kappa'(0)}{-\kappa'(-\gamma\_i)} e^{-\gamma\_i \mu}. \tag{12}$$

*Similar formulas hold for the scale function and related quantities (derivatives and integrals of the scale function, etc).*

*This refocuses the question of ruin probabilities and similar quantities to the harder cases of*


For non-matrix exponential jumps, one can start by a Padé approximation of the Laplace transform, which is perhaps the oldest method of Laplace transform inversion. For heavy-tailed claims, whose moment generating function is not analytic at the origin, one may apply, as it will be shown in Section 7, a "shifted" Padé approximation before Laplace transform inversion—see, for example, Avram et al. (2018). In this paper, we also compare the precision of the classic Padé and two-point Padé Laplace transform inversion methods with the Laguerre–Tricomi–Weeks inversion, as applied to the optimal dividends barrier problem. As test cases, we consider mixed exponentials, for which exact answers are available for the case of rational Wiener–Hopf factorization roots described in Section 8.

### **3. Padé/matrix exponential Approximations of the Scale Function**

The essence of classic ruin theory is the availability of explicit "output Laplace transforms" expressed in terms of "input Laplace transforms", Equations see (3), (10) and (16). Equivalently, "output moments" (i.e., Taylor coefficients of functions of interest) are expressed in terms of "input data moments". In the case of insufficient data, Padé approximations of the output function seem to deserve special attention, because they involve only a few input moments, and seem to extract at low orders the essence of the data. Note, as pointed out in Avram et al. (2011, 2018); Avram and Pistorius (2014), that the classical ruin theory approximations of Cramér–Lundberg, De Vylder and Renyi are all first order Padé approximations. Slightly increasing the order yields more sophisticated moments based approximations Avram et al. (2018); Avram and Pistorius (2014); Ramsay (1992). Similar approximations, which could be useful in dividends optimization, are presented in Section 4—see also Hu et al. (2017) for a related application to reinsurance.

With more reliable data, higher order Padé approximations are just as easy to obtain, using computer systems such as Maple, Mathematica, Sage, Matlab, etc. Let *νk* = *λmk*, *k* ≥ 1 denote the moments of the Lévy measure2, and let

$$
\nu\_{2,\sigma} = \nu\_2 + \sigma^2. \tag{13}
$$

The simplest Padé approximation of the scale function and optimal dividend barrier requires implementing the following algorithm, which is valid for claim distributions having 2*n* moments.

1. Obtain the power series expansion of the Laplace exponent in terms of the moments of the Lévy measure

$$\kappa(s) = s\left(c - \hat{\overline{\nu}}(s) + \frac{\sigma^2}{2}s\right) = \left(c - \nu\_1\right)s + \nu\_{2,\sigma}\frac{s^2}{2} + \sum\_{k=3}^{\infty} \nu\_k \frac{(-s)^k}{k!}.\tag{14}$$

2. Construct a Padé approximation

$$
\hat{\mathcal{W}}\_q(\mathbf{s}) = \frac{1}{(\mathbf{c} - \nu\_1)\mathbf{s} + \nu\_{2,\sigma}\mathbf{s}^2/2 - \nu\_3\mathbf{s}^3/6 + \dots - q} \sim \frac{P\_{n-1}(\mathbf{s})}{Q\_n(\mathbf{s})} = \frac{\sum\_{i=0}^{n-1} a\_i \mathbf{s}^i}{\mathbf{c}\mathbf{s}^n + \sum\_{i=0}^{n-1} b\_i \mathbf{s}^i},
$$

i.e., find *ai*, *bi*, *i* = 0, ... , *n* − 1, by solving the linear system *Pn*−<sup>1</sup>(*s*)(*κ*(*s*) − *q*) − *Qn*(*s*) = *<sup>O</sup>*(*sn*+<sup>1</sup>). As emphasized in Figure 2, the series expansion of *κ*(*s*) is a good approximation only for small *s*; since the most important part of the approximation is the root <sup>Φ</sup>*q* in Equation (25), simple Padé approximation will only work for small *q*.

<sup>2</sup> Note that only moments starting from 1 are required, so this may be applied to processes whose subordinator part has infinite activity as well.

**Figure 2.** The series expansion of *κ*(*s*) has multiple contact with *κ*(*s*) at *s* = 0, but increases asymptotically at a smaller rate. Despite that, it yields quite reasonable approximations of <sup>Φ</sup>*q*, for *q* small enough.

3. Factor the denominator as

$$\mathbf{c}s^n + \sum\_{i=0}^{n-1} b\_i \mathbf{s}^i = \mathbf{c}(\mathbf{s} - \Phi\_{\overline{q}}) \prod\_{i=1}^n (\mathbf{s} + \gamma\_i). \tag{15}$$

This operation is the essential numerical difficulty, but may be achieved nowadays with arbitrarily high precision.

4. Then, partial fractions plus inversion quickly yield an approximate Laplace transform inverse

$$\mathcal{W}\_{\emptyset}(\mathbf{x}) \sim \mathbb{C}e^{\Phi\_{\emptyset}\mathbf{x}} + \sum\_{i=1}^{n} \mathbb{C}\_{i} e^{-\gamma\_{i}\mathbf{x}}.$$

5. The dividend barrier is obtained by finding the largest nonnegative root of *<sup>W</sup>*(*x*) = 0.

**Remark 2.** *Comparing Equation* (15) *with the fundamental relation of Equation* (44)*, we see that the Padé approximation of the output Laplace transform may be viewed as an approximate Wiener–Hopf factorization.*

**Remark 3.** *Note that, as already mentioned in Remark 1, Padé inversion cannot be improved upon when the density of the claims is matrix exponential of some order n, (or, equivalently, when its Laplace transform is rational), since it results into the exact decomposition of the scale function into exponentials*

$$W\_{\emptyset}(\mathbf{x}) = \Phi\_{\emptyset}^{\prime} e^{\Phi\_{\emptyset} \mathbf{x}} + \sum\_{i=1}^{n} \mathbb{C}\_{i} e^{-\gamma\_{i} \mathbf{x}} \lambda\_{i}$$

*where the Ci are partial fraction decomposition coefficients (and where the first coefficient* Φ*q follows from the normalization constant chosen in Equation* (3)*).*

*The coefficients Ci are less important than the roots, and in fact one may easily improve on Padé, for example, by recomputing them so that the approximation is optimal in L*2 *sense. Positivity of the Ci's, or, better, of their "Coxian linear combinations" can be further ensured by linear programming, if desired.*

The preceding observations make Laplace inversion by rational approximation of the Laplace transform a quite attractive tool in ruin theory and related fields, from both a pedagogical and numerical point of view Avram et al. (2011, 2018); Avram and Pistorius (2014). One drawback is the appearance of non-admissible roots (for example, with (*<sup>γ</sup>i*) < 0), which will imply non-admissible ruin probabilities for very large initial reserves *x*. It may be argued, however, that if the admissibility range is large enough, in practice this may not cause problems. Moreover, this problem may be fixed by convenient mathematical tools, implemented in packages such as BUTools and SOPE Bobbio et al. (2005); Dumitrescu et al. (2016); Horváth and Telek (2016).

For the sake of simplicity, we will mostly take *σ* = 0 from now on. In this case, the Pollaczek–Khinchine formula, Equation (10), may be simply expressed in terms of the so-called "excess distribution *fe*" of the jumps, and of the "intensity" *ρ*:

$$\begin{split} \hat{\overline{\Psi}}(s) &= \frac{c(1-\rho)}{\kappa(s)} = \frac{1-\rho}{s(1-\rho\hat{f}\_{\varepsilon}(s))} \Leftrightarrow\\ \hat{\Psi}(s) &= \int\_{0}^{\infty} e^{-sx} \Psi(x) dx = \frac{1}{s} \left( 1 - \frac{1-\rho}{1-\rho\hat{f}\_{\varepsilon}(s)} \right) = \rho \frac{\hat{\overline{F}}\_{\varepsilon}(s)}{1-\rho\hat{f}\_{\varepsilon}(s)} = \rho \frac{\hat{\overline{F}}\_{\varepsilon}(s)}{1-\rho\frac{\hat{\overline{F}}(s)}{m\_{1}}}, \end{split} \tag{16}$$

where *f* (*s*) = ∞ 0 *<sup>e</sup>*<sup>−</sup>*sxdF*(*x*), *<sup>F</sup>*(*x*) = 1 − *<sup>F</sup>*(*x*), *m*1 = ∞ 0 *xdF*(*x*), *fe*(*x*) = *<sup>F</sup>*(*x*)/*<sup>m</sup>*1, *f e*(*s*) = ∞ 0 *e*<sup>−</sup>*sx fe*(*x*)*dx* = (1 − *f* (*s*))/(*<sup>m</sup>*1*<sup>s</sup>*), *ρ* = *λ <sup>m</sup>*1/*<sup>c</sup>*, *Fe*(*s*)=(<sup>1</sup> − *f <sup>e</sup>*(*s*))/*<sup>s</sup>*. The last formula is especially convenient for cases where the excess distribution is simply related to the original one, as in the case of the Pareto distribution.

We would like to note however that the Padé approach works with no problems in the presence of Brownian motion, the only modifications being the form of *κ*(*s*) and the second moment *ν*2 of the associated Lévy measure—see Equation (13).

### **4. Two-Point Padé Approximations, with Low Order Examples**

One may obtain better results by incorporating into the Padé approximation the following initial values, which can be derived easily from the Laplace transform:

$$\mathcal{W}\_q(0) = \lim\_{s \to \infty} s \hat{\mathcal{W}}\_q(s) = \begin{cases} \frac{1}{c'} & \text{if } X \text{ is of bounded variation/Cramér-Lundberg} \\ 0, & \text{if } X \text{ is of unbounded variation} \end{cases},\tag{17}$$

$$\mathcal{W}'\_q(0) = \lim\_{s \to \infty} s \left( \frac{s}{\kappa(s) - q} - \mathcal{W}\_q(0) \right) = \begin{cases} \frac{q + \nu(0, \infty)}{\epsilon^2}, & \text{if } X \text{ is of bounded variation} \\\frac{2}{\sigma^2}, & \text{if } \sigma > 0, \\\infty, & \text{if } \sigma = 0 \text{ and } \nu(0, \infty) = \infty \end{cases} \tag{18}$$

Furthermore, when the jump distribution has a density *f* , it holds that :3

$$\begin{split} \mathcal{W}\_{q}^{\prime\prime}(0\_{+}) &= \lim\_{s \to \infty} \mathrm{s} \left( \mathrm{s} \left( \frac{\mathrm{s}}{\kappa(s) - q} - \mathcal{W}\_{q}(0) \right) - \mathcal{W}\_{q}^{\prime}(0\_{+}) \right) \\ &= \begin{cases} \frac{1}{\varepsilon} \Big( (\frac{\lambda + q}{\varepsilon})^{2} - \frac{\lambda}{\varepsilon} f(0) \Big), & \text{if } X \text{ is of bounded variation} \\ -\varepsilon (\frac{2}{\sigma^{2}})^{2}, & \text{if } \sigma > 0 \end{cases} . \end{split} \tag{19}$$

Further derivatives at 0 could be computed, but we stop at order 2, since *W q* (0) already requires estimating *fC*(0), which is a rather delicate task starting from real data.

We will provide in Proposition 1 below a couple of two-point Padé approximations, when *n* = 2. Before that, it is worth recalling the case of exponential claims.

**Example 1.** *The Cramér–Lundberg model with exponential jumps. Consider the Cramér–Lundberg model with exponential jump sizes with mean* 1/*μ, jump rate λ, premium rate c* > 0*, and Laplace exponent κ*(*s*) = *s c* − *λ μ*+*s , assuming κ* (0) = *c* − *λ μ* > 0*. Let γ* = *μ* − *λ*/*c denote the adjustment coefficient, and let ρ* = *λ cμ .*

<sup>3</sup> This equation is important in establishing the nonnegativity of the optimal dividends barrier.

*Solving κ*(*s*) − *q* = 0 ⇔ *cs*<sup>2</sup> + *s*(*cμ* − *λ* − *q*) − *qμ* = 0 *for s yields two distinct solutions γ*2 ≤ 0 ≤ *γ*1 = <sup>Φ</sup>*q given by*

$$\begin{aligned} \gamma\_1 &= \frac{1}{2c} \left( -\left( \mu c - \lambda - q \right) + \sqrt{\left( \mu c - \lambda - q \right)^2 + 4\mu qc} \right), \\\gamma\_2 &= \frac{1}{2c} \left( -\left( \mu c - \lambda - q \right) - \sqrt{\left( \mu c - \lambda - q \right)^2 + 4\mu qc} \right). \end{aligned}$$

*The W scale function is:*

$$\mathcal{W}\_{\boldsymbol{q}}(\mathbf{x}) = \frac{A\_1 \boldsymbol{\varepsilon}^{\gamma\_1 \mathbf{x}} - A\_2 \boldsymbol{\varepsilon}^{\gamma\_2 \mathbf{x}}}{\boldsymbol{c}(\gamma\_1 - \gamma\_2)} \Leftrightarrow \hat{\mathcal{W}}\_{\boldsymbol{q}}(\mathbf{s}) = \frac{\mathbf{s} + \boldsymbol{\mu}}{\boldsymbol{c} \mathbf{s}^2 + \boldsymbol{s}(\boldsymbol{c} \boldsymbol{\mu} - \boldsymbol{\lambda} - \boldsymbol{q}) - q\boldsymbol{\mu}} \tag{20}$$

*where A*1 = *μ* + *γ*1, *A*2 = *μ* + *γ*2*.*

*Furthermore, it is well-known and easy to check that the function <sup>W</sup>q*(*x*) = *HD*(*x*)−<sup>1</sup> *is in this case unimodal with global minimum at*

$$b^\* = \frac{1}{\gamma\_1 - \gamma\_2} \begin{cases} \log \frac{(\gamma\_2)^2 A\_2}{(\gamma\_1)^2 A\_1} = \log \frac{(\gamma\_2)^2 (\mu + \gamma\_2)}{(\gamma\_1)^2 (\mu + \gamma\_1)} & \text{if } W\_q''(0) < 0 \Leftrightarrow (q + \lambda)^2 - c\lambda\mu < 0\\ 0 & \text{if } W\_q'''(0) \ge 0 \Leftrightarrow (q + \lambda)^2 - c\lambda\mu \ge 0 \end{cases} \tag{21}$$

*since W q* (0) = (*<sup>γ</sup>*1)<sup>2</sup>(*μ*+*γ*1)−(*<sup>γ</sup>*2)<sup>2</sup>(*μ*+*γ*2) *<sup>c</sup>*(*<sup>γ</sup>*1−*γ*2) ∼ (*q* + *λ*)<sup>2</sup> − *cλμ and that the optimal strategy for the de Finetti problem is the barrier strategy at level b*∗ *Avram et al. (2007).*

**Proposition 1.** *1. To secure both the values of Wq*(0) *and <sup>W</sup>q*(0)*, take into account Equations* (17) *and* (18)*, i.e., use the Padé approximation*

$$
\hat{\mathcal{W}}\_{\emptyset}(\mathbf{s}) \sim \frac{\sum\_{i=0}^{n-1} a\_i \mathbf{s}^i}{\mathbf{c} \mathbf{s}^n + \sum\_{i=0}^{n-1} b\_i \mathbf{s}^i}, a\_{n-1} = 1, b\_{n-1} = c a\_{n-2} - \lambda - q.
$$

*This yields*

$$
\hat{\mathcal{W}}\_{\boldsymbol{q}}(\boldsymbol{s}) \sim \frac{\frac{1}{m\_1} + \boldsymbol{s}}{c s^2 + s \left(\frac{\boldsymbol{c}}{m\_1} - \lambda - \boldsymbol{q}\right) - \frac{\boldsymbol{q}}{m\_1}}.\tag{22}
$$

*In view of Equation* (20)*, this yields the same result as approximating the claims by exponential claims, with μ* = 1*m*1*.*

*2. To ensure Wq*(0) = 1*c , we must only impose the behavior specified in Equation* (17)*, i.e., use the Padé approximation*

$$
\widehat{\mathcal{W}}\_{\mathfrak{q}}(\mathbf{s}) \sim \frac{\sum\_{i=0}^{n-1} a\_i \mathbf{s}^i}{\mathrm{cs}^n + \sum\_{i=0}^{n-1} b\_i \mathbf{s}^i}, a\_{n-1} = 1.
$$

*For n* = 2*, this yields*

$$\widehat{\mathcal{W}}\_{q}(s) \sim \frac{\frac{2w\_{1}}{m\_{2}} + s}{cs^{2} + \frac{s\left(2cm\_{1} - 2\lambda m\_{1}^{2} - m\_{2}q\right)}{m\_{2}} - \frac{2w\_{1}q}{m\_{2}}} = \frac{\frac{1}{\overline{m}\_{1}} + s}{cs^{2} + s\left(\frac{c}{\overline{m}\_{1}} - \lambda\frac{m\_{1}}{\overline{m}\_{1}} - q\right) - \frac{q}{\overline{m}\_{1}}},\tag{23}$$

*where we denoted by m*˜ 1 = *m*2 2*m*1 *the first moment of the excess density fe*(*x*)*. For exponential claims this coincides with Equation* (22) *(since fe*(*x*) = *f*(*x*)*). This is the De Vylder B) method (Gerber et al. 2008, (5.6–5.7)), derived therein by assuming exponential claims, with μ* = 2*m*1 *m*2 *, and simultaneously modifying λ to fit the first two moments of the risk process.*

### *3. When the pure Padé approximation is applied, the first step yields*

$$\widehat{\mathcal{W}}\_{q}(\mathbf{s}) \sim \frac{\mathbf{s} + \frac{3m\_{2}}{m\_{3}}}{\mathbf{s}^{2}\left(\mathbf{c} + \lambda\left(\frac{3m\_{2}}{2m\_{3}} - m\_{1}\right)\right) + \mathbf{s}\left(\mathbf{c}\frac{3m\_{2}}{m\_{3}} - \frac{3m\_{1}m\_{2}}{m\_{3}}\lambda - q\right) - \frac{3m\_{2}}{m\_{3}}q}$$

$$\mathbf{s} = \frac{\mathbf{s} + \frac{1}{\tilde{m}\_{3}}}{\mathbf{s}^{2}\left(\mathbf{c} + \lambda\left(\frac{\tilde{m}\_{2}}{\tilde{m}\_{3}} - 1\right)\right) + \mathbf{s}\left(\mathbf{c}\frac{1}{\tilde{m}\_{3}} - \frac{m\_{1}}{\tilde{m}\_{3}}\lambda - q\right) - \frac{1}{\tilde{m}\_{3}}q'}\tag{24}$$

*where <sup>m</sup>i* = *mi i mi*−1 *is a so-called "normalized moment" Bobbio et al. (2005).*

*This is the De Vylder A) method (Gerber et al. 2008, (5.2–5.4)), derived therein by assuming exponential claims, with μ* = 3*<sup>m</sup>*2 *m*3 *, and simultaneously modifying both λ and c to fit the first three moments of the risk process.*

**Lemma 1.** *In the case of exponential claims, the three approximations given above are exact.*

**Proof.** It suffices to check that for exponential claims all the normalized moments are equal to *m*1 = *μ* −1 .

In particular, the optimal barrier for exponential claims obtained by the explicit Equation (21) is the same. For example, *μ* = 2/5, *λ* = 9/10, *c* = 1, *q* = 1/10 yields

> *Wq*(*x*) ∼ 0.652989 2.718280.0659646*<sup>x</sup>* − 0.152989 2.71828−1.51596*<sup>x</sup>*

and *b*∗ = 3.04576.
