**7. Padé Approximation of Heavy-Tailed Claims**

In this section we experiment with heavy-tailed claims in the context of dividends barrier. To this end we consider hereby the risk process with Pareto II/Lomax claim size distribution (also known as American Pareto) with survival function given by *<sup>F</sup>*(*y*) = 1 + *yβ* −*α*, *y* > 0, *α*, *β* > 0, mean *m*1 = *βα*−1 and *Fe*(*y*) = 1 + *yβ* −(*<sup>α</sup>*−<sup>1</sup>), *y* > 0. The Laplace transform of the density is Nadarajah and Kotz (2006)

$$\hat{f}(\mathbf{s}) = \int\_0^\infty e^{-\mathbf{s}x} dF(\mathbf{x}) = a(\beta \mathbf{s})^a e^{\beta \mathbf{s}} \Gamma(-a, \beta \mathbf{s}) = a! I(1, 1 - a, \beta \mathbf{s}) = 1 - (\beta \mathbf{s})^a e^{\beta \mathbf{s}} \Gamma(1 - a, \beta \mathbf{s}), \mathbf{s} \ge 0$$

where <sup>Γ</sup>(*<sup>α</sup>*,*<sup>s</sup>*) = ∞*s t<sup>α</sup>*−1*e*<sup>−</sup>*tdt* is the incomplete gamma function and

$$\mathcal{U}[a, a+c, z] = \frac{1}{\Gamma[a]} \int\_0^\infty e^{-zt} t^{a-1} (t+1)^{c-1} dt, \text{Re}[z] > 0, \text{Re}[a] > 0$$

is Tricomi's Hypergeometric U function (Abramowitz and Stegun 1965, 13.2.5), where we used *U*1, 1 + *λ*˜ , *v*) = *<sup>e</sup>vv*−*λ*Γ(*<sup>λ</sup>*, *v*) (see, (Abramowitz and Stegun 1965, 6.5.22)). It follows that the Laplace transform of the survival function is

$$
\widehat{\vec{F}}(s) = \beta l I(1, 2 - \alpha, \beta s). \tag{43}
$$

As discussed in Albrecher et al. (2010), the model is suitable for obtaining precise quantities by numerical inverse Laplace transform using the fixed Talbot (FT) algorithm with parameter *M* = 200. We will refer to results obtained by the FT algorithm as "exact".

We apply *α* = 3/2, *β* = 1 with which the distribution is of infinite variance. One can still obtain moment based Padé approximations by "shifting" the claim distribution as described in Avram et al. (2018). The resulting approximate claim density function approximations of order *K* = 4, 6, and 8 are the following (the shift factor we applied is 1.6):

$$\begin{aligned} f\_{4}(\mathbf{x}) &= 0.3674e^{-4.8697\mathbf{x}} + 0.7439e^{-2.1699\mathbf{x}} + 0.3412e^{-0.8795\mathbf{x}} + 0.0465e^{-0.2883\mathbf{x}} \\ f\_{6}(\mathbf{x}) &= 0.1053e^{-6.6998\mathbf{x}} + 0.4679e^{-3.5877\mathbf{x}} + 0.5442e^{-1.9191\mathbf{x}} + 0.2898e^{-0.9714\mathbf{x}} + \\ &0.08279e^{-0.4298\mathbf{x}} + 0.0099e^{-0.1339\mathbf{x}} \\ f\_{8}(\mathbf{x}) &= 0.02944e^{-8.3899\mathbf{x}} + 0.2211e^{-4.9512\mathbf{x}} + 0.4433e^{-2.9899\mathbf{x}} + 0.424e^{-1.7882\mathbf{x}} + \\ &0.2512e^{-1.0331\mathbf{x}} + 0.1008e^{-0.555\mathbf{x}} + 0.02605e^{-0.2566\mathbf{x}} + 0.003042e^{-0.08227\mathbf{x}} \end{aligned}$$

which match 2 *K* moments (including *m*0) of the shifted density.

Based on the Laplace transform of the above approximations it is straightforward to obtain the Laplace transform of *Wq*(*x*) and then by symbolic inverse Laplace transform *Wq*(*x*) itself. We studied the model with *λ* = 1, *c* = 9/4, *q* = 1/10. We considered both the perturbed model with *σ* = 1 and the non-perturbed one. In the following we present figures where dots represent values obtained by numerical Laplace inversion while solid lines are given by the order 4 approximation. Figures 6 and 7 gives the result for the perturbed case and Figures 8 and 9 for the non-perturbed model. Numerical inverse Laplace transform requires about 60 seconds on a modern laptop to obtain a single point of *Wq*(*x*) while the Padé approximation is immediate. Order 4 Padé approximation provides results that are not distinguishable from the exact values. The right panels of Figures 7 and 9 show that increasing the order of the approximation improves accuracy.

**Figure 6.** The scale function *Wq*(*x*) (on the left) and its derivative (on the right) in the case of the perturbed model.

**Figure 7.** On the left the expected value of total discounted dividends, *<sup>V</sup><sup>b</sup>*](*x*), as function of *b* with *x* = 1 in the case of the perturbed model. On the right the absolute error of approximations of *<sup>V</sup><sup>b</sup>*](*x*) of order 4 (blue), 6 (orange) and 8 (green).

**Figure 8.** The scale function *Wq*(*x*) (on the left) and its derivative (on the right) in the case of the unperturbed model.

**Figure 9.** On the left the expected value of total discounted dividends, *<sup>V</sup><sup>b</sup>*](*x*), as function of *b* with *x* = 1 in the case of the unperturbed model. On the right the absolute error of approximations of *<sup>V</sup><sup>b</sup>*](*x*) of order 4 (blue), 6 (orange) and 8 (green).

### **8. The Wiener–Hopf Factorization and Risk Theory with Rational Roots**

Lévy processes are naturally parameterized via the additive decomposition of the symbol Equation (14). The roots and poles of the symbol play an important role and this has lead eventually to the discovery of the well-known Wiener–Hopf factorization Kyprianou (2014)

$$
\eta \hat{\mathcal{W}}\_q(\mathbf{s}) = \frac{q}{\kappa(\mathbf{s}) - q} = -\phi\_q^+(\mathbf{s})\phi\_q^-(\mathbf{s}), q \ge 0. \tag{44}
$$

With mixed exponential claims of rates *βi* for example, this holds with *φ*+*q* (*s*) = 1 <sup>1</sup>−*<sup>s</sup>*/Φ*q*and

$$\phi^{-}(s) = \frac{\prod\_{i=1}^{n} (1 + s/\beta\_i)}{\prod\_{i=1}^{n} (1 + s/\gamma\_i)},\tag{45}$$

where <sup>Φ</sup>*q* and −*γi* are, respectively, the nonnegative and negative roots of the Cramér–Lundberg equation.

One may bypass the numerical difficulties of finding roots by parameterizing a risk process by the poles and roots of its Wiener–Hopf factor *φ*−(*s*). To further simplify matters, one may choose them as integers; then, the parameters of the process will be rational.

The simplest illustration is provided in the tables of (Dufresne and Gerber 1991a, p. 24): The input parameters therein are the interspersed non-zero integer roots *γi* (denoted there by *ri*) and the poles *βi* of the symbol *κ*(*s*); essentially, this is equivalent to specifying Equation (45). The output parameter is the limit of Equation (45) when *s* → <sup>∞</sup>, yielding

$$\prod \frac{\gamma\_i}{\beta\_i} = 1 - \rho\_{\prime}$$

where *ρ* := Ψ(0) = *λm*1/*c* (Dufresne and Gerber (1991a) gives instead the safety load *θ* := *ρ*<sup>−</sup><sup>1</sup> − 1). The other output parameters are partial fraction decomposition coefficients, given by the well-known partial fractions decomposition formulas, see Equation (12) for an example.<sup>9</sup> Once these coefficients are known, it is easy to show that

$$\begin{cases} m\_1^{-1} = f\_\ell(0) = \lim\_{s \to \infty} s \, f\_\ell^\*(s) = \rho^{-1} (\sum \mu\_i - \sum \gamma\_i) \Leftrightarrow \\ \lambda := \frac{\lambda}{\varepsilon} = \frac{\rho}{m\_1} = \sum \mu\_i - \sum \gamma\_i \end{cases}$$

Note that the Levy model of Equation (1) is overdetermined, since by scaling the time, the symbol will be multiplied by an arbitrary factor; thus, one of the parameters, for example, *c*, may be specified at will.

*A Computer Program "Rat" That Outputs a Spectrally Negative Lévy Process with Given Roots and Poles of Its Symbol*

The first program inputs consists of


The algorithm must distinguish between four cases, depending on whether or not *q* = 0 and *σ* = 0, and this is achieved by different inputs. The case *q* > 0 is specified via the presence of a third input parameter <sup>Φ</sup>*q* (which is 0 if *q* = 0), and the cases *σ* > 0/*σ* = 0 are distinguished by having *N* = *n* + 1 and *N* = *n*, respectively.

We will now introduce an artificial normalization constant *qc*,*σ*, defined by

$$q\_{\mathbb{C},\mathcal{I}} = \begin{cases} \frac{q}{\Phi\_{\mathbb{P}}(1-\rho\_{\mathbb{P}})} := \frac{q}{\Phi\_{\mathbb{P}}} \prod \frac{\mathcal{G}\_{i}}{\gamma\_{i}}, & q > 0\\ \frac{1}{\Phi\_{\mathbb{P}}'(0)(1-\rho\_{\mathbb{P}})} = \frac{1}{\Phi\_{\mathbb{P}}'(0)} \prod \frac{\mathcal{G}\_{i}}{\gamma\_{i}}, & q = 0 \end{cases} \tag{46}$$

.

(the factor 1 − *ρq* being incorporated just for convenience, see Equation (48)). Note that this definition allows one to deal with the case *q* = 0 = <sup>Φ</sup>*q* as a limiting case of *q* > 0. *qc*,*<sup>σ</sup>* may be thought of at

<sup>9</sup> Nowadays, the simplest way to obtain them is by solving linear systems with Mathematica, Maple, Sage, etc.

first as the leading coefficient of *κ*(*s*) (which is *c* when *σ* = 0 and *σ*22 when *σ* > 0); for a more precise statement, see Equation (48).

When *q* > 0, we will write the symbol as

$$\kappa(\mathbf{s}) - q \quad = \quad q\_{\varepsilon,\sigma} \left( \frac{\vartheta^2}{2} \mathbf{s}^2 + \mathbb{E}\mathbf{s} - \Phi\_q(1-\rho\_q) + \sum\_{j=1}^n \bar{A}\_j \left( \frac{\beta\_j}{\beta\_j + \mathbf{s}} - 1 \right) \right), \tag{47}$$

where *c*˜ = *cqc*,*σ* , *σ*˜ 2 = *σ*<sup>2</sup> *qc*,*<sup>σ</sup>* , *A*˜*j* = *Aj qc*,*<sup>σ</sup>* .

With this parameterization, results when *q* → 0 will have a limit, provided that *qc*,*<sup>σ</sup>* is kept constant. If we assume (w.l.o.g.) that the polynomial part inside the parentheses in Equation (47) is monic (i.e., *c*˜ = 1, if *σ* = 0, and *σ*˜ 2/2 = 1, else), then the parameter *qc*,*<sup>σ</sup>* equals either the original parameter *c* or *σ*22.

**Remark 8.** *In terms of qc*,*σ, we may rewrite the Wiener–Hopf factorization as*

$$\begin{split} \hat{\mathcal{W}}\_{q}(\mathbf{s}) &= \frac{1}{\kappa(\mathbf{s}) - q} = \frac{\Phi\_{q}}{q(\mathbf{s} - \Phi\_{q})} \Phi\_{q}^{-}(\mathbf{s}) := \frac{\Phi\_{q}}{q} (1 - \rho\_{q}) \frac{\tilde{\phi}\_{q}^{-}(\mathbf{s})}{\mathbf{s} - \Phi\_{q}} = \frac{1}{q\_{\mathbf{c}, \sigma}} \frac{\tilde{\phi}\_{q}^{-}(\mathbf{s})}{\mathbf{s} - \Phi\_{q}} \\ &\Leftrightarrow \kappa(\mathbf{s}) - q = q\_{\mathbf{c}, \sigma} \frac{\mathbf{s} - \Phi\_{q}}{\tilde{\Phi}\_{q}^{-}(\mathbf{s})} \end{split} \tag{48}$$

*where the last factor in Equation* (48) *may be represented as a quotient of monic polynomials. This equation motivated the introduction of qc*,*σ.*

**Remark 9.** *The various components of the cumulant generating function are readily obtained: σ*, *c and q by extracting the polynomial part; the remaining jump part j*(*s*) *then yields λ* = − lim*s*→∞ *j*(*s*)*.*

**Remark 10.** *Note that if the condition*

$$
\gamma\_{n+1} > \beta\_n > \gamma\_n > \dots > \gamma\_2 > \beta\_1 > \gamma\_1 \ge 0 \tag{49}
$$

*is satisfied, then the nonnegativity of the density of the jumps is ensured—see, for example, Kuznetsov et al. (2012) (this condition is also necessary with hyperexponential jumps).*

*If this condition is not satisfied, then the user must deal somehow with the nonnegativity by restricting the position of the Cramér–Lundberg roots (not an easy task).*

The output of our program consists of several first passage characteristics:


They are most easily described as combinations of exponentials.

### **9. Further Background on the Laguerre–Tricomi–Weeks Method**

We now proceed to review an extension of the Laguerre expansion, the "Laguerre–Tricomi–Weeks" method of inverting the Laplace transform, initially proposed by Tricomi in 1935 and McCully (see McCully (1960)), which is considered as one of the most efficient methods of inverting the Laplace transform Weideman (1999). We mention also that an exact explicit closed-form Laguerre series expansion formula was proposed recently in Zhang and Cui (2019).

Consider a given function

$$f(\mathbf{x}) = \mathfrak{e}^{\beta \mathbf{x}} \mathfrak{g}(\mathbf{a} \mathbf{x}), \quad \mathfrak{g} \in L^2 \tag{50}$$

(note that the transformation Equation (50) corresponds to a linear transformation for the Laplace transform *f* (*s*) = *g*(*<sup>s</sup>*−*βα* ) and its singularities). First, one attempts to construct a Laguerre expansion

$$f(\mathbf{x}) = e^{\hat{\beta}\mathbf{x}} \sum\_{n=0}^{\infty} B\_n e^{-\mathbf{a}\mathbf{x}/2} L\_n(\mathbf{a}\mathbf{x}) \Leftrightarrow f(\mathbf{x}/a) e^{-\hat{\beta}\mathbf{x}/\mathbf{a}} = \sum\_{n=0}^{\infty} B\_n e^{-\mathbf{x}/2} L\_n(\mathbf{x})$$

$$\Leftrightarrow \hat{f}(\mathbf{a}\mathbf{s} - \mathbf{a}/2 + \beta) = \mathbf{s}^{-1} \sum\_{n=0}^{\infty} B\_n (1 - \mathbf{a}/\mathbf{s})^n \tag{51}$$

with two judiciously chosen parameters *α* > 0 and *β* > *sf* , where *sf* is the maximum real part of the singular points of *f* <sup>∗</sup>(*s*).

After the change of variables 1 − *α*/*s* = *z*, the coefficients *Bn* may be found from the Taylor expansion

$$\frac{\alpha}{1-z}f^\*\left(\frac{\alpha}{1-z} - \alpha/2 + \beta\right) = \sum\_{n=0}^{\infty} B\_n z^n, |z| < R\_{\Phi} \tag{52}$$

where *R*Φ is the radius of convergence of <sup>Φ</sup>(*z*) = *α*1−*<sup>z</sup> f* ∗ *α*1−*<sup>z</sup>* − *α*/2 + *β*, or, equivalently,

$$a f^\* \left( \frac{n}{1-z} - n/2 + \beta \right) = \sum\_{n=0}^{\infty} b\_n z^n, \quad \text{with } b\_0 = B\_0, b\_n = B\_n - B\_{n-1}, \ n \ge 1.$$

Judicious choices of *β* (related to maximizing *R*Φ) and *α* (related to finding a minimal circle including the singular points of *f* (*s*)) may turn this into one of the most effective inversion methods Weideman (1999), Abate et al. (1998).

Some particular choices are *β* = *α*/2, *β* = 0 and *β* = −*<sup>α</sup>*/2, in which case the Taylor expansion for the coefficients *Bn* becomes, respectively,

$$\frac{n}{1-z}f^\*\left(\frac{n}{1-z}\right) \quad = \sum\_{n=0}^{\infty} B\_n z^n \Leftrightarrow f(\mathbf{x}/n) = \sum\_{n=0}^{\infty} B\_n L\_n(\mathbf{x}) \tag{53}$$

$$\frac{n}{1-z}f^\*\left(\frac{n}{2}\cdot\frac{1+z}{1-z}\right)^{-1}=\sum\_{n=0}^{\infty}B\_nz^n\Leftrightarrow f(\mathbf{x}/n)=\sum\_{n=0}^{\infty}B\_nz^{-\mathbf{x}/2}L\_n(\mathbf{x})\tag{54}$$

$$\frac{n}{1-z}f^\*\left(n\frac{z}{1-z}\right) \quad = \quad \sum\_{n=0}^{\infty}B\_n z^n \Leftrightarrow f(\mathbf{x}/\mathbf{a}) = \sum\_{n=0}^{\infty}B\_n e^{-\mathbf{x}}L\_n(\mathbf{x}).\tag{55}$$

As noted by Keilson and Nunn Keilson and Nunn (1979) (see also Keilson et al. (1980, 1981)), and Abate et al. (1998), respectively, in the second and third cases, the expansion in the of powers Equation (32) implies simple relations between the Laguerre transform and another interesting technique, the Erlang transform, which is a useful tool for the solving integral equations of the convolution type. We refer to these papers for a further discussion of the relative advantages of the Erlang and Laguerre transforms, and the spaces of the functions within which they converge.

**Author Contributions:** F.A., A.H. and S.P. contributed equally to Conceptualization, Methodology, Formal Analysis, Investigation, Writing-Original Draft Preparation, Writing-Review, Project Administration, Supervision and Editing, and U.S. contributed to Software, Validation and Data Curation.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.
