**5. Two Numerical Methods for Computing** *Wq*(*x*)**: The Tijms–Padé and Laguerre–Tricomi–Weeks Approximations**

Computing *Wq*(*x*) reduces in principle to the case *q* = 0, if <sup>Φ</sup>*q* may be accurately computed, via the well-known relation

$$\mathcal{W}\_{\emptyset}(\mathbf{x}) = \boldsymbol{\varepsilon}^{\mathbf{x}} \boldsymbol{\Phi}\_{\emptyset} \, \, \mathcal{W}\_{0}^{(\Phi\_{\emptyset})}(\mathbf{x}), \tag{25}$$

where *W*(<sup>Φ</sup>*q*) 0 (*x*) denotes the 0-scale function with respect to the "Esscher transformed" measure *P*(<sup>Φ</sup>*q*) (in general, the transform *P*(*r*) of the measure *P* of a Lévy process with Laplace exponent *κ*(*s*) is the measure of the Lévy process with Laplace exponent *κ*(*s* + *r*) − *<sup>κ</sup>*(*r*), with *r* in the domain of *<sup>κ</sup>*(·) (Albrecher and Asmussen 2010, Proposition 4.2), (Kyprianou 2014, 3.3 pg.83)).<sup>4</sup>

The Esscher transform removes the exponential growth, or, equivalently, the unique positive pole of *W <sup>q</sup>*(*s*). When working with the Esscher transform, let

$$G(\mathfrak{x}) := \mathcal{W}\_0^{(\Phi\_q)}(\infty) - \mathcal{W}\_0^{(\Phi\_q)}(\mathfrak{x})$$

denote the function after the behavior at <sup>Φ</sup>*q* and ∞ has been exploited<sup>5</sup> We construct then a mixed exponential/rational approximation

$$RG(\mathbf{x}) := \sum\_{i} A\_{i} \mathfrak{e}^{-\gamma\_{i} \mathbf{x}} \Leftrightarrow \widehat{RG}(\mathbf{s}) := \frac{a(\mathbf{s})}{\prod(\mathbf{s} + \gamma\_{i})} \sim \widehat{G}(\mathbf{s}).\tag{26}$$

<sup>4</sup> Equation (25) is easy to check by taking the Laplace transform, but quite important, numerically, since *W*(<sup>Φ</sup>*<sup>q</sup>* ) 0 (*x*) is a monotone bounded function (with values in the interval (lim*s*→∞ *s κ*(*s*) , 1 *<sup>κ</sup>*(<sup>Φ</sup>*q*) )).

<sup>5</sup> Since *W*(<sup>Φ</sup>*q* )(*x*) still converges to a constant when *x* → <sup>∞</sup>, this non-zero limit at ∞ must be removed first.

For general rational approximations, we will use the name of "Esscher–Tijms type algorithms", and for the Padé case we will use the name of "Esscher–Tijms–Padé approximation ", in reference to the fact that the special positive pole, supposed to be known somehow exactly, has been removed.

**Remark 4.** *Clearly, a very precise approximation of the unique positive pole* <sup>Φ</sup>*q and the "moments" <sup>κ</sup>*(*n*)(<sup>Φ</sup>*q*) *are an essential for determining Wq*, *q* > 0*, which was missing in the ruin probability problem.* <sup>Φ</sup>*q is quite easy to compute under parametric models where the Laplace transform of the jumps is explicit, as the unique positive root of a convex function. It may also be estimated empirically by the plug-in method, but this seems not yet studied in this context. We took therefore the shortcut of assuming that* <sup>Φ</sup>*q is available with infinite precision. Conceptually, this amount to studying the Esscher transform.*

The Padé approximation of the Esscher transform is the first method we have investigated. An interesting alternative to the Esscher transform is to consider the simplified Laplace transform:

$$
\hat{\mathcal{W}}(s) := (s - \Phi\_q)\hat{\mathcal{W}}\_q(s) = \frac{\Phi\_q}{q}\phi^-(s), \tag{27}
$$

whose positive pole has been removed and thus coincides with the negative Wiener–Hopf factor *φ*−(*s*) (see Section 8, (44)), up to the constant <sup>Φ</sup>*q q*.

Subsequently, we may apply a Padé (or more sophisticated rational approximation Nakatsukasa et al. (2018)), with the goal of approximating the dominant poles of the presumably imprecisely known transform.

Note that a Padé of *W*˜ (*s*) will then yield immediately a rational approximation of *<sup>W</sup>q*(*s*) and a mixed exponential approximation of *Wq*(*x*). The Padé approximation requires solving

$$
\tilde{\mathcal{W}}\_{\emptyset}(\mathbf{s}) := (\mathbf{s} - \boldsymbol{\Phi}\_{\emptyset}) \hat{\mathcal{W}}\_{\emptyset}(\mathbf{s}) = \frac{\mathbf{s} - \boldsymbol{\Phi}\_{\emptyset}}{(\mathbf{c} - \boldsymbol{\nu}\_{1})\mathbf{s} + \boldsymbol{\nu}\_{2,\boldsymbol{\sigma}}\mathbf{s}^{2}/2 - \boldsymbol{\nu}\_{3}\mathbf{s}^{3}/6 + \dots - q} \sim \frac{P\_{\boldsymbol{n}}(\mathbf{s})}{Q\_{\boldsymbol{n}}(\mathbf{s})} + o(\mathbf{s}^{\boldsymbol{n}}).
$$

For the Cramér–Lundberg case, also fixing

$$w\_0 = \dot{W}\_q(0) = \frac{1}{c}, w\_0' = \dot{W}\_q'(0) = \frac{1}{c}(\frac{q+\lambda}{c} - \Phi\_q) \tag{28}$$

leads to the linear system

$$\mathfrak{c}\left(\sum\_{i=0}^{n-1} a\_i \mathbf{s}^i + \mathbf{s}^n\right) \mathfrak{c}\mathfrak{x}(\mathbf{s}) - q\mathfrak{y} = \left(\mathbf{s} - \Phi\_{\emptyset}\right) \mathfrak{c}\left(\mathbf{s}^n + (a\_{n-1} + \Phi\_{\emptyset} - \frac{\lambda + q}{\mathfrak{c}}) \mathbf{s}^{n-1} + \sum\_{i=0}^{n-2} b\_i \mathbf{s}^i\right) + o(\mathfrak{s}^n).$$

**Remark 5.** *The "Tijms-two-point Padé approximation " with n* = 1 *and w*0 *satisfied by Equation (28) is obtained by solving:*

$$\frac{s-\Phi\_{q}}{c(s-\nu\_{1})s+\cdots-q} \sim \frac{s+a}{c(s+\gamma)} \Leftrightarrow (s-\Phi\_{q})c(s+\gamma) \sim (ps-q)(s+a\_{0}) \implies \gamma = a\frac{q}{c\Phi\_{q}},\\ a = \frac{q-c\Phi\_{q}}{p-q/\Phi\_{q}}.\tag{29}$$

*This is exact for the Cramér–Lundberg process with exponential jumps of rate μ, since then a* = *μ and* − *γ reduces to the second root of the Cramér–Lundberg equation.*

*Applying Equation* (21) *to Equation* (29) *with A*1 = *a* <sup>Φ</sup>*q*+*<sup>γ</sup>* , *A*2 = 1 − *A*1 *yields then an approximation for the optimal barrier.*

*Starting with order two, we may also satisfy w* 0*. We must now solve*

$$\begin{split} &\frac{s-\Phi\_{q}}{(c-\nu\_{1})s+\cdots-q} \sim \frac{s^{2}+a\_{1}s+a\_{0}}{c(s^{2}+(a\_{0}+\Phi\_{q}-\frac{\lambda+q}{c})s+b\_{0})} \\ &\Leftrightarrow (s-\Phi\_{q})c(s^{2}+(a\_{0}+\Phi\_{q}-\frac{\lambda+q}{c})s+b\_{0}) \sim (ps+\nu\_{2}s^{2}/2-q)(s^{2}+a\_{1}s+a\_{0}) \\ &\Longrightarrow b\_{0}=a\_{0}\frac{q}{c\Phi\_{q}},a\_{1}q=a\_{0}(p-q/\Phi\_{q}+c\Phi\_{q})+c\Phi\_{q}(\Phi\_{q}-\frac{\lambda+q}{c}),a\_{0}=\dots \end{split}$$

*Tedious computations or the help of Mathematica reveal that this is exact for the Cramér–Lundberg process with mixed exponential jumps of order* 2*.*

In the next subsection we propose a method similar in spirit to the above: Obtaining first the simplest exponential approximation, and refining this subsequently (analogously to Step 3) above by means of applying the Laguerre–Tricomi–Weeks expansion.

### *5.1. Laguerre–Tricomi–Weeks Laplace Transform Inversion with Prescribed Exponent*

In this section we obtain first an exponential approximation of the transformed scale function of the form *W*(<sup>Φ</sup>*q*) (*x*) ∼ *W*(<sup>Φ</sup>*q*) (∞) − *Cα* 2 *e*<sup>−</sup> *α* 2 *x*, by a first order Padé approximation<sup>6</sup> of its Laplace transform

$$\widehat{\mathcal{W}^{(\Phi\_q)}}(s) = \frac{1}{\kappa(s + \Phi\_q) - q} = \frac{1}{\kappa(s + \Phi\_q) - \kappa(\Phi\_q)} := \frac{1}{\kappa^{(\Phi\_q)}(s)}.\tag{30}$$

The exponent *α*/2 thus obtained is the "second step output" form the previous section.

Next, apply the Tricomi–Weeks method for Laplace transform inversion (see, for example, Abate et al. (1998)), i.e., search for a Laguerre series expansion

$$G(\mathbf{x}) := \mathcal{W}^{(\Phi\_{\emptyset})}(\infty) - \mathcal{W}^{(\Phi\_{\emptyset})}(\mathbf{x}) \sim \mathbb{C} \sum\_{n=0}^{\infty} B\_n \, \mathrm{e}^{-a\mathbf{x}/2} L(n, a\mathbf{x}) \tag{31}$$

where, following tradition, *C* normalizes the sum following it to be a pdf7, and where

$$L\_n(\mathbf{x}) = \frac{\varepsilon^\mathbf{x}}{n!} \frac{d^n}{d\mathbf{x}^n} (e^{-\mathbf{x}}\mathbf{x}^n) = \sum\_{k=0}^n \binom{n}{k} \frac{(-\mathbf{x})^k}{k!}, \mathbf{x} \ge 0, n = 0, 1, \dots, \tag{32}$$

denote the Laguerre polynomials, which are orthogonal with respect to the weight e<sup>−</sup>*<sup>x</sup>*/2. The Laplace transform of the Laguerre polynomial is

$$\hat{L}(n,s) = \frac{(s-1)^n}{s^{n+1}}, n = 0, 1, \dots,\tag{33}$$

The Laguerre–Tricomi–Weeks method is based on the fact that Equation (31) is equivalent to

$$
\hat{G}(s) := \frac{W^{(\Phi\_q)}(\infty)}{s} - \hat{W}^{(\Phi\_q)}(s) \sim \mathbb{C} \sum\_{n=0}^{\infty} B\_n \frac{(s - n/2)^n}{(s + n/2)^{n+1}} \Leftrightarrow
$$

$$
\hat{H}(s) := (s + n/2)G(s) \sim \mathbb{C} \sum\_{n=0}^{\infty} B\_n \frac{(s - n/2)^n}{(s + n/2)^n}.\tag{34}
$$

<sup>6</sup> Higher-order rational approximations may be considered as well.

<sup>7</sup> C can also be included in the coefficients *Bn*, but introducing it does render the computation of Equation (40) more convenient.

Now the last RHS may be obtained using the "collocation transformation" *z* = *<sup>s</sup>*−*<sup>α</sup>*/2 *s*+*α*/2 ⇔ *s* = *α*2 1+*z* 1−*<sup>z</sup>* and a power series expansion of the LHS, yielding

$$H\left(\frac{a}{2}\frac{1+z}{1-z}\right) = \frac{a}{1-z}G\left(\frac{a}{2}\frac{1+z}{1-z}\right) \sim \mathbb{C}\sum\_{n=0}^{\infty}B\_n\ z^n.\tag{35}$$

In conclusion, extracting the Taylor coefficients *CBn* of *H*( *α*2 1+*z* 1−*<sup>z</sup>* ) yields approximations to *W*(<sup>Φ</sup>*q*) by substituting it into Equation (31). After multiplication by *<sup>e</sup>*Φ*qx*, we obtain approximations of *Wq*(*x*).

**Determining** *α*. Previous work on choosing *α* involved the radius of convergence of a Laguerre–Tricomi–Weeks expansion which is slightly more general than Equation (35) Weideman (1999)—see also Section 9. We now introduce a different method.

Recall that our Laplace transform *W* (<sup>Φ</sup>*q*)(*s*) has been shifted by <sup>Φ</sup>*q*, so that 0 is the largest pole, and that the growth at ∞ has been removed by subtracting


$$\mathcal{W}^{(\Phi\_q)}(\infty) = \lim\_{s \to 0} \widehat{\mathcal{W}^{(\Phi\_q)}}(s) = \frac{s}{\mathfrak{x}(s + \Phi\_q) - q} = \frac{1}{\mathfrak{x}'(\Phi\_q)} = \Phi'\_q. \tag{36}$$

Now, let

$$\mathcal{C} := \lim\_{s \to 0} \hat{G}(s) = \lim\_{s \to 0} \frac{\kappa(s + \Phi\_q) - q - s\kappa'(\Phi\_q)}{s\kappa'(\Phi\_q)(\kappa(s + \Phi\_q) - q)} = \lim\_{s \to 0} \frac{s^2 \kappa''(\Phi\_q)/2}{s\kappa'(\Phi\_q)(s\kappa'(\Phi\_q))} = \frac{\kappa''(\Phi\_q)}{2(\kappa'(\Phi\_q))^2} \tag{37}$$

denote the mass of our measure.

A reasonable *α* may be obtained by approximating *<sup>G</sup>*(*x*) via a first order exponential approximation of *C*

$$G(x) \sim \mathbb{C}\frac{\mathfrak{a}}{2} \mathfrak{e}^{-\frac{\mathfrak{a}}{2}x} \Leftrightarrow \hat{G}(s) \sim \mathbb{C}\frac{\mathfrak{a}/2}{s + \mathfrak{a}/2} \tag{38}$$

(or, equivalently, by fitting the first two moments of *G* (*s*)), so that the Laguerre exponent in Equation (31) is fitted at order *n* = 0.

Fitting moments may be achieved by Padé's method. At order 1, this yields

$$\begin{split} \hat{G}(s) &= \frac{\kappa(s + \Phi\_{q}) - q - \kappa\kappa'(\Phi\_{q})}{\kappa\kappa'(\Phi\_{q})(\kappa(s + \Phi\_{q}) - q)} \sim \frac{\mathsf{C}\kappa/2}{\kappa/2 + s} \Leftrightarrow \\ \kappa/2 &= \lim\_{s \to 0} \frac{s(\kappa(s + \Phi\_{q}) - q - \kappa\kappa'(\Phi\_{q}))}{(\kappa(s + \Phi\_{q}) - q)(\mathsf{C}\kappa\kappa'(\Phi\_{q}) - 1) + \kappa\kappa'(\Phi\_{q})} = \\ \lim\_{s \to 0} \frac{s(\kappa(s + \Phi\_{q}) - q - \kappa\kappa'(\Phi\_{q}))}{(\kappa\kappa'(\Phi\_{q}) + s^{2}\kappa''(\Phi\_{q})/2 + s^{3}\kappa'''(\Phi\_{q})/6 + \dots)(\mathsf{C}\kappa\kappa'(\Phi\_{q}) - 1) + \kappa\kappa'(\Phi\_{q})} = \\ \lim\_{s \to 0} \frac{s(\kappa(s + \Phi\_{q}) - q - \kappa\kappa'(\Phi\_{q}))}{s^{3}\frac{\kappa''(\Phi\_{q})}{2}\mathsf{C}\kappa'(\Phi\_{q}) - \kappa\kappa''(\Phi\_{q})/6} = \frac{\frac{\kappa''(\Phi\_{q})}{2}}{\mathsf{C}\kappa'(\Phi\_{q})\frac{\kappa''(\Phi\_{q})}{2} - \frac{\kappa'''(\Phi\_{q})}{6}} = \frac{\frac{\kappa''(\Phi\_{q})}{2}}{\frac{\kappa''(\Phi\_{q})}{2\kappa'(\Phi\_{q})}\frac{\kappa'''(\Phi\_{q})}{2} - \frac{\kappa'''(\Phi\_{q})}{6}} \end{split} \tag{39}$$

Finally, a bit of algebra yields

$$\alpha/2 = \frac{\frac{\mathbf{x}^{\prime\prime}(\Phi\_{q})}{2}}{\frac{(\mathbf{x}^{\prime\prime}(\Phi\_{q})/2)^{2}}{\mathbf{x}^{\prime}(\Phi\_{q})} - \frac{\mathbf{x}^{\prime\prime}(\Phi\_{q})}{6}} = \frac{6\mathbf{x}^{\prime}(\Phi\_{q})\mathbf{x}^{\prime\prime}(\Phi\_{q})}{3(\mathbf{x}^{\prime\prime}(\Phi\_{q}))^{2} - 2\mathbf{x}^{\prime}(\Phi\_{q})\mathbf{x}^{\prime\prime\prime}(\Phi\_{q})}.\tag{40}$$

**Remark 6.** *It may also be checked that*

$$G(0) = \lim\_{s \to \infty} \frac{1 - s \frac{\mathbf{x}'(\Phi\_q)}{\mathbf{x}(s + \Phi\_q) - q}}{\mathbf{x}'(\Phi\_q)} = \begin{cases} \frac{1}{\mathbf{x}'(\Phi\_q)} & \text{unbounded variation} \\\frac{1}{\mathbf{x}'(\Phi\_q)} - \frac{1}{\varepsilon} & \text{bounded variation} \end{cases}$$

*and*

$$G'(0) = \lim\_{s \to \infty} s \hat{G}'(s) = \begin{cases} \frac{2}{3\mathbf{k}'(\Phi\_\emptyset)} & \text{unbounded variation} \\ \frac{1}{2\mathbf{k}'(\Phi\_\emptyset)} - \frac{1}{2\mathbf{c}} & \text{bounded variation} \end{cases}$$

*This will be useful for developing higher order two-point Padé approximations of <sup>G</sup>*(*s*)*.*

**Remark 7.** *It is possible that the choice in Equation* (40) *for the Laguerre–Tricomi–Weeks exponential parameter (often left to the user) has previously been proposed, but we have been unable to find this result in the literature. Note the natural generalization to the case of Laplace transforms with a finite number of positive poles (instead of a single one, as in our case).*

We show in the next sections that the approximations obtained by the Tijms–Padé and Laguerre–Tricomi–Weeks methods are accurate enough so that the resulting inverse *Wq*, *W q* retain the features of the intricate Azcue and Loeffen examples Azcue and Muler (2005); Loeffen (2008).

### *5.2. Combining the Tijms–Padé and Laguerre–Tricomi–Weeks Approximations*

We end this section by proposing a more general strategy for inverting *<sup>W</sup>q*(*s*) and other Laplace transforms of interest in probability, given the uncertainty inherent in real world data. This consists in refining the rational approximation of Equation (26) by constructing the quotient of the original Laplace transform and its approximation

$$\hat{G}(s) := \frac{\hat{G}(s)}{\widehat{R}G(s)} = \hat{G}(s)\frac{\prod(s+\gamma\_i)}{a(s)}\tag{41}$$

.

and by applying to it an inversion method with good convergence properties like Laguerre–Tricomi–Weeks. The final approximation of *Wq*(*x*) will be then a convolution of the inverse Laplace transforms *RG*(*x*) and *<sup>G</sup>*˜(*x*). This may be viewed as a to a two stage filter, aimed at removing uncertainties of different orders of magnitude.
