# Optimal Dividend Bands Revisited: A gradient-based method and Evolutionary
Algorithms

Hansjörg Albrecher  Department of Actuarial Science, Faculty of Business and
Economics, University of Lausanne and Swiss Finance Institute, Bâtiment
Extranef, Quartier UNIL-Dorigny, 1015 Lausanne, Switzerland
<EMAIL_ADDRESS>and  Brandon Garcia Flores  Department of
Actuarial Science, Faculty of Business and Economics, University of Lausanne,
Bâtiment Extranef, Quartier UNIL-Dorigny, 1015 Lausanne, Switzerland
<EMAIL_ADDRESS>
###### Abstract.

We reconsider the study of optimal dividend strategies in the Cramér-Lundberg
risk model. It is well-known that the solution of the classical dividend
problem is in general a band strategy. However, the numerical techniques for
the identification of the optimal bands available in the literature are very
hard to implement and explicit numerical results are known for very few cases
only. In this paper we put a gradient-based method into place which allows to
determine optimal bands in more general situations. In addition, we adapt an
evolutionary algorithm to this dividend problem, which is not as fast, but
applicable in considerable generality, and can serve for providing a
competitive benchmark. We illustrate the proposed methods in concrete
examples, reproducing earlier results in the literature as well as
establishing new ones for claim size distributions that could not be studied
before.

###### Key words and phrases:

Optimal dividend strategies; stochastic control; genetic algorithms; Cramér-
Lundberg risk model

##  1\. Introduction

Consider the optimal dividend problem for an insurance company whose surplus
process evolves according to the Cramér-Lundberg model (see e.g. [2]). The
company pays dividends to shareholders in continuous time according to some
admissible strategy $\pi$, and the objective is to identify the strategy that
maximizes the expected sum of discounted dividend payments until the event of
ruin. If there are no constraints on the size of the payments, Gerber [11]
showed that such a strategy always exists and is given by a band strategy that
partitions the interval $[0,\infty)$ into three sets
$\mathscr{B}_{0},\mathscr{B}_{c}$ and $\mathscr{B}_{\infty}$: whenever the
current surplus level is in $\mathscr{B}_{0}$, no dividends are paid, when the
current surplus level is in $\mathscr{B}_{c}$, all incoming premium is paid as
dividends so that the same surplus level is maintained until the next claim
arrives, and when the current surplus level is in $\mathscr{B}_{\infty}$, a
lump sum payment to the first surplus level outside of $\mathscr{B}_{\infty}$
is carried out. This leads to a cascading strategy towards ruin (see e.g. the
sample path illustration in Figure 1 later). Since Gerber’s result, the
optimal dividend problem and variants have been studied under many different
and more general model assumptions and under various constraints (see for
instance Avanzi [3] and Albrecher & Thonhauser [1] for surveys).

For many claim size distributions of practical relevance, the above band
strategy turns out to collapse to a barrier strategy (i.e.,
$|\mathscr{B}_{c}|=1$), see e.g. Gerber & Shiu [13], and Avram et al. [4] and
Loeffen [20] for sufficient conditions on the Lévy measure under which a
barrier strategy is optimal in the general framework of spectrally negative
Lévy processes. In such a case the determination of the respective optimal
barrier is rather straight-forward, if the scale function corresponding to the
underlying surplus process is explicitly available (see e.g. Hubalek &
Kyprianou [14]). In contrast, already for the classical Cramér-Lundberg model
it turns out to be surprisingly challenging to (even numerically) identify the
sets $\mathscr{B}_{0},\mathscr{B}_{c},\mathscr{B}_{\infty}$ for more involved
claim size distributions.

The value function of the optimal dividend problem has been identified as a
viscosity solution of the corresponding Hamilton-Jacobi-Bellman equation by
Azcue & Muler [6], which led to an iterative procedure for the numerical
determination of the optimal bands (cf. [6, 7] as well as Schmidli [25] and
Berdel [8]). However, when following the respective algorithm, one faces two
main difficulties:

  1. (1)

There is no complete understanding of the internal mechanism which furnishes
the sets $\mathscr{B}_{0},\mathscr{B}_{c}$ and $\mathscr{B}_{\infty}$ given an
arbitrary set of parameters of the Cramér-Lundberg model, so that, as of now,
concrete numerical solutions are known only for very few concrete and simple
claim distributions.

  2. (2)

The numerical approaches suggested in the literature require the solution of
several differential equations, which can be very expensive in terms of
computation time.

While the numerical algorithms suggested so far were a by-product of the
meticulous study of the existence and uniqueness of the solution of the
optimal dividend problem and its characterization, in this paper we would like
to take a different route. Relying on the fact that a band strategy is optimal
and assuming that there are only finitely many such bands, we are interested
to see if there are numerical alternatives to determine these optimal bands
more efficiently and/or generally. In particular, we propose two respective
numerical algorithms that differ from previous approaches.
The first one exploits the ’cascading’ nature of band strategies and
establishes a method based on gradients, when the value function is considered
as a function of all band levels rather than the initial surplus level. This
will lead to a rather fast numerical routine that can also be tailored to work
for cases in which the scale function is not explicitly available, but its
Laplace transform is (which is the case in general, since that Laplace
transform is defined as the reciprocal of the (shifted) Laplace exponent of
the underlying spectrally negative Lévy process). For every fixed number of
bands, one then obtains an iterative procedure and can finally compare whether
increasing the number of bands still improves the solution.
The second method also uses the explicit iterative formula for the expected
discounted dividend payments for given band levels and uses it as the
objective function in a general evolutionary algorithm. Evolutionary
strategies (ES) have been applied in the past with some success in reinsurance
problems where the evaluation of the function to optimize is only possible
through numerical procedures due to the non-existence of explicit algebraic
expressions (see, for example, Salcedo-Sanz et al. [23] and Roman et al.
[21]). In our context, we propose the use of an evolutionary self-adaptive
strategy based on the algorithm originally proposed in the survey paper by
Beyer & Schwefel [10] which does not use the derivatives of involved functions
and which can be easily implemented in common programming languages. This
genetic algorithm is rather flexible and works particularly well in high-
dimensional problems. We will hence adapt such an algorithm to the needs of
our dividend problem, and indeed get to the same solutions as the other
methods do. While in reasonably low dimension (like the dividend problem
typically is, as there are only a few bands to consider) the computation time
using this algorithm is not favorable when compared to the gradient-based
approach, it is applicable in very general setups and can nicely serve as a
benchmark for numerical procedures. Furthermore, it can also be useful in
other application areas in risk theory, and since the idea and implementation
of evolutionary algorithms may not be so commonly familiar in the risk theory
community, we present the underlying principle and implementation variants in
some detail in a separate section.

The rest of the paper is organized as follows. Section 2 contains some
definitions and preliminaries on the model assumptions and the nature of the
dividend band strategies. In Section 3 we summarize the previously available
numerical procedures for the determination of optimal dividend bands. Section
4 then provides the expressions for the expected discounted payments that will
be used in the numerical algorithms, in particular with respect to the
’cascading’ view. Section 5 develops the gradient-based algorithm and
discusses issues of its implementation. In Section 6 we give a general account
of the idea behind evolutionary algorithms, and the necessary adaptations to
the optimal dividend problem are discussed in Section 7. Finally, in Section 8
we provide numerical illustrations. We first re-derive the known optimal bands
for the well-known example of an Erlang(2) claim size distribution derived in
Azcue & Muler [6] as well as the mixture of Erlang claim size distribution
established in Berdel [8]. We then establish a new instance of a mixure of
Erlang distributions for which a 4-band strategy is optimal. Subsequently we
use our algorithms to derive the optimal barrier level in a risk model with
Pareto claim sizes (where a barrier is known to be optimal due to Loeffen
[20]). We also implement an example with a mixture of Erlang and Pareto
claims, which could not be handled with previously existing techniques and for
which also two bands turn out to be optimal. Section 9 concludes.

##  2\. Definitions and preliminaries

Consider the surplus process of an insurance portfolio in a Cramer-Lundberg
model

$C_{t}=u+pt-\sum_{k=1}^{N_{t}}Y_{k},\quad t\geq 0,$

with $(N_{t})_{t\geq 0}$ a homogeneous Poisson process with rate $\lambda$
representing the arrival of claims, $(Y_{k})_{k\geq 0}$ a sequence of
absolutely continuous i.i.d. claim size random variables with density $f_{Y}$
and finite mean $\mu$, and the premium rate $p$ satisfying the positive safety
loading condition $p=(1+\eta)\lambda\mu$ for some $\eta>0$. Let
$(\mathscr{F}_{t})_{t\geq 0}$ be the usual augmentation of the filtration
generated by $(C_{t})_{t\geq 0}$.
The dividend strategy $\pi$ is represented by the process $(U_{t})_{\geq 0}$,
where $U_{t}$ are the dividends paid up to time $t$. A dividend strategy is
called admissible if the associated process $(U_{t})_{t\geq 0}$ is adapted to
$(\mathscr{F}_{t})_{t\geq 0}$, non-decreasing and càglàd. Denote by $\Pi$ the
set of all admissible strategies.
For an admissible strategy $\pi\in\Pi$ we denote by $X_{t}=C_{t}-U_{t}$ the
surplus after dividend payments. Let

$\tau_{D}=\inf\\{t\geq 0\mid X_{t}<0\\}$

be the time until ruin, then

$V_{\pi}(u)=\mathbb{E}\left[\int_{0}^{\tau_{D}}e^{-\delta
t}dU_{t}{{\mid}}X_{0}=u\right]$

is the expected value of the aggregated discounted dividends paid until ruin,
where $\delta>0$ is the force of interest. The objective is then to identify
the strategy that maximizes $V$ over all admissible strategies, that is, to
find a strategy $\pi^{\ast}$ such that

(1) $V_{\pi^{\ast}}=\sup\limits_{\pi\in\Pi}V_{\pi}(x).$

As pointed out in the introduction, the class of band strategies is known to
be optimal in this case (cf. [11]). A band strategy is defined by a partition
of the positive half-line
${\mathbb{R}}=\mathscr{B}_{0}\cup\mathscr{B}_{c}\cup\mathscr{B}_{\infty}$ with
the following properties:

  * •

If $x\in\mathscr{B}_{0}$, there exists $\varepsilon>0$ such that
$[x,x+\varepsilon)\subset\mathscr{B}_{0}$.

  * •

$\mathscr{B}_{c}$ is compact.

  * •

$\mathscr{B}_{\infty}$ is open in $[0,\infty)$.

  * •

If $x\not\in\mathscr{B}_{\infty}$ and there is a sequence
$(x_{n})\subset\mathscr{B}_{\infty}$ such that $x_{n}\to x$, then
$x\in\mathscr{B}_{c}$.

  * •

$(\sup\mathscr{B}_{c},\infty)\subset\mathscr{B}_{\infty}$.

$\mathscr{B}_{0}$ corresponds to all surplus levels for which $dU_{t}=0$ (no
dividends are being paid), $\mathscr{B}_{c}$ is the set of surplus levels for
which $dU_{t}=p\,dt$ (all incoming premium is paid as dividends) and
$\mathscr{B}_{\infty}$ is the set of surplus levels at which
$U_{t+}-U_{t}=X_{t}-\sup\\{b\in\mathscr{B}_{c}\mid X_{t}>b\\}$ is applied (the
smallest possible lump sum is paid with which one leaves the set
$\mathscr{B}_{\infty}$). As the focus of this paper is to provide alternatives
to numerically identifying the optimal bands, we restrict ourselves to
finitely many bands, i.e., for given levels $a_{0}=0\leq b_{0}\leq
a_{1}\leq\cdots\leq b_{m-1}$, the band strategy is given by

$\mathscr{B}_{0}=[0,b_{0})\cup\bigcup_{k=1}^{m-1}[a_{k},b_{k}),\;\mathscr{B}_{c}=\bigcup_{k=1}^{m-1}\\{b_{k}\\},\;\mathscr{B}_{\infty}=\bigcup_{k=0}^{m-2}(b_{k},a_{k+1})\cup(b_{m-1},\infty).$

We refer to this strategy as an m-band strategy (see Figure 1 for an
illustration of a sample path with a 2-band strategy).

Figure 1. A sample path with a 2-band strategy

We conclude this section by making some remarks about notation: for a set
$A\subset\mathbb{R}^{n}$, a function $f:A\to\mathbb{R}$, and a limit point
$x\in A$, we denote by $f(x_{1},\ldots,x_{j}-,\ldots,x_{n})$ the limit
$\lim_{y\to x}f(y)$ through points $y$ for which $y_{j}<x_{j}$, given that
this limit exists and similarly for $f(x_{1},\ldots,x_{j}+,\ldots,x_{n})$. The
function $f$ is then continuous at $x$ if and only if
$f(x_{1},\ldots,x_{j}-,\ldots,x_{n})=f(x_{1},\ldots,x_{j}+,\ldots,x_{n})$ for
each $j=1,\ldots,n$, provided all right-hand and left-hand limits exist.
Finally, in order to avoid cumbersome notation, we denote the partial
derivative in the $i$-th variable by $D_{i}f$.

##  3\. The identification of bands in previous literature

As mentioned in the introduction, the explicit identification of optimal bands
has proved challenging even in the classical Cramér-Lundberg risk model. In
the following we briefly summarize the available approaches in the literature.
The typical approach is to derive the Hamilton-Jacobi-Bellman (HJB) equation

(2) $\max\\{1-V^{\prime}(x),\mathcal{L}(V)(x)\\}=0,$

related to the stochastic control problem (1), where

(3)
$\mathcal{L}(f)(x)=pf^{\prime}(x)-(\lambda+\delta)f(x)+\lambda\int_{0}^{x}f(x-y)f_{Y}(y)\;dy,\quad
x>0$

is the infinitesimal generator of the discounted surplus process (see e.g.
Azcue and Muler [6, 7]). Since $V$ is typically not sufficiently regular to
satisfy the needs of (2) as a classical solution, one needs to look for
viscosity solutions, and it can be shown that $V(x)$ is the unique viscosity
solution of (2) satisfying a growth condition and a particular initial
condition [6]. The numerical approach to find this solution is then an
iterative procedure. In particular, Schmidli [25] proposed an algorithm for
finding the levels of the optimal bands. This algorithm was formalized by
Berdel [8], who considered the problem for the general case of phase-type
claim distributions (cf. Algorithm 1).

1

2

Input : Scale function $W_{\delta}$ and infinitesimal generator $\mathcal{L}$.

Output : Levels $B^{*}=(b_{0}^{*},a_{1}^{*},\ldots,b_{M-1}^{*})$ of the best
band strategy.

3

4begin

5 $m:=0$;

6 $b_{0}:=\sup\left\\{x\geq 0\mid W_{\delta}^{\prime}(x)=\inf\limits_{y\geq
0}W_{\delta}^{\prime}(y)\right\\}$;

7 $V_{0}(x):=\begin{cases}W_{\delta}(x)/W_{\delta}^{\prime}(b_{0})&x\leq
b_{0}\\\ V_{0}(b_{0})+x-b_{0}&x>b_{0}\end{cases}$;

8 while  _$\mathcal{L}(V_{m})(x) >0$ for some $x>b_{m}$_ do

9 $\mathscr{G}_{m}:=\\{f^{a}_{m+1},a>b_{m}\mid f^{a}_{m+1}(x)=V_{m}(x),\;x\leq
a\text{ and}$ $\mathcal{L}(f)^{a}_{m+1}(x)=0,\;x>a\\}$;

10
$a_{m+1}:=\inf\left\\{a>b_{m}\mid\inf\limits_{x>a}f^{a\prime}_{m+1}(x)=1\right\\}$;

11 $b_{m+1}:=\sup\left\\{x>a_{m+1}\mid f^{a_{m+1}\prime}_{m+1}(x)=1\right\\}$;

12 $V_{m+1}(x):=\begin{cases}f^{a_{m+1}}_{m+1}(x)&x\leq b_{m+1}\\\
f^{a_{m+1}}_{m+1}(b_{m+1})+x-b_{m+1}&x>b_{m+1}\end{cases}$;

13 $m:=m+1$

14 end while

15

16 end

17

Algorithm 1 Schmidli’s algorithm

Here, $W_{\delta}$ is the scale function (see e.g. [2, Ch.XI]). When the
optimal strategy is in fact a finite band strategy, Algorithm 1 is guaranteed
to converge. However, depending on the particular distribution of the claims,
the procedure can be computationally complex, as can be seen from Lines 1 to
1. In each iteration of the algorithm, a family of functions $\mathscr{G}_{m}$
parametrized by the interval $(b_{m},\infty)$ is defined in such a way that
for each $a>b_{m}$, the function $f^{a}_{m+1}$ solves (2) on $x>a$ with
boundary condition $f^{a}_{m+1}(a)=V_{m}(a)$. Apart from some cases where this
can be done explicitly, (2) has in general to be solved numerically. While
this might not represent a problem for a couple of values of $a>b_{m}$, the
difficulty arises when we consider Lines 1 and 1, since they presuppose a full
knowledge of the solutions in the entire interval $(b_{m},\infty)$ in order to
compute the necessary extrema. As stated in [8], one can define
$\bar{a}_{m+1}=\inf\left\\{x>b_{m}\mid\mathcal{L}(V_{m})(x)>0\right\\}$ and
restrict the family to the interval $(b_{m},\bar{a}_{m+1})$. However, this
introduces another extremum to be computed and one still has to consider the
trade-off that arises at each step of the procedure when choosing a grid fine
enough to discretize this new interval.

An alternative iterative algorithm for finding the optimal bands is proposed
in Avram et al. [5], using stochastic sub- and super-solutions of (2) (their
approach is formulated for general spectrally negative Lévy processes and the
inclusion of fixed transaction costs with every dividend payment). Similarly
to Algorithm 1, optimal levels are found in a sort of “upwards” approach
finding higher band levels at each step of the procedure. However, instead of
solving integro-differential equations, each step consists of solving a
stochastic control problem expressed through Gerber-Shiu functions. The
advantage of this is that the problem is then reduced to finding the extrema
of a low-dimensional function at each iteration, and there is no need for a
full set of solutions of (2) as seen in, for example, Lines 1 and 1 of
Algorithm 1. The Algorithm 2 presented later in this paper sets out from a
top-down approach, and then also leads to a bottom-up procedure that is
formulated via discounted deficit densities explicitly (rather than general
Gerber-Shiu functions), so that it eventually can be interpreted as a
particular customization and implementation of the algorithm by Avram et al.,
see the details below.

##  4\. Properties of the Value Function

In this section we recollect some properties of $V_{\pi}$, which will form the
basis for the implementation of the numerical algorithms presented later.
For a fixed set of levels $a_{0}=0\leq b_{0}\leq a_{1}\leq\cdots\leq b_{m-1}$
of an $m$-band strategy $\pi$, we observe the following: for any $0\leq k\leq
m-1$ and initial capital $u$ in $[a_{k},a_{k+1})$ (here $a_{0}=$ and
$a_{m}=\infty$), the amount of dividends paid in a realization of the process
will be the same as in a process with a barrier strategy with initial capital
$u-a_{k}$ and barrier $b_{k}-a_{k}$, up until a claim makes the original
process go below $a_{k}$. Denoting by $V_{b}$ the value function of a barrier
strategy with barrier $b$, space-homogeneity and the Markov property imply
that

(4)
$V_{\pi}(u)=V_{b_{k}-a_{k}}(u-a_{k})+\mathbb{E}\left[V_{\pi}(a_{k}-D(u-a_{k},b_{k}-a_{k}))\right],\quad
a_{k}\leq u<a_{k+1},$

where $D(u-a_{k},b_{k}-a_{k})$ denotes the deficit at ruin of a process with
initial capital $u-a_{k}$, for which a barrier strategy with barrier
$b_{k}-a_{k}$ is applied. In many cases, the density of the deficit at ruin
can be computed by means of Gerber-Shiu functions (see [13, 17]) and the
dividends-penalty identity (see [19, 12]). Assume henceforth that
$D(u-a_{k},b_{k}-a_{k})$ has a density, which we denote by
$f_{D}(\cdot,u-a_{k},b_{k}-a_{k})$. We observe that the variable inside the
expectation in (4) is non-zero only when the deficit is at most $a_{k}$. We
can therefore rewrite (4), for any $k=0,\ldots,m-1$, as

(5)
$V_{\pi}(u)=V_{b_{k}-a_{k}}(u-a_{k})+\int_{0}^{a_{k}}V_{\pi}(a_{k}-y)f_{D}(y,u-a_{k},b_{k}-a_{k})dy,\quad
a_{k}\leq u<a_{k+1}.$

This set of equations provides the central formulas for computing the value of
$V_{\pi}$ given a fixed set of levels: For $0\leq u<a_{1}$, the value of
$V_{\pi}(u)$ is equal to $V_{b_{0}}(u)$, which is given in terms of the scale
function of the process. We can then plug in these values in the integral in
Equation (5) to obtain the value of $V_{\pi}(u)$ for $a_{1}\leq u<a_{2}$ and
repeat this procedure in an iterative way to obtain the value of $V_{\pi}(u)$
for every $u$. The problem of evaluating $V_{\pi}(u)$ is therefore reduced to
the computation of the scale function $W_{\delta}$ and the density $f_{D}$.
However, with knowledge of the scale function, the latter can be computed by
means of the formula

(6)
$\begin{split}f_{D}(y,u,b)=\lambda\int_{0}^{\infty}\left(\frac{W_{\delta}(u)W_{\delta}^{\prime}(b-z)}{W_{\delta}^{\prime}(b)}-W_{\delta}(u-z)\right)f_{Y}(y+z)\;dz.\end{split}$

see, e.g., [17, Ch.X]. The setting in this last reference is that of general
Lévy processes. A more basic approach is to consider first the density
$f_{D^{0}}(y,u)$ of the deficit at ruin with initial capital $u$ in the
absence of a dividend strategy. Let $\hat{f_{Y}}(s)$ denote the Laplace
transform of the claim size density $f_{Y}$. Following e.g. [2, Ch.XII], we
know that $f_{D^{0}}(y,u)$ can be obtained as the inverse Laplace transform of

(7)
$\begin{split}\int_{0}^{\infty}e^{-su}f_{D^{0}}(y,u)\,du=\frac{\lambda(\hat{w}(y,\rho)-\hat{w}(y,s))}{ps-\delta-\lambda+\lambda\hat{f_{Y}}(s)},\end{split}$

where $\hat{w}:[0,\infty)^{2}\to\mathbb{R}$ is given by

(8) $\hat{w}(y,s)=\int_{0}^{\infty}e^{-su}f_{Y}(y+u)\,du.$

From the dividends-penalty identity [12] we then have

(9)
$f_{D}(y,u,b)=f_{D^{0}}(y,u)-\frac{W_{\delta}(u)}{W_{\delta}^{\prime}(b)}D_{2}f_{D^{0}}(y,b).$

Note that we require $W_{\delta}$ to be differentiable in order to be able to
use these formulas. Nonetheless, if for some $\alpha<1$ and $C>0$, we have
$f_{Y}(x)\leq Cx^{-1-\alpha}$ for $x$ in some neighbourhood of the origin,
then $W_{\delta}\in C^{q+2}(0,\infty)$ whenever $f_{Y}\in C^{q}(0,\infty)$
(see [16]), where $C^{q}(0,\infty)$ refers to the set of $q$-times
continuously differentiable functions on the real positive line. A formula
similar to (6) shows that the statement is also valid whenever we replace
$W_{\delta}$ by $f_{D^{0}}$.
Equation (5) also reveals further properties of $V_{\pi}$ when we shift the
focus from the initial capital $u$ to the band limits: for $m\geq 1$, we can
identify the set of $m$-band strategies with the set

$\mathcal{B}_{m}=\\{x\in\mathbb{R}^{2m-1}\mid 0\leq x_{1}\leq\cdots\leq
x_{2m-1}\\},$

and, for fixed $u>0$, we can consider the function
$V^{m}:\mathcal{B}_{m}\to[0,\infty)$ given by $x\mapsto V_{x}(u)$. We have the
following:

######  Proposition 1.

If $f_{Y}\in C^{q}(0,\infty)$, $q\geq 2$, the function $V^{m}$ is continuously
differentiable in the interior of the set

$\mathcal{C}_{m}=\mathcal{B}_{m}\cap\\{x\in\mathbb{R}^{2m-1}\mid x_{2j-2}\neq
u,j=1,\ldots,m\\}.$

When $m=1$, we take $\\{x\in\mathbb{R}^{2m-1}\mid x_{2m-2}\neq u\\}$ to be
equal to $\mathbb{R}$.

###### Proof.

We proceed by induction on $m$. For $m=1$, $\mathcal{C}_{1}=[0,\infty)$ and
$\mathcal{D}_{1}=[0,u)\cup(u,\infty)$. For $0\leq u<b$, we have

$V^{1}(b)=\frac{W_{\delta}(u)}{W_{\delta}^{\prime}(b)},\quad
V^{1\prime}(b)=-\frac{W_{\delta}(u)W_{\delta}^{\prime\prime}(b)}{W_{\delta}^{\prime}(b)^{2}},$

while for $0\leq b<u$,

$V^{1}(b)=u-b+V_{b}(b)=u-b+\frac{W_{\delta}(b)}{W_{\delta}^{\prime}(b)},\quad
V^{1\prime}(b)=-\frac{W_{\delta}(b)W_{\delta}^{\prime\prime}(b)}{W_{\delta}^{\prime}(b)^{2}}.$

Since $V^{1}(u-)=V^{1}(u+)=W_{\delta}(u)/W_{\delta}^{\prime}(u)$ and
$V^{1\prime}(u-)=V^{1\prime}(u+)=W_{\delta}(u)W_{\delta}^{\prime\prime}(u)/W_{\delta}^{\prime}(u)^{2}$,
the claim follows. Now, assume the claim is true for some $m\in\mathbb{N}$. We
can write $\mathcal{C}_{m+1}$ as $\mathcal{C}_{m+1}=A\cup B\cup C$, where

$\displaystyle A$ $\displaystyle=\\{x\in\mathcal{C}_{m+1}\mid u<x_{2m-2}\\},$
$\displaystyle B$ $\displaystyle=\\{x\in\mathcal{C}_{m+1}\mid x_{2m-2}<u\leq
x_{2m-1}\\},$ $\displaystyle C$ $\displaystyle=\\{x\in\mathcal{C}_{m+1}\mid
x_{2m-1}<u\\}.$

We observe the following: on $A$, $V^{m+1}=V^{m}\circ\pi$, where
$\pi:\mathbb{R}^{2m+1}\to\mathbb{R}^{2m-1}$ is the projection onto the first
$2m-1$ coordinates. Since $\pi$ maps the interior of $A$ into the interior of
$\mathcal{C}_{m}$, by the induction hypothesis, it follows that $V^{m+1}$ is
continuously differentiable in that set. Now, from (5) we have, for
$b=(b_{0},a_{1},\ldots,b_{m})$ in $B$

$\begin{split}V(b)=\frac{W_{\delta}(u-a_{m})}{W_{\delta}^{\prime}(b_{m}-a_{m})}+\int_{0}^{a_{m}}V^{m}_{y}(b)f_{D}(a_{m}-y,b_{m}-a_{m},b_{m}-a_{m})dy,\end{split}$

while for $b$ in $C$,

$\begin{split}V(b)=u-b_{m}+\frac{W_{\delta}(b_{m}-a_{m})}{W_{\delta}^{\prime}(b_{m}-a_{m})}+\int_{0}^{a_{m}}V^{m}_{y}(b)f_{D}(a_{m}-y,b_{m}-a_{m},b_{m}-a_{m})dy.\end{split}$

From (7) and (9) we see that under the assumptions of the proposition, $f_{D}$
is twice continuously differentiable. Hence, using the induction hypothesis
once again, it follows that $V$ is continuously differentiable in the
interiors of $B$ and $C$. Moreover, since

$\displaystyle V(b_{0},a_{1},\ldots,u-)$
$\displaystyle=V(b_{0},a_{1},\ldots,u+),$ $\displaystyle
D_{2m}V(b_{0},a_{1},\ldots,u-)$
$\displaystyle=D_{2m}V(b_{0},a_{1},\ldots,u+),$

and

$D_{2m+1}V(b_{0},a_{1},\ldots,u-)=D_{2m+1}V(b_{0},a_{1},\ldots,u+),$

we conclude that $V$ is continuously differentiable in the interior of
$\mathcal{C}_{m+1}$, concluding the proof. ∎

######  Remark 4.1.

By considering instead the set

$\mathcal{D}_{m}=\mathcal{B}_{m}\cap\\{x\in\mathbb{R}^{2m-1}\mid x_{j}\neq
u,j=1,\ldots,2m-1\\},$

we can, in a similar manner, conclude that if $f_{Y},W_{\delta},f_{D^{0}}\in
C^{q}(0,\infty)$, $q\geq 2$, then $V^{m}$ is $q-1$ times differentiable in the
interior of $\mathcal{D}_{m}$, since in this case one does not have to
consider the “pasting” at the points where $b_{j}=u$.

##  5\. A gradient-based method

From Proposition 1 and its proof, we can compute the gradient of the value
function when we fix the initial capital $u$ and we look at it as a function
of the levels. Assuming that $W_{\delta}$ and $f_{D^{0}}$ are twice
differentiable and setting

$\mathcal{E}_{m}=\mathcal{B}_{m}\cap\\{x\in\mathbb{R}^{2m-1}\mid
x_{2m-2}<u\\},$

we have for $b\in\mathcal{E}_{m}$ and $u<b_{m-1}$,

(10)
$\displaystyle\begin{split}D_{2m-1}V^{m}(b)&=-\frac{W_{\delta}(u-a_{m-1})W_{\delta}^{\prime\prime}(b_{m-1}-a_{m-1})}{W_{\delta}^{\prime}(b_{m-1}-a_{m-1})^{2}}\\\
&\quad+\int_{0}^{a_{m-1}}D_{3}f_{D}(a_{m-1}-y,b_{m-1}-a_{m-1},b_{m-1}-a_{m-1})V^{m-1}_{y}(b)dy,\end{split}$
(11)
$\displaystyle\begin{split}D_{2m-2}V^{m}(b)&=\frac{W_{\delta}(u-a_{m-1})W_{\delta}^{\prime\prime}(b_{m-1}-a_{m-1})-W_{\delta}^{\prime}(u-a_{m-1})W_{\delta}^{\prime}(b_{m-1}-a_{m-1})}{W_{\delta}^{\prime}(b_{m-1}-a_{m-1})^{2}}\\\
&\quad+\int_{0}^{a_{m-1}}D_{1}f_{D}(a_{m-1}-y,b_{m-1}-a_{m-1},b_{m-1}-a_{m-1})V^{m-1}_{y}(b)dy\\\
&\quad\quad-\sum_{i=2}^{3}\int_{0}^{a_{m-1}}D_{i}f_{D}(a_{m-1}-y,b_{m-1}-a_{m-1},b_{m-1}-a_{m-1})V^{m-1}_{y}(b)dy\\\
&\quad\quad\quad+f_{D}(0,b_{m-1}-a_{m-1},b_{m-1}-a_{m-1})V^{m-1}_{y}(a_{m-1}),\end{split}$

while for $u>b_{m-1}$

(12)
$\displaystyle\begin{split}D_{2m-1}V^{m}(b)&=-\frac{W_{\delta}(b_{m-1}-a_{m-1})W_{\delta}^{\prime\prime}(b_{m-1}-a_{m-1})}{W_{\delta}^{\prime}(b_{m-1}-a_{m-1})^{2}}\\\
&\quad+\int_{0}^{a_{m-1}}D_{3}f_{D}(a_{m-1}-y,b_{m-1}-a_{m-1},b_{m-1}-a_{m-1})V^{m-1}_{y}(b)dy,\end{split}$
(13)
$\displaystyle\begin{split}D_{2m-2}V^{m}(b)&=\frac{W_{\delta}(b_{m-1}-a_{m-1})W_{\delta}^{\prime\prime}(b_{m-1}-a_{m-1})}{W_{\delta}^{\prime}(b_{m-1}-a_{m-1})^{2}}-1\\\
&\quad+\int_{0}^{a_{m-1}}D_{1}f_{D}(a_{m-1}-y,b_{m-1}-a_{m-1},b_{m-1}-a_{m-1})V^{m-1}_{y}(b)dy\\\
&\quad\quad-\sum_{i=2}^{3}\int_{0}^{a_{m-1}}D_{i}f_{D}(a_{m-1}-y,b_{m-1}-a_{m-1},b_{m-1}-a_{m-1})V^{m-1}_{y}(b)dy\\\
&\quad\quad\quad+f_{D}(0,b_{m-1}-a_{m-1},b_{m-1}-a_{m-1})V^{m-1}_{y}(a_{m-1}),\end{split}$

and, in both cases, for $1\leq i<2m-2$,

(14)
$D_{i}V^{m}(b)=\int_{0}^{a_{m-1}}f_{D}(a_{m-1}-y,b_{m-1}-a_{m-1},b_{m-1}-a_{m-1})D_{i}V^{m-1}_{y}(b)dy.$

We can solve these equations in an iterative manner to find candidate levels
for the optimal band strategy whenever it is finite: call
$b_{0}^{*}<a_{1}^{*}<\ldots<b_{M-1}^{*}$ the levels of the optimal strategy
and assume for the moment that $b_{0}^{*}>0$. Since the first level has to
occur at the largest global minimum of $W_{\delta}^{\prime}$, we have
$b_{0}^{*}=\sup\\{x\geq 0\mid W_{\delta}^{\prime}(x)=\inf_{y\geq
0}W_{\delta}^{\prime}(y)\\}$. We observe that
$W_{\delta}^{\prime\prime}(b_{0}^{*})=0$, so $D_{1}V^{1}(b_{0}^{*})=0$,
independently of $u$.
If the barrier strategy at $b_{0}^{*}$ is not optimal and $u$ is such that
$b_{1}^{*}<u<a_{2}^{*}$, $b_{0}^{*}$, $a_{1}^{*}$ and $b_{1}^{*}$ solve (12)
to (14) when $m=2$, since the global maximum is attained at the 2-band
strategy with levels $b_{0}^{*}$, $a_{1}^{*}$ and $b_{1}^{*}$. Hence, since
(14) is zero regardless of $a_{1}$ and $b_{1}$ when $b_{0}=b_{0}^{*}$, we see
that $a_{1}^{*}$ and $b_{1}^{*}$ are within the solutions to (12) and (13).
Moreover, these equations can be solved without regards to $u$, and, if we
optimally choose the solution (so that we end up obtaining $a_{1}^{*}$ and
$b_{1}^{*}$) we will see that $b_{0}^{*}$, $a_{1}^{*}$ and $b_{1}^{*}$ solve
(10) to (14) regardless of the value of $u>a_{1}^{*}$.
Continuing in this fashion, if a two-band strategy is not optimal, for $m=3$
and $i=2,3$, we have,

$\begin{split}D_{i}V(b)&=\int_{0}^{a_{2}}f_{D}(a_{2}-y,b_{2}-a_{2},b_{2}-a_{2})D_{i}V^{2}_{y}(b)dy\\\
&=\int_{0}^{a_{1}^{*}}f_{D}(a_{2}-y,b_{2}-a_{2},b_{2}-a_{2})D_{i}V^{2}_{y}(b)dy\\\
&\quad+\int_{a_{1}^{*}}^{a_{2}}f_{D}(a_{2}-y,b_{2}-a_{2},b_{2}-a_{2})D_{i}V^{2}_{y}(b)dy.\end{split}$

On the interval $(0,a_{1}^{*})$, $V^{2}(b)=V^{1}(b)$ as functions of the
initial capital, so $D_{i}V^{2}(b)=0$. By the remarks of the previous
paragraph, we also have $V^{2}_{y}(b)=0$ for all $y>a_{1}^{*}$, so we see that
(14) is always zero regardless of the value of $a_{2}$ and $b_{2}$. Hence, we
can proceed again by solving (12) and (13), choosing an optimal solution and
test whether we proceed further with another band. At the $m+1$-th step, the
equations that we need to solve can be written in a simpler form as

(15) $\displaystyle 0$
$\displaystyle=\frac{W_{\delta}(b_{m}-a_{m})W_{\delta}^{\prime\prime}(b_{m}-a_{m})}{W_{\delta}^{\prime}(b_{m}-a_{m})^{2}}-\int_{0}^{a_{m}}D_{3}f_{D}(a_{m}-y,b_{m}-a_{m},b_{m}-a_{m})V^{m}_{y}(b^{*})dy,$
(16) $\displaystyle 1$
$\displaystyle=f_{D}(0,b_{m}-a_{m},b_{m}-a_{m})V^{m}_{a_{m}}(b)+\int_{0}^{a_{m}}D_{1}f_{D}(a_{m}-y,b_{m}-a_{m},b_{m}-a_{m})V^{m}_{y}(b^{*})dy$

where $b^{*}=(b_{0}^{*},a_{1}^{*},\ldots,b_{m-1}^{*})$.
If $b_{0}^{*}=0$, we can discard the equation for $D_{1}V^{m}$ and work
instead on the interior of the set

$\mathcal{E}^{\prime}_{m}=\\{x\in\mathbb{R}^{2m-2}\mid b_{0}^{*}\leq
x_{1}\leq\cdots\leq x_{2m-3}\leq\min(u,x_{2m-2})\\}$

by realizing that (10) to (14) for $2\leq i<2m-2$ carry over verbatim to this
situation. The same argument then shows that we can use the same procedure for
obtaining the optimal levels. The procedure is described in Algorithm 2.

1

2

3

4

Input : Scale function $W_{\delta}$ and density of deficit $f_{D}$.

Output : Levels $B^{*}=(b_{0}^{*},a_{1}^{*},\ldots,b_{M-1}^{*})$ of the best
band strategy.

5

6begin

7

8 $m:=0$;

9 $b_{0}:=\sup\left\\{x\geq 0\mid W_{\delta}^{\prime}(x)=\inf\limits_{y\geq
0}W_{\delta}^{\prime}(y)\right\\}$;

10 $V^{0}(x):=\begin{cases}W_{\delta}(x)/W_{\delta}^{\prime}(b_{0})&x\leq
b_{0}\\\ V^{0}(b_{0})+x-b_{0}&x>b_{0}\end{cases}$;

11 while  _$\mathcal{L}(V^{m})_{x} >0$ for some $x>b_{m}$_ do

12 $m:=m+1$;

13
$B^{*}=\,$solve(_$D_{2m}V(b_{0}^{*},a_{1}^{*},\ldots,b_{m-1}^{*},a_{m},b_{m})=0,$
$D_{2m+1}V(b_{0}^{*},a_{1}^{*},\ldots,b_{m-1}^{*},a_{m},b_{m})=0$_);

14 $a_{m}^{*},b_{m}^{*}:=\;$select(_$B^{*}$_);

15 $V^{m}:=\,$defineV(_$b_{0}^{*},a_{1}^{*},\ldots,b_{m}^{*}$_);

16

17 end while

18

19 end

20

Algorithm 2 Gradient-based algorithm for optimal dividends

Algorithm 2 starts by obtaining the first level $b_{0}^{*}$ of the optimal
band strategy. If the barrier strategy at $b_{0}^{*}$ is optimal, the
algorithm finishes. Otherwise, the algorithm enters into its main loop. After
updating the number of bands, the loop proceeds to an application of an
abstract solve function in Line 2, which is used to solve simultaneously
Equations (12) and (13) (or equivalently, (15) and (16)) when the first $2m-1$
variables of $D_{2m}V$ and $D_{2m+1}V$ are fixed. In Line 2, the best levels
$a_{m-1}^{*}$ and $b_{m-1}^{*}$ are chosen by selecting the couple that
produces the best value for the value function when the initial capital is set
to $b_{m-1}^{*}$. Line 2 uses the function defineV for creating the value
function of the strategy with the levels
$(b_{0}^{*},a_{1}^{*},\ldots,b_{m-1}^{*})$ found so far, finalizing the loop.
As explained before, Algorithm 2 can be considered as a particular
implementation of the algorithm proposed by Avram et al. [5] but obtained
after trying to solve the gradient equations in a sort of ”backward” way. It
is similar to Algorithm 1 and Algorithm 4 in the sense that it is not possible
to determine beforehand whether a finite band strategy is optimal or not. The
advantage is, however, that one avoids having to fully specify the solutions
to the HJB equation as in Algorithm 1 while also avoiding the randomness
involved in Algorithm 4.

######  Remark 5.1.

Note that, in principle, one could derive explicit expressions that the
optimal band levels should satisfy by means of equations (15) and (16).
However, even in the simple case of the claims following an Erlang(2)
distribution, one already arrives at rather complicated expressions which
involve combinations and products of exponentials and polynomials, and the
resulting levels cannot be given in terms of elementary functions. However,
these equations can still be solved through numerical methods, which is the
basis of the gradient-based technique that is going to be introduced in the
sequel.

##  6\. Evolutionary Strategies

Evolutionary strategies belong to a class of nature-inspired optimization
algorithms which intend to mimic biological evolution by means of procedures
roughly categorized as mutation, recombination and selection procedures which
incorporate tasks that resemble the way evolution is carried out in nature.
Starting with a set of candidate solutions (called the parental population),
one produces a new set of candidate solutions by means of recombination and,
through mutation, randomly alters it to form a second set of candidate
solutions (called the offspring population). One then uses selection to filter
out the best candidate solutions from these two populations and iterates the
process, replacing the previous parental population with the new population
thus obtained. In general, recombination, mutation and selection tasks are
problem-dependent and adjusted according to a diverse set of criteria. ES are
classically referred by the way the offspring population is generated, and the
notation for expressing it is the
$(\mu_{\mathrm{ES}}/\rho_{ES}\overset{+}{,}\lambda_{ES})$-notation (the
subscripts ES are used here to differentiate these symbols from the previously
defined $\mu$ and $\lambda$ in earlier sections). In this notation,
$\mu_{\mathrm{ES}}$ refers to the size of the parental population at the
beginning and end of each iteration, $\rho_{ES}$ refers to the amount of
parents involved in the creation of one single offspring, randomly chosen
without replacement, and $\lambda_{\mathrm{ES}}$ refers to the number of
offsprings created in each iteration. The symbols “+” and “,” refer to the way
selection is carried out: the first one indicates that the $\mu_{\mathrm{ES}}$
members of the new parental population are going to be extracted from a set
obtained by merging the parental population together with the offspring
population, while the symbol “,” indicates that the parental population is
discarded after the creation of the $\lambda_{\mathrm{ES}}$ offsprings (so
that in a $(\mu_{\mathrm{ES}}/\rho_{ES},\lambda_{ES})$ strategy we necessarily
require $\lambda_{\mathrm{ES}}\geq\mu_{\mathrm{ES}}$).

The pseudo-code for the algorithm from Beyer [10] is presented below in
Algorithm 3, formulated in terms of a maximization problem. As suggested
before, the basic objects handled by ES are populations, which in this
strategy are modeled by tuples of the form $(x,s,f(x))$. In this
representation, $x$ is simply a candidate solution belonging to the search
space $X$. The element $s$ is used as a set of parameters aiding in the
mutation procedure of the members of the population and leading the self-
adaptive properties of the strategy. The last element is the value of the
function to optimize at $x$, which needs to be stored in order to select
elements in each iteration.

1

2

3

Input : Function $f$ to maximize in unconstrained object space $E$.

Output : Solution to problem $x^{\ast}=\operatorname*{argmax}\limits_{x\in
E}f(x)$

4 begin

5 $g:=0$;

6 initialize(_$\mathscr{P}^{(0)}:=\\{(x_{0,k},s_{0,k},f(x_{0,k}))\mid
k=1,\ldots,\mu_{\mathrm{ES}}\\}$_);

7 repeat

8 for  _$l:=1$ to $\lambda_{\mathrm{ES}}$_ do

9 $\mathscr{S}_{l}:=\,$sample(_$\mathscr{P}^{(g)},\rho_{\mathrm{ES}}$_);

10 $\tilde{s}_{l}:=\,$s_recombination(_$\mathscr{S}_{l}$_);

11 $\tilde{x}_{l}:=\,$x_recombination(_$\mathscr{S}_{l}$_);

12 $\tilde{s}_{l}^{\prime}:=\,$s_mutation(_$\tilde{s}_{l}$_);

13 $\tilde{x}_{l}^{\prime}:=\,$x_mutation(_$\tilde{x}_{l}$
,$\tilde{s}_{l}^{\prime}$_);

14 $\tilde{F}_{l}:=F(\tilde{x}_{l}^{\prime})$;

15

16 end for

17
$\mathscr{O}^{(g)}:=\\{(\tilde{x}_{l}^{\prime},\tilde{s}_{l}^{\prime},\tilde{F}_{l})\mid
l=1,\ldots,\lambda_{\mathrm{ES}}\\}$;

18 $\mathscr{P}^{(g+1)}:=\,$Selection(_$\mathscr{P}^{(g)}$
,$\mathscr{O}^{(g)}$,$\mu_{\mathrm{ES}}$_);

19 $g:=g+1$

20 until  _$\mathrm{terminal\\_condition}$_ ;

21 end

22

Algorithm 3 The basic ES-algorithm

After initialization of the algorithm, which is usually carried out randomly,
the algorithm enters into the main loop of the strategy for generating
subsequent populations. This loop can roughly be described as an alternation
of creating offsprings out of the parental population and selecting the
replacing parental population out of these offsprings.
The process for creating offsprings is carried out from Line 3 to 3. Lines 3
to 3 carry out the recombination procedure by first extracting a subsample
from $\mathscr{P}^{(g)}$ of size $\rho_{\mathrm{ES}}$ without replacement. In
real-valued spaces, a common recombination operator is the arithmetic mean,
so, for example,
$\texttt{s\\_recombination}(\mathscr{S}_{l})=\rho_{\mathrm{ES}}^{-1}\sum_{\mathscr{S}_{l}^{s}}s_{k}^{g}$
where $\mathscr{S}_{l}^{s}$ is the set of $s_{k}^{g}$ that belong to some
tuple in $\mathscr{S}_{l}$. The mutation operator is then applied in lines 3
to 3 by first mutating the strategy parameters and then the candidate
solutions afterwards. While there is no established methodology for choosing
the mutation operator, in [9], Beyer suggests that any operator should satisfy
three requirements for a successful implementation of ES: scalability (the
ability to tune the strength of the mutation), reachability (the ability to
reach any other state $(x,s)$ within a finite number of steps) and
unbiasedness. Scalability is achieved by allowing the mutation of the object
states, $x$, to be dependent on $s$. For $N$-dimensional real-valued search
spaces, the parameter $s$ is generally used for controlling the variance of
the mutation and in this regard, theoretical and practical considerations lead
to a common mutation operator given by

(17) $\texttt{s\\_mutation}(s_{l})^{j}=s_{l}^{j}\exp\left(\tau
N_{j}\right),\quad j=1,\ldots,N,$

where $\tau$ is the learning-rate parameter and $N_{j}$ is a standard normally
distributed random variable. Given the current parental state, the
unbiasedness requirement simply means that the mutation procedure should not
introduce any bias, and following the so-called maximum entropy principle,
this requirement immediately leads to mutation operators given by

(18)
$\texttt{x\\_mutation}(x_{l},\tilde{s}_{l})^{j}=x_{l}^{j}+\tilde{s}_{l}^{j}Z_{j},\quad
j=1,\ldots,N$

with $Z_{j}$ a standard normal random variable independent from the variables
used to mutate $s$. However, Yao et al. suggest in [26] and [18] that, more
generally, allowing $Z_{j}$ to have other kinds of stable distributions
improves convergence speed and deals better with problems where several local
extrema exist, dealing at once as well with a better handling of the
reachability requirement.
The creation of the offspring is finalized in Line 3 by evaluating the
objective function in the mutated objects $x$ and the offspring population is
gathered in Line 3. The last step of the main loop is achieved in Line 3,
where the desired selection (plus or comma) takes place and the new parental
population is created.
Figure 2 illustrates one iteration of a $(10/5+5)$ evolutionary strategy in a
two-dimensional real space for the function
$f(x,y)=\sqrt{\max\\{10-(x-10)^{2}/2-(y-10)^{2}/2,0\\}}$. Mutations for the
parameters occur as in equations (17) and (18) with $N_{j}$ independent
standard normal variables. Figure 2a represents the state of the population at
the beginning of the iteration, where, in the notation of Algorithm 3, the
$x_{k}^{g}$ are shown as the center of the ellipses, the $s_{k}^{g}$ as their
axes and, using a blue-black-red scale, each point and ellipse is colored
according to the value of $f$ at $x_{k}^{g}$. Figure 2b shows the first step
in the creation of a single offspring: after randomly selecting 5 individuals
from the original population (marked by the 5 darkest ellipses), the olive-
colored ellipse is created after applying the recombination operator. In this
case, recombination is given by the arithmetic mean, so that the olive point
and the axes of the olive ellipse are the arithmetic means of the other 5
points and the ellipses’ axes respectively. After recombination, mutation
takes places, which is represented in Figure 2c. Here, the axes of the olive
ellipse are mutated according to equation (17), generating the green ellipse.
Sampling from the normal distribution centered in the olive point and variance
given by the green ellipse, the green dot is created. Conclusion of the
offspring individual’s creation is depicted in Figure 2d, where the olive
point and ellipse are deleted and the green ellipse is “associated” with the
green point. After all offspring individuals have been created and their
values according to $f$ have been computed, the parent and offspring
population are merged, which is shown in Figure 2e. Finally, the plus
selection operator is used to discard the individuals with the lowest values
of $f$, which finishes the iteration.

a Initial Population

b Recombination

c Mutation 1

d Mutation 2

e Evaluation

f Selection

Figure 2. Illustration of a $(10/5+5)$ ES.

The presentation of the evolutionary strategy given so far raises two
questions. First, in the case of real search spaces with the mutations defined
in (18), Algorithm 3 does not consider possible constraints imposed on the
search space. A solution to this is the incorporation of restriction-handling
techniques, like the inclusion of penalty functions, reparation of offspring,
multiobjective optimization, etc., and the algorithm has to be adapted
accordingly to the technique used (see [15] for an overview of several
constraint-handling techniques in the context of ES).
The second question concerns the convergence of the algorithm. While
theoretical results exist ensuring the almost sure convergence of the
iterations (see, e.g., [22]), the assumptions used in the statements of such
results are usually quite restrictive or require a deep knowledge of the
explicit form of the objective function, which leads to difficulties at the
moment of the implementation. Despite this, evolutionary algorithms have been
tested in a wide set of scenarios, proving to be effective tools for solving
optimization problems.

##  7\. ES for the optimal dividend bands problem

It is now of interest to see how competitive evolutionary strategies for the
numerical determination of optimal band levels are in the present context.
Following the considerations from the previous section, in Section 8 we will
use a $(\mu_{\mathrm{ES}}/\rho_{\mathrm{ES}}+\lambda_{\mathrm{ES}})$
evolutionary algorithm to find the optimal band strategy for three distinct
claim distributions: mixtures of Erlang distributions, a shifted-Pareto and a
mixture of shifted-Pareto and Erlang distributions. As discussed earlier, only
very few instances of explicit non-barrier optimal band strategies have been
identified. In a what has become a classical example by now, Azcue and Muler
[6] identified a 2-band strategy for a case with Erlang(2,1) claims. Adding to
this, in [8] Berdel managed to expand this work by developing an algorithm for
identifying non-barrier band strategies in the case of a mixture of Erlang
distributions and some more general phase-type distributions. Our selection of
mixtures of Erlang distributions for testing the ES was therefore made to
compare its efficacy against an established baseline. Further, as can be seen
from Section 8.1, the lack of explicit formulas for the scale function in the
case of Pareto claims imposes the need of numerical approximations to the
evaluation of $V_{\pi}$. As shown in Loeffen [20], for any choice of
parameters, a barrier strategy is the optimal one for a Cramér-Lundberg model
with shifted-Pareto claims. The second choice of claim distribution for the
present work was then made to test the ES in a numerically-driven situation
and test its respective efficacy. Finally, the mixture of Erlang and shifted-
Pareto claim distribution was used as a means of testing the algorithm in
uncharted territory.
Algorithm 4 displays the ES-algorithm adapted for the dividend-bands
optimization, where $\odot$ stands for the element-wise multiplication
operator, and $\bar{1}_{k}$ and $\bar{0}_{k}$ are the $k$-dimensional vectors
of ones and zeros, respectively.

1

2

3

4

Input : Initial capital $u_{0}$, upper bound on number of bands $M$, number of
generations $G$ and variance bound $\varepsilon$ and value function
$V_{\cdot}$.

Output : Levels $B^{*}=(b_{0}^{*},a_{1}^{*},\ldots,b_{M-1}^{*})$ of the best
$M$-band strategy.

5 begin

6 for  _$m:=1$ to $M$_ do

7 $g:=0$;

8 initialize(_$\\{B_{0,k}\mid k=1,\ldots,\mu_{\mathrm{ES}}\\}$_);

9 $s_{0,k}:=1.0,\;k=1,\ldots,\mu_{\mathrm{ES}}$;

10 $\mathscr{P}^{(0)}:=\\{(B_{0,k},s_{0,k},V_{B_{0,k}}(u_{0}))\mid
k=1,\ldots,\mu_{\mathrm{ES}}\\}$;

11 repeat

12
$\overline{s}=\frac{1}{\mu_{\mathrm{ES}}}\sum_{i=1}^{\mu_{\mathrm{ES}}}s_{g,i}$;

13
$\overline{B}=\frac{1}{\mu_{\mathrm{ES}}}\sum_{i=1}^{\mu_{\mathrm{ES}}}B_{g,i}$;

14 $R:=\;$random_normal(_$1$_);

15 for  _$l:=1$ to $\lambda_{\mathrm{ES}}$_ do

16 $s_{R}:=\;$random_normal(_$2m-1$_);

17 $B_{R}:=\;$random_normal(_$2m-1$_);

18
$\tilde{s}_{l}^{\prime}:=\,\overline{s}\odot\exp\left(s_{R}/\sqrt{4m-2}+R\cdot\bar{1}_{2m-1}/\sqrt{2\sqrt{4m-2}}\right)$;

19 $\tilde{B}_{l}^{\prime}:=\,\max(\tilde{s}_{l}^{\prime}\odot
B_{R}+\overline{B},0_{2m-1})$;

20 $\tilde{s}_{l}^{\prime\prime}:=\,$sort(_$\tilde{s}_{l}^{\prime}$
,order(_$\tilde{B}_{l}^{\prime}$_)_);

21 $\tilde{B}_{l}^{\prime\prime}:=\,$sort(_$\tilde{B}_{l}^{\prime}$_);

22 $\tilde{V}_{l}:=V_{\tilde{B}_{l}^{\prime\prime}}(u_{0})$;

23

24 end for

25
$\mathscr{O}^{(g)}:=\\{(\tilde{B}_{l}^{\prime\prime},\tilde{s}_{l}^{\prime\prime},\tilde{V}_{l})\mid
l=1,\ldots,\lambda_{\mathrm{ES}}\\}$;

26 $\mathscr{P}^{(g+1)}:=\,$selection(_$\mathscr{P}^{(g)}$
,$\mathscr{O}^{(g)}$,$\mu_{\mathrm{ES}}$_);

27 $g:=g+1$

28 until  _$g=G$ or $\max(s_{g,0})<\varepsilon$_;

29

30 end for

31

32 end

33

Algorithm 4 ES-algorithm for optimal dividends

Algorithm 4 is a
$(\mu_{\mathrm{ES}}/\mu_{\mathrm{ES}}+\lambda_{\mathrm{ES}})$-ES based on the
basic strategy described in [10] for a search in a real unconstrained object
space. Initialization is carried out in Lines 4 to 4, where the function
initialize stands for random initialization of the candidate levels.
Recombination of the parental population is done in Lines 4 and 4 using the
arithmetic mean. In Lines 4, 4 and 4 the call of the function
$\texttt{random\\_normal}(k)$ represents the creation of an independent vector
of dimension $k$ of standard normal random variables. Lines 4 to 4 show the
implementation of the mutation operator, where each parent individual produces
one offspring individual using log-normal multiplicative mutations for the
exogenous parameters and normal mutations for the object parameters. The
coefficients $1/\sqrt{4m-2}$ and $R\cdot\bar{1}_{2m-1}/\sqrt{2\sqrt{4m-2}}$
are learning rates which depend on the dimension of the search space and are
based on both theoretical and empirical investigations. After the offspring
has been created, repairing is carried out to ensure that the levels satisfy
the condition $0\leq b_{0}\leq a_{1}<\cdots\leq b_{m-1}$. Line 4 sorts the $s$
parameters according to the increasing order of $B$, while Line 4 sorts the
object parameters in increasing order. Finally, the function in Line 4
performs plus selection and outputs the population
$\mathscr{P}^{(g+1)}:=\\{(B_{g+1,k},s_{g+1,k},V_{B_{g+1,k}}(u_{0}))\mid
k=1,\ldots,\mu_{\mathrm{ES}}\\}$ ordered in decreasing order according to the
value of $V$, so when the terminal condition in Line 4 is evaluated, $s_{g,0}$
holds the variances of the levels with the best fit.
Notice that the algorithm requires a value for the initial capital. While in
principle this is a technical condition for the evaluation of $V_{\pi}$,
caution should be taken: in case the optimal band strategy $\pi^{\ast}$ is
finite with levels $b_{0}^{*}\leq a_{1}^{*}<\cdots<b_{m-1}^{*}$, for
$b_{i}^{*}\leq u\leq a_{i+1}$, any other band strategy $\pi$ with first $i+1$
bands given by $b_{0}^{*}\leq a_{1}^{*}<\cdots<b_{i}^{*}$ and $u_{0}<a_{i+1}$
will satisfy $V_{\pi^{*}}(u_{0})=V_{\pi}(u_{0})$. Therefore, unless we can
ensure $b_{m-1}^{*}\leq u_{0}$, any such $\pi$ will be the output of any
optimization algorithm for which the initial capital is fixed. Following Lemma
3.3.1 in [24], the inequality $b_{m-1}^{*}\leq u_{0}$ can be guaranteed by
taking $u_{0}=p\lambda/(\delta(\lambda+\delta))$, which is the value that we
use for all the iterations of the algorithm.
At this point, it is worthwhile to mention a key difference between the search
method employed by this algorithm and the iterative algorithm discussed in,
for example, [6] or [7]: given that the dimension of the search space has to
be kept constant during the procedure, one has to fix in advance the number of
bands for which the ES will try to identify the optimal levels. By observing
that one can “collapse” levels in a band strategy, the $n$-bands strategies
can be thought of as $m$-band strategies for $n\leq m$ and hence the algorithm
would identify at once the best levels for all $n$-band strategies for $n\leq
m$. If the optimal strategy is finite, one could then set $m$ large enough and
use the ES to find the optimal levels. However, this approach requires the
evaluation of $V_{\pi}$ for several bands and as explained in Section 8.1
below, this is not efficient. Hence, a more efficient approach is instead to
consecutively compute the the best 1-band, 2-band, 3-band, etc. strategies
until collapsing of the levels is observed and then verify, by means of the
HJB equation, that the proposed solution is in fact the optimal band strategy.
Finally, efficiency is improved by skipping the search for the optimal 1-band
strategy and set $b_{0}:=\sup\left\\{x\geq 0\mid
W_{\delta}^{\prime}(x)=\inf_{y\geq 0}W_{\delta}^{\prime}(y)\right\\}$ in all
searches.

##  8\. Numerical results

We evaluate the performance of the procedures shown in the two previous
sections by finding optimal band strategies for three study cases: claims
distributed as a mixture of Erlang distributions, the case for a pure
(shifted) Pareto distribution and a mixture between Erlang and Pareto
distributions. The mixture of Erlang distributions is chosen because there are
already explicit results available (see [6, 8]) so that we can benchmark our
algorithms. Given that no explicit expressions exist for the scale function
when the claims follow a Pareto distribution, the second case was chosen to
test the algorithms in a purely numerical situation (and in the case of a
Pareto distribution, it is known that the optimal strategy is a barrier
strategy, see [20]). Finally, the mixture of Erlang and Pareto distributions
was carried out to study the problem in a new context.

###  8.1. Objective function evaluation

In the case where $f_{Y}$ comes from a mixture of Erlang distributions, the
Laplace transform $\widehat{W_{\delta}}$ of the scale function is given in
terms of a rational function, so explicit expressions in terms of the roots of
the Lundberg equation can be found for $W_{\delta}$ and $f_{D}$. These
expressions are then used for computing the value of $V_{\pi}$.
In the other two cases, numerical inversion of $\widehat{W_{\delta}}$ and (7)
have to be carried out to find the values of $f_{D}$. Since the evaluation of
this function is needed at several points, we opted for using a piece-wise
linear approximation for $W_{\delta}$, $W_{\delta}^{\prime}$, $f_{D^{0}}$ and
its partial derivatives. The approximation was carried out in the following
way: from the remarks of Section 7, it can easily be seen that in the case
where the optimal band strategy is finite, it is enough to restrict the domain
of $f_{D}$ and $W_{\delta}$ to $[0,p\lambda/(\delta(\lambda+\delta))]^{3}$ and
$[0,p\lambda/(\delta(\lambda+\delta))]$ respectively in order to find the
optimal band levels. Hence, the approximation was done by evaluating 10,000
equidistant points in the interval $[0,p\lambda/(\delta(\lambda+\delta))]$
(including boundaries) and linear interpolation in between. The functions
$f_{D^{0}}$ and $f_{D}$ were then computed using Equations (7) and (9).
Finally, except for the mixture of Erlang distributions appearing in Case I,
the integrals appearing in (5), (12), (13) and (14) were evaluated using
numerical methods by using the numpy, scipy and mpmath libraries for Python 3.
Equations (15) and (16) were solved by means of the MINPACK’s hybrd and hybrj
algorithms implemented in the scipy library through the fsolve function.
Finally, although there are no explicit results supporting the fact, the
examples found in the literature indicate that the values for the different
band values $b_{1},b_{2},\ldots$ are close to local minima of
$W_{\delta}^{\prime}$ found after the optimal $b_{0}$.
We now apply both numerical procedures introduced in this paper to the
concrete examples. Below, the reporting times mean the clock time used to
produce the results and do not include the time used for verification of the
solution in an interval through the HJB equation, which is the same for both
procedures.

###  8.2. Case I: Erlang mixture claims

The following three examples are considered (for which we know the explicit
result already from [6] and [8] for the first two):

  * •

An $\mathrm{Erlang}(2,1)$ distribution with parameters $\lambda=10$,
$\delta=0.1$ and $\eta=0.07$.

  * •

A mixture of the distributions $\mathrm{Erlang}(2,10)$, $\mathrm{Erlang}(3,1)$
and $\mathrm{Erlang}(4,0.1)$ with weights $0.025,0.225$ and $0.75$
respectively, and parameters $\lambda=1$, $\delta=0.1$ and $\eta=0.405$.

  * •

A mixture of the distributions $\mathrm{Erlang}(2,10)$,
$\mathrm{Erlang}(3,1.06775)$, $\mathrm{Erlang}(4,0.2325)$ and
$\mathrm{Erlang}(5,0.05)$ with weights $0.005,0.045$, $0.225$ and $0.725$
respectively, and parameters $\lambda=1$, $\delta=0.1$ and $\eta=0.4$.

For the first two distributions, a $(30/30+60)$-ES was used in both cases,
with bound on the variance equal to 0.01. For reasons that will be explained
later, a combination of a $(30/30+60)$-ES and a $(1/1+1)$-ES was used for the
third distribution, with same bound in the variance. The number of iterations
vary from distribution to distribution.
For the Erlang(2,1) distribution, Table 1 shows that we indeed find the
optimal two-band strategy established in [6].

| Time | Iterations | $b_{0}$ | $a_{1}$ | $b_{1}$
---|---|---|---|---|---
Evolutionary Strategy |  $\sim 1000$s | 1000 | 0 | 1.8064 | 10.2158
Gradient-Based |  $<1$s | – | 0 | 1.8030 | 10.2161
Table 1. Results for the $\mathrm{Erlang}(2,1)$ distribution.

One can observe that the gradient-based approach is very fast, while the ES
algorithm takes considerably longer time, but also arrives at the correct
solution.
Table 2 shows the results for the first mixture of distributions, where a
3-band strategy is optimal.

| Time | Iterations | $b_{0}$ | $a_{1}$ | $b_{1}$ | $a_{2}$ | $b_{2}$
---|---|---|---|---|---|---|---
Evolutionary Strategy |  $\sim 1$h | $5000$ | $0.2617$ | $0.4668$ | $3.5249$ | $25.7390$ | $34.7857$
Gradient |  $\sim 1$h | – | $0.2615$ | $1.5230$ | $3.5246$ | $25.5763$ | $34.7696$
Table 2. Results for the first mixture of Erlang distributions.

Several remarks for this example are in order. First, there is a clear
discrepancy between the $a_{1}$-values obtained by the two methods. Comparing
the value functions of both strategies shows that the strategy found by the ES
provides higher values. However, the difference between the value functions is
of the order of $10^{-4}$ while the norm of the gradient at that point is of
the order $10^{-5}$ and the iterations do not reduce this value significantly,
which shows why the gradient method has an early stop. Note that for this
example, Berdel [8] already studied the optimal dividend strategy, and her
results are very similar to the ones in Table 2, with only the values of
$a_{1}$ and $a_{2}$ differing. The ES parameters above provide a larger value
function, but again the difference is only of order $10^{-4}$. The explanation
is that the step size in [8] for solving the inf and sup in Lines 7 and 8 in
Algorithm 1 was set to be $10^{-4}$, whereas a smaller step size would have
been needed to arrive at the above result.
As stated before, for the evolutionary strategy, the integrals in equations
(5), (12), (13) and (14) were not evaluated numerically but instead were
computed exactly, by means of symbolical calculus in Mathematica. Now, for
$k=2$, and $u\geq b_{2}$, (5) can be more explicitly written as

$\displaystyle V_{\pi}(u)$
$\displaystyle=u-b_{2}+\frac{W_{\delta}(h_{2})}{W_{\delta}^{\prime}(h_{2})}$
$\displaystyle+\int_{0}^{b_{0}}f_{D}(a_{2}-y,h_{2},h_{2})\frac{W_{\delta}(y)}{W_{\delta}^{\prime}(b_{0})}\;dy$
$\displaystyle+\int_{b_{0}}^{a_{1}}f_{D}(a_{2}-y,h_{2},h_{2})\left(y-b_{0}+\frac{W_{\delta}(b_{0})}{W_{\delta}^{\prime}(b_{0})}\right)\;dy$
$\displaystyle+\int_{a_{1}}^{b_{1}}f_{D}(a_{2}-y,h_{2},h_{2})\frac{W_{\delta}(y-a_{1})}{W_{\delta}^{\prime}(h_{1})}\;dy$
$\displaystyle+\int_{a_{1}}^{b_{1}}f_{D}(a_{2}-y,h_{2},h_{2})\int_{0}^{b_{0}}f_{D}(a_{1}-z,y-a_{1},h_{1})\frac{W_{\delta}(z)}{W_{\delta}^{\prime}(b_{0})}\;dz\;dy$
$\displaystyle+\int_{a_{1}}^{b_{1}}f_{D}(a_{2}-y,h_{2},h_{2})\int_{0}^{b_{0}}f_{D}(a_{1}-z,y-a_{1},h_{1})\left(y-b_{0}+\frac{W_{\delta}(b_{0})}{W_{\delta}^{\prime}(b_{0})}\right)\;dz\;dy$
$\displaystyle+\int_{b_{1}}^{a_{2}}f_{D}(a_{2}-y,h_{2},h_{2})\left(y-b_{1}+V_{\pi^{1}}(b_{1})\right)\;dy$

where $h_{i}=b_{i}-a_{i}$ and $V_{\pi^{1}}$ is the 2-band strategy obtained
after deleting the last band from $V_{\pi}$. For mixtures of Erlang
distributions, the scale function can be written as a linear combination of
complex exponential functions with as many terms as roots of the Lundberg
equation, assuming all of them are different. For the present case, there are
10 different roots. By means of formula (6), it follows that $f_{D}(y,u,b)$
can be written as a sum of approximately 90 different terms involving $y$ with
coefficients dependent on $u$ and $b$. Following this line of thought, the
single integrals from the paragraph above have, in rough terms, 900 terms,
while the first double integral has around 81,000 (in theory, further
reductions that decrease these numbers considerably could in principle be
possible, but the amount of terms implies that the human or computational
effort for carrying out such operations is beyond reason). The computational
effort for explicitly computing the integrals above was of around 1 hour,
which implies that the total time for the ES was of around 2 hours, which
doubled the computational time of the gradient method, but provided a slightly
more accurate result. Although this procedure could also be carried out to
test the computational time of the gradient method (which by virtue of the
other cases would be expected to be smaller), we observe that the
computational time for obtaining explicit expressions for (12) and (13) would
be at least doubled, matching the current time of the ES.

Finally, using the intuition that the number of bands in the optimal strategy
is related to the number of modes of the claim distribution, we were
interested to establish a case where a 4-band strategy is optimal, and the
second mixture of Erlang distributions indeed leads to such an optimal 4-band
strategy. The resulting optimal bands are given in Table 3.

Time | Iterations | $b_{0}$ | $a_{1}$ | $b_{1}$ | $a_{2}$ | $b_{2}$ | $a_{3}$ | $b_{3}$
---|---|---|---|---|---|---|---|---
$\sim 3$h | $6000$ |  $0.2562$, | $1.0543$ | $3.1988$ | $10.6647$ | $19.5499$ | $127.9288$ | $171.6044$
Table 3. Results for the second mixture of Erlang distributions.

These values were computed using only the ES technique. As before, the
procedure was carried out in two steps, first using a $(30/30+60)$-ES for
computing the values of the first three bands and consecutively using these
values to reduce the problem to a two-dimensional optimization exercise, where
a $(1/1+1)$-ES was used for computing the final two values. The reported time
is for the combination of both procedures. The reason for proceeding in a two-
step fashion was due to the observation that, as shown by these examples, the
optimal values of of the $b_{i}$’s are usually located in the vicinity of the
local minima of $W_{\delta}^{\prime}$. The three smallest local minima are
found in the interval $(0,20)$, which does not present any numerical
complication. However, the last one is found at $172.7545$, which due to the
nature of the scale function, produces exponentials with very large exponents
at the moment of evaluating the value function, creating considerable
numerical instabilities. To solve this issue, the value function and
evolutionary strategy were re-implemented using arbitrary-precision floating-
point-arithmetic, which decreased the speed at the moment of evaluating the
value function. Since the $(1/1+1)$-ES is the strategy that requires the least
evaluations of the objective function, it was chosen for obtaining the final
values of $a_{3}$ and $b_{3}$. Figure 3 illustrates that this 4-band strategy
is indeed the optimal strategy. Concretely, Figures 3a,3b and 3c show that
none of the first three strategies (with 1, 2 or 3 bands) is optimal, as the
HJB equation attains positive values. Figure 3d shows that this does not
happen for the 4-band strategy and, moreover, Figure 3e reveals that whenever
the derivative exists, it is at least 1, so the solution is optimal.

a HJB equation for the 1-band strategy

b HJB equation for the 2-band strategy

c HJB equation for the 3-band strategy

d HJB equation for the 4-band strategy

e Derivative of the value function

Figure 3. Plots of the l.h.s. of the HJB equation for the strategy with 1, 2,
3, and 4 bands based on the $a_{i}$’s and $b_{i}$’s from Table 3 as well as
the derivative of the value function in the points where it exists.

###  8.3. Case II: Pure Pareto claims

Following [20], the optimal strategy will be a barrier strategy when claims
have a shifted Pareto distribution with density function

$f_{Y}(y)=\alpha x_{0}^{-1}(1+x_{0}^{-1}y)^{-\alpha-1},\quad y>0,$

and Laplace transform

$\hat{f_{Y}}(s)=\alpha
x_{0}^{\alpha}s^{\alpha}e^{sx_{0}}\Gamma(-\alpha,sx_{0}),\quad s>0,$

with $\Gamma$ the upper incomplete Gamma function and $\alpha,x_{0}>0$. For
the case at hand, we considered $x_{0}=1$ and $\alpha=1.5$, so that the claims
have finite expectation and infinite variance. Moreover, the parameters of the
Cramér-Lundberg process were taken to be $\lambda=10$, $\delta=0.1$ and
$\eta=0.1$. The derivatives of the scale function were computed through their
Laplace transforms and all of these were inverted using the de Hoog, Knight
and Stones algorithm implemented in the library mpmath. The results are given
in Table 4.

| Time | Iterations | $b_{0}$
---|---|---|---
Evolutionary Strategy |  $\sim 1000$s | 100 | 2.71036
Gradient-Based |  $<1$s | - | 2.71036
Table 4. Results for the Pareto distribution.

Indeed, one arrives at an optimal barrier strategy, where for the evolutionary
algorithm we only used 100 iterations to arrive at a running time that is
comparable to the ones of the Erlang case, and the result is already well-
aligned with the one of the gradient-based method.

###  8.4. Case III: Erlang and Pareto mixture claims

Finally, let us consider a mixture of an $\mathrm{Erlang}(2,1)$ distribution
and a shifted Pareto distribution ($\alpha=1.5$, $x_{0}=1$) with weights $0.8$
and $0.2$ respectively. The parameters of the Pareto distribution were chosen
to match the mean of the Erlang component, while the weights were chosen to
avoid a monotonicity of $W_{\delta}^{\prime}$. The other parameters are again
$\lambda=1$, $\delta=0.1$ and $\eta=0.1$. A $(150/150+100)$-ES was used and a
2-band strategy was found to be optimal. The results are shown in Table 5.

| Time | Iterations | $b_{0}$ | $a_{1}$ | $b_{1}$
---|---|---|---|---|---
Evolutionary Strategy |  $\sim 8$h | 1000 | $0$ | $0.1524$ | $3.5115$
Gradient-Based |  $\sim 1$h | - | $0$ | $0.0053$ | $3.8877$
Table 5. Results for a mixture of an Erlang and a Pareto distribution.

In this case, the discrepancy between the results is more significant than in
the other cases, with the gradient method providing a better solution. Due to
the time it took for each evaluation, the ES was stopped after 1000
iterations, which meant that convergence was not fully achieved. The results
from the gradient method also help to explain the poor performance of the ES:
observe that the value of $a_{1}$ is rather close to zero, with an order of
magnitude of $10^{-3}$. At the initialization of the algorithm, there is no
knowledge of what the final order of magnitude will be, and during this
experiment the initial values for the exogenous parameters were set to be 1,
so, due to the projection into the zero for negative values, it takes many
evaluations before the desired order of magnitude is achieved.

##  9\. Conclusions

In this paper we added two numerical alternatives to identify optimal dividend
bands in the classical optimal dividend problem of risk theory. We illustrate
that both of them are efficient, and their scope and applicability goes beyond
the one of the previously discussed methods in the literature. The gradient-
based method can be particularly efficient. The second algorithm based on
evolutionary strategies is satisfactory as well, and whereas in terms of
computation times it can not compete with the gradient-based method for the
complexity of this concrete problem, its range of applicability is even wider.
In fact, ES algorithms can be an interesting competitor whenever an objective
function can be efficiently evaluated, and it is known to work particularly
well for higher-dimensional optimization problems, in which case the gradient
alternative can be hard to explicitly compute or implement. We rederived
optimal bands for some known cases, established new ones and also derived
results for cases that were beyond the scope of previously available methods.
The focus of this paper was on optimal dividend strategies in the Cramér-
Lundberg model. However, since the equations used to derive the necessary
functions for these two algorithms were obtained by means of Gerber-Shiu
functions, one can in principle easily extend the range of applications to the
case where a diffusion is added to the surplus process or even to the case
where the surplus process is modelled by a generally spectrally-negative Lévy
process satisfying the safety loading condition. Since evolutionary algorithms
can be applied in rather general settings, it will also be interesting to see
in future research other applications of this method in risk theory,
particularly also in optimization problems with constraints, which may be
handled with an introduction of a penalty term in the objective function (see
e.g. [15]).

## References

  * [1] H. Albrecher and S. Thonhauser.  Optimality results for dividend problems in insurance.  RACSAM-Revista de la Real Academia de Ciencias Exactas, Fisicas y Naturales. Serie A. Matematicas, 103(2):295–320, 2009.
  * [2] S. Asmussen and H. Albrecher.  Ruin probabilities, volume 14.  World scientific, 2010.
  * [3] B. Avanzi.  Strategies for dividend distribution: A review.  North American Actuarial Journal, 13(2):217–251, 2009.
  * [4] F. Avram, Z. Palmowski, and M. R. Pistorius.  On the optimal dividend problem for a spectrally negative lévy process.  The Annals of Applied Probability, 17(1):156–180, 2007.
  * [5] F. Avram, Z. Palmowski, and M. R. Pistorius.  On gerber–shiu functions and optimal dividend distribution for a lévy risk process in the presence of a penalty function.  The Annals of Applied Probability, 25(4):1868–1935, 2015.
  * [6] P. Azcue and N. Muler.  Optimal reinsurance and dividend distribution policies in the cramér-lundberg model.  Mathematical Finance: An International Journal of Mathematics, Statistics and Financial Economics, 15(2):261–308, 2005.
  * [7] P. Azcue and N. Muler.  Optimal dividend policies for compound poisson processes: The case of bounded dividend rates.  Insurance: Mathematics and Economics, 51(1):26–42, 2012.
  * [8] J. Berdel.  Optimale Bandstrategien für Dividendenzahlungen im Cramér-Lundberg-Modell.  PhD thesis, Karlsruher Instituts für Technologie, 2014.
  * [9] H.-G. Beyer.  The theory of evolution strategies.  Springer Science & Business Media, 2001.
  * [10] H.-G. Beyer and H.-P. Schwefel.  Evolution strategies–a comprehensive introduction.  Natural computing, 1(1):3–52, 2002.
  * [11] H. U. Gerber.  Entscheidungskriterien für den zusammengesetzten Poisson-Prozess.  PhD thesis, ETH Zurich, 1969.
  * [12] H. U. Gerber, X. S. Lin, and H. Yang.  A note on the dividends-penalty identity and the optimal dividend barrier.  ASTIN Bulletin, 36(2):489–503, 2006.
  * [13] H. U. Gerber and E. S. Shiu.  On the time value of ruin.  North American Actuarial Journal, 2(1):48–72, 1998.
  * [14] F. Hubalek and E. Kyprianou.  Old and new examples of scale functions for spectrally negative lévy processes.  In Seminar on Stochastic Aanalysis, Random Fields and Applications VI, pages 119–145. Springer, 2011.
  * [15] O. Kramer.  A review of constraint-handling techniques for evolution strategies.  Applied Computational Intelligence and Soft Computing, 2010, 2010\.
  * [16] A. Kuznetsov, A. E. Kyprianou, and V. Rivero.  The theory of scale functions for spectrally negative lévy processes.  Lévy matters II, pages 97–186, 2012.
  * [17] A. E. Kyprianou.  Fluctuations of Lévy processes with applications: Introductory Lectures.  Springer Science & Business Media, 2014.
  * [18] C.-Y. Lee and X. Yao.  Evolutionary programming using mutations based on the lévy probability distribution.  IEEE Transactions on Evolutionary Computation, 8(1):1–13, 2004\.
  * [19] X. S. Lin, G. E. Willmot, and S. Drekic.  The classical risk model with a constant dividend barrier: analysis of the gerber–shiu discounted penalty function.  Insurance: Mathematics and Economics, 33(3):551–566, 2003.
  * [20] R. L. Loeffen.  On optimality of the barrier strategy in de finetti’s dividend problem for spectrally negative lévy processes.  The Annals of Applied Probability, pages 1669–1680, 2008.
  * [21] S. Román, A. M. Villegas, and J. G. Villegas.  An evolutionary strategy for multiobjective reinsurance optimization.  Journal of the Operational Research Society, 69(10):1661–1677, 2018\.
  * [22] G. Rudolph.  Convergence of evolutionary algorithms in general search spaces.  In Proceedings of IEEE international conference on evolutionary computation, pages 50–54. IEEE, 1996.
  * [23] S. Salcedo-Sanz, L. Carro-Calvo, M. Claramunt, A. Castañer, and M. Mármol.  Effectively tackling reinsurance problems by using evolutionary and swarm intelligence algorithms.  Risks, 2(2):132–145, 2014.
  * [24] H. Schmidli.  Optimisation in non-life insurance.  Stochastic models, 22(4):689–722, 2006.
  * [25] H. Schmidli.  Stochastic control in insurance.  Springer Science & Business Media, 2007.
  * [26] X. Yao, Y. Liu, and G. Lin.  Evolutionary programming made faster.  IEEE Transactions on Evolutionary computation, 3(2):82–102, 1999\.


# A comparative study of back propagation and its alternatives on multilayer
perceptrons

###### Abstract

The de facto algorithm for training the back pass of a feedforward neural
network is backpropagation (BP). The use of almost-everywhere differentiable
activation functions made it efficient and effective to propagate the gradient
backwards through layers of deep neural networks. However, in recent years,
there has been much research in alternatives to backpropagation. This analysis
has largely focused on reaching state-of-the-art accuracy in multilayer
perceptrons (MLPs) and convolutional neural networks (CNNs). In this paper, we
analyze the stability and similarity of predictions and neurons in MLPs and
propose a new variation of one of the algorithms.

##  1 Introduction

When training neural networks there are many decisions that are made either
manually or automatically through any of brute force, trial and error, or
theoretical assumptions. These decisions include the overall architecture of
the network, the learning rate, the initialization, and the learning algorithm
among others. Each of these decisions results in a different model.
Occasionally, as can be the case with different weight initialization
techniques, the trained network converges to the same local extrema.

Backpropagation is the most widely used algorithm to train a neural
network[1]. The development of differential, non-linear, activation made it
possible to use the chain rule to train deep neural networks.

However, it is now believed that due to the weight transpose in the backward
pass, this is not a biologically plausible algorithm [2]. This led to a
significant increase in algorithms such as feedback alignment (FA) and direct
feedback alignment (DFA) [3]. There is much research that many of these
algorithms perform nearly as well in the MNIST and CIFAR10 tests when trained
on multilayer perceptrons (MLPs) [2,3,4,6]. In this paper, we will examine not
only the accuracy and stability of these algorithms on these tasks, but also
the similarity of the predictions and neurons of models trained using these
algorithms.

##  2 Algorithms

We will be comparing 5 different training algorithms. The first,
backpropagation, is the benchmark we will be comparing everything to. The
other four are variations of feedback alignment.

###  2.1 Backpropagation

By far the mostly widely used algorithm used to train neural networks. By
backpropagating the gradient through each layer we can efficiently and
effectively train deeper networks. The backprop algorithm is as follows,

$\displaystyle\delta W_{i}$ $\displaystyle=-((W^{T}_{i+1}\delta
z_{i+1})\odot\phi^{{}^{\prime}}(z_{i}))y^{T}_{i-1}$ $\displaystyle\delta
z_{i+1}$ $\displaystyle=\diffp{\mathcal{L}}$ $\displaystyle y_{i}$
$\displaystyle=\phi(z_{i})$ $\displaystyle z_{i}$
$\displaystyle=W_{i}y_{i-1}+b_{i}$

Backpropagation tend to set the highest accuracy in MLPs, RNNS, CNNs, and
ViTs. Newer algorithms attempt to achieve this accuracy with caveats of
holding less information, not having a weight transport, or converging faster.

###  2.2 Feedback Alignment

In feedback alignment algorithms, the weight matrix W information in the back
pass is replaced with a random matrix B. This allows us to not need a weight
transport in the back pass. Shown in figure 1 below is a comparison of
training in backpropagation and feedback alignment methods.

Figure 1: Comparison of training algorithms

Notice that in all feedback alignment algorithms a random matrix B is used in
the backpass instead of the transport of the W matrix. Below are the formal
definitions of the feedback alignment algorithms we will test,

####  2.2.1 Vanilla Feedback Alignment

Feedback alignment is an algorithm designed to bypass the $W^{T}$ weight
transpose. Instead, a random matrix B with the same dimensions as $W^{T}$ for
each layer in the network is created. Although this algorithm tends to learn a
little slower in the early epochs, information is able to be learned through
this random transformation. The algorithm is as follows,

$\displaystyle\delta W_{i}$ $\displaystyle=-((B_{i+1}\delta
z_{i+1})\odot\phi^{{}^{\prime}}(z_{i}))y^{T}_{i-1}$ $\displaystyle\delta
z_{i+1}$ $\displaystyle=\diffp{\mathcal{L}}{{z_{i+1}}}$ $\displaystyle y_{i}$
$\displaystyle=\phi(z_{i})$ $\displaystyle z_{i}$
$\displaystyle=W_{i}y_{i-1}+b_{i}$

####  2.2.2 Uniform Signed Feedback Alignment

Uniform signed feedback alignment is very similar to vanilla FA. However,
unlike FA it updates the B matrix as follows, using some information from
$W^{T}$. Specifically, the algorithm is,

$\displaystyle\delta W_{i}$ $\displaystyle=-((B_{i+1}\delta
z_{i+1})\odot\phi^{{}^{\prime}}(z_{i}))y^{T}_{i-1}$ $\displaystyle\delta
z_{i+1}$ $\displaystyle=\diffp{\mathcal{L}}{{z_{i+1}}}$ $\displaystyle y_{i}$
$\displaystyle=\phi(z_{i})$ $\displaystyle z_{i}$
$\displaystyle=W_{i}y_{i-1}+b_{i}$

And after every iteration we update the B matrix as follows,

$\displaystyle B_{i}$ $\displaystyle=\text{sign}(W_{i}^{T})$

Therefore, this method is somewhat of a hybrid between the between
backpropagation and vanilla FA. Unlike in vanilla FA, the B matrix is now a
function of $W^{T}$. However, unlike backpropagation it is not a perfect
representation on the weight transport as it only contains partial information
on whether each element is non negative.

####  2.2.3 Direct Feedback Alignment

For direct feedback alignment (DFA) we directly propagate the error to each
layer. Instead of $B$ having the same dimensions as $W^{T}$ , $B_{i}$ has the
same input dimension as $W$ but the output dimension equal to the output layer

$\displaystyle\delta W_{i}$ $\displaystyle=-((B_{i}\delta
z_{n})\odot\phi^{{}^{\prime}}(z_{i}))y^{T}_{i-1}$ $\displaystyle\delta z_{n}$
$\displaystyle=\diffp{\mathcal{L}}{{z_{n}}}$ $\displaystyle y_{i}$
$\displaystyle=\phi(z_{i})$ $\displaystyle z_{i}$
$\displaystyle=W_{i}y_{i-1}+b_{i}$

DFA was created as a means to effectively use feedback alignment in larger and
deeper networks [5]. Whereas FA gradients may get lost during backwards pass
through many randomized B layers DFA attempts to resolve that noise.

####  2.2.4 Weighted Direct Feedback Alignment

We propose a variation of DFA that slows down the learning rate of in the
earlier layers by the 59 following logic,

$\displaystyle lr_{j}=lr\frac{n\sqrt{j}}{\sum_{i=1}^{n}\sqrt{i}}$

Notice that the average learning rate over all layers i to n is the starting
learning rate $lr$. Therefore, since DFA is projecting the error propagation
directly across all layers, we can slow down this change and avoid rapid
changes to earlier layers using this method. DFA is expected to rival BP in
deeper networks [6]. We propose that this change could increase stability in
the training of these deeper networks.

##  3 Procedure

In order to compare these algorithms we build models with architecture as
follows:

###  3.1 MNIST Architecture

We flattened the 28*28 images to have the following model architecture,

$\displaystyle\text{Layer}_{\text{input}}$
$\displaystyle:\text{Size}:(786,768)...\text{Activation: tanh}$
$\displaystyle\text{Layer}_{\text{hidden 1}}$
$\displaystyle:\text{Size}:(768,256)...\text{Activation: tanh}$
$\displaystyle\text{Layer}_{\text{hidden 2}}$
$\displaystyle:\text{Size}:(256,128)...\text{Activation: tanh}$
$\displaystyle\text{Layer}_{\text{output}}$
$\displaystyle:\text{Size}:(128,10)...\text{Activation: sigmoid}$

###  3.2 CIFAR10 Architecture

We flattened the 3*32*32 images to have the following model architecture,

$\displaystyle\text{Layer}_{\text{input}}$
$\displaystyle:\text{Size}:(3072,768)...\text{Activation: tanh}$
$\displaystyle\text{Layer}_{\text{hidden 1}}$
$\displaystyle:\text{Size}:(768,256)...\text{Activation: tanh}$
$\displaystyle\text{Layer}_{\text{hidden 2}}$
$\displaystyle:\text{Size}:(256,128)...\text{Activation: tanh}$
$\displaystyle\text{Layer}_{\text{output}}$
$\displaystyle:\text{Size}:(128,10)...\text{Activation: sigmoid}$

For each dataset we created 10 models of each type. For the $i$th model
$i\in[1,10]$, of each type we set the random seed equal to $i$ before
initializing the weights of the model. That way, for each type of model, the
$i$th model started as the exact same position. Therefore, any differences in
the model are solely due to the differences in the training algorithm. All
models were trained for 60 epochs with a learning rate of 1e-4. We used cross-
entropy loss for all our models.

##  4 Experimental Results

###  4.1 Accuracy Comparison

79 For a similar comparison of the training process of these algorithms, all
models were trained with a simple minor weight decay but without an optimizer.
For accuracy purposes and for bench marking we included BP with an Adam
optimizer.

Table 1: Accuracy Comparison Algorithm | MNIST | CIFAR 10
---|---|---
BP | 0.973 | 0.491
DFA | 0.961 | 0.480
DFA Weighted | 0.960 | 0.478
FA | 0.959 | 0.454
FA Uni Sign | 0.955 | 0.471

Notice that while FA and Uniform Sign FA preform on par with BP, and DFA in
MNIST, it suffers once more parameters and complexity is added. On the other
hand, DFA tends to perform very well and stays in line, albeit a little below
BP, for both tasks.

###  4.2 Prediction Similarity

####  4.2.1 Stability of the algorithms

First we will examine the similarity of the output of each model. We will look
at the cosine similarity between the top predicted class from each of the 10
iterations for each model.

Table 2: Cosine Similarity of Predictions for each Algorithm Algorithm | MNIST | CIFAR 10
---|---|---
BP | 0.98 | 0.851
DFA | 0.99 | 0.873
DFA Weighted | 0.99 | 0.877
FA | 0.99 | 0.794
FA Uni Sign | 0.99 | 0.825

From Table 2, we notice that from the MNIST results that when the model has
few parameters the algorithms all tend to converge to similar results. This
may be largely due to the fact that the MNIST accuracy is so high that the the
few pictures the models are getting wrong are ambiguous and all the models
fail to accurately classify them.

Comparing the CIFAR10 results, BP, DFA, and Weighted DFA tend to produce
relatively stable models. Especially weighted DFA when the similarity scores
among the predicted output are highest. This would suggest that Weighted DFA
tends to be less dependent on the randomization of the B matrix then the other
algorithms.

Lastly, we note that using the Uniform Sign FA, we tend to get back some of
the stability lost using vanilla FA. Vanilla FA has the most randomized
parameters and subsequently takes the lost to train and has the least
stability.

####  4.2.2 Comparison of Predictions Across Algorithms

Table 3: Cosine Similarity Across Algorithm Outputs for CIFAR10 | BP | DFA | DFA Weighted | FA | FA Uni Sign
---|---|---|---|---|---
BP | - | 0.844 | 0.853 | 0.806 | 0.831
DFA | 0.844 | - | 0.891 | 0.791 | 0.826
DFA Weighted | 0.853 | 0.891 | - | 0.791 | 0.828
FA | 0.806 | 0.791 | 0.791 | - | 0.806
FA Uni Sign | 0.831 | 0.826 | 0.828 | 0.806 | -

In Table 3, we record the average similarity scores of each instance of each
model to the other instances. We see that our weighted DFA method most
resembles the output of backpropagation followed by regular DFA. Much like,
Table 2 we see that the more stable algorithms tend to have similar outputs
compared to each other. BP is much closer to DFA methods than FA methods.

As expected methods in the same type tend to be more similar to eachother.
Noticeably, there is much less difference in Weighted DFA and DFA than there
is between FA and Uniform Sign FA.

###  4.3 Layer Similarity

Lastly, we will measure the similarity among the 4 layers of the trained
models. To do so we will once again measure the cosine similarity of the each
of the layers both within the same algorithm and across algorithms.

Given the initial layer weights produced given the random seed, we then train
a model starting with these weights using each of the five training
algorithms. Then, for each seed and for each layer we check the cosine
similarity between than layer for each of the 5 training algorithms. The
results are as follows,

Table 4: Cosine Similarity of First Layer Weights Across Algorithms | BP | DFA | DFA Weighted | FA | FA Uni Sign
---|---|---|---|---|---
BP | - | 0.993 | 0.997 | 0.989 | 0.982
DFA | 0.993 | - | 0.997 | 0.982 | 0.976
DFA Weighted | 0.997 | 0.997 | - | 0.986 | 0.980
FA | 0.989 | 0.982 | 0.986 | - | 0.973
FA Uni Sign | 0.982 | 0.976 | 0.980 | 0.973 | -

In the first layer we can see that the most algorithms, are highly correlated.
Weighted DFA tends to be the most similar to BP.

Table 5: Cosine Similarity of Third Layer Weights Across Algorithms | BP | DFA | DFA Weighted | FA | FA Uni Sign
---|---|---|---|---|---
BP | - | 0.855 | 0.863 | 0.926 | 0.971
DFA | 0.855 | - | 0.904 | 0.799 | 0.835
DFA Weighted | 0.863 | 0.904 | - | 0.807 | 0.843
FA | 0.926 | 0.799 | 0.807 | - | 0.905
FA Uni Sign | 0.971 | 0.835 | 0.843 | 0.905 | -

In the middle layers we can see that DFA starts to have some significant
differences from BP and FA. Due to the direct propagation of the error term
these methods are updated in a more distinct method than FA and BP.

Table 6: Cosine Similarity of Output Layer Weights Across Algorithms | BP | DFA | DFA Weighted | FA | FA Uni Sign
---|---|---|---|---|---
BP | - | 0.104 | 0.098 | 0.225 | 0.782
DFA | 0.104 | - | 0.857 | 0.062 | 0.107
DFA Weighted | 0.098 | 0.857 | - | 0.057 | 0.112
FA | 0.225 | 0.062 | 0.057 | - | 0.257
FA Uni Sign | 0.782 | 0.107 | 0.112 | 0.257 | -

But in the output layer shown in Table 6, the algorithms really start to have
high cosine distances to each other. For most algorithms, this layer is pretty
distinct. This layer is the only layer we see significant gaps between the DFA
and Weighted DFA. Note that for this layer, the FA methods tends to be far
more similar to BP than the DFA methods. This is especially true for Uni Sign
FA.

##  5 Conclusion

After studying various training algorithms for MLPs, we see that in terms of
similarity of weights across layers particularly in the earlier hidden layers,
stability of algorithm, and similarity of final predictions, the Weighted DFA
tends to best mimic the backpropagation algorithm. However, in terms of
accuracy it tends to slightly underperform DFA on average.

In contrast, FA algorithms tend to update the first layers faster than DFA and
BP and tend to be less stable. We also saw that FA fell off faster when we
added more complexity to the problem when we shifted from MNIST to CIFAR10.

It remains to be seen if these similarities carry over to training of
CNNs.Recent studies show that DFA performs well on CNNs in terms of accuracy
[7].

It is clear that FA algorithms cannot handle increased complexity or deeper
networks as efficiently as BP. However, both the Weighted DFA and regular DFA
algorithms seem to be plausible candidates for moderate sized MLPs.

## References

[1] Rumelhart, David E., Geoffrey E. Hinton, and Ronald J. Williams. Learning
internal representations by error propagation. California Univ San Diego La
Jolla Inst for Cognitive Science, 1985.

[2] Song, Yuhang, et al. "Can the Brain Do Backpropagation?—Exact
Implementation of Backpropagation in Predictive Coding Networks." Advances in
neural information processing systems 33 (2020): 22566-22579.

[3] Nøkland, Arild. "Direct feedback alignment provides learning in deep
neural networks." Advances in neural information processing systems 29 (2016).

[4] Akrout, Mohamed, et al. "Deep learning without weight transport." Advances
in neural information processing systems 32 (2019).

[5] Moskovitz, Theodore H., Ashok Litwin-Kumar, and L. F. Abbott. "Feedback
alignment in deep convolutional networks." arXiv preprint arXiv:1812.06488
(2018).

[6] Launay, Julien, et al. "Direct feedback alignment scales to modern deep
learning tasks and architectures." Advances in neural information processing
systems 33 (2020): 9346-9360.

[7] Han, Donghyeon, and Hoi-jun Yoo. "Efficient convolutional neural network
training with direct feedback alignment." arXiv preprint arXiv:1901.01986
(2019).