# The bilateral generalized inverse Gaussian process with applications to
financial modeling

Gaetano Agazzotti  Société Générale & École des Mines de Nancy
<EMAIL_ADDRESS>and  Jean-Philippe Aguilar  Corresponding
Author, Société Générale<EMAIL_ADDRESS>
(Date: This version: July 15, 2024)

###### Abstract.

We introduce and document a class of probability distributions, called
bilateral generalized inverse Gaussian (BGIG) distributions, that are obtained
by convolution of two generalized inverse Gaussian distributions supported by
the positive and negative semi-axis. We prove several results regarding their
analyticity, shapes and asymptotics, and we introduce the associated Lévy
processes as well as their main properties. We study the behaviour of these
processes under change of measure, their simulations and the structure of
their sample paths, and we introduce a stock market model constructed by means
of exponential BGIG processes. Based on real market data, we show that this
model is easy to calibrate thanks notably to idiosyncratic properties of BGIG
distributions, and that it is well suited to Monte Carlo and Fourier option
pricing.

###### Key words and phrases:

Bilateral generalized inverse Gaussian distribution, Bilateral generalized
inverse Gaussian process, Jaeger integral, Lévy process, Esscher transform,
Calibration, Option pricing

###### 2020 Mathematics Subject Classification:

60G51, 60G46, 60E07, 60E10, 91B28

##  1\. Introduction

The rich interplay between mathematics, physics and finance has given birth to
a vast collection of stochastic processes whose study remains an active field
of research across academia and industry. In mathematical finance in
particular, two classes of processes appear to be of peculiar importance:
first, the tempered stable (TS) processes, that generalize the so-called
truncated Lévy flights first introduced in [43, 32] in the fully symmetric
case and then extended to single sided and asymmetric processes (notably via
KoBoL distributions, see [13, 14]), and to even more general configurations in
[35]. TS processes are particularly important because they recover (or are
limiting cases) of many popular processes previously introduced in finance to
model asset returns, such as the Carr-Geman-Madan-Yor (CGMY) process
introduced in [15] and sometimes called classical tempered sable process, or
the bilateral Gamma (see [36, 37, 34]) and Variance Gamma (see [41])
processes. Being Lévy processes, they allow to extend the classical Black-
Scholes-Merton framework of [11] by capturing the heavy tails that are
ubiquitous in financial data, while preserving the stationarity of returns.

The second important class is arguably the class of generalized hyperbolic
(GH) processes; they were introduced in [6] to formalize mathematically the
transport of sand by wind, following the initial observations by Bagnold in
the 1940s regarding the formation of dunes in the Lybian desert (paving the
way to what is now known as Aeolian processes). GH processes are typically
constructed as mean-variance mixtures of normal distributions, the mixing
distribution being the generalized inverse Gaussian (GIG) distribution (see
[19, 10] for GH distributions, and [19, 8, 30] for GIG distributions), and,
interestingly, have also found applications in financial modeling and risk
management, as discussed for instance in [21]. They recover other processes
that are well known in quantitative finance, such as the normal inverse
Gaussian process (see [7]) but also the Variance Gamma processes, which
demonstrates that the intersection of TS and GH classes is actually nonempty.

In this paper, we would like to follow a somewhat different approach and,
instead of considering GIG distributions only from a mixing distribution
perspective, to construct a bilateral distribution by convolution of two
independent GIG distributions defined on the positive and negative semi-axis
(note that a similar path had already been taken notably by Küchler and Tappe
when introducing bilateral Gamma processes in [36], but the generating
variables were Gamma distributed in this case). The resulting distributions we
obtain are double sided, and we name them bilateral generalized inverse
Gaussian (BGIG) distributions. They have six parameters, like TS distributions
(GH have five) and, indeed, have a nonempty intersection with the TS class;
but BGIG distributions also form a very a rich set on their own, giving birth
to a wide range of random variables that are not TS distributed. Moreover they
are infinitely divisible, allowing us to construct associated Lévy processes;
even if BGIG processes are not convolution-closed (because GIG processes
themselves are not, see e.g. [47]), they have a tractable characteristic
function, which makes them ideal for financial applications.

The present article has therefore multiple purposes. First, we wish to define
and document BGIG distributions as well as their immediate properties in
section 2 and to study their particular and limiting cases in section 3.
Thanks to recent results obtained in [5, 24, 27] on a specific family of
integrals (namely Jaeger integrals, originally introduced in [29] and that
happen to be involved in the Lévy measure of BGIG distributions), we will be
able to prove several results regarding the smoothness, shapes and asymptotics
of BGIG distributions in section 4. Then, we will define the BGIG Lévy process
and provide a detailed study of its behavior, paths and simulations in section
5, allowing us to properly define a stock market model based on the
exponential BGIG process in section 6. As an application, we will calibrate
this model over real market data and price European style options (via Monte
Carlo and Fourier methods) in section 7. Last, section 8 will be dedicated to
conclusions.

##  2\. Bilateral generalized inverse Gaussian distributions

We start our study by recalling some properties of GIG distributions, before
introducing BGIG distributions as well as some of their first properties
(representation of densities, cumulants).

###  2.1. Notations

In all of the following, we will consider a probability space
$(\Omega,\mathcal{F},\mathbb{P})$ and, unless otherwise mentioned,
expectations will be taken under $\mathbb{P}$, that is
$\mathbb{E}[.]=\mathbb{E}^{\mathbb{P}}[.]$. Given a random variable
$X:\Omega\rightarrow\mathbb{R}$, we define its characteristic function
(Fourier transform) as

(1) $\forall
u\in{\mathbb{R}},\quad\Phi_{X}(u)=\mathbb{E}\left[e^{\mathrm{i}uX}\right]=\int_{\mathbb{R}}e^{\mathrm{i}ux}f_{X}(x)\mathrm{d}x,$

where the second equality holds if we assume that $X$ possesses a
$\mathbb{P}$-density $f_{X}$.

###  2.2. Generalized inverse Gaussian distributions

Let us recall some essential facts about the generalized inverse Gaussian
(GIG) distribution; details and proofs can be found in [10, 9, 19].

####  2.2.1. Density

The GIG distribution is a three parameter probability distribution whose
density function is defined for all $x\in\mathbb{R}$ by :

(2)
$f_{\mathrm{GIG_{+}}}(x)=\frac{(a/b)^{p/2}}{2K_{p}\left(\sqrt{ab}\right)}x^{p-1}e^{-\left(\frac{ax+b/x}{2}\right)}\mathbbm{1}_{\\{x>0\\}}(x),\quad
a,b>0,\quad p\in\mathbb{R},$

where $K_{p}(.)$ denotes the modified Bessel function of the second kind (see
[1]). We will denote this distribution as $\mathrm{GIG}_{+}(a,b,p)$. It can
easily be translated to the negative semi-axis to define the “negative-sided”
GIG distribution, that we denote by $\mathrm{GIG}_{-}(a,b,p)$, and whose
density is:

(3)
$f_{\mathrm{GIG_{-}}}(x)=\frac{(a/b)^{p/2}}{2K_{p}\left(\sqrt{ab}\right)}|x|^{p-1}e^{-\left(\frac{a|x|+b/|x|}{2}\right)}\mathbbm{1}_{\\{x<0\\}}(x)$

under the same conditions on the parameters. Let us now provide some details
(Characteristic function, Lévy measure) for the $\mathrm{GIG}_{+}(a,b,p)$
distribution only, the extension to the $\mathrm{GIG}_{-}(a,b,p)$ case being
straightforward.

####  2.2.2. Characteristic function

The characteristic function of a random variable
$X\sim\allowbreak\mathrm{GIG}_{+}(a,b,$
$p)$ is known under closed form as

(4) $\forall
u\in{\mathbb{R}},\quad\Phi_{\mathrm{GIG}_{+}}(u)=\left(\frac{a}{a-\mathrm{i}2u}\right)^{p/2}\frac{K_{p}\left(\sqrt{b(a-\mathrm{i}2u)}\right)}{K_{p}\left(\sqrt{ab}\right)}.$

####  2.2.3. Lévy measure

Let us introduce the Hankel function of the first kind
$H_{\lambda}(z)=J_{\lambda}(z)+\mathrm{i}Y_{\lambda}(z)$ where $J_{\lambda}$
(resp. $Y_{\lambda}$) denotes the Bessel function of the first (resp. second)
kind (see [1] or any monograph on special functions), as well as the Jaeger
integral (see [29, 40]):

(5) $\forall
x\in(0,+\infty),\quad\mathscr{I}(x,p):=\frac{2}{\pi^{2}}\int_{0}^{\infty}\frac{e^{-xy^{2}}}{y\left|H_{|p|}(y)\right|^{2}}\mathrm{d}y.$

Then, defining the Lévy measure as

(6) $\forall
x\in{\mathbb{R}},\quad\pi_{\mathrm{GIG}_{+}}(\mathrm{d}x)=\frac{e^{-ax/2}}{x}\left(\int_{0}^{\infty}\frac{e^{-xy}}{\pi^{2}y\left|H_{|p|}(\sqrt{2by})\right|^{2}}\mathrm{d}y+\max(p,0)\right)\mathbbm{1}_{(0,\infty)}(x)\mathrm{d}x$

leads to the Lévy-Khintchine representation for the characteristic function
(4):

(7) $\forall
u\in{\mathbb{R}},\quad\Phi_{\mathrm{GIG}_{+}}(u)=\exp\left(\int_{\mathbb{R}}\left(e^{\mathrm{i}ux}-1\right)\pi_{\mathrm{GIG}_{+}}(\mathrm{d}x)\right).$

As noted in [27], the change of measure $\tilde{y}=\sqrt{2by}$ allows to
obtain an equivalent representation for the Lévy measure in terms of the
Jaeger integral:

(8)
$\pi_{\mathrm{GIG}_{+}}(\mathrm{d}x)=\frac{e^{-ax/2}}{x}\left(\mathscr{I}(x/2b,p)+\max(p,0)\right)\mathbbm{1}_{(0,\infty)}(x)\mathrm{d}x.$

The following lemma regarding asymptotic behaviours of the Jaeger integral
will be useful in many occurrences throughout the paper.

######  Lemma 2.1.

Let $p\in[0,+\infty)$, we have the following properties:

  1. (1)

Asymptotic behaviour in $0$:

(9) $\quad\mathscr{I}(x,p)\underset{x\to
0}{=}\frac{1}{2\pi^{1/2}\sqrt{x}}+\frac{1-2p}{4}+\mathcal{O}\left(\sqrt{x}\right),$

  2. (2)

Asymptotic behaviour in $+\infty$:

(10)
$\mathscr{I}(x,p)\underset{x\to\infty}{=}\frac{x^{-|p|}}{2^{2|p|}\Gamma(|p|)}+o\left(x^{-|p|}\right).$

###### Proof.

Lets us distinguish two cases:

  * •

for $p=0$, [24, Theorem D] and [24, Theorem I] prove respectively the
asymptotic (9) and (10),

  * •

for $p>0$, (9) and (10) are given by [24, Theorem I].

∎

###  2.3. Bilateral generalized inverse Gaussian distributions

We define the bilateral generalized inverse Gaussian (BGIG) distribution as
the convolution

(11)
$\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-}):=\mathrm{GIG}_{+}(a_{+},b_{+},p_{+})\star\mathrm{GIG}_{-}(a_{-},b_{-},p_{-}),$

where $\star$ stands for the convolution operator and
$\mathrm{GIG}_{+}(a_{+},b_{+},p_{+})$ (resp.
$\mathrm{GIG}_{-}(a_{-},b_{-},p_{-})$) is the $\mathrm{GIG}$ distribution
defined on the positive (resp. negative) semi-axis. The parameters satisfy
$(a_{+},b_{+},a_{-},b_{-})\in(0,+\infty)^{4}$ and
$(p_{+},p_{-})\in{\mathbb{R}}^{2}$. An example of a such distribution can be
found in Figure 1.

######  Remark 1.

Note that by construction, if $X_{+}\sim\mathrm{GIG}_{+}(a_{+},b_{+},p_{+})$
and $X_{-}\sim\mathrm{GIG}_{+}(a_{-},b_{-},p_{-})$, the difference is BGIG
distributed, i.e.
$X_{+}-X_{-}\sim\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$.

Let us give an immediate property of BGIG distributions (that will allow us to
construct the associated BGIG process in section 5):

######  Proposition 2.1.

BGIG distributions are infinitely divisible.

###### Proof.

GIG distributions are infinitely divisible (see [10]). From definition 11, it
follows that BGIG distributions are infinitely divisible as well. ∎

####  2.3.1. Characteristic function

An immediate consequence of definition 11 and of (4) is that the
characteristic function of a random variable
$X\sim\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$ can be expressed as:

(12) $\forall u\in{\mathbb{R}},~{}\Phi(u)=\\\
\left(\frac{a_{+}}{a_{+}-\mathrm{i}2u}\right)^{p_{+}/2}\frac{K_{p_{+}}\left(\sqrt{b_{+}(a_{+}-\mathrm{i}2u)}\right)}{K_{p_{+}}\left(\sqrt{a_{+}b_{+}}\right)}\times\left(\frac{a_{-}}{a_{-}+\mathrm{i}2u}\right)^{p_{-}/2}\frac{K_{p_{-}}\left(\sqrt{b_{-}(a_{-}+\mathrm{i}2u)}\right)}{K_{p_{-}}\left(\sqrt{a_{-}b_{-}}\right)}.$

####  2.3.2. Lévy measure

From definition (11) and from (6), we obtain the Lévy measure of the BGIG
distribution:

(13) $\forall
x\in{\mathbb{R}},\quad\pi(\mathrm{d}x)=\left(\frac{e^{-a_{+}x/2}}{x}\left(\mathscr{I}(x/2b_{+},p_{+})+\max(p_{+},0)\right)\mathbbm{1}_{(0,\infty)}(x)\right.\\\
\left.+\frac{e^{-a_{-}|x|/2}}{|x|}\left(\mathscr{I}(|x|/2b_{-},p_{-})+\max(p_{-},0)\right)\mathbbm{1}_{(-\infty,0)}(x)\right)\mathrm{d}x.$

The BGIG characteristic function (12) can therefore be written in the Lévy-
Khintchine form:

(14) $\forall
u\in{\mathbb{R}},\quad\Phi(u)=\exp\left(\int_{\mathbb{R}}\left(e^{\mathrm{i}ux}-1\right)\pi(\mathrm{d}x)\right).$

####  2.3.3. Cumulants

From definition 11, the cumulants $\kappa_{n}$ of the
$\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$ distribution can be
expressed as the sum of the cumulants of $\mathrm{GIG}_{+}(a_{+},b_{+},p_{+})$
and $\mathrm{GIG}_{-}(a_{-},b_{-},\\\ p_{-})$. From the characteristic
function (12), it is immediate to see that the cumulants of
$\mathrm{GIG}_{-}(a_{-},b_{-},p_{-})$ are equal to the cumulants of
$\mathrm{GIG}_{+}(a_{-},b_{-},p_{-})$ for even order and opposite otherwise.
If we denote $\kappa_{n}^{+}$ (resp. $\kappa_{n}^{-}$) the cumulants of
$\mathrm{GIG}_{+}(a_{+},b_{+},p_{+})$ (resp. $\mathrm{GIG}_{+}(a_{-},b_{-},\\\
p_{-})$), we can therefore write:

(15) $\kappa_{n}=\kappa_{n}^{+}+(-1)^{n}\kappa_{n}^{-}.$

The expression of $\kappa_{n}^{+}$ and $\kappa_{n}^{-}$ can be found up to
order four in [30, Eq. 2.25]; we therefore have:

(16) $\kappa_{n}^{\pm}=W_{n}(\omega_{\pm})\eta_{\pm}^{n}$

where $\omega_{\pm}:=\sqrt{a_{\pm}b_{\pm}}$, $\eta:=\sqrt{a_{\pm}/b_{\pm}}$
and:

(17) $\left\\{\begin{aligned} W_{1}(\omega_{\pm})&=R(\omega_{\pm}),\\\
W_{2}(\omega_{\pm})&=-R(\omega_{\pm})^{2}+\frac{2(p_{\pm}+1)}{\omega_{\pm}}R(\omega_{\pm})+1,\\\
W_{3}(\omega_{\pm})&=2R(\omega_{\pm})^{3}-\frac{6(p_{\pm}+1)}{\omega_{\pm}}R(\omega_{\pm})^{2}+\left(\frac{4(p_{\pm}+1)(p_{\pm}+2)}{\omega_{\pm}^{2}}-2\right)R(\omega_{\pm})+\frac{2(p_{\pm}+2)}{\omega_{\pm}},\\\
W_{4}(\omega_{\pm})&=-6R(\omega_{\pm})^{4}+\frac{24(p_{\pm}+1)}{\omega_{\pm}}+\left(8-\frac{4(p_{\pm}+1)(7p_{\pm}+11)}{\omega_{\pm}^{2}}\right)R(\omega_{\pm})^{2}\\\
&+\left(\frac{8(p_{\pm}+1)(p_{\pm}+2)(p_{\pm}+3)}{\omega_{\pm}^{3}}-\frac{4(4p_{\pm}+5)}{\omega_{\pm}}\right)R(\omega_{\pm})+\frac{4(p_{\pm}+1)(p_{\pm}+2)}{\omega_{\pm}^{2}}-2.\end{aligned}\right.$

In (17),
$R(\omega_{\pm})=K_{p_{\pm}+1}(\omega_{\pm})/K_{p_{\pm}}(\omega_{\pm})$, where
the function $K_{\nu}$ stands for the modified Bessel function of second kind
(MacDonald function) [1].

######  Proposition 2.2.

Let $X\sim\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$, we have:

(18) $\forall
n\in\mathbb{N},\quad\mathbb{E}\left[X^{n}\right]=\frac{1}{K_{p_{+}}\left(\sqrt{a_{+}b_{+}}\right)K_{p_{-}}\left(\sqrt{a_{-}b_{-}}\right)}\\\
\times\sum_{k=0}^{n}\binom{n}{k}(-1)^{n-k}K_{p_{+}+k}\left(\sqrt{a_{+}b_{+}}\right)K_{p_{-}+n-k}\left(\sqrt{a_{-}b_{-}}\right).$

###### Proof.

Using definition 11 and the expression of the moment for GIG distributions
(see [30]), the proof is straightforward. ∎

####  2.3.4. Mellin transform

In this subsection, we provide a representation of the BGIG density in the
Mellin space, that can be used as an alternative to the representation in the
Fourier space (12). We start by a lemma on GIG densities:

######  Lemma 2.2.

For every $x>0$, the GIG density admits the representation:

(19)
$f_{\mathrm{GIG}_{+}}(x)=\int\limits_{c-\mathrm{i}\infty}^{c+\mathrm{i}\infty}\frac{K_{s+p-1}(\sqrt{ab})}{K_{p}(\sqrt{ab})}\left(\frac{a}{b}\right)^{\frac{1-s}{2}}x^{-s}\frac{\mathrm{d}s}{2\mathrm{i}\pi},\quad
c\in\mathbb{R}.$

###### Proof.

By definition of the Mellin transform [23, Eq. 1], the Mellin transform of the
GIG density is:

(20) $\displaystyle f_{\mathrm{GIG}_{+}}^{*}(s)$
$\displaystyle=\frac{(a/b)^{p/2}}{2K_{p}\left(\sqrt{ab}\right)}\int_{0}^{\infty}x^{s+p-2}e^{-\left(\frac{ax+b/x}{2}\right)}\mathrm{d}x$
(21)
$\displaystyle=\left(\frac{a}{b}\right)^{\frac{1-s}{2}}\frac{K_{s+p-1}(\sqrt{ab})}{K_{p}(\sqrt{ab})}$

where the integral converges for all $s\in{\mathbb{C}}$. Applying the Mellin
inversion theorem [23, Eq. 7] yields (19). ∎

######  Proposition 2.3.

For every $x>0$, the BGIG density admits the representation:

(22)
$f(x)=\sqrt{\frac{a_{+}a_{-}}{b_{+}b_{-}}}\int\limits_{c+\mathrm{i}\mathbb{R}^{2}}\frac{\Gamma(s_{1}+s_{2}-1)\Gamma(1-s_{2})}{\Gamma(s_{1})}\times\frac{K_{s_{1}+p_{+}-1}(\sqrt{a_{+}b_{+}})K_{s_{2}+p_{-}-1}(\sqrt{a_{-}b_{-}})}{K_{p_{+}}(\sqrt{a_{+}b_{+}})K_{p_{-}}(\sqrt{a_{-}b_{-}})}\\\
\times\left(\frac{a_{+}}{b_{+}}\right)^{-\frac{s_{1}}{2}}\left(\frac{a_{-}}{b_{-}}\right)^{-\frac{s_{2}}{2}}x^{1-s_{1}-s_{2}}\frac{\mathrm{d}s_{1}\mathrm{d}s_{2}}{(2\mathrm{i}\pi)^{2}},$

where $c$ belongs to the polyhedron
$P=\\{(s_{1},s_{2})\in{\mathbb{C}}^{2},\mathrm{Re}(s_{1})>0,\mathrm{Re}(s_{2})<1,\mathrm{Re}(s_{1}+s_{2}-1)>0\\}$.

###### Proof.

From definition (11), the BGIG density can be written as the convolution of
two GIG densities as follows:

(23)
$f(x)=\int_{-\infty}^{\infty}f_{\mathrm{GIG}_{+}}(y,a_{+},b_{+},p_{+})f_{\mathrm{GIG}_{-}}(x-y,a_{-},b_{-},p_{-})\mathrm{d}y,$

where we have used a slight abuse of notations to indicate the parameters
involved in each density function. Using definitions (2) and (3) we are left
with, for $x>0$:

(24)
$f(x)=\int\limits_{0}^{\infty}f_{\mathrm{GIG}_{+}}(x+y,a_{+},b_{+},p_{+})f_{\mathrm{GIG}_{+}}(y,a_{-},b_{-},p_{-})\mathrm{d}y$

From lemma 2.2, we can therefore write:

(25)
$f(x)=\int\limits_{c_{1}-\mathrm{i}\infty}^{c_{1}+\mathrm{i}\infty}\int\limits_{c_{2}-\mathrm{i}\infty}^{c_{2}+\mathrm{i}\infty}\frac{K_{s_{1}+p_{+}-1}(\sqrt{a_{+}b_{+}})K_{s_{2}+p_{-}-1}(\sqrt{a_{-}b_{-}})}{K_{p_{+}}(\sqrt{a_{+}b_{+}})K_{p_{-}}(\sqrt{a_{-}b_{-}})}\left(\frac{a_{+}}{b_{+}}\right)^{\frac{1-s_{1}}{2}}\left(\frac{a_{-}}{b_{-}}\right)^{\frac{1-s_{2}}{2}}\\\
\times\int_{0}^{\infty}(x+y)^{-s_{1}}y^{-s_{2}}\mathrm{d}y\frac{\mathrm{d}s_{1}\mathrm{d}s_{2}}{(2\mathrm{i}\pi)^{2}}$

where the $y$ integral is a particular case of a Beta integral:

(26)
$\int_{0}^{\infty}(x+y)^{-s_{1}}y^{-s_{2}}\mathrm{d}y=\frac{\Gamma(s_{1}+s_{2}-1)\Gamma(1-s_{2})}{\Gamma(s_{1})}x^{1-s_{1}-s_{2}},$

which converges if $\mathrm{Re}(s_{1})>0$, $\mathrm{Re}(1-s_{2})>0$ and
$\mathrm{Re}(-1+s_{1}+s_{2})>0$. Applying the Mellin inversion theorem yields
(22). ∎

######  Remark 2.

It is straightforward to extend proposition 2.3 to the case $x<0$, as it
follows from the characteristic function (12) that under the transformation
$x\rightarrow-x$ the BGIG parameters are transformed as:

(27)
$(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})\rightarrow(a_{-},b_{-},p_{-},a_{+},b_{+},p_{+}).$

See also proposition 4.5 in section 4.

##  3\. Particular cases

In this section, we study some particular cases of BGIG distributions, as well
as some of their limiting cases.

###  3.1. Decomposition of BGIG distributions

Let $D$ be the distribution defined by its characteristic function:

(28) $\forall
u\in{\mathbb{R}},\quad\Phi_{\mathrm{D}}(u)=\frac{K_{p_{+}}(\sqrt{b_{+}(a_{+}-\mathrm{i}2u)})K_{p_{-}}(\sqrt{b_{-}(a_{-}+\mathrm{i}2u)})}{K_{p_{+}}(\sqrt{a_{+}b_{+}})K_{p_{-}}(\sqrt{a_{-}b_{-}})}.$

From equation 13, it is immediate to see that any BGIG distribution with
parameters $p_{\pm}>0$ can be written as

(29)
$\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})=\mathrm{D}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})\star\mathrm{BG}(p_{+}/2,a_{-}/2,p_{-}/2,a_{-}/2),$

where $\mathrm{BG}(\alpha_{+},\lambda_{+},\alpha_{-},\lambda_{-})$ denotes the
bilateral Gamma (BG) distribution whose parameters are strictly positive and
whose characteristic function is (see details in [36, 37, 3])

(30)
$\Phi_{\mathrm{BG}}(u)=\left(\frac{\lambda_{+}}{\lambda_{+}-\mathrm{i}u}\right)^{\alpha_{+}}\left(\frac{\lambda_{-}}{\lambda_{-}+\mathrm{i}u}\right)^{\alpha_{-}}$

for all $u\in\mathbb{R}$.

######  Remark 3.

Let us notice that, even if they can be represented by a convolution with a BG
density (in the case $p_{\pm}>0$), BGIG distributions are not a subclass of
tempered stable (TS) distributions [35]; indeed from lemma 2.1 and definition
13, it is immediate to see that for $p_{\pm}\leq 0$ the BGIG Lévy measure has
the following behaviours:

(31) $\left\\{\begin{aligned}
\pi(\mathrm{d}x)&=\mathcal{O}\left(\frac{e^{-a_{\pm}|x|/2}}{|x|^{1+|p_{\pm}|}}\right)\mathrm{d}x\quad\mathrm{as}\quad
x\to\pm\infty,\\\
\pi(\mathrm{d}x)&=\mathcal{O}\left(x^{-3/2}\right)\mathrm{d}x\quad\mathrm{as}\quad
x\to 0.\end{aligned}\right.$

Such difference between asymptotics in $\pm\infty$ and 0 is impossible to
obtain with TS distributions whose Lévy measure is given by

(32) $\forall
x\in{\mathbb{R}},\quad\pi_{\mathrm{TS}}(\mathrm{d}x)=\left(\frac{\alpha_{+}}{x^{1+\beta_{+}}}e^{-\lambda_{+}x}\mathbbm{1}_{(0,+\infty)}(x)+\frac{\alpha_{-}}{|x|^{1+\beta_{-}}}e^{-\lambda_{-}|x|}\mathbbm{1}_{(-\infty,0)}(x)\right)\mathrm{d}x$

where $\alpha_{\pm},\lambda_{\pm}>0$ (like in the BG case) and
$\beta_{\pm}\in(0,1)$. The only exception occurs for $p_{\pm}=-1/2$: in this
case both behaviours in $\eqref{eq:pi_behav}$ coincide and the BGIG
distribution degenerates into a bilateral inverse Gaussian distribution (see
section 3.4), which is a TS distribution with $\beta_{\pm}=1/2$. For
$p_{\pm}>0$, the closure of BGIG distributions recovers other TS
distributions, as we will see in the next sections.

###  3.2. Bilateral Gamma and Gamma distributions

For $p_{\pm}>0$, we have the following convergence in law:

(33)
$\mathrm{BGIG}\left(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-}\right)\overset{\mathcal{L}}{\underset{b_{+},b_{-}\to
0}{\longrightarrow}}\mathrm{BG}(p_{+},a_{+}/2,p_{-},a_{-}/2)$

This is an immediate consequence of the weak convergence of each of the GIG
subordinators to Gamma subordinators (see [36, Eq. 3.3]). Letting $p_{-}\to 0$
recovers the usual Gamma distribution, which is supported by the positive
semi-axis only.

######  Remark 4.

BG distributions belong then to the closure of BGIG class of distributions,
however, as proved in [37], these distributions are not smooth whereas BGIG
distributions are (see proposition 4.1).

###  3.3. Variance Gamma and symmetric Variance Gamma distributions

If $p_{\pm}=p>0$, the following convergence holds:

(34)
$\mathrm{BGIG}\left(a_{+},b_{+},p,a,b_{-},p\right)\overset{\mathcal{L}}{\underset{b_{+},b_{-}\to
0}{\longrightarrow}}\mathrm{VG}(\sigma^{2},\theta,\nu),$

where $\mathrm{VG}$ denotes the Variance Gamma distribution (see [41, 4]); it
is usually parametrized via a scale (volatility) parameter $\sigma$, an
asymmetry (or skewness) parameter $\theta$ and a kurtosis parameter $\nu$, and
has the following characteristic function:

(35) $\Phi_{\mathrm{VG}}(u)=\left(1-i\theta\nu
u+\frac{\sigma^{2}\nu}{2}u^{2}\right)^{-\frac{1}{\nu}}$

Comparing (35) and (12) in the limits $b\pm\to 0$ and $p_{\pm}=p$ allows to
express the VG parameters in terms of the BGIG parameters as (see also [3]):

(36) $\left\\{\begin{aligned} &\sigma^{2}=8p(a_{+}a_{-})^{-1},\\\
&\theta=2p(a_{+}^{-1}-a_{-}^{-1}),\\\ &\nu=p^{-1}.\end{aligned}\right.$

The VG distribution is itself a particular case of the CGMY distribution [15]
with $Y=1$ (which are themselves particular cases of the more general KoBoL
distributions [14]). If furthermore $a_{\pm}=a$, then the BGIG distribution
becomes the symmetric Variance Gamma distribution (see [42, 2]) for which the
asymmetry parameter $\theta$ is equal to zero.

###  3.4. Half-integer values of $p$

In this subsection, we focus on negative parameters $p_{\pm}$ only since, for
$p_{\pm}$ positive, the distribution is simply the convolution with a
bilateral Gamma distribution (see section 3.1). For simplicity we restrict
ourselves to the case $p_{+}=p_{-}=p$.

####  3.4.1. Case $p=-1/2$: the “bilateral inverse Gaussian” distribution

For $p=-1/2$, the Jaeger integral, which is symmetric in $p$, is available in
closed form (see [5]):

(37) $\mathscr{I}(x,1/2)=\frac{1}{2\sqrt{\pi x}}$

and the Lévy measure simplifies to:

(38) $\forall
x\in{\mathbb{R}},\quad\pi(\mathrm{d}x)=\left(\sqrt{\frac{b_{+}}{2\pi
x^{3}}}e^{-a_{+}x/2}\mathbbm{1}_{(0,\infty)}(x)+\sqrt{\frac{b_{-}}{2\pi|x|^{3}}}e^{-a_{-}|x|/2}\mathbbm{1}_{(-\infty,0)}(x)\right)\mathrm{d}x.$

We recognize the distribution of the difference of two inverse Gaussian
distributions and we can write:

(39)
$\mathrm{BGIG}(-1/2,a_{+},b_{+},-1/2,a_{-},b_{-})=\mathrm{IG}_{+}(a_{+},b_{+})*\mathrm{IG}_{-}(a_{-},b_{-}).$

This distribution is a special case of Tempered Stable (TS) distribution with,
under the parametrization of (32), TS parameters equal to
$\lambda_{\pm}=a_{\pm}/2$, $\beta_{\pm}=1/2$ and
$\alpha_{\pm}=\sqrt{b_{\pm}/2\pi}$; it corresponds to the distribution of the
difference of two independent inverse Gaussian variables.

####  3.4.2. Case $p=-3/2$.

In general, BGIG distributions are not a particular case of TS distributions
(the case $p=-1/2$ being an exception). Indeed, let us assume $p=-3/2$, in
that case we have (see [5]):

(40) $\mathscr{I}(x,3/2)=\frac{1}{2}\left(\frac{1}{\sqrt{\pi
x}}-e^{x}\mathrm{erfc}(\sqrt{x})\right).$

The corresponding Lévy measure is then:

(41) $\forall
x\in{\mathbb{R}},\quad\pi(\mathrm{d}x)=\left(\sqrt{\frac{b_{+}}{2\pi
x^{3}}}e^{-a_{+}x/2}\mathbbm{1}_{(0,\infty)}(x)+\sqrt{\frac{b_{-}}{2\pi|x|^{3}}}e^{-a_{-}|x|/2}\mathbbm{1}_{(-\infty,0)}(x)\right.\\\
\left.+x^{-1}\mathrm{erfc}(\sqrt{x/2b_{+}})e^{x/2b_{+}-a_{+}x/2}\mathbbm{1}_{(0,\infty)}(x)+|x|^{-1}\mathrm{erfc}(\sqrt{|x|/2b_{-}})e^{|x|/2b_{-}-a_{-}|x|/2}\mathbbm{1}_{(-\infty,0)}(x)\right)\mathrm{d}x.\\\
$

where erfc stands for the complementary error-function. The first line of
equation 41 is the Lévy measure of a TS distribution with $\beta_{\pm}=1/2$,
but the second line is a more complex component showing that BGIG
distributions are not a subclass of Tempered-Stable distributions.

####  3.4.3. Case $p=-(n+1/2)$.

More generally, it is shown in [5] that for $p=n+1/2$ with
$n\in\mathbb{N}_{>1}$, we have:

(42)
$\mathscr{I}(x,n+1/2)=\frac{1}{\pi}\int_{0}^{\infty}\frac{u^{2n}e^{-xu^{2}}}{\prod_{i=1}^{n}(u^{2}+\alpha_{j})}\mathrm{d}u$

and therefore:

(43)
$\mathscr{I}(x,n+1/2)=\frac{B_{0}}{\sqrt{x}}+\sum_{j=1}^{n}B_{j}e^{\alpha_{j}^{2}x}\left(\frac{\alpha_{j}}{|\alpha_{j}|}-1+\mathrm{erfc}(\alpha_{j}\sqrt{x})\right)$

where the $B_{j}$ are constants and the $\alpha_{j}$ are the zeros of the
modified Bessel function of second kind. The particular values of the Jaeger
integrals generate a component involving the erfc function in the Lévy
measure, demonstrating yet again that these distributions do not belong to the
TS class.

##  4\. Shape and asymptotics of BGIG distributions

Let us now study the smoothness and analyticity of BGIG distributions, as well
as some asymptotic behaviours for large arguments of BGIG random variables as
well as their maximum and minimum.

###  4.1. Smoothness

In this section, we prove that BGIG distributions are strictly unimodal and
smooth, which can be seen in the example of a BGIG distribution $f$ displayed
in Figure 1. First, remark that $x\mapsto\mathscr{I}(x,p)$ is decreasing on
$(0,+\infty)$ (which is immediate by differentiating with respect to $x$ and
noticing that the integrand of (5) is negative). As a consequence, if we re-
write the BGIG Lévy measure (13) as:

(44) $\forall
x\in{\mathbb{R}},\quad\pi(\mathrm{d}x)=\frac{k(x)}{x}\mathrm{d}x$

where $k$ is defined as:

(45) $\forall x\in{\mathbb{R}},\quad
k(x):=e^{-a_{+}x/2}\left(\mathscr{I}(x/2b_{+},p_{+})+\max(p_{+},0)\right)\mathbbm{1}_{(0,\infty)}(x)\\\
-e^{-a_{-}|x|/2}\left(\mathscr{I}(|x|/2b_{-},p_{-})+\max(p_{-},0)\right)\mathbbm{1}_{(-\infty,0)}(x),$

then it follows from lemma 2.1, that $x\rightarrow k(x)$ is decreasing on both
$(-\infty,0)$ and $(0+\infty)$. From [44, Corollary 15.11], this implies that
BGIG distributions are self-decomposable and hence of class $L$ in the
terminology of [45].

Figure 1. $\mathrm{BGIG}(1,2,1,3,4,5)$ distribution and its two underlying GIG
distributions.

######  Proposition 4.1.

Densities of BGIG distributions are smooth, i.e., they belong to
$\mathscr{C}^{\infty}({\mathbb{R}},{\mathbb{R}})$.

###### Proof.

From lemma 2.1, we have that $k(0^{+})=|k(0^{-})|=+\infty$ and thus, neither
the density nor its reflection is of type $I_{4}$ in the sense of [45]. The
use of [45, Theorem 1.9] concludes the proof. ∎

###  4.2. Analyticity

Before discussing the analyticity of the density, let us remark that, for
$a,b>0$, we have:

(46) $\left\\{\begin{aligned} &\arg\sqrt{b(a\pm
2\mathrm{i}u)}\underset{u\to\infty}{\longrightarrow}\pm\frac{\pi}{4},\\\
&\left|\sqrt{b(a\pm
2\mathrm{i}u)}\right|=\sqrt{b}(a^{2}+4u^{2})^{1/4}=\sqrt{2bu}(1+\mathcal{O}(u^{-2})).\end{aligned}\right.$

The following Lemma will also be useful.

######  Lemma 4.1.

If $X\sim\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$, then its
characteristic function $\Phi$ satisfies:

(47)
$\mathrm{Re}\left[\Phi(u)\right]\underset{u\to\infty}{\sim}Ku^{-(1+\bar{p})/2}e^{-\sqrt{u}\bar{b}}\cos\left(-\underaccent{\bar}{b}\sqrt{u}+\frac{\pi}{4}\underaccent{\bar}{p}\right),$

where:

(48) $\left\\{\begin{aligned} K&:=\frac{\pi
a_{+}^{p_{+}/2}a_{-}^{p_{-}/2}}{2^{(\bar{p}+3)/2}(b_{+}b_{-})^{1/4}K_{p_{+}}(\sqrt{a_{+}b_{+}})K_{p_{-}}(\sqrt{a_{-}b_{-}})},\\\
\bar{p}&:=p_{+}+p_{-},\\\ \underaccent{\bar}{p}&:=p_{+}-p_{-},\\\
\bar{b}&:=\sqrt{b_{+}}+\sqrt{b_{-}},\\\
\underaccent{\bar}{b}&:=\sqrt{b_{-}}-\sqrt{b_{+}}.\end{aligned}\right.$

###### Proof.

The BGIG characteristic function 12 can be written as:

(49)
$\Phi(u)=M\frac{K_{p_{+}}\left(\sqrt{b_{+}(a_{+}-\mathrm{i}2u)}\right)K_{p_{-}}\left(\sqrt{b_{-}(a_{-}+\mathrm{i}2u)}\right)}{(a_{+}-\mathrm{i}2u)^{p_{+}/2}(a_{-}+\mathrm{i}2u)^{p_{-}/2}}.$

where
$M=a_{+}^{p_{+}/2}a_{-}^{p_{-}/2}/(K_{p_{+}}(\sqrt{a_{+}b_{+}})K_{p_{-}}(\sqrt{a_{-}b_{-}}))$.
We first note that, with notations (48), equation LABEL:eq:arg-module gives:

(50)
$(a_{+}-\mathrm{i}2u)^{-p_{+}/2}(a_{-}+\mathrm{i}2u)^{-p_{-}/2}=(2u)^{-\bar{p}/2}e^{\mathrm{i}\frac{\pi}{4}\underaccent{\bar}{p}}(1+\mathcal{O}(u^{-2})).$

Since $\lim_{u\to+\infty}|\arg{\sqrt{b(a\pm 2\mathrm{i}u)}}|<3\pi/2$, we have
by [18, Eq. 10.40.2]:

(51)
$K_{p_{\pm}}\left(\sqrt{b_{\pm}(a_{\pm}\mp\mathrm{i}2u)}\right)\underset{u\to\infty}{\sim}\sqrt{\frac{\pi}{2\sqrt{b_{\pm}}(a_{\pm}^{2}+4u^{2})^{1/4}}e^{\pm\mathrm{i}\pi/4}}e^{-\sqrt{b_{\pm}}(a_{\pm}^{2}+4u^{2})^{1/4}e^{\mp\mathrm{i}\pi/4}}.$

Inserting this asymptotic in equation LABEL:eq:arg-module, and re-arranging
the terms, we get:

(52)
$K_{p_{+}}\left(\sqrt{b_{+}(a_{+}-\mathrm{i}2u)}\right)K_{p_{-}}\left(\sqrt{b_{-}(a_{-}+\mathrm{i}2u)}\right)\underset{u\to\infty}{\sim}\frac{\pi}{2\sqrt{2}(b_{+}b_{-})^{1/4}\sqrt{u}}e^{-\sqrt{u}(\bar{b}+\mathrm{i}\underaccent{\bar}{b})}.$

Combining 52 and 50 yields:

(53)
$\Phi(u)\underset{u\to\infty}{\sim}Ku^{-(1+\bar{p})/2}e^{-\sqrt{u}(\bar{b}+\mathrm{i}\underaccent{\bar}{b})+\mathrm{i}\frac{\pi}{4}\underaccent{\bar}{p}}.$

Taking the real part concludes the proof. ∎

######  Proposition 4.2.

Let $f$ be the density function of the
$\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$-distribution. The
following holds:

  1. (1)

$f$ is not real-analytic in 0,

  2. (2)

$f$ is real-analytic on ${\mathbb{R}}\backslash\\{0\\}$.

###### Proof.

To prove these results, we use the Fourier–Bros–Iagolnitzer (FBI) transform
defined as [28, 33]:

(54) $T_{t}[\varphi](t,\xi):=\int_{\mathbb{R}}\varphi(s)e^{-\mathrm{i}2\pi
s\xi}e^{-\pi t(s-x)^{2}}\mathrm{d}s$

where $(t,\xi,x)\in(0,\infty)\times{\mathbb{R}}^{2}$ and $\varphi\in
L^{1}({\mathbb{R}})$. In the case of the BGIG density (i.e. taking $\varphi=f$
in (54)), we therefore have:

(55) $T_{t}[f](t,t\xi)=e^{-\pi x^{2}t}\left(\Phi\star
g_{t}\right)(t\xi+\mathrm{i}tx)$

where $g_{t}(u):=(\pi t)^{-1/2}\exp(-u^{2}/\pi t)$. This equality is obtained
by expanding the square in the exponential term, using the convolution
property of the Fourier transform and noticing that $\mathcal{F}[f]=\Phi$.

  * $\bullet$

When $x\neq 0$, it is easy to see that $t\mapsto\Phi(t\xi+\mathrm{i}tx)$ is
bounded on $(0,\infty)$, it is then immediate that $\Phi\star g_{t}$ is also
bounded. As a consequence, there exists $K>0$ such that:

(56) $|T_{t}[f](t,t\xi)|\leq Ke^{-\pi x^{2}t}$

which, in virtue of [33, Thm. 5.3.2], ensures analyticity of $f$ for $x\neq
0$.

  * $\bullet$

For $x=0$, we have:

(57) $T_{t}[f](t,t\xi)=\left(\Phi\star g_{t}\right)(t\xi).$

The convolution can be approximated by the introducing a truncation centered
around $t\xi$ and of size $2C\sqrt{t}$ for $C>0$:

(58) $\left(\Phi\star
g_{t}\right)(t\xi)=\int_{\mathbb{R}}g_{t}(t\xi-u)\Phi(u)\mathrm{d}u=\int_{t\xi-C\sqrt{t}}^{t\xi+C\sqrt{t}}g_{t}(t\xi-u)\Phi(u)\mathrm{d}u+\mathcal{O}\left(C^{-1}e^{-C^{2}/\pi}\right)$

where the $\mathcal{O}(C^{-1}e^{-C^{2}/\pi})$ terms comes from the asymptotic
behavior of the complementary error function $\mathrm{erfc}$ (see [18, Eq.
7.12.1]). Applying the Mean Value theorem to the real part of the truncated
integral in 58, there exists $\eta_{t}\in[t\xi-C\sqrt{t},t\xi+C\sqrt{t}]$ such
that:

(59)
$\mathrm{Re}\left[\int_{t\xi-C\sqrt{t}}^{t\xi+C\sqrt{t}}g_{t}(t\xi-u)\Phi(u)\mathrm{d}u\right]=2C\sqrt{t}g_{t}(t\xi-\eta_{t})\mathrm{Re}\left[\Phi(\eta_{t})\right].$

Since $\eta_{t}\in[t\xi-C\sqrt{t},t\xi+C\sqrt{t}]$, we have
$g_{t}(t\xi-\eta_{t})=\mathcal{O}(t^{-1/2}C^{-1}e^{-C^{2}/\pi t}))$. If we
assume $\xi\neq 0$, then $\eta_{t}\to\infty$ when $t\to\infty$ and thus, from
the asymptotic of Lemma 4.1, we have:

(60)
$\mathrm{Re}\left[\int_{t\xi-C\sqrt{t}}^{t\xi+C\sqrt{t}}g_{t}(t\xi-u)\Phi(u)\mathrm{d}u\right]=\\\
\mathcal{O}\left(e^{-C^{2}/(\pi
t)}(t\xi)^{-(1+\bar{p})/2}e^{-\sqrt{t\xi}\bar{b}}\cos\left(-\underaccent{\bar}{b}\sqrt{t\xi}+\frac{\pi}{4}\underaccent{\bar}{p}\right)\right).$

Noting that:

(61) $\left|\left(\Phi\star
g_{t}\right)(t\xi)\right|\geq\left|\mathrm{Re}\left[\left(\Phi\star
g_{t}\right)(t\xi)\right]\right|$

shows that there does not exist $M>0$ such that for $|\xi|>M$,
$\left(\Phi\star g_{t}\right)(t\xi)$ decays as $e^{-at}$ for some positive $a$
(here (60) shows that the decay is at best in $e^{-a\sqrt{t}}$). Thus, $f$
does not satisfy the RA(0) condition of [33, Thm. 5.3.2] and therefore is not
analytic in $x=0$.

∎

###  4.3. Unimodality

######  Proposition 4.3.

BGIG distributions are strictly unimodal, i.e., there exists a unique
$x_{0}\in{\mathbb{R}}$ such that their density function is strictly increasing
on $(-\infty,x_{0})$ and strictly decreasing on $(x_{0},+\infty)$.

###### Proof.

The same argument exposed in the proof of proposition 4.1 and the use of [45,
Theorem 1.2] instead of [45, Theorem 1.9] concludes the proof. ∎

######  Proposition 4.4.

The mode $x_{0}$ of a BIG distribution is located on $(0,\infty)$ if and only
if:

(62)
$(a_{-}b_{+}-a_{+}b_{-})K_{1-p_{+}-p_{-}}\left(\sqrt{\Lambda}\right)>\frac{2(\Lambda-(a_{+}+a_{-})(b_{+}p_{-}+b_{-}p_{+}))K_{2-p_{-}-p_{+}}\left(\sqrt{\Lambda}\right)}{\sqrt{b_{-}+b_{+}}}$

where $\Lambda:=(a_{+}+a_{-})(b_{+}+b_{-})$. The mode is located $(-\infty,0)$
if and only if the inequality is strictly reversed. Finally, $x_{0}=0$ if and
only if the inequality becomes an equality.

###### Proof.

Using proposition 4.1 and simplifying positive multiplicative terms, we get:

(63)
$\mathrm{sgn}[f_{\mathrm{BGIG}}^{\prime}(0)]=\mathrm{sgn}\bigg{[}[(a_{-}b_{+}-a_{+}b_{-})K_{1-p_{+}-p_{-}}\left(\sqrt{(a_{+}+a_{-})(b_{+}+b_{-})}\right)\\\
+\frac{2((a_{+}+a_{-})(b_{+}(1-p_{-})+b_{-}(1-p_{+})))K_{2-p_{-}-p_{+}}\left(\sqrt{(a_{+}+a_{-})(b_{+}+b_{-})}\right)}{\sqrt{b_{-}+b_{+}}}\bigg{]}$

and the result follows immediately from proposition 4.3 and proposition 4.1
that ensures continuity of the first derivative. ∎

###  4.4. Asymptotics

######  Proposition 4.5.

Let $a_{\pm},b_{\pm}>0$, $p_{\pm}\in{\mathbb{R}}$, and denote by $f$ the
density function of the distribution
$\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$ and by $\widetilde{f}$
the density function of the distribution $\mathrm{BGIG}(a_{-},b_{-},$
$p_{-},a_{+},b_{+},p_{+})$. We have the symmetry relation:

(64) $\forall x\in{\mathbb{R}},\quad f(x)=\widetilde{f}(-x).$

###### Proof.

This is an immediate consequence of the subsitution $u\rightarrow-u$ in the
BGIG characteristic function 12. ∎

######  Proposition 4.6.

Let $f$ be the density function of a
$\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$ distribution. We have the
asymptotic behaviour:

(65) $\displaystyle
f(x)\underset{x\to\pm\infty}{\sim}Z_{\pm}|x|^{p_{\pm}-1}\exp\left(-\frac{a_{\pm}|x|}{2}\right)$

where:

(66)
$Z_{\pm}:=\left(\frac{a_{\pm}}{b_{\pm}}\right)^{p_{\pm}/2}\left(\frac{a_{\pm}}{a_{\pm}+a_{\mp}}\right)^{p_{\mp}/2}\frac{K_{p_{\mp}}\left(\sqrt{(a_{\mp}+a_{\pm})b_{\mp}}\right)}{2K_{p_{\mp}}\left(\sqrt{a_{\mp}b_{\mp}}\right)K_{p_{\pm}}(\sqrt{a_{\pm}b_{\pm}})}.$

###### Proof.

By definition of the BGIG distribution as a convolution of GIG distributions
(11), we have

(67) $\forall x\geq 0,\quad
f(x)=\int_{0}^{\infty}f_{\mathrm{GIG}_{+}}(y+x)f_{\mathrm{GIG}_{-}}(-y)\mathrm{d}y.$

It is immediate to see that:

(68) $\quad
f_{\mathrm{GIG}_{+}}(y)\underset{y\to+\infty}{\sim}\alpha_{+}y^{p_{+}-1}e^{-\frac{a_{+}y}{2}}$

where
$\alpha_{+}:=\left(a_{+}/b_{+}\right)^{p_{+}/2}(2K_{p_{+}}(\sqrt{a_{+}b_{+}}))^{-1}$.
Inserting this asymptotic in (67) gives:

(69)
$f(x)\underset{x\to+\infty}{\sim}\alpha_{+}x^{p_{+}-1}e^{-\frac{a_{+}x}{2}}\int_{0}^{\infty}(1+y/x)^{p_{+}-1}e^{-\frac{a_{+}y}{2}}f_{\mathrm{GIG}_{-}}(-y)\mathrm{d}y.$

By dominated convergence theorem and the definition of the
$\mathrm{GIG}_{-}(a,b,p)$ distribution, we can write:

(70)
$\int_{0}^{\infty}(1+y/x)^{p_{+}-1}e^{-\frac{a_{+}y}{2}}f_{\mathrm{GIG}_{-}}(-y)\mathrm{d}y\underset{x\to+\infty}{\longrightarrow}\left(\frac{a_{+}}{a_{+}+a_{-}}\right)^{p_{-}/2}\frac{K_{p_{-}}\left(\sqrt{(a_{-}+a_{+})b_{-}}\right)}{K_{p_{-}}\left(\sqrt{a_{-}b_{-}}\right)}.$

Using the symmetry relation in proposition 4.5 leads to the analog result for
$x\to-\infty$. ∎

The next proposition is an application of proposition 4.6 to the maximum and
minimum of BGIG processes; it will prove to be particulary useful in section 7
for the calibration of BGIG parameters from empirical data.

######  Proposition 4.7.

For $n\in\mathbb{N}$, let $(X_{i})_{1\leq i\leq n}$ be a sequence of i.i.d.
BGIG variables and denote $\displaystyle M_{n}:=\max_{1\leq i\leq n}X_{i}$ and
$\displaystyle m_{n}:=\min_{1\leq i\leq n}X_{i}$. We have the following
$\mathbb{P}$-almost sure asymptotic behaviours:

(71) $\begin{cases}M_{n}\underset{n\to+\infty}{\sim}\frac{2}{a_{+}}\ln n,\\\
m_{n}\underset{n\to+\infty}{\sim}-\frac{2}{a_{-}}\ln n.\end{cases}$

###### Proof.

From proposition 4.6 and the asymptotic behavior of the upper incomplete gamma
function (see [18, Eq. 8.11.2]), we have

(72)
$\mathbb{P}\left(X_{i}>x\right)\underset{x\to+\infty}{\sim}\Gamma(p_{+},x)\underset{x\to+\infty}{\sim}Z_{+}x^{p_{+}-1}\exp\left(-\frac{a_{+}x}{2}\right).$

Let $\widetilde{M}_{n}:=\frac{2}{a_{+}}\ln n$. For $\alpha>0$, we have:

(73)
$\mathbb{P}\left(M_{n}>\alpha\widetilde{M}_{n}\right)=1-\mathbb{P}\left(X_{i}<\alpha\widetilde{M}_{n}\quad\forall
i=1\dots n\right)=1-(1-\mathbb{P}(X_{i}>\alpha M_{n}))^{n}.$

By equation 72, we have (using $\ln(1+x)=x+\mathcal{O}(x^{2})$ for small $x$):

(74) $\displaystyle\mathbb{P}\left(M_{n}>\alpha\widetilde{M}_{n}\right)$
$\displaystyle=1-\exp\left(-n^{1-\alpha}\ln(n)^{p_{+}-1}\frac{(Z_{+}+o(1))a_{+}(\alpha\widetilde{M}_{n})^{p_{+}-1}}{2}\right)$
$\displaystyle\underset{n\to+\infty}{\longrightarrow}\begin{cases}1\quad\mathrm{if}\quad\alpha<1,\\\
0\quad\mathrm{if}\quad\alpha>1.\\\ \end{cases}$

The symmetry relation in proposition 4.5 leads to the analog result for the
minimum. ∎

##  5\. Bilateral generalized inverse Gaussian processes

In this section, we define the BGIG process and provide some of its
properties. We also study its Blumenthal-Getoor index and its stability under
a specific change of measure, as well as the jump structure of its sample
paths.

###  5.1. Definition and first properties

As the BGIG distribution is infinitely divisible (proposition 2.1), there
exists an associated Lévy process, said to be generated by the distribution
$\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$. Denoting such process by
$X=(X_{t})_{t\geq 0}$, we will say, with a slight abuse of language, that $X$
is a $\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$ process, or simply a
BGIG process; as before $(a_{+},b_{+},a_{-},b_{-})\in(0,+\infty)^{4}$ and
$(p_{+},p_{-})\in{\mathbb{R}}^{2}$.

######  Remark 5.

From remark 1, we see that the difference of two GIG processes is a BGIG
process (where the GIG process is the Lévy process generated by the GIG
distribution $\mathrm{GIG}_{+}(a,b,p)$).

As for all Lévy processes, the characteristic function of the BGIG process be
simply expressed as:

(75) $\forall(u,t)\in{\mathbb{R}}\times(0,+\infty),\quad\Phi(u,t)=\Phi(u)^{t}$

where $\Phi$ is the characteristic function of the generating distribution 12.
Of course, it follows that
$X_{1}\sim\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$. The Lévy-
Khintchine representation for the characteristic function (75) follows
immediately from (14) and reads

(76) $\forall u\in{\mathbb{R}},\forall
t\in(0,+\infty),\quad\Phi(u,t)=\exp\left(t\int_{\mathbb{R}}\left(e^{\mathrm{i}ux}-1\right)\pi(\mathrm{d}x)\right).$

where $\pi$ is the Lévy measure defined in (13). This shows that the BGIG
process is a pure jump Lévy process, whose characteristic triplet is
$(0,0,\pi)$ . We display an example of its transition densities in Figure 2
for different values of $t$, obtained by inverting the characteristic function
in eq. 75.

Figure 2. Transition densities for a $\mathrm{BGIG}(1,2,1,3,4,5)$-process and
different values of $t$.

######  Proposition 5.1.

If $X$ is a BGIG process, then:

  1. (1)

$X$ is a finite variation process,

  2. (2)

$X$ has infinite activity,

  3. (3)

$X$ is equal to the sum of its jumps, i.e.,

(77) $\forall t>0,\quad X_{t}=\sum_{s\leq t}\Delta X_{s}.$

###### Proof.

The Lévy measure $\pi$ satisfies $\pi({\mathbb{R}})=\infty$ and
$\int_{|x|<1}|x|\pi(\mathrm{d}x)<\infty$ because of the asymptotic for the
Jaeger integral obtained in lemma 2.1. The process $X$ is therefore of type B
in the sense of [44]. Then, points (1), (2) and (3) follow from [44, Th. 21.9,
21.3, 19.3] respectively. ∎

###  5.2. Blumenthal-Getoor index

The Blumenthal-Getoor index $\beta(X)$ (see [12]) of a Lévy process $X$
characterizes its small jumps. Formally, it is defined as:

(78)
$\beta(X):=\inf\left\\{p>0,\quad\int_{|x|<1}|x|^{p}\pi_{X}(\mathrm{d}x)<\infty\right\\},$

where $\pi_{X}$ denotes the Lévy measure of $X$. It is kown that a Lévy
process is of finite variation if and only if $\beta(X)<1$. For example, a
Tempered stable process has Blumenthal-Getoor index equal to
$\max(\beta_{+},\beta_{-})$ (see [35]). As consequence, a bilateral Gamma
process has a Blumenthal-Getoor index of 0.

######  Proposition 5.2.

If $X$ a BGIG process, we have $\beta(X)=1/2$.

###### Proof.

Let $p>0$; by definition of the BGIG Lévy measure (13), we can write:

(79)
$\int_{-1}^{1}|x|^{p}\pi(\mathrm{d}x)=\int_{0}^{1}\frac{e^{-a_{+}x/2}}{x^{1-p}}\left(\mathscr{I}(x/2b_{+},p_{+})+\max(p_{+},0)\right)\mathrm{d}x\\\
+\int_{0}^{1}\frac{e^{-a_{-}x/2}}{x^{1-p}}\left(\mathscr{I}(x/2b_{-},p_{-})+\max(p_{-},0)\right)\mathrm{d}x.$

By the asymptotic given in lemma 2.1, the two integrals in the right hand side
of (79) are finite if and only if $p>\max(1/2,0)=1/2$ and the desired results
follows. ∎

######  Proposition 5.3.

If X is a BGIG process, then we have the following almost-sure convergences:

(80) $\left\\{\begin{aligned} &\forall p>\frac{1}{2},\quad\lim_{t\to
0}t^{-1/p}X_{t}=0,\\\ &\forall p<\frac{1}{2},\quad\limsup_{t\to
0}t^{-1/p}X_{t}=\infty.\end{aligned}\right.$

###### Proof.

This is a direct consequence of proposition 5.2 and [12, Thm. 3.1, Thm. 3.3].
∎

###  5.3. Change of measure

We provide here some properties of BGIG distributions under change of
measures, which, as we will discuss in section 6, are important for financial
applications. Let us first define $\nu$, the Radon-Nykodym derivative of the
two Lévy measure $\pi$ and $\widetilde{\pi}$ of two BGIG distributions. We
have:

(81) $\forall
x\in{\mathbb{R}},\quad\nu(\mathrm{d}x):=e^{-(\widetilde{a}_{+}-a_{+})x/2}\left(\frac{\mathscr{I}(x/2\widetilde{b}_{+};\widetilde{p}_{+})+\max(\widetilde{p}_{+},0)}{\mathscr{I}(x/2b_{+};p_{+})+\max(p_{+},0)}\right)\mathbbm{1}_{(0,\infty)}(x)\mathrm{d}x\\\
+e^{-(\widetilde{a}_{-}-a_{-})|x|/2}\left(\frac{\mathscr{I}(|x|/2\widetilde{b}_{-};\widetilde{p}_{-})+\max(\widetilde{p}_{-},0)}{\mathscr{I}(|x|/2b_{-};p_{-})+\max(p_{-},0)}\right)\mathbbm{1}_{(0,\infty)}(x)\mathrm{d}x.\\\
$

######  Proposition 5.4.

Let $X$ be a $\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$-process and
let
$(\widetilde{a}_{+},\widetilde{b}_{+},\widetilde{p}_{+},\widetilde{a}_{-},\widetilde{b}_{-},$
$\widetilde{p}_{-})$ be six real numbers in the acceptable parameter space of
BGIG distributions. The following statements are equivalent:

  1. (1)

There exists a measure $\mathbb{Q}\overset{\mathrm{loc}}{\sim}\mathbb{P}$
under which $X$ is a BGIG process with parameters
$(\widetilde{a}_{+},\widetilde{b}_{+},\widetilde{p}_{+},\widetilde{a}_{-},\widetilde{b}_{-},\widetilde{p}_{-})$.

  2. (2)

$b_{\pm}=\widetilde{b}_{\pm}$ and $p_{\pm}=\widetilde{p}_{\pm}$.

###### Proof.

We use Theorem 33.1 of [44]. The only condition (the others are immediately
satisfied) to fullfill is:

(82)
$\int_{\mathbb{R}}\left(1-\sqrt{\nu(x)})\right)^{2}\pi(\mathrm{d}x)<\infty.$

Knowing $\pi$ and $\nu$ (see equations 13 and 81), this integral can be
rewritten as:

(83)
$\int_{0}^{\infty}x^{-1}\left(h(x;a_{+},b_{+},p_{+})-h(x;\widetilde{a}_{+},\widetilde{b}_{+},\widetilde{p}_{+}\right)^{2}\mathrm{d}x\\\
+\int_{0}^{\infty}x^{-1}\left(h(x;a_{-},b_{-},p_{-})-h(x;\widetilde{a}_{-},\widetilde{b}_{-},\widetilde{p}_{-}\right)^{2}\mathrm{d}x\\\
$

where $h(x;a,b,p):=e^{-ax/4}\sqrt{\mathscr{I}(x/2b;p)+\max(p,0)}$. The
integrability condition (82) is satisfied if and only if:

(84)
$\mathscr{I}(x/2b_{\pm},p_{\pm})+\max(p_{\pm},0)=\mathscr{I}(x/2\widetilde{b}_{\pm},\widetilde{p}_{\pm})+\max(\widetilde{p}_{\pm},0).$

for $x$ in a neighborhood of 0. The asymptotic for $x\to 0$ given in
proposition 2.1 gives that this is true if and only if:

(85) $b_{\pm}=\widetilde{b}_{\pm}\quad\mathrm{and}\quad
p_{\pm}=\widetilde{p}_{\pm}.$

∎

######  Remark 6.

As we will see in Section 6, this proposition is verified for the Esscher
transform; this will allow us to construct a risk-neutral version of the
(exponential) BGIG process, opening the way to the valuation of financial
instruments.

###  5.4. Jump structure of the sample paths

In this subsection we prove a result that can be useful for estimating the
parameters $b_{\pm}$ of a BGIG process from its sample path; this result is
particularly interesting if one wishes to infer BGIG parameters from
historical data featuring a very short time step (for instance, tick by tick
data). First, let us consider a process $X\sim\mathrm{GIG}(a,b,p)$ and let us
define, for $n\geq 1$ and $T>0$,

(86) $U_{n}:=\frac{1}{nT}\\#\left\\{t\leq T,\quad\Delta X_{t}\geq
n^{-2}\right\\}.$

######  Proposition 5.5.

The sequence of random variables $(U_{n})_{n\geq 1}$ converges to
$\sqrt{\frac{b}{2\pi}}$ almost surely.

###### Proof.

The random measure $\mu_{X}$ of the jumps of the BGIG process is a Poisson
random measure, whose intensity measure $\nu$ is given by (see [44]):

(87)
$\nu(\mathrm{d}t,\mathrm{d}x)=\frac{e^{-ax/2}}{x}\mathscr{I}(x/2b)\mathbbm{1}_{(0,+\infty)}(x)\mathrm{d}x\mathrm{d}t+\max(0,p)\frac{e^{-ax/2}}{x}\mathbbm{1}_{(0,+\infty)}(x)\mathrm{d}x\mathrm{d}t.$

We define $\widetilde{U}_{n}$ as:

(88) $\widetilde{U}_{n}=\frac{1}{T}\mu_{X}([0,T]\times[(n+1)^{-2},n^{-2}]).$

By definition, the $U_{n}$ are independent. Let us prove that the first two
moments of the $U_{n}$ are finite. By a change of variable $x=u^{-1/2}-n$ in
the expectation integral, we have:

(89)
$\mathbb{E}(\widetilde{U}_{n})=2\int_{0}^{1}\frac{e^{-a(x+n)^{-2}/2}}{(x+n)^{-2}}\mathscr{I}((x+n)^{-2}/2b)(x+n)^{-3}\mathrm{d}x\\\
\quad+2\max(0,p)\int_{0}^{1}\frac{e^{-a(x+n)^{-2}/2}}{(x+n)^{-2}}(x+n)^{-3}\mathrm{d}x.$

The second integral tends to $0$ when $n\to\infty$ by the dominated
convergence theorem. Furthermore, using the asymptotic of lemma 2.1, we have:

(90)
$\frac{\mathscr{I}((x+n)^{-2}/2b)}{(x+n)^{-1}}\underset{n\to\infty}{\longrightarrow}\sqrt{\frac{b}{2\pi}}.$

Again, by the dominated convergence theorem, we have:

(91)
$\int_{0}^{1}\frac{e^{-a(x+n)^{-2}/2}}{(x+n)^{-2}}\mathscr{I}((x+n)^{-2}/2b)(x+n)^{-3}\mathrm{d}x\underset{n\to\infty}{\longrightarrow}\sqrt{\frac{b}{2\pi}}.$

Thus, we have $\mathbb{E}(\widetilde{U}_{n})\to\sqrt{b/2\pi}$ and similarly
$\mathrm{Var}(\widetilde{U}_{n})\to(1/T)\sqrt{b/2\pi}$ when $n\to\infty$.
Consequently,

(92)
$\sum_{n=1}^{\infty}\frac{\mathrm{Var}(\widetilde{U}_{n})}{n^{2}}<\infty,$

This allows to apply Kolmogorov’s strong law of large numbers:

(93)
$\lim_{n\to\infty}U_{n}=\lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}\widetilde{U}_{n}+\lim_{n\to\infty}\frac{1}{nT}\mu_{X}([0,T]\times[1,\infty))=\sqrt{\frac{b}{2\pi}}$

since $\mu_{X}([0,T]\times[1,\infty))<\infty$. ∎

Now, let $X$ be a $\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$ process
and let us define the random variables:

(94) $\left\\{\begin{aligned} U_{n}^{+}&:=\frac{1}{nT}\\#\\{t\leq
T,\quad\Delta X_{t}\geq n^{-2}\\},\\\ U_{n}^{-}&:=\frac{1}{nT}\\#\\{t\leq
T,\quad\Delta X_{t}\leq-n^{-2}\\}.\end{aligned}\right.$

######  Proposition 5.6.

The following $\mathbb{P}$-almost sure convergence holds:

(95)
$\lim_{n\to\infty}(U_{n}^{+},U_{n}^{-})=\left(\sqrt{\frac{b_{+}}{2\pi}},\sqrt{\frac{b_{-}}{2\pi}}\right).$

###### Proof.

The proof is straightforward with proposition 5.5 and since a BGIG process can
be decomposed as $X=X^{+}-X^{-}$ with $X^{\pm}$ two independent GIG processes.
∎

###  5.5. Simulation of BGIG processes

For a stochastic process to gain popularity in quantitative finance (notably
from the point of view of market practitioners), it is important that it can
be efficiently simulated; the reason is that, in order to value very exotic
instruments (that can not be reached via closed or semi-closed formulas), one
typically rely on Monte-Carlo simulations (see Section 6.2.1), hence the
necessity of generating a possibly large number of sample paths. Let us
therefore provide two different methods for simulating BGIG processes.

####  5.5.1. Direct simulation

First, a simulation of the BGIG process at a collection of integer times
$\\{t_{i}\\}_{i=1\dots N}\in\mathbb{N}^{N}$ can be achieved very simply.
Indeed, since BGIG process are Lévy process and therefore have independent
increments, in this case we can write:

(96) $\forall i=1\dots N,\quad
X_{t_{i}}=\sum_{k=0}^{t_{i}}\left(X_{k}-X_{k-1}\right)\overset{\mathcal{L}}{=}\sum_{k=1}^{t_{i}}X_{1}^{(k)}$

where the $X_{1}^{(k)}$ are i.i.d. random variables following the generating
BGIG distribution. We illustrate this method in Figure 3; the theoretical
density are obtained by inverting the characteristic function in equation 75.

(a) $t=5$

(b) $t=20$

Figure 3. Histogram of 10 000 simulations of a BGIG(1,2,1,3,4,5)-process for
$t=5$ and $t=20$, and the corresponding theoretical densities.

####  5.5.2. Difference of GIG processes

Since BGIG distributions are not convolution-closed, this method cannot be
applied in general for $\\{t_{i}\\}_{i=1\dots N}\in{\mathbb{R}}_{+}^{N}$. In
this case, recalling that a BGIG process $X$ can be written $X:=X_{+}-X_{-}$
where $X_{\pm}$ are two independent GIG processes, we can use the simulation
algorithm developed in [27] for each GIG process to obtain a simulation for
the BGIG process $X$. In Figure 4, we display two set of simulations of BGIG
processes using this methodology.

(a) $\mathrm{BGIG}(2,1,-1,4,1,-3)$

(b) $\mathrm{BGIG}(2,1,-3,10,8,6)$

Figure 4. 20 simulations of two BGIG-processes for $t\in(0,1)$.

##  6\. The exp-BGIG stock market model

Let us now demonstrate how BGIG processes can be used in mathematical finance,
in particular for the purpose of asset returns modeling and derivatives
pricing.

###  6.1. Definition and change of measure

In continuous time finance, a popular approach to modeling asset prices is via
the introduction of exponential Lévy processes, because they allow to capture
realistic features (most notably the occurrence of jumps) while remaining
computationally tractable; see, among many other references, [8, 46, 17].
Here, we will assume that the exponential Lévy process is an exponential BGIG
(exp-BGIG) process; More precisely, we introduce the stochastic process
$S=\\{S_{t}\\}_{t\geq 0}$ describing the price of a financial asset (a stock
price, an index, an intereste rate…), that is defined as:

(97) $\forall t\geq 0,\quad S_{t}=S_{0}e^{X_{t}},\quad S_{0}>0,$

where $X=\\{X_{t}\\}_{t\geq 0}$ is a BGIG process. One essential point when
one wishes to conduct pricing of financial instruments written on $S$ is to
find a measure (called risk-neutral measure) under which the discounted asset
price process $\widetilde{S}:=\left(e^{-rt}S_{t}\right)_{t\geq 0}$ is a
martingale [20, 46]; here $r$ is the risk-free rate that is assumed to be
deterministic and continuously compounded. In order to find such measure, let
us start by mentioning the following equivalence:

######  Proposition 6.1.

Assume $a_{+}>2$ and let $S$ be an exp-BGIG process, we have the equivalence:

  1. (1)

$S$ is a martingale,

  2. (2)

the following equality holds:

(98)
$\left(\frac{a_{+}}{a_{+}-2}\right)^{p_{+}/2}\frac{K_{p_{+}}\left(\sqrt{b_{+}(a_{+}-2)}\right)}{K_{p_{+}}\left(\sqrt{b_{+}a_{+}}\right)}\times\left(\frac{a_{-}}{a_{-}+2}\right)^{p_{-}/2}\frac{K_{p_{-}}\left(\sqrt{b_{-}(a_{-}+2)}\right)}{K_{p_{-}}\left(\sqrt{b_{-}a_{-}}\right)}=1.$

###### Proof.

For exponential Lévy models, we have by [17, Prop. 3.17] that the martingale
condition is equivalent to have $\mathbb{E}[e^{X_{1}}]=1$. By noticing that
$\mathbb{E}[e^{X_{1}}]=\Phi(-\mathrm{i},1)$, the result follows. ∎

######  Remark 7.

Of course for the discounted price process $\tilde{S}$, the right hand side of
(98) is replaced by $e^{r}$.

Since Lévy markets are incomplete, there exists a potentially infinite number
of measures under which $\widetilde{S}$ is a martingale [46]. One classical
way to achieve martingality is via the so-called Esscher transform (see [25]),
that provides a straightforward procedure to build a risk-neutral measure.
Indeed, let us define $\theta^{*}$ to be the solution of the equation:

(99)
$\log\Phi_{X}(-\mathrm{i}(\theta^{*}+1))-\log\Phi_{X}(-\mathrm{i}\theta^{*})=r.$

where $\Phi_{X}$ is the characteristic function of $X$, and introduce the
Esscher transform $\mathbb{P}^{*}$ of the $\mathbb{P}$ measure defined by the
following Radon-Nikodym derivative

(100)
$\frac{d\mathbb{P}^{*}_{\mathcal{F}_{t}}}{d\mathbb{P}_{\mathcal{F}_{t}}}=\frac{e^{\theta^{*}X_{t}}}{\mathbb{E}\left[e^{\theta^{*}X_{t}}\right]}\,=\,e^{\theta^{*}X_{t}-\ln\Phi_{X}(-\mathrm{i}\theta^{*})t}.$

Then, following remark 7, $\tilde{S}$ is a martingale under $\mathbb{P}^{*}$.
We note that, as a direct consequence of definition 100, the characteristic
function of the process $X$ under $\mathbb{P}^{*}$ can be explicitly written
as

(101)
$\Phi_{X}^{*}(u,t)=\mathbb{E}^{\mathbb{P}^{*}}\left[e^{\mathrm{i}uX_{t}}\right]=\frac{\Phi_{X}(u-\mathrm{i}\theta^{*},t)}{\Phi_{X}(-\mathrm{i}\theta^{*},t)}.$

######  Proposition 6.2.

Let $S$ be an exp-BGIG process process with parameters
$(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$. Under the Esscher measure
$\mathbb{P}^{*}$, $S$ is still an exp-BGIG process with parameters
$(a_{+}-2\theta^{*},b_{+},p_{+},$
$a_{-}+2\theta^{*},b_{-},p_{-})$.

###### Proof.

By the properties of the convolution, it is sufficient to prove the result for
a one sided GIG distribution. By the time-scaling property of the Lévy
characteristic functions, it is sufficient to prove it for $t=1$. For $a,b>0$
and $p\in{\mathbb{R}}$, we have:

(102) $\forall
u\in{\mathbb{R}},\quad\Phi^{*}_{\mathrm{GIG}_{+}}(u)=\frac{\Phi_{\mathrm{GIG}_{+}}(u-\mathrm{i}\theta^{*},t)}{\Phi_{\mathrm{GIG}_{+}}(-\mathrm{i}\theta^{*},t)}=\left(\frac{a-2\theta^{*}}{a-2\theta^{*}-\mathrm{i}u}\right)^{p/2}\frac{K_{p}\left(\sqrt{b(a-2\theta^{*}-\mathrm{i}2u)}\right)}{K_{p}\left(\sqrt{b(a-2\theta^{*})}\right)}.$

Thus, under the Esscher measure, a $\mathrm{GIG}_{+}(a,b,p)$ distribution is a
$\mathrm{GIG}(a-2\theta^{*},b,p)$ distribution. The desired result immediately
follows. ∎

This proposition obviously satisfies the parameters condition of proposition
5.4 and the process thus obtained satisfies martingale condition of
proposition 6.1 if $a_{+}-2\theta^{*}>2$. We will see in Section 7.2 that this
condition is satisfied in practice.

###  6.2. Pricing Methods

Once a risk neutral measure $\mathbb{P}^{*}$ has been constructed, it is
possible to apply the risk-neutral pricing theory (see [20] or any classical
monograph on quantitative finance) to compute the value of derivatives, that
is, of financial instruments written on the the underlying asset price $S$.
Such instruments (typically, options) are characterized by the payoff $g(S)$
they deliver, that can be a function of either the final underlying price
$S_{T}$ or of the whole realization of the process on $[0,T]$; their price at
time $t\in[0,T]$ corresponds to the risk-neutral expectation of the discounted
payoff, that is:

(103) $\forall t\geq 0,\quad
V_{t}=\mathbb{E}^{\mathbb{P}^{*}}[e^{-r(T-t)}g(S)|\mathcal{F}_{t}].$

In the following, we will conduct valuations on $t=0$, and we will focus on
European call options, which are among the most widely traded derivatives and
whose payoff is given by

(104) $g_{\mathrm{eur}}(S)=\max(S_{T}-K,0)=(S_{T}-K)^{+}$

for some predetermined $K>0$, called strike price. We will denote their price
by $C(S_{0},K,T)$ to highlight the dependence on current (spot) asset value
$S_{0}$, strike $K$ and maturity (or expiry) $T$. Note that the reciprocal of
a European call is a European put, which pays $(K-S_{T})^{+}$ at maturity.

There exists several methods to evaluate the expectation (103); we briefly
describe two popular methods that we will implement, namely the Monte-Carlo
approach and the Fourier pricing approach.

####  6.2.1. Monte-Carlo pricing

We can approximate option prices by introducing

(105)
$C^{(n)}(S_{0},K,T)=\frac{e^{-rT}}{n}\sum_{k=1}^{n}g_{\mathrm{eur}}(S_{k})$

where the $S_{k}$ are i.i.d. BGIG processes that are simulated as difference
of two risk-neutral GIG processes (see section 5.5). By the strong law of
large numbers,

(106) $C^{(n)}(S_{0},K,T)\underset{n\to\infty}{\longrightarrow}C(S_{0},K,T)$

alsmot surely. For a complete introduction to Monte-Carlo methods in finance,
see [26].

####  6.2.2. Fourier pricing

While (103) can be carried out via a numerical integration of the payoff
function over the density function, this method is generally not practical in
the case of exponential Lévy processes because the densities involve special
functions or because the processes are not convolution closed. However, the
characteristic functions of such processes are generally known in a simpler
form, and several methods exist that allow to directly express European option
prices in terms of integrals over the characteristic function. For instance,
one can use the Lewis-Lipton formula [38]:

(107)
$C(S_{0},K,T)=\frac{-Ke^{-rT}}{2\pi}\int_{\mathrm{i}v_{1}-\infty}^{\mathrm{i}v_{1}+\infty}\frac{e^{-\mathrm{i}zk}\Phi^{*}(-z,T)}{z^{2}-\mathrm{i}z}\mathrm{d}z,$

where $k=\ln(S_{0}/K)+rT$, $v_{1}>1$ and $\Phi^{*}$ is the BGIG characteristic
function under $\mathbb{P}^{*}$. Note that many extensions of this idea exist
that improve its accuracy and performance, such as the Carr-Madan formula [16]
or the COS [22], CONV [39] or PROJ [31] methods.

##  7\. Calibration procedure and numerical examples

We terminate the paper by performing a calibration of the exp-BGIG model (97)
over real market data, and by pricing European call and put options under both
Monte Carlo and Fourier methods.

###  7.1. Calibration methodology

In this subsection, we provide a procedure to calibrate the exp-BGIG model
over empirical data, that is to infer the parameters
$(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-}$) from a given data set composed of
daily close prices, for a set of trading dates $(t_{1},t_{2},\dots,t_{n})$
with for simplicity constant spacing $t_{k+1}-t_{k}=\delta>0$. We denote the
corresponding underlying prices by $S_{k}=S_{t_{k}}$ and we define the series
of log-returns as $L=(\ln(S_{k+1}/S_{k}))_{k=1\dots n}$. Let us remark that

(108)
$\ln\left(\frac{S_{k+1}}{S_{k}}\right)=X_{t_{k+1}}-X_{t_{k}}\overset{\mathcal{L}}{=}X_{\delta}$

where the second equality (in law) is a consequence of the stationarity of
increments of Lévy processes. In particular when $\delta=1$, $L$ is a
collection of $n$ independent
$\mathrm{BGIG}(a_{+},b_{+},p_{+},a_{-},b_{-},p_{-})$ random variables.

Our calibration procedure can be divided into two steps: first we estimate
$a_{+}$ and $a_{-}$ thanks to the asymptotic properties of BGIG distributions
(proposition 4.7), then we estimate $b_{+}$, $b_{-}$, $p_{+}$ and $p_{-}$ via
the method of moments.

####  7.1.1. Estimating $a_{\pm}$

Let us denote $m_{n}=\min_{1\leq i\leq n}L_{i}$ and $M_{n}:=\max_{1\leq i\leq
n}L_{i}$. Taking $\delta=1$ and using proposition 4.7 we have:

(109) $m_{n}\underset{n\to\infty}{\sim}-\frac{2}{a_{-}}\ln
n\quad\mathrm{and}\quad M_{n}\underset{n\to\infty}{\sim}\frac{2}{a_{+}}\ln n.$

Our estimators $\widehat{a}_{+}$ andd $\widehat{a}_{-}$ will therefore be
simply defined as:

(110) $\widehat{a}_{+}:=\frac{2\ln
n}{M_{n}}\quad\mathrm{and}\quad\widehat{a}_{-}:=-\frac{2\ln n}{m_{n}}.$

These estimators are directly linked to the nature of the BGIG distribution
and thus have the advantage to be robust. Moreover they allow to restrict the
subsequent method of moments to the estimation of four parameters instead of
six, which reduces the complexity of the system to be solved.

####  7.1.2. Estimating $b_{\pm}$ and $p_{\pm}$

The parameters $b_{\pm}$ and $p_{\pm}$ are estimated thanks to the relations
between moments and cumulants. For $\delta>0$, all the returns $L_{k}$ are
i.i.d and we have:

(111) $\left\\{\begin{aligned}
&\displaystyle\mathbb{E}(L_{k})=\delta\kappa_{1},\\\
&\displaystyle\mathrm{Var}(L_{k})=\delta\kappa_{2},\\\
&\displaystyle\gamma_{1}(L_{k})=\frac{\kappa_{3}}{\kappa_{2}^{3/2}\sqrt{\delta}},\\\
&\displaystyle\gamma_{2}(L_{k})=3+\frac{\kappa_{4}}{\kappa_{2}^{2}\delta}.\end{aligned}\right.$

where $\gamma_{1}(L_{k})$ (rep. $\gamma_{2}(L_{k})$) denotes the skewness
(resp. kurtosis) of $L_{k}$, and where the BGIG cumulants $\kappa_{j}$ are the
ones given in section 2.3.3. The estimators $\mathfrak{m}_{j}$ of these four
central and standardized moments are defined as:

(112) $\left\\{\begin{aligned}
\displaystyle\mathfrak{m}_{1}&=\frac{1}{n}\sum_{k=0}^{n}L_{k},\\\
\displaystyle\mathfrak{m}_{2}&=\frac{1}{n}\sum_{k=0}^{n}(L_{k}-\mathfrak{m}_{1})^{2},\\\
\displaystyle\mathfrak{m}_{3}&=\frac{1}{n\mathfrak{m}_{2}^{3/2}}\sum_{k=0}^{n}(L_{k}-\mathfrak{m}_{1})^{3},\\\
\displaystyle\mathfrak{m}_{4}&=\frac{1}{n\mathfrak{m}_{2}^{2}}\sum_{k=0}^{n}(L_{k}-\mathfrak{m}_{1})^{4}.\end{aligned}\right.$

Cumulants in equation 111 are functions of the parameters $a_{\pm}$, $b_{\pm}$
and $p_{\pm}$: we thus look for a vector
$(b_{+},p_{+},b_{-},p_{-})\in((0,+\infty)\times{\mathbb{R}})^{2}$ solving:

(113) $\left\\{\begin{aligned}
\displaystyle\mathfrak{m}_{1}&=\delta\kappa_{1}(\widehat{a}_{+},b_{+},p_{+},\widehat{a}_{-},b_{-},p_{-}),\\\
\displaystyle\mathfrak{m}_{2}&=\delta\kappa_{2}(\widehat{a}_{+},b_{+},p_{+},\widehat{a}_{-},b_{-},p_{-}),\\\
\displaystyle\mathfrak{m}_{3}&=\frac{\kappa_{3}(\widehat{a}_{+},b_{+},p_{+},\widehat{a}_{-},b_{-},p_{-})}{\kappa_{2}(\widehat{a}_{+},b_{+},p_{+},\widehat{a}_{-},b_{-},p_{-})^{3/2}\sqrt{\delta}},\\\
\displaystyle\mathfrak{m}_{4}&=3+\frac{\kappa_{4}(\widehat{a}_{+},b_{+},p_{+},\widehat{a}_{-},b_{-},p_{-})}{\kappa_{2}(\widehat{a}_{+},b_{+},p_{+},\widehat{a}_{-},b_{-},p_{-})^{2}\delta},\end{aligned}\right.$

where $\widehat{a}_{+}$ and $\widehat{a}_{-}$ are the estimated parameters
already obtained in Section 7.1.1. In practice, we approach a solution with a
least-squares minimization and we denote
$(\widehat{b}_{+},\widehat{p}_{+},\widehat{b}_{-},\widehat{p}_{-})$ the
approximated solution.

####  7.1.3. Risk-neutral parameters

Sections 7.1.1 and 7.1.2 allow to estimate the parameters of an exp-BGIG
process; however these parameters are historical parameters and therefore
cannot be used directly for pricing purposes. The last step is to transform
them into risk-neutral parameters via the Esscher transform. First, we have
determine $\theta^{*}$ the solution of the Esscher equation 99; this can be
achieved by any root-finding method on ${\mathbb{R}}$. Once $\theta^{*}$ is
found, the risk-neutral parameters
$(\widehat{a}_{+}^{*},\widehat{b}_{+}^{*},\widehat{p}_{+}^{*},\widehat{a}_{-}^{*},\widehat{b}_{-}^{*},\widehat{p}_{-}^{*})$
are easily obtained thanks to proposition 6.2:

(114) $\widehat{a}_{\pm}^{*}=\widehat{a}_{\pm}\mp
2\theta^{*},\quad\widehat{b}_{\pm}^{*}=\widehat{b}_{\pm}\quad\textrm{and}\quad\widehat{p}_{\pm}^{*}=\widehat{p}_{\pm}.$

###  7.2. Calibration on S&P data

We apply the calibration methodology described in subsection 7.1 to daily
close prices of the S&P500 index, for a 2 years trading period (2021-01-01 to
2024-01-01) with $\delta=1$ (daily returns). To avoid conidering outliers, we
remove the 1% highest and lowest log-returns from our dataset. The obtained
parameters for the daily BGIG distributions are:

(115) $\left\\{\begin{aligned} \widehat{a}_{+}&=5.58753&&\times 10^{2},\\\
\widehat{b}_{+}&=4.43139&&\times 10^{-2},\\\ \widehat{p}_{+}&=2.53084&&\times
10^{0},\\\ \widehat{a}_{-}&=4.39902&&\times 10^{2},\\\
\widehat{b}_{-}&=2.42973&&\times 10^{-2},\\\ \widehat{p}_{-}&=2.26669&&\times
10^{0}.\end{aligned}\right.$

In Figure 5, we display in a same graph the histograms of the log-returns of
S&P 500 index as well as the calibrated density. The corresponding moments and
the calibration errors are given in table 1. Lets us observe that the kurtosis
is strictly bigger than $3$, confirming the presence of fat tails, and that
the distribution is skewed to the left (negative third moment), indicating a
prevalence of downward moves over upward ones.

Moments Estimation and Errors
---
| Moments | Errors
Order |  Empirical |  Estimated |  Absolute |  Percentage
1 |  1.38234e-04 |  1.38194e-04 |  4.02246e-08 |  2.91074e-04
2 |  9.29722e-05 |  9.29885e-05 |  1.62870e-08 |  1.75151e-04
3 |  -2.46379e-01 |  -2.46271e-01 |  1.07698e-04 |  4.37315e-04
4 |  3.16266e+00 |  3.16231e+00 |  3.58132e-04 |  1.13250e-04
Table 1. Estimated and empirical moments (up to order 4), and the associated
errors. Figure 5. Historical log-returns of the S&P 500 index and the
calibrated BGIG density (parameters from (115)).

###  7.3. Option pricing

Let us terminate by pricing some financial instruments using the calibration
results from section 7.2. We solve the Esscher equation (99) numerically and
we find

(116) $\theta^{*}=-1.98436$

and, using (115) and (114), we immediately deduce the risk-neutral parameters.
We can then price European call and put options with either Monte-Carlo (105)
or Lewis-Lipton (107) methods. Results can be found in Table 2, for different
strike prices (from deep out-of-the money to deep in-the-money configurations)
and a maturity of one year (252 trading days).

European Option Prices
---
Option | Call | Put
$K$ |  Lewis-Lipton |  Monte-Carlo |  Lewis-Lipton |  Monte-Carlo
$0.5\times S_{0}$ |  5.00000e-01 |  5.00098e-01 |  9.25269e-08 |  4.00395e-07
$0.8\times S_{0}$ |  2.04505e-01 |  2.04489e-01 |  4.50541e-03 |  4.61485e-03
$1\times S_{0}$ |  6.10931e-02 |  6.10825e-02 |  6.10931e-02 |  6.10295e-02
$1.2\times S_{0}$ |  9.52017e-03 |  9.45246e-03 |  2.09520e-01 |  2.08716e-01
$1.5\times S_{0}$ |  2.25386e-04 |  2.40956e-04 |  5.00225e-01 |  5.00270e-01
Table 2. Prices of European calls/puts with $S_{0}=1$, $T=252$ days, $r=0\%$,
and $n=50~{}000$ paths for the Monte-Carlo (MC) method. MC prices are more
accurate for put options since the payoff is bounded.

##  8\. Conclusion

We have introduced and studied BGIG distributions and processes. These
distributions, constructed by convolution of GIG random variables, are proved
to be smooth, real-analytic on ${\mathbb{R}}\backslash\\{0\\}$ and unimodal.
Since they are also infinitely divisible, they generate associated Lévy
processes, that we have shown have finite variations, infinite activity and a
Blumenthal-Getoor index of $1/2$. With the help of known algorithms, these
processes can also be efficiently simulated. Moreover they are stable under
Esscher transform, which makes them suitable for financial modeling and option
pricing.

As an application of this theory, we have defined a stock market model based
on exponential BGIG processes and provided a calibration methodology relying
on the asymptotics for extrema of BGIG distributions and on the method of
moments; we have tested this methodology over historical stock market returns,
and, combined with the estimation of the Esscher parameter, we have applied
the model to the pricing of S&P 500 index options via Monte-Carlo and Fourier
based methods.

Future work should include calibration of the model over option data,
investigation of the BGIG implied volatility smile and comparison with the
market smile, and pricing of various generations of exotic options using
state-of-the-art pricing techniques. The existence of a closed or a semi-
closed (e.g. via series expansion) pricing formula for European-style options
would also be an interesting question to address.

##  9\. Acknowledgement

We thank Iosif Pinelis for his help proving proposition 4.7.

## References

  * [1] M. Abramowitz and I. Stegun.  Handbook of Mathematical Functions.  Dover Publications, Mineola, NY, 1972.
  * [2] J. Ph. Aguilar.  Some pricing tools for the Variance Gamma model.  International Journal of Theoretical and Applied Finance, 23(4):2050025, 2020.
  * [3] J. Ph. Aguilar and J. L. Kirkby.  Robust and nearly exact option pricing with bilateral Gamma processes.  The Journal of Derivatives, 30(1):8–31, 2022.
  * [4] J. Ph. Aguilar and J. L. Kirkby.  Closed form option pricing for exponential Lévy models: a residue approach.  Quantitative Finance, 23(2):251–278, 2023.
  * [5] A. Baricz, T. Pogány, S. Ponnusamy, and I. Rudas.  Bounds for Jaeger integrals.  Journal of Mathematical Chemistry, 53, 2015.
  * [6] O. E Barndorff-Nielsen.  Exponentially decreasing distributions for the logarithm of particle size.  Proceedings of the Royal Society of London, 353:401–419, 1977.
  * [7] O. E. Barndorff-Nielsen.  Normal inverse Gaussian distributions and stochastic volatility models.  Scandinavian Journal of Statistics, 24(1):1–133, 1997.
  * [8] O. E. Barndorff-Nielsen, T. Mikoschand, and Resnick S.I.  Lévy Processes: Theory and Applications.  Birkhäuser Boston, 2001.
  * [9] O. E. Barndorff-Nielsen and N. Shephard.  Non-Gaussian Ornstein-Uhlenbeck-based models and some of their uses in financial economics.  J. R. Stat. Soc. Ser. B Stat. Methodol., 63(2):167–241, 2001.
  * [10] O.E. Barndorff-Nielsen and C. Halgreen.  Infinite divisibility of the hyperbolic and generalized inverse gaussian distributions.  Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 38(4):309–311, 1977.
  * [11] F. Black and M. Scholes.  The pricing of options and corporate liabilities.  Journal of Political Economy, 81(3):637–654, 1973.
  * [12] R. M. Blumenthal and R. K. Getoor.  Sample functions of stochastic processes with stationary independent increments.  Journal of Mathematics and Mechanics, 10(3):493–516, 1961.
  * [13] S. I. Boyarchenko and S. Z. Levendorskii.  Option pricing for truncated Lévy processes.  International Journal of Theoretical and Applied Finance, 3(3):549–552, 2000.
  * [14] S. I. Boyarchenko and S. Z. Levendorskii.  Non-Gaussian Merton-Black-Scholes Theory.  World Scientific Publishing Co., River Edge, NJ, 2002.
  * [15] P. Carr, H. Geman, D. Madan, and M. Yor.  The fine structure of asset returns: An empirical investigation.  The Journal of Business, 75:305–332, 2002.
  * [16] P. Carr, H. Geman, D. Madan, and M. Yor.  The fine structure of asset returns: An empirical investigation.  Journal of Business, 75(2):305–332, 2002.
  * [17] R. Cont and P. Tankov.  Financial Modelling with Jump Processes.  Chapman & Hall, New York, 2004.
  * [18] NIST Digital Library of Mathematical Functions, 2024.  F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * [19] E. Eberlein and E.-A. v. Hammerstein.  Generalized hyperbolic and inverse Gaussian distributions: limiting cases and approximation of processes.  Progress in Probability, 58:221–264, 2004.
  * [20] E. Eberlein and J. Kallsen.  Mathematical Finance.  Springer Cham, Switzerland, 2019.
  * [21] E. Eberlein and K. Prause.  The Generalized Hyperbolic Model: Financial Derivatives and Risk Measures. In: Geman, H., Madan, D., Pliska, S.R., Vorst, T. (eds) Mathematical Finance — Bachelier Congress 2000.  Springer, Berlin, Heidelberg, 2002.
  * [22] F. Fang and C. W. Oosterlee.  A novel pricing method for European options based on Fourier-cosine series expansions.  SIAM Journal on Scientific Computing, 31(2):826–848, 2009.
  * [23] P. Flajolet, X. Gourdon, and P. Dumas.  Mellin transforms and asymptotics: Harmonic sums.  Theoretical Computer Science, 144:3–58, 1995.
  * [24] P. Freitas.  Sharp bounds for the modulus and phase of Hankel functions with applications to Jaeger integrals.  Mathematics of Computation, 87(309):289–308, 2018.
  * [25] H. U. Gerber and E. S. W. Shiu.  Option pricing by Esscher transforms.  Insurance Mathematics & Economics, 3:287, 1995.
  * [26] P. Glasserman.  Monte Carlo Methods in Financial Engineering.  Springer New York, NY, 2003.
  * [27] S.J. Godsill and Y. Kındap.  Point process simulation of generalised inverse Gaussian processes and estimation of the Jaeger integral.  Statistics and Computing, 32:13, 2022.
  * [28] L. Hörmander.  The Analysis of Linear Partial Differential Operators I: Distribution Theory and Fourier Analysis.  Springer Berlin, Heidelberg, 2003.
  * [29] J.C. Jaeger and M. Clarke.  A short table of I (O, I; x).  Proc. R. Soc. Edinburgh A, 61:229–230, 1942.
  * [30] B. Jorgensen.  Statistical Properties of the Generalized Inverse Gaussian Distribution.  Lecture Notes in Statistics. Springer New York, 2012.
  * [31] J.L. Kirkby.  Efficient option pricing by frame duality with the fast Fourier transform.  SIAM Journal on Financial Mathematics, 6(1):713–747, 2015.
  * [32] I. T. Koponen.  Analytic approach to the problem of convergence of truncated Lévy flights towards the Gaussian stochastic process.  Physical review. E, Statistical physics, plasmas, fluids, and related interdisciplinary topics, 52 1:1197–1199, 1995.
  * [33] S. G. Krantz and H. R. Parks.  A Primer of Real Analytic Functions.  Birkhäuser Boston, MA, 2002.
  * [34] U. Küchler and S. Tappe.  Option pricing in bilateral Gamma stock models.  Statistics & Decisions, 27(4):281–307, 2009.
  * [35] U. Küchler and S. Tappe.  Tempered stable distributions and processes.  Stochastic Processes and their Applications, 123(12):4256–4293, 2013.
  * [36] U. Küchler and S. Tappe.  Bilateral Gamma distributions and processes in financial mathematics.  Stochastic Processes and their Applications, 118(2):261–283, 2008.
  * [37] U. Küchler and S. Tappe.  On the shapes of bilateral Gamma densities.  Statistics & Probability Letters, 78(15):2478–2484, 2008.
  * [38] A.L. Lewis.  A simple option formula for general jump-diffusion and other exponential Lévy processes.  Available at SSRN: https://ssrn.com/abstract=282110, 2001.
  * [39] R. Lord, F. Fang, F. Bervoets, and C.W. Oosterlee.  A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes.  SIAM J. Sci. Comput., 10:1678–1705, 2008.
  * [40] Y. L. Luke.  Integrals of Bessel functions.  Dover Publications, Mineola, NY, 2014.
  * [41] D. Madan, P. Carr, and E. Chang.  The Variance Gamma process and option pricing.  European Finance Review, 2:79–105, 1998.
  * [42] D. Madan and E. Seneta.  The Variance Gamma (VG) model for share market returns.  The Journal of Business, 63:511–24, 1990.
  * [43] R. N. Mantegna and H. E. Stanley.  Stochastic process with ultraslow convergence to a Gaussian: The truncated Lévy flight.  Phys. Rev. Lett., 73:2946–2949, 1994.
  * [44] K. I. Satō.  Lévy Processes and Infinitely Divisible Distributions.  Cambridge studies in advanced mathematics. Cambridge University Press, 1999.
  * [45] K. I. Satō and M. Yamazato.  On distribution functions of class L.  Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, 43:273–308, 1978.
  * [46] W. Schoutens.  Lévy processes in finance: pricing financial derivatives.  John Wiley & Sons, Hoboken, NJ, 2003.
  * [47] M. Zhang, M. Revie, and J. Quigley.  Saddlepoint approximation for the generalized inverse gaussian lévy process.  Journal of Computational and Applied Mathematics, 411:114275, 2022.