# Exchange option pricing under variance gamma-like models

Matteo Gardini Department of Mathematics, University of Genoa, Via Dodecaneso
16146, Genoa, Italy, email<EMAIL_ADDRESS>Piergiacomo Sabino
Quantitative Risk Management, E.ON SE, Brüsseler Platz 1, 45131 Essen,
Germany, email<EMAIL_ADDRESS>of Mathematics and
Statistics, University of Helsinki, P.O. Box 68 FI-00014 Finland, email
<EMAIL_ADDRESS>
###### Abstract

In this article we focus on the pricing of exchange options when the dynamic
of log-prices follows either the well-known variance gamma or the recent
variance gamma++ process introduced in Gardini et al. [19]. In particular, for
the former model we can derive a Margrabe’s type formula whereas, for the
latter one we can write an “integral free” formula. Furthermore, we show how
to construct a general multidimensional versions of the variance gamma++
processes preserving both the mathematical and numerical tractability.

Finally we apply the derived models to German and French energy power markets:
we calibrate their parameters using real market data and we accordingly
evaluate exchange options with the derived closed formulas, Fourier based
methods and Monte Carlo techniques. Keywords: Lévy Processes; Exchange
Options; Energy Markets; Derivative Pricing

##  1 Introduction

Spread options are widely used in many fields of finance, in particular in
energy markets, and their payoff depends on the difference of the value of two
risky underlying assets instead of one. The spread might be between spot and
futures prices, between currencies, interest rates or commodities, among
others. As mentioned, this type of options is particularly frequent in energy
markets and the name changes depending on the commodity one is dealing with.
The most notable examples are the crack, crush, and spark spreads, which
measure profits in the oil, soybean, and gas markets, respectively. The
interested reader can refer to Carmona and Durrleman [10] for a detailed
discussion on this topic. By denoting $S_{i}=\left\\{S_{i}(t);t\geq
0\right\\},i=1,2$, the evolution of the value of two risky assets
respectively, and labelling the strike price with $K$, the payoff of a spread
option at maturity $T$ is given by

$\Phi(T)=\left(S_{2}(T)-S_{1}(T)-K\right)^{+}.$ (1)

The spread option valuation has been widely investigated: the first attempt to
properly compute the fair price of such a contract can be found in Margrabe
[33] where the author finds a closed formula assuming a Black and Scholes [6]
market and $K=0$. Unfortunately, so far, no closed form for the more general
case $K\neq 0$ is known. Nevertheless, many approximate formulas have been
proposed over the years. One which is commonly used has been proposed by Kirk
[25] and only provides a good approximation for small values of $K$. Many
other alternatives are available in the literature under the assumption of a
the Black-Scholes (BS) market, see for instance Bjerksund and Stensland [5]
and Carmona and Durrleman [10]. They can be also extended taking a multi-asset
spread (see for instance S. J. Deng and Zhou [39], Pellegrino and Sabino [37,
38]) or a non-Gaussian market (see for instance Pellegrino [36]). On the other
hand, it is well know that price dynamics present jumps, fat tails, volatility
clusters or even stochastic volatility and many other stylized facts which
cannot be properly modeled by Gaussian processes. One approach to address some
of the issues mentioned above is to rely on Lévy processes which have gained
popularity because of their mathematical tractability also from the numerical
point of view (see Sato [41] and Applebaum [2] for an extensive discussion and
Cont and Tankov [14] for applications to financial markets). On the order
hand, if one leaves the Gaussian framework things might be hard to handle.
Nevertheless, in some cases, such as the variance gamma (VG) process (Madan
and Seneta [32]), the jump diffusion models proposed by Merton [34] and Kou
[26] and the variance gamma++ (VG++) process recently investigated by Gardini
et al. [19], closed formula expressions for European options valuation can be
derived. For more general models or complicated derivative structures one has
to resort to numerical methods, for instance based on Monte Carlo simulations
(Glasserman [20]), PDE’s (Seydel [44] and Cont and Tankov [14]), Fourier
methods (see Carr and Madan [11], Lewis [27], Attari [3], Lord et al. [28]),
or on the knowledge of the expression of the characteristic function (see Fang
[16]). Therefore, the toolkit is rich to properly model single asset dynamics.
Unfortunately, this is not always the case for the multidimensional framework.
Modeling, calibration and pricing with multivariate Lévy processes are
challenging tasks, several attempts are available in the literature: one
example is the copula approach proposed by Cherubini et al. [12] and Cont and
Tankov [14] which sometimes is difficult to handle. Moreover, in a
multidimensional framework, numerical methods based on PDE or on the Fourier
transform gain computational complexity. Finally, the calibration step, which
is essential, might be very unstable. To address these issues Schoutens [42],
Semeraro [43], Luciano and Semeraro [30] and Ballotta and Bonfiglioli [4]
proposed a framework, based on multivariate Brownian subordination, which
leads to a tractable class of multidimensional Lévy models both from a
numerical and a theoretical point of view. However, in such frameworks, closed
formula expressions for spread options are rare (see for instance Cufaro
Petroni and Sabino [15]) and numerical methods, such as those proposed by Choi
[13], Hurd and Zhou [23],Van Belle et al. [46] and Caldana and Fusai [9],
should be used. Nevertheless, in some very particular cases, closed form
expressions for very simple exchange options can be obtained giving a
computational advantage and a useful benchmark to test numerical algorithms.
In addition, they can be used in numerical inversion to guess the so called
“implied parameters” as observed by Carmona and Durrleman [10].

The first contribution of this article is the derivation of an analytic
expression for the exchange options under a bi-dimensional version of the VG
model which improves the results of Hürlimann [24] and Van Belle et al. [46].
The second contribution is the finding of an integral free Margrabe’s style
formula for the VG++ model introduced by Gardini et al. [19] which can be
considered an extension of the VG model. To this end, we derive a
multidimensional versions of the VG++ process, we briefly study its
properties, we provide closed form expressions of their characteristic
function and that of the linear correlation coefficient. Accordingly, we
compare the performance of the new formulas.

The article is organized as follows: in Sections 2 and 3 we recall the
Margrabe formula and extend it to the VG setting. In Section 4 we recall the
notion of self-decomposability and we briefly introduce the VG++ process.
These concepts are preparatory for the derivation of the integral free pricing
formula for exchange options in the VG++ framework presented in Section 5.
Then, in Section 6 we derive multidimensional versions of the VG++ process, we
briefly study its properties, provide closed form expressions of their
characteristic function and that of the linear correlation coefficient.
Section 7 contains the applications to concrete problems: we consider
electricity markets and calibrate the model parameters to real data and price
exchange options written on futures power German and France calendar contracts
using the new closed formulas, Monte Carlo methods and Fourier techniques.
Section 8 concludes the paper and discusses possible extensions along with
further inquires.

##  2 The Margrabe’s formula

In this section we briefly recall some well known results which are
preparatory for the sequel. Consider a market with two risky assets, $S_{1}$
and $S_{2}$ and the money market account $M$, with the following dynamics
under the risk neutral measure $\mathbb{Q}$

$\begin{split}dM(t)=rM(t)dt,\quad
dS_{i}(t)=\mu_{i}S_{i}(t)dt+\sigma_{i}S_{i}(t)dW_{i}^{\mathbb{Q}}(t),\quad
i=1,2,\\\ \end{split}$ (2)

where $W_{1}^{\mathbb{Q}}$ and $W_{2}^{\mathbb{Q}}$ are Brownian motions such
that $\mathbb{E}\left[W_{1}^{\mathbb{Q}}(t)W_{2}^{\mathbb{Q}}(t)\right]=\rho
t$. Consider the pricing problem of an exchange option, whose payoff at
maturity $T$ is given by

$\Phi(T)=\left(S_{2}(T)-S_{1}(T)\right)^{+},$ (3)

namely a spread option with a zero strike price. By using a change of
numéraire technique, presented in Shreve [45], it is easy to show that the
value of this contract $V(0)$ at time $t=0$ is given by

$V(0)=e^{(\mu_{2}-r)T}S_{2}(0)\mathcal{N}(d_{1})-e^{(\mu_{1}-r)T}S_{1}(0)\mathcal{N}(d_{2}),$
(4)

where

$d_{1}=\frac{\log\left(\frac{S_{2}(0)}{S_{1}(0)}\right)+\frac{1}{2}\bar{\sigma}^{2}T}{\bar{\sigma}\sqrt{T}},\quad
d_{2}=d_{1}-\bar{\sigma}\sqrt{T},\quad\bar{\sigma}^{2}=\sigma_{1}^{2}+\sigma_{2}^{2}-2\rho\sigma_{1}\sigma_{2}.$

Note that if $\mu_{1}=\mu_{2}=r$ one gets the well known formula which was
derived by Margrabe [33]. If the value of the strike price $K$ is different
from zero the payoff at time $T$ is given by Equation (1) and no closed
pricing formula is known therefore, numerical techniques or approximations
must be adopted to properly evaluate the derivative.

It is customary to model the asset prices by the so called Lévy exponentiation
framework, presented in Cont and Tankov [14]. The idea is to model the asset
price process $S=\left\\{S(t);t\geq 0\right\\}$ as $S(t)=e^{rt+X(t)}$, where
$X=\left\\{X(t);t\geq 0\right\\}$ is a Lévy Process. If $X$ is a Brownian
motion with drift, we talk about Brownian exponentiation and we reduce to the
original BS framework. In particular, for $T>0$, we assume that

$\displaystyle M(T)=M(0)\exp\left(rT\right),\quad
S_{i}(T)=S_{i}(0)\exp\left(\omega_{i}T+rT+\theta_{i}T+\sigma_{i}W_{i}^{\mathbb{Q}}(T)\right),\quad
i=1,2,$

where $\mathbb{E}\left[W_{1}^{\mathbb{Q}}(t)W_{2}^{\mathbb{Q}}(t)\right]=\rho
t$ and $\omega_{1}$ and $\omega_{2}$ are the drift corrector parameters111In
BS model we have that
$\omega_{i}=-\left(\theta_{i}+\frac{\sigma_{i}^{2}}{2}\right)$, for $i=1,2$.:
they are chosen so that the discounted price processes are martingales.
Applying Itō’s Lemma we obtain the dynamic of $S_{1}$ and $S_{2}$ which is
given by

$\displaystyle\frac{dS_{i}(t)}{S_{i}(t)}$
$\displaystyle=\left(r+\omega_{i}+\theta_{i}+\frac{\sigma_{i}^{2}}{2}\right)dt+\sigma_{i}dW_{i}^{\mathbb{Q}}(t),\quad
i=1,2.$

In order to price exchange options we can use Equation (4) where we take:

$\mu_{i}=r+\omega_{i}+\theta_{i}+\frac{\sigma_{i}^{2}}{2},\quad i=1,2.$

This formula allows us to easily cope with the more general problem of pricing
an exchange option under a time changed stochastic process like, for example,
the VG model introduced by Madan and Seneta [32], which overcomes some limits
of the standard BS approach.

##  3 Margrabe’s formula under the VG model

In this section, we derive a Margrabe style option pricing formula assuming
that the asset log-price dynamic follows a VG process. Let
$G=\left\\{G(t);t\geq 0\right\\}$ be a gamma subordinator with law
$\Gamma\left(\alpha t,\beta\right)$, we model the evolution of the price of
the risky assets under the risk neutral measure as follows

$S_{i}(T)=e^{\omega_{i}T+rT+\theta_{i}G(T)+\sigma_{i}W_{i}(G(T))},\quad
i=1,2,$

where $\mathbb{E}\left[W_{1}(t)W_{2}(t)\right]=\rho t$ and:

$\omega_{i}=\alpha
t\cdot\log\left(1-\frac{\theta_{i}}{\beta}-\frac{\sigma_{i}^{2}}{2\beta}\right).$

Note that the two dynamics share the same subordinator $G$. This choice
ensures that the couple
$\left(\log\left(S_{1}(t)\right),\log\left(S_{2}(t)\right)\right)$ is a Lévy
process. If we use two different subordinators, $G_{1}$ and $G_{2}$, there is
no guarantee that the couple is a Lévy process, as observed in Buchmann et al.
[7], Buchmann et al. [8] and Michaelsen and Szimayer [35].

######  Remark 1.

In financial modeling it is customary to impose that
$\mathbb{E}\left[G(t)\right]=t$, which means that, on average, the stochastic
time runs as fast as the deterministic one. This condition entails some
parameter constrains: for example, if $G(t)\sim\Gamma\left(\alpha
t,\beta\right)$, we have $\mathbb{E}\left[G(t)\right]=\alpha t/\beta$ and
hence $\alpha=\beta$. Nevertheless, we keep our setting as general as possible
but adopt such an assumption in the numerical section.

Under the usual risk-neutral argument, the price of an exchange option at time
$t=0$ is

$\begin{split}V(0)&=e^{-rT}\mathbb{E}^{\mathbb{Q}}\left[\left(S_{2}(T)-S_{1}(T)\right)^{+}\right]\\\
&=\mathbb{E}^{\mathbb{Q}}\left[\mathbb{E}^{\mathbb{Q}}\left[\left.S_{2}(0)e^{\omega_{2}T+\theta_{2}g+\sigma_{2}W_{2}(g)}-S_{1}(0)e^{\omega_{1}T+\theta_{1}g+\sigma_{1}W_{1}(g)}\right]\right|G(T)=g]\right]\\\
&=\mathbb{E}^{\mathbb{Q}}\left[S_{2}(0)e^{\omega_{2}T+\left(\theta_{2}+\frac{\sigma_{2}^{2}}{2}\right)G(T)}\mathcal{N}\left(d_{1}(G(T))\right)-S_{1}(0)e^{\omega_{1}T+\left(\theta_{1}+\frac{\sigma_{1}^{2}}{2}\right)G(T)}\mathcal{N}\left(d_{2}(G(T))\right)\right],\end{split}$

where in the last step, we used the Margrabe’s formula derived in Equation (4)
and where

$\displaystyle d_{1}(g)$
$\displaystyle=\frac{\log\left(\frac{S_{2}(0)}{S_{1}(0)}+\left(\omega_{2}-\omega_{2}\right)T+\left(\theta_{2}-\theta_{1}\right)g+\frac{1}{2}\left(\sigma_{2}^{2}-\sigma_{1}^{2}\right)g+\frac{1}{2}\bar{\sigma}g\right)}{\bar{\sigma}\sqrt{g}},$
(5) $\displaystyle d_{2}(g)$ $\displaystyle=d_{1}(g)-\bar{\sigma}\sqrt{g}.$
(6)

Since the probability density function of $\Gamma\left(\alpha T,\beta\right)$
is given by

$f\left(x;\alpha T,\beta\right)=\frac{\beta^{\alpha T}}{\Gamma\left(\alpha
T\right)}x^{\alpha T-1}e^{-\beta x}\mathbbm{1}_{x>0},$

we get:

$\begin{split}V(0)&=\underbrace{\int_{0}^{\infty}S_{2}(0)e^{\omega_{2}T+\left(\theta_{2}+\frac{\sigma_{2}^{2}}{2}\right)g}\mathcal{N}\left(d_{1}(g)\right)\frac{\beta^{\alpha
T}}{\Gamma\left(\alpha T\right)}g^{\alpha T-1}e^{-\beta g}dg}_{I_{1}}\\\
&-\underbrace{\int_{0}^{\infty}S_{1}(0)e^{\omega_{1}T+\left(\theta_{1}+\frac{\sigma_{1}^{2}}{2}\right)g}\mathcal{N}\left(d_{2}(g)\right)\frac{\beta^{\alpha
T}}{\Gamma\left(\alpha T\right)}g^{\alpha T-1}e^{-\beta
g}dg}_{I_{2}}.\end{split}$ (7)

defining

$\psi\left(a,b;\gamma\right)=\int_{0}^{\infty}\mathcal{N}\left(\frac{a}{\sqrt{u}}+b\sqrt{u}\right)\frac{u^{\gamma-1}}{\Gamma\left(\gamma\right)}e^{-u}du,$

and focusing on the computation of $I_{1}$ we have that

$\begin{split}I_{1}&=S_{2}(0)\beta^{\alpha
T}e^{\omega_{2}T}\int_{0}^{\infty}\frac{1}{\Gamma\left(\alpha
T\right)}g^{\alpha
T-1}e^{-\left(\beta-\left(\theta_{2}+\frac{\sigma_{2}^{2}}{2}\right)\right)g}\\\
&\mathcal{N}\left(\frac{\log\left(\frac{S_{2}(0)}{S_{1}(0)}\right)+\left(\omega_{2}-\omega_{2}\right)T}{\bar{\sigma}\sqrt{g}}+\frac{\theta_{2}-\theta_{1}+\frac{1}{2}\left(\sigma_{2}^{2}-\sigma_{1}^{2}+\bar{\sigma}^{2}\right)}{\bar{\sigma}}\sqrt{g}\right)dg.\end{split}$

Now defining

$A=\beta-\left(\theta_{2}+\frac{\sigma_{2}^{2}}{2}\right),\qquad u=Ag,$

we get

$I_{1}=S_{2}(0)\beta^{\alpha
T}e^{\omega_{2}T}\int_{0}^{\infty}\frac{1}{\Gamma\left(\alpha
T\right)}\frac{1}{A^{\alpha}}u^{\alpha-1}e^{-u}\mathcal{N}\left(\frac{\tilde{a}}{u}+\tilde{b}\sqrt{u}\right),$

where

$\tilde{a}=\frac{\log\left(\frac{S_{2}(0)}{S_{1}(0)}+\left(\omega_{2}-\omega_{1}\right)T\right)\sqrt{A}}{\bar{\sigma}},\qquad\tilde{b}=\frac{\theta_{2}-\theta_{1}+\frac{1}{2}\left(\sigma_{2}^{2}-\sigma_{1}^{2}+\bar{\sigma^{2}}\right)}{\bar{\sigma}\sqrt{A}}.$
(8)

Finally, we obtain:

$I_{1}=S_{2}(0)\frac{\beta^{\alpha T}e^{\omega_{2}T}}{A^{\alpha
T}}\psi\left(\tilde{a},\tilde{b},\alpha T\right).$

The second integral $I_{2}$ can be computed in a similar way. defining

$C=\beta-\left(\theta_{1}+\frac{\sigma_{1}^{2}}{2}\right),$

and

$\tilde{c}=\frac{\log\left(\frac{S_{2}(0)}{S_{1}(0)}+\left(\omega_{2}-\omega_{1}\right)T\right)\sqrt{C}}{\bar{\sigma}},\qquad\tilde{d}=\frac{\theta_{2}-\theta_{1}+\frac{1}{2}\left(\sigma_{2}^{2}-\sigma_{1}^{2}+\bar{\sigma^{2}}\right)}{\bar{\sigma}\sqrt{C}},$

we obtain that

$I_{2}=S_{1}(0)\frac{\beta^{\alpha T}e^{\omega_{1}T}}{C^{\alpha
T}}\psi\left(\tilde{c},\tilde{d},\alpha T\right).$

Finally, the price of the exchange option in the VG model is given by

$V(0)=S_{2}(0)\frac{\beta^{\alpha
T}e^{\omega_{1}T}}{A^{\alpha}}\psi\left(\tilde{a},\tilde{b};\alpha
T\right)-S_{1}(0)\frac{\beta^{\alpha
T}e^{\omega_{1}T}}{C^{\alpha}}\psi\left(\tilde{c},\tilde{d};\alpha T\right).$
(9)

The function $\psi\left(a,b;\gamma\right)$ can be written in terms of the
modified Bessel function of the second type (see Abramowitz and Stegun [1,
pag. 374]) and of the confluent hypergeometric function of two variables
introduced by Humbert [22] hence, it can be efficiently computed, as shown in
Madan et al. [31]. Under the additional hypothesis $S_{1}(0)=S_{2}(0)$ and
$\mu_{1}=\mu_{2}$, Van Belle et al. [46] have shown that a closed form option
pricing formula for exchange options can be obtained in terms of the gamma and
of the hypergeometric functions (see Abramowitz and Stegun [1, pag 556]). It
is easy to check that such a result is a particular case of the more general
formula presented in Equation (9).

In Figure 1 we compare three different methods to compute the exchange option
price under the VG model: a standard Monte Carlo method (see Cont and Tankov
[14]), a Gauss-Kronrod quadrature method to directly integrate Equation (7)
and finally the closed formula in Equation (9). We observe that all methods
produce the same results, as expected. Of course, the direct integration
approach should be implemented carefully since the integrand function is quite
complex. In Table 1 we compare the computational time of both the Monte Carlo
approach with $10^{6}$ simulations and the closed formula. We observe that the
closed formula is approximately one hundred times faster than the Monte Carlo
method and therefore such an approach is the preferable one222The
computational performance has been measured using MATLAB on a PC with an Intel
Core i5-10210U 2.11 GHz processor..

Figure 1: Price of an exchange option with $S_{1}(0)=50$ priced by using quadrature, Monte Carlo simulation and the closed formula. Algorithm | Computational time (s)
---|---
Monte Carlo | $2.60\cdot 10^{-1}$
Closed formula | $3.91\cdot 10^{-3}$
Table 1: Computational time to price an exchange option with the closed
formula and $N_{sim}=10^{6}$ Monte Carlo simulations: the resulting time of
the Monte Carlo is the average of $100$ runs.

##  4 Self-decomposability and the VG++ process

Gardini et al. [19] have studied a new Lévy process called the VG++ process.
Its construction relies upon the notion of self-decomposability which on the
other hand can be also used to built multidimensional Lévy processes including
stochastic delay in the propagation of market news, as discussed in Gardini et
al. [17, 18].

In defining the VG++ process the main idea is to replace the gamma
subordinator by another subordinator strictly related to it which in turn, has
finite activity and the step-wise increasing trajectories. Therefore, seeing
the subordinator as the evolution of the number of transactions, the
probability of having no transaction in a finite period of time will not be in
general zero.

We recall that a random variable $X$ is said to have a self-decomposable law
if for all $a\in\left(0,1\right)$ we can always find two independent random
variables $Y$, with the same law as $X$, and $Z_{a}$ such that in distribution
we have

$X\stackrel{{\scriptstyle d}}{{=}}aY+Z_{a}.$ (10)

We refer to $Z_{a}$ as the $a$-remainder of the sd law. Accordingly, denoting
$\phi_{X}\left(u\right)$ the characteristic function of $X$ and by
$\phi_{Z_{a}}\left(u\right)$ that of $Z_{a}$ we have that

$\phi_{X}\left(u\right)=\phi_{X}\left(au\right)\phi_{Z_{a}}\left(u\right).$
(11)

The law of $Z_{a}$ turns out to be infinitely divisible (see Sato [41,
Proposition 15.5]) and one can construct the associated Lévy process
$Z_{a}=\left\\{Z_{a}(t);t\geq 0\right\\}$ which is actually a subordinator
that can be used for Brownian subordination.

The gamma distribution is a well-known self-decomposable law (see Grigelionis
[21]) and the distribution $\Gamma^{++}\left(a,\alpha,\beta\right)$ of its
$a$-remainder $Z_{a}$ is infinitely divisible. The VG++ process is then
defined as follows

######  Definition 4.1.

Consider a Brownian motion $W=\left\\{W(t);t\geq 0\right\\}$, with drift
$\theta\in\mathbb{R}$, diffusion $\sigma\in\mathbb{R}^{+}$ independent of
$Z_{a}^{++}$. We call the process $X=\left\\{X(t);t\geq 0\right\\}$ defined as

$X(t)=\theta Z_{a}(t)+\sigma W\left(Z_{a}(t)\right),\quad t\geq 0$ (12)

VG++ process.

Such a process has been studied in details in Gardini et al. [19] and it has
several nice properties. It can be written as difference of two $Z_{a}^{++}$
processes with suitable parameters, it is of finite activity and therefore of
finite variation. In particular, the probability density function of $X$ at
time $t>0$ and its characteristic function can be written in closed form. In
addition, the path-simulation is very efficient. Accordingly, all the common
numerical techniques for calibration and pricing can be adopted. Finally, a
numerically convenient integral-free closed form option pricing formula for
European call (and put) options is available.

From an economical point, the VG++ process can be used to model illiquid
markets. Indeed, the value of the parameters $a$ and $\alpha$ are related to
the liquidity activity of the market. Taking the change $\Delta X=X(t)-X(t-1)$
of the log-price over the time interval $\Delta t$, it can be shown that the
probability that the increment equals zero over the time interval $\Delta t$
is strictly larger than zero and, more precisely, it is given by

$P\left(\Delta X=0\right)=a^{\alpha\Delta t},$

since the density of the VG++ process has an atom in zero. This is the main
financial difference from the standard VG process which does imply that non-
zero trading activity takes place in every time interval. On the other hand,
the VG++ model inherits the mathematical tractability of the standard VG
process which is, in any case, recovered when $a$ tends to zero.

It turns out that an integral free closed formula for exchange options can
obtained for the VG++ model as well, which will be detailed in the next
section.

##  5 An integral free formula for exchange options under the VG ++ process

In this section we derive the explicit formula for exchange options under the
VG++ process but first recall some useful relations.

Following Humbert [22], we define the confluent hypergeometric function of two
variables $x$ and $y$ as

$\Phi\left(\alpha,\beta,\gamma;x,y\right)=\frac{\Gamma(\gamma)}{\Gamma(\alpha)\Gamma(\gamma-\alpha)}\int_{0}^{1}u^{\alpha-1}\left(1-u\right)^{\gamma-\alpha-1}\left(1-ux\right)^{-\beta}e^{uy}du,$
(13)

and observe that if $\alpha=n,\beta=1-n$ and $\gamma=1+n$ for
$n\in\mathbb{N}$, Equation (13) becomes:

$\Phi\left(n,1-n,1+n;x,y\right)=n\int_{0}^{1}u^{n-1}\left(1-ux\right)^{n-1}e^{uy}du.$
(14)

It turns out that the integral in Equation (14) can be computed exactly, as it
will be shown later. $\Phi\left(1+n,1-n,2+n;x,y\right)$ can be computed in a
similar way.

Another important formula relative to the Bessel function of the second kind
$K_{\nu}$ when $\nu=n+\frac{1}{2}$ (see Abramowitz and Stegun [1, pag. 444])
is

$\sqrt{\frac{\pi}{2x}}K_{n+\frac{1}{2}}(x)=\left(\frac{\pi}{2}\right)e^{-x}\sum_{k=0}^{n}\left(n+\frac{1}{2},k\right)\left(2x\right)^{-k},$
(15)

where

$\left(n+\frac{1}{2},k\right)=\frac{\left(n+k\right)!}{k!\Gamma\left(n-k+1\right)}.$

Finally, the following trivial identity will be used

$K_{n-\frac{1}{2}}(x)=K_{(n-1)+\frac{1}{2}}(x),\quad n\in\mathbb{N}.$

###  5.1 Computation of $\Phi\left(n,1-n,1+n;x,y\right)$

There are several possible ways to compute Equation (14) analytically.
Observing that

$\begin{split}I_{2}&=\int_{0}^{1}\left(u(1-ux)\right)^{n-1}e^{uy}du=\int_{0}^{1}\left(u(1-ux)\right)^{n-1}\sum_{m=0}^{\infty}\frac{(uy)^{m}}{m!}du\\\
&=\sum_{m=0}^{\infty}\frac{y^{m}}{m!}\int_{0}^{1}u^{m+n-1}\left(1-ux\right)^{n-1}du,\end{split}$
(16)

the problem reduces to compute an integral of the form

$I_{3}=\int u^{p}\left(1-ux\right)^{q}du=I_{3,1}+I_{3,2}+I_{3,3}+C.$

Simple computations show that

$\displaystyle I_{3,1}$
$\displaystyle=\frac{u^{p+1}}{p+1}\left(1-ux\right)^{q},$ $\displaystyle
I_{3,2}$
$\displaystyle=\sum_{j=2}^{q}x^{j-1}\left(1-ux\right)^{q+1-j}\frac{\prod_{k=2}^{j}(q+2-k)}{\prod_{k=1}^{j}(p+j)},$
$\displaystyle I_{3,3}$
$\displaystyle=x^{q}q!\frac{u^{p+q+1}}{(p+1)(p+2)\dots(p+q)(p+q+1)}.$

and hence we obtain an analytical expression for Equation
$\Phi\left(n,1-n,1+n;x,y\right)$.

On the other hand, relying upon the binomial theorem

$\left(x+y\right)^{n}=\sum_{k=0}^{n}\binom{n}{k}x^{n-k}y^{k},$

it follows that

$\begin{split}I_{2}&=\int\left(u-u^{2}x\right)^{n-1}e^{uy}du=\sum_{k=0}^{n-1}\binom{n-1}{k}\int
u^{n-1-k}\left(-u^{2}x\right)^{k}e^{uy}du\\\
&=\sum_{k=0}^{n-1}\left(-1\right)^{k}x^{k}\binom{n-1}{k}\underbrace{\int
u^{n-1+k}e^{uy}du}_{I_{4}}.\end{split}$ (17)

Since:

$\int
u^{m}e^{uy}du=e^{uy}\cdot\left\\{\sum_{k=0}^{m}\left(-1\right)^{k}\binom{m}{k}k!\frac{u^{m-k}}{y^{k+1}}\right\\},$

we can conclude that

$\Phi\left(n,1-n,1+n;x,y\right)=n\sum_{k=0}^{n-1}\left\\{\left(-1\right)^{k}x^{k}\binom{n-1}{k}\cdot
e^{uy}\cdot\left(\sum_{j=0}^{m}\left(-1\right)^{j}\binom{m}{j}j\frac{u^{m-j}}{y^{j+1}}\right)\right\\}.$

###  5.2 Exchange option Pricing

In this section we derive a Margrabe’s style formula for exchange option for
the VG++ model. For $T>0$, we model the two risky asset $S_{1}$,$S_{2}$ under
measure $\mathbb{Q}$ as follows

$\begin{split}S_{i}(T)=S_{i}(0)e^{\omega_{i}T+rT+\theta_{i}Z_{a}(T)+\sigma_{i}W_{i}\left(Z_{a}(T)\right)},\quad
i=1,2,\end{split}$ (18)

where $W_{1}=\left\\{W_{1}(t);t\geq 0\right\\}$ and
$W_{2}=\left\\{W_{2}(t);t\geq 0\right\\}$ are two standard Brownian motions,
$\mathbb{E}\left[W_{1}(t)W_{2}(t)\right]=\rho t$ and
$Z_{a}=\left\\{Z_{a}(t);t\geq 0\right\\}$ is a $\Gamma^{++}\left(a,\alpha
t,\beta\right)$ subordinator, as discussed in Section 4. Observe that we use
the subordinator $Z_{a}$ to time change both Brownian motions $W_{1}$ and
$W_{2}$ as proposed by Schoutens [42]: this guarantees us to remain in a Lévy
framework also when $W_{1}$ and $W_{2}$ are correlated. From Gardini et al.
[19] we recall that the density of
$Z_{a}\sim\Gamma^{++}\left(a,\alpha,\beta\right)$ is given by

$\ h_{a}(x)=a^{\alpha}\delta_{0}(x)+\sum_{n\geq
1}\binom{\alpha+n-1}{n}a^{\alpha}\left(1-a\right)^{n}f_{n,\beta/a}(x)\mathbbm{1}_{x\geq
0},$

where, for $\alpha>0$,

$\binom{\alpha+n-1}{n}=\frac{\alpha\left(\alpha-1\right)\dots\left(\alpha-k+1\right)}{k!},\quad\binom{\alpha}{0}=1,$

and

$f_{n,\beta/a}=\left(\frac{\beta}{a}\right)^{n}x^{n-1}\frac{e^{-\beta
x/a}}{\Gamma\left(n\right)},$

which is the density of an Erlang distribution with parameters
$n\in\mathbb{N}$ and $\beta/a$. Proceeding by conditioning on $Z_{a}(T)$, as
done in the VG case, we have that the value of the exchange option $V(0)$ at
time $t=0$ is given by

$\begin{split}V(0)&=e^{-rT}\mathbb{E}^{\mathbb{Q}}\left[\left(S_{2}(T)-S_{1}(T)\right)^{+}\right]\\\
&=\mathbb{E}^{\mathbb{Q}}\left[S_{2}(0)e^{\omega_{2}T+\left(\theta_{2}+\frac{\sigma_{2}^{2}}{2}\right)Z_{a}(T)}\mathcal{N}\left(d_{1}\left(Z_{a}(T)\right)\right)-S_{1}(0)e^{\omega_{1}T+\left(\theta_{1}+\frac{\sigma_{1}^{2}}{2}\right)Z_{a}(T)}\mathcal{N}\left(d_{2}\left(Z_{a}(T)\right)\right)\right],\end{split}$

where $d_{1}(x)$ and $d_{2}(x)$ are given by Equation (6).

We use the linearity of expected value and observe that we can focus only on
the computation of the first expectation, since the one of the second term is
similar. Omitting the the subscripts we have that

$\begin{split}&\mathbb{E}^{\mathbb{Q}}\left[S(0)e^{\omega T+\theta
Z_{a}(T)+\frac{\sigma^{2}}{2}Z(T)}\mathcal{N}\left(d_{1}\left(Z_{a}(T)\right)\right)\right]=S(0)e^{\omega
T}\int_{0}^{\infty}e^{\theta
x+\frac{\sigma^{2}}{2}x}\mathcal{N}\left(d_{1}(x)\right)h_{a}(x)dx\\\
&=\psi\cdot a^{\alpha T}S(0)e^{\omega T}+S(0)e^{\omega T}\sum_{n\geq
1}\binom{\alpha T+n-1}{n}a^{\alpha
T}\left(1-a\right)^{n}\underbrace{\int_{0}^{\infty}e^{\theta
x+\frac{\sigma^{2}}{2}x}f_{n,\beta/a}(x)\mathcal{N}\left(d_{1}(x)\right)dx}_{I_{1}},\end{split}$
(19)

where

$\psi=\begin{cases}1&\text{if
}\log\left(S_{2}(0)/S_{1}(0)\right)+\left(\omega_{2}-\omega_{1}\right)T>0,\\\
0&\text{otherwise.}\end{cases}$

Focusing on the computation of the term $I_{1}$

$I_{1}=\int_{0}^{\infty}\left(\frac{\beta}{a}\right)^{n}\frac{1}{\Gamma(n)}e^{-\left(\beta/a-\left(\theta+\frac{\sigma^{2}}{2}\right)x\right)}x^{n-1}\mathcal{N}\left(d_{1}(x)\right)dx,$

and defining $A_{j}=\frac{\beta}{a}-\left(\theta+\frac{\sigma^{2}}{2}\right)$
and $u=A_{j}x$ we have that

$I_{1}=\left(\frac{\beta}{aA_{j}}\right)^{n}\int_{0}^{\infty}\mathcal{N}\left(\frac{\tilde{a}}{u}+\tilde{b}\sqrt{u}\right)\frac{u^{n-1}}{\Gamma(n)}e^{-u}du=\left(\frac{\beta}{aA_{j}}\right)^{n}\Psi\left(\tilde{a},\tilde{b};n\right).$

for properly defined values of $\tilde{a}$ and $\tilde{b}$, given by Equations
(8).

In particular, as shown in Madan et al. [31],
$\Psi\left(\tilde{a},\tilde{b};n\right)$ can be expressed in terms of the
modified Bessel function of the second kind $K_{\nu}(z)$ and
$\Phi\left(\alpha,\beta,\gamma;x,y\right)$ where $\alpha$, $\beta$ and
$\gamma$ assume integer values. $\Psi\left(a,b;\gamma\right)$ has the
following expression

$\begin{split}\Psi\left(a,b;\gamma\right)&=\frac{c^{\gamma+1/2}\exp\left[sign\left(a\right)c\right]\left(1+u\right)^{\gamma}}{\sqrt{2\pi}\Gamma\left(\gamma\right)\gamma}\cdot
K_{\gamma+1/2}\left(c\right)\\\
&\Phi\left(\gamma,1-\gamma,1+\gamma,\frac{1+u}{2},-sign\left(a\right)c\left(1+u\right)\right)\\\
&-sign\left(a\right)\frac{c^{\gamma+1/2}\exp\left[sign\left(a\right)c\right]\left(1+u\right)^{\gamma+1}}{\sqrt{2\pi}\Gamma\left(\gamma\right)\left(1+\gamma\right)}\cdot
K_{\gamma-1/2}\left(c\right)\\\
&\Phi\left(1+\gamma,1-\gamma,2+\gamma,\frac{1+u}{2},-sign\left(a\right)c\left(1+u\right)\right)\\\
&+sign\left(a\right)\frac{c^{\gamma+1/2}\exp\left[sign\left(a\right)c\right]\left(1+u\right)^{\gamma}}{\sqrt{2\pi}\Gamma\left(\gamma\right)\gamma}\cdot
K_{\gamma-1/2}\left(c\right)\\\
&\Phi\left(\gamma,1-\gamma,1+\gamma,\frac{1+u}{2},-sign\left(a\right)c\left(1+u\right)\right),\end{split}$

where $c=|a|\sqrt{2+b^{2}}$ and $u=\frac{b}{\sqrt{2+b^{2}}}$. $K_{\nu}(z)$ can
be computed using Equation (15) when $\nu=n\pm 1/2$, whereas, the confluent
hypergeometric function in Equation (13) has an explicit expression as
observed in Section 5.1. Summarizing all the previous results we have the
following proposition.

######  Proposition 5.1.

The price $V(0)$ of an exchange option at time $t=0$ with payoff at maturity
$T$ given by $\Phi\left(T\right)=\left(S_{2}(T)-S_{1}(T)\right)^{+}$ under the
dynamics in Equation (18) is given by

$\begin{split}V(0)&=S_{2}(0)e^{\omega_{2}T}\left(\psi\cdot a^{\alpha
T}+\sum_{n\geq 1}\binom{\alpha T+n-1}{n}a^{\alpha
T}\left(1-a\right)^{n}\Psi\left(\tilde{a},\tilde{b};n\right)\frac{1}{A^{n}}\right)\\\
&-S_{1}(0)e^{\omega_{1}T}\left(\psi\cdot a^{\alpha T}+\sum_{n\geq
1}\binom{\alpha T+n-1}{n}a^{\alpha
T}\left(1-a\right)^{n}\Psi\left(\tilde{c},\tilde{d};n\right)\frac{1}{C^{n}}\right),\end{split}$

where

$A=\frac{a}{\beta}\frac{1}{A_{2}},\qquad
C=\frac{a}{\beta}\frac{1}{A_{1}},\qquad\psi=\begin{cases}1&\text{if
}\log\left(S_{2}(0)/S_{1}(0)\right)+\left(\omega_{2}-\omega_{1}\right)T>0,\\\
0&\text{otherwise.}\end{cases}$

and

$\displaystyle\tilde{a}=\frac{\log\left(\frac{S_{2}(0)}{S_{1}(0)}+\left(\omega_{2}-\omega_{1}\right)T\right)\sqrt{A_{2}}}{\bar{\sigma}},\qquad$
$\displaystyle\tilde{b}=\frac{\theta_{2}-\theta_{1}+\frac{1}{2}\left(\sigma_{2}^{2}-\sigma_{1}^{2}+\bar{\sigma^{2}}\right)}{\bar{\sigma}\sqrt{A_{2}}},$
$\displaystyle\tilde{c}=\frac{\log\left(\frac{S_{2}(0)}{S_{1}(0)}+\left(\omega_{2}-\omega_{1}\right)T\right)\sqrt{A_{1}}}{\bar{\sigma}},\qquad$
$\displaystyle\tilde{d}=\frac{\theta_{2}-\theta_{1}+\frac{1}{2}\left(\sigma_{2}^{2}-\sigma_{1}^{2}+\bar{\sigma^{2}}\right)}{\bar{\sigma}\sqrt{A_{1}}},$

and

$A_{1}=\frac{\beta}{a}-\left(\theta_{1}+\frac{1}{2}\sigma_{1}^{2}\right),\quad
A_{2}=\frac{\beta}{a}-\left(\theta_{2}+\frac{1}{2}\sigma_{2}^{2}\right),\quad\sigma^{2}=\sigma_{1}^{2}+\sigma_{2}^{2}-2\rho\sigma_{1}\sigma_{2}.$

######  Remark 2.

We stress that the pricing formula for the exchange option $V(0)$ in
Proposition 5.1 does not require numerical integration since
$\Psi\left(x,y;n\right)$ can be computed analytically for $n\in\mathbb{N}$ as
shown in Section 5.1.

###  5.3 Numerical test of the option pricing formula

In this section we compare the price of an exchange option calculated with the
exact formula of the previous section to those ones obtained with Monte Carlo
and Fourier-based techniques.

To this end, we choose the following set of parameters: $r=0.01$,
$\sigma_{1}=0.2$, $\sigma_{2}=0.3$, $\theta_{1}=-0.2012$,
$\theta_{2}=-0.1712$, $S_{1}(0)=100$, $S_{2}(0)=105$, $T=1$, $\alpha=2$,
$\rho=0.8$.

We focus on the bi-variate VG++ process for $S_{1}$ and $S_{2}$ defined in
Equation (18) and we define

$Y_{i}(t)=\theta_{i}Z_{a}(T)+\sigma_{i}W_{i}(Z_{a}(T)),\quad i=1,2.\\\ $

The generation of the skeleton of each component is straightforward given the
simulation of $Z_{a}$ which can be easily obtained relying on the method
illustrated in Sabino and Petroni [40]. Namely,

$Z_{a}\stackrel{{\scriptstyle
d}}{{=}}\left\\{\begin{array}[]{ll}\sum_{i=1}^{S}X_{i},&\mbox{when $S>0$}\\\
0,&\mbox{when $S=0$}\end{array}\right.$

where $X_{i}\sim\mathcal{E}(\beta/a)$ iid and
$S\sim\overline{\mathfrak{B}}\left({\alpha,1-a}\right)$. In particular
${Z_{a}}_{|S=s}\sim\Gamma(s,\beta/a)$, when $s>0$.

Several FFT-based methods are available in the literature, among other for
instance the one proposed in Hurd and Zhou [23] or Caldana and Fusai [9] or
that based on the cosine expansions proposed by Pellegrino [36]. All these
approaches require the characteristic function of the log-distribution of
prices under the risk-neutral measure in a closed form. To this end, it is
possible to show that the characteristic function of
$\boldsymbol{Y}(t)=\left(Y_{1}(t),Y_{2}(t))\right)$ is given by

$\phi_{\boldsymbol{Y}(t)}=\phi_{Z(t)}\left(\boldsymbol{\theta}^{T}\boldsymbol{u}+\frac{i}{2}\boldsymbol{u}^{T}\Sigma\boldsymbol{u}\right),$

where $\boldsymbol{\theta}=\left[\theta_{1},\theta_{2}\right]$,
$\boldsymbol{u}=\left[u_{1},u_{2}\right]$ and

$\Sigma=\begin{bmatrix}\sigma_{1}^{2}&\rho\sigma_{1}\sigma_{2}\\\
\rho\sigma_{1}\sigma_{2}&\sigma_{2}^{2}\end{bmatrix}.$

The prices of the exchange options for different values of $a$ are shown in
Figure 2. Here we compare the analytic formula with the Monte Carlo scheme and
the approximate formula for spread options proposed by Caldana and Fusai [9],
which is simpler than the one proposed by Hurd and Zhou [23] since it requires
a single Fourier transform inversion. In Figure 2 we observe that for small
values of $a$ the computation of the exact formula might become problematic.
This is because the exact formula requires the computation of several terms
involving many factorials, which might be numerically unstable in some cases:
some attention should be payed in this case to avoid memory overflow.
Furthermore, even with a smart truncation strategy, the convergence of the
series in Equation (19) can be slow and consequently the computation of the
formula might require some time. For these reasons, in practical purposes, the
pricing techniques based on Fourier transform or on the Monte Carlo method are
preferable. Nevertheless, the exact formula provides a benchmark to test
numerical algorithms and more importantly, could be inverted to infer the so
called “implied parameters”, as commonly done to get the option implied
volatility.

Figure 2: Exchange option pricing in the VG++ model computed by exact formula,
Monte Carlo simulations and FFT method.

##  6 Multidimensional modeling

In this section we go beyond the modeling approach used in the previous
section which relies upon a single subordinator (see Luciano and Schoutens
[29]). Such an approach can be extended using a multidimensional subordinator
while still remaining in a Lévy framework. Successful approaches have been
proposed by Semeraro [43], Luciano and Semeraro [30], Ballotta and Bonfiglioli
[4] and taken up by Gardini et al. [17, 18] to introduce delays in market news
propagation. In order to generalize the results of Section 5 we recall some
properties of the $\Gamma^{++}\left(a,\alpha,\beta\right)$ distribution which
are directly inherited from the gamma law (see Sabino and Petroni [40]).

$\displaystyle Z_{a}\sim\Gamma^{++}\left(a,\alpha,\beta\right),\quad$
$\displaystyle
cZ_{a}\sim\Gamma^{++}\left(a,\alpha,\frac{\beta}{c}\right),\;c>0,$
$\displaystyle Z_{a,i}\sim\Gamma^{++}\left(a,\alpha_{i},\beta\right),\quad$
$\displaystyle\sum_{i=1}^{n}Z_{a,i}\sim\Gamma^{++}\left(a,\sum_{i=1}^{n}\alpha_{i},\beta\right).$

Furthermore, the characteristic function of
$Z_{a}\sim\Gamma^{++}\left(a,\alpha,\beta\right)$ is given by

$\phi_{Z_{a}}(u)=\left(\frac{\beta-iua}{\beta-iu}\right)^{\alpha},\quad
u\in\mathbb{R},.$ (20)

Moreover, for $c\in\mathbb{R}^{+}$, we have that
$\phi_{cZ_{a}}(u)=\phi_{Z_{a}}(cu)$.

###  6.1 Semeraro’s approach

We adapt the approach illustrated in Semeraro [43], in the context of the
VG++process. We consider independent subordinators
$I_{j}=\left\\{I_{j}(t);t\geq 0\right\\},\;j=1,\dots,n$ also independent from
$Z_{a}=\left\\{Z_{a}(t);t\geq 0\right\\}$ and, for $\alpha_{j}\geq
0,\;j=1,\dots,n$, we consider $G_{j}=\left\\{G_{j}(t);t\geq 0\right\\}$
defined as

$G_{j}(t)=I_{j}(t)+\alpha_{j}Z_{a}(t),\quad j=1,\dots,n.$ (21)

By choosing:

$\begin{split}I_{j}(t)&\sim\Gamma^{++}\left(a,A_{j}t,\frac{B}{\alpha_{j}}\right),\quad
A_{j},B>0,\\\ Z_{a}(t)&\sim\Gamma^{++}\left(a,At,B\right),\quad
A,B>0,\end{split}$ (22)

we have that

$G_{j}(t)\sim\Gamma^{++}\left(a,\left(A_{j}+A\right)t,\tfrac{B}{\alpha_{j}}\right).$

Then, given a standard Brownian motion with drift
$\boldsymbol{\mu}=\left(\mu_{1},\dots,\mu_{n}\right)$ and volatility
$\boldsymbol{\sigma}=\left(\sigma_{1},\dots,\sigma_{n}\right)$ we can define a
subordinated Brownian motion
$\boldsymbol{Y}=\left\\{\left(Y_{1}(t),\dots,Y_{n}(t)\right),t\geq 0\right\\}$
as

$Y_{j}(t)=\mu_{j}G_{j}(t)+\sigma_{j}W_{j}\left(G_{j}(t)\right),\quad
j=1,\dots,n.$

The characteristic function of $\boldsymbol{Y}$ at time $t$ can be easily
computed and is given by

$\mathbb{E}\left[e^{i\langle\boldsymbol{u},\boldsymbol{Y}\rangle}\right]=\prod_{j=1}^{n}\phi_{I_{j}(t)}\left(u_{j}\mu_{j}+i\frac{\sigma_{j}^{2}}{2}u_{j}^{2}\right)\cdot\phi_{Z_{a}(t)}\left(\alpha_{j}\left(u_{j}\mu_{j}+i\frac{\sigma_{j}^{2}}{2}u_{j}^{2}\right)\right),$

whereas, the linear correlation coefficient at time $t$ has the following form

$\rho_{Y_{i}(t),Y_{i}(t)}=\frac{\mu_{i}\mu_{j}\alpha_{i},\alpha_{j}Var\left[Z_{a}(t)\right]}{\sqrt{Var\left[Y_{i}(t)\right]Var\left[Y_{j}(t)\right]}}.$
(23)

Linear correlation coefficient is important for the model calibration while
instead, the knowledge of the characteristic function allows us to use pricing
methods based on Fourier techniques. It is worth noting that the model we have
just proposed is not always able to properly catch the correlation observed in
the market, especially if the drift $\mu_{j}$ is small. For this reason,
Luciano and Semeraro [30] introduced an extended version of this model which
aims at better fitting the market correlation even if the drift $\mu_{j}$
vanishes.

###  6.2 Luciano-Semeraro’s approach

In this section we adapt the approach of Luciano and Semeraro [30] which in
contrast, introduces correlated Brownian motions (hereafter labeled LS).
Consider mutually independent subordinators $I_{i}=\left\\{I_{i}(t);t\geq
0\right\\},\;i=1,\dots,n$ and another independent subordinator
$Z_{a}=\left\\{Z_{a}(t);t\geq 0\right\\}$. Define the multidimensional process
$\boldsymbol{Y}=\left\\{Y_{1},\dots,Y_{n}\right\\}$ as follows

$\boldsymbol{Y}\left(t\right)=\left(\begin{array}[]{ll}\mu_{1}I_{1}\left(t\right)+\sigma_{1}W_{1}\left(I_{1}\left(t\right)\right)+\alpha_{1}\mu_{1}Z_{a}\left(t\right)+\sqrt{\alpha_{1}}\sigma_{1}W_{1}^{\rho}\left(Z_{a}\left(t\right)\right)\\\
\dots\\\
\mu_{n}I_{n}\left(t\right)+\sigma_{n}W_{n}\left(I_{n}\left(t\right)\right)+\alpha_{n}\mu_{n}Z_{a}\left(t\right)+\sqrt{\alpha_{n}}\sigma_{n}W_{n}^{\rho}\left(Z_{a}\left(t\right)\right)\end{array}\right),$
(24)

where $I_{j}(t)\;j=1,\dots,n$ and $Z_{a}(t)$ are distributed as in Equation
(22). As shown by Luciano and Semeraro [30], the process $\boldsymbol{Y}$ is a
Lévy process. The expression of its characteristic function and that of its
linear correlation coefficient are known in closed form and they are
respectively given by

$\phi_{\boldsymbol{Y}(t)\left(\boldsymbol{u}\right)}=\prod_{i=1}^{n}\phi_{I_{i}(t)}\left(u_{i}\mu_{i}+\frac{i}{2}\sigma_{i}^{2}u_{i}^{2}\right)\phi_{Z_{a}(t)}\left(\boldsymbol{u}^{T}\boldsymbol{\mu}+\frac{i}{2}\boldsymbol{u}^{T}\Sigma\boldsymbol{u}\right),$

where
$\boldsymbol{\mu}=\left[\alpha_{1}\mu_{1},\dots,\alpha_{n}\mu_{n}\right]$ and

$\Sigma=\begin{bmatrix}\alpha_{1}\sigma_{1}^{2}&\sqrt{\alpha_{1}\alpha_{2}}\sigma_{1}\sigma_{2}\rho_{1,2}&\dots&\sqrt{\alpha_{1}\alpha_{n}}\sigma_{1}\sigma_{n}\rho_{1,n}&\\\
\vdots&\ddots&\ddots&\vdots\\\
\alpha_{n}\sigma_{n}^{2}&\sqrt{\alpha_{n}\alpha_{1}}\sigma_{n}\sigma_{1}\rho_{n,1}&\dots&\sqrt{\alpha_{n}\alpha_{1}}\sigma_{n}\sigma_{1}\rho_{n,1}\end{bmatrix},$

and

$\rho_{Y_{i}\left(t\right),Y_{j}\left(t\right)}=\frac{\mu_{i}\mu_{j}\alpha_{i}\alpha_{j}Var\left[Z_{a}\left(t\right)\right]+\rho\sigma_{i}\sigma_{j}\sqrt{\alpha_{i}\alpha_{j}}\mathbb{E}\left[Z_{a}\left(t\right)\right]}{\sqrt{Var\left[Y_{i}\left(t\right)\right]Var\left[Y_{j}\left(t\right)\right]}}.$

In this case we observe that, even if the drift $\mu_{j}$ vanishes the linear
correlation coefficient can still be different from zero. Furthermore, if the
correlation $\rho=0$ we obtain the same expression of the original model
proposed by Semeraro [43].

###  6.3 Ballotta-Bonfiglioli’s approach

Alternatively to the previous approaches, the dependence can be also
introduced at the level of subordinated Brownian motions as proposed by
Ballotta and Bonfiglioli [4] (hereafter labeled BB) . This approach is
intuitive, but some attention is required in order to guarantee that the sum
of two subordinated Brownian motions is still a subordinated Brownian motion
of the same family. In order to ensure that, one has to impose some
convolutions conditions. Unfortunately, in some case, they can be difficult to
handle, especially from a numerical point of view.

####  6.3.1 Convolution conditions

In this section first we investigate the convolution conditions which
guarantee that a linear combination of two random variables with a VG++ keeps
the same law. Accordingly, we use such conditions in order to build a
multidimensional version of the VG++ process.

Consider three independent subordinators $G_{Y},G_{X}$ and $G_{Z}$ with
$\Gamma^{++}$ law such that at time $t\geq 0$:

$\displaystyle G_{Y}(t)\sim\Gamma^{++}\left(a,A_{y},B_{y}\right),\quad
A_{y},B_{y}>0,$ $\displaystyle
G_{X}(t)\sim\Gamma^{++}\left(a,A_{x},B_{x}\right),\quad A_{x},B_{x}>0,$
$\displaystyle G_{Z}(t)\sim\Gamma^{++}\left(a,A_{z},B_{z}\right),\quad
A_{z},B_{z}>0,.$

and three independent Brownian motions $W_{x},W_{y}$ and $W_{z}$. We define
then, the processes $X,Y$ and $Z$ as follows

$\displaystyle Y(t)$ $\displaystyle=\theta G_{Y}(t)+\sigma
W_{y}\left(G_{Y}(t)\right),\quad\theta\in\mathbb{R},\;\sigma>0,$
$\displaystyle X(t)$ $\displaystyle=\beta G_{X}(t)+\gamma
W_{x}\left(G_{X}(t)\right),\quad\beta\in\mathbb{R},\;\gamma>0,$ $\displaystyle
Z(t)$
$\displaystyle=\beta_{z}G_{Z}(t)+\gamma_{z}W_{z}\left(G_{Z}(t)\right),\quad\beta_{z}\in\mathbb{R},\;\gamma_{z}>0.$

Finally, we define the process $\hat{Y}=X+a_{1}Z$, with $a_{1}\in\mathbb{R}$.
The goal is to find some relations between the parameters such that
$Y(t)\stackrel{{\scriptstyle d}}{{=}}\hat{Y}(t),\;\forall t\geq 0$. Following
the approach of Ballotta and Bonfiglioli [4] and using (20) we equate the
characteristic exponents, obtaining

$A_{y}\log\left(\frac{1-iya/B_{y}}{1-iy/B_{y}}\right)=A_{x}\log\left(\frac{1-ixa/B_{x}}{1-ix/B_{x}}\right)+A_{z}\log\left(\frac{1-iza/B_{z}}{1-iz/A_{z}}\right),$

where

$\displaystyle y$ $\displaystyle=u\theta+iu^{2}\frac{\sigma^{2}}{2},$
$\displaystyle x$ $\displaystyle=u\beta+iu^{2}\frac{\gamma^{2}}{2},$
$\displaystyle z$
$\displaystyle=u\beta_{z}a_{1}+iu^{2}a_{1}^{2}\frac{\gamma_{z}^{2}}{2}.$

The first condition we need to impose is $A_{y}=A_{x}+A_{z}$ which leads to

$\begin{split}A_{y}\log\left(\frac{1-iya/B_{y}}{1-ix/B_{y}}\right)&=A_{y}\log\left(\frac{1-ixa/B_{x}}{1-ix/B_{x}}\right)-A_{z}\log\left(\frac{1-ixa/B_{x}}{1-ix/B_{x}}\right)\\\
&+A_{z}\log\left(\frac{1-iza/B_{z}}{1-iz/A_{z}}\right).\end{split}$

The second one is

$\displaystyle\frac{y}{B_{y}}$ $\displaystyle=\frac{x}{B_{x}},$
$\displaystyle\frac{x}{B_{x}}$ $\displaystyle=a_{1}\frac{z}{B_{z}}.$

From these equations we obtain the following conditions

$\displaystyle\beta$
$\displaystyle=\frac{B_{x}}{B_{z}}a_{1}\beta_{z},\qquad\gamma^{2}=\frac{B_{x}}{B_{z}}\gamma_{z}a_{1}^{2},$
$\displaystyle\theta$
$\displaystyle=\frac{B_{y}}{B_{x}}\gamma^{2}\beta,\qquad\sigma^{2}=\frac{B_{y}}{B_{x}}\gamma^{2},$

and therefore

$\mathbb{E}\left[Y(t)\right]=\frac{A_{y}}{B_{y}}\theta=\beta\frac{A_{x}}{B_{x}}+a_{1}\beta_{z}\frac{A_{z}}{B_{z}}.$

Finally, it is easy to check that

$Var\left[Y(t)\right]=\theta^{2}\left(1-a^{2}\right)\frac{A_{y}}{B_{y}^{2}}+\sigma^{2}\left(1-a^{2}\right)\frac{A_{y}}{B_{y}}.$

Consequently, we have $Y(t)\stackrel{{\scriptstyle
d}}{{=}}\tilde{Y}(t),\;\forall t\geq 0$ and a tool to construct
multidimensional Lévy processes as proposed by Ballotta and Bonfiglioli [4].

####  6.3.2 A multidimensional VG++ process

According to the convolution conditions of the previous section we can
introduce a multi-dimensional version of the VG++ as follows

######  Definition 6.1 (VG++ Ballotta and Bonfiglioli’s model).

Let $Z_{a}$ be a sd subordinator and define the process
$\boldsymbol{Y}=\left\\{\boldsymbol{Y}\left(t\right);t\geq 0\right\\}$ as

$\boldsymbol{Y}\left(t\right)=\left(Y_{1}\left(t\right),\dots,Y_{n}\left(t\right)\right)=\left(X_{1}\left(t\right)+a_{1}Z\left(t\right),\dots,X_{n}\left(t\right)+a_{n}Z\left(t\right)\right),$
(25)

where

  * •

$X_{j}=\left\\{X_{j}\left(t\right);t\geq 0\right\\}$ is a subordinated
Borwnian motion with parameters $\left(\beta_{j},\gamma_{j},\nu_{j}\right)$,
$j=1,\dots,n$, where $\beta_{j}\in\mathbb{R}$ is the drift,
$\gamma_{j}\in\mathbb{R}^{+}$ is the diffusion and $\nu_{j}\in\mathbb{R}^{+}$
is the variance of the subordinator. Let be
$G_{X_{j}}=\left\\{G_{X_{j}}\left(t\right);t\geq 0\right\\}$ the subordinator
of $X_{j}$ and let be $G_{X_{j}}$ be mutually independent for $j=1,\dots,n$.
We define

$X_{j}\left(t\right)=\beta_{j}G_{X_{j}}\left(t\right)+\gamma_{j}W_{x_{j}}\left(G_{X_{j}}\left(t\right)\right),\quad
j=1,\dots,n.$

  * •

Define $Z=\left\\{Z\left(t\right);t\geq 0\right\\}$:

$\displaystyle Z\left(t\right)$
$\displaystyle=\beta_{z}Z_{a}\left(t\right)+\gamma_{z}W_{z}\left(Z_{a}\left(t\right)\right),$

where $W_{z}=\left\\{W_{z}\left(t\right);t\geq 0\right\\}$ is a _BM_ and
$\beta_{z}\in\mathbb{R}$ and $\gamma_{z}\in\mathbb{R}^{+}$.

  * •

$\displaystyle G_{X_{j}}(t)$
$\displaystyle\sim\Gamma^{++}\left(a,A_{x_{j}}t,B_{x_{j}}\right),$
$\displaystyle Z_{a}(t)$
$\displaystyle\sim\Gamma^{++}\left(a,A_{z}t,B_{z}\right).$

  * •

Let $a_{1},\dots,a_{n}\in\mathbb{R}$.

The resulting process is such that $Y_{j}\stackrel{{\scriptstyle
d}}{{=}}\theta_{j}G_{Y_{j}}(t)+\sigma_{j}W_{y}(G_{Y_{j}}(t))$ where
$G_{Y_{j}}(t)\sim\Gamma^{++}\left(a,A_{y_{j}}t,B_{y_{j}}\right)$. In financial
modeling it is customary to require that
$\mathbb{E}\left[G_{Y_{j}}(t)\right]=t$ which can be obtained by imposing
$B_{y_{j}}=\left(1-a\right)A_{y_{j}}$.

After straightforward computations one can show that the characteristic
function of the process $\boldsymbol{Y}$ at time $t$ is given by

$\phi_{\boldsymbol{Y}(t)}\left(\boldsymbol{u}\right)=\prod_{k=1}^{n}\phi_{G_{k}(t)}\left(\beta_{k}u_{k}+\frac{i}{2}u_{k}^{2}\gamma_{k}^{2}\right)\cdot\phi_{Z_{a}(t)}\left(\sum_{k=1}^{n}a_{k}u_{k}\beta_{z}+\frac{i}{2}\left(\sum_{k=1}^{n}a_{k}u_{k}\gamma_{z}\right)^{2}\right).$

Finally, the linear correlation coefficient at time $t$ has the following
expression

$\rho_{Y_{i}(t),Y_{j}(t)}=\frac{a_{i}a_{j}\left(\beta_{z}^{2}\left(1-a^{2}\right)A_{z}/B_{z}^{2}+\left(1-a\right)\gamma_{z}^{2}A_{z}/B_{z}\right)}{\sqrt{\theta_{i}^{2}\left(1-a^{2}\right)\frac{A_{y_{i}}}{B_{y_{i}}^{2}}+\sigma_{i}^{2}\left(1-a^{2}\right)\frac{A_{y_{i}}}{B_{y_{i}}}}\cdot\sqrt{\theta_{j}^{2}\left(1-a^{2}\right)\frac{A_{x_{j}}}{B_{x_{j}}^{2}}+\sigma_{j}^{2}\left(1-a^{2}\right)\frac{A_{y_{j}}}{B_{y_{j}}}}}.$

##  7 Numerical Experiment

In this section we apply the models we have studied so far to real data and
investigate their performance. We focus on the multidimensional versions of
the VG and VG++ process presented in Section 3 and 5 and on those discussed in
Section 6. For each of them we have calibrated the model on market data and
priced exchange options using different numerical techniques.

In order to calibrate our models we need the quoted prices of the derivatives
contracts written on each of the two forward products and their joint
historical time series. The data-set333Data Source: www.eex.com. we rely upon
consists of

  * •

Forward quotations from 4 January 2021 to 19 May 2022 of calendar 2023 power
forward. A forward calendar 2023 contract is a contract to buy or sell a
specific volume of energy in MWh at fixed price for all the hours of 2023.
Calendar power forward in Germany and France are denoted with DEBY and F7BY,
respectively.

  * •

Quotations of European call options on power forward 2023 in Germany and
France with settlement date 19 May 2022. We use strikes in a range of $\pm
15\,[EUR/MWh]$ around the settlement price of the forward contract, i.e. we
excluded deep ITM and OTM options.

  * •

We assume a risk-free rate $r=0.015$.

  * •

The estimated historical correlation between markets log-returns is
$\rho_{mkt}=0.96$.

We have adopted the same two-steps calibration procedure of Luciano and
Semeraro [30]: to this end, it is worthwhile noticing that the marginal
distributions do not depend on the parameters required to model the dependence
structure. The vector of the marginal parameters $\Theta^{*}$ is obtained
solving the following optimization problem:

$\Theta^{*}=\operatorname*{arg\,min}_{\Theta}\sum_{i=1}^{n}\left(C_{i}^{\Theta}\left(K,T\right)-C_{i}\right)^{2},$
(27)

where $C_{i},i=1,\dots,n$ are the values of $n$ quoted vanilla products and
$C_{i}^{\Theta}\left(K,T\right),i=1,\dots,n$ are the relative model prices.
Once $\Theta^{*}$ is found we have to calibrate the remaining parameters for
the dependence structure. Generally derivative contracts written on multiple
underlying assets are not very liquid and market quotes are rarely available
therefore, the vector $\boldsymbol{\eta}^{*}$, that encompasses the dependence
parameters, has been estimated fitting the correlation matrix on historical
data. The expression of the theoretical correlation matrix for the extension
of Luciano and Semeraro [30] and Ballotta and Bonfiglioli [4] in the VG++
framework has been derived in Section 6: the parameters have been chosen so
that the correlation is fitted accordingly to the market.

For the first calibration step we have combined the NLLS approach with the FFT
method proposed by Carr and Madan [11] (the version proposed by Lewis [27]
returns similar results) whereas for the second one, we have used the plain
NLLS method for the minimization of the distance between the theoretical and
the observed correlation coefficient. The calibrated parameters are reported
in Tables 6, 6, 6, 6 and 6 , whereas Table 7 reports the values of the
exchange option computed with Monte Carlo, the FFT method and the closed
formula expression, when available.

Parameter | Value
---|---
$\sigma_{1}$ | 0.84
$\sigma_{2}$ | 0.91
$\rho$ | 0.96
Table 2: _BS_

Parameter | Value
---|---
$\mu_{1}$ | -0.27
$\mu_{2}$ | -0.24
$\sigma_{1}$ | 0.98
$\sigma_{2}$ | 0.92
$\alpha$ | 2.04
$\rho$ | 0.96
Table 3: _VG_

Parameter | Value
---|---
$\mu_{1}$ | -0.28
$\mu_{2}$ | -0.26
$\sigma_{1}$ | 0.98
$\sigma_{2}$ | 0.92
$\alpha$ | 2.43
$a$ | 0.04
$\rho$ | 0.96
Table 4: _VG++_

Parameter | Value
---|---
$\mu_{1}$ | -0.47
$\mu_{2}$ | -0.35
$\sigma_{1}$ | 1.06
$\sigma_{2}$ | 0.96
$\alpha_{1}$ | 1.43
$\alpha_{2}$ | 1.64
$a$ | 0.01
$A$ | 1.43
$B$ | 1224.83
$\rho$ | 0.99
Table 5: _LS-VG++_

Parameter | Value | Parameter | Value
---|---|---|---
$\mu_{1}$ | -0.38 | $\mu_{2}$ | -0.33
$\sigma_{1}$ | 1.025 | $\sigma_{2}$ | 0.96
$A_{y_{1}}$ | 1.5839 | $A_{y_{2}}$ | 1.4439
$A_{x_{1}}$ | 0.141 | $A_{x_{2}}$ | 0.0001
$B_{x_{1}}$ | 12.20 | $B_{x_{2}}$ | 0.006
$\beta_{1}$ | -2.95 | $\beta_{2}$ | -0.00014
$\gamma_{1}$ | 2.84 | $\gamma_{2}$ | 0.062
$a_{1}$ | 4.80 | $a_{2}$ | 4.74
$A_{z}$ | 1.4438 | $B_{z}$ | 62.77
$\beta_{z}$ | -3.15 | $\gamma_{z}$ | 1.34
$a$ | 0.001 |  |
Table 6: _BB-VG++_

Model | Monte Carlo | Exact | FFT
---|---|---|---
BS | 80.04 | 79.92 | 79.91
VG | 81.69 | 81.68 | 81.65
VG++ | 81.62 | 81.28 | 81.15
VG++ LS | 85.18 | - | 85.19
VG++ BB | 83.15 | - | 83.29
Table 7: _Exchange option prices_

From Table 7 we can observe that all models return very similar prices and
that all of them can replicate the observed market log-returns correlation
level with the exception of the LS-VG++ model which attains a slightly lower
level, $\rho_{mod}=0.93$. Since the value of the spread option is very
sensitive to the correlation, in order to avoid mispricing, one should avoid
selecting a model which does not fit it properly.

We remark that we have constructed the VG and VG++ models with a single
subordinator whereas, in the extended models LS-VG++ and BB-VG++ we have
introduced different dependent subordinators. We observe that the option
prices obtained with the different constructions are very similar. This fact
is somehow intuitive from an economical point of view. Indeed, the German and
French power futures markets are extremely intertwined because of the physical
electrical grid connections between the two related countries. We could say
that the stochastic time “runs in the same way” in both markets and hence,
using a single subordinator for price modeling is not a restrictive choice. On
the other hand, there may be other market situations in which multi-
subordination might be the right approach (see Luciano and Semeraro [30] and
Ballotta and Bonfiglioli [4]). In general, when different models lead to
comparable results, it is reasonable to select the simpler one and avoid
unnecessary complexity.

In Section 4 we have briefly mentioned that the parameter
$a\in\left(0,1\right)$ can be interpreted as an index of liquidity. Here we
observe that $a$, and hence the probability of having a zero increment in the
underlying price over the time interval $\Delta t$, is very close to zero in
all cases. During the previous years the options power market has lacked of
liquidity in many European countries but recently such derivative instruments
have become more and more popular leading to an increase in the number of
exchanged trades.

Possible limitations in the construction of the multidimensional version of
the VG++ processes might arise if we consider markets with different level of
liquidities. Indeed, in order to use the summing property of the $\Gamma^{++}$
law, the value of $a\in\left(0,1\right)$ must be the same for all the process
we consider. For this reason, such a construction is reasonable if the
considered markets have approximately the same level of liquidity. This fact
might seem restrictive but, in the most important European power options
markets, almost the same level of liquidity is often observed and hence, the
proposed approach results reasonable.

In this section we have focused only on exchange options, namely spread
options where $K$ is equal to zero. If the strike is not zero, no closed form
expression is know, even in the most simple case of the BS model.
Nevertheless, one can use the approximate formula proposed by Kirk [25] or the
numerical methods proposed by Hurd and Zhou [23], Caldana and Fusai [9] or
Pellegrino and Sabino [37].

##  8 Conclusions

In this article we have firstly discussed the problem of spread option pricing
when the dynamic of the log-prices follows a VG or a VG++ process. Focusing on
exchange options, we have derived closed form pricing formulas when all the
involved Brownian motions are subordinated by a single subordinator. In
particular, for the multidimensional version of the VG++ process of Section 5
we have derived an integral free formula. Accordingly, we have investigated
the numerical performance of the derived formulas, comparing the results with
the ones obtained by Monte Carlo and Fourier techniques. The explicit option
pricing formulas are accurate and efficient, leading to computational benefits
with respect to the Monte Carlo method, especially in the variance gamma case.

As a second contribution, we have proposed more general versions of the
multidimensional VG++ model. Instead of using a common subordinator, these
versions have been obtained by time changing each Brownian motion with
different dependent subordinator relying upon the summing and scaling
properties of the $\Gamma^{++}$ distribution. This new Lévy framework leads to
tractable models, both from a theoretical and numerical point of view. For
instance, the relative expressions of the characteristic function and of the
linear correlation coefficient can be easily derived. Such results are
instrumental for the model calibration and for the option pricing based on FFT
techniques.

Finally, we have applied all these models to the German and French power
futures market. We have calibrated them on vanilla contracts and on historical
futures quotations and accordingly have priced exchange options. Under these
assumptions we have shown that the models with a single subordinator are
sufficiently accurate. This fact has a clear economic interpretation since the
German and French electrical grids are strongly connected and therefore, one
can assume that the business time modeled by the subordinator “flows equally
fast” in both markets. Furthermore, we have observed that the self-
decomposability parameter $a$, that is related to the market liquidity, is
very small and reflects the fact that the option market is sufficiently
liquid, as has been recently observed also from an empirical point of view. In
this work we have mainly focused on the pricing of exchange options
nevertheless, all numerical methods could be easily adapted for the pricing of
spread options. On the other hand, the investigation of approximate formulas
for spread options, for instance, in Kirk’s spirit, under the VG and the VG++
models could be of some interest and this will be the subject of a future
research.

## References

  * Abramowitz and Stegun [1964] M. Abramowitz and I. A. Stegun.  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_.  Dover, New York, ninth dover printing, tenth gpo printing edition, 1964\.
  * Applebaum [2009] D. Applebaum.  _Lévy Processes and Stochastic Calculus_.  Cambridge Studies in Advanced Mathematics. Cambridge University Press, 2 edition, 2009.  doi: 10.1017/CBO9780511809781.
  * Attari [2004] M. Attari.  Option Pricing Using Fourier Transforms: A Numerically Efficient Simplification.  _Derivatives_ , 2004.
  * Ballotta and Bonfiglioli [2013] L. Ballotta and E. Bonfiglioli.  Multivariate Asset Models Using Lévy Processes and Applications.  _The European Journal of Finance_ , 13(22):1320–1350, 2013.
  * Bjerksund and Stensland [2006] P. Bjerksund and G. Stensland.  Closed Form Spread Option Valuation.  Discussion Papers 2006/20, Norwegian School of Economics, Department of Business and Management Science, December 2006.  URL https://ideas.repec.org/p/hhs/nhhfms/2006_020.html.
  * Black and Scholes [1973] F. Black and M. Scholes.  The Pricing of Options and Corporate Liabilities.  _Journal of Political Economy_ , 81(3):637–654, 1973.
  * Buchmann et al. [2017] B. Buchmann, B. Kaehler, R. Maller, and A. Szimayer.  Multivariate Subordination Using Generalised Gamma Convolutions with Applications to Variance Gamma Processes and Option Pricing.  _Stochastic Processes and their Applications_ , 127(7):2208–2242, 2017.
  * Buchmann et al. [2020] B. Buchmann, K. Lu, and D. Madan.  Self-Decomposability of Variance Generalised Gamma Convolutions.  _Stochastic Processes and their Applications_ , 130(2):630–655, 2020.  doi: 10.1016/j.spa.2019.02.012.
  * Caldana and Fusai [2016] R. Caldana and G. Fusai.  A General Closed-Form Spread Option Pricing Formula.  _Journal of Banking & Finance_, 12(37):4863–4906, 2016.
  * Carmona and Durrleman [2003] R. Carmona and V. Durrleman.  Pricing and Hedging Spread Options.  _SIAM Review_ , 45:627–685, 2003.
  * Carr and Madan [1999] P. Carr and D.B. Madan.  Option Valuation Using the Fast Fourier Transform.  _Journal of Computational Finance_ , 2:61–73, 1999.
  * Cherubini et al. [2013] U. Cherubini, E. Luciano, and Vecchiato V.  _Copula Methods in Finance_.  Wiley Finance, 2013.
  * Choi [2018] J. Choi.  Sum of all Black–Scholes–Merton Models: An Efficient Pricing Method for Spread, Basket, and Asian Options.  _Journal of Futures Markets_ , 38(6):627–644, 2018.  doi: https://doi.org/10.1002/fut.21909.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/fut.21909.
  * Cont and Tankov [2003] R. Cont and P. Tankov.  _Financial Modelling with Jump Processes_.  Chapman and Hall, 2003.
  * Cufaro Petroni and Sabino [2018] N. Cufaro Petroni and P. Sabino.  Pricing Exchange Options with Correlated Jump Diffusion Processes.  _Quantitative Finance_ , 0(0):1–13, 2018.  doi: 10.1080/14697688.2017.1423371.
  * Fang [2010] F. Fang.  _The COS Method: An Efficient Fourier Method for Pricing Financial Derivatives_.  PhD thesis, Tecnical University of Delft, 2010.
  * Gardini et al. [2021a] M. Gardini, P. Sabino, and E. Sasso.  Correlating Lévy Processes with Self-Decomposability: Applications to Energy Markets.  _Decision in Economics and Finance_ , 2021a.
  * Gardini et al. [2021b] M. Gardini, P. Sabino, and E. Sasso.  A Bivariate Normal Inverse Gaussian Process with Stochastic Delay: Efficient Simulations and Applications to Energy Markets.  _Applied Mathematical Finance_ , 28(2):178–199, 2021b.  doi: 10.1080/1350486X.2021.2010106.  URL https://doi.org/10.1080/1350486X.2021.2010106.
  * Gardini et al. [2022] M. Gardini, P. Sabino, and E. Sasso.  The Variance Gamma++ Process and Applications to Energy Markets.  _Applied Stochastic Models in Business and Industry_ , 38(2):391–418, 2022.  doi: https://doi.org/10.1002/asmb.2666.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/asmb.2666.
  * Glasserman [2004] P. Glasserman.  _Monte Carlo Methods in Financial Engineering_.  Springer-Verlag New York, 2004.
  * Grigelionis [2003] B. Grigelionis.  On the Self-Decomposability of Euler’s Gamma Function.  _Lithuanian Mathematical Journal_ , 43(3):295–305, 2003.
  * Humbert [1922] P. Humbert.  Ix.—the Confluent Hypergeometric Functions of Two Variables.  _Proceedings of the Royal Society of Edinburgh_ , 41:73–96, 1922.  doi: 10.1017/S0370164600009810.
  * Hurd and Zhou [2009] T.R. Hurd and Z. Zhou.  A Fourier Transform Method for Spread Option Pricing.  https://arxiv.org/pdf/0902.3643.pdf, 2009.
  * Hürlimann [2013] W. Hürlimann.  Margrabe Formulas for a Simple Bivariate Exponential Variance-Gamma Price Process (i) Theory.  _International Journal of Scientific and Innovative Mathematical Research_ , 1(1):33–44, 2013.
  * Kirk [1995] E. Kirk.  Correlations in the Energy Markets.  _In Managing Energy Price Risk, Risk Publications and Enron_ , pages 71–78, 1995.
  * Kou [2002] S. G. Kou.  A Jump-Diffusion Model for Option Pricing.  _Manage. Sci._ , 48(8):1086–1101, August 2002\.  ISSN 0025-1909.
  * Lewis [2001] A. Lewis.  A Simple Option Formula for General Jump-Diffusion and Other Exponential Lévy Processes.  available at http://optioncity.net/pubs/ExpLevy.pdf, 2001.
  * Lord et al. [2007] R. Lord, F. Fang, F. Bervoets, and C.W. Oosterlee.  The CONV Method for pricing options.  _Proceedings in Applied Mathematics & Mechanics_, 7, 2007.
  * Luciano and Schoutens [2006] E. Luciano and W. Schoutens.  A Multivariate Jump-driven Financial Asset Model.  _Quantitative Finance_ , 6(5):385–402, 2006.  URL https://doi.org/10.1080/14697680600806275.
  * Luciano and Semeraro [2010] E. Luciano and P. Semeraro.  Multivariate Time Changes for Lévy Asset Models: Characterization and Calibration.  _Journal of Computational and Applied Mathematics_ , 233(1):1937–1953, 2010.
  * Madan et al. [1999] D. Madan, P. Carr, M. Stanley, and E. Chang.  The Variance Gamma Process and Option Pricing.  _Review of Finance_ , 2, 11 1999.  doi: 10.1023/A:1009703431535.
  * Madan and Seneta [1990] D. B. Madan and E. Seneta.  The Variance Gamma (V.G.) Model for Share Market Returns.  _The Journal of Business_ , 63(4):511–524, 1990\.
  * Margrabe [1978] W. Margrabe.  The Value of an Option to Exchange One Asset for Another.  _Journal of Finance_ , 33(1):177–186, 1978.
  * Merton [1976] R.C. Merton.  Options Pricing when Underlying Shocks are Discontinuous.  _Journal of Financial Economics_ , 3:125–144, 1976.
  * Michaelsen and Szimayer [2018] M. Michaelsen and A. Szimayer.  Marginal Consistent Dependence Modelling using Weak Subordination for Brownian Motions.  _Quantitative Finance_ , 18(11):1909–1925, 2018\.
  * Pellegrino [2016] T. Pellegrino.  A General Closed Form Approximation Pricing Formula for Basket and Multi-Asset Spread Options.  _The Journal of Mathematical Finance_ , 6(5):944–974, 2016.
  * Pellegrino and Sabino [2014a] T. Pellegrino and P. Sabino.  On the Use of the Moment-matching Technique for Pricing and Hedging Multi-asset Spread Options.  _Energy Economics_ , 45:172–185, 2014a.
  * Pellegrino and Sabino [2014b] T. Pellegrino and P. Sabino.  Pricing and Hedging Multiasset Spread Options Using a Three-dimensional Fourier Cosine Series Expansion Method.  _The Journal of Energy Markets_ , 7(2):71–92, 2014b.
  * S. J. Deng and Zhou [2008] M. Li S. J. Deng and J. Zhou.  Closed-form Approximations for Spread Option Prices and Greeks.  _Journal of Derivatives_ , 15(3):58–80, 2008\.  doi: 10.3905/jod.2008.702506.
  * Sabino and Petroni [2020] P. Sabino and N. Cufaro Petroni.  Gamma-related Ornstein–Uhlenbeck Processes and their Simulation.  _Journal of Statistical Computation and Simulation_ , 0(0):1–26, 2020.  doi: 10.1080/00949655.2020.1842408.
  * Sato [1999] K. Sato.  _Lévy Processes and Infinitely Divisible Distributions_.  Cambridge U.P., Cambridge, 1999.
  * Schoutens [2003] W. Schoutens.  _Lévy Processes in Finance: Pricing Financial Derivatives_.  John Wiley and Sons Inc, 2003.
  * Semeraro [2008] P. Semeraro.  A Multivariate Variance Gamma Model For Financial Applications.  _International Journal of Theoretical and Applied Finance_ , 11(1):1–18, 2008.
  * Seydel [2004] R. Seydel.  _Tools for Computational Finance_.  Universitext (1979). Springer, 2004.  ISBN 9783540406044.
  * Shreve [2004] S.E. Shreve.  _Stochastic Calculus for Finance Volume 2_.  Springer, 2004.
  * Van Belle et al. [2019] J. Van Belle, S. Vanduffel, and J. Yao.  Closed-form Approximations for Spread Options in Lévy Markets.  _Applied Stochastic Models in Business and Industry_ , 35(3):732–746, 2019.  doi: https://doi.org/10.1002/asmb.2391.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/asmb.2391.