Kyoto University,
Graduate School of Economics
Research Project Center Discussion Paper Series

Solving Optimal Dividend Problems via Phase-Type Fitting
Approximation of Scale Functions

Masahiko Egami and Kazutoshi Yamazaki

Discussion Paper No. E-10-011

Research Project Center
Graduate School of Economics
Kyoto University
Yoshida-Hommachi, Sakyo-ku
Kyoto City, 606-8501, Japan

December 2010

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTINGAPPROXIMATION OF SCALE FUNCTIONSMASAHIKO EGAMI AND KAZUTOSHI YAMAZAKI
A BSTRACT. The optimal dividend problem by De Finetti (1957) has been recently generalized to the spectrally negative Lévy model where the implementation of optimal strategies draws upon the computation ofscale functions and their derivatives. This paper proposes a phase-type fitting approximation of the optimal
strategy. We consider spectrally negative Lévy processes with phase-type jumps as well as meromorphic
Lévy processes (Kuznetsov et al., 2010a), and use their scale functions to approximate the scale function
for a general spectrally negative Lévy process. We obtain analytically the convergence results and illustratenumerically the effectiveness of the approximation methods using examples with the spectrally negativeLévy process with i.i.d. Weibull-distributed jumps, the β-family and CGMY process.

Key words: De Finetti’s dividend problem; phase-type models; Meromorphic Lévy processes; Spectrally negativeLévy processes; Scale functions
JEL Classification: G22, D81, C61
Mathematics Subject Classification (2000) : Primary: 93E20 Secondary: 60G51

1. I NTRODUCTION
In the optimal dividend problem, an insurance company wants to maximize the cumulative amount of dividends
paid out to its beneficiaries until the time of ruin. The original problem by De Finetti (1957) considered a discretetime case by modeling the surplus by a random walk. It has been later extended to the continuous-time diffusion
model (Jeanblanc and Shiryaev, 1995; Asmussen and Taksar, 1997; Gerber and Shiu, 2004) and to the CramérLundberg model (Gerber, 1969; Azcue and Muler, 2005). It was recently generalized to the spectrally negativeLévy model by Avram et al. (2007) and Kyprianou and Palmowski (2007) where the surplus is a general Lévy
process with only negative jumps.
Date: November 22, 2010.
M. Egami is in part supported by Grant-in-Aid for Scientific Research (C) No. 20530340, Japan Society for the Promotionof Science. K. Yamazaki is supported by Grant-in-Aid for Young Scientists (B) No. 22710143, the Ministry of Education,
Culture, Sports, Science and Technology.
1

2

M. EGAMI AND K. YAMAZAKI

The implementation of the optimal strategies in the spectrally negative Lévy model draws upon the computationof the so-called scale function. Avram et al. (2007) first obtained the expected value under a barrier strategy interms of the scale function. Loeffen (2008) then showed that a barrier strategy is indeed optimal under a suitablecondition. The scale function also plays great roles in its extensions with transaction costs and additional terminal
values at ruin. See Loeffen (2009a,b) and Loeffen and Renaud (2010).
Despite these advances, a major obstacle still remains in putting in practice the above-mentioned results because
scale functions are in general known only up to their Laplace transforms. The implementation is even harder forthe optimal dividend problem because it requires the first derivative of the scale function. One can in principle
approximate the scale function and its derivative via the numerical Laplace inversion as discussed in Surya (2008).
However, the method is only heuristic and one cannot determine the accuracy of the approximation. Moreover, thederivative of the scale function tends to explode in the neighborhood of zero, and therefore the error is expected
to be large near zero. For these reasons, there is a clear need of new approximation methods that work for any
spectrally negative Lévy process.
As a new tool to overcome these problems, this paper proposes a phase-type fitting approach to approximate thescale functions as well as the solutions to the optimal dividend problem and its extensions. We obtain the scale
functions of spectrally negative Lévy processes with phase-type jumps and those in the M-class (meromorphic
Lévy processes), and apply these for a general spectrally negative Lévy model.
We first consider the class of spectrally negative Lévy processes with phase-type jumps. Consider a continuoustime Markov chain with some initial distribution and state space consisting of a single absorbing state and a finite
number of transient states. The phase-type distribution is the distribution of the time to absorption. It is known
that the class of phase-type distributions is dense in the class of all positive-valued distributions. We obtain thescale functions for these processes and show that they can approximate the scale function of a general spectrally
negative Lévy process arbitrarily closely.
The phase-type fitting approach has mainly three advantages. First, thanks to the smoothness and monotonicity
properties of scale functions as proved by, for example, Chan et al. (2009) and Loeffen (2008), the approximation
can be applied also to its derivative. Second, the Laplace transform of the phase-type distribution has an explicit
expression and hence can avoid the error caused while approximating the Laplace transform for a general jump
distribution. This type of errors tends to occur in other approximation methods such as Surya (2008). Third, thephase-type fitting approach enjoys a variety of fitting algorithms. See, for example, Asmussen (1996), Bladt et al.
(2003) and Feldmann and Whitt (1998). The fitting can be applied also to empirical data and this is another majoradvantage. In the first half of our numerical results, we use the results by Feldmann and Whitt (1998) and consider
the case with a Brownian motion plus a compound Poisson with Weibull-distributed jumps.

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS

3

We also consider meromorphic Lévy processes (Kuznetsov et al., 2010a), which generalizes a number ofrecently-discovered Lévy processes such as Lamperti-stable processes (Caballero and Chaumont, 2006; Chaumont et al., 2009), hypergeometric processes (Kyprianou, 2010; Kuznetsov et al., 2010b) and processes in the βand Θ-families (Kuznetsov, 2009a,b). The spectrally negative versions of these processes commonly have Lévy
measures in the form
(1.1)

ν(dz) ="∞
X#αi ηi e−ηi |z|1{z<0} dz,z ∈ R,

i=1

which can be seen as an extension of the hyperexponential density and as a “discrete version” of the completely
monotone density.
The reason we consider the Lévy measure above is that it can approximate efficiently any Lévy measure witha completely monotone density. Furthermore, if the Lévy measure is completely monotone, the barrier-strategy is
guaranteed to be optimal for the dividend problem (Loeffen, 2008). The class of Lévy measures with completely
monotone densities is rich. It enables us to model compound-Poisson-type jumps with long-tailed distributions
such as the Pareto, Weibull and gamma distributions; see Feldmann and Whitt (1998). It further allows us toconstruct many of recently-introduced Lévy processes such as variance gamma processes (Madan and Milne,
1991; Madan et al., 1998), CGMY processes (Carr et al., 2002), generalized hyperbolic processes (Eberlein et al.,1998) and normal inverse Gaussian processes (Barndorff-Nielsen, 1998). These processes can be approximatedefficiently by meromorphic Lévy processes.
Given a Lévy measure in the form (1.1), the corresponding scale function and its derivative can be expressed
explicitly as infinite sums of exponential functions. We first show that these can be approximated by finite sums
with some analytical error bounds, concluding that the error bounds for the solutions to the optimal dividend
problem can also be obtained. We then show numerically the effectiveness of the approximation procedure using
examples with β-processes (Kuznetsov, 2009a) and CGMY processes. We obtain bounds on scale functions andsolutions for the former and use them to approximate for the latter.
It should be emphasized here that the approximation procedure discussed in this paper can be applied outside
the class of optimal dividend problems. The spectrally negative Lévy model has been recently introduced widely
and the scale function plays a great role in characterizing the solutions. We refer the reader to Avram et al. (2004)
and Alili and Kyprianou (2005) for derivative pricing, Kyprianou and Surya (2007) for optimal capital structure,
Baurdoux and Kyprianou (2008, 2009) for stochastic games. For a comprehensive account, see Kyprianou (2006).
The rest of the paper is organized as follows. Section 2 summarizes the results on the classical dividend problem
and its extensions. In Section 3, we obtain the scale functions for spectrally negative Lévy processes with phasetype jumps and show that it can approximate the scale function of a general spectrally negative Lévy process

4

M. EGAMI AND K. YAMAZAKI

arbitrarily closely. We then obtain in Section 4 the scale functions for those in the M-class including its upperand lower bounds. We conclude this paper by giving numerical results in Section 5. All proofs are given in the
appendix.
2. O PTIMAL S OLUTIONS TO D IVIDEND P ROBLEM VIA S CALE F UNCTIONSThis section reviews the classical dividend problem and its extensions focusing on the spectrally negative Lévy
model. Let (Ω, F, P) be a probability space hosting a spectrally negative Lévy process X = {Xt ; t ≥ 0} with its
Laplace exponent
(2.1)



1
ψ(s) := log E esX1 = µ̂s + σ 2 s2 +2

Z 0(esz − 1 − sz1{z>−1} )ν(dz),

s∈C

−∞

R
where ν is a Lévy measure with the support (−∞, 0) that satisfies the integrability condition (−∞,0) (1∧z 2 )ν(dz) <∞. Moreover, let Px be the conditional probability under which X0 = x (also let P ≡ P0 ), and F := {Ft : t ≥ 0}
be the filtration generated by X. The process X is called the Lévy insurance risk process (or the risk process in
short), and models the surplus of an insurance company before dividends are deducted.
The classical dividend problem is a control problem where the cumulative amount of dividends prior to ruin is
maximized. A (dividend) strategy π := {Lπt ; t ≥ 0} is given by a nondecreasing, left-continuous and F-adapted
process starting at zero. Corresponding to every strategy π, the remaining amount of surplus after dividends are
deducted is given by U π = {Utπ : t ≥ 0} whereUtπ := Xt − Lπt ,t ≥ 0,and its ruin time is the first time it goes below zero:σ π = inf {t > 0 : Utπ < 0} .
A lump-sum payment must be smaller than the available fund and hence it is required that
(2.2)

Lπt+ − Lπt ≤ Utπ ,t < σ π a.s.

Let Π be the set of all admissible strategies satisfying (2.2). The problem concerns the expected sum of total
discounted dividends until ruinxZ σ πvπ (x) := E

−qt

e

dLπt


,

π ∈ Π,

0and wants to obtain an admissible strategy that maximizes it. Hence the classical dividend problem is written as(2.3)v(x) := sup vπ (x).
π∈ΠSOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS

5

2.1. Scale functions. The solutions to the dividend problems can be written in terms of the scale function. Here,
we describe the scale function and summarize its properties that will be used in this paper.
Fix q ≥ 0 and any spectrally negative Lévy process with its Laplace exponent ψ as defined in (2.1). The scale
function W (q) : R 7→ R is a function whose Laplace transform is given byZ ∞1e−sx W (q) (x)dx =(2.4),
s > ζqψ(s) − q0
whereζq := sup{s ≥ 0 : ψ(s) = q},(2.5)

q ≥ 0.

We assume W (q) (x) = 0 on (−∞, 0).
Let us define the first down- and up-crossing times, respectively, by(2.6)

τa := inf {t ≥ 0 : Xt < a}and

τb+ := inf {t ≥ 0 : Xt > b} ,for every 0 ≤ a ≤ x ≤ b. Then we have
(2.7)

i W (q) (x)h+Ex e−qτb 1{τ + <τ0 } =b
W (q) (b)andihW (q) (x)Ex e−qτ0 1{τ + >τ0 } = Z (q) (x) − Z (q) (b) (q)
b
W (b)

where
Z

(2.8)

(q)

Z x
(x) := 1 + q

W (q) (y)dy,x ∈ R.

0Here, we disregard the case when X is a negative subordinator (or decreasing a.s.).
We also consider a version of the scale function Wζq : R 7→ R that satisfies
W (q) (x) = eζq x Wζq (x),

(2.9)

x∈R

with its Laplace transform
Z ∞(2.10)

e−sx Wζq (x)dx =01,ψ(s + ζq ) − q

s > 0.

Suppose Pc , for any given c > 0, is the probability measure defined by the Esscher transformdPc= ecXt −ψ(c)t ,dP Ft

t ≥ 0;

see page 78 of Kyprianou (2006). Then Wζq under Pζq is analogous to W (0) under P. Furthermore, it is known
that Wζq is monotonically increasing and(2.11)Wζq (x) % (ψ 0 (ζq ))−1

as x → ∞,6

M. EGAMI AND K. YAMAZAKI

which also implies that the scale function W (q) increases exponentially in x;
W (q) (x) ∼(2.12)eζq x
ψ 0 (ζq )as x → ∞.

Due to the fact that Wζq does not explode for large x as opposed to W (q) , it is often convenient to deal with Wζq
and convert it to W (q) using (2.9), especially when numerical computations are involved; see Surya (2008).
Recall that a spectrally negative Lévy process has paths of bounded variation if and only if
Z(1 ∧ |x|)ν(dx) < ∞;σ = 0 and(−∞,0)

see, for example, Kyprianou (2006), Lemma 2.12. In this case, we can rewrite the Laplace exponent (2.1) by
Z 0
1 2 2
(esz − 1)ν(dz), s ∈ C,
ψ(s) = µs + σ s +
2
−∞
with
Z 0
µ := µ̂ −zν(dz).
−1

Regarding the smoothness of the scale function, it has been shown by Chan et al. (2009) that if a Lévy process
has a Gaussian component (σ > 0), we have W (q) ∈ C 2 (0, ∞). When it does not have a Gaussian component
and if its jump distribution has no atoms, we have W (q) ∈ C 1 (0, ∞). In particular, a stronger result holds forthe completely monotone jump case. Recall that a density function f is called completely monotone if all thederivatives exist and, for every n ≥ 1,(−1)n f (n) (x) ≥ 0,x ≥ 0,

where f (n) denotes the nth derivative of f .
Lemma 2.1 (Loeffen (2008)). If the Lévy measure has a completely monotone density, Wζq is again completely
monotone.
Finally, the behavior in the neighborhood of zero is given as follows. See Lemmas 4.3 and 4.4 of Kyprianou
and Surya (2007).
Lemma 2.2. For every q ≥ 0, we have

W (q) (0) =


 0,


unbounded variation  1 , bounded variation
µ



and

0

W (q) (0+) =


2


,

 σ2

σ>0







.
∞,
σ = 0 and ν(−∞, 0) = ∞






 q+ν(−∞,0) , compound Poissonµ2

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS

7

2.2. Solutions in terms of scale functions. In a spectrally negative Lévy model, the expected value under thebarrier strategy can be expressed in terms of the scale function as shown by Avram et al. (2007). A barrier strategy
at level a ≥ 0 is denoted by πa := {Lat ; t ≤ σa } where
Lat := sup(Xs − a) ∨ 0,t ≥ 0.

s≤t

We further let σa := inf {t > 0 : Utπa < 0} denote the corresponding ruin time.
Theorem 2.1 (Avram et al. (2007), (5.1)). For every a ≥ 0, we have
 u (x),x ≤ a,
a(2.13)vπa (x) = x − a + ua (a), x > a,whereua (x) :=(2.14)

W (q) (x).
W (q)0 (a)Furthermore, the barrier strategy attains optimality under a suitable condition.
Theorem 2.2 (Loeffen (2008)). Suppose a∗ satisfies00

W (q) (a) ≤ W (q) (b),a∗ ≤ a ≤ b.

Then the barrier strategy πa∗ is an optimal strategy.
In view of the above, an optimal barrier strategy exists, for example, when W 0 is convex, which holds wheneverLévy measure has a completely monotone density as in Lemma 2.1.
Extension with bail-out. A variant called the bail-out problem is discussed in Avram et al. (2007). Here thebeneficiary of the dividends must inject capital to keep the risk process from going below zero. A strategy is now
a pair π = Lπ , Rπ where Lπ is the cumulative amount of dividends as in the classical model and Rπ is aright-continuous process representing the cumulative amount of injected capital satisfying
Z ∞(2.15)
e−qt dRtπ < ∞, a.s.
0Assume that ϕ > 1 is the cost per unit injected capital, the problem is to maximize
Z ∞

Z ∞
x
−qt
π
−qt
π
v π (x) = E
e dLt − ϕ
e dRt
0

0

among all strategies that satisfy (2.2) and (2.15). In this model, the optimal strategy reduces to the double-barrier
strategy that regulates the risk process Vt := Xt − Lπt + Rtπ inside the interval [0, d∗ ] whereno0(2.16)d∗ := inf a > 0 : G(a) := (ϕZ (q) (a) − 1)W (q) (a) − ϕqW (q) (a)2 ≤ 0 .

8

M. EGAMI AND K. YAMAZAKI

The value function becomes v d∗ (x) with
 R
i

h
 ϕ x Z (q) (y)dy + ψ0 (0+) + Z (q) (x) 1−ϕZ (q) (d) , 0 ≤ x ≤ d,(q)q
0qW (d)(2.17)v d (x) = x − d + v (d),x > d.
dExtension with terminal values at ruin. Loeffen and Renaud (2010) recently considered the case with additional terminal costs at ruin; the objective function is
Z σπ

x
−qt
π
−qσ π
π
max E
e dLt + e
P (Uσπ )1{σπ <∞}π∈Π

0for some affine function P (y) := S + Ky. Its special case with constant terminal value was studied by Loeffen(2009a). With the assumption that the tail of the Lévy measure is log-convex, the optimal solution is either thebarrier strategy or the take-the-money-and-run strategy where the latter immediately pays out all the dividends andforces the ruin to occur immediately. For the former case, the optimal barrier level is given byb∗ := sup {b ≥ 0 : F (b) ≥ F (x) for all x ≥ 0}
where for every x ≥ 0
F (x) :=

(2.18)

1 − A(x),W (q)0 (x)andA(x) := K Z (q) (x) − ψ 0 (0+)W (q) (x) + SqW (q) (x).
The value function is given by veb∗ (x) withR
 S + x A(y)dy + 1−A(b)W (q) (x), 0 ≤ x ≤ b,00
W (q) (b)(2.19)veb (x) = x − b + v (b),x > b.
b
Extension with transaction costs. An extension allowing transaction costs is discussed in Loeffen (2009b)
where the objective is to maximize



Z σπ
X
Ex 
e−qt d Lπt −δ1{∆Lπ >0}  .
s

0

0≤s<t

Here δ > 0 is the unit transaction cost and the strategy is assumed to be defined by a pure jump process in the form
Lπt =

X

∆Lπs ,t ≥ 0.

0≤s<t

For this impulse control problem, the role of the barrier strategy in the classical model is now replaced by the socalled (c1 , c2 )-policy which is commonly known in inventory control. The (c1 , c2 )-policy brings the risk process

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS

9down to the level c1 whenever the risk process goes above the level c2 . Loeffen (2009b) showed that (c∗1 , c∗2 )-policy
is optimal if it satisfies c∗2 − c∗1 − δ ≥ 0, minimizes the functiong(c1 , c2 ) =0

W (q) (c2 ) − W (q) (c1 ),c2 − c1 − δ0and satisfies W (q) (a) ≤ W (q) (b) for every c∗2 ≤ a ≤ b. The value function has the form (2.13) where the optimal
threshold level is replaced with c∗2 .
3. S CALE FUNCTIONS FOR SPECTRALLY NEGATIVE L ÉVY PROCESSES WITH PHASE - TYPE JUMPSAs we have seen in the last section, the implementation of the optimal strategies in the optimal dividend problem
draws upon the computation of scale functions. This section obtains the scale function of the spectrally negativeLévy process with phase-type jumps and shows that it can be used to approximate the scale function of any
spectrally negative Lévy process.
3.1. Spectrally negative Lévy processes with phase-type jumps. Consider a continuous-time Markov chain
Y = {Yt ; t ≥ 0} with finite state space {1, . . . , m} ∪ {∆} where 1, . . . , m are transient and ∆ is absorbing. Its
initial distribution is given by a simplex α = [α1 , . . . , αm ] such that αi = P {Y0 = i} for every i = 1, . . . , m. The
intensity matrix Q is partitioned into the m transient states and the absorbing state ∆, and is given by


T t
.
Q := 
0 0Here T is an m × m-matrix called the phase-type generator, and t = −T 1 where 1 = [1, . . . , 1]0 . A distribution iscalled phase-type with representation (m, α, T ) if it is the distribution of the absorption time to ∆ in the Markov
chain described above. It is known that T is non-singular and thus invertible; see Asmussen (1996). Its distribution
and density functions are given, respectively, by
F (z) = 1 − αeT z 1

andf (z) = αeT z t,

z ≥ 0.

Let X = {Xt ; t ≥ 0} be a spectrally negative Lévy process of the form
(3.1)

Xt − X0 = µt + σBt −

Nt
X

Zn ,

0 ≤ t < ∞,n=1

for some µ ∈ R and σ ≥ 0. Here B = {Bt ; t ≥ 0} is a standard Brownian motion, N = {Nt ; t ≥ 0} is a Poissonprocess with arrival rate λ, and Z = {Zn ; n = 1, 2, . . .} is an i.i.d. sequence of phase-type distributed randomvariables with representation (m, α, T ). These processes are assumed independent. Its Laplace exponent is then
(3.2)


1
ψ(s) = µs + σ 2 s2 + λ α(sI − T )−1 t − 1 ,2

10

M. EGAMI AND K. YAMAZAKI

which is analytic for every s ∈ C except for the eigenvalues of T .
Disregarding the case when X is a negative subordinator, we consider the following two cases:
Case 1: when σ > 0 (i.e. X has unbounded variation),
Case 2: when σ = 0 and µ > 0 (i.e. X is a compound Poisson process).
P t
Notice, in Case 2, that we can write Xt = Ut − N
n=1 Zn where Ut = x + µt is a (positive) subordinator. Thisimplies that down-crossing of a threshold can occur only by jumps; see, for example, Chapter III of Bertoin (1996).
On the other hand, in Case 1, down-crossing can occur also by creeping downward (by the diffusion components).
Due to this difference, the form of the scale function differs as we shall see.
Fix q > 0. Consider the Cramér-Lundberg equation
ψ(s) = q,(3.3)

and define the set of (the absolute values of) negative roots and the set of poles:
Iq := {i : ψ(−ξi,q ) = q and R(ξi,q ) > 0} ,qJq := j := 0 and R(ηj ) > 0 .
q − ψ(−ηj )The elements in Iq and Jq may not be distinct, and, in this case, we take each as many times as its multiplicity. By
Lemma 1 of Asmussen et al. (2004), we have
 |J | + 1, for Case 1,q|Iq | =
 |Jq |,
for Case 2.
In particular, if the representation is minimal (see Asmussen et al. (2004)), we have |Jq | = m.
Let κq be an independent exponential random variable with parameter q and denote the running maximum andminimum, respectively, by
X t = sup Xs

and

0≤s≤t

X t = inf Xs ,0≤s≤tt ≥ 0.

−The Wiener-Hopf factorization states that q/(q − ψ(s)) = ϕ+q (s)ϕq (s) for every s ∈ C such that R(s) = 0, with

the Wiener-Hopf factors
(3.4)

h
i
ϕ−
q (s) := E exp(sX κq )

and



ϕ+
q (s) := E exp(sX κq )that are analytic for s with R(s) > 0 and R(s) < 0, respectively. By Lemma 1 of Asmussen et al. (2004), wehave, for every s such that R(s) > 0,
Q
(3.5)

ϕ−
q (s) =Q

j∈Jq (s + ηj )Q

j∈Jq ηj

Q

i∈Iq ξi,qi∈Iq (s + ξi,q ),SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 11

from which we can obtain the distribution of X κq by the Laplace inverse via partial fraction expansion.
As in Remark 4 of Asmussen et al. (2004), let n denote the number of different roots in Iq and mi denote themultiplicity of a root ξi,q for i = 1, . . . , n. Then we have
mi
n X
n
o X
(ξi,q x)k−1 −ξi,q x(k)
P −X κq ∈ dx =Ai,q ξi,qedx,(k − 1)!

(3.6)

x>0

i=1 k=1

where(k)

Ai,q :=mi∂ mi −k ϕ−1
q (s)(s + ξi,q )k(mi − k)! ∂smi −kξi,q.
s=−ξi,qNotice that this can be simplified significantly when all the roots in Iq are distinct.
3.2. Scale functions for spectrally negative Lévy processes with phase-type jumps. Here we obtain the scale
function. We focus on the case q > 0 because the scale function when q = 0 (and X drifts to infinity) can bederived by using W (0) (x) = Px {X ∞ ≥ 0} /ψ 0 (0) and the ruin probability (19) of Asmussen et al. (2004) bytaking q → 0. Kyprianou and Palmowski (2007) briefly stated the scale function when q = 0 and all the roots in
Iq are distinct.
Before obtaining the scale function, we shall first represent the positive root ζq (2.5) in terms of the negativeroots {ξi,q ; i ∈ Iq }. Let us define
%q :=

(3.7)

n
X

(1)

Ai,q ξi,q ,q > 0,

i=1

and by Lemma 2.2
(3.8)




2
,
σ2
(q)
(q)0
θ := −ζq W (0) + W (0+) = − ζq + q+λ ,µµ2

for Case 1 .
for Case 2 Lemma 3.1. For every q > 0, we have
ζq
θ
= .
q
%q
We now obtain the version of the scale function Wζq (·). In the lemma below, Wζq (0) = W (q) (0) is either 0 or
1
µ depending on if it is Case 1 or Case 2; see Lemma 2.2.

Lemma 3.2. For every q > 0, we haveWζq (x) − Wζq (0) =n mi
ζq X Xq

i=1 k=1

(k)

Ai,q



ξi,q
ζq + ξi,qk
1 − e−(ζq +ξi,q )xk−1
X((ζq + ξi,q )x)jj=0

Lemma 3.2 together with (2.12) and Lemmas 2.2 and 3.1 shows the following.

j!


,

x ≥ 0.

12

M. EGAMI AND K. YAMAZAKI

Proposition 3.1. For every q > 0 and x ≥ 0, we have the following.
(1) For Case 1, we have
W (q) (x) =

2

mi
n X
X

σ 2 %q
i=1

(k)

Ai,q

k=1



ξi,q
ζq + ξi,qkeζq x − e−ξi,q x

k−1
X((ζq + ξi,q )x)j
j=0

j!


.

(2) For Case 2, we have
W (q) (x) =

1
%q


−

ζq
q+λ
+
µ
µ2

X
mi
n X
i=1 k=1

(k)

Ai,q



ξi,q
ζq + ξi,qkeζq x − e−ξi,q x

k−1
X((ζq + ξi,q )x)j
j=0

j!


 + 1 eζq x .
µRecall that the solution to the dividend problem requires the derivative. The scale functions obtained above areinfinitely differentiable. In particular, the first derivative becomes



k
mi
n X
k−2
j
k−1
X
X
ξi,q
((ζq + ξi,q )x) 20(k)ζq eζq x + ξi,q e−ξi,q x ((ζq + ξi,q )x)
Ai,q− ζq e−ξi,q x
W (q) (x) = 2
σ %q
ζq + ξi,q(k − 1)!
j!
j=0

i=1 k=1

for Case 1 and
W

(q)0

1
(x) =
%q


 n mikζq
ξi,qq + λ XX
(k)
− +
Ai,q
µ
µ2
ζq + ξi,q
i=1 k=1


k−2
j
k−1
X
((ζq + ξi,q )x)  1 ζq x((ζq + ξi,q )x)− ζq e−ξi,q x+ ζq e× ζq eζq x + ξi,q e−ξi,q x(k − 1)!
j!
µj=0

for Case 2.
When all the roots in Iq are distinct, the scale functions above can be simplified and have nice properties as
discussed in the following corollary.
Corollary 3.1. If all the roots in Iq are distinct, we have the followings.
(1) The scale function can be simplified to

h
n
i
ξi,q
2 X (1)(q)
W (x) = 2
Ai,qeζq x − e−ξi,q x ,σ %q
ζq + ξi,q
i=1
 n


h
i 1
ζq
ξi,q
q + λ X (1)
1
(q)
ζq x
−ξi,q x
W (x) =− +Ae−e
+ eζq x ,
i,q
%q
µ
µ2
ζq + ξi,qµi=1

for Case 1 and Case 2, respectively.
0(2) W (q) is convex.
(3) Wζ0q is completely monotone.

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 13

This guarantees the optimality of the barrier-strategy in view of Theorem 2.2; the optimal strategy can be00

obtained by finding a unique a∗ such that W (q) (a∗ ) = 0.
Example 3.1 (Hyperexponential Case). As an important example where all the roots in Iq are distinct, we consider
the case where Z has a hyperexponential distribution with density function
f (z) =mX

αi ηi e−ηi z ,z ≥ 0,

i=1

for some 0 < η1 < · · · < ηm < ∞. Its Laplace exponent (2.1) is then
m

X
s
1
.
αi
ψ(s) = µs + σ 2 s2 − λ2ηi + s(3.9)

i=1

Notice in this case that −η1 , . . . , −ηm are the poles of the Laplace exponent. Furthermore, all the roots in Iq are
distinct and satisfy the following interlacing condition for every q > 0:(1) when σ > 0, there are m + 1 roots −ξ1,q , . . . , −ξm+1,q such that
0 < ξ1,q < η1 < ξ2,q < · · · < ηm < ξm+1,q < ∞;(2) when σ = 0 and µ > 0, there are m roots −ξ1,q , . . . , −ξm,q such that
0 < ξ1,q < η1 < ξ2,q < · · · < ξm,q < ηm < ∞.
The class of hyperexponential distributions is important as it is dense in the class of all positive-valued distributions
with completely monotone densities.
3.3. Approximation of the scale function of a general spectrally negative Lévy process. The scale function
obtained in Proposition 3.1 can be used to approximate the scale function of a general spectrally negative Lévy
process. By Proposition 1 of Asmussen et al. (2004), there exists, for any spectrally negative Lévy process X, asequence of spectrally negative Lévy processes with phase-type jumps X (n) converging to X in D[0, ∞). This(n)is equivalent to saying that X1→ X1 in distribution by Jacod and Shiryaev (2003), Corollary VII 3.6; see also(q)

Pistorius (2006). Suppose ψn (ψ), ζq,n (ζq ) and Wn /Wζq ,n (W (q) /Wζq ) are the Laplace exponent, the positiveroot (2.5) and the scale function of X (n) (X), respectively. Because these processes are spectrally negative andψ is continuous, we have, by the continuity theorem, ψn (β + ζq,n ) → ψ(β + ζq ) for every β > 0. Now inview of (2.10), the convergence of the scale function holds by the continuity theorem; see Feller (1971), Theorem
R
R
R
R
(q)
2a, XIII.1. More precisely, we have I Wζq ,n (y)dy → I Wζq (y)dy and I Wn (y)dy → I W (q) (y)dy for any
interval I.

14

M. EGAMI AND K. YAMAZAKI

The smoothness and monotonicity properties of the scale function can be additionally used to obtain stronger
results. The scale functions in Proposition 3.1 are in C ∞ (0, ∞). In addition, when all the roots of Iq are different,
its first derivative Wζ0q is completely monotone as discussed in Corollary 3.1.
If the target scale function is in C 1 (0, ∞) (which holds whenever the jump distribution has no atoms), notingthat Wζq (x) ≤ (ψ 0 (ζq ))−1 for every x and hence e−βx Wζq (x) vanishes in the limit for any β > 0, we have by(2.10)
Z ∞0e−βx Wζ0q (x)dx =Because Wζ0q (x) is nonnegative and F (x) :=β− Wζq (0),
ψ(β + ζq ) − q

β > 0.

Rx00−1 − W (0)) is a probability distribution,ζq
0 Wζq (y)dy/((ψ (ζq ))n↑∞Wζq ,n (x) −−−→ Wζq (x) andn↑∞Wn(q) (x) −−−→ W (q) (x),x ≥ 0.

Furthermore, suppose that it is in C 2 (0, ∞) (which holds, for example, when σ > 0 by Chan et al. (2009)),x↑∞

Wζ0q (0+) < ∞ (i.e., σ > 0 or Π(−∞, 0) < ∞) and Wζ00q (x) ≤ 0 for every x ≥ 0, because Wζ0q (x) −−−→ 0, weRxhave F (x) := (Wζ0q (0+))−1 0 |Wζ00q (y)|dy is a probability distribution and
Z ∞
e
0

−βx

F (dx) = (Wζ0q (0+))−1




β2
0
−
+ βWζq (0) + Wζq (0+) ,ψ(β + ζq ) − q

β > 0.

0

Therefore, noting that W (q) (x) = ζq W (q) (x) + eζq x Wζ0q (x) and assuming that the convergent sequence Wζq ,n (x)has the same property, we can obtain by the continuity theoremn↑∞Wζ0q ,n (x) −−−→ Wζ0q (x) and

0

n↑∞

0

Wn(q) (x) −−−→ W (q) (x),x ≥ 0.

The negativity of Wζ00q holds, for example, for the completely monotone jump case because Wζ0q is completely
monotone by Loeffen (2008). We can also choose the sequence Wζ0q ,n completely monotone in view of Corollary
3.1 because approximation can be done via hyperexponential distributions. In fact, it also means that Wζq is
C ∞ (0, ∞) and the convergence of higher derivatives can be pursued. Even for a general jump distribution, thenegativity of Wζ00q is a reasonable assumption in view of the numerical plots given by Surya (2008).
This phase-type fitting approach complements the approach by Surya (2008) where scale functions are approximated by numerical Laplace inversion. A major disadvantage of using this inversion method is the fact that it
requires the exact value of the right-hand side of (2.4). However, the Laplace transform of a jump distribution does
not in general have an explicit closed-form expression. Surya (2008)’s approach, therefore, contains two types oferrors: 1) the approximation error caused while computing ψ and 2) the error caused while inverting the Laplacetransform. On the other hand, the phase-type fitting approach only contains the phase-type fitting error thanks tothe closed-form Laplace transform of the phase-type distribution.

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 15

The phase-type fitting approach enjoys a variety of fitting algorithms typically developed in queueing analysis.
Well-known examples are the moment-matching approach (e.g. MEFIT and MEDA) and the maximum-likelihood
approach (e.g. MLAPH and EMPHT), and a thorough study of pros and cons of each fitting techniques has been
conducted in, for example, Horváth and Telek (2000) and Lang and Arthur (1996). The fitting can be applied also
to empirical data and this is another major advantage over the Laplace inversion approach.
4. S CALE FUNCTIONS FOR M EROMORPHIC L ÉVY PROCESSESIn this section, we consider another class of spectrally negative Lévy processes called meromorphic Lévy processes. We obtain their scale functions and use these as approximation tools for a general spectrally negative Lévy
process with a completely monotone Lévy measure. Similarly to the approach applied in the last section, we obtainthe scale function using its Wiener-Hopf factorization. It has a form expressed as an infinite sum of exponentialfunctions which can be bounded efficiently by finite sums.
4.1. Meromorphic Lévy processes. The following is due to Kuznetsov et al. (2010a), Definition 1.
Definition 4.1 (spectrally negative meromorphic Lévy process). A spectrally negative Lévy process X is said tobelong to the M-class if the following conditions hold.
(1) The Laplace exponent ψ(s) (2.1) has a countable set of real negative poles.
(2) For every q ≥ 0, the Cramér-Lundberg equation (3.3) has a countable set of real negative roots.
(3) Let {ηk ; k ≥ 1} and {ξk,q ; k ≥ 1}, respectively, be the sets of the absolute values of the poles and thenegative roots of (3.3) for fixed q ≥ 0. Then it satisfies the following interlacing conditions:· · · < −ηk < −ξk,q < · · · < −η2 < −ξ2,q < −η1 < −ξ1,q < 0.
(4) There exists α > 21 such that ηk ∼ ck α as k → ∞.
(5) The Wiener-Hopf factor (3.4) is expressed as convergent infinite products
(4.1)

ϕ−
q (s) =∞
Y(s + ηk )
k=1

ηk

ξk,q
.
(s + ξk,q )The M-class complements the class of Lévy processes with phase-type jumps described in the previous section
because it also contains those of infinite activity. As noted by Corollary 3 of Kuznetsov et al. (2010a), the property(3) in Definition 4.1 is equivalent to the condition that the Lévy measure has the form (1.1). This can be seen as anextension to the hyperexponential case as described in Example 3.1.
We consider the M-class alternatively to the “hyperexponential fitting”. When a Lévy measure is completely
monotone, approximation via hyperexponential distributions is in principle possible. However, as in, for example,

16

M. EGAMI AND K. YAMAZAKI

Asmussen et al. (2007), special care is needed for the infinitesimal jumps, and one needs to approximate separately
the process with Lévy measure ν(−ε, 0) for small ε > 0. Fitting via the M-class is more tractable in the sense that
this procedure is not necessary. Although property (4) requires one to choose η’s in a certain way, the approximation for the Lévy process with a completely monotone density is still effective by choosing the value of c and α
sufficiently small. For more details, see Kuznetsov et al. (2010a).
The Wiener-Hopf factor (4.1) is again a rational function as in (3.5) for the phase-type case. Therefore, this canbe inverted again by partial fraction decomposition, and we have∞o X
nP −X κq ∈ dx =Ak,q ξk,q e−ξk,q x dx,

(4.2)

x>0

k=1

where
s + ξk,q −ϕq (s)=
Ak,q :=
ξk,qs=−ξk,q



ξk,q
1−
ηk

Y

ξ

1 − ηk,qiξk,qi6=k 1 − ξi,q,k ≥ 1.

Notice by the interlacing condition that Ak,q > 0 for every k ≥ 1.
4.2. Scale functions for meromorphic Lévy processes. We now obtain the scale function for the M-class. We
omit the proof because it is similar to the phase-type case; see Appendix A.1.
Lemma 4.1. For every q > 0, we haveWζq (x) − Wζq (0) =

∞
X

i
h
Ck,q 1 − e−(ζq +ξk,q )x ,

x≥0

k=1

where
Ck,q :=(4.3)

ζq ξk,q Ak,q,q ζq + ξk,qk ≥ 1.

By (2.11) and Lemma 4.1, we have, by taking the limit,
κq :=

(4.4)

∞
X

Ck,q = (ψ 0 (ζq ))−1 − Wζq (0) < ∞.

k=1

The scale function can be therefore obtained by Lemma 4.1 and (4.4).
Proposition 4.1. For every q > 0, we have
(4.5)

W (q) (x) =

∞
X
i=1

∞
h
i
X
Ci,q eζq x − e−ξi,q x + Wζq (0)eζq x = (ψ 0 (ζq ))−1 eζq x −Ci,q e−ξi,q x ,i=1

x ≥ 0.

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 17

By straightforward differentiation, we have, for every q > 0 and x ≥ 0,0

W (q) (x) = (ψ 0 (ζq ))−1 ζq eζq x +∞XCi,q ξi,q e−ξi,q x ,

i=1

W

(q)00

0

−1

(x) = (ψ (ζq ))2 ζq x(ζq ) e−∞X

Ci,q (ξi,q )2 e−ξi,q x .

i=10

Remark 4.1. From the derivatives above, it can be easily verified that W (q) (·) is indeed convex. In view of
00

Theorem 2.1, the optimal solution to the classical dividend problem (2.3) is the unique point x such that W (q) (x)vanishes or(ψ 0 (ζq ))−1 (ζq )2 eζq x −∞XCi,q (ξi,q )2 e−ξi,q x = 0.

i=1

As an extension to Lemma 3.1, we have the following.
Lemma 4.2.

(1) The following two statements are equivalent:(a) σ = 0 and ν(−∞, 0) = ∞,P∞(b)
k=1 Ak,q ξk,q = ∞.
(2) Suppose σ > 0 or ν(−∞, 0) < ∞. Then, for every q > 0, we have
∞
X
ζq
=θ
Ak,q ξk,q
q

(4.6)

!−1

k=1

where



2
,
0
σ2
θ := −ζq W (q) (0) + W (q) (0+) =ζ − q + q+ν(−∞,0) ,
µ
µ2

(4.7)


when σ > 0 .

when σ = 0 

4.3. Approximation of the scale functions via finite sum. The scale function obtained in Proposition 4.1 is an
infinite sum of exponential functions and in reality its exact value cannot be computed. Here, we obtain bounds
0for W (q) (·), W (q) (·) and Z (q) (·) in terms of finite sums.
For every m ≥ 1, let(m)
Ak,q := 1{k≤m}



ξk,q
1−
ηk

(m)

ξ



1 − ηk,qiY

ξk,q1≤i≤m,i6=k 1 − ξi,q

and

(m)
ζq ξk,q Ak,q(m)
Ck,q :=,q ζq + ξk,q(m)

By the interlacing condition, Ak,q and Ck,q are all positive and, for every k ≥ 1,(m)

Ak,q ↑ Ak,q

and

(m)

Ck,q ↑ Ck,q as m → ∞.

k ≥ 1.

18

M. EGAMI AND K. YAMAZAKI

Now we define candidates for upper and lower bounds of Wζq respectively by(m)W ζq (x) := (ψ 0 (ζq ))−1 −

m
X

(m)

Ci,q e−(ζq +ξi,q )x ,

i=1

h
i
(m)
(m)
W ζq (x) := W ζq (x) − δm e−ζq x + e−(ζq +ξm+1,q )x ,for every m ≥ 1 and x ≥ 0, whereδm := κq −

m
X

(m)

Ci,q > 0,i=1

which vanishes in the limit as m → ∞ by (4.4). As candidates for upper and lower bounds of W (q) , we also define
W

(q,m)

(m)

(x) := eζq x W ζq (x) and(m)

W (q,m) (x) := eζq x W ζq (x),x ≥ 0.

The following proposition shows that the scale functions are bounded and approximated by these functions.
Proposition 4.2. For every m ≥ 1 and x ≥ 0, we have
(m)

(m)

W ζq (x) ≤ Wζq (x) ≤ W ζq (x)

(4.8)

and

W (q,m) (x) ≤ W (q) (x) ≤ W

(q,m)

(x).

Furthermore, we have
(m)

m↑∞

(m)

m↑∞

W ζq (x) −−−→ Wζq (x) andW ζq (x) −−−→ Wζq (x) and

W

(q,m)

m↑∞

(x) −−−→ W (q) (x),m↑∞W (q,m) (x) −−−→ W (q) (x),uniformly on x ∈ [0, ∞).
By straightforward calculation, we can bound Z (q) in (2.8). Let, for every m ≥ 1,
Z x
Z x
(q,m)
(m,q)
Z
(x) := 1 + q
W(y)dy and Z (q,m) (x) := 1 + q
W (m,q) (y)dy,0

Then by Proposition 4.2, we have Z (q,m) (x) ≤ Z (q,m) (x) ≤ Z
0≤Z

(q,m)

(x) − Z (q,m) (x) = qx ≥ 0.

0

Z x

W

(q,m)

(q,m)

(x) and(y) − W (q,m) (y) dy

0


Z xh
i
−ξm+1,q y= qδm1+edy = qδm x +01
ξm+1,q(1 − e−ξm+1,q x) .

We therefore have the following.
Corollary 4.1 (Bounds on Z (q) ). We have Z
for every x ≥ 0.

(q,m)(x) → Z (q) (x) and Z (q,m) (x) → Z (q) (x) as m → ∞ pointwiseSOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 19

We now obtain bounds for the derivative. Define, for every x > 0,
w

(m)

0

−1

(x) := (ψ (ζq ))

ζq e

ζq x

+

m
X

(m)

Ci,q ξi,q e−ξi,q x ,i=1

w(m) (x) := w(m) (x) +max (ξk,q e−ξk,q x ) + max (ξk,q e−ξk,q x ) δm .

1≤k≤m

k≥m+1

Here notice thatmax (ξk,q e−ξk,q x ) =

k≥m+1


 1,

ξm+1,q ≤ x1 ,x ξm+1,q e−ξm+1,q x , ξm+1,q > 1 .
x

Proposition 4.3. For every m ≥ 1, we have0w(m) (x) ≤ W (q) (x) ≤ w(m) (x),x ≥ 0.

00Furthermore, we have w(m) (x) → W (q) (x) and w(m) (x) → W (q) (x) uniformly on x ≥ x0 for any x0 > 0.
A stronger result holds when σ > 0 or ν(−∞, 0) < ∞. Recall in this case that θ < ∞ by Lemma 4.2 (2) andhence we can definemm := θ −ζq X
(m)ξi,q Ai,q > 0,qm ≥ 1,

i=1

which vanishes in the limit as m → ∞ by Lemma 4.2 (2).
Corollary 4.2. When σ > 0 or ν(−∞, 0) < ∞, we have0e(m) (x),w(m) (x) ≤ W (q) (x) ≤ w(m) (x) ∧ w

x>0

where
w
e(m) (x) := w(m) (x) + max (ξk,q e−ξk,q x )δm + e−ξm+1,q x m .
1≤k≤m

The bounds obtained above on the derivative of the scale function can be used to obtain bounds on the value0functions and optimal barriers. For simplicity, let w and w be the lower and upper bounds for W (q) obtained above
and letw∗ = min w(x)x≥0and w∗ = min w(x).
x≥0Clearly, w and w are convex because w(m) , w(m) and we(m) are for every fixed m ≥ 1.

20

M. EGAMI AND K. YAMAZAKI
0Remark 4.2 (classical dividend problem). Because w and w are convex, we have w∗ ≤ W (q) (a∗ ) ≤ w∗ and we
have
W

(q,m)

(q,m)

(x) − δm (1 + e−ξm,q x )W(x)≤ ua∗ (x) ≤.
w∗
w∗

Furthermore, the optimal barrier must lie in the following interval:a∗ ∈ {x ≥ 0 : w(x) ≤ w∗ } .
Remark 4.3 (bail-out problem). For the bail-out problem, recall that the optimal barrier is the minimum value ofa such that G(a) becomes negative (see (2.16)). We can obtain its bounds simply by
G(a) := [ϕZ (q,m) (a) − 1]w(a) − ϕqW

(q,m)

(a)2

and

G(a) := [ϕZ(q,m)(a) − 1]w(a) − ϕqW (q,m) (a)2 .
∗We can therefore obtain the bounds on the optimal barrier level d∗ ; we have d∗ ≤ d∗ ≤ d whered∗ := inf {a > 0 : G(a) ≤ 0}

and


∗
d := inf a > 0 : G(a) ≤ 0 .

5. N UMERICAL E XAMPLES
We conclude this paper by illustrating numerically the effectiveness of the phase-type fitting approximation fora general spectrally negative Lévy process. First, we use the classical hyperexponential fitting algorithm for a completely monotone density function by Feldmann and Whitt (1998) and approximate, as an example, the scale function and solutions for the case with a (Brownian motion plus) compound Poisson with Weibull-distributed jumps.
Second, we consider, as an example of the meromorphic Lévy process, the β-family introduced by Kuznetsov(2009a) and extend the results to the spectrally negative version of the CGMY process.
5.1. Brownian motion plus compound Poisson process with Weibull-distributed jumps. As noted earlier,
any spectrally negative Lévy process with a completely monotone Lévy measure can be approximated arbitrarilyclosely by fitting hyperexponential distributions. Here, we use the fitted data computed by Feldmann and Whitt(1998) to approximate the scale function when it is a Brownian motion plus a compound Poisson process with i.i.d.
Weibull-distributed jumps. Recall that the Weibull distribution with parameters c and a (Weibull(c,a)) is give bycF (t) = 1 − e−(t/a) ,t ≥ 0.

If c < 1, it has long-tails, or eδt (1 − F (t)) → ∞ as t → ∞ for any δ > 0, and has a completely monotone density.
Feldmann and Whitt (1998) constructed a recursive algorithm to approximate completely monotone densities
in terms of hyperexponential densities. We use their results and compute the scale functions of spectrally negativeLévy processes with Weibull-distributed jumps.

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 21

i

αi

ηi

i

αi

ηi

1 0.029931 676.178 4 0.476233 0.76100
2 0.093283 38.7090 5 0.068340 0.24800
3 0.332195 4.27400 6 0.000018 0.09700
TABLE 1. Parameters of the hyperexponential distribution fitted to Weibull(0.6,0.665) (takenfrom Tables 3 of Feldmann and Whitt (1998)).

Table 1 shows the parameters of the hyperexponential distribution fitted to Weibull(0.6,0.665) when m = 6. As
can be seen in Figure 4 of Feldmann and Whitt (1998), this fitting is very accurate. We consider the Lévy process
in the form (3.1) where Z is hyperexponential specified in Table 1 as an approximation to Weibull(0.6,0.665). We
use various values of σ with the common values of λ = 1, µ = 0.1 and q = 0.03. The roots ξ·,q ’s and ζq are
calculated via the bisection method with error bound 1.0E − 10.
00Figure 1 shows the scale function W (q) and its derivative W (q) . The optimal barrier levels that minimize W (q)are given by a∗ = 0.05, 0.481 and 0.643 for the cases σ = 0, 0.2 and 0.4, respectively. The results are consistent
with Lemma 2.2; with the existence of a diffusion component, the scale function is forced to converge to 0 as x
goes to 0. Using these barrier levels, optimal value functions for the classical dividend problem can be computedby Theorem 2.1. Figure 2 shows the value functions vπa∗ . Notice that they are monotonically decreasing in σ.

45

55
σ=0
σ = 0.2
σ = 0.4

35

45

30

40

25

35

20
15

30
25

10

20

5

15

0

10

-5

0

0.2

0.4

0.6
x

scale function W (q)

0.8

σ=0
σ = 0.2
σ = 0.4

50

derivative

scale function

40

1

5

0

0.2

0.4

0.6

0.8

x

0

derivative W (q)F IGURE 1. Scale functions and their derivatives for the case with Weibull-distributed jumps.

122

M. EGAMI AND K. YAMAZAKI

1.6
σ=0
σ = 0.2
σ = 0.4

1.4
1.2

value function

1
0.8
0.6
0.4
0.2
0
-0.2

0

0.2

0.4

0.6

0.8

1

x

F IGURE 2. Value function vπa∗ of the classical dividend problem.
We now consider the extensions described in the end of Section 2. Here we use the same parameters as in the
results above. Figure 3 shows the results on the bail-out problem (Avram et al., 2007) when ϕ = 1.3. It plots G
in (2.16) as well as the value function v d∗ in (2.17). Here the optimal barrier level is obtained by computing theunique level d∗ that satisfies G(d∗ ) = 0. Figure 4 shows the results on the extension with terminal values at ruin
(Loeffen and Renaud, 2010) with the plots of F in (2.18) and the value function veb∗ in (2.19). We consider the casewith constant terminal value (i) S = −1 and K = 0 and (ii) S = 1 and K = 0. The maximizer of F becomes thebarrier level b∗ . Figure 5 shows the value functions vc∗2 on the extension with transaction costs (Loeffen, 2009b)
when (i) δ = 0.5 and (ii) δ = 0.1. In order to obtain the optimal impulse contol (c∗1 , c∗2 ), we use the technique
discussed in Section 4 of Loeffen (2009b). Unlike the other results, the value functions are no longer monotone
in σ unless δ is sufficiently small. However, as δ decreases to zero, c∗1 converges to c∗2 and the value function
converges to that of the classical model as shown in Figure 2.
5.2. Numerical results on the β-class and CGMY process. We now consider, as an example of meromorphic
Lévy processes, the β-class introduced by Kuznetsov (2009a). The following definition is due to Kuznetsov(2009a), Definition 4.
Definition 5.1. A spectrally negative Lévy process is said to be in the β-class if its Lévy measure is in the form
(5.1)

ν(dx) = c

eαβx
1
dx,
(1 − eβx )λ {x<0}x ∈ R,

for some α > 0, β > 0, c ≥ 0 and λ ∈ (0, 3). It is equivalent to saying that its Laplace exponent is
1 2 2 cz
ψ(z) = µ̂z + σ z +B(α + , 1 − λ) − B(α, 1 − λ)2ββSOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 23

20-448.5σ=0
σ = 0.2
σ = 0.410

σ=0
σ = 0.2
σ = 0.4-449

0-449.5Gvalue function

-10

-20

-450

-450.5
-30
-451

-40

-50

0

0.5

1

1.5

2

-451.5

0

0.5

x

1

1.5

2

x

G(x)

value function v d∗

F IGURE 3. Bail-out problem: G in (2.16) and the value function v d∗ when ϕ = 1.3.
Optimal barriers are 0.38, 0.775 and 1.495 for σ = 0, 0.2 and 0.4, respectively.
where B is the beta function B(x, y) := Γ(x)Γ(y)/Γ(x + y).
The special case σ = 0 and β = 1 reduces to the class of Lamperti-stable processes, which are obtained by
the Lamperti transformation (Lamperti (1972)) from the stable processes conditioned to stay positive; see Bertoin
and Yor (2001) and Caballero et al. (2008) and references therein. For the scale function of a related process, see
Kyprianou and Rivero (2008).
It can be also seen that this is a “discrete-version” of the (spectrally negative) CGMY process, whose Lévy
measure is given by
ν(dx) = c

(5.2)

eαx
1
dx,
|x|λ {x<0}x ∈ R.

Indeed, if we set c = ecβ λ and α = αeβ −1 in (5.1), we have
(5.3)

c

eαβx
eαex
β↓0
1{x<0} −−→ ec λ 1{x<0} ,βx
λ(1 − e )|x|x ∈ R.

See Asmussen et al. (2007) for approximation of (double-sided) CGMY processes using hyperexponential distributions.
We shall use the results in Section 4 to obtain the bounds on the scale functions and the solutions to the classicaloptimal dividend problem. Figure 6 shows the approximation results when q = 0.03, σ = 0.2, µ̂ = 0.1, λ = 1.5,
α = 3, β = 1 and c = 0.1 in (5.1). We plot, for m = 15 and m = 150, the upper and lower bounds on the scale

24

M. EGAMI AND K. YAMAZAKI

0.14

2
σ=0
σ = 0.2
σ = 0.4

0.12

σ=0
σ = 0.2
σ = 0.4

1.5

1

0.5

value function

F

0.1

0.08

0

0.06
-0.5
0.04

0.02

-1

0

0.5

1

1.5

2

-1.5

0

0.5

x

F (x) when (i) S = −1, K = 0

2

3.5
σ=0
σ = 0.2
σ = 0.4

0.08

σ=0
σ = 0.2
σ = 0.4

3

2.5
value function

0.06

F

1.5

value function veb∗ when (i) S = −1, K = 0

0.1

0.04

2

0.02

1.5

0

1

-0.02

1
x

0

0.5

1

1.5

2

x

F (x) when (ii) S = 1, K = 0

0.5

0

0.5

1

1.5

2

x

value function veb∗ when (ii) S = 1, K = 0

F IGURE 4. Extension with terminal values at ruin: F in (2.18) and the value function veb∗
when (i) S = −1, K = 0 and (ii) S = 1, K = 0. The optimal barriers are (i) 0.0628,
0.0793 and 0.1384 and (ii) 0.0317, 0.0383 and 0.0955, respectively, for σ = 0, 0.2 and
0.4.

function, its derivative and the function ua∗ defined in (2.14). As shown in the previous section, the difference
between the upper and lower bounds indeed converges to zero.

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 25

1.8

2.5
σ=0
σ = 0.2
σ = 0.4

1.6

2

1.4
1.2

1.5

1

value function

value function

σ=0
σ = 0.2
σ = 0.4

0.8
0.6

1

0.5

0.4
0.2

0

0
-0.2

0

0.5

1

1.5

2

-0.5

0

0.5

1

x

1.5

2

x

value function vc∗2 when (i) β = 0.5

value function vc∗2 when (ii) β = 0.1

F IGURE 5. Extension with transaction costs: value function when (i) β = 0.5 and (ii)β = 0.1. The optimal impulse controls (c∗1 , c∗2 ) are (i) (0, 1.173), (0.069, 1.527) and(0, 1.885) and (ii) (0, 0.05), (0.222, 0.481) and (0.197, 0.643) for σ = 0, 0.2 and 0.4,
respectively.
We now take β in (5.3) to zero and see how the approximation for the CGMY process works. Here we set α̃ = 3
and c̃ = 0.1 and use the same values as the above for the other parameters. Figure 7 shows the upper and lower
bounds of scale function and its derivative for various values of β. Figure 8 shows the mean value of the upperand lower bounds on the function u. Here we can indeed observe the convergence as β → 0. This implies that it
effectively approximates the scale function and the solution for the CGMY case.
A PPENDIX A. P ROOFS
A.1. Proof of Lemmas 3.1 and 3.2. By (3.6), it is easy to verify that
Z ∞
mi
n X

 X
(ξi,q y)k−1 −ξi,q y(k)x −qτaedy,
E e1{τa <∞} =Ai,q ξi,qx−a (k − 1)!

0 ≤ a < x,i=1 k=1

and hence, because W (q) ∈ C 1 (0, ∞),
n

(A.1)

(A.2)

m

i
 XX
(ξi,q (x − a))k−1 −ξi,q (x−a)∂ x  −qτa
(k)
E e
e
,
1{τa <∞} =Ai,q ξi,q∂a(k − 1)!

i=1 k=1nX∂ x  −qτ0
E e1{τ0 <∞}
=−
∂x
x=0+

(1)

Ai,q ξi,q .

i=1

0 ≤ a < x,

26

M. EGAMI AND K. YAMAZAKI

250

250
lower boundupper bound

200

200

150

150
scale function

scale function

lower boundupper bound

100

100

50

50

0

0

-50

0

0.5

1

1.5

-50

2

0

0.5

1

x

Bounds on W (q) (x) when m = 15

100lower boundupper bound

90
80

80

70

70

60

60

50
40

50
40

30

30

20

20

10

10

0

0.2

0.4

0.6

0.8

1

1.2

lower boundupper bound90

derivative

derivative

2

Bounds on W (q) (x) when m = 150

100

0

1.5

x

0

1.4

0

0.2

0.4

0.6

x

0.8

1

1.2

1.4

x

0

0

Bounds on W (q) (x) when m = 15

Bounds on W (q) (x) when m = 150

76lower boundupper bound6

lower boundupper bound
5

5

4

4

u

u

3
3

2
2
1

1

0

0

-1

0

0.5

1

1.5

x

Bounds on u(x) when m = 15

2

-1

0

0.5

1

1.5

2

x

Bounds on u(x) when m = 150

F IGURE 6. Approximation of the scale function and its derivative for the β-class.

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 27

400

1800
2
1
0.5
0.2
0.1

350

1400
1200

250

derivative (upper value)

scale function (upper value)

300

2
1
0.5
0.2
0.1

1600

200

150

1000

100

600
400

50

0

800

200

0

0.5

1

1.5

0

2

0

0.5

x

1

1.5

2

x

upper bounds for the scale functions

upper bounds for the derivatives

400

700
2
1
0.5
0.2
0.1

350
300

2
1
0.5
0.2
0.1

600

derivative (lower value)

scale function (lower value)

500
250
200
150
100

400

300

200

50
100

0
-50

0

0.5

1

1.5

2

0

0

0.5

x

1

1.5

2

x

lower bounds for the scale functions

lower bounds for the derivatives

F IGURE 7. Convergence of scale functions to the CGMY model.

In fact, different representations of (A.1) and (A.2) can be pursued. By Theorem 8.1 of Kyprianou (2006) and(2.9),

q
E e−qτa 1{τa <∞} = Z (q) (x − a) − W (q) (x − a) = 1 + q
ζq
x

Z x−a
0

W (q) (y)dy −q ζq (x−a)
eWζq (x − a)ζq28

M. EGAMI AND K. YAMAZAKI

10
2
1
0.5
0.2
0.1

9
8

u (mean value)

7
6
5
4
3
2
1
0

0

0.5

1

1.5

2

x

F IGURE 8. Convergence of u to the CGMY model. Here we use M = 150.
for every 0 ≤ a < x. Its derivative with respect to a becomes

(A.3)


q
∂ x  −qτaE e1{τa <∞} = −qW (q) (x − a) + qeζq (x−a) Wζq (x − a) + eζq (x−a) Wζ0q (x − a)∂aζqq ζq (x−a) 0
Wζq (x − a).
= eζq

In particular, when a = 0, the derivative with respect to x and its limit as x → 0 are(A.4)
i x↓0+i∂ x  −qτ0
q h
q h
q
0
0
E e
−ζq W (q) (x) + W (q) (x) −−−→ −−ζq W (q) (0) + W (q) (0+) = − θ.
1{τ0 <∞} = −
∂x
ζq
ζq
ζq
By matching (A.2) and (A.4), Lemma 3.1 is immediate.
For the proof of Lemma 3.2, by matching (A.1) and (A.3) and using Lemma 3.1, we havenmi
(ξi,q y)k−1 −(ζq +ξi,q )yθ XX(k)Ai,q ξi,qe%q(k − 1)!
i=1 k=1

k
mi
n X
X
ξi,q
((ζq + ξi,q )y)k−1 −(ζq +ξi,q )yθ(k)=Ai,q (ζq + ξi,q )e,%qζq + ξi,q(k − 1)!

Wζ0q (y) =i=1 k=1

y ≥ 0.

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 29

Integrating the above and changing variables, we havek
Z (ζq +ξi,q )x
n mi
ξi,q
θ XX
1
(k)
z k−1 e−z dz,Wζq (x) − Wζq (0) =Ai,q%qζq + ξi,q(k − 1)! 0

x ≥ 0.

i=1 k=1

Lemma 3.2 is now immediate because the integral part is a lower incomplete gamma function.
A.2. Proof of Proposition 4.2. Notice, for every m ≥ 1, that
0≤
(A.5)
0≤

m 
X

(m)
Ck,q − Ck,q

k=1
∞
X



≤ κq −mX(m)

Ck,q = δm ,k=1

Ck,q = κq −

k=m+1

m
X

m
X

Ck,q ≤ κq −

k=1

(m)

Ck,q = δm ,k=1

and hence by (A.5)
m 
∞

X
X
(m)
(m)
−(ζq +ξi,q )x
Ck,q − Ck,q e+
Ck,q e−(ζq +ξk,q )x0 ≤ W ζq (x) − Wζq (x) =
k=1
k=m+1
−ζq x

≤e

m 
X

(m)
Ck,q − Ck,q+e−(ζq +ξm+1,q )x

k=1

∞
X

i
h
Ck,q ≤ δm e−ζq x + e−(ζq +ξm+1,q )x .

k=m+1

Therefore we have the bounds for Wζq in (4.8). The bounds for W (q) are immediate by multiplying eζq x . Finally,
the convergence results hold because δm → 0 and e−ζq x +e−(ζq +ξm+1,q )x and 1+e−ξm+1,q x are bounded uniformly
in x ≥ 0.
(m)

A.3. Proof of Proposition 4.3. The lower bound is immediate by the fact that 0 ≤ Ci,q ≤ Ci,q . For every x ≥ 0,
we have by (A.5)
0

W (q) (x) = (ψ 0 (ζq ))−1 ζq eζq x +∞XCi,q ξi,q e−ξi,q x

i=1

=w

(m)

(x) +

m
X

(m)
(Ci,q − Ci,q )ξi,q e−ξi,q x +

i=1

≤w
≤w

(m)

−ξk,q x(x) + max (ξk,q e(x) +max (ξk,q e

1≤k≤m

Ci,q ξi,q e−ξi,q x

i=m+1

1≤k≤m

(m)

∞
X

)

m
X

(m)
(Ci,q − Ci,q ) +i=1−ξk,q x−ξk,q xmax (ξk,q ek≥m+1−ξk,q x

) + max (ξk,q e
k≥m+1

)

∞
X

Ci,q

i=m+1


) δm = w(m) (x),which shows the first claim. For the second claim, notice that for the given x0 > 0,
0 ≤ w(m) (x) − w(m) (x) ≤ max (ξk,q e−ξk,q x0 ) + max (ξk,q e−ξk,q x0 ) δm
1≤k≤m

k≥m+1

uniformly on [x0 , ∞). Because δm vanishes as m → 0 and supλ≥0 λe−λx < ∞, the convergence is immediate.

30

M. EGAMI AND K. YAMAZAKI

A.4. Proof of Corollary 4.2. We haveW(q)0(x) = w

(m)

(x) +

m
X

(m)
(Ci,q − Ci,q )ξi,q e−ξi,q x +

i=1

≤w

(m)

∞
X

Ci,q ξi,q e−ξi,q x

i=m+1(x) + max (ξk,q e−ξk,q x

1≤k≤m

)

m
X

ζq
(m)
(Ci,q − Ci,q ) + e−ξm+1,q x
q

i=1

∞
X

Ai,q ξi,qi=m+1≤ w(m) (x) + max (ξk,q e−ξk,q x )δm + m e−ξm+1,q x
1≤k≤mwhere the last inequality holds because
∞
m
m
ζq X
ζq X
ζq X
(m)
ξi,q Ai,q = θ −ξi,q Ai,q ≤ θ −ξi,q Ai,q = m ,
q
q
q
k=m+1

k=1

m ≥ 1.

k=1

This together with Proposition 4.3 shows the claim.

R EFERENCES
Alili, L. and Kyprianou, A. E. (2005). Some remarks on first passage of Lévy processes, the American put andpasting principles. Ann. Appl. Probab., 15(3):2062–2080.
Asmussen, S. (1996). Fitting phase-type distributions via the em algorithm. Scand. J. Statist, 23:419–441.
Asmussen, S., Avram, F., and Pistorius, M. R. (2004). Russian and American put options under exponential
phase-type Lévy models. Stochastic Process. Appl., 109(1):79–111.
Asmussen, S., Madan, D., and Pistorius, M. (2007). Pricing equity default swaps under an approximation to the
CGMY levy model. Journal of Computational Finance, 11(2):79–93.
Asmussen, S. and Taksar, M. (1997). Controlled diffusion models for optimal dividend pay-out. Insurance: Math.
Econom., 20:1–15.
Avram, F., Kyprianou, A. E., and Pistorius, M. R. (2004). Exit problems for spectrally negative Lévy processes
and applications to (Canadized) Russian options. Ann. Appl. Probab., 14(1):215–238.
Avram, F., Palmowski, Z., and Pistorius, M. R. (2007). On the optimal dividend problem for a spectrally negativeLévy process. Ann. Appl. Probab., 17(1):156–180.
Azcue, P. and Muler, N. (2005). Optimal reinsurance and dividend distribution policies in the Cramér-Lundberg
model. Math. Finance, 15(2):261–308.
Barndorff-Nielsen, O. E. (1998). Processes of normal inverse Gaussian type. Finance Stoch., 2(1):41–68.
Baurdoux, E. and Kyprianou, A. E. (2008). The McKean stochastic game driven by a spectrally negative Lévy
process. Electron. J. Probab., 13:no. 8, 173–197.

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 31

Baurdoux, E. and Kyprianou, A. E. (2009). The shepp-shiryaev stochastic game driven by a spectrally negative
lévy process. Theory of Probability and Its Applications, 53.
Bertoin, J. (1996). Lévy processes, volume 121 of Cambridge Tracts in Mathematics. Cambridge University Press,
Cambridge.
Bertoin, J. and Yor, M. (2001). On subordinators, self-similar Markov processes and some factorizations of the
exponential variable. Electron. Comm. Probab., 6:95–106 (electronic).
Bladt, M., Gonzalez, A., and Lauritzen, S. L. (2003). The estimation of phase-type related functionals using
Markov chain Monte Carlo methods. Scand. Actuar. J., (4):280–300.
Caballero, M., Pardo, J., and Pérez, J. (2008). On lamperti stable processes. Preprint.
Caballero, M. E. and Chaumont, L. (2006). Conditioned stable Lévy processes and the Lamperti representation. J.
Appl. Probab., 43(4):967–983.
Carr, P., Hélyette, G., Madan, D. B., and Yor, M. (2002). The structure of asset returns: an emperical investigation.
Journal of Business, 75:305–332.
Chan, T., Kyprianou, A., and Savov, M. (2009). Smoothness of scale functions for spectrally negative Lévy
processes. Probability Theory and Related Fields.
Chaumont, L., Kyprianou, A. E., and Pardo, J. C. (2009). Some explicit identities associated with positive selfsimilar Markov processes. Stochastic Process. Appl., 119(3):980–1000.
De Finetti, B. (1957). Su un’impostazion alternativa dell teoria collecttiva del rischio. Trans. XVth Internat. Congr.
Actuaries, 2:433–443.
Eberlein, E., Keller, U., and Prause, K. (1998). New insights into smile, mispricing and value at risk: the hyperbolicmodel. Journal of Business, 71:371–405.
Feldmann, A. and Whitt, W. (1998). Fitting mixtures of exponentials to long-tail distributions to analyze network
performance models. Performance Evaluation, (31):245–279.
Feller, W. (1971). An introduction to probability theory and its applications. Vol. II. Second edition. John Wiley& Sons Inc., New York.
Gerber, H. U. (1969). Entscheidungskriterien fr den zusammengesetzten poisson-prozess. Schweiz. Verein. Versicherungsmath. Mitt., 69:185–228.
Gerber, H. U. and Shiu, E. S. W. (2004). Optimal dividends: analysis with Brownian motion. N. Am. Actuar. J.,
8(1):1–20.
Horváth, A. and Telek, M. (2000). Approximating heavy tailed behaviour with phase type distributions. Advancesin Algorithmic Methods for Stochastic Models.

32

M. EGAMI AND K. YAMAZAKI

Jacod, J. and Shiryaev, A. N. (2003). Limit theorems for stochastic processes, volume 288 of Grundlehren derMathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]. Springer-Verlag, Berlin,second edition.
Jeanblanc, M. and Shiryaev, A. N. (1995).

Optimization of the flow of dividends.

Uspekhi Mat. Nauk,50(2(302)):25–46.
Kuznetsov, A. (2009a). Wiener-hopf factorization and distribution of extrema for a family of Lévy processes. Ann.
Appl. Probab.
Kuznetsov, A. (2009b). Wiener-hopf factorization for a faimily of Lévy processes related to theta functions.
preprint.
Kuznetsov, A., Kyprianou, A., and Pardo, J. (2010a). Meromorphic Lévy processes and their fluctuation identities.
Preprint.
Kuznetsov, A., Kyprianou, A., Pardo, J., and Schaik, K. v. (2010b). A wiener-hopf monte-carlo simulation technique for Lévy processes. Preprint.
Kyprianou, A. E. (2006). Introductory lectures on fluctuations of Lévy processes with applications. Universitext.
Springer-Verlag, Berlin.
Kyprianou, A. E. (2010). Exact and asymptotic n-tuple lwas at first and last passage. Ann. Appl. Probab.,20(2):522–564.
Kyprianou, A. E. and Palmowski, Z. (2007). Distributional study of de Finetti’s dividend problem for a generalLévy insurance risk process. J. Appl. Probab., 44(2):428–443.
Kyprianou, A. E. and Rivero, V. (2008). Special, conjugate and complete scale functions for spectrally negativeLévy processes. Electron. J. Probab., 13:no. 57, 1672–1701.
Kyprianou, A. E. and Surya, B. A. (2007). Principles of smooth and continuous fit in the determination of endogenous bankruptcy levels. Finance Stoch., 11(1):131–152.
Lamperti, J. (1972). Semi-stable Markov processes. I. Z. Wahrscheinlichkeitstheorie und Verw. Gebiete, 22:205–
225.
Lang, A. and Arthur, J. (1996). Parameter Approximation for Phase-Type Distributions. Matrix-analytic Methods
in Stochastic Models (Lecture Notes in Pure and Applied Mathematics). CRC Press.
Loeffen, R. (2009a). An optimal dividends problem with a terminal value for spectrally negative Lévy processes
with a completely monotone jump density. Journal of Applied Probability, 46(1):85–98.
Loeffen, R. L. (2008). On optimality of the barrier strategy in de Finetti’s dividend problem for spectrally negativeLévy processes. Ann. Appl. Probab., 18(5):1669–1680.

SOLVING OPTIMAL DIVIDEND PROBLEMS VIA PHASE-TYPE FITTING APPROXIMATION OF SCALE FUNCTIONS 33

Loeffen, R. L. (2009b). An optimal dividends problem with transaction costs for spectrally negative Lévy processes. Insurance Math. Econom., 45(1):41–48.
Loeffen, R. L. and Renaud, J.-F. (2010). De Finetti’s optimal dividends problem with an affine penalty function at
ruin. Insurance Math. Econom., 46(1):98–108.
Madan, D., P.P., C., and E.C., C. (1998). The variance gamma processes and option pricing. European Finance
Review, 2:79–105.
Madan, D. B. and Milne, F. (1991). Option pricing with vg martingale components. Mathematical Finance,
1(4):39–55.
Pistorius, M. (2006). On maxima and ladder processes for a dense class of Lévy process. J. Appl. Probab.,43(1):208–220.
Surya, B. A. (2008). Evaluating scale functions of spectrally negative Lévy processes. J. Appl. Probab., 45(1):135–
149.
(M. Egami) G RADUATE S CHOOL OF E CONOMICS , K YOTO U NIVERSITY, S AKYO -K U , K YOTO , 606-8501, JAPANE-mail address: egami@econ.kyoto-u.ac.jp(K. Yamazaki) C ENTER FOR THE S TUDY OF F INANCE AND I NSURANCE , O SAKA U NIVERSITY, 1-3 M ACHIKANEYAMA CHO , T OYONAKA C ITY, O SAKA 560-8531, JAPANE-mail address: k-yamazaki@sigmath.es.osaka-u.ac.jp