# Optimizing dividends and capital injections limited by bankruptcy, and
practical approximations for the Cramér-Lundberg process

F. Avram 111Corresponding author. E-mail<EMAIL_ADDRESS>Laboratoire
de Mathématiques Appliquées, Université de Pau, France 64000  D. Goreac  LAMA,
Univ Gustave Eiffel, UPEM, Univ Paris Est Creteil, CNRS, F-77447 Marne-la-
Valleée, France  School of Mathematics and Statistics, Shandong University,
Weihai, Weihai 264209, PR China  R. Adenane  Département des Mathématiques,
Université Ibn-Tofail, Kenitra, Maroc 14000  U. Solon  Laboratoire de
Mathématiques Appliquées, Université de Pau, France 64000  Institute of
Mathematics, University of the Philippines Diliman, Philippines

###### Abstract

The recent papers Gajek-Kucinsky(2017) and Avram-Goreac-Li-Wu(2020)
investigated the control problem of optimizing dividends when limiting capital
injections stopped upon bankruptcy. The first paper works under the spectrally
negative Lévy model; the second works under the Cramér-Lundberg model with
exponential jumps, where the results are considerably more explicit. The
current paper has three purposes. First, it illustrates the fact that quite
reasonable approximations of the general problem may be obtained using the
particular exponential case studied in Avram-Goreac-Li-Wu(2020). Secondly, it
extends the results to the case when a final penalty $P$ is taken into
consideration as well besides a proportional cost $k>1$ for capital
injections. This requires amending the “scale and Gerber-Shiu functions"
already introduced in Gajek-Kucinsky(2017). Thirdly, in the exponential case,
the results will be made even more explicit by employing the Lambert-W
function. This tool has particular importance in computational aspects and can
be employed in theoretical aspects such as asymptotics.

Keywords: dividend problem, capital injections, penalty at default, scale
functions, Lambert-W function, De Vylder-type approximations, rational Laplace
transform

###### Contents

  1. 1 Introduction
  2. 2 The cost function of $(-a,0,b)$ policies, for the spectrally negative Lévy case
    1. 2.1 The HJB System
  3. 3 Explicit determination of $a^{*},b^{*}$ when $F(x):=1-e^{-\mu x},P>-\frac{c}{q}$
    1. 3.1 The simplified cost function and optimality equations
    2. 3.2 Existence of the roots of the equations $\eta(b)=0,\delta_{k,P}=0$
  4. 4 Which exponential approximation?
    1. 4.1 Three de Vylder-type exponential approximations for the ruin probability
    2. 4.2 Three two point Padé approximations of the Laplace transform $\widehat{W}_{q}$ of scale function
  5. 5 Examples of computations involving scale function and dividend value approximations
    1. 5.1 A Cramér-Lundberg process with hyperexponential claims of order 2
    2. 5.2 A Cramér-Lundberg process with hyperexponential claims of order 3
    3. 5.3 A Cramér-Lundberg process with oscillating density and scale function
  6. 6 The maximal error of exponential approximations $J_{0}$ along one parameter families of Cramér-Lundberg processes
  7. 7 The profit function when the claims are distributed according to a matrix exponential jumps density

##  1 Introduction

This paper concerns the approximate optimization of a new type of boundary
mechanism, which emerged recently in the actuarial literature [APY18, GK17,
AGLW20], in the context of the optimal control of dividends and capital
injections.

The model. Consider the spectrally negative Lévy risk model:

$X_{t}=x+ct-\bar{\xi}_{t},c\geq 0,$

where $\bar{\xi}_{t}$ is a spectrally positive Lévy process, with Lévy measure
$\nu(x)dx$. The classic example is that of the perturbed Cramèr-Lundberg risk
model with

$\bar{\xi}_{t}=\sum_{i=1}^{N_{t}}\xi_{i}+\sigma B_{t},$

where $B_{t}$ is a Brownian motion, where $N_{t}$ is an independent Poisson
process of intensity $\lambda>0$, and $\left(\xi_{i}\right)_{i\geq 1}$ is a
family of i.i.d.r.v. whose distribution, density and moments are denoted
respectively by $F,f,m_{i},i\in\\{1,...\\}$. Furthermore,

  * •

the process is modified by dividends and capital injection:

$\pi:=\left(D,I\right)\Rightarrow X_{t}^{\pi}:=X_{t}-D_{t}+I_{t},$

where $D,I$ are adapted, non-decreasing and càdlàg processes, with
$D_{0-}=I_{0-}=0$;

  * •

the first time when we do not bail-out to positive reserves
$\sigma^{\pi}_{0}:=\inf\left\\{t>0:\
X_{t-}^{\pi}-\bigtriangleup\bar{\xi}_{t}+\bigtriangleup I_{t}<0\right\\}$ is
called bankruptcy/absolute ruin;

  * •

prior to bankruptcy, dividends are limited by the available reserves:
$\bigtriangleup D_{t}:=D_{t}-D_{t-}\leq
X_{t-}^{\pi}-\bigtriangleup\bar{\xi}_{t}+\bigtriangleup I_{t},\textnormal{
where }\bar{\xi}_{t}:=\sum_{i=1}^{N_{t}}\xi_{i}$. The set of “admissible"
policies satisfying this constraint will be denoted by $\tilde{\Pi}(x)$.

The objective is to maximize the profit:

$J_{x}^{\pi}:=\mathbb{E}_{x}\left[\int_{\left[0,\sigma^{\pi}_{0}\right]}e^{-qt}\left(dD_{s}-kdI_{s}\right)-Pe^{-q\sigma^{\pi}_{0}}\right],k\geq
1,P\geq 0.$

The value function is

$V\left(x\right):=\sup_{\pi\in\tilde{\Pi}(x)}J_{x}^{\pi},\ x\in\mathbb{R}.$

Motivation. The recent papers [GK17, AGLW20] investigated the above control
problem of optimizing dividends and capital injections for processes with
jumps, when bankruptcy is allowed as well. The second paper works under the
Cramér-Lundberg model with exponential jumps, while the first works under the
spectrally negative Lévy model, allowing also for the presence of Brownian
motion and infinite activity jumps. It turns out that the optimal policy
belongs to the class of $(-a,0,b),a>0,b\geq 0$ “bounded buffer policies",
which consist in allowing only capital injections smaller than a given $a$ and
declaring bankruptcy at the first time when the size of the overshoot below 0
exceeds $a,$ and in paying dividends when the reserve reaches an upper barrier
$b$. These will briefly be described as $(-a,0,b)$ policies from now on.
Furthermore, the optimal $(a^{*},b^{*})$ are the roots of one variable
equations with explicit solutions related to the Lambert-W(right) function
(ProductLog in Mathematica).

Below, our goal is to show numerically that exponential approximations provide
quite reasonable results (as the de Vylder approximation provides for the ruin
problem). We will focus in our examples on the case of matrix exponential
jumps(which are known to be dense in the class of general nonnegative jumps,
with even error bounds for completely monotone jumps being available
[VAVZ14]), for two reasons. One is in order to highlight certain exact
equations which are similar to their exponential versions, and which may at
their turn be used to produce even more accurate approximations in the future,
and, secondly, since numerical Laplace inversion for this class may easily
tuned to have arbitrarily small errors.

History of the problem: The case of no capital injections (also characterized
by $k=\infty$ or absorption below 0) is the dividend problem posed by De
Finetti [DF57, Ger69] where dividends are paid above barrier $b^{*}$ and
$a^{*}=0$ is imposed. “The challenge is to find the right compromise between
paying early in view of the discounting or paying late in order not to reach
ruin too early and thus profit from the positive safety loading for a longer
time" [AAM20].

Forced injections and no bankruptcy at $0$ (also characterized by a reflection
at $0$) is studied in Shreve [SLG84] where dividends are paid above barrier
$b^{*}$ and $a^{*}=\infty$ is imposed.

From Lokka, Zervos, [LZ08] we know that in the Brownian motion case, it is
optimal to either always inject, if $k\leq k_{c}$, for some critical cost
$k_{c}$ (i.e. use Shreve), or, stop at $0$ (use De Finetti). We propose to
call this the Lokka, Zervos alternative. The “proof" of this alternative
starts by largely assuming it via a heuristically justified border Ansatz
[LZ08, (5.2)]: $\max\left\\{-V(0),\
V^{\prime}(0)-k\right\\}=0\Longrightarrow\text{ either }V(0)=0\text{ or
}V^{\prime}(0)=k.$

Extensive literature on SLG forced bailouts (no bankruptcy) can be found at
Avram et al., (2007) [APP07], Kulenko and Schmidli, (2008) [KS08], Eisenberg
and Schmidli, (2011) [ES11], Pérez et al., (2018) [PYB18], Lindensjo, Lindskog
(2019) [LL19], Noba et al., 2020 [NPY20].

Articles [GK17, AGLW20] are the only papers which relate declaring bankruptcy
to the size of jumps, with general and exponential jumps, respectively. [GK17]
deals also with the presence of Brownian motion and infinite activity jumps,
by conditioning at the first draw-down time; the optimality proof is quite
involved.

In [AGLW20], it is also shown that neither $V(0)=0$ nor  $V^{\prime}(0)=k$ are
possible: the Lokka-Zervos alternative disappears, but another interesting
alternative holds. Above a certain critical $k_{c}$ the optimal dividends
barrier switches from strictly positive to $0$, and $k_{c}$ is related in (26)
to the Lambert-W function .

The results of [GK17, AGLW20] may be divided in three parts:

  1. 1.

Compute the value of bounded buffer policies. The key result is (5) below, an
explicit determination of the objective $J_{0}=J_{0}^{a,b}$, which allows
optimizing it numerically.

######  REMARK 1.

Computing the value function is considerably simplified by the use of first
passage recipes available for spectrally negative Lévy processes [AKP04,
Kyp14, KKR13, AGVA19], which are built around two ingredients: the $W_{q}$ and
$Z_{q}$ $q$-scale functions, defined respectively for $x\geq 0,q\geq 0$ as:

    1. (a)

the inverse Laplace transform of $\frac{1}{\kappa(s)-q}$, where $\kappa(s)$ is
the Laplace exponent (which characterizes a Lévy process) and

    2. (b)

$Z_{q}(x)=1+q\int_{0}^{x}W_{q}(y)dy$

– see the papers [Sup76, Ber98, AKP04] for the first appearance of these
functions. The name $q$-scale/harmonic functions is justified by the fact that
these functions are harmonic for the process $X$ killed upon entering
$(-\infty,0)$, in the sense that

$\\{e^{-q\min[t,T_{0}]}\;W_{q}(X_{\min[t,T_{0}]}),e^{-q\min[t,T_{0}]}\;Z_{q}(X_{\min[t,T_{0}]})\\},t\geq
0$

are martingales, as shown in [Pis04, Prop. 3] (in the case of $Z_{q}$, there
is also a penalty of $1$ at ruin, generalizing to other penalties produces the
so-called Gerber-Shiu function).

  2. 2.

Equations determining candidates for the optimal $a^{*},b^{*}$ are obtained by
differentiating the objective (which is expressed in terms of the scale
functions $W_{q},Z_{q}$), and the optimal pair $(a^{*},b^{*})$ is identified.
As a result, the critical $k_{c}$ is related in (26) to the Lambert-W
function.

  3. 3.

The optimality of the $(-a^{*},0,b^{*})$ policy is established.

Note that the last step is quite non-trivial and is achieved by different
methods in [GK17] and [AGLW20]. The latter paper starts by formulating a (new)
HJB equation associated to this stochastic control problem – see (8).

######  REMARK 2.

The objective may be optimized numerically using the first step only (the
equation (10) for $J_{0}^{a,b}$).

Exponential approximations may also be used, which are similar in spirit with
the de Vylder-type approximations. Recall that the philosophy of the de Vylder
approximation is to approximate a Cramèr-Lundberg process by a simpler process
with exponential jumps, with cleverly chosen exponential rate $\mu$, and the
parameters $\lambda,c$ may also be modified, if one desires to make the
approximations exact at $x=0$ –see for example [AHPS19] for more details)

The efficiency of the de Vylder approximation for approximating ruin
probabilities is well documented [DV78]. The natural question of whether this
type of techniques may work for other objectives, like for example for
optimizing dividends and/or reinsurance was already discussed in [Høj02, DD05,
BDW07, GSS08, AHPS19]. In this paper, following on previous works [ACFH11,
AP14, ABH18], we draw first the attention to the fact that we have not one,
but three de Vylder-type approximations for $W_{q}(x)$ (as for the ruin
probability). The best approximation in our experiments when the loading
coefficient $\theta$ is large turn out to be the classic de Vylder
approximation. However, for approximating near the origin, the two point Padé
approximation which fixes both the values
$W_{q}(0)=\frac{1}{c},W_{q}^{\prime}(0)=\frac{q+\lambda}{c^{2}}$ works better.
The end result here is simply replacing the inverse rate $\mu^{-1}$ by
$m_{1},$ in the formula for the scale function of the Cramèr-Lundberg process
with exponential jumps. In between $x=0$ and $x\to\infty$, the winner is
sometimes the “Renyi approximation" which replaces the inverse rate by
$\frac{m_{2}}{2m_{1}}$, and modifies $\lambda$ as well (for the de Vylder
approximation, $\mu^{-1}$ is replaced by $\frac{m_{3}}{3m_{2}}$, and both
$\lambda,c$ are modified).

We end this introduction by highlighting in figure 1 the fact that for
exponential jumps, the limited capital injections objective function $J_{0}$
given by (14) for arbitrary $b$ but optimal $a=s(b)$ (via a complicated
formula) improves the value function with respect to de Finetti and Shreve,
Lehoczky and Gaver, for any $b$.

Figure 1:  The value function $J_{0}$ given by (14), for arbitrary $b$ but
optimal $a$. The inequality observed is a consequence of the properties of the
Lambert function. The improvement with respect to de Finetti is considerable,
of $0.382292\%$ (the SLG approach is not competitive in this case). Note also
that the optimal barrier $b=0.109023$ is smaller than the de Finetti and SLG
optima of $0.626672,1.82726$ respectively.

Contents and contributions. Section 2 offers a conjectured profit formula for
$(-a,0,b)$ policies, where we include also a final penalty $P$. The
theoretical result of the section 1 revisits [GK17, AGLW20] by linking the two
formulations together and emphasizing the impact of the bankruptcy penalty $P$
(via the scale function $G$). Its proof is beyond the scope of the present
(already lengthy enough) paper and it can be inferred from either one of
[GK17] and [AGLW20] through a three step argument:

  1. 1.

express the cost by conditioning on the reserve ($J_{x}$) starting from $0\leq
x\leq b$ hitting either 0 or b;

  2. 2.

get a further relationship on costs $J_{b}$ and $J_{0}$ by conditioning on the
first claim;

  3. 3.

finally, mix these conditions together in order to obtain the explicit formula
for $J_{x}$.

We also wish to point out (in section 2.1) the link to an appropriate HJB
variational inequality equation8 and the definition 1 specifying the action
regions and their computation starting from the regimes in the HJB system. To
the best of our knowledge, this is new and the preliminary studies conducted
on more complicated problems (involving reinsurance and reserve-dependent
premium) seem to reinforce the relevance of this tool.

In section 7 we provide an alternative matrix exponential form of the exact
cost, in the case of matrix exponential jumps.

An explicit determination of $a^{*},b^{*}$ and an equity cost dichotomy when
dealing with exponential jumps are given in section 3, taking also advantage
of properties of the Lambert-W function, which were not exploited before. The
two main novelties of the section are:

  * •

emphasizing the computations of the optimal buffer/barrier (from [AGLW20]) in
relation with the scale-like quantities appearing in [GK17];

  * •

making explicit use of the (computation-ready) Lambert-W function to describe
the dependency of optimal $a^{*}b$ (in equation 14) and of the dichotomy-
triggering cost $k_{c}$ in equation 26.

Again, a further novelty is the presence of the bankruptcy cost $P$.

Section 4 reviews, for completeness, the de Vylder approximation-type
approximations. Section 4.1 recalls, for warm-up, some of the oldest
exponential approximations for ruin probabilities. Section 4.2 recalls in
Proposition 5, following [AP14, AHPS19] three approximations of the scale
function $W_{q}(x)$ 444essentially, this is the “dividend function with fixed
barrier”, which had been also extensively studied in previous literature
before the introduction of $W_{q}(x)$ , obtained by approximating its Laplace
transform. These amount finally to replacing our process by one with
exponential jumps and cleverly crafted parameters based on the first three
moments of the claims.

In section 5, we consider particular examples and obtain very good
approximations for two fundamental objects of interest: the growth exponent
$\Phi_{q}$ of the scale function $W_{q}(x)$, and the (last) global minimum of
$W_{q}^{\prime}(x)$, which is fundamental in the de Finetti barrier problem.
Proceeding afterwards to the problem of dividends and limited capital
injections, concepts in section 3 are used to compute a straightforward
exponential approximation based on an exponential approximation of the claim
density, and a new “correct ingredients approximation" which consists of
plugging into the objective function (10) for exponential claims the exact
"non-exponential ingredients" (scale functions and, survival and mean
functions) of the non-exponential densities. Both methods are observed to
yield reasonable values in approximating the objective.

This leads us to our conclusion that from a practical point of view,
exponential approximations are typically sufficient in the problems discussed
in this paper.

##  2 The cost function of $(-a,0,b)$ policies, for the spectrally negative
Lévy case

In this section, we allow $\bar{\xi}_{t}$ to be a spectrally positive Lévy
process, with a Lévy measure admitting a density $\nu(dy)=\nu^{\prime}(y)dy$.
The simplest example is that of the perturbed Cramèr-Lundberg risk model with

$\bar{\xi}_{t}=\sum_{i=1}^{N_{t}}\xi_{i}+\sigma B_{t},$

where $N_{t}$ is a Poisson process of intensity $\lambda>0$,
$\left(\xi_{i}\right)_{i\geq 1}$ is an independent family of i.i.d.r.v. with
density $f(y)$, and $B_{t}$ is an independent Brownian motion.

We revisit here the problem of optimizing the value of "bounded buffer
$(-a,0,b)$ policies", following [GK17, AGLW20] (in order to relate the
results, one needs to replace $\gamma$ in the objective of [GK17] by $1/k$),
while taking into account also the bankruptcy penalty $P$.

An important role in the results will be played by the expected scale after a
jump

$C(x)=\int_{0}^{x}W_{q}(x-y)\overline{\nu}(y)dy=c{W_{q}(x)}-Z_{q}(x)+\frac{\sigma^{2}}{2}W_{q}^{\prime}(x),$
(1)

where $\overline{\nu}(y)=\int_{y}^{\infty}\nu(u)du$ is the tail of the Lévy
measure and $\sigma$ is the Brownian volatility (the identity above follows
easily from the $q$-harmonicity of $Z_{q}$, after an integration by parts of
the convolution term and a division by $q$).

The problem of limited reflection requires introducing a new "scale function
$S_{a}(x)$ and Gerber-Shiu function $G_{a,\sigma}(x)$"– see Remark 4 for
further comments on this terminology:

$\begin{cases}S_{a}(x)=Z_{q}(x)+C_{a}(x),C_{a}(x)=\int_{0}^{x}W_{q}(x-y)\;\overline{\nu}(a+y)\;dy\\\
G_{a,\sigma}(x)=G_{a}(x)+k\frac{\sigma^{2}}{2}W_{q}(x)\end{cases}$ (2)

where

$\displaystyle
G_{a}(x)=\int_{0}^{x}W_{q}(x-y)\left(k\;m_{a}(y)+P\overline{\nu}(a+y)\right)dy:=kM_{a}(x)+PC_{a}(x),$
$\displaystyle m_{a}(y)=\int_{0}^{a}z\nu(y+z)dz.$

######  EXAMPLE 1.

With exponential jumps and possibly $\sigma>0$, using the identities

$\overline{\nu}(y+a)=e^{-\mu a}\overline{\nu}(y),m_{a}(y)=\lambda e^{-\mu
y}m(a),m(a)=\int_{0}^{a}y\;\mu e^{-\mu y}dy=\frac{1-e^{-\mu a}}{\mu}-ae^{-\mu
a},$

we find that the functions (2) are expressible as products of $C(x)$ and the
survival or mean function of the jumps:

$\begin{cases}C_{a}(x)=C(x)e^{-\mu
a}=C(x)\overline{F}(a),\;\overline{F}(a)=1-F(a)\\\ S_{a}(x)=Z_{q}(x)+e^{-\mu
a}C(x)\\\ G_{a}(x)=\left(km(a)+P\overline{F}(a)\right)C(x)\end{cases}$ (3)

([GK17] use $s_{c},r_{c}$, instead of
$M_{a}(x):=\int_{0}^{x}W_{q}(x-y)\;m_{a}(y)\;\;dy,C_{a}(x)$, respectively).
When $P=0=\sigma$, these reduce to quantities in [AGLW20].

The formulas above will be used below as a heuristic approximation in non-
exponential cases.

######  REMARK 3.

Note that

$C_{a}(0)=0,G_{a}(0)=0,S_{a}(0)=1,C(0)=0,C^{\prime}(0)=\begin{cases}\frac{\lambda}{c}&\sigma=0\\\
0&\sigma>0\end{cases},$ (4)

and that $C(x),G_{a}(x),S_{a}(x)$ are increasing functions in $x$.

We state now a generalization of [GK17, Thm. 4] for the value function
$J_{0}^{a,b}$ of $(-a,0,b)$ policies, in terms of $S_{a}(x),G_{a}(x)$. In the
Cramèr-Lundberg case illustrated below, the proof is straightforward,
following [AGLW20]. In the other case, one needs to adapt the proof of [GK17].

######  THEOREM 1.

Cost function for $(a,b)$ policies For a spectrally negative Lévy processes,
let

$J_{x}=J^{a,b}(x):=\mathbb{E}_{x}\left[\int_{0}^{T_{-a}}e^{-qt}\left(\mathrm{d}D_{t}-k\;\mathrm{d}I_{t}\right)-Pe^{-qT_{-a}}\right]$

denote the expected discounted dividends minus capital injections associated
to policies consisting in paying capital injections with proportional cost
$k\geq 1$, provided that the severity of ruin is smaller than $a>0$, and
paying dividends as soon as the process reaches some upper level $b$. Put

$\displaystyle G_{a,\sigma}(x)=G_{a}(x)+k\frac{\sigma^{2}}{2}W_{q}(x).$

Then, it holds that

$J_{x}=\begin{cases}G_{a,\sigma}(x)+J_{0}^{a,b}S_{a}(x)=G_{a,\sigma}(x)+\frac{1-G_{a,\sigma}^{\prime}(b)}{S_{a}^{\prime}(b)}S_{a}(x),&x\in[0,b]\\\
kx+J_{0}^{a,b}&x\in[-a,0]\\\ 0&x\leq-a\end{cases}.$ (5)

######  REMARK 4.

The first equality in (5) will be easily obtained by applying the strong
Markov property at the stopping time $T=\min[T_{0,-},T_{b,+}]$, but it still
contains the unknown $J_{0}$.

This relation suggests a definition of the scale $S_{a}$ and the Gerber-Shiu
function $G_{a,\sigma}$, as the coefficient of $J_{0}$ and the part
independent of $J_{0}$, respectively.

This equality is also equivalent to

$J_{0}^{a,b}=\frac{J_{x}-G_{a,\sigma}(x)}{S_{a}(x)}=\frac{1-G_{a,\sigma}^{\prime}(b)}{S_{a}^{\prime}(b)},$
(6)

which suggests another analytic definition of the scale and Gerber-Shiu
function corresponding to an objective $J_{x}$ which involves reflection at
$b$.

The functions $S_{a}(x),G_{a}(x)$ may be shown to stay the same for problems
which require only modifying the boundary condition at $b$, like the problem
of capital injections for the process reflected at $b$, or the problem of
dividends for the process reflected at $b$, with proportional retention
$k_{D}$ (this is in coherence with previously studied problems).

######  COROLLARY 2.

Let us consider the Cramèr-Lundberg setting without diffusion (i.e.
$\sigma=0$), For fixed $k\geq 1,b\geq 0$, the optimality equation
$\frac{\partial}{\partial a}J_{0}^{a,b}=0$ may be written as

$J_{0}^{a,b}=ka-P\Leftrightarrow J_{-a}^{a,b}=-P.$ (7)

######  REMARK 5.

The first equality in (7) provides a relation between the objective $J_{0}$
and the variable $a$; the second recognizes this as the smooth fit equation
$J_{-a}=0$.

Proof: Recalling the expressions of $J_{0}^{a,b}$, $G_{a}(x)$, in (6), in (2),
and from [GK17, Lem. A.4]

$M_{a^{\prime}}(x)=-aC_{a^{\prime}}(x),$

where $C_{a^{\prime}}(x),M_{a^{\prime}}(x)$ denote derivatives with respect to
the subscript $a$, Whenever $b>0$, if $a$ achieves the maximum in
$J^{a,b}_{0}$, it is straightforward (think of the economic interpretation)
that $a$ achieves the maximum of $a\mapsto J^{a,b}_{x}$ for every $x\in[0,b]$.
Therefore, we find

$\displaystyle\begin{aligned} &\frac{\partial}{\partial
a}J_{0}^{a,b}=0\Leftrightarrow
J_{0}^{a,b}=\frac{-G_{a^{\prime}}(x)}{C_{a^{\prime}}(x)}=\frac{-kM_{a^{\prime}}(x)-PC_{a^{\prime}}(x)}{C_{a^{\prime}}(x)}=ka-P\\\
&\Leftrightarrow J_{-a}^{a,b}=J_{0}^{a,b}-ka\Leftrightarrow
J_{-a}^{a,b}=-P.\end{aligned}$

∎

###  2.1 The HJB System

The optimality proof in [AGLW20] is based on showing that the function $J_{x}$
(5) with $a^{*},b^{*}$ defined in (18), (19) is the minimal AC-supersolution
of the HJB system

$\left\\{\begin{split}&\max\left\\{{H}\left(x,{V},{V}^{\prime}(x)\right),1-{V}^{\prime}(x),{V}^{\prime}(x)-k\right\\}=0,\
\forall x\in\mathbb{R}_{+}\\\
&\max\left\\{{V}^{\prime}(x)-k,-P-V(x)\right\\}=0,\ \forall
x\in\mathbb{R}_{-}\end{split}\right.,$ (8)

where the Hamiltonian ${H}$ is given by

${H}\left(x,\phi,v\right):=cv+\lambda\int_{\mathbb{R}_{+}}\phi(x-y)\mu e^{-\mu
y}-\left(q+\lambda\right)\phi(x).$ (9)

To discuss (8), it is useful to introduce the concept of dividend-limited
injections strategies and barrier strategies. The following are also valid for
its generalizations to mixed singular/continuous controls taking into account
reinsurance:

######  DEFINITION 1.

Dividend-limited injections strategies are stationary strategies where the
dividends are paid according to a partition of the state space ${\mathbb{R}}$
in five sets $\mathcal{A},\mathcal{B},\mathcal{C},\mathcal{C}_{0},\mathcal{D}$
as follows:

  1. 1.

If the surplus is in $\mathcal{A}$ (absolute ruin), bankruptcy is declared and
a penalty $P$ is paid;

  2. 2.

If the surplus is in $\mathcal{B}$ bailouts/capital injections are used for
bringing the surplus to the closest point of $\mathcal{C}$.

  3. 3.

If the surplus is in the open set $\mathcal{C}$ (continuation/no action set),
no controls are used.

  4. 4.

If the current surplus is in $\mathcal{C}_{0}\subset\mathcal{D}$ (these are
upper-accumulation points of $\mathcal{C}$), dividends are paid at a positive
rate, in order to keep the surplus process from moving.

  5. 5.

If the current surplus is in $\mathcal{D}$, a positive amount of money is paid
as dividends in order to bring the surplus process to $\mathcal{C}_{0}$.

Barrier strategies are stationary strategies for which
$\mathcal{A},\mathcal{B},\mathcal{C},\mathcal{D}$ are four consecutive
intervals.

######  REMARK 6.

The four sets $\mathcal{A},\mathcal{B},\mathcal{C},\mathcal{D}$ correspond to
the cases when equality in the HJB equation (8) is attained by at least one of
the operators $-V-P,V^{\prime}-k,({\mathcal{G}}-q)V,$ and $1-V^{\prime},$
respectively. Note this generalizes [AM14, Ch. 5.3], where only the last two
operators are considered.

######  REMARK 7.

One may conjecture that dividend-limited injections strategies are of a
(recursive) multi-band nature. In the case of exponential jumps, [AGLW20] show
that the four sets $\mathcal{A},\mathcal{B},\mathcal{C},\mathcal{D}$ in the
optimal solution are intervals, denoted respectively by
$(-\infty,-a),[-a,0],(0,b),[b,\infty)$.

Cheap equity corresponds the case when $\mathcal{C}=\emptyset$, and the
partition reduces to three sets.

##  3 Explicit determination of $a^{*},b^{*}$ when $F(x):=1-e^{-\mu
x},P>-\frac{c}{q}$

In this section we turn to the exponential case, where explicit formulas for
the optimizers $a^{*},b^{*}$ are available. In particular, we will take
advantage of properties of the Lambert-W function, which were not exploited in
[AGLW20]. Subsequently, in sections 5, 6 we will show that exponential
approximations work typically excellently in the general case. Although these
results have already been established in [AGLW20], the present formulations
have two achievements:

  1. 1.

allow an unified formulation of [AGLW20] and [GK17] (via the previously
introduced scale functions);

  2. 2.

make use of a numerical tool (Lambert-W function) to express the optimal
quantities of interest $a^{*},b^{*}$.

###  3.1 The simplified cost function and optimality equations

######  PROPOSITION 3.

Cost function and optimality equations in the exponential case

  1. 1.

$J_{0}^{a,b}=\frac{1-C^{\prime}(b)\left(k\;m(a)+P\overline{F}(a)\right)}{(\overline{F}(a))C^{\prime}(b)+qW_{q}(b)}=\frac{{\gamma(b)}-k\;m(a)-P\overline{F}(a)}{q\theta(b)+\overline{F}(a)},$
(10)

where we put

$\gamma(b)=\frac{1}{C^{\prime}(b)},\theta(b)=\frac{W_{q}(b)}{C^{\prime}(b)}.$

  2. 2.

Put

$j(b):=\frac{\gamma^{\prime}(b)}{q\theta^{\prime}(b)}.$ (11)

For fixed $a\geq 0$, the optimality equation $\frac{\partial}{\partial
b}J_{0}^{a,b}=0$ may be written as

$J_{0}^{a,b}=j(b).$ (12)

  3. 3.

For fixed $k\geq 1$ and $b\geq 0$, at critical points with
$a(b)=a^{(k,P)}(b)\neq 0$ satisfies $\frac{\partial}{\partial
a}J_{0}^{a(b),b}=0$ we must have

$\left[J_{0}^{a,b}-(ka-P)\right]_{a=a(b)}=0.$

Explicitly,

$0=\eta(b,a):=\frac{\gamma(b)}{\theta(b)}-\frac{k}{\mu\theta(b)}F\left(a\right)-q\left(ka-P\right).$
(13)

  4. 4.

When $P\geq-\frac{c}{q}$ and $b\geq 0$ is fixed, the solution of (13) may be
expressed in terms of the principal value of the “Lambert-W(right)" function
(an inverse of $L(z)=ze^{z}$)

$[-e^{-1},\infty)\ni L_{0}(z),z\in[-1,\infty)$

[CGH+96, Boy98, BFS08, Pak15, VLSHGG+19] (this observation is missing in
[AGLW20]).

$\displaystyle(0,\infty)\ni
a(b)=\mu^{-1}\left(-h(b)+L_{0}\left(\frac{e^{h(b)}}{q\theta(b)}\right)\right),\quad
h(b)=h(b,P)=\frac{1}{q\theta(b)}-\frac{\mu}{k}\left(\frac{\gamma(b)}{q\theta(b)}+P\right)$
(14)

It follows that

$J_{0}^{a(b),b}=\frac{k}{\mu}\left(-h(b)+L_{0}\left(\frac{1}{q\theta(b)}e^{h(b)}\right)\right)-P.$
(15)

  5. 5.

In the special case $b=0$, (13) implies that $a=a^{(k,P)}=a^{(k,P)}(0)$
satisfies the simpler equation

$\displaystyle
0=\delta_{k,P}(a):=\lambda\eta(0,a)=\widetilde{c}-k\left(aq+\frac{\lambda}{\mu}(1-e^{-\mu
a})\right),\;\widetilde{c}=c+qP>0,$ (16)

with solution

$\displaystyle\mu\;a^{(k,P)}=-g+L_{0}\left(\frac{\lambda}{q}e^{g}\right)>0,\;g=h(0)=\frac{\lambda}{q}-\frac{\mu}{kq}\widetilde{c}.$
(17)

  6. 6.

At a critical point $(a^{*},b^{*}),a^{*}>0,b^{*}>0$, we must have both
$J_{0}^{a^{*},b^{*}}=j(b^{*})=ka^{*}-P\Longrightarrow$

$a^{*}=s(b^{*}),s(b):=\frac{j(b)+P}{k},$ (18)

and

$0=\eta(b^{*}),\quad\eta(b):=\eta(b,s(b))=\frac{\gamma(b)}{\theta(b)}-qj(b)-\frac{k}{\mu\theta(b)}F\left(\frac{j(b)+P}{k}\right)=0.$
(19)

  7. 7.

The equation $0=\eta(b)$ may be solved explicitly for $P$, yielding

$P=-\frac{k}{\mu}\log\left[1+\frac{q\theta(b)j(b)-\gamma(b)}{\frac{k}{\mu}}\right]-j(b).$
(20)

Proof: 1\. follows from Theorem 1.1.

2\. Let $M(b),N(b)$ denote the numerator an denominator of
$J_{0}^{a,b}:=\frac{M(b)}{N(b)}$ in (10). The optimality equation
$\frac{\partial}{\partial
b}J_{0}^{a,b}=\frac{N^{\prime}(b)}{N(b)}\left(\frac{M^{\prime}(b)}{N^{\prime}(b)}-J_{0}^{a,b}\right)=0$
simplifies to

$J_{0}^{a,b}=\frac{M^{\prime}(b)}{N^{\prime}(b)}=\frac{{\gamma^{\prime}(b)}}{q\theta^{\prime}(b)}=j(b).$

3\. (13) is a consequence of 1 and of the smooth fit result Corollary 2.

4\. See the proof of the particular case 5; $a\in(0,\infty)$ holds since
$P\geq-\frac{c}{q}\Longrightarrow h(b)<\frac{1}{q\theta(b)}.$

5\. (16) follows from $W_{q}(0)=\frac{1}{c},\theta(0)=\lambda^{-1}.$ To get
(17), rewrite the equation (16) as $ze^{z}=\frac{\lambda}{q}e^{g},z=\mu a+g$;
$a\in(0,\infty)$ holds since $P\geq-\frac{c}{q}\Longrightarrow
g<\frac{\lambda}{q}.$

6\. follows from 2. and .3.

7\. is straightforward.

∎

######  REMARK 8.

Note that the de Finetti and Shreve, Lehoczky and Gaver solutions
$a^{*}=a(b^{*})=\begin{cases}0\\\ \infty\end{cases}$ are always non-optimal,
when $P\geq-\frac{c}{q}$ (see (14)).

However, as $k\rightarrow\infty$, $h(b)\rightarrow\frac{1}{q\theta(b)}\neq 0$
and, $a(b)=\mu^{-1}\left(-h(b)+L_{0}\left(h(b)e^{h(b)}\right)\right)=0$. This
suffices to infer that you get de Finetti case.

On the other hand,

${P\to\infty}\Longrightarrow h(b)\to-\infty\Longrightarrow
a(b)\to\infty\Longrightarrow\begin{cases}J_{0}^{a(b),b}\to\frac{\gamma(b)-k/\mu}{q\theta(b)}=\frac{1-kC^{\prime}(b)/\mu}{qW_{q}(b)}=J_{0,SLG}(b)\\\
\eta(b)\to q\left(J_{0,k}^{SLG}(b)-j(b)\right)\end{cases},\forall b>0.$ (21)

Thus, these regimes can be recovered asymptotically. Let now
$b_{k}^{*,S},b_{P}^{*,D}$ denote the unique roots of $\eta(b)=0$ in the two
asymptotic cases, which coincide with the classic Shreve, Lehoczky and Gaver
and de Finetti barriers.

Then, it may be checked that $b^{*}\leq\min[b_{k}^{*,S},b_{P}^{*,D}]$.

###  3.2 Existence of the roots of the equations $\eta(b)=0,\delta_{k,P}=0$

The following (new) result discusses the existence of the roots of the
equations $\eta(b)=0,\delta_{k,P}=0$ introduced in proposition (3) and relates
them to the Lambert-W function.

######  PROPOSITION 4.

  1. 1.

$\theta$ increases from $\theta(0)=\frac{1/c}{\lambda/c}=\frac{1}{\lambda}$ to
$\theta(\infty)=\frac{1}{c\Phi_{q}-q}$, as we see it in the figure below.

Figure 2: Plot of $\theta$ with $\theta(0)=2$ and $\theta(\infty)=22.8743$,
for $\mu=2,\ c=3/4,\ \lambda=1/2,\ q=1/10$, $P=1$ and $k=3/2$.

$\gamma$ is increasing-decreasing (from $\frac{c}{\lambda}$ to $0$), with a
maximum at the unique root of $C^{\prime\prime}(x)=0$ given by

$\bar{b}:=\frac{1}{\Phi_{q}-\rho_{-}}\log\left(\frac{\rho_{-}^{2}}{\Phi_{q}^{2}}\right),$
(22)

where $\Phi_{q},\rho_{-}$ denote the positive and negative roots of the
Cramèr-Lundberg equation $\kappa(s)=0$.

The figure below illustrates the plot of the function $\gamma$ and $j({b})$ in
which the $\bar{b}$ is represented by the black point.

Figure 3: Plots of $j({b})$ and $\gamma({b})$ with $\bar{b}=2.5046$ and
$j(0)=4.5$, for $\mu=2,\ c=3/4,\ \lambda=1/2,\ q=1/10$, $P=1$ and $k=3/2$.

If $c\mu-\left(q+\lambda\right)>0$, then $\bar{b}>0$ defined in (22) is the
unique positive root of $j({b})$ and
$\eta\left(\bar{b}\right)=\frac{1}{W_{q}\left(\bar{b}\right)}>0$. See the
figure (4)

Figure 4: For $\mu=2,\ c=3/4,\ \lambda=1/2,\ q=1/10$, $P=1$ and $k=3/2$, the
root of $\eta({b})=0$ is at $b=0.469843$

The function $j(b)=\frac{\gamma^{\prime}(b)}{q\theta^{\prime}(b)}$ is
nonnegative and decreasing to $0$ on $[0,\bar{b}]$, with

$j(0)=\frac{\lambda}{\mu
q}\frac{-C^{\prime\prime}(0)}{(C^{\prime}(0))^{2}}=\frac{c\mu-\left(q+\lambda\right)}{\mu
q}.$ (23)

See the figure (3).

  2. 2.

Put

$\delta_{k,P}:=\delta_{k,P}(a^{(k,P)})=\delta_{k,P}(\frac{j(0)+P}{k})=\delta_{k,P}(\frac{\widetilde{c}\mu-\left(q+\lambda\right)}{k\mu
q})=\frac{\lambda+q-\lambda
k\left(1-e^{-\frac{\widetilde{c}\mu-\lambda-q}{qk}}\right)}{\mu},$ (24)

and assume

${\underset{k\rightarrow\infty}{\lim}\delta_{k,P}=\frac{\lambda+q}{\mu}-\frac{\lambda\left(\mu\widetilde{c}-(\lambda+q)\right)}{q\mu}=\frac{\left(\lambda+q\right)^{2}-\lambda\mu\widetilde{c}}{q\mu}<0}\Leftrightarrow\widetilde{c}\mu>\lambda^{-1}\left(\lambda+q\right)^{2}.$
(25)

Then, $\forall P>-\frac{c}{q}$, the function $\delta_{k,P}$ is decreasing in
$k$ with $\delta_{1,P}>0$, and has a unique root

$k_{c}=k_{c}(P):=\frac{q+\lambda}{\lambda}\frac{f}{f+L_{0}\left(-fe^{-f}\right)}>\frac{q+\lambda}{\lambda},$
(26)

where

$f:=\frac{\lambda}{q+\lambda}\frac{\widetilde{c}\mu-\left(\lambda+q\right)}{q}>1\Leftrightarrow\widetilde{c}\mu>\lambda^{-1}\left(\lambda+q\right)^{2}\Leftrightarrow
P>P_{l}:=q^{-1}\left(\mu^{-1}\lambda^{-1}\left(\lambda+q\right)^{2}-c\right)$
(27)

(note that the denominator $f+L_{0}\left(-fe^{-f}\right)$ does not equal $0$
since $f>1$ and $L_{0}$ takes always values bigger than $-1$; or, note that
$-f=L_{-1}\left(L(-f)\right)$, where $L_{-1}$ is the other real branch of the
Lambert function).

Furthermore,

$\delta_{k,P}<0\Leftrightarrow k>k_{c}(P).$ (28)

  3. 3.

It follows that $\eta(b)=0$ has at least one solution of in $(0,\bar{b}]$ iff

$\eta(0)=\frac{c}{\lambda}-\frac{1}{\lambda}\left(c-\frac{q+\lambda}{\mu}\right)-\frac{k}{\mu}F\left(a^{(k,P)}\right)=\frac{1}{\lambda\mu}\left(\lambda+{q}-\lambda
kF\left(a^{(k,P)}\right)\right)=\frac{1}{\lambda}\delta_{k,P}<0\Leftrightarrow
k>k_{c}.$ (29)

The first such solution will be denoted by $b^{*}$.

Proof: For 1. see [AGLW20, Proof of Theorem 11, A2].

2\. By using the assumption $\widetilde{c}\mu>\lambda^{-1}(\lambda+q)^{2}$ we
get $\widetilde{c}\mu\geq\lambda+q$, and
$k\in[1,\infty)\rightarrow\delta_{k,P}$ is decreasing.

Put $d=\frac{\widetilde{c}\mu-\left(\lambda+q\right)}{q}$. The inequality
$\delta_{k,P}<0$ (see (29)) may be reduced to

$\displaystyle e^{-\frac{d}{k}}<1-\frac{q+\lambda}{\lambda k}\Leftrightarrow
1<e^{\frac{d}{k}}\left(1-\frac{(q+\lambda)}{d\lambda}\frac{d}{k}\right):=e^{z}\left(1-z/f\right).$

Rewriting the latter as $-f>e^{z}\left(z-f\right)$ we recognize, by putting
$z=y+f$, an inequality reducible to $ye^{y}<-fe^{-f}.$ The solution is

$y<L_{0}\left(-fe^{-f}\right),$

where $L_{0}$ is the principal branch of the Lambert-W function.

The final solution is (28), where we may note that the variables $k,P$ have
been separated.

3\. is straightforward.

∎

######  REMARK 9.

The function $[1,\infty)\ni
f\mapsto\frac{f}{f+L_{0}\left(-fe^{-f}\right)}\in(1,\infty)$ blows up at
$f=1$, and converges to $1$ when $f\to\infty$ (or when either $\mu$ or
$\widetilde{c}=c+qP$ are large enough) as may be noticed in the figure below,
which blows up at the value $P_{l}:=-4/5.$ Note also that when $f$ (or one of
$c,P,\mu$ are large enough), $k_{c}$ given by (26) stabilizes to the
equilibrium $\frac{\lambda+q}{\lambda}=6/5$; this is related to [APP07, Lemma
2], [KS08, Lemma 7], who obtain the same condition for $b^{*}=0$ (without
buffering capital injections). Intuitively, under these conditions, buffering
is not crucial.

At the other end, as $f$ tends to its lower limit and to the regime B, the
notion of equity expensiveness vanishes, and $k_{c}\to\infty$.

Figure 5:  $k_{c}$ as function of $P$, for several values of $c$, with the
vertical asymptote at $P_{l}$ fixed.

The next two figures illustrate how $k_{c}$ blows up at the critical values
$q_{l}:=(1/2)(-2\lambda+P\lambda\mu+\sqrt{\lambda}\sqrt{\mu}\sqrt{4}c-4P\lambda+P^{2}\lambda\mu)$
and $\lambda_{l}:=\left(c\mu-\sqrt{(-c\mu+\mu(-P)q+2q)^{2}-4q^{2}}+\mu
Pq-2q\right)$ (represented by red points in the figures below). The dark
(blue) parts correspond to the regime $A$.

(a) $k_{c}$ defined in (26)as function of $q$, for $\mu=2,\ c=3/2,\ \lambda=1$
and $P=1$; $q_{l}=\sqrt{2}$.

(b) $k_{c}$ as function of $\lambda$, for $\mu=2,\ c=3/2,\ q=1/10$ and $P=1$.

Figure 6: $k_{c}$ as a function of $q$, $\lambda$

##  4 Which exponential approximation?

###  4.1 Three de Vylder-type exponential approximations for the ruin
probability

In the simplest case of exponential jumps of rate $\mu$ and $\sigma=0$, the
formula for the ruin probability is

$\Psi(x)=P_{x}[\exists t\geq 0:X_{t}<0]=\frac{1}{1+{\theta}}\
exp\left(-\frac{x{\theta}\mu}{1+{\theta}}\right)=\frac{1}{1+{\theta}}\
exp\left(-\frac{x{\theta}m_{1}^{-1}}{1+{\theta}}\right),$ (30)

where $\theta=\frac{c-\lambda m_{1}}{\lambda m_{1}}$ is the loading
coefficient. By plugging the correct mean of the claims in the second formula
yields the simplest approximation for processes with finite mean claims.

More sophisticated is the Renyi exponential approximation

$\Psi_{R}(x)=\frac{1}{1+{\theta}}\
exp\left(-\frac{x{\theta}\widehat{m}_{2}^{-1}}{1+{\theta}}\right),\widehat{m}_{2}=\frac{m_{2}}{2m_{1}};$
(31)

This formula can be obtained as a two point Padé approximation of the Laplace
transform, which conserves also the value $\Psi(0)=(1+\theta)^{-1}$ [AP14]. It
may be also derived heuristically from the first formula in (30), via
replacing $\mu$ by the correct “excess mean" of the excess/severity density

$f_{e}(x)=\frac{\overline{F}(x)}{m_{1}}=\frac{1-F(x)}{m_{1}},$

which is known to be $\widehat{m}_{2}$. Heuristically, it makes more sense to
approximate $f_{e}(x)$ instead of the original density $f(x)$, since
$f_{e}(x)$ is a monotone function, and also an important component of the
Pollaczek-Khinchine formula for the Laplace transform
$\widehat{\Psi}(s)=\int_{0}^{\infty}e^{-sx}F(dx)$ – see [Ram92, AP14].

More moments are put to work in the de Vylder approximation

$\Psi_{DV}(x)=\frac{1}{1+\widetilde{\theta}}\
exp\left(-\frac{x\widetilde{\theta}\widehat{m}_{3}^{-1}}{1+\widetilde{\theta}}\right),\;\widehat{m}_{3}:=\frac{m_{3}}{3m_{2}},\;\widetilde{\lambda}=\frac{9m_{2}^{3}}{2m_{3}^{2}}\lambda,\;\;\widetilde{c}=c-\lambda
m_{1}+\widetilde{\lambda}\widehat{m}_{3},\;\widetilde{\theta}=\frac{2m_{1}m_{3}}{3m_{2}^{2}}\theta=\frac{\widehat{m}_{3}}{\widehat{m}_{2}}\theta.$
(32)

Interestingly, the result may be expressed in terms of the so-called
"normalized moments"

${\widehat{m}_{i}}=\frac{m_{i}}{i\;m_{i-1}}$ (33)

introduced in [BHT05].

The de Vylder approximation parameters above may be obtained either from

  1. 1.

equating the first three cumulants of our process to those of a process with
exponentially distributed claim sizes of mean $\widehat{m}_{3}$, and modified
$\lambda,c$ [DV78] (however $p=c-\lambda m_{1}=E_{0}[X_{1}]$ must be
conserved, since this is the first cumulant), or

  2. 2.

a Padé approximation of the Laplace transform of the ruin probabilities
[ACFH11].

The second derivation via Padé shows that higher order approximations may be
easily obtained as well. They might not be admissible, due to negative values,
but packages for “repairing" the non-admissibility are available – see for
example [DcSA16].

The first derivation of the de Vylder approximation is a process approximation
(i.e., independent of the problem considered); as such, it may be applied to
other functionals of interest besides ruin probabilities($W_{q}(x)$, dividend
barriers, etc), simply by plugging the modified parameters in the exact
formula for the ruin probability of the simpler process.

###  4.2 Three two point Padé approximations of the Laplace transform
$\widehat{W}_{q}$ of scale function

The simplest approximations for the scale function $W_{q}(x)$ will now be
derived heuristically from the following example.

######  EXAMPLE 2.

The Cramér-Lundberg model with exponential jumps  Consider the Cramér-Lundberg
model with exponential jump sizes with mean $1/\mu$, jump rate $\lambda$,
premium rate $c>0$, and Laplace exponent
$\kappa(s)=s\left(c-\frac{\lambda}{\mu+s}\right)$. Solving
$\kappa(s)-q=0\Leftrightarrow cs^{2}+s(c\mu-\lambda-q)-q\mu=0$ for $s$ yields
two distinct solutions $\gamma_{2}\leq 0\leq\gamma_{1}=\Phi_{q}$ given by

$\displaystyle\gamma_{1}(\mu,\lambda,c)=\gamma_{1}=$
$\displaystyle\frac{1}{2c}\left(-\left(\mu c-\lambda-q\right)+\sqrt{\left(\mu
c-\lambda-q\right)^{2}+4\mu qc}\right),$
$\displaystyle\gamma_{2}(\mu,\lambda,c)=\gamma_{2}=$
$\displaystyle\frac{1}{2c}\left(-\left(\mu c-\lambda-q\right)-\sqrt{\left(\mu
c-\lambda-q\right)^{2}+4\mu qc}\right).$

The $W$ scale function is:

$W_{q}(x)=\frac{A_{1}e^{\gamma_{1}x}-A_{2}e^{\gamma_{2}x}}{c(\gamma_{1}-\gamma_{2})}\Leftrightarrow\widehat{W}_{q}(s)=\frac{s+\mu}{cs^{2}+s(c\mu-\lambda-q)-q\mu},$
(34)

where $A_{1}=\mu+\gamma_{1},A_{2}=\mu+\gamma_{2}$.

Furthermore, it is well-known and easy to check that the function
$W_{q}^{\prime}(x)$ is in this case unimodal with global minimum at

$b_{DeF}=\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})}\quad&\text{if
$W_{q}^{\prime\prime}(0)<0\Leftrightarrow(q+\lambda)^{2}$}-c\lambda\mu<0\\\
0&\text{if $W_{q}^{\prime\prime}(0)\geq
0\Leftrightarrow(q+\lambda)^{2}-c\lambda\mu\geq 0$}\end{cases},$ (35)

since
$W_{q}^{\prime\prime}(0)=\frac{(\gamma_{1})^{2}(\mu+\gamma_{1})-(\gamma_{2})^{2}(\mu+\gamma_{2})}{c(\gamma_{1}-\gamma_{2})}=\frac{(q+\lambda)^{2}-c\lambda\mu}{c^{3}}$
and that the optimal strategy for the de Finetti problem is the barrier
strategy at level $b_{DeF}$ (see for example [APP07], [AGVA19, Sec. 3]).

Plugging now the respective parameters of the de Vylder type approximations in
the exact formula (34) for the Cramèr-Lundberg process with exponential
claims, we obtain three approximations for $\widehat{W}_{q}$:

  1. 1.

“Naive exponential" approximation obtained by plugging $\mu^{-1}\rightarrow
m_{1}$ in (34) (as was done, for a different purpose) in (30)

  2. 2.

Renyi 444This is called DeVylder B) method in [GSS08, (5.6-5.7)], since it is
the result of fitting the first two cumulants of the risk process. , obtained
by plugging
$\mu^{-1}\rightarrow\widehat{m}_{2},\lambda_{R}\rightarrow\lambda\frac{m_{1}}{\widehat{m}_{2}}$
(since $c$ is unchanged, the latter equation is equivalent to the conservation
of $\rho=\frac{\lambda m_{1}}{c},$ and to the conservation of $\theta$, so
this coincides with the Renyi ruin approximation used in (31).)

  3. 3.

De Vylder, obtained by plugging
$\mu^{-1}\rightarrow\widehat{m}_{3},\widetilde{\lambda}\rightarrow\lambda\frac{9m_{2}^{3}}{2m_{3}^{2}},\widetilde{c}=c-\lambda
m_{1}+\widetilde{\lambda}\widehat{m}_{3}$.

######  REMARK 10.

In the case of exponential claims, these three approximations are exact, by
definition (or check that for exponential claims all the normalized moments
are equal to $\mu^{-1}$).

######  REMARK 11.

The conditions for the non-negativity of the barrier is
$W_{q}^{\prime\prime}(0_{+})<0\Leftrightarrow(\frac{\lambda+q}{c})^{2}<\frac{\lambda}{c}f(0)$.
Here, this condition is satisfied for the exact when
$\theta>\frac{(\lambda+q)^{2}(1-\rho)}{\lambda^{2}f(0)m_{1}}$.

It is shown in [AHPS19, Prop. 1] that the three de Vylder type approximations
are two-point Padé approximations of the Laplace transform (hence higher order
generalizations are immediately available).

We recall that two-point Padé approximations incorporate into the Padé
approximation two initial values of the function (which can be derived easily
via the initial value theorem, from the Pollaczek-Khinchine Laplace
transform):

$\displaystyle W_{q}(0_{+})=\lim_{s\to\infty}s\widehat{W}_{q}(s)=\frac{1}{c},$
(36) $\displaystyle
W_{q}^{\prime}(0_{+})=\lim_{s\to\infty}s\left(\frac{s}{\kappa(s)-q}-W_{q}(0_{+})\right)=\frac{q+\lambda}{c^{2}}.$
(37)

In our case, incorporating both $W_{q}(0_{+}),W_{q}^{\prime}(0_{+})$ leads to
the natural exponential approximation which is therefore the best near $x=0$.
Incorporating none of them yields the de Vylder approximation, which is the
best asymptotically. Incorporating only $W_{q}(0_{+})$ leads to Renyi, which
is expected to be the best in an intermediate regime.

Note that when the jump distribution has a density $f$, it holds that :
666This equation is important in establishing the nonnegativity of the optimal
dividends barrier.

$\displaystyle W_{q}^{\prime\prime}(0_{+})$
$\displaystyle=\lim_{s\to\infty}s\left(s\left(\frac{s}{\kappa(s)-q}-W_{q}(0_{+})\right)-W_{q}^{\prime}(0_{+})\right)=\frac{1}{c}\Big{(}(\frac{\lambda+q}{c})^{2}-\frac{\lambda}{c}f(0)\Big{)}.$
(38)

Thus, $W_{q}^{\prime\prime}(0)$ already requires knowing $f_{C}(0)$ (which is
a rather delicate task starting from real data); therefore we will not
incorporate into the Padé approximation more than two initial values of the
function.

We recall below in Proposition 5 three types of two-point Padé approximations
[AHPS19, Prop. 1], and particularize them to the case when the denominator
degree is $n=2$ (which are further illustrated below).

######  PROPOSITION 5.

Three matrix exponential approximations for the scale function.

  1. 1.

To secure both the values of $W_{q}(0)$ and $W_{q}^{\prime}(0)$, take into
account (36) and (37), i.e. use the Padé approximation

$\widehat{W}_{q}(s)\sim\frac{\sum_{i=0}^{n-1}a_{i}s^{i}}{cs^{n}+\sum_{i=0}^{n-1}b_{i}s^{i}},a_{n-1}=1,b_{n-1}=ca_{n-2}-\lambda-q.$

For $n=2$ we recover the “natural exponential " approximation of plugging
$\mu\to\frac{1}{m_{1}}$ in (34):

$\widehat{W}_{q}(s)\sim\frac{\frac{1}{m_{1}}+s}{cs^{2}+s\left(\frac{c}{m_{1}}-\lambda-q\right)-\frac{q}{m_{1}}},$
(39)

used also (for a different purpose) in (30).

  2. 2.

To ensure only $W_{q}(0)=\frac{1}{c}$, we must use the Padé approximation

$\widehat{W}_{q}(s)\sim\frac{\sum_{i=0}^{n-1}a_{i}s^{i}}{cs^{n}+\sum_{i=0}^{n-1}b_{i}s^{i}},a_{n-1}=1.$

For $n=2$, we find

$\widehat{W}_{q}(s)\sim\frac{\frac{2m_{1}}{m_{2}}+s}{cs^{2}+\frac{s\left(2cm_{1}-2\lambda
m_{1}^{2}-m_{2}q\right)}{m_{2}}-\frac{2m_{1}q}{m_{2}}}=\frac{\frac{1}{\widehat{m}_{2}}+s}{cs^{2}+{s\left(\frac{c}{\widehat{m}_{2}}-\lambda\frac{m_{1}}{\widehat{m}_{2}}-q\right)}-\frac{q}{\widehat{m}_{2}}},$
(40)

where $\widehat{m}_{2}=\widehat{m}_{2}=\frac{m_{2}}{2m_{1}}$ is the first
moment of the excess density $f_{e}(x)$. Note that it equals the scale
function of a process with exponential claims of rate $\widehat{m}_{2}^{-1}$
and with $\lambda$ modified to
$\lambda_{R}=\lambda\frac{m_{1}}{\widehat{m}_{2}}$. Since $c$ is unchanged,
the latter equation is equivalent to the conservation of $\rho=\frac{\lambda
m_{1}}{c},$ and to the conservation of $\theta$, so this coincides with the
Renyi approximation 444This is called DeVylder B) method in [GSS08,
(5.6-5.7)], since it is the result of fitting the first two cumulants of the
risk process. used in (31).

  3. 3.

The pure Padé approximation yields for $n=2$

$\displaystyle\widehat{W}_{q}(s)\sim\frac{s+\frac{3m_{2}}{m_{3}}}{s^{2}\left(c-\lambda
m_{1}+\lambda\frac{3m_{2}^{2}}{2m_{3}}\right)+s\left(c\frac{3m_{2}}{m_{3}}-\frac{3m_{1}m_{2}}{m_{3}}\lambda-q\right)-\frac{3m_{2}}{m_{3}}q}$
$\displaystyle=\frac{s+\frac{1}{\widehat{m}_{3}}}{\widetilde{c}s^{2}+s\left(\widetilde{c}\frac{1}{\widehat{m}_{3}}-\widetilde{\lambda}-q\right)-\frac{1}{\widehat{m}_{3}}q},\;\widetilde{c}=c-\lambda
m_{1}+\widetilde{\lambda}\widehat{m}_{3},\widetilde{\lambda}=\lambda\frac{9m_{2}^{3}}{2m_{3}^{2}}.$

Note that both the coefficient of $s^{2}$ in the denominator coincides with
the coefficient $\widetilde{c}$ in the classic de Vylder approximation, since
$\widetilde{\lambda}\widehat{m}_{3}=\lambda\frac{9m_{2}^{3}}{2m_{3}^{2}}\frac{m_{3}}{3m_{2}}=\lambda\frac{3m_{2}^{2}}{2m_{3}}$,
and so does the coefficient of $s$, since

$c\frac{3m_{2}}{m_{3}}-\frac{3m_{1}m_{2}}{m_{3}}\lambda=\widetilde{c}\frac{1}{\widehat{m}_{3}}-\widetilde{\lambda}=\left(c-\lambda
m_{1}+\widetilde{\lambda}\widehat{m}_{3}\right)\frac{1}{\widehat{m}_{3}}-\widetilde{\lambda}.$

##  5 Examples of computations involving scale function and dividend value
approximations

Our goal in this section is to investigate whether exponential approximations
are precise enough to yield reasonable estimates for quantities important in
control like

  1. 1.

the dominant exponent $\Phi_{q}$ of $W_{q}(x)$

  2. 2.

the last local minimum of $W_{q}^{\prime}(x)$, $b_{DeF}$, which yields, when
being the global minimum, the optimal De Finetti barrier

  3. 3.

$W^{\prime\prime}_{q}(0)$, which determines if $b_{DeF}=0$

  4. 4.

the functional $J_{0}$ yielding the maximum dividends with capital injections.

All the examples considered involve a Cramèr-Lundberg model with rational
Laplace transform $\widehat{W}_{q}(s)$ (since in this case, the computation of
$W_{q},Z_{q}$ is fast and in principle arbitrarily large precision may be
achieved with symbolic algebra systems).

  1. 1.

For the first three problems, we will use de Vylder type approximations.
Graphs of $W^{\prime}_{q}$, $W^{\prime\prime}_{q}$ and some tables summarizing
the simulation results will be presented. We note that in most of the cases
that we observed, the de Vylder approximation of $\Phi_{q}$ deviates from the
exact value the least – see for example Table 2. For the De Finetti barrier,
the "winner" depends on the size of $b_{DeF}$. Unsurprisingly, when near $0$,
the natural exponential approximation wins, and as $b_{DeF}$ increases, Renyi
and subsequently the de Vylder approximation take the upper hand – see for
example Table 3.

  2. 2.

For the computation of $J_{0}$, we provide, besides the exact value, also two
approximations:

    1. (a)

For a given density of claims $f$ one computes an exponential density
approximation $f_{e}(x)=\frac{1}{m_{1}}\text{exp}(-\frac{x}{m_{1}})$ where
$m_{1}$ is the first moment of $f$. Subsequently, $W$, $Z$, $J_{0}$ and $a,b$
are obtained using the exponential approximation $f_{e}$. Quantities obtained
by this method would be referred to with an affix ‘expo pure’.

    2. (b)

For a given density of claims $f$, the value function is computed via the
formula which assumes exponential claims in equation 3, but the "ingredients"
$W$, $Z,\overline{F}$ and the mean function $m$ are the correct ingredients
corresponding to our original density $f$. Quantities obtained by this method
would be referred to with an affix ‘expo CI’.

It turns out that the pure expo approximation works better for large $\theta$,
and the correct ingredients approximation works better for small $\theta$.

Note that we only included tables illustrating approximating $J_{0}$ for the
first two examples, to keep the length of the paper under control, but similar
results were obtained for the other examples.

###  5.1 A Cramér-Lundberg process with hyperexponential claims of order 2

We take a look at a Cramér-Lundberg process with density function
$f(x)=\frac{2}{3}e^{-x}+\frac{2}{3}e^{-2x}$ with $\lambda=1$, $\theta=1$ and
$q=\frac{1}{10}$.

(a) $W^{\prime}_{q}(x)$

(b) $W^{\prime\prime}_{q}(x)$

Figure 7: Exact and approximate plots of $W^{\prime}_{q}(x)$ and $W^{\prime\prime}_{q}(x)$ for $f(x)=\frac{2}{3}e^{-x}+\frac{2}{3}e^{-2x}$, $\theta=1$, $q=\frac{1}{10}$. |  | Dominant exponent
---
$\Phi_{q}$
| Percent relative error
---
($\Phi_{q}$)
| Optimal barrier
---
$b_{DeF}$
| Percent relative error
---
($b_{DeF}$)
Exact | 0.110113 | 0 | 3.45398 | 0
Expo | 0.110657 | 0.494313 | 3.51173 | 1.67191
Dev | 0.110115 | 0.00195933 | 3.48756 | 0.972251
Renyi | 0.110078 | 0.0321413 | 3.5323 | 2.26744
Table 1: Exact and approximate values of $\Phi_{q}$ and $b_{DeF}$ for $f(x)=\frac{2}{3}e^{-x}+\frac{2}{3}e^{-2x}$, $theta=1$, $q=\frac{1}{10}$, as well as percent relative errors, computed as the absolute value of the difference between the approximation and the exact, divided by the exact, times 100. Relative errors for $\Phi_{q}$ are less than $0.5\%$, with the pure exponential approximation proving to be the worst and the DeVylder the best approximations, respectively. The optimal barrier $b_{DeF}$ is also best approximated by DeVylder, with Renyi being the worst at $2.26\%$ . $\theta$ | Closest approximation |  $\Phi_{q}$ exact |  $\Phi_{q}$ approximation | % error $\Phi_{q}$
---|---|---|---|---
1 | Dev | 0.110113 | 0.110115 | 0.00195933
0.9 | Dev | 0.120328 | 0.120331 | 0.00269878
0.8 | Dev | 0.132452 | 0.132457 | 0.00380056
0.7 | Dev | 0.147017 | 0.147025 | 0.00548411
0.6 | Dev | 0.16475 | 0.164763 | 0.00812643
0.5 | Dev | 0.186652 | 0.186675 | 0.0123901
0.4 | Dev | 0.214122 | 0.214163 | 0.0194631
0.3 | Dev | 0.249118 | 0.249196 | 0.0315039
0.2 | Dev | 0.294396 | 0.294551 | 0.0524528
0.1 | Dev | 0.353829 | 0.354145 | 0.0894466
Table 2: Exact and the winning DeVylder approximate values of $\Phi_{q}$ for $f(x)=\frac{2}{3}e^{-x}+\frac{2}{3}e^{-2x}$ when $\theta$ varies. $\theta$ | Closest approximation | Barrier exact | Barrier approx | % error Barrier
---|---|---|---|---
1 | Dev | 3.45398 | 3.48756 | 0.972251
0.9 | Dev | 3.20191 | 3.23103 | 0.909487
0.8 | Dev | 2.90951 | 2.93074 | 0.729628
0.7 | Dev | 2.57043 | 2.57742 | 0.272088
0.6 | Dev | 2.1804 | 2.16054 | 0.910666
0.5 | Ren | 1.74216 | 1.75266 | 0.60278
0.4 | Expo | 1.2735 | 1.29456 | 1.65378
0.3 | Expo | 0.81068 | 0.652264 | 19.5412
Table 3: Exact and approximate values of $b_{DeF}$ for for $f(x)=\frac{2}{3}e^{-x}+\frac{2}{3}e^{-2x}$ when $\theta$ varies. As $b_{DeF}$ approaches $0$, errors of all the approximations increase dramatically, with the pure exponential approximation performing better than the rest. Meanwhile, as $b_{DeF}$ increases, Renyi and subsequently the de Vylder approximation take the upper hand. For $\theta=0.2$ and $\theta=0.1$, all the approximations yield a 0 barrier approximate for exact barrier values of $0.392105$ and $0.0354538$ respectively, hence failing to predict the non-zero barrier. | J0
---|---
$\theta$ | J0 exact | J0 expo pure | J0 expo pure error | J0 expo CI | J0 expo CI error
1 | 5.95034 | 5.99151 | 0.691856 | 6.26009 | 5.20551
0.9 | 5.15579 | 5.17573 | 0.386663 | 5.45269 | 5.7584
0.8 | 4.39383 | 4.38494 | 0.202205 | 4.67042 | 6.29489
0.7 | 3.68299 | 3.63933 | 1.18555 | 3.92937 | 6.68958
0.6 | 3.04577 | 2.96728 | 2.57704 | 3.25112 | 6.7423
0.5 | 2.50331 | 2.39942 | 4.15022 | 2.65901 | 6.21974
0.4 | 2.06833 | 1.9585 | 5.31006 | 2.17044 | 4.93725
0.3 | 1.74095 | 1.65616 | 4.86984 | 1.78984 | 2.80878
0.2 | 1.50439 | 1.44242 | 4.11969 | 1.50871 | 0.286672
0.1 | 1.30271 | 1.25324 | 3.79748 | 1.30271 | 0
Table 4: Values of $J_{0}$ compared with approximations using all exponential inputs ($J_{0}$ expo pure) and actual inputs but computed using the exponential formula ($J_{0}$ expo CI). The pure exponential approximation does a good job of approximating $J_{0}$ for higher values of $\theta$ considered, while the exponential CI approximation seemed to fair better for lower $\theta$ values | a
---|---
$\theta$ | a exact | a expo pure | a expo pure error | a expo CI | a expo CI error
1 | 3.9669 | 3.99434 | 0.691861 | 4.17339 | 5.20551
0.9 | 3.4372 | 3.45049 | 0.386665 | 3.63512 | 5.7584
0.8 | 2.92922 | 2.9233 | 0.202204 | 3.11361 | 6.29489
0.7 | 2.45533 | 2.42622 | 1.18555 | 2.61958 | 6.68958
0.6 | 2.03051 | 1.97818 | 2.57704 | 2.16741 | 6.7423
0.5 | 1.66888 | 1.59961 | 4.15022 | 1.77268 | 6.21974
0.4 | 1.37888 | 1.30566 | 5.31006 | 1.44696 | 4.93725
0.3 | 1.16063 | 1.10411 | 4.86983 | 1.19323 | 2.80878
0.2 | 1.00293 | 0.961612 | 4.11969 | 1.0058 | 0.286672
0.1 | 0.868476 | 0.835496 | 3.79748 | 0.868476 | 0
Table 5: Values of $a$ compared with approximations using all exponential inputs ($a$ expo pure) and actual inputs but computed using the exponential formula ($a$ expo CI) | b
---|---
$\theta$ | b exact | b expo pure | b expo pure error | b expo CI | b expo CI error
1 | 1.41036 | 1.46188 | 3.65293 | 1.25374 | 11.1045
0.9 | 1.37645 | 1.44439 | 4.93621 | 1.23362 | 10.3761
0.8 | 1.31492 | 1.40417 | 6.78809 | 1.19529 | 9.09781
0.7 | 1.21057 | 1.32258 | 9.25207 | 1.12775 | 6.84178
0.6 | 1.04634 | 1.17215 | 12.0245 | 1.01753 | 2.7529
0.5 | 0.810767 | 0.920406 | 13.5229 | 0.853397 | 5.25805
0.4 | 0.510085 | 0.538725 | 5.61475 | 0.634716 | 24.4335
0.3 | 0.17425 | 0.0105496 | 93.9457 | 0.376872 | 116.282
0.2 | 0 | 0 | 0 | 0.105322 | 100
0.1 | 0 | 0 | 0 | 0 | 0
Table 6: Values of $b$ compared with approximations using all exponential
inputs ($b$ expo pure) and actual inputs but computed using the exponential
formula ($b$ expo CI)

###  5.2 A Cramér-Lundberg process with hyperexponential claims of order 3

Consider a Cramér-Lundberg process with density function
$f(x)=\frac{12}{83}e^{-x}+\frac{42}{83}e^{-2x}+\frac{150}{83}e^{-3x}$, and
$c=1$, $\lambda=\frac{83}{48}$, $\theta=\frac{263}{235}$, $p=\frac{263}{498}$,
$q=\frac{5}{48}$.

The Laplace exponent of this process is
$\kappa(s)=s-\frac{12s}{83(s+1)}-\frac{21s}{83(s+2)}-\frac{50s}{83(s+3)}$ and
from this one can invert $\frac{1}{\kappa(s)-q}=\widehat{W}_{q}(s)$ to obtain
the scale function 222Laplace inversion done via Mathematica; coefficients and
exponents are decimal approximations of the real values.

$\displaystyle W_{q}(x)$
$\displaystyle=-0.0813294e^{(-2.60997x)}-0.179472e^{(-1.68854x)}-0.373887e^{(-0.779311x)}+1.63469e^{(0.18198x)}.$

From this, we see that the dominant exponent is $\Phi_{q}=0.18198$.

Figure 8 shows the exact and approximate plots of the first two derivatives of
$W_{q}$. The exact plots are labelled Wxexact, and coloured as the darkest.
The plots of $W_{q}^{\prime}$ exhibit noticeable unique minima around $x=2$,
with the exact one being at $b_{DeF}=1.89732$, which is the optimal barrier
that maximizes dividends here. Note that the approximations are practically
indistinguishable from the exact around this point (which is our main object
of interest here).

(a) $W^{\prime}_{q}(x)$

(b) $W^{\prime\prime}_{q}(x)$

Figure 8: Exact and approximate plots of $W^{\prime}_{q}(x)$ and $W^{\prime\prime}_{q}(x)$ for $f(x)=\frac{12}{83}e^{-x}+\frac{42}{83}e^{-2x}+\frac{150}{83}e^{-3x}$, $c=1$, $q=\frac{5}{48}$. |  | Dominant exponent
---
$\Phi_{q}$
| Percent relative error
---
($\Phi_{q}$)
| Optimal barrier
---
$b_{DeF}$
| Percent relative error
---
($b_{DeF}$)
Exact | 0.18198 | 0 | 1.89732 | 0
Expo | 0.184095 | 1.162215628 | 2.04608 | 7.840532962
Renyi | 0.181708 | 0.149466974 | 2.08136 | 9.699997892
Dev | 0.182011 | 0.017034839 | 1.91233 | 0.79111589
Table 7: Exact and approximate values of $\Phi_{q}$ and $b_{DeF}$ for $f(x)=\frac{12}{83}e^{-x}+\frac{42}{83}e^{-2x}+\frac{150}{83}e^{-3x}$, $c=1$, $q=\frac{5}{48}$, as well as percent relative errors, computed as the absolute value of the difference between the approximation and the exact, divided by the exact, times 100. Relative errors for the $\Phi_{q}$ value are less than $1.25\%$, even for the worse natural exponential approximation, and the DeVylder approximation is the winner. The optimal barrier $b_{DeF}$ is well approximated only by DeVylder. $\theta$ | Closest approximation |  $\Phi_{q}$ exact |  $\Phi_{q}$ approximation | % error $\Phi_{q}$
---|---|---|---|---
263/235 | Dev | 0.18198 | 0.182011 | 0.0168217
243/235 | Dev | 0.194712 | 0.194754 | 0.0213671
223/235 | Dev | 0.209221 | 0.209279 | 0.0274827
203/235 | Dev | 0.225876 | 0.225957 | 0.0358309
183/235 | Dev | 0.245146 | 0.245262 | 0.0474032
163/235 | Dev | 0.267635 | 0.267806 | 0.0637063
143/235 | Dev | 0.294126 | 0.294382 | 0.0870647
123/235 | Dev | 0.325643 | 0.326038 | 0.121115
103/235 | Dev | 0.363539 | 0.364163 | 0.171618
83/235 | Dev | 0.40961 | 0.410625 | 0.247788
63/235 | Dev | 0.466261 | 0.46796 | 0.364457
43/235 | Dev | 0.536719 | 0.539647 | 0.545532
23/235 | Dev | 0.62533 | 0.630516 | 0.829419
3/235 | Dev | 0.737962 | 0.747389 | 1.27736
Table 8: Exact and the winning DeVylder approximate values of $\Phi_{q}$, for $f(x)=\frac{12}{83}e^{-x}+\frac{42}{83}e^{-2x}+\frac{150}{83}e^{-3x}$, $\theta$ varies. $\theta$ | Closest approximation | Barrier exact | Barrier approx | % error Barrier
---|---|---|---|---
263/235 | Dev | 1.89732 | 1.91233 | 0.791183
243/235 | Dev | 1.79954 | 1.78002 | 1.08482
183/235 | Ren | 1.45224 | 1.52484 | 4.9989
163/235 | Ren | 1.31579 | 1.33691 | 1.60463
143/235 | Ren | 1.16804 | 1.12368 | 3.79796
123/235 | Expo | 1.00898 | 1.04123 | 3.19653
103/235 | Expo | 0.839228 | 0.794964 | 5.27444
83/235 | Expo | 0.660338 | 0.513179 | 22.2854
63/235 | Expo | 0.474896 | 0.196234 | 58.6785
43/235 | Expo | 0.286563 | 0 | 100
23/235 | Expo | 0.0998863 | 0 | 100
3/235 | All | 0 | 0 | 0
Table 9: Exact and approximate values of $b_{DeF}$ for
$f(x)=\frac{12}{83}e^{-x}+\frac{42}{83}e^{-2x}+\frac{150}{83}e^{-3x}$, when
$\theta$ varies. Unsurprisingly, when $b_{DeF}$ is near $0$, the natural
exponential approximation wins, but has poor performance, and as $b_{DeF}$
increases, Renyi and subsequently the de Vylder approximation take the upper
hand. For the smallest three values of $\theta$, all of the approximations
yielded $\widetilde{b}_{DeF}=0$ as the optimal barrier, with this being true
only for $\theta=3/235$.

We move now to the dividend problem with capital injections with cost $k\geq
1$ as in Theorem 1. One can compute the value function $J_{0}$ at $x=0$ in
terms of $W$, $Z$,$C$, $S$, and $G$ – see equation 5.

To provide a more concrete example, fixing $q=\frac{5}{48}$,$P=0$, $k=3/2$ as
input parameters we compute for values of $J_{0}$ as a function of $\theta$,
with results summarized in the tables below. The tables provide comparisons of
the computed optimal quantities $J_{0}$, $a$, and $b$ to an approximation
using all exponential inputs (referred to as $J_{0}$, $a$, and $b$ expo pure)
and to an approximation which uses actual inputs but computed using the
exponential formula as described in equation 3 (referred to as $J_{0}$, $a$,
and $b$ expo CI).

| J0
---|---
$\theta$ | $J_{0}$ |  $J_{0}$ expo pure |  $J_{0}$ expo pure error |  $J_{0}$ expo CI |  $J_{0}$ expo CI error
263/235 | 3.7747 | 3.76883 | 0.155556 | 4.11784 | 9.09041
243/235 | 3.41491 | 3.38603 | 0.845606 | 3.74156 | 9.5654
223/235 | 3.0636 | 3.00802 | 1.81444 | 3.36985 | 9.99637
203/235 | 2.72335 | 2.63828 | 3.1238 | 3.00466 | 10.3296
183/235 | 2.39737 | 2.28225 | 4.80185 | 2.64879 | 10.4871
163/235 | 2.08958 | 1.94765 | 6.79226 | 2.3062 | 10.3665
143/235 | 1.80446 | 1.64396 | 8.8946 | 1.9823 | 9.85516
123/235 | 1.54668 | 1.38072 | 10.73 | 1.68379 | 8.86472
103/235 | 1.32041 | 1.16526 | 11.7499 | 1.4178 | 7.37587
83/235 | 1.12864 | 1.00194 | 11.2258 | 1.19022 | 5.45555
63/235 | 0.972835 | 0.88785 | 8.73585 | 1.00404 | 3.20798
43/235 | 0.852739 | 0.789923 | 7.36635 | 0.859039 | 0.738837
23/235 | 0.751597 | 0.701299 | 6.69218 | 0.751597 | 0
3/235 | 0.660372 | 0.620567 | 6.02761 | 0.660372 | 0
Table 10: Values of $J_{0}$ compared with approximations using all exponential inputs ($J_{0}$ expo pure) and actual inputs but computed using the exponential formula ($J_{0}$ expo CI). The pure exponential approximation does a good job of approximating $J_{0}$ for higher values of $\theta$ considered, while the exponential CI approximation seemed to fair better for lower $\theta$ values | a
---|---
$\theta$ | $a$ |  $a$ expo pure |  $a$ expo pure error |  $a$ expo CI |  $a$ expo CI error
263/235 | 2.51647 | 2.74523 | 0.155536 | 2.51255 | 9.09042
243/235 | 2.27661 | 2.49437 | 0.845596 | 2.25735 | 9.56541
223/235 | 2.0424 | 2.24657 | 1.81443 | 2.00535 | 9.99638
203/235 | 1.81557 | 2.00311 | 3.1238 | 1.75885 | 10.3296
183/235 | 1.59825 | 1.76586 | 4.80185 | 1.5215 | 10.4871
163/235 | 1.39306 | 1.53747 | 6.79226 | 1.29844 | 10.3665
143/235 | 1.20298 | 1.32153 | 8.8946 | 1.09598 | 9.85516
123/235 | 1.03112 | 1.12252 | 10.73 | 0.920479 | 8.86472
103/235 | 0.880271 | 0.945198 | 11.7499 | 0.77684 | 7.37587
83/235 | 0.752428 | 0.793477 | 11.2258 | 0.667962 | 5.45555
63/235 | 0.648557 | 0.669362 | 8.73585 | 0.5919 | 3.20798
43/235 | 0.568493 | 0.572693 | 7.36635 | 0.526616 | 0.738838
23/235 | 0.501065 | 0.501065 | 6.69218 | 0.467532 | 0
3/235 | 0.440248 | 0.440248 | 6.02761 | 0.413711 | 0
Table 11: Values of $a$ compared with approximations using all exponential inputs ($a$ expo pure) and actual inputs but computed using the exponential formula ($a$ expo CI) | b
---|---
$\theta$ | $b$ |  $b$ expo pure |  $b$ expo pure error |  $b$ expo CI |  $b$ expo CI error
263/235 | 0.709355 | 0.805116 | 13.4997 | 0.677918 | 4.43179
243/235 | 0.695874 | 0.801936 | 15.2416 | 0.671779 | 3.46259
223/235 | 0.677601 | 0.794377 | 17.2337 | 0.662801 | 2.18425
203/235 | 0.653005 | 0.779265 | 19.3352 | 0.649805 | 0.490147
183/235 | 0.620126 | 0.751601 | 21.2012 | 0.631097 | 1.76912
163/235 | 0.576553 | 0.704104 | 22.123 | 0.604293 | 4.81128
143/235 | 0.519526 | 0.627369 | 20.7579 | 0.566198 | 8.98366
123/235 | 0.446259 | 0.511076 | 14.5246 | 0.512961 | 14.947
103/235 | 0.354524 | 0.346046 | 2.39143 | 0.440755 | 24.3231
83/235 | 0.243362 | 0.126054 | 48.2032 | 0.347059 | 42.6099
63/235 | 0.113593 | 0 | 100 | 0.231975 | 104.216
43/235 | 0 | 0 | 0 | 0.0987484 | 0
23/235 | 0 | 0 | 0 | 0 | 0
3/235 | 0 | 0 | 0 | 0 | 0
Table 12: Values of $b$ compared with approximations using all exponential
inputs ($b$ expo pure) and actual inputs but computed using the exponential
formula ($b$ expo CI)

To provide a point of comparison, we fix $q=\frac{5}{48}$, and compute the de
Finetti barrier to be $b_{DeF}=1.89732$ and the corresponding dividend value
function when starting at $x=0$ to be $J_{DeF}=1.99847$.

$k$ |  $J_{0}$ % deviation | $J_{0}-J_{DeF}$ | $J_{0}$ |  $b$ % deviation | $b$
---|---|---|---|---|---
1 | 154.123 | 3.0801 | 5.07857 | 100 | 0
2 | 65.714 | 1.31327 | 3.31174 | 43.0208 | 1.08108
3 | 42.863 | 0.856604 | 2.85507 | 25.4792 | 1.4139
4 | 31.4465 | 0.628448 | 2.62692 | 17.6853 | 1.56178
5 | 24.6995 | 0.493611 | 2.49208 | 13.4234 | 1.64264
6 | 20.2855 | 0.4054 | 2.40387 | 10.7759 | 1.69287
7 | 17.1884 | 0.343504 | 2.34197 | 8.98441 | 1.72686
8 | 14.9014 | 0.2978 | 2.29627 | 7.69619 | 1.7513
9 | 13.1463 | 0.262724 | 2.26119 | 6.72732 | 1.76968
10 | 11.758 | 0.23498 | 2.23345 | 5.97306 | 1.78399
100 | 1.11095 | 0.0222019 | 2.02067 | 0.533448 | 1.8872
1000 | 0.110409 | 0.00220648 | 2.00067 | 0.0527257 | 1.89632
10000 | 0.0109536 | 0.000218903 | 1.99869 | 0.00533526 | 1.89722
Table 13: Values of $J_{0}$ and $b$ in presence of capital injections compared
to the case where capital injections are non-existent, $J_{DeF}=1.99847$ and
$b_{DeF}=1.89732$. As $k$ is increased one can see that $J_{0}$ and $b$
approaches $J_{DeF}$ and $b_{DeF}$. This is expected since higher costs of
injecting capital makes it less viable, hence it is treated like the concept
does not exist.

###  5.3 A Cramér-Lundberg process with oscillating density and scale function

In the following example, we study a Cramèr-Lundberg model with density of
claims given by

$\displaystyle f(x)=$ $\displaystyle ue^{-ax}2\cos^{2}\left(\frac{\omega
x+\phi}{2}\right)=ue^{-ax}(1+\cos(\omega x+\phi))=$ $\displaystyle=$
$\displaystyle e^{-ax}(u+u\cos(\phi)\cos(\omega x)-u\sin(\phi)\sin(\omega x))$

where

$u=\frac{a\left(a^{2}+\omega^{2}\right)}{a^{2}+\omega^{2}+a^{2}\cos(\phi)-a\omega\sin(\phi)}.$

Assuming further that $a=1$, $\phi=2$, $\omega=20$, and that $\theta=1$,
$q=1/10$, the Laplace exponent for this process is
$\kappa(s)=\frac{s\left(2.09898s^{3}+5.29695s^{2}+843.502s+420.846\right)}{(s+1.)\left(s^{2}+2.s+401.\right)}$
and the scale function is

$\displaystyle W_{q}(x)$
$\displaystyle=0.824723e^{0.0881484x}-0.348141e^{-0.540677x}$
$\displaystyle+e^{-1.0117x}\cos(19.9957x)~{}\Big{(}-(0.000285494\,+0.0000804151i)\sin(39.9914x)$
$\displaystyle-(0.0000804151\,+0.000285494i)+(-0.0000804151+0.000285494i)\cos(39.9914x)\Big{)}$
$\displaystyle+e^{-1.0117x}\sin(19.9957x)~{}\Big{(}-(0.0000804151\,-0.000285494i)\sin(39.9914x)$
$\displaystyle+(0.000285494\,+0.0000804151i)\cos(39.9914x)-(0.000285494\,-0.0000804151i)\Big{)}.$

(a) $W^{\prime}_{q}(x)$

(b) $W^{\prime\prime}_{q}(x)$

Figure 9: Plots of $W^{\prime}_{q}(x)$, and $W^{\prime\prime}_{q}(x)$ of the exact solution and the approximations for $f(x)=ue^{-ax}2\cos^{2}\left(\frac{\omega x+\phi}{2}\right)$, $\theta=1$, $q=\frac{1}{10}$. |  | Dominant exponent
---
$\Phi_{q}$
| Percent relative error
---
($\Phi_{q}$)
| Optimal barrier
---
$b_{DeF}$
| Percent relative error
---
($b_{DeF}$)
Exact | 0.0881484 | 0 | 4.38201 | 0
Expo | 0.0878658 | 0.32053 | 4.42263 | 0.927122
Renyi | 0.0881481 | 0.000314617 | 4.39788 | 0.362284
Dev | 0.0881484 | 6.11743*10^-6 | 4.39745 | 0.352331
Table 14: Exact and approximate values of $\Phi_{q}$ and $b_{DeF}$ for
$f(x)=ue^{-ax}2\cos^{2}\left(\frac{\omega x+\phi}{2}\right)$, $\theta=1$,
$q=\frac{1}{10}$. The DeVylder approximation wins on both fronts.

Clearly, our completely monotone approximation cannot fully reproduce
functions like $W_{q}^{\prime}(x),W_{q}^{\prime\prime}(x)$ in examples like
this where oscillations occur (note however that the de Finetti optimal
barrier is well approximated here). If a more exact reproduction is necessary,
higher order approximations should be used.

##  6 The maximal error of exponential approximations $J_{0}$ along one
parameter families of Cramér-Lundberg processes

In this section, we provide the two approximations for the dividend value with
capital injections $J_{0}$, and the dividend barrier $b$, for two one
parameter families of Cramér-Lundberg processes, with densities given
respectively by:

$\displaystyle f(x)=k_{\epsilon}\left[e^{-x}+\epsilon e^{-2x}\right]$ (41)
$\displaystyle
f(x)=k_{\epsilon}\left[\frac{12}{83}e^{-x}+\epsilon\left(\frac{42}{83}e^{-2x}+\frac{150}{83}e^{-3x}\right)\right],$
(42)

where $k_{\epsilon}$ is the normalization constant, and compute the maximal
error of approximation when $\epsilon\in(0,\infty)$ and $\theta\approx 1$. For
this choice, the pure exponential approximation works considerably better,
$\forall\epsilon$.

| J0
---|---
| J0 exact | J0 expo pure | J0 expo pure error | J0 expo CI | J0 expo CI error
0.001 | 7.1879 | 7.18802 | 0.0016603 | 7.18849 | 0.00827967
0.01 | 7.17075 | 7.17193 | 0.0164663 | 7.17666 | 0.0824782
0.1 | 7.008 | 7.01863 | 0.151653 | 7.06358 | 0.793041
1 | 5.95034 | 5.99151 | 0.691856 | 6.26009 | 5.20551
10 | 4.20175 | 4.19406 | 0.183122 | 4.40089 | 4.73941
100 | 3.66909 | 3.6654 | 0.100631 | 3.69555 | 0.721228
1000 | 3.6025 | 3.60208 | 0.0117065 | 3.60523 | 0.075585
Table 15: $\lambda=1$, $\theta=1$, $q=\frac{1}{10}$ $k=3/2$ and $P=0$. As
expected, the errors decrease both as $\epsilon$ goes to zero and infinity
since the densities approach an exponential density.

Figure 10: $J_{0}$ values and errors plotted against $\epsilon$. Errors peak
at $\epsilon=1$.

We do the same thing for the family of densities given by
$f(x)=k_{\epsilon}\left[\frac{12}{83}e^{-x}+\epsilon\left(\frac{42}{83}e^{-2x}+\frac{150}{83}e^{-3x}\right)\right]$.

| J0
---|---
$\epsilon$ | J0 exact | J0 expo pure | J0 expo pure error | J0 expo CI | J0 expo CI error
0.001 | 7.95508 | 7.95771 | 0.0330111 | 7.96565 | 0.132796
0.01 | 7.68765 | 7.71127 | 0.307292 | 7.78772 | 1.30166
0.1 | 6.06176 | 6.15381 | 1.5186 | 6.62641 | 9.31507
1 | 3.7747 | 3.76883 | 0.155556 | 4.11784 | 9.09041
10 | 3.1382 | 3.1379 | 0.00959692 | 3.23354 | 3.03813
100 | 3.05894 | 3.06427 | 0.174306 | 3.12284 | 2.08921
1000 | 3.0508 | 3.05678 | 0.196219 | 3.11149 | 1.98956
Table 16: $\lambda=1$, $c=1$, $q=\frac{5}{48}$ $k=3/2$ and $P=0$. As
$\epsilon$ goes to zero, the density becomes exponential hence the decrease in
errors. As $\epsilon$ goes to infinity, the density approaches a hyper
exponential density of order 2, but still both methods of approximating
$J_{0}$ yield reasonable results.

Figure 11: $J_{0}$ values and errors plotted against $\epsilon$. Errors peak
at $\epsilon=0.1$.

##  7 The profit function when the claims are distributed according to a
matrix exponential jumps density

Consider now the more general case when the claims are distributed according
to a matrix exponential density generated by a row vector $\vec{\beta}$ and by
an invertible matrix $B$ of order $n$, which are such that the vector
$\vec{\beta}e^{xB}$ is decreasing componentwise to $0$, and
$\vec{\beta}.{\mbox{\boldmath$1$}}\neq 0$, with $1$ a column vector. As
customary, we restrict w.l.o.g. to the case when $\vec{\beta}$ is a
probability vector, and $\vec{\beta}.{\mbox{\boldmath$1$}}=1$, so that

$\overline{F}(x)=\vec{\beta}e^{xB}{\mbox{\boldmath$1$}}$

is a valid survival function.

The matrix versions of our functions are:

$\begin{cases}C_{a}(x)=\lambda\int_{0}^{x}W_{q}(x-y)\;\overline{F}(y+a)\;dy=\lambda\vec{\beta}\int_{0}^{x}W_{q}(x-y)\;e^{yB}\;dy\;e^{aB}{\mbox{\boldmath$1$}}=\vec{C}(x)e^{aB}{\mbox{\boldmath$1$}}\\\
m_{a}(y)=\int_{0}^{a}zf(y+z)dz=\vec{\beta}\;e^{yB}\;\int_{0}^{a}z\;e^{zB}(-B)\;dz\;{\mbox{\boldmath$1$}}=\vec{\beta}\;e^{yB}M(a){\mbox{\boldmath$1$}}\\\
G_{a}(x)=\lambda\int_{0}^{x}W_{q}(x-y)\;m_{a}(y)\;\;dy=\vec{C}(x)M(a){\mbox{\boldmath$1$}}\;\;\;\;\end{cases},$
(43)

where

$\begin{cases}C(x)=\lambda\int_{0}^{x}W_{q}(x-y)\;e^{yB}\;dy\\\
\vec{C}(x)=\lambda\vec{\beta}\int_{0}^{x}W_{q}(x-y)\;e^{yB}\;dy\end{cases}.$
(44)

The product formulas (43) may also be established directly in the phase-type
case, using the conditional independence of the ruin probability of the
overshoot size.

We derive first these extensions from scratch for $(\vec{\beta},B)$ phase-type
densities, in order to highlight their probabilistic interpretation. Later, we
will show that the matrix exponential jumps case follows as a particular case
of [GK17].

Recall first [AA10] that $\Psi_{q}(x)=\vec{\Psi}_{q}(x){\mbox{\boldmath$1$}},$
where $\vec{\Psi}_{q}(x)$ is a vector whose components represent the
probability that ruin occurs during a certain phase, and that the conditional
independence of ruin and overshoots translates into the product formula

$\Psi_{q}(x,y):=P_{x}[T_{0,-}<\infty,X_{T_{0,-}}<-y]=\vec{\Psi}_{q}(x)e^{yB}{\mbox{\boldmath$1$}}.$
(45)

To take advantage of this, it is convenient to replace from the beginning
$Z_{q}(x)$ by $\Psi_{q}(x)$, taking advantage of the formula [AKP04, Kyp14]

$Z_{q}(x)=\Psi_{q}(x)+W_{q}(x)\frac{q}{\Phi_{q}}\Longrightarrow
C(x)=(c-\frac{q}{\Phi_{q}})W_{q}(x)-\Psi_{q}(x).$ (46)

Alternatively, one may introduce a vector function

$\vec{Z}_{q}(x):=\vec{\Psi}_{q}(x)+W_{q}(x)\frac{q}{\Phi_{q}}{\mbox{\boldmath$1$}}.$
(47)

On the other hand, the mean function may be written as

$m_{a}=\int_{0}^{a}y\;F(\mathrm{d}y)=-a\overline{F}(a)+\int_{0}^{a}\overline{F}(x)dx=\vec{\beta}M(a){\mbox{\boldmath$1$}},M(a)=-B^{-1}-e^{aB}\left(aI_{n}-B^{-1}\right).$

The following result follows in the phase-type case just as in the exponential
case [AGLW20]; in the matrix exponential jumps case, it may be obtained from
[GK17]:

######  PROPOSITION 6.

For a Cramèr-Lundberg process (compound Poisson ) with matrix exponential
jumps of type $(\vec{\beta},B)$, it holds that

  1. 1.

$J_{x}=\begin{cases}kG_{a}(x)+J_{0}S_{a}(x)=kG_{a}(x)+\frac{1-kG^{\prime}(b)}{S^{\prime}(b)}S_{a}(x),&x\in[0,b]\\\
kx+J_{0}&x\in[-a,0]\\\ 0&x\leq-a\end{cases},$ (48)

where

$\begin{cases}C(x)=\lambda\int_{0}^{x}W_{q}(x-y)\;e^{yB}\;dy\\\
\vec{C}(x)=\lambda\vec{\beta}\int_{0}^{x}W_{q}(x-y)\;e^{yB}\;dy=cW_{q}(x)\vec{1}-\vec{Z}_{q}(x)=(c-\frac{q}{\Phi_{q}})W_{q}(x)\vec{1}-\vec{\Psi}_{q}(x)\\\
G_{a}(x)=\vec{C}_{q}(x)M(a){\mbox{\boldmath$1$}}\\\
R_{a}(x)=S_{a}(x)-Z_{q}(x)=\vec{C}_{q}(x)e^{aB}{\mbox{\boldmath$1$}}\end{cases},$
(49)

and

$J_{0}=\frac{1-k\vec{C}_{q}^{\prime}(b)\;M(a){\mbox{\boldmath$1$}}}{qW_{q}(b)+\vec{C}_{q}^{\prime}(b)e^{aB}{\mbox{\boldmath$1$}}}.$
(50)

  2. 2.

For fixed $a$, the optimality equation $\frac{\partial}{\partial
b}J_{0}^{a,b}=0$ simplifies to

$J_{0}=\frac{k\vec{C}_{q}^{\prime\prime}(b)\;M(a){\mbox{\boldmath$1$}}}{qW_{q}^{\prime}(b)+\vec{C}_{q}^{\prime\prime}(b)e^{aB}{\mbox{\boldmath$1$}}}.$
(51)

######  REMARK 12.

The additive separation of $a,b$ which was the basis of proving optimality in
the exponential case does not seem possible anymore, but (50) allows the
numeric computation of the optimum.

## References

  * [AA10] Hansjörg Albrecher and Sören Asmussen, _Ruin probabilities_ , vol. 14, World Scientific, 2010.
  * [AAM20] Hansjoerg Albrecher, Pablo Azcue, and Nora Muler, _Optimal ratcheting of dividends in insurance_ , SIAM Journal on Control and Optimization 58 (2020), no. 4, 1822–1845.
  * [ABH18] Florin Avram, Abhijit Datta Banik, and András Horváth, _Ruin probabilities by Padé’s method: simple moments based mixed exponential approximations (Renyi, De Vylder, Cramér–Lundberg), and high precision approximations with both light and heavy tails_ , European Actuarial Journal (2018), 1–27.
  * [ACFH11] Florin Avram, Donatien Chedom-Fotso, and András Horváth, _On moments based Padé approximations of ruin probabilities_ , Journal of Computational and Applied mathematics 235 (2011), no. 10, 3215–3228.
  * [AGLW20] Florin Avram, D. Goreac, J Li, and X Wu, _Equity cost induced dichotomy for optimal dividends in the cramér-lundberg model_ , IME (2020).
  * [AGVA19] Florin Avram, D. Grahovac, and C. Vardar-Acar, _The $W,Z$ scale functions kit for first passage problems of spectrally negative Lévy processes, and applications to the optimization of dividends_, arXiv preprint (2019).
  * [AHPS19] Florin Avram, Andras Horváth, Serge Provost, and Ulyses Solon, _On the padé and laguerre–tricomi–weeks moments based approximations of the scale function w and of the optimal dividends barrier for spectrally negative lévy risk processes_ , Risks 7 (2019), no. 4, 121.
  * [AKP04] Florin Avram, Andreas Kyprianou, and Martijn Pistorius, _Exit problems for spectrally negative Lévy processes and applications to (Canadized) Russian options_ , The Annals of Applied Probability 14 (2004), no. 1, 215–238.
  * [AM14] Pablo Azcue and Nora Muler, _Stochastic optimization in insurance: A dynamic programming approach_ , Springer, 2014.
  * [AP14] Florin Avram and Martijn Pistorius, _On matrix exponential approximations of ruin probabilities for the classic and brownian perturbed cramér-lundberg processes_ , Insurance: Mathematics and Economics 59 (2014), 57–64.
  * [APP07] Florin Avram, Zbigniew Palmowski, and Martijn R Pistorius, _On the optimal dividend problem for a spectrally negative Lévy process_ , The Annals of Applied Probability (2007), 156–180.
  * [APY18] Florin Avram, José Luis Pérez, and Kazutoshi Yamazaki, _Spectrally negative Lévy processes with Parisian reflection below and classical reflection above_ , Stochastic processes and Applications 128 (2018), 255–290.
  * [BDW07] Christopher J Beveridge, David CM Dickson, and Xueyuan Wu, _Optimal dividends under reinsurance_ , Centre for Actuarial Studies, Department of Economics, University of …, 2007.
  * [Ber98] Jean Bertoin, _Lévy processes_ , vol. 121, Cambridge university press, 1998.
  * [BFS08] PB Brito, F Fabiao, and A Staubyn, _Euler, lambert, and the lambert w-function today._ , Mathematical Scientist 33 (2008), no. 2.
  * [BHT05] Andrea Bobbio, András Horváth, and Miklós Telek, _Matching three moments with minimal acyclic phase type distributions_ , Stochastic models 21 (2005), no. 2-3, 303–326.
  * [Boy98] JP Boyd, _Global approximations to the principal real-valued branch of the lambert w-function_ , Applied mathematics letters 11 (1998), no. 6, 27–31.
  * [CGH+96] Robert M Corless, Gaston H Gonnet, David EG Hare, David J Jeffrey, and Donald E Knuth, _On the lambertw function_ , Advances in Computational mathematics 5 (1996), no. 1, 329–359.
  * [DD05] David Dickson and Steve Drekic, _Optimal dividends under a ruin probabilityconstraint_ , 2005.
  * [DF57] Bruno De Finetti, _Su un’impostazione alternativa della teoria collettiva del rischio_ , Transactions of the XVth international congress of Actuaries, vol. 2, 1957, pp. 433–443.
  * [DcSA16] Bogdan Dumitrescu, Bogdan C cSicleru, and Florin Avram, _Modeling probability densities with sums of exponentials via polynomial approximation_ , Journal of Computational and Applied Mathematics 292 (2016), 513–525.
  * [DV78] Fl De Vylder, _A practical solution to the problem of ultimate ruin probability_ , Scandinavian Actuarial Journal 1978 (1978), no. 2, 114–119.
  * [ES11] Julia Eisenberg and Hanspeter Schmidli, _Minimising expected discounted capital injections by reinsurance in a classical risk model_ , Scandinavian Actuarial Journal 2011 (2011), no. 3, 155–176.
  * [Ger69] Hans U Gerber, _Entscheidungskriterien für den zusammengesetzten poisson-prozess_ , Ph.D. thesis, ETH Zurich, 1969.
  * [GK17] Leslaw Gajek and Lukasz Kuciński, _Complete discounted cash flow valuation_ , Insurance: Mathematics and Economics 73 (2017), 1–19.
  * [GSS08] Hans U Gerber, Elias SW Shiu, and Nathaniel Smith, _Methods for estimating the optimal dividend barrier and the probability of ruin_ , Insurance: Mathematics and Economics 42 (2008), no. 1, 243–254.
  * [Høj02] Bjarne Højgaard, _Optimal dynamic premium control in non-life insurance. maximizing dividend pay-outs_ , Scandinavian Actuarial Journal 2002 (2002), no. 4, 225–245.
  * [KKR13] Alexey Kuznetsov, Andreas E Kyprianou, and Victor Rivero, _The theory of scale functions for spectrally negative Lévy processes_ , Lévy Matters II, Springer, 2013, pp. 97–186.
  * [KS08] Natalie Kulenko and Hanspeter Schmidli, _Optimal dividend strategies in a cramér–lundberg model with capital injections_ , Insurance: Mathematics and Economics 43 (2008), no. 2, 270–278.
  * [Kyp14] Andreas Kyprianou, _Fluctuations of Lévy processes with applications: Introductory lectures_ , Springer Science & Business Media, 2014\.
  * [LL19] Kristoffer Lindensjö and Filip Lindskog, _Optimal dividends and capital injection under dividend restrictions_ , arXiv preprint arXiv:1902.06294 (2019).
  * [LZ08] Arne Løkka and Mihail Zervos, _Optimal dividend and issuance of equity policies in the presence of proportional costs_ , Insurance: Mathematics and Economics 42 (2008), no. 3, 954–961.
  * [NPY20] Kei Noba, José-Luis Pérez, and Xiang Yu, _On the bailout dividend problem for spectrally negative markov additive models_ , SIAM Journal on Control and Optimization 58 (2020), no. 2, 1049–1076.
  * [Pak15] Anthony G Pakes, _Lambert’s w meets kermack–mckendrick epidemics_ , IMA Journal of Applied Mathematics 80 (2015), no. 5, 1368–1386.
  * [Pis04] Martijn R Pistorius, _On exit and ergodicity of the spectrally one-sided Lévy process reflected at its infimum_ , Journal of Theoretical Probability 17 (2004), no. 1, 183–220.
  * [PYB18] José-Luis Pérez, Kazutoshi Yamazaki, and Alain Bensoussan, _Optimal periodic replenishment policies for spectrally positive Lévy demand processes_ , arXiv preprint arXiv:1806.09216 (2018).
  * [Ram92] Colin M Ramsay, _A practical algorithm for approximating the probability of ruin_ , Transactions of the Society of Actuaries 44 (1992), 443–461.
  * [SLG84] Steven E Shreve, John P Lehoczky, and Donald P Gaver, _Optimal consumption for general diffusions with absorbing and reflecting barriers_ , SIAM Journal on Control and Optimization 22 (1984), no. 1, 55–75.
  * [Sup76] VN Suprun, _Problem of destruction and resolvent of a terminating process with independent increments_ , Ukrainian Mathematical Journal 28 (1976), no. 1, 39–51.
  * [VAVZ14] Eleni Vatamidou, Ivo Jean Baptiste Franccois Adan, Maria Vlasiou, and Bert Zwart, _On the accuracy of phase-type approximations of heavy-tailed risk models_ , Scandinavian Actuarial Journal 2014 (2014), no. 6, 510–534.
  * [VLSHGG+19] H Vazquez-Leal, MA Sandoval-Hernandez, JL Garcia-Gervacio, AL Herrera-May, and UA Filobello-Nino, _Psem approximations for both branches of lambert function with applications_ , Discrete Dynamics in Nature and Society 2019 (2019).