# Pricing European Options under Stochastic Volatility Models:
Case of five-Parameter Variance-Gamma Process

A.H.Nzokem1

###### Abstract.

The paper builds a Variance-Gamma (VG) model with five parameters: location
($\mu$), symmetry ($\delta$), volatility ($\sigma$), shape ($\alpha$), and
scale ($\theta$); and studies its application to the pricing of European
options. The results of our analysis show that the five-parameter VG model is
a stochastic volatility model with a $\Gamma(\alpha,\theta)$
Ornstein–Uhlenbeck type process; the associated Lévy density of the VG model
is a KoBoL family of order $\nu=0$, intensity $\alpha$, and steepness
parameters
$\frac{\delta}{\sigma^{2}}-\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}$
and
$\frac{\delta}{\sigma^{2}}+\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}$;
and the VG process converges asymptotically in distribution to a Lévy process
driven by a normal distribution with mean $(\mu+\alpha\theta\delta)$ and
variance $\alpha(\theta^{2}\delta^{2}+\sigma^{2}\theta)$. The data used for
empirical analysis were obtained by fitting the five-parameter Variance-Gamma
(VG) model to the underlying distribution of the daily SPY ETF data. Regarding
the application of the five-parameter VG model, the twelve-point rule
Composite Newton–Cotes Quadrature and Fractional Fast Fourier (FRFT)
algorithms were implemented to compute the European option price. Compared to
the Black–Scholes (BS) model, empirical evidence shows that the VG option
price is underpriced for out-of-the-money (OTM) options and overpriced for in-
the-money (ITM) options. Both models produce almost the same option pricing
results for deep out-of-the-money (OTM) and deep-in-the-money (ITM) options..

###### Key words and phrases:

stochastic volatility, Lévy process, Ornstein-Uhlenbeck process, infinitely
divisible distributions, Variance-Gamma (VG) model, function characteristic,
Esscher transform

1<EMAIL_ADDRESS>
##  1\. Introduction

Black-Scholes (BS) model [12] is considered the cornerstone of option pricing
theory. The model relies on the fundamental assumption that the asset returns
have a normal distribution with a known mean and variance. However, based on
empirical studies, the Black-Scholes (BS) model is inconsistent with a set of
well-established stylized features [15]. Due to the subsequent development of
the option pricing theory, a new class of models has emerged in the literature
to address the stylized characteristics of the markets. The probabilistic
property of infinity divisibility is the main characteristic of these new
models, and they belong to the family of Levy processes [23].
The new class of models can be divided into two subclasses of Jump-Diffusion
models and Stochastic Volatility models. The Jump-Diffusion process is modeled
as an independent Brownian motion plus a Compound Poisson Process. The popular
models in the literature are Merton’s jump-diffusion model [27] and Kou’s
jump-diffusion model [22], where logarithmic jump size follows a normal
distribution and an asymmetric double exponential distribution respectively.
The stochastic volatility (SV) model is another extension of the standard
geometric Brownian motion (GBM) model, where the observed volatility is
modeled as a stochastic process. In a stochastic volatility framework [2], the
constant volatility ($\sigma$) in a standard geometric Brownian motion (GBM)
model is replaced by a deterministic function of a stochastic process
($\sigma(Y_{t})$) where $Y_{t}$ represents the solution of the stochastic
differential equation (SDE). We have two main types of SV models in the
literature: Diffusion based SV models and non-Gaussian Ornstein-Uhlenbeck-
based SV models. In the popular diffusion-based SV models, $Y_{t}$ follows a
Cox-Ingersoll-Ross(CIR) process [17] or a Log-normal process [18] and the
deterministic function is a squared root of the stochastic process
($\sigma(Y_{t})=\sqrt{Y_{t}}$). The non-Gaussian Ornstein-Uhlenbeck-based SV
models have been introduced and thoroughly studied in [10, 11, 6, 7]. The SV
model with the Ornstein-Uhlenbeck type process is mathematically tractable and
has many appealing features.
From the perspective of derivative asset analysis, we will build a five-
parameter VG model as a Stochastic volatility model with
$\Gamma(\alpha,\theta)$ Ornstein-Uhlenbeck type process. While there is a
great number of studies on option pricing under the VG Model, most of the VG
Model in the literature has three parameters [24, 28, 25, 1], certainly due to
technical issues of fitting a high parametric model to the marginal
distribution of asset returns. The amount of literature considering the VG
Model with five parameters is rather limited. Using the five-parameter
Variance-Gamma model as an underlying distribution of the European option will
allow controlling both the excess kurtosis and the skewness of the underlying
data. In the option pricing theory, the challenge is often the existence of
the Equivalent Martingale Measure (EMM) and whether it preserves the structure
of the Variance-Gamma measure. The Variance-Gamma (VG) process is not a
Gaussian process, and the market is incomplete; therefore, the Equivalent
Martingale Measure is not unique. The Esscher transform of the historical
measure is considered optimal with respect to some optimization criterion
[13]. The Esscher martingale measure was shown [3] to coincide with the
minimal entropy martingale measure for Lévy processes.
The remainder of this paper is organized as follows. Section 2 is devoted to
building a five-parameter VG process and presenting parameter estimations and
simulations of the VG process. Section 3 investigates the Lévy density and the
asymptotic distribution of the VG process. And section 4 extends the Black-
Scholes framework, provides the integral representation for the option price,
and computes the VG option price numerically.

##  2\. Variance - Gamma Process: Stochastic Volatility Model

###  2.1. Lévy Framework and Asset Pricing

Let $(\Omega,\mathscr{F},\left\\{\mathscr{F}_{t}\right\\}_{t\geq
0},\mathbf{P})$ be a filtered probability space, with
$\mathscr{F}=\vee_{t>0}\mathscr{F}_{t}$ and
$\left\\{\mathscr{F}_{t}\right\\}_{t\geq 0}$ is a filtration.
$\mathscr{F}_{t}$ is a $\sigma$-algebra included in $\mathscr{F}$ and for
$\tau<t$, $\mathscr{F}_{\tau}\subseteq\mathscr{F}_{t}$.
A stochastic process $Y=\\{Y_{t}\\}_{t\geq 0}$ is a Lévy process, if it has
the following properties
(L1): $Y_{0}=0$ a.s;
(L2): $Y_{t}$ has independent increments, that is, for
$0<t_{1}<t_{2}<\dots<t_{n}$, the random variables $Y_{t_{1}}$,
$Y_{t_{1}}-Y_{t_{2}}$, …, $Y_{t_{n}}-Y_{t_{n-1}}$ are independent;
(L3): $Y_{t}$ has stationary increments, that is, for any
$t_{1}<t_{2}<+\infty$, the probability distribution of $Y_{t_{1}}-Y_{t_{2}}$
depends only on $t_{1}-t_{2}$ ;
(L4): $Y_{t}$ is stochastically continuous: for any $t$ and $\epsilon>0$,
$\lim_{s\to t}P\left(|Y_{s}-Y_{t}|>\epsilon\right)=0$;
(L5): $c\grave{a}dl\grave{a}g$ paths, that is, $t\mapsto Y_{t}$ is a.s right
continuous with left limits

Given a Lévy process $Y=\\{Y_{t}\\}_{t\geq 0}$ on the filtered probability
space $(\Omega,\mathscr{F},\left\\{\mathscr{F}_{t}\right\\}_{t\geq
0},\mathbf{P})$, we define the asset value process $S=\\{S_{t}\\}_{t\geq 0}$
such as $S_{t}=S_{0}e^{Y_{t}}$.

###### Theorem 2.1

(Lévy-Khintchine representation)
Let $Y=\\{Y_{t}\\}_{t\geq 0}$ be a Lévy process on $\mathbb{R}$. Then the
characteristic exponent admits the following representation.

$\displaystyle\varphi(\xi)=
-Log\left(Ee^{iY_{1}\xi}\right)=-i\gamma\xi+\frac{1}{2}\sigma^{2}\xi^{2}+\int_{0}^{t}\left(e^{i\xi
y}-1-y\xi 1_{|y|\leq 1}\right)\Pi(dy) $ (2.1)

where $\gamma\in\mathbb{R}$, $\sigma\geq 0$ and $\Pi$ is a $\sigma$-finite
measure called the Lévy measure of Y, satisfying the property

$\int_{-\infty}^{+\infty}Min(1,|y|^{2})\Pi(dy)<+\infty$

For the Theorem-proof, see [4, 20, 36]
Each Lévy process is uniquely determined by the Lévy–Khintchine triplet
$(\gamma,\sigma^{2},\Pi)$. The terms of this triplet suggest that a Lévy
process can be seen as having three independent components: a linear drift, a
Brownian motion, and a Lévy jump process. When the diffusion term $\sigma=0$,
we have a Lévy jump process; in addition, if $\gamma=0$, we have a pure jump
process.

###  2.2. $\Gamma(\alpha,\theta)$ Ornstein-Uhlenbeck proces

The Ornstein-Uhlenbeck process is a diffusion process introduced by Ornstein
and Uhlenbeck [37] to model the stochastic behavior of the velocity of a
particle undergoing Brownian motion. The Ornstein-Uhlenbeck diffusion
$\sigma^{2}=\\{\sigma^{2}(t),t\geq 0\\}$ is the solution of the Langevin
Stochastic Differential Equation (SDE) (2.2)

$\displaystyle d\sigma^{2}(t)=-\lambda\sigma^{2}(t)dt+dB(\lambda t). $ (2.2)

where $\lambda>0$ and $B=\\{B_{t},t\geq 0\\}$ is a Brownian motion. In recent
years, the Ornstein-Uhlenbeck process has been used in finance to capture
important distributional deviations from Gaussianity and to model dependence
structures. The extension of the Ornstein-Uhlenbeck processes was obtained by
replacing the Brownian motion in (2.2) by z(t), a background driving Lévy
process (BDLP) [11, 8, 10]. The SDE (2.2) becomes

$\displaystyle d\sigma^{2}(t)=-\lambda\sigma^{2}(t)dt+dz(\lambda t)
\quad\lambda>0\. $ (2.3)

where the process $z(t)=\\{z(t),t\geq 0,z(0)=0\\}$ is subordinator; a process
with non-negative, independent and stationary increments, which implies
$\sigma^{2}(t)\geq 0$. Correspondingly z(t) moves up entirely by jumps and
then tails off exponentially [7].

###### Lemma 2.2

The general form of the stationary process $\sigma^{2}(t)$ , solution of (2.3)
is given by :

$\displaystyle\sigma^{2}(t)=-\int_{0}^{+\infty}e^{-s}dz(\lambda t-s)
\quad\lambda>0\. $ (2.4)

###### Proof.

$\displaystyle\sigma^{2}(t)=-\int_{0}^{+\infty}e^{-s}dz(\lambda
t-s)=-\int_{0}^{+\infty}e^{-\lambda
s}dz(\lambda(t-s))=\int_{-\infty}^{t}e^{-\lambda(t-s)}dz(\lambda s) $ (2.5)

By using the variable changing method, we can have different expressions of
(2.4).

$\displaystyle\sigma^{2}(t)=\int_{-\infty}^{t}e^{-\lambda(t-s)}dz(\lambda s)
\quad\Longrightarrow\quad d\sigma^{2}(t)=-\lambda\sigma^{2}(t)dt+dz(\lambda t)
$ (2.6)

∎

Expression (2.5) can be written as follows:

$\displaystyle\sigma^{2}(t)=e^{-\lambda
t}\sigma^{2}(0)+\int_{0}^{t}e^{-\lambda(t-s)}dz(\lambda
s)\quad\quad\sigma^{2}(0)=\int_{-\infty}^{0}e^{\lambda s}dz(\lambda s) $ (2.7)

###### Theorem 2.3

Assume $z(t)=\sum_{k=1}^{N(t)}\xi_{k}$ is a compound poison process, that is,
N(t) is Poisson process with the instantaneous rate $\alpha$, and $\xi_{k}$
follows an exponential distribution with the rate $\theta$.
The stationary marginal distribution of $\sigma^{2}(t)$ is Gamma distribution
$\Gamma(\alpha,\theta)$

###### Proof.

$\displaystyle\sigma^{2}(t+u)=\int_{-\infty}^{t+u}e^{-\lambda(t+u-s)}dz(\lambda
s)=e^{-\lambda u}\sigma^{2}(t)+\int_{0}^{u}e^{-\lambda(u-s)}dz(\lambda s)\quad
u\geq 0$ (2.8)

The stationary solution $\sigma^{2}(t)$ of (2.3) can be written as in (2.8).
Because of the stationarity, we have

$\displaystyle\vartheta(\xi)$ $\displaystyle=\vartheta(\xi e^{-\lambda
u})\Phi(u,\xi)$ (2.9)

$\vartheta(\xi)$ is the characteristic function of the stationary distribution
of $\sigma^{2}(t)$ and $\Phi(u,\xi)$ is the characteristic function of
$\int_{0}^{u}e^{-\lambda(u-s)}dz(\lambda s)$. We have $0\leq e^{-\lambda
u}\leq 1$ for $u\geq 0$, and the relation (2.8) shows that $\sigma^{2}(t)$ is
self-decomposable.
$z(t)$ is a compound poison process with the function characteristic.

$\displaystyle g(\xi)=E(e^{i\xi z(1)})=exp\left\\{\int_{0}^{\infty}(e^{i\xi
x}-1)\alpha
f(x)dx\right\\}=exp({\rho(\xi)})\quad\quad\rho(\xi)=\frac{i\xi\alpha}{\theta-i\xi}$
(2.10)

It was shown in [6] that $\Phi(u,\xi)$ can be expressed as follows

$\displaystyle\Phi(u,\xi)=exp\left\\{\lambda\int_{0}^{u}\rho(\xi
e^{-\lambda(u-s)})ds\right\\}=exp\left\\{\int_{\xi e^{-\lambda
u}}^{\xi}\frac{\rho(w)}{w}dw\right\\}$ (2.11)

By replacing, $\frac{\rho(w)}{w}=\frac{i\alpha}{\theta-iw}$, we have

$\displaystyle\Phi(u,\xi)=\left(\frac{\theta-i\xi e^{-\lambda
u}}{\theta-i\xi}\right)^{\alpha}$ (2.12)

$\vartheta(\xi)$ is continuous at zero, and we have:

$\displaystyle\vartheta(\xi)=\lim_{u\to\infty}\vartheta(\xi e^{-\lambda
u})\Phi(u,\xi)=\left(\frac{1}{1-i\frac{1}{\theta}\xi}\right)^{\alpha}=\left(1-i{\theta}^{-1}\xi\right)^{-\alpha}$
(2.13)

From (2.13), $\vartheta(\xi)$ is the function characteristics of the gamma
distribution; and the stationary marginal distribution of $\sigma^{2}(t)$ is
the $\Gamma(\alpha,\theta)$ Gamma distribution.
Another method developed in [7, 9, 11, 8, 10] uses the relationship between
the $z(t)$ Lévy density $w(x)$ and the Lévy density $u(x)$ of $\sigma^{2}(t)$.

$\displaystyle u(x)=\int_{1}^{\infty}w(xr)dr $ (2.14)

From (2.10), we have the Lévy density $W(x)=\alpha f(x)=\alpha\theta
e^{-\theta x}$ and the Lévy density $u(x)$ of $\sigma^{2}(t)$ can be deduced
as follows.

$\displaystyle u(x)=\int_{1}^{\infty}\alpha\theta e^{-\theta
xr}dr=\frac{\alpha}{x}e^{-\theta x}\quad x>0$ (2.15)

u(x) is the Lévy density of Gamma distribution $\Gamma(\alpha,\theta)$. ∎

We can integrate the stationary non-negative process $\sigma^{2}(t)$.

$\displaystyle\sigma^{2*}(t)=\int_{0}^{t}\sigma^{2}(s)ds$ (2.16)

By integration by part method, (2.16) becomes

$\displaystyle\sigma^{2*}(t)$
$\displaystyle=\lambda^{-1}\sigma^{2}(0)(1-e^{-\lambda
t})+\lambda^{-1}\int_{0}^{t}\left(1-e^{-\lambda(t-s)}\right)dz(\lambda s)$
(2.17) $\displaystyle=\lambda^{-1}\left(-\sigma^{2}(t)+z(\lambda
t)+\sigma^{2}(0)\right)$ (2.18)

It results from (2.18) that the process $\sigma^{2*}(t)$ is continuous as
$z(\lambda t)$ and $\sigma^{2}(t)$ co-break [10, 9]. In addition, the shape of
$\sigma^{2*}(t)$ is determined by $z(\lambda t)$. In fact, $\sigma^{2*}(t)$
and $z(\lambda t)$ co-integrate. The co-integration can be shown by
transforming the equation (2.18) into (2.19). $\lambda\sigma^{2*}(t)-z(\lambda
t)$ is a stationary process such that.

$\displaystyle\lambda\sigma^{2*}(t)-z(\lambda t)=-\sigma^{2}(t)+\sigma^{2}(0)$
(2.19)

For $\lambda=1$ and $\sigma^{2}(0)=0$, the compound poison process ($z(t)$),
the $\Gamma(\alpha,\theta)$ Ornstein-Uhlenbeck process in ($\sigma^{2}(t)$),
and $\sigma^{2*}(t)$ in (2.20) were simulated and the results are in Fig 1(a),
Fig 1(b), and Fig 1(c) respectively.

$\displaystyle
z(t)=\sum_{k=1}^{N(t)}\xi_{k}\quad\quad\sigma^{2}(t)=\sigma^{2}(0)e^{\lambda
t}+\sum_{k=1}^{N(t)}exp(-\lambda(t-a_{k}))\xi_{k}\quad\quad\sigma^{2*}(t)=\int_{0}^{t}\sigma^{2}(s)ds$
(2.20)

(a) Compound Poison process: $\hat{z}(t)$

(b) Gamma process: $\hat{\sigma}^{2}(t)$

(c) Subordinator: $\hat{\sigma}^{2*}(t)$

Figure 1. Simulations: $\hat{\alpha}=0.8845$, $\hat{\theta}=0.9378$

The estimations of the gamma distribution parameter $\Gamma(\alpha,\theta)$
were performed by the FRFT Maximum likelihood on the daily SPY prices [31].

###  2.3. Variance - Gamma Process: Semi-Martingale

Let $Y^{*}=\\{Y_{t}^{*}\\}$, a stochastic process used to model the log of an
asset price.

$\displaystyle Y_{t}^{*}=A_{t}+M_{t}\quad$ $\displaystyle\quad A_{t}=\beta
t+\delta\sigma^{2*}(t)$ (2.21) $\displaystyle M_{t}$
$\displaystyle=\sigma\int_{0}^{t}\sigma(s)dW(s)$ (2.22)

where $\beta$ and $\delta$ are the drift parameters, $t$ represents the
continuous time clock, and $W(t)$ is the standard Brownian motion and
independent of $\sigma^{2}(t)$.

$\displaystyle\sigma(t)=\sqrt{\sigma^{2}(t)}\quad\quad\sigma^{2*}(t)$
$\displaystyle=\int_{0}^{t}\sigma^{2}(s)ds$ (2.23)

$\sigma(t)$ is the spot or instantaneous volatility, and $\sigma^{2*}(t)$ is
the chronometer or the integrated variance of the process. As shown in Fig
1(c), the Gamma process ($\sigma^{2*}(t)$) is a strictly increasing process of
the stationary process ($\sigma^{2}(t)$).
The mean process $A_{t}$ is a predictable process with locally bounded
variation. In fact, $A_{t}$ is continuous and differentiable because of
$\sigma^{2*}(t)$.
$M_{t}$ is a local martingale. The derivative of $M_{t}$ in (2.22) can be
written as a Stochastic Differential Equation (SDE) (2.24)

$\displaystyle dM_{t}$ $\displaystyle=\sigma\sigma(t)dW(t)$ (2.24)

$Y_{t}^{*}$ is a special semi-martingale [35, 10] and the decomposition
$Y_{t}^{*}=A_{t}+M_{t}$ is unique.

(a) Simulations: $Y^{*}=\\{Y_{t}^{*}\\}$

(b) Simulations versus SPY ETF data

Figure 2. VG Model: $\hat{\mu}=0.0848$, $\hat{\delta}=-0.0577$,
$\hat{\sigma}=1.0295$, $\hat{\alpha}=0.8845$, $\hat{\theta}=0.9378$

Fig 2(a) and Fig 2(b) display, in blue color, the simulation of the
logarithmic of the asset price ($Y^{*}$) in (2.21). The simulation is compared
with the daily SPY ETF historical return data from January 4, 2010, to
December 30, 2020. The daily SPY ETF data is displayed in red color in Fig
2(a).

###  2.4. Variance - Gamma Process: Parameter Estimations

The stochastic process in (2.21) is the solution of the following Stochastic
Differential Equation (SDE):

$\displaystyle dY_{t}^{*}=(\beta+\delta\sigma^{2}(t))dt+\sigma\sigma(t)dW(t)$
(2.25)

Given an interval of length $\Delta$, we define $\sigma^{2}_{n}$ and $Y_{n}$
over the interval [$(n-1)\Delta$; $n\Delta$].

$\displaystyle\sigma^{2}_{n}=\int_{(n-1)\Delta}^{n\Delta}d\sigma^{2*}(s)=\sigma_{n\Delta}^{2*}-\sigma_{(n-1)\Delta}^{2*}\quad
Y_{n}=\int_{(n-1)\Delta}^{n\Delta}dY_{s}^{*}=Y_{n\Delta}^{*}-Y_{(n-1)\Delta}^{*}$
(2.26)

The volatility component can be transformed into a normally distributed
variable $X(1)$ as follows:

$\displaystyle\int_{(n-1)\Delta}^{n\Delta}\sigma(t)dW(t)$
$\displaystyle\overset{\text{d}}{=}N\left(0,\int_{(n-1)\Delta}^{n\Delta}\sigma^{2}(s)ds\right)=N\left(0,\sigma_{n\Delta}^{2*}-\sigma_{(n-1)\Delta}^{2*}\right)=N\left(0,\sigma^{2}_{n}\right)$
(2.27)
$\displaystyle\overset{\text{d}}{=}{\sigma_{n}}N\left(0,1\right)\overset{\text{d}}{=}{\sigma_{n}}X(1).$

where $X(1)\overset{\text{d}}{=}N(0,1)$ and $N(0,1)$ denotes a standard normal
distribution.

By integrating the instantaneous return rate (2.25) per component, we have:

$\displaystyle\int_{(n-1)\Delta}^{n\Delta}dY_{s}^{*}=\beta\Delta+\delta\int_{(n-1)\Delta}^{n\Delta}d\sigma^{2*}(s)+\sigma\int_{(n-1)\Delta}^{n\Delta}\sigma(t)dW(t)$

Based on (2.26) and (2.27), we have the following equation over the interval
[$(n-1)\Delta$; $n\Delta$]

$\displaystyle
Y_{n}=\mu+\delta\sigma^{2}_{n}+\sigma{\sigma_{n}}X(1)\quad\quad\quad\mu=\beta\Delta\quad\quad\sigma^{2}_{n}\overset{\text{d}}{=}\Gamma(\alpha,\theta)$
(2.28)

In case $\Delta$ is a daily length, $Y_{n}$ becomes the daily return rate. The
equation (2.28) was analyzed in [31, 32] as a daily return rate, and the
parameters were estimated. The data came from the daily SPY ETF historical
data and the period spans from January 4, 2010, to December 30, 2020. See [31,
32, 34] for more details on the methodology and the results.

Table 1 presents the estimation results of the five parameters
$(\mu,\delta,\alpha,\theta,\sigma)$ of $Y_{n}$ in (2.28) along with some
statistical indicators.

Table 1. FRFT Maximum Likelihood VG Parameters Estimations Model | Parameters | Statistics
---|---|---
| $\hat{\mu}=0.0848$ | $\hat{E(Y)}=0.0369$
$\hat{\delta}=-0.0577$ | $\hat{Var(Y)}=0.8817$
VG | $\hat{\sigma}=1.0295$ | $\hat{Skew(Y)}=-0.173$
$\hat{\alpha}=0.8845$ | $\hat{Kurt(Y)}=6.412$
$\hat{\theta}=0.9378$
Source: Nzokem(2021) [31, 32]

As shown in Table 2, with initial parameter values ($\sigma=\alpha=\theta=1$,
$\delta=\mu=0$), the maximization procedure convergences after 21 iterations.
$\log(ML)$ is the function to maximize and $||\frac{d\log(ML)}{dV}||$ is the
norm of the partial derivative function ($\log(ML)$). During the maximization
process, both quantities converge respectively to $-3549.692$ and $0$; where
the parameter vector is stable. The location parameter $\mu$ is positive, the
symmetric parameter $\delta$ is negative, and other parameters have the
expected sign.

Table 2. Results of VG Model Parameter Estimations

Iterations | $\mu$ | $\delta$ | $\sigma$ | $\alpha$ | $\theta$ | $Log(ML)$ | $||\frac{dLog(ML)}{dV}||$
---|---|---|---|---|---|---|---
1 | 0 | 0 | 1 | 1 | 1 | -3582.8388 | 598.743231
2 | 0.05905599 | -0.0009445 | 1.03195903 | 0.9130208 | 1.03208412 | -3561.5099 | 833.530396
3 | 0.06949925 | 0.00400035 | 1.04101444 | 0.88478895 | 1.05131996 | -3559.5656 | 447.807305
4 | 0.07514039 | 0.00055771 | 1.17577397 | 0.67326429 | 1.17778666 | -3569.6221 | 211.365781
5 | 0.08928373 | -0.0263716 | 1.03756321 | 0.83842661 | 0.94304967 | -3554.4434 | 498.289445
6 | 0.08676498 | -0.0521887 | 1.03337015 | 0.85591875 | 0.95066351 | -3550.6419 | 204.467192
7 | 0.086995 | -0.0608517 | 1.02788937 | 0.87382621 | 0.95054954 | -3549.8465 | 66.8039738
8 | 0.08542912 | -0.058547 | 1.02705241 | 0.88258411 | 0.94321299 | -3549.7023 | 15.3209117
9 | 0.08478622 | -0.0576654 | 1.02995166 | 0.88447791 | 0.93670036 | -3549.6921 | 1.14764198
10 | 0.08477798 | -0.0577736 | 1.02922308 | 0.88449072 | 0.93831041 | -3549.692 | 0.17287708
11 | 0.08476475 | -0.0577271 | 1.02960343 | 0.88450434 | 0.93755549 | -3549.692 | 0.07850459
12 | 0.08477094 | -0.0577488 | 1.02942608 | 0.8844984 | 0.93790784 | -3549.692 | 0.03723941
13 | 0.08476804 | -0.0577386 | 1.02950937 | 0.88450117 | 0.93774266 | -3549.692 | 0.01732146
14 | 0.0847694 | -0.0577434 | 1.02947043 | 0.88449987 | 0.93781995 | -3549.692 | 0.00813465
15 | 0.08476876 | -0.0577411 | 1.02948868 | 0.88450048 | 0.93778375 | -3549.692 | 0.00380345
16 | 0.08476906 | -0.0577422 | 1.02948014 | 0.88450019 | 0.9378007 | -3549.692 | 0.00178206
17 | 0.08476892 | -0.0577417 | 1.02948414 | 0.88450033 | 0.93779276 | -3549.692 | 0.00083415
18 | 0.08476898 | -0.0577419 | 1.02948226 | 0.88450026 | 0.93779648 | -3549.692 | 0.00039063
19 | 0.08476895 | -0.0577418 | 1.02948314 | 0.88450029 | 0.93779474 | -3549.692 | 0.00018289
20 | 0.08476897 | -0.0577419 | 1.02948273 | 0.88450028 | 0.93779555 | -3549.692 | 8.56E-05
21 | 0.08476896 | -0.0577418 | 1.02948292 | 0.88450029 | 0.93779517 | -3549.692 | 4.01E-05

##  3\. Variance - Gamma Process: Probability versus Lévy Density

Based on (2.21) and (2.22), the VG Process $Y=\\{Y_{t}\\}_{t\geq 0}$ with five
parameters $(\mu,\delta,\sigma,\alpha,\theta)$ can be written as follows:

$\displaystyle Y_{t}=\mu
t+\delta\sigma^{2*}(t)+\sigma\int_{0}^{t}\sigma(s)dW(s).$ (3.1)

where $\mu,\delta\in R$ , $\sigma>0$, $\alpha>0$, $\theta>0$, $t$ represents
the continuous time clock, and $W(t)$ is the standard Brownian motion and
independent of $\sigma^{2}(t)$.

$\displaystyle\sigma(t)=\sqrt{\sigma^{2}(t)}\quad\quad\sigma^{2*}(t)$
$\displaystyle=\int_{0}^{t}\sigma^{2}(s)ds.$ (3.2)

where $\sigma(t)$ is the spot or instantaneous volatility, $\sigma^{2}(t)$ is
the spot or instantaneous variance, and $\sigma^{2*}(t)$ is the chronometer or
the integrated variance of the process.

We consider the characteristic function of the VG process $Y=\\{Y_{t}\\}$

$\displaystyle E\left[e^{i\xi Y_{t}}\right]=E\left[e^{i\xi(\mu
t+\delta\sigma^{2*}(t)+\sigma\int_{0}^{t}\sigma(s)dW(s))}\right]=e^{it\mu\xi}E\left[e^{i\xi(\delta\sigma^{2*}(t)+\sigma\int_{0}^{t}\sigma(s)dW(s))}\right].
$ (3.3)

$\int_{0}^{t}\sigma(s)dW(s)$ is the It$\hat{o}$ integral with respect to the
Brownian motion, and we have :

$\displaystyle\int_{0}^{t}\sigma(s)dW(s)$
$\displaystyle\overset{\text{d}}{=}N\left(0,\int_{0}^{t}\sigma^{2}(s)ds\right)=N\left(0,\sigma^{2*}(t)\right).
$ (3.4)

where $N(0,1)$ is a standard normal distribution.

From expressions (3.3) and (3.4), we have

$\displaystyle
E\left[e^{i\xi(\delta\sigma^{2*}(t)+\sigma\int_{0}^{t}\sigma(s)dW(s))}\right]$
$\displaystyle=E\left[e^{i\xi
N\left(\delta\sigma^{2*}(t),\sigma^{2}\sigma^{2*}(t)\right)}\right]=E\left[E\left[e^{i\xi
N\left(\delta\sigma^{2*}(t),\sigma^{2}\sigma^{2*}(t)\right)}|\sigma^{2*}(t)\right]\right]
$ (3.5)
$\displaystyle=E\left[e^{(i\delta\xi-\frac{1}{2}\sigma^{2}\xi^{2})\sigma^{2*}(t)}\right]
$

$\sigma^{2*}(t)$ is a Lévy process generated by the Gamma distribution
$\Gamma(\alpha,\theta)$ and we have

$\displaystyle
E\left[e^{(i\delta\xi-\frac{1}{2}\sigma^{2}\xi^{2})\sigma^{2*}(t)}\right]$
$\displaystyle=\frac{1}{(1+\frac{1}{2}\theta\sigma^{2}\xi^{2})^{t\alpha}}E\left[e^{i\delta\xi
W}\right] \quad
\hbox{$\sigma^{2*}(t)\overset{\text{d}}{=}\Gamma(t\alpha,\theta)$} $ (3.6)
$\displaystyle=\frac{1}{\left(1-i\delta\theta\xi+\frac{1}{2}\sigma^{2}\theta\xi^{2}\right)^{t\alpha}}\quad
\hbox{$W\overset{\text{d}}{=}\Gamma(t\alpha,\frac{1}{2}\sigma^{2}\xi^{2}+\frac{1}{\theta})$}
$

From expressions (3.3), (3.5) and (3.6), we have

$\displaystyle E\left[e^{iY_{t}\xi}\right]$
$\displaystyle=\frac{e^{it\mu\xi}}{\left(1-i\delta\theta\xi+\frac{1}{2}\sigma^{2}\theta\xi^{2}\right)^{t\alpha}}$
(3.7)

We define two related functions $\phi(\xi)$ and $\varphi(\xi,t)$ such that

$\displaystyle\phi(\xi)=\frac{e^{i\mu\xi}}{\left(1-i\delta\theta\xi+\frac{1}{2}\sigma^{2}\theta\xi^{2}\right)^{\alpha}}\quad$
$\displaystyle\quad\quad\varphi(\xi,t)=-Log(E\left[e^{iY_{t}\xi}\right])=-tLog(\phi(\xi))$
(3.8)

The characteristic function can be written as follows.

$\displaystyle
E\left[e^{iY_{t}\xi}\right]=\left(\phi(\xi)\right)^{t}=E[e^{-tLog(\phi(\xi))}]$
(3.9)

####  3.0.1. Lévy measure and the structure of the jumps

###### Lemma 3.1

(Frullani integral) $\forall\alpha,\beta>0$ and $\forall z\in\mathscr{C}$ with
$\Re(z)\leq 0$.
we have

$\displaystyle
\frac{1}{\left(1-\frac{z}{\alpha}\right)^{\beta}}=e^{-\int_{0}^{\infty}(1-e^{zx})\beta
x^{-1}e^{-\alpha x}dx} $

For lemma proof, see [5]

###### Theorem 3.2

(Variance-Gamma Model representation)
Let $Y=\\{Y_{t}\\}_{t\geq 0}$ be a Lévy process on $\mathbb{R}$ generated by
the VG model with parameter $(\mu,\delta,\sigma,\alpha,\theta)$. The
characteristic exponent of the Lévy process has the following representation.

$\displaystyle\varphi(\xi,1)$ $\displaystyle=-
Log\left(Ee^{iY_{1}\xi}\right)=i\mu\xi+\int_{-\infty}^{+\infty}\left(1-e^{-i\xi
u}\right)\Pi(u)du $ (3.10)

$\Pi(u)$ is the Lévy density of $Y$ :

$\displaystyle \Pi(u)=\alpha\left(\frac{1_{\\{u>0\\}}}{u}e^{-x_{1}u}
+\frac{1_{\\{u<0\\}}}{\lvert u\rvert}e^{{-x_{2}}u}\right)$ (3.11)

with

$\displaystyle
x_{1}=\frac{\delta}{\sigma^{2}}+\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}
\quad\quad
x_{2}=\frac{\delta}{\sigma^{2}}-\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}$
(3.12)

and $\Pi(u)$ satisfies the properties

$\displaystyle\int_{-\infty}^{+\infty}\Pi(u)du=+\infty\quad\hbox{and}\quad\int_{-\infty}^{+\infty}Min(1,|u|)\Pi(u)du<+\infty$
(3.13)

###### Proof.

We consider the characteristic function $\phi(\xi)$ in (3.8) of the VG model
with parameter $(\mu,\delta,\sigma,\alpha,\theta)$ developed previously

$\displaystyle\phi(\xi)=\frac{e^{i\mu\xi}}{\left(1-i\delta\theta\xi+\frac{1}{2}\sigma^{2}\theta\xi^{2}\right)^{\alpha}}$

We factor the quadratic function in the denominator of $\phi(\xi)$.

$\displaystyle\left(1+\frac{1}{2}\theta\sigma^{2}x^{2}-i\delta\theta
x\right)^{\alpha}=\left(\frac{1}{2}\theta\sigma^{2}\right)^{\alpha}\left(x-ix_{1}\right)^{\alpha}\left(x-ix_{2}\right)^{\alpha}.$
(3.14)

with

$\displaystyle
x_{1}=\frac{\delta}{\sigma^{2}}+\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}\quad\quad\quad
x_{2}=\frac{\delta}{\sigma^{2}}-\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}$

We apply the lemma 3.1 on each factor of the quadratic function (3.14).

$\displaystyle\left(\frac{1}{2}\theta\sigma^{2}\right)^{\alpha}{\left(x-ix_{1}\right)^{\alpha}}{\left(x-ix_{2}\right)^{\alpha}}$
$\displaystyle={\left(1+\frac{ix}{x_{1}}\right)^{\alpha}}{\left(1+\frac{ix}{x_{2}}\right)^{\alpha}}$
$\displaystyle=\left(e^{\int_{0}^{\infty}{(1-e^{-ixu})\frac{\alpha}{u}e^{-x_{1}u}du}}\right)\left(e^{\int_{0}^{\infty}{(1-e^{ixu})\frac{\alpha}{u}e^{x_{2}u}du}}\right)$
$\displaystyle=e^{\int_{0}^{\infty}{(1-e^{-ixu})\frac{\alpha}{u}e^{-x_{1}u}du}+\int_{-\infty}^{0}{(1-e^{-ixv})\frac{\alpha}{|v|}e^{-x_{2}v}dv}}$
$\displaystyle=e^{\int_{-\infty}^{+\infty}{(1-e^{-ixu})\Pi(u)du}}$

we take into account the expression (3.14) and have

$\displaystyle\left(1+\frac{1}{2}\theta\sigma^{2}x^{2}-i\delta\theta
x\right)^{\alpha}=e^{\int_{-\infty}^{+\infty}{(1-e^{-ixu})\Pi(u)du}}.$ (3.15)

where $\Pi(u)=\alpha\left(\frac{1_{\\{u>0\\}}}{u}e^{-x_{1}u}
+\frac{1_{\\{u<0\\}}}{\lvert u\rvert}e^{{-x_{2}}u}\right)$
From expression (3.8), we have:

$\displaystyle\varphi(\xi,t)=-tLog(\phi(\xi))$
$\displaystyle=-it\mu\xi+tLog\left(1+\frac{1}{2}\theta\sigma^{2}x^{2}-i\delta\theta
x\right)^{\alpha}$
$\displaystyle=-it\mu\xi+tLog\left(1+\frac{1}{2}\theta\sigma^{2}x^{2}-i\delta\theta
x\right)^{\alpha}$
$\displaystyle=-it\mu\xi+\int_{-\infty}^{+\infty}{(1-e^{-ixu})t\Pi(u)du}$

We have

$\displaystyle\varphi(\xi,t)=-tLog(\phi(\xi))=-it\mu\xi+\int_{-\infty}^{+\infty}{(1-e^{-ixu})t\Pi(u)du}$
(3.16)

For $t=1$, we have the expression (3.10)
We can check some properties of $\Pi(u)$

$\displaystyle\int_{-\infty}^{+\infty}\Pi(u)du$
$\displaystyle=+\infty\quad\quad \hbox{ in fact \quad$\lim_{|u|\to
0}\Pi(u)=+\infty$ }$ (3.17)
$\displaystyle\int_{-\infty}^{+\infty}Min(1,|u|)\Pi(du)$
$\displaystyle=\int_{-1}^{1}Min(1,|u|)\Pi(du)+\int_{1}^{+\infty}Min(1,|u|)\Pi(du)$
$\displaystyle+\int_{-\infty}^{-1}Min(1,|u|)\Pi(du)$
$\displaystyle=\alpha\left(\frac{1-e^{-x_{1}}}{x_{1}}+\frac{1-e^{x_{2}}}{-x_{2}}+\Gamma(0,x_{1})+\Gamma(0,-x_{2})\right).$

with $\Gamma(s,u)=\int_{u}^{+\infty}{y^{s-1}}e^{-y}dy$
And we have:

$\displaystyle\int_{-\infty}^{+\infty}Min(1,|u|)\Pi(du)<+\infty$ (3.18)

∎

The results in (3.17) show that the VG process is not a finite activities
process and can not be written as a Compound Poisson process [10]. The VG
process is an infinite activity process with an infinite number of jumps in
any given time interval. The arrival rate of jumps of all sizes in the VG
process is defined by the Lévy density (3.19)

$\displaystyle\Pi(u)=\begin{cases}\frac{\alpha}{\lvert
u\rvert}{e^{-x_{2}u}}&\text{if }u<0\\\ \frac{\alpha}{u}{e^{-x_{1}u}}&\text{if
}u>0.\end{cases}$ (3.19)

As shown in Fig 3(a), the high arrival rates of jumps are concentrated around
the origin $0$. The smaller the jump size, the higher the arrival rate for the
VG model. The steepness parameters [13],$-x_{2}$ and $x_{1}$, defined the rate
of exponential decay of the tails on each side. As shown in Fig 3(a) and
(3.19), the Lévy density is asymmetric, and the left tail is heavier as
$-x_{2}<x_{1}$. On the other hand, the result in (3.18) proves that the VG
process is a finite variation process, which is contrary to the Brownian
motion process. The Gamma distribution Parameter ($\alpha$), called the
process intensity [13], plays an important role in the Lévy density. The
intensity of the process $(\alpha)$ has a similar role as the variance
parameter in the Brownian motion process. The Lévy density function (3.19) is
different for negative and positive jump size. The difference has led [25] to
see the VG process as the difference between two increasing processes, with
one process providing the upward movement and another the downward movement in
the market.

(a) Lévy density of the VG model

(b) Probability density of the VG model

Figure 3. VG Model: $\hat{\mu}=0.0848$, $\hat{\delta}=-0.0577$,
$\hat{\sigma}=1.0295$, $\hat{\alpha}=0.8845$, $\hat{\theta}=0.9378$

Using the VG parameter estimations in table 1, we have $x_{1}=1.4775$ and
$x_{2}=-1.3640$. Fig 3(a) and Fig 3(b) display the Lévy and the probability
densities. As shown in Fig 3, the shape of the density functions are
different; even-though, both densities are linked by the same characteristic
function.
Variance-Gamma (VG) Process can be described as a subfamily of the KoBoL
family, which is the extension of Koponen’s family by Boyarchenko and
Levendorskii [13]. The KoBoL family is also called CGMY- model (named after
Carr, German, Madan, and Yor) [14]. Under the KoBoL family, the Lévy density
has the following general form. See [13] for more details

$\displaystyle\Pi(u)=\begin{cases}C_{-}{|u|^{-\nu-1}}{e^{\lambda_{-}u}}&\text{if
}u<0\\\ C_{+}{u^{-\nu-1}}{e^{-\lambda_{+}u}}&\text{if }u>0\end{cases}$ (3.20)

where $C_{+}>0$, $C_{-}>0$, $\nu>0$ and $\lambda_{-}<0<\lambda_{+}$
As a subfamily of the KoBoL family, the VG process belongs to the process
class of order $\nu=0$, intensity $C_{+}=C_{-}=\alpha$ and steepness
parameters
$\lambda_{-}=-x_{2}=-\frac{\delta}{\sigma^{2}}+\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}$
and
$\lambda_{+}=x_{1}=\frac{\delta}{\sigma^{2}}+\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}$.
For $0<\nu<1$, see [30] for a general case of tempered stable distribution.

###  3.1. Variance - Gamma Process: Asymptotic distribution

###### Theorem 3.3

(Variance-Gamma process probability density)
Let $Y=\\{Y_{t}\\}_{t\geq 0}$ be a Lévy process on $\mathbb{R}$ generated by
the VG model with parameter $(\mu,\delta,\sigma,\alpha,\theta)$. The
probability density function can be written as follows

$\displaystyle f(y,t)$
$\displaystyle=\frac{1}{\sigma\Gamma(t\alpha)\theta^{t\alpha}}\int_{0}^{+\infty}\frac{1}{\sqrt{2\pi
v}}e^{-\frac{(y-t\mu-\delta
v)^{2}}{2v\sigma^{2}}}v^{t\alpha-1}e^{-\frac{v}{\theta}}\,dv\quad\quad t\geq
0\quad y\in\mathbb{R} $ (3.21)

###### Proof:.

$\varphi(\xi,t)$ in (3.16) provides the relation between the characteristic
exponent and the Lévy density. the expression is used as follows:

$\displaystyle\varphi(\xi,t)= -Log\left(Ee^{iY_{t}\xi}\right)$
$\displaystyle=-it\mu\xi+\int_{-\infty}^{+\infty}{(1-e^{-ixu})t\Pi(u)du}$
$\displaystyle t \Pi(u)$
$\displaystyle=t\alpha\left(\frac{1_{\\{u>0\\}}}{u}e^{-x_{1}u}
+\frac{1_{\\{u<0\\}}}{\lvert u\rvert}e^{{\lvert x_{2}\rvert}u}\right) $
$\displaystyle\mu_{t}=t\mu \quad\quad\alpha_{t}=t\alpha $ (3.22)

It was shown in [31] that the probability density of a VG model with parameter
$(\mu,\delta,\sigma,\alpha,\theta)$ can be written

$\displaystyle f(y)$
$\displaystyle=\frac{1}{\sigma\Gamma(\alpha)\theta^{\alpha}}\int_{0}^{+\infty}\frac{1}{\sqrt{2\pi
v}}e^{-\frac{(y-\mu-\delta
v)^{2}}{2v\sigma^{2}}}v^{\alpha-1}e^{-\frac{v}{\theta}}\,dv $

By replacing the parameters in (3.22), we have the result in Theorem 3.3. ∎

###### Theorem 3.4

(Asymptotic distribution of Variance-Gamma process)
Let $Y=\\{Y_{t}\\}_{t\geq 0}$ be a Lévy process on $\mathbb{R}$ generated by
the VG model with parameter $(\mu,\delta,\sigma,\alpha,\theta)$.
Then $Y_{t}$ converges in distribution to a Lévy process driving by a Normal
distribution with mean $a=\mu+\alpha\theta\delta$ and variance
$\sigma^{2}=\alpha(\theta^{2}\delta^{2}+\sigma^{2}\theta)$.

$\displaystyle Y_{t}\overset{\text{d}}{\sim}N(ta,t\sigma^{2})\quad$
$\displaystyle\text{as}\quad t\to+\infty $ (3.23)

###### Proof:.

Let us have

$\displaystyle b_{t}$ $\displaystyle=\sqrt{t}b\ $ $\displaystyle\ a_{t}$
$\displaystyle=ta$ $\displaystyle b$
$\displaystyle=\sqrt{\alpha(\theta^{2}\delta^{2}+\sigma^{2}\theta)}\ $
$\displaystyle\ a$ $\displaystyle=\mu+\alpha\theta\delta $

$\phi(\xi,t)$ is the characteristic function of the process
$Y=\\{Y_{t}\\}_{t\geq 0}$ , we use the expression (3.7).

$\displaystyle\phi(\xi,t)=E\left[e^{iY_{t}\xi}\right]$
$\displaystyle=\frac{e^{it\mu\xi}}{\left(1-i\delta\theta\xi+\frac{1}{2}\sigma^{2}\theta\xi^{2}\right)^{t\alpha}}$

$\phi^{T}(\xi,t)$ is the characteristic function of the stochastic process
$\\{\frac{{Y_{t}}-a_{t}}{b_{t}}\\}_{t\geq 0}$ and we have

$\displaystyle\phi^{T}(\xi,t)=
E\\{e^{i\frac{{Y_{t}}-a_{t}}{b_{t}}\xi}\\}=e^{-i\frac{a_{t}}{b_{t}}\xi}E\\{e^{i\frac{\xi}{b_{t}}Y_{t}}\\}$
$\displaystyle=e^{-i\frac{a_{t}}{b_{t}}\xi}\phi(\frac{\xi}{b_{t}},t)=\frac{e^{it\alpha\theta\delta\frac{\xi}{b_{t}}}}{\left(1+\frac{1}{2}\theta\sigma^{2}\frac{\xi^{2}}{b_{t}^{2}}-i\delta\theta\frac{\xi}{b_{t}}\right)^{t\alpha}}$
$\displaystyle=e^{it\alpha\theta\delta\frac{\xi}{b_{t}}}{\left(1+\frac{1}{2}\theta\sigma^{2}\frac{\xi^{2}}{tb^{2}}-i\delta\theta\frac{\xi}{\sqrt{t}b}\right)^{-t\alpha}}
$

Let us have

$\displaystyle
u(t)=\frac{1}{2}\theta\sigma^{2}\frac{\xi^{2}}{tb^{2}}-i\delta\theta\frac{\xi}{\sqrt{t}b}\quad\quad\quad\lim_{t\to+\infty}u(t)=0
$

We can use the Taylor expansions of $\ln(1+u)$

$\displaystyle\ln(1+\frac{1}{2}\theta\sigma^{2}\frac{\xi^{2}}{tb^{2}}-i\delta\theta\frac{\xi}{\sqrt{t}b})$
$\displaystyle=\frac{1}{2}(\theta\sigma^{2}+\delta^{2}\theta^{2})\frac{\xi^{2}}{tb^{2}}-i\delta\theta\frac{\xi}{\sqrt{t}b}+o\left(\frac{1}{t\sqrt{t}}\right)$
$\displaystyle\lim_{t\to+\infty}o\left(\frac{1}{t\sqrt{t}}\right)$
$\displaystyle=0 $

The characteristic function, $\phi^{T}(\xi,t)$, developed previously becomes:

$\displaystyle\phi^{T}(\xi,t)=e^{it\alpha\theta\delta\frac{\xi}{b_{t}}}{\left(1+\frac{1}{2}\theta\sigma^{2}\frac{\xi^{2}}{tb^{2}}-i\delta\theta\frac{\xi}{\sqrt{t}b}\right)^{-t\alpha}}$
$\displaystyle=e^{it\alpha\theta\delta\frac{\xi}{b_{t}}}e^{-t\alpha\ln(1+\frac{1}{2}\theta\sigma^{2}\frac{\xi^{2}}{tb^{2}}-i\delta\theta\frac{\xi}{\sqrt{t}b})}$
$\displaystyle=e^{-\frac{1}{2}\alpha(\theta\sigma^{2}+\delta^{2}\theta^{2})\frac{\xi^{2}}{b^{2}}+o\left(\frac{1}{\sqrt{t}}\right)}$
$\displaystyle=e^{-\frac{1}{2}\xi^{2}+o\left(\frac{1}{\sqrt{t}}\right)} $

We have

$\displaystyle\lim_{t\to+\infty}\phi^{T}(\xi,t)=\lim_{t\to+\infty}E\\{e^{i\frac{{Y_{t}}-a_{t}}{b_{t}}\xi}\\}=e^{-\frac{1}{2}\xi^{2}}$
(3.24)

By applying the limit in (3.24), we produce the cumulant-generating function
[21] of the Normal distribution. We have the following convergence in
distribution

$\displaystyle\frac{{Y_{t}}-a_{t}}{b_{t}}\overset{\text{d}}{\sim}N(0,1)\quad$
$\displaystyle\text{as}\quad t\to+\infty $

∎

As shown in (3.22), the dynamic of the probability density $f(y,t)$ is carried
by two parameters: $t\mu$ and $t\alpha$. $f(y,t)$ can be compared to the
histogram of the daily SPY ETF return data, as shown in Fig 4(a). Fig 4(b)
shows the shape of the probability densities (3.21) adjusted at different
timeframes: Quarterly ($\tau=0.25$), Semi-Annual ($\tau=0.5$), Third-Quarterly
($\tau=0.75$), and Annual ($\tau=1$).

(a) Daily SPY ETF Return in (%)

(b) $f(y,\tau)$: $\tau$ in years

Figure 4. $f(y,t)$: $\hat{\mu}=0.0848$, $\hat{\delta}=-0.0577$,
$\hat{\sigma}=1.0295$, $\hat{\alpha}=0.8845$, $\hat{\theta}=0.9378$

The discrepancy between the shape of the probability densities (3.21) can be
explained by the Asymptotic distribution in Theorem 3.4. When the time frame
becomes large, The VG probability density generated by the Lévy Process
changes its nature and becomes a Normal distribution Process. Empirically, the
convergence is illustrated in Fig 4(b).

##  4\. Variance - Gamma Process: Pricing European Options

###  4.1. Variance - Gamma Process: Risk Neutral Esscher Measure

The method of Esscher transforms introduced by [16] as an efficient technique
for pricing derivative securities if a Lévy process models the logarithms of
the underlying asset prices. An Esscher transform of a stock-price process
provides an equivalent martingale measure; under such measure, the price of
any derivative security is calculated as the expectation of the discounted
payoffs. In some cases, the Esscher transform of a distribution [16] remains
in the family of the original distributions. In particular, Gamma,
Exponential, Normal, Inverse Gaussian, Negative Binomial, Geometric, Poisson,
and Compound Poisson distributions are examples of such conservative
distributions. The existence of the equivalent Esscher transform measure is
not always guaranteed, and the issue of the unicity of the equivalent
martingale measure remains recurrent when pricing an option with a Lévy
process.
From the characteristic function $\phi(\xi)$ in (3.8), We have the Moment
generating function of the VG model.

$\displaystyle M(h,t)$ $\displaystyle= \phi(-ih)^{t}=\frac{e^{t\mu
h}}{\left(1-\frac{1}{2}\theta\sigma^{2}h^{2}-\delta\theta
h\right)^{t\alpha}}={M(h,1)}^{t}\hskip 14.22636pt \hbox{with $h_{1}<h<h_{2}$
}$ (4.1) $\displaystyle M(h,1)$ $\displaystyle=\frac{e^{\mu
h}}{\left(1-\frac{1}{2}\theta\sigma^{2}h^{2}-\delta\theta
h\right)^{\alpha}}\hskip 14.22636pt \hbox{with $h_{1}<h<h_{2}$ }$
$\displaystyle h_{1}$
$\displaystyle=-\frac{\delta}{\sigma^{2}}-\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}
\quad
h_{2}=-\frac{\delta}{\sigma^{2}}+\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}$

Under the Esscher transform with parameter h, the probability density of
$Y={Y_{t}}$ becomes:

$\displaystyle\hat{f}(x,t,h)= \frac{e^{hx}{f(x,t)}}{M(h,t)}\quad\quad
\hbox{with $h_{1}<h<h_{2}$ }$ (4.2)

The Moment generating function of the Esscher transform VG model with
$h_{1}<h<h_{2}-z$

$\displaystyle M(z,t,h)=
E^{h}\left[e^{zX_{t}}\right]=\int_{0}^{+\infty}e^{zx}\hat{f}(x,t,h)dx$
$\displaystyle= \int_{0}^{+\infty}\frac{e^{(h+z)x}{f(x,t)}}{M(x,t)}dx$ (4.3)
$\displaystyle=\frac{M(h+z,t)}{M(h,t)}\quad \hbox{with $h_{1}<h<h_{2}-z$ }$
$\displaystyle=\left(\frac{M(h+z,1)}{M(h,1)}\right)^{t}=M(z,1,h)^{t}$

with

$\displaystyle M(z,1,h)=\frac{M(h+z,1)}{M(h,1)}=e^{\mu
z}{\left(M_{**}(z,1,h)\right)}^{\alpha}$ (4.4) $\displaystyle
M_{**}(z,1,h)=\frac{1-\frac{1}{2}\theta\sigma^{2}h^{2}-\delta\theta
h}{1-\frac{1}{2}\theta\sigma^{2}(h+z)^{2}-\delta\theta(h+z)}$

$\hat{f}(x,t,h)$ is the modified probability density of $f(x,t)$ defined in
(3.21). The function $exp(x)$ is a strictly increasing function, and the
probability measure generated by $\hat{f}(x,t,h)$ is equivalent to the
original probability measure generated by $f(x,t)$. Both probability measures
have the same null sets [16](sets of probability measure zero).
Given the process $\left\\{e^{-r\tau}S(\tau)\right\\}_{\tau\geq 0}$ with $r$
the constant risk-free rate of interest. We look into the conditions to have
$h=h^{*}$ such that

$\displaystyle E^{h^{*}}\left[e^{-r\tau}S(\tau)\right]=S(0)$ (4.5)

We have $S(\tau)=S(0)e^{Y_{\tau}}$, with $Y_{\tau}$ is the Variance - Gamma
process. The equation (4.5) becomes

$\displaystyle
e^{r\tau}=E^{h^{*}}\left[e^{Y_{\tau}}\right]=M(1,1,h^{*})^{\tau}=\left(\frac{M(h^{*}+1,1)}{M(h^{*},1)}\right)^{\tau}\quad
\hbox{with $h_{1}<h^{*}<h_{2}-1$ }$ (4.6)

The first condition is that

$\displaystyle
h_{2}-h_{1}>1\quad\quad\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}>\frac{1}{4}$

The equation (4.6) is equivalent to (4.7).

$\displaystyle e^{\frac{r-\mu}{\alpha}}$
$\displaystyle=M_{**}(1,1,h^{*})=\frac{M(h^{*}+1,1)}{M(h^{*},1)}=\frac{1-\frac{1}{2}\theta\sigma^{2}h^{*2}-\delta\theta
h^{*}}{1-\frac{1}{2}\theta\sigma^{2}(h^{*}+1)^{2}-\delta\theta(h^{*}+1)}$
(4.7)

We consider the function $g(h)$ defined as follows

$g(h)=\frac{1-\frac{1}{2}\theta\sigma^{2}h^{2}-\delta\theta
h}{1-\frac{1}{2}\theta\sigma^{2}(h+1)^{2}-\delta\theta(h+1)}$
$\displaystyle\frac{dg}{dh}(h)=\frac{\frac{1}{2}\theta^{2}\sigma^{4}h^{2}+\delta\theta^{2}(\frac{1}{2}\sigma^{2}+\delta)h+\delta\theta^{2}(\frac{1}{2}\sigma^{2}+\delta)+\theta\sigma^{2}}{(1-\frac{1}{2}\theta\sigma^{2}(h+1)^{2}-\delta\theta(h+1))^{2}}$

we have

$\displaystyle\frac{dg}{dh}(h)>0\quad\hbox{$h_{1}<h<h_{2}-1$}\quad\quad\quad\lim_{h\to
h_{1}}g(h)$
$\displaystyle=0\quad\quad\quad\lim_{h\to{h_{2}-1}^{-}}g(h)=+\infty$ (4.8)

(4.8) shows the existence and the unicity of $h^{*}$ in $[h_{1},h_{2}-1[$ such
that

$\displaystyle e^{\frac{r-\mu}{\alpha}}=g(h^{*})$

For VG Model in Table 1 [31, 32], the existence and unicity of $h^{*}$ can be
studied empirically as shown in Fig 5. Over the interval $[h_{1};h_{2}-1]$,
g(h) is an increasing function, as shown in Fig 5(a). Fig 5(b) provides the
solution $h^{*}$ of equation (4.7) for free interest rate less than 10%. The
solution $h^{*}$ increases with the free interest rate $r$.

(a) $g(h)$

(b) $e^{\frac{r-\mu}{\alpha}}=g(h)$

Figure 5. VG Model: $\hat{\mu}=0.0848$, $\hat{\delta}=-0.0577$,
$\hat{\sigma}=1.0295$, $\hat{\alpha}=0.8845$, $\hat{\theta}=0.9378$,
$h_{1}=-13.6511$ and $h_{2}=14.7399$

From the Esscher transform, we have the Equivalent Martingale Measure (EMM)
$\mathbf{Q}$, which can be written as the Radon-Nikodym derivative:

$\displaystyle\frac{d\mathbf{Q}}{d\mathbf{P}}=\frac{e^{h^{*}Y_{\tau}}}{M(h^{*},\tau)}=e^{h^{*}Y_{\tau}-log(M(h^{*},\tau))}$
(4.9)

and $E^{\mathbf{Q}}$ for expectation with respect to $\mathbf{Q}$

$\displaystyle E^{\mathbf{Q}}\left[e^{-r\tau}S(\tau)\right]$
$\displaystyle=E^{\mathbf{P}}\left[e^{-r\tau}S(\tau)\frac{d\mathbf{Q}}{d\mathbf{P}}\right]=S(0)E^{P}\left[e^{(1+h^{*})Y_{\tau}-log(M(h^{*},\tau))-r\tau}\right]$
$\displaystyle=S(0)E^{\mathbf{P}}\left[e^{(1+h^{*})Y_{\tau}}\right]e^{-log(M(h^{*},\tau))-r\tau}$
$\displaystyle=S(0)\frac{M(1+h^{*},\tau)}{M(h^{*},\tau)}e^{-r\tau}\quad\quad(\ref{eq:l34})\quad
e^{r\tau}=\left(\frac{M(h^{*}+1,\tau)}{M(h^{*},\tau)}\right)$
$\displaystyle=S(0)$

we have

$\displaystyle E^{\mathbf{Q}}\left[e^{-r\tau}S(\tau)\right]=S(0)$ (4.10)

###### Theorem 4.1

(Variance - Gamma Esscher transform distribution)
Esscher transform of Variance - Gamma process $Y=\\{Y_{t}\\}_{t\geq 0}$ with
parameter $(t\mu,\delta,\sigma,t\alpha,\theta)$ is also a Variance - Gamma
process with parameter $(t\mu,\tilde{\delta},\sigma,t\alpha,\tilde{\theta})$

$\displaystyle\tilde{\delta}=\delta+h\sigma^{2}\quad\quad\tilde{\theta}=\frac{\theta}{1-\frac{1}{2}\theta\sigma^{2}{h}^{2}-\delta\theta{h}}$
(4.11)

###### Proof:.

From (4.4), we have

$\displaystyle
M_{**}(z,1,h)=\frac{1-\frac{1}{2}\theta\sigma^{2}h^{2}-\delta\theta
h}{1-\frac{1}{2}\theta\sigma^{2}(h+z)^{2}-\delta\theta(h+z)}$ (4.12)

We can divide the denominator by the numerator of the function
$M_{**}(z,1,{h})$ in (4.12) and rearrange the resulting expression.

$\displaystyle
M_{**}(z,1,h)=\frac{1}{1-\frac{1}{2}\tilde{\theta}\sigma^{2}z^{2}-\tilde{\delta}\tilde{\theta}z}\quad\tilde{\theta}=\frac{\theta}{1-\frac{1}{2}\theta\sigma^{2}h^{2}-\delta\theta
h}\quad\tilde{\delta}=\delta+h\sigma^{2}$ (4.13)

$M(z,1,h)$ in (4.4) becomes

$\displaystyle M(z,1,{h})=\frac{e^{\mu
z}}{\left(1-\frac{1}{2}\tilde{\theta}\sigma^{2}z^{2}-\tilde{\delta}\tilde{\theta}z\right)^{\alpha}}$
(4.14)

Using the Esscher transform method, the moment generating function for
Variance - Gamma process $Y=\\{Y_{t}\\}_{t\geq 0}$ becomes:

$\displaystyle M(z,t,{h})=
E^{{h}}\left[e^{zY_{t}}\right]=M(z,1,{h})^{t}=\frac{e^{t\mu
z}}{\left(1-\frac{1}{2}\tilde{\theta}\sigma^{2}z^{2}-\tilde{\delta}\tilde{\theta}z\right)^{t\alpha}}\quad
\hbox{with $\tilde{h}_{1}<z<\tilde{h}_{2}$ }$ (4.15)

We have a new Variance - Gamma process with parameter
$(t\mu,\tilde{\delta},\sigma,t\alpha,\tilde{\theta})$ ∎

The Esscher transform method preserves the structure of the five-parameter VG
process; it introduces an addition symmetric parameter ($h\sigma^{2}$) and
inflates the Gamma scale parameter by
$\frac{1}{1-\frac{1}{2}\theta\sigma^{2}h^{2}-\delta\theta h}$ factor.

###  4.2. Variance - Gamma Model: Extended Black-Scholes Formula

###### Corollary 4.2

(Extended Black-Scholes)
Let $r$ a continuously compounded risk-free rate of interest;
$Y=\\{Y_{t}\\}_{t\geq 0}$, a VG Process with parameter $(\mu
t,\delta,\sigma,\alpha t,\theta)$; and $(S(0)e^{X_{T}}-K)^{+}$ the terminal
payoff for a contingent claim with the expiry date $T$.
Then at time $t<T$, the arbitrage price of the European call option with the
strike price $K$ can be written as follows.

$\displaystyle F^{GV}_{call}(S_{t},t)$
$\displaystyle=S(t)\left[1-\hat{F}(log(\frac{K}{S(t)}),\tau,h^{*}+1)\right]-Ke^{-r\tau}\left[1-\hat{F}(log(\frac{K}{S(t)}),\tau,h^{*})\right]$
(4.16)
$\displaystyle\hat{F}(log(\frac{K}{S(t)}),\tau,h^{*})=\int_{-\infty}^{log(\frac{K}{S(t)})}\hat{f}(\xi,\tau,h^{*})d\xi\quad
\hbox{$\hat{f}(\xi,\tau,h^{*})$ in (\ref{eq:l3021})}$ (4.17)

where $\tau=T-t$, and $\hat{F}(k,\tau,h^{*})$ and $\hat{F}(k,\tau,h^{*}+1)$
are the cumulative distribution of VG Model with parameter
$(\tau\mu,\tilde{\delta},\sigma,\tau\alpha,\tilde{\theta})$ and parameter
$(\tau\mu,\tilde{\delta}+\sigma^{2},\sigma,\tau\alpha,\tilde{\theta}e^{\frac{r-\mu}{\alpha}})$
respectively

###### Proof:.

$\displaystyle f(Y_{T},K)= (S(0)e^{Y_{T}}-K)^{+}$
$\displaystyle=S(t)(e^{Y_{\tau}}-k)^{+}\quad\quad
\hbox{$Y_{T}=Y_{\tau}+Y_{t}$}$ $\displaystyle=S(t)g(Y_{\tau})\quad\quad
\hbox{$S(t)=S(0)e^{Y_{t}}$ and $k=\frac{K}{S(t)}$} $

Under the Equivalent Martingale Measure (EMM), $\hat{f}(\xi,\tau,h^{*})$ is
the probability density of VG model with parameter
$(\tau\mu,\tilde{\delta},\sigma,\tau\alpha,\tilde{\theta})$. We note
$k=\frac{K}{S(t)}$

$\displaystyle
S(t)e^{-r\tau}\int_{log(k)}^{+\infty}e^{\xi}\hat{f}(\xi,\tau,h^{*})d\xi$
$\displaystyle=S(t)e^{-r\tau}\int_{log(k)}^{+\infty}e^{\xi}\frac{e^{h^{*}\xi}{f(\xi,\tau)}}{M(h,t)}d\xi\quad
\hbox{$\hat{f}(\xi,\tau,h^{*})$ in (\ref{eq:l3021}) }$
$\displaystyle=S(t)e^{-r\tau}\int_{log(k)}^{+\infty}\frac{e^{(1+h^{*})\xi}{f(\xi,\tau)}}{M(h^{*},t)}d\xi\quad
\hbox{$e^{r\tau}=\frac{M(h^{*}+1,\tau)}{M(h^{*},\tau)}$ in (\ref{eq:l34})}$
$\displaystyle=S(t)\int_{log(k)}^{+\infty}\frac{e^{(1+h^{*})\xi}{f(\xi,\tau)}}{M(1+h^{*},t)}d\xi=S(t)\int_{log(k)}^{+\infty}\hat{f}(\xi,\tau,h^{*}+1)d\xi
$

We can now show the relation in (4.16)

$\displaystyle F^{GV}_{call}(S_{t},t)$
$\displaystyle=S(t)e^{-r\tau}E^{Q}\left[g(X_{\tau})\right]=S(t)e^{-r\tau}\int_{-\infty}^{+\infty}\hat{f}(\xi,\tau,h^{*})(e^{\xi}-k)^{+}dy$
$\displaystyle=S(t)e^{-r\tau}\int_{log(k)}^{+\infty}e^{\xi}\hat{f}(\xi,\tau,h^{*})d\xi-
Ke^{-r\tau}\int_{log(k)}^{+\infty}\hat{f}(\xi,\tau,h^{*})d\xi$
$\displaystyle=S(t)\int_{log(k)}^{+\infty}\hat{f}(\xi,\tau,h^{*}+1)d\xi-
Ke^{-r\tau}\int_{log(k)}^{+\infty}\hat{f}(\xi,\tau,h^{*})d\xi$
$\displaystyle=S(t)\left[1-\hat{F}(log(k),\tau,h^{*}+1)\right]-Ke^{-r\tau}\left[1-\hat{F}(log(k),\tau,h^{*})\right]
$

with

$\displaystyle\hat{F}(log(k),\tau,h^{*})=\int_{-\infty}^{log(k)}\hat{f}(\xi,\tau,h^{*})d\xi\quad\quad\hat{F}(log(k),\tau,h^{*}+1)=\int_{-\infty}^{log(k)}\hat{f}(\xi,\tau,h^{*}+1)d\xi
$

∎

From Theorem 3.3 and Theorem 4.1, $\hat{f}(\xi,\tau,h^{*})$ is the probability
density of the VG model with parameter
$(\tau\mu,\tilde{\delta},\sigma,\tau\alpha,\tilde{\theta})$.

$\displaystyle\hat{f}(\xi,\tau,h^{*})$
$\displaystyle=\frac{1}{\sigma\Gamma(\tau\alpha)\tilde{\theta}^{\tau\alpha}}\int_{0}^{+\infty}\frac{1}{\sqrt{2\pi
v}}e^{-\frac{(y-\mu-\tilde{\delta}v)^{2}}{2v\sigma^{2}}}v^{\tau\alpha-1}e^{-\frac{v}{\tilde{\theta}}}\,dv
\quad\quad \hbox{($\tilde{\delta}$, $\tilde{\theta}$) in ($\ref{eq:l377a}$) }$
(4.18)

Following the same methodology, $\hat{f}(\xi,\tau,h^{*}+1)$ is the probability
density of the VG model with parameter
$(\tau\mu,\tilde{\delta}^{\prime},\sigma,\tau\alpha,\tilde{\theta}^{\prime})$.
we have

$\displaystyle\tilde{\delta}^{\prime}$
$\displaystyle=\tilde{\delta}+\sigma^{2}\quad\quad\tilde{\theta}^{\prime}=\tilde{\theta}e^{\frac{r-\mu}{\alpha}}$
(4.19)

In fact, as in (4.18), we have :

$\displaystyle\tilde{\delta}^{\prime}$
$\displaystyle=\delta+(h^{*}+1)\sigma^{2}=\tilde{\delta}+\sigma^{2}$

And

$\displaystyle\tilde{\theta}^{\prime}$
$\displaystyle=\frac{\theta}{1-\frac{1}{2}\theta\sigma^{2}{(h^{*}+1)}^{2}-\delta\theta{(h^{*}+1)}}=\frac{\theta}{1-\frac{1}{2}\theta\sigma^{2}{h^{*}}^{2}-\delta\theta{h^{*}}}e^{\frac{r-\mu}{\alpha}}=\tilde{\theta}e^{\frac{r-\mu}{\alpha}}$

We have the expression of the probability density,

$\displaystyle\hat{f}(\xi,\tau,h^{*}+1)$
$\displaystyle=\frac{1}{\sigma\Gamma(\tau\alpha)\tilde{\theta}^{\prime\tau\alpha}}\int_{0}^{+\infty}\frac{1}{\sqrt{2\pi
v}}e^{-\frac{(y-\mu-\tilde{\delta}^{\prime}v)^{2}}{2v\sigma^{2}}}v^{\tau\alpha-1}e^{-\frac{v}{\tilde{\theta}^{\prime}}}\,dv
$ (4.20)

####  4.2.1. Equivalent Martingale Measure (EMM) Computation

Under the Equivalent Martingale Measure (EMM), $\hat{f}(\xi,\tau,h^{*})$ is
the probability density of VG model with parameter
$(\tau\mu,\tilde{\delta},\sigma,\tau\alpha,\tilde{\theta})$. The Fourier
transform is:

$\displaystyle\mathscr{F}[\hat{f}](y,\tau,h^{*})$
$\displaystyle=E\left[e^{-iyX_{\tau}}\right]=\left(\frac{e^{-i\mu
y}}{\left(1+\frac{1}{2}\tilde{\theta}\sigma^{2}y^{2}+i\tilde{\delta}\tilde{\theta}y\right)^{\alpha}}\right)^{\tau}\quad
\hbox{see (\ref{eq:l3312})}$ $\displaystyle=\phi(y)^{\tau}=e^{\tau
log(\phi(y))}=e^{-\tau\varphi(-y)}$
$\displaystyle\mathscr{F}[\hat{f}](y,\tau,h^{*})=e^{-\tau\varphi(-y)}\quad\quad\varphi(y)=-i\mu
y+\alpha
log(1+\frac{1}{2}\tilde{\theta}\sigma^{2}y^{2}-i\tilde{\delta}\tilde{\theta}y)
$ (4.21)

$\phi(y)$ and $\varphi(y)$ are, respectively, the Fourier Transform of the
probability density and the characteristic exponent of the VG model with
parameter $(\mu,\tilde{\delta},\sigma,\alpha,\tilde{\theta})$

$\hat{f}(\xi,\tau,h^{*})$ can be written as the inverse Fourier Transform from
(4.21)

$\displaystyle\hat{f}(\xi,\tau,h^{*})=
\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{i\xi
z}\mathscr{F}[\hat{f}](z,\tau,h^{*})dz$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{i\xi
z-\tau\varphi(-z)}dz$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{-i\xi
z-\tau\varphi(z)}dz$

It was shown in [31] that we can have

$\displaystyle\mathscr{F}[\hat{F}](\xi,\tau,h^{*})=\frac{\mathscr{F}[\hat{f}](\xi,\tau,h^{*})}{i\xi}+\pi\mathscr{F}[\hat{f}](0)\delta(\xi)$
(4.22)

Based on (4.22), we deduce

$\displaystyle\hat{F}(\xi,\tau,h^{*})=
\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{i\xi
z}\mathscr{F}[\hat{F}](z,\tau,h^{*})dz=\frac{-1}{2\pi}\int_{-\infty}^{+\infty}\frac{e^{-i\xi
z-\tau\varphi(z)}}{iz}dz+\frac{1}{2} $

We have the probability density and cumulative functions

$\displaystyle\hat{f}(\xi,\tau,h^{*})=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{-i\xi
z-\tau\varphi(z)}dz\quad\quad\hat{F}(\xi,\tau,h^{*})=\frac{-1}{2\pi}\int_{-\infty}^{+\infty}\frac{e^{-i\xi
z-\tau\varphi(z)}}{iz}dz+\frac{1}{2} $ (4.23)

For VG Model in Table 1 [31, 32]: : $\hat{\mu}=0.0848$,
$\hat{\delta}=-0.0577$, $\hat{\sigma}=1.0295$, $\hat{\alpha}=0.8845$,
$\hat{\theta}=0.9378$; and a risk-free rate of interest $6\%$, we have a
Esscher transform parameter $h^{*}=-2.6997$ in (4.6).
$\hat{f}(\xi,\tau,h^{*})$ and $\hat{f}(\xi,\tau,h^{*}+1)$ were computed by the
Fractional Fast Fourier (FRFT) [31] as shown in Fig 6(b) and Fig 7(b).
Using the numerical integration technique, $\hat{f}(\xi,\tau,h^{*})$ and
$\hat{f}(\xi,\tau,h^{*}+1)$ were computed by implementing the following
12-point rule Composite Newton-Cotes Quadrature Formulas [33, 29].

$\displaystyle\hat{f}(\xi,\tau,h)$
$\displaystyle\approx\frac{b}{n}\sum_{p=0}^{\frac{n}{Q}-1}\sum_{j=0}^{Q}W_{j}g(x_{Qp+j},\tau,h)$
(4.24) $\displaystyle g(x,\tau,h)$
$\displaystyle=\frac{1}{\sigma\tilde{\theta}^{\tau\alpha}\Gamma(\tau\alpha)\sqrt{2\pi}}e^{-\frac{(x-\tau\mu-\tilde{\delta}v)^{2}}{2v\sigma^{2}}}v^{\tau\alpha-\frac{1}{2}}e^{-\frac{v}{\tilde{\theta}}}$

In order to compute (4.24) for $h=h^{*}$ and $h=h^{*}+1$, the following
parameter values are used: $a=0$, $b=20$, $Q=12$, $n=5000Q$, $n_{0}=5000$ and
the weights $\\{W_{j}\\}_{0\leq j\leq Q}$ values come from Table 1 [33] and
table 4.1 [29]. The results are shown in Fig 6(a) and Fig 7(a).

(a) Newton-Cote Martingale

(b) FRFT Martingale Measure

(c) Error: FRFT versus Newton

Figure 6. Estimations: $\hat{f}(\xi,\tau,h^{*})$ versus
$\hat{f}(\xi,\tau,h^{*}+1)$, $\tau=0.25$

Both methods produce smooth density functions as shown in Fig 6 and Fig 7. Fig
6(c) and Fig 7(c) provide the estimation error of $\hat{f}(\xi,\tau,h^{*}+1)$.
The Fractional Fast Fourier (FRFT) underestimates the peakedness of the
density function when the timeframe is small ($\tau=0.25$) in Fig 6(c). The
estimation error decreases significantly when the timeframe increases. see Fig
7(c) when $\tau=0.5$.

(a) Newton-Cote Martingale

(b) FRFT Martingale Measure

(c) Error: FRFT versus Newton

Figure 7. Estimations: $\hat{f}(\xi,\tau,h^{*})$ versus
$\hat{f}(\xi,\tau,h^{*}+1)$, $\tau=0.5$

Both methods will be used to compute the arbitrage price of the European call
option.

###  4.3.  Variance - Gamma Model: Generalized Black-Scholes Formula

###### Theorem 4.3

Let $r$ a continuously compounded risk-free rate of interest,
$Y=\\{Y_{t}\\}_{t\geq 0}$, a VG Process with parameter $(\mu
t,\delta,\sigma,\alpha t,\theta)$, and $(S(0)e^{X_{T}}-K)^{+}$, the terminal
payoff for a contingent claim with the expiry date $T$. Then at time $t<T$,
the arbitrage price of the European call option with the strike price $K$ can
be written as follows.

$\displaystyle
F^{GV}_{call}(S_{t},t)=\frac{K}{2\pi}\int_{-\infty+iq}^{+\infty+iq}{\frac{e^{\left(i\xi
log(\frac{S(t)}{K})-\tau(r+\varphi(\xi))\right)}}{i\xi(i\xi-1)}d\xi}$ (4.25)

where $\varphi(z)$ is the characteristic exponent of VG model with parameter
$(\mu,\tilde{\delta},\sigma,\alpha,\tilde{\theta})$ in (4.11), $\tau=T-t$ and
$q<-1$.

###### Proof:.

$\displaystyle (S(0)e^{Y_{T}}-K)^{+}$
$\displaystyle=S(t)(e^{Y_{\tau}}-k)^{+}\quad\quad \hbox{$Y_{T}=Y_{\tau}+Y_{t}$
}$ (4.26) $\displaystyle=S(t)g(Y_{\tau},k)\quad\quad
\hbox{$S(t)=S(0)e^{Y_{t}}$ and $k=\frac{K}{S(t)}$} $ (4.27)

$(S(0)e^{Y_{T}}-K)^{+}$ is the payoff of the call option. The Fourier
transform can be written

$\displaystyle\mathscr{F}[g](y,k)$
$\displaystyle=\int_{0}^{+\infty}e^{-iyx}g(x,k)dx=\int_{0}^{+\infty}e^{-iyx}{(e^{x}-k)^{+}}dx=\int_{log(k)}^{+\infty}e^{-iyx}{(e^{x}-k)}dx$
$\displaystyle=\int_{log(k)}^{+\infty}e^{(1-iy)x}dx-k\int_{log(k)}^{+\infty}e^{-iyx}dx=\frac{1}{1-iy}\left[e^{(1-iy)x}\right]_{ln(k)}^{+\infty}+\frac{k}{iy}\left[e^{-iyx}\right]_{ln(k)}^{+\infty}$
$\displaystyle=\frac{ke^{-iylog(k)}}{iy(iy-1)}\hskip 14.22636pt \hbox{for
$\Im(y)<-1$}$

We have the Fourier transform of call payoff

$\displaystyle\hat{g}(y,k)=\mathscr{F}[g](y,k)=\frac{ke^{-iylog(k)}}{iy(iy-1)}\hskip
14.22636pt \hbox{for $\Im(y)<-1$}$ (4.28)

It is shown in (4.23) and in (4.21) that $\hat{f}(\xi,\tau,h^{*})$ and
$\varphi(y)$ can be written as follows

$\displaystyle\hat{f}(\xi,\tau,h^{*})=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{-i\xi
z-\tau\varphi(z)}dz\quad\quad\varphi(y)=-i\mu y+\alpha
log(1+\frac{1}{2}\tilde{\theta}\sigma^{2}y^{2}-i\tilde{\delta}\tilde{\theta}y)
$ (4.29)

with ($\tilde{\delta}$, $\tilde{\theta}$) defines in (4.11)

$F^{GV}_{call}(S_{t},t)$ is the call function under the Equivalent Martingale
Measure(EMM), and we can now find a good expression of the function

$\displaystyle F^{GV}_{call}(S_{t},t)$
$\displaystyle=e^{-r\tau}E^{h^{*}}\left[(S(0)e^{X_{T}}-K)^{+}\right]=S(t)e^{-r\tau}E^{h^{*}}\left[g(X_{\tau})\right]\quad\quad\quad
\hbox{recall (\ref{eq:l92})}$
$\displaystyle=S(t)e^{-r\tau}\int_{-\infty}^{+\infty}\hat{f}(y,\tau,h^{*})g(y,k)dy=\frac{S(t)}{2\pi}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}e^{-iyz-\tau(r+\varphi(z))}g(y,k)dydz$
$\displaystyle=\frac{S(t)}{2\pi}\int_{-\infty}^{+\infty}e^{-\tau(r+\varphi(z))}\hat{g}(z,k)dz\quad\quad\quad
\hbox{recall (\ref{eq:l101})}$
$\displaystyle=\frac{K}{2\pi}\int_{-\infty+iq}^{+\infty+iq}\frac{exp\left[{-izlog(k)-\tau(r+\varphi(z)))}\right]}{iz(iz-1)}dz\quad\quad\hbox{$\Im(z)\leq
q<-1$ and $k=\frac{K}{S(t)}$}$
$\displaystyle=\frac{K}{2\pi}\int_{-\infty+iq}^{+\infty+iq}\frac{exp\left[{izlog(\frac{S(t)}{K})-\tau(r+\varphi(z))}\right]}{iz(iz-1)}dz$

We have the formula (4.25)

$\displaystyle
F^{GV}_{call}(S_{t},t)=\frac{K}{2\pi}\int_{-\infty+iq}^{+\infty+iq}\frac{exp\left[{izlog(\frac{S(t)}{K})-\tau(r+\varphi(z))}\right]}{iz(iz-1)}dz$
(4.30)

∎

###  4.4. European Option Pricing by Fractional Fast Fourier (FRFT)

####  4.4.1. Parameter $q$ Evaluation

Let us considerate the stock or index price $S=S_{0}e^{Y}$ and the strike
price $K$; it was shown in (4.28) that the Fourier transform of the call
payoff can be written as follows.

$\displaystyle\hat{g}(y,k)=\mathscr{F}[g](y,k)=\frac{ke^{-iylog(k)}}{iy(iy-1)}\hskip
14.22636pt \hbox{for $\Im(y)<-1$}$

We can recover the call payoff from the inverse of Fourier in (4.28)

$\displaystyle\check{g}(x,k)=\frac{1}{2\pi}\int_{-\infty+iq}^{+\infty+iq}e^{iyx}\mathscr{F}[g](y,k)dy\hskip
14.22636pt \hbox{for $q<-1$}$ (4.31)

The payoff in (4.31) depends on the parameter $q$. As shown in Fig 8(a), for
$q=-2$, the inverse of Fourier in (4.31) produces poor results; in fact, the
curve in red fluctuates around real call payoff $(e^{Y}-k)^{+}$. For
$q=-1.002$, the inverse of Fourier over-estimates the call payoff.

(a) $(e^{x}-k)^{+}$ versus $\check{g}(x,k)$ (ATM)

(b) $ER(k,q)$ and $q$ (ATM)

Figure 8. Optimal $q$ parameter

To find $q$ value with a high level of accuracy, we define the error function
($ER(k,q)$) between the real call payoff and the inverse Fourier payoff, with
$k$ (strike price) and the parameter $q$ as inputs.

$\displaystyle
ER(k,q)=\sqrt{\frac{1}{m}\sum_{j=1}^{m}\left[(e^{x_{j}}-k)^{+}-\check{g}(x_{j},k)\right]^{2}}\hskip
14.22636pt \hbox{with $-M\leq x_{j}\leq M$}$ (4.32)

At the money (ATM) option, the strike price $k=1$ and the $ER(k,q)$ can be
analyzed as a function of one variable $q$. Fig 8(b) displays the error
function ($ER$) graph as a function of $q$. $ER$ is a convex function, which
decreases and increases over the interval $]-\infty,-1[$. The section method
was applied to determine $q=-1.0086$, which minimizes $ER(1,q)$.
Fig 9(a) displays the $ER(k,q)$ minimum value as a function of the strike
price $k$; and the correspondent optimal parameter $q$ as a function of the
strike price $k$ (Fig 9(b)). Both graphs display almost a constant function
with respect to the strike price.

(a) $ER(k,q)$ minimum value

(b) optimal parameter $q$

Figure 9. Optimal $q$ parameter and Error $ER(k,q)$ value

####  4.4.2. Calculating the Fourier Integral by FRFT

The method of Fourier transform [24] provides valuable and powerful tools for
option pricing under a class of Lévy processes when the characteristic
function is much simpler than their density function. We compute the call
option’s value on the SPY ETF with the Fractional Fast Fourier (FRFT).
For $x=log(\frac{S(t)}{K})$, $\frac{F^{GV}_{call}(S(t),t)}{K}$ is the price
per one dollar of the strike price. We have

$\displaystyle\frac{F^{GV}_{call}(S_{t},\tau)}{K}$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty+iq}^{+\infty+iq}\frac{exp\left[{izlog(\frac{S(t)}{K})-\tau(r+\varphi(z))}\right]}{iz(iz-1)}dz$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty+iq}^{+\infty+iq}e^{izx}\frac{exp\left[{-\tau(r+\varphi(z))}\right]}{iz(iz-1)}dz$

we assume

$\displaystyle
f(\xi)=\frac{exp\left[{-\tau(r+\varphi(z))}\right]}{i\xi(i\xi-1)}\quad\quad
F(x)=\frac{1}{2\pi}\int_{-\infty+iq}^{+\infty+iq}e^{i\xi x}f(\xi)d\xi$ (4.33)

Following the notation of the Fractional Fast Fourier (FRFT) developed in [31]
(Section 2 and Appendix A.1.)

$\displaystyle F(x_{k})$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty+iq}^{+\infty+iq}e^{i\xi
x_{k}}f(\xi)d\xi=\frac{e^{-qx_{k}}}{2\pi}\int_{-\infty}^{+\infty}e^{i\xi
x_{k}}f(\xi+iq)d\xi$
$\displaystyle\approx\frac{e^{-qx_{k}}}{2\pi}\int_{-\frac{a}{2}}^{\frac{a}{2}}e^{i\xi
x_{k}}f(\xi+iq)d\xi\approx\frac{\gamma}{2\pi}e^{-(q+\frac{m}{2}\beta)x_{k}}G_{k}(f(\xi_{j}+iq)e^{\pi
ijn\delta},-\delta)$ $\displaystyle
F(x_{k})\approx\frac{\gamma}{2\pi}e^{-(q+\frac{m}{2}\beta)x_{k}}G_{k}(f(\xi_{j}+iq)e^{\pi
ijn\delta},-\delta)$ (4.34)

###  4.5. Empirical Analysis

Based on parameter data from VG Model [31, 32]: $\hat{\mu}=0.0848$,
$\hat{\delta}=-0.0577$, $\hat{\sigma}=1.0295$, $\hat{\alpha}=0.8845$,
$\hat{\theta}=0.9378$; $6\%$ risk-free interest rate is added and the Esscher
transform parameter ($h^{*}=-2.6997$) was computed. The VG option pricing will
be calculated across maturity and option moneyness using Extended and
Generalized Black-Scholes formulas. The closed-form Black-Scholes model [19]
was added to the analysis as a benchmark.

$\displaystyle F^{BS}_{call}(S_{t},\tau)$
$\displaystyle=S_{t}N(d_{1})-Ke^{-r\tau}N(d_{2})$ (4.35) $\displaystyle
d_{1}=\frac{Ln(\frac{S_{t}}{K})+(r+\frac{1}{2}\sigma^{*2})\tau}{\sigma^{*}\sqrt{\tau}}\quad
d_{2}$ $\displaystyle=d_{1}-\sigma^{*}\sqrt{\tau}\quad
N(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}\exp\left(-\frac{t^{2}}{2}\
\right)\,dt$

The variance $\sigma^{*2}=0.1848$ is the annualized variance computed from the
daily SPY ETF return variance in [31].
Option Moneyness describes the intrinsic value of an option in its current
state. It indicates whether the option would make money if exercised
immediately. Option moneyness can be classified into three categories: At-The-
Money (ATM) option ($k=\frac{S_{t}}{K}=1$), Out-of-The Money (OTM) option
($k=\frac{S_{t}}{K}<1$), and In-The-Money (ITM) option
($k=\frac{S_{t}}{K}>1$). On August 04, 2021, the SPY ETF market price closed
at $438.98$. We compute the VG call option price on SPY ETF with the spot
price ($S_{0}$) $438.98$. The results are summarised in Table 3.

Table 3. The price of a European call on SPY ETF

Strike Price | Moneyness | BSM | VGM(4.16) | VGM(4.25) | BSM | VGM(4.16) | VGM(4.25) | BSM | VGM(4.16) | VGM(4.25) | BSM | VGM(4.16) | VGM(4.25) | BSM | VGM(4.16) | VGM(4.25) | BSM | VGM(4.16) | VGM(4.25)
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Period ( in year) | 0.0625 | 0.125 | 0.25 | 0.5 | 0.75 | 1
219.49 | 2.00 | 220.31 | 220.28 | 219.86 | 221.13 | 221.10 | 220.48 | 222.76 | 222.72 | 221.71 | 225.98 | 225.93 | 224.14 | 229.15 | 229.13 | 226.53 | 232.27 | 232.26 | 228.88
225.12 | 1.95 | 214.70 | 214.74 | 214.10 | 215.54 | 215.58 | 214.74 | 217.21 | 217.25 | 216.01 | 220.52 | 220.54 | 218.53 | 223.77 | 223.82 | 221.01 | 226.97 | 227.04 | 223.45
231.04 | 1.90 | 208.80 | 208.83 | 208.18 | 209.66 | 209.69 | 208.84 | 211.38 | 211.41 | 210.16 | 214.77 | 214.80 | 212.77 | 218.10 | 218.16 | 215.34 | 221.39 | 221.47 | 217.88
237.29 | 1.85 | 202.58 | 202.53 | 202.10 | 203.47 | 203.42 | 202.79 | 205.23 | 205.18 | 204.16 | 208.71 | 208.67 | 206.86 | 212.13 | 212.13 | 209.53 | 215.51 | 215.53 | 212.16
243.88 | 1.80 | 196.02 | 196.06 | 195.38 | 196.92 | 196.97 | 196.10 | 198.73 | 198.78 | 197.51 | 202.31 | 202.37 | 200.32 | 205.83 | 205.93 | 203.10 | 209.31 | 209.43 | 205.84
250.85 | 1.75 | 189.07 | 189.16 | 188.47 | 190.01 | 190.10 | 189.21 | 191.87 | 191.96 | 190.68 | 195.55 | 195.66 | 193.61 | 199.17 | 199.33 | 196.50 | 202.75 | 202.94 | 199.36
258.22 | 1.70 | 181.72 | 181.81 | 181.36 | 182.69 | 182.78 | 182.13 | 184.60 | 184.69 | 183.66 | 188.39 | 188.52 | 186.70 | 192.12 | 192.30 | 189.71 | 195.81 | 196.03 | 192.70
266.05 | 1.65 | 173.93 | 173.98 | 173.53 | 174.92 | 174.97 | 174.33 | 176.89 | 176.95 | 175.92 | 180.79 | 180.91 | 179.09 | 184.64 | 184.83 | 182.24 | 188.45 | 188.69 | 185.37
274.36 | 1.60 | 165.64 | 165.64 | 165.45 | 166.67 | 166.66 | 166.28 | 168.70 | 168.70 | 167.94 | 172.73 | 172.82 | 171.26 | 176.70 | 176.87 | 174.56 | 180.63 | 180.88 | 177.84
283.21 | 1.55 | 156.83 | 156.75 | 156.56 | 157.88 | 157.81 | 157.43 | 159.98 | 159.91 | 159.18 | 164.14 | 164.21 | 162.66 | 168.25 | 168.42 | 166.14 | 172.33 | 172.59 | 169.60
292.65 | 1.50 | 147.42 | 147.28 | 146.81 | 148.51 | 148.38 | 147.73 | 150.68 | 150.56 | 149.56 | 154.98 | 155.05 | 153.24 | 159.24 | 159.45 | 156.93 | 163.49 | 163.79 | 160.59
302.75 | 1.45 | 137.37 | 137.50 | 136.72 | 138.50 | 138.63 | 137.69 | 140.74 | 140.89 | 139.62 | 145.20 | 145.61 | 143.53 | 149.64 | 150.21 | 147.45 | 154.08 | 154.76 | 151.34
313.56 | 1.40 | 126.60 | 126.45 | 126.29 | 127.77 | 127.64 | 127.31 | 130.09 | 129.99 | 129.36 | 134.73 | 135.00 | 133.53 | 139.39 | 139.85 | 137.71 | 144.06 | 144.63 | 141.86
325.17 | 1.35 | 115.03 | 115.01 | 114.85 | 116.24 | 116.25 | 115.94 | 118.65 | 118.71 | 118.15 | 123.51 | 124.07 | 122.64 | 128.44 | 129.20 | 127.14 | 133.40 | 134.25 | 131.59
337.68 | 1.30 | 102.57 | 102.48 | 102.35 | 103.83 | 103.80 | 103.53 | 106.34 | 106.39 | 105.94 | 111.50 | 112.20 | 110.84 | 116.79 | 117.68 | 115.72 | 122.10 | 123.05 | 120.53
351.18 | 1.25 | 89.11 | 89.15 | 88.69 | 90.42 | 90.55 | 90.00 | 93.08 | 93.33 | 92.69 | 98.68 | 99.71 | 98.12 | 104.44 | 105.61 | 103.48 | 110.16 | 111.34 | 108.71
365.82 | 1.20 | 74.53 | 74.60 | 74.52 | 75.91 | 76.15 | 76.02 | 78.82 | 79.18 | 79.08 | 85.10 | 86.32 | 85.17 | 91.45 | 92.73 | 91.09 | 97.65 | 98.88 | 96.78
381.72 | 1.15 | 58.69 | 59.15 | 58.38 | 60.22 | 60.94 | 60.20 | 63.67 | 64.34 | 63.82 | 70.94 | 72.47 | 70.85 | 77.99 | 79.48 | 77.46 | 84.70 | 86.09 | 83.68
399.07 | 1.10 | 41.51 | 42.09 | 41.76 | 43.56 | 44.29 | 44.08 | 48.03 | 48.30 | 48.53 | 56.55 | 57.74 | 56.68 | 64.32 | 65.45 | 64.03 | 71.52 | 72.55 | 70.80
418.08 | 1.05 | 23.73 | 24.37 | 24.15 | 27.03 | 27.33 | 27.36 | 32.83 | 32.25 | 33.00 | 42.51 | 43.25 | 42.47 | 50.85 | 51.63 | 50.57 | 58.42 | 59.16 | 57.83
438.98 | 1.00 | 8.92 | 6.45 | 6.76 | 13.11 | 11.13 | 11.40 | 19.53 | 18.35 | 18.43 | 29.61 | 29.17 | 29.01 | 38.11 | 38.01 | 37.61 | 45.79 | 45.79 | 45.20
462.08 | 0.95 | 1.64 | 1.19 | 1.27 | 4.40 | 2.94 | 3.02 | 9.62 | 7.41 | 7.50 | 18.70 | 17.17 | 17.09 | 26.72 | 25.85 | 25.50 | 34.12 | 33.67 | 33.01
487.76 | 0.90 | 0.10 | 0.35 | 0.40 | 0.89 | 0.96 | 1.02 | 3.68 | 2.82 | 2.94 | 10.42 | 8.69 | 8.83 | 17.24 | 15.73 | 15.71 | 23.87 | 22.75 | 22.48
516.45 | 0.85 | 0.00 | 0.10 | 0.13 | 0.09 | 0.31 | 0.34 | 1.02 | 1.03 | 1.11 | 4.96 | 3.92 | 4.09 | 10.02 | 8.48 | 8.63 | 15.46 | 13.90 | 13.93
548.73 | 0.80 | 0.00 | 0.03 | 0.04 | 0.00 | 0.10 | 0.12 | 0.19 | 0.37 | 0.41 | 1.94 | 1.66 | 1.78 | 5.12 | 4.19 | 4.37 | 9.10 | 7.80 | 7.96
585.31 | 0.75 | 0.00 | 0.01 | 0.01 | 0.00 | 0.03 | 0.04 | 0.02 | 0.12 | 0.14 | 0.60 | 0.64 | 0.72 | 2.23 | 1.85 | 2.02 | 4.76 | 3.91 | 4.14
627.11 | 0.70 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 0.01 | 0.00 | 0.04 | 0.04 | 0.14 | 0.23 | 0.26 | 0.80 | 0.75 | 0.83 | 2.15 | 1.77 | 1.91
675.35 | 0.65 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 0.01 | 0.02 | 0.07 | 0.09 | 0.22 | 0.27 | 0.31 | 0.81 | 0.72 | 0.81
731.63 | 0.60 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.02 | 0.03 | 0.05 | 0.09 | 0.11 | 0.24 | 0.26 | 0.31
798.15 | 0.55 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 0.01 | 0.01 | 0.02 | 0.03 | 0.06 | 0.08 | 0.10
877.96 | 0.50 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 0.01 | 0.01 | 0.02 | 0.03

(4.16): 12-point rule Composite Newton-Cotes Quadrature (4.25): Fractional
Fourier Transform (FRFT)

The Fractional Fourier Transform (FRFT) algorithm performs poorly for the in-
the-money (ITM) option. The FRFT underprices the VG option for in-the-money
(ITM), whereas the 12-point rule Composite Newton-Cotes Quadrature produces
consistent option pricing results with the Black-Scholes model. Both
algorithms yield consistent results for at-the-money and out-of-the-money
options.
To generalize the analysis and account for a large range of option moneyness
and maturity, The error (4.36) was computed as the difference between VG
option and BS option prices.

$\displaystyle
Error(k,\tau)=\frac{F^{GV}_{call}(S_{t},\tau)}{K}-\frac{F^{BS}_{call}(S_{t},\tau)}{K}\quad\quad
\hbox{$(k=\frac{S_{t}}{K})$}$ (4.36)

Fig 10 graphs the error as a function of the time to maturity ($\tau$) and the
option moneyness ($k$). The spot price ($S_{t}$) is a constant, and the option
moneyness depends on the strike price.

(a) Error(k,$\tau$)

(b) Error(k,$\tau$) (top view)

Figure 10. Combined Effects of time to Maturity ($\tau$) and Option Moneyness
($k=\frac{S_{t}}{K}$)

The Black-Scholes (BS) and VG models produce different option pricing results.
The Black-Scholes model is overpriced for the out-of-the-money (OTM) option
(blue color in Fig 10) and underpriced for the in-the-money (ITM) option (red
color).
The results in Fig 10 are consistent with [28], where the VG pricing was
performed on S&P500 index data. The shape in Fig 10(a) looks similar to Fig 6
in [28] when the Option moneyness variable replaces the strike price. However,
the overpriced Black-Scholes model in bleu color (Fig 10) does not support the
findings [26] that the VG option prices are typically higher than the Black-
Scholes model prices, with the percentage bias rising when the stock gets out-
of-the-money (OTM). One of the limitations of these studies is that the VG
model is symmetric and uses three parameters. The five-parameter VG model
controls the excess kurtosis and the skewness of the daily SPY ETF return
data.

##  5\. Conclusion

In the paper, a $\Gamma(\alpha,\theta)$ Ornstein-Uhlenbeck type process was
used to build a continuous sample path of a five-parameter Variance-Gamma (VG)
process ($\mu,\delta,\sigma,\alpha,\theta$): location ($\mu$), symmetric
($\delta$), volatility ($\sigma$), and shape ($\alpha$) and scale ($\theta$).
The data parameters [31, 32] were used to simulate the gamma process
($\sigma^{2}(t)$) and the continuous sample path of the subordinator process
($\sigma^{2*}(t)$). Both simulations were used as inputs to simulate the VG
process’s continuous sample path. The Lévy density of the VG process was
derived and shown to belong to a KoPoL family of order $\nu=0$, intensity
$\alpha$ and steepness parameters
$\frac{\delta}{\sigma^{2}}-\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}$
and
$\frac{\delta}{\sigma^{2}}+\sqrt{\frac{\delta^{2}}{\sigma^{4}}+\frac{2}{\theta\sigma^{2}}}$.
It was shown that the VG process converges asymptotically in distribution to a
Lévy process driven by a Normal distribution with mean
$(\mu+\alpha\theta\delta)$ and variance
$\alpha(\theta^{2}\delta^{2}+\sigma^{2}\theta)$. The existence of the
Equivalent Martingale Measure (EMM) of the five-parameter VG process was
shown. The EMM preserves the structure of the five-parameter VG process with
an inflated Gamma scale parameter and a constant term adjustment symmetric
parameter. The extended Black-Scholes formula provides the closed form of the
VG option price. The Lévy process generated by the VG model provides the
Generalized Black-Scholes Formula. The daily SPY ETF return data illustrates
the computation of the European option price under the five-parameter VG
process. The 12-point rule Composite Newton-Cotes Quadrature and the
Fractional Fast Fourier (FRFT) algorithms were implemented to compute the
European option price. It results that the FRFT yields inconsistent European
option prices for in-the-money options. The Black-Scholes (BS) and VG models
produce different option pricing results. The Black-Scholes model is
overpriced for out-of-the-money (OTM) options and underpriced for in-the-money
(ITM) options. However, for deep out-of-the-money (OTM) and deep-in-the-money
(ITM) options, Black-Scholes (BS) and VG models yield almost the same option
price.

## References

  * [1] ME Adeosun, SO Edeki, and OO Ugbebor.  On a variance gamma model (vgm) in option pricing: A difference of two gamma processes.  Journal of Informatics and Mathematical Sciences, 8(1):1–16, 2016\.
  * [2] Mohammed Alhagyan, Masnita Misiran, and Zurni Omar.  Discussions on continuous stochastic volatility models.  Global and Stochastic Analysis, 7(1):55–64, 2020.
  * [3] Andrii Andrusiv and Hans-Jürgen Engelbert.  On the minimal entropy martingale measure for lévy processes.  Stochastics, 92(8):1223–1243, 2020.
  * [4] David Applebaum.  Lévy processes and stochastic calculus.  Cambridge university press, 2009.
  * [5] Juan Arias-de Reyna.  On the theorem of frullani.  Proceedings of the American Mathematical Society, 109(1):165–175, 1990.
  * [6] Ole E Barndorff-Nielsen, Jens Ledet Jensen, and Michael Sørensen.  Some stationary processes in discrete and continuous time.  Advances in Applied Probability, 30(4):989–1007, 1998.
  * [7] Ole E Barndorff-Nielsen and Neil Shephard.  Non-Gaussian OU based models and some of their uses in financial economics.  Nuffield College Oxford, 1999.
  * [8] Ole E Barndorff-Nielsen and Neil Shephard.  Modelling by lévy processess for financial econometrics.  Lévy processes, pages 283–318, 2001.
  * [9] Ole E Barndorff-Nielsen and Neil Shephard.  Non-gaussian ornstein–uhlenbeck-based models and some of their uses in financial economics.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(2):167–241, 2001.
  * [10] Ole E Barndorff-Nielsen and Neil Shephard.  Financial volatility, lévy processes and power variation.  Unpublished book, Nuffield College, 2002.
  * [11] Ole E Barndorff-Nielsen and Neil Shephard.  Integrated ou processes and non-gaussian ou-based stochastic volatility models.  Scandinavian Journal of statistics, 30(2):277–295, 2003.
  * [12] Fischer Black and Myron Scholes.  The pricing of options and corporate liabilities.  Journal of Political Economy, 81(3):637–654, 1973.
  * [13] Svetlana Boyarchenko and Sergei Z Levendorskii.  Non-Gaussian Merton-Black-Scholes Theory, volume 9.  World Scientific, 2002.
  * [14] Peter Carr, Hélyette Geman, Dilip B Madan, and Marc Yor.  Stochastic volatility for lévy processes.  Mathematical finance, 13(3):345–382, 2003.
  * [15] R. Cont.  Empirical properties of asset returns: stylized facts and statistical issues.  Quantitative Finance, 1(2):223–236, feb 2001.
  * [16] Hans U Gerber, Elias SW Shiu, et al.  Option pricing by Esscher transforms.  HEC Ecole des hautes études commerciales, 1993.
  * [17] Steven L Heston.  A closed-form solution for options with stochastic volatility with applications to bond and currency options.  The review of financial studies, 6(2):327–343, 1993.
  * [18] John Hull and Alan White.  The pricing of options on assets with stochastic volatilities.  The journal of finance, 42(2):281–300, 1987.
  * [19] John C Hull.  Options futures and other derivatives.  Pearson Education India, 2003.
  * [20] Sato Ken-Iti.  Lévy processes and infinitely divisible distributions.  Cambridge university press, 1999.
  * [21] Maurice George Kendall et al.  The advanced theory of statistics.  The advanced theory of statistics., (2nd Ed), 1946.
  * [22] Steven G Kou.  A jump-diffusion model for option pricing.  Management science, 48(8):1086–1101, 2002.
  * [23] Andreas E Kyprianou.  Fluctuations of Lévy processes with applications: Introductory Lectures.  Springer Science & Business Media, 2014.
  * [24] Cuixiang Li, Huili Liu, Mengna Wang, and Wenhan Li.  The pricing of compound option under variance gamma process by fft.  Communications in Statistics-Theory and Methods, pages 1–15, 2020\.
  * [25] Dilip B Madan, Peter P Carr, and Eric C Chang.  The variance gamma process and option pricing.  Review of Finance, 2(1):79–105, 1998.
  * [26] Dilip B Madan and Frank Milne.  Option pricing with vg martingale components 1.  Mathematical finance, 1(4):39–55, 1991.
  * [27] Kazuhisa Matsuda.  Introduction to merton jump diffusion model.  Department of Economics. The Graduate Center, The City University of New York, 2004.
  * [28] Sharif Mozumder, Ghulam Sorwar, and Kevin Dowd.  Revisiting variance gamma pricing: An application to s&p500 index options.  International Journal of Financial Engineering, 2(02):1550022, 2015\.
  * [29] A. H. Nzokem.  Stochastic and Renewal Methods Applied to Epidemic Models.  PhD thesis, York University , YorkSpace institutional repository, 2020\.
  * [30] A. H. Nzokem and V. T. Montshiwa.  Fitting generalized tempered stable distribution: Fractional fourier transform (frft) approach.  ARXIV.2205.00586[q-fin.ST], 2022.
  * [31] A.H. Nzokem.  Fitting infinitely divisible distribution: Case of gamma-variance model, 2021.
  * [32] A.H. Nzokem.  Gamma variance model: Fractional fourier transform (FRFT).  Journal of Physics: Conference Series, 2090(1):012094, nov 2021\.
  * [33] A.H. Nzokem.  Numerical solution of a gamma - integral equation using a higher order composite newton-cotes formulas.  Journal of Physics: Conference Series, 2084(1):012019, nov 2021\.
  * [34] A.H. Nzokem.  Sis epidemic model: Birth-and-death markov chain approach.  International Journal of Statistics and Probability, 10(4):10–20, July 2021.
  * [35] Philip Protter.  Stochastic integration and differential equations. A new approach, volume 21.  01 2005.
  * [36] Peter Tankov.  Financial modelling with jump processes.  Chapman and Hall/CRC, 2003.
  * [37] GE Uhlenback and LS Ornstein.  On the theory of the brownian motion.  Phys. Rev, 36:823–841, 1930.