# Pricing Energy Derivatives in Markets Driven by Tempered Stable and CGMY
Processes of Ornstein-Uhlenbeck Type

Piergiacomo<EMAIL_ADDRESS>222The views, opinions,
positions or strategies expressed in this article are those of the author and
do not necessarily represent the views, opinions, positions or strategies of,
and should not be attributed to E.ON SE.
Quantitative Risk Management
E.ON SE
Brüsseler Platz 1, 45131 Essen, Germany

###### Abstract

In this study we consider the pricing of energy derivatives when the evolution
of spot prices follows a tempered stable or a CGMY driven Ornstein-Uhlenbeck
process. To this end, we first calculate the characteristic function of the
transition law of such processes in closed form. This result is instrumental
for the derivation of non-arbitrage conditions such that the spot dynamics is
consistent with the forward curve. Moreover, based on the results of Cufaro
Petroni and Sabino [21], we also conceive efficient algorithms for the exact
simulation of the skeleton of such processes and propose a novel procedure
when they coincide with compound Poisson processes of Ornstein-Uhlenbeck type.
We illustrate the applicability of the theoretical findings and the simulation
algorithms in the context of the pricing different contracts namely, strips of
daily call options, Asian options with European style and swing options.
Finally, we present an extension to future markets.

_Keywords_ : Lévy-driven Ornstein-Uhlenbeck Processes; CGMY process; Tempered
Stable Distributions; Exact simulation; Energy Markets; Derivative Pricing

##  1 Introduction

Most energy and commodity markets exhibit seasonality, mean-reversion high
volatilities and occasional distinctive price spikes, which results in demand
for derivative products which protect the holder against high prices. In
equity markets there is a clear evidence that asset returns are not Gaussian
and it is common practice to rely on Lévy processes, other than the Brownian
motion, in order to capture heavy-tails and jumps of the log-prices. Several
empirical studies (see for instance Carr and Crosby [14]) have shown that CGMY
Lévy processes introduced by Carr et al. [15], named after its authors, are a
valuable alternative. Moreover, the class of such processes is quite flexible
and also encompasses Variance Gamma processes introduced in Madan and Seneta
[39] and is on the other hand, a special case of the wider class of bilateral
tempered stable processes (see Küchler and Tappe [33]).

Commodity and energy markets however, exhibit mean-reversion which cannot be
described by plain Lévy processes but rather by Lévy-driven Ornstein-Uhlenbeck
(OU) processes. Based on these observations, in this study we consider the
pricing of energy derivatives assuming that the spot price is driven by CGMY
and tempered stable processes of OU type. The first contribution consists in
the derivation of the closed formula of the characteristic function of the
transition law of these processes that is instrumental to find non-arbitrage
conditions. It also gives the fundamental ingredient to calculate the price of
financial derivatives with FFT-based methods.

Based on the results of Cufaro Petroni and Sabino [21], the second
contribution is the derivation of the exact methods for the simulation of the
skeleton of CGMY and bilateral tempered stable processes processes of OU type
with finite variation. Particular emphasis is given to the case when such
processes consists of compound Poisson processes of OU type.

The calibration and the model selection is not the focus of this study,
instead we rather illustrate the theoretical applicability of our findings and
the proposed simulation algorithms to the pricing of a few energy derivative
contracts. As a first application, we consider the pricing of a daily strip of
call options on the day-ahead spot price driven by tempered stable OU
processes using the FFT technique of Carr and Madan [16]. As mentioned, the
parameters calibration is not the focus of this paper, nevertheless such a
task can be easily accomplished combining the FFT pricing with an optimization
to minimize the difference between the available option prices and the
theoretical ones. Secondly, we consider the pricing of Asian options depending
on the day-ahead spot price described by a CGMY-driven OU process via Monte
Carlo simulations where we also highlight the differences between our exact
simulation schemes and the standard approximation procedures. The last example
consists in pricing swing options with the modified version of the Least-
Squares Monte Carlo method detailed in Boogert and C. de Jong [12, 13] using
market models based on CGMY-driven OU process that coincide with compound
Poisson processes of OU type. Finally, we show that our results are not
restricted to OU processes and to the modeling of spot prices. Indeed, in the
spirit of Benth et al. [9], Latini et al. [34] and Piccirilli et al. [41],
they can be adapted to capture the Samuelson effect and different implied
volatility profiles displayed by options futures.

The paper is organized as follows. Section 2 introduces tempered stable and
CGMY processes and the general results relatively to Lévy-driven OU processes.
In Section 3 we focus on the classical tempered stable and CGMY processes of
OU type with finite variation and derive the characteristic function of their
transition law. In Section 4 we present the algorithms for the simulation of
the skeleton of the processes under study and we focus on the case of compound
Poisson processes of OU type. We also present numerical experiments
demonstrating their efficiency. The application of these results is
illustrated in Section 5 in the context of the pricing of energy derivative
contracts, namely daily strips of call options, Asian options with European
exercise, swing options written on the day-ahead spot price and futures.
Finally Section 6 concludes the paper with an overview of future inquiries and
possible further applications.

###  1.1 Notation

Before proceeding, we introduce some notation and shortcuts that will be used
throughout the paper. We write $\Gamma(\alpha,\beta)$ to denote the gamma
distribution with shape parameter $\alpha>0$ and rate parameter $\beta>0$.
Moreover, we write $\mathcal{U}([0,1])$ to denote the uniform distribution in
$[0,1]$ and $\mathcal{P}(\lambda)$ to denote the Poisson distribution with
parameter $\lambda>0$. We use the shortcuts id and sd for _infinitely
divisible_ and _self-decomposable_ distributions, respectively. We use the
shortcut rv for _random variable_ and iid for _independently and identically
distributed_ whereas we use chf, lch cgf and pdf as shortcuts for
_characteristic function_ , _logarithmic characteristic_ , _cumulant
generating function_ and _density function_ , respectively.

##  2 Preliminaries

Take a Lévy process $L(\cdot)$ of classic tempered stable type namely, with
Lévy measure having density

$\nu(x)=\nu_{p}(x)+\nu_{x}(x)=c_{p}\,\frac{e^{-\beta_{p}\,x}}{x^{1+\alpha_{p}}}\mathds{1}_{x\geq
0}+c_{n}\,\frac{e^{\beta_{n}\,x}}{|x|^{1+\alpha_{n}}}\mathds{1}_{x<0}$ (1)

where $c_{p},c_{n},\beta_{p},\beta_{n}$ are all positive numbers,
$\alpha_{p}<2$ and $\alpha_{n}<2$. Hereafter we will denote with
$\mathcal{BCTS}(\alpha_{p},\alpha_{n},\beta_{p},\beta_{n},c_{p}\,t,c_{n}\,t)$
the law of $L(t)$.

Different applications of such a process can be found among other in Koponen
[32], Carr et al. [15] Poirot and Tankov [42] and Ballotta and Kyriakou [2].
In particular, the model introduced in Carr et al. is named CGMY and assumes
$C=c_{p}=c_{n}$, $Y=\alpha_{p}=\alpha_{n}$, $G=\beta_{n}$ and $M=\beta_{p}$
from the names of the authors. It can also be proven that a classic tempered
stable process is a time-changed Brownian motion provided $c_{p}=c_{n}$ and
$\alpha_{p}=\alpha_{n}=\alpha>-1$ (see Cont and Tankov [18] Proposition 4.1).
For sake of completeness, it is worthwhile mentioning that we are referring to
classic tempered stable processes because different processes can be
constructed applying an alternative tempering function rather than the
exponential function used in (1) (see Rosiński [44]). An overview of such
processes, named general tempered stable processes, can be found in Grabchak
[24].

In the following we consider the subset of classic tempered stable processes
with finite variation for which it holds $\alpha_{p}<1$ and $\alpha_{n}<1$ and
in particular when $\alpha_{p}<0$ and $\alpha_{n}<0$ the subset consists of
Poisson processes (see Cont and Tankov [18]). Due to the fact that any process
of finite variation can be seen as the difference of two independent
subordinators, the process $L(\cdot)$ can be written as
$L(t)=L_{p}(t)-L_{n}(t)$ where $L_{p}(\cdot)$ and $L_{n}(\cdot)$ are two
classic tempered stable subordinators with Lévy densities $\nu_{p}(x)$ and
$\nu_{n}(-x)$, respectively. In the following we will dub classic tempered
stable subordinators with CTS, whereas the full bilateral case will be denoted
with BCTS. Moreover, we will denote the law of a CTS subordinator at time $t$
with $\mathcal{CTS}(\alpha,\beta,c\,t)$.

Consider now an Ornstein-Uhlenbeck (OU) process $X(\cdot)$ solution of the
stochastic differential equation

$\displaystyle dX(t)$ $\displaystyle=$ $\displaystyle-bX(t)dt+dL(t)\quad\qquad
X(0)=X_{0}\quad\bm{P}\hbox{-\emph{a.s.}}\qquad b>0$ (2)

namely,

$\displaystyle X(t)$ $\displaystyle=$ $\displaystyle
X_{0}\,e^{-bt}+Z(t)\qquad\quad Z(t)=\int_{0}^{t}e^{-b\,(t-s)}dL(s).$ (3)
$\displaystyle Z(t)$ $\displaystyle=$ $\displaystyle
Z_{p}(t)-Z_{n}(t)\qquad\quad
Z_{d}(t)=\int_{0}^{t}e^{-b\,(t-s)}dL_{d}(s),d\in\\{p,n\\}.$ (4)

Following the convention in Barndorff-Nielsen and Shephard [5], $X(\cdot)$ is
then named OU-BCTS process or, if the above parameter constrain holds, OU-CGMY
process.

There is a close relation between the concept of self-decomposability and the
theory of Lévy driven OU processes, indeed as observed in Barndorff-Nielsen et
al. [4], the solution process (3) is stationary if and only if its
chf$\varphi_{X}(u,t)$ is constant in time and steadily coincides with the chf
$\overline{\varphi}_{X}(u)$ of the sd invariant initial distribution that
turns out to be decomposable according to

$\overline{\varphi}_{X}(u)=\overline{\varphi}_{X}(u\,e^{-b\,t})\varphi_{Z}(u,t)$

where now, at every given $t$, $\varphi_{Z}(u,t)=e^{\psi_{Z}(u,t)}$ denotes
the idchf of the rv $Z(t)$ in (3) and $\psi_{Z}(u,t)$ its lch. We remark that
the process $Z(\cdot)$ is not a Lévy process but rather an additive process.

We recall here that a law with chf $\eta(u)$ is said to be sd (see Sato [48],
Cufaro Petroni [19]) when for every $0<a<1$ we can find another law with chf
$\chi_{a}(u)$ such that

$\eta(u)=\eta(au)\chi_{a}(u).$ (5)

Of course a rv$X$ with chf $\eta(u)$ is also said to be sd when its law is sd,
and looking at the definitions this means that for every $0<a<1$ we can always
find two _independent_ rv’s – a $Y$ with the same law of $X$, and a $Z_{a}$
with chf $\chi_{a}(u)$ – such that in distribution

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

Hereafter the rv $Z_{a}$ will be called the _$a$ -remainder_ of $X$ and in
general has an id (see Sato [48])

This last statement apparently means that the law of $Z(t)$ in the solution
(3) coincides with that of the $a$-remainder of the sd, stationary law
$\overline{\varphi}_{X}$ provided that $a=e^{-b\,t}$. It is easy indeed to see
from (3) that the chf of the time homogeneous transition law with a degenerate
initial condition $X(0)=x_{0},\;\bm{P}\hbox{-\emph{a.s.}}$ is

$\varphi_{X}(u,t|x_{0})=e^{\,ix_{0}ue^{-bt}}\varphi_{Z}(u,t)=\frac{\overline{\varphi}_{X}(u)\,e^{\,ix_{0}ue^{-bt}}}{\overline{\varphi}_{X}(u\,e^{-b\,t})}$
(7)

moreover we have

$\psi_{Z}(u,t)=\psi_{Z_{p}}(u)+\psi_{Z_{n}}(-u)$ (8)

and

$\psi_{X}(u,t|x_{0})=iux_{0}e^{-bt}+\psi_{Z}(u,t)=iux_{0}e^{-bt}+\psi_{Z_{p}}(u,t)+\psi_{Z_{n}}(-u,t).$
(9)

where $\psi_{Z_{d}}(u,t)=\ln\varphi_{Z_{d}}(u,t),d\in\\{p,n\\}$ is the lch of
$Z_{d}(t),d\in\\{p,n\\}$. Moreover, the transition lch of a OU process can
also be written in terms of the corresponding $\psi_{L}(u)$ in the form

$\psi_{X}(u,t|x_{0})=iux_{0}e^{-bt}+\psi_{Z}(u,t)=iux_{0}e^{-b\,t}+\int_{0}^{t}\psi_{L}\left(ue^{-b\,s}\right)ds.$
(10)

Finally in virtue of the results of Cufaro Petroni and Sabino [21], one can
relate the Lévy density $\nu_{Z}(x,t)$ of $Z(t)$ to that of the BDLP
$L(\cdot)$ at $t=1$ denoted with $\nu_{L}(x)$

$\displaystyle\nu_{Z}(x,t)$ $\displaystyle=$
$\displaystyle\frac{1}{b\,|x|}\left\\{\begin{array}[]{ll}\int_{\,{}^{x}\\!/_{a}}^{x}\nu_{L}(y)dy&\;\hbox{$x<0$}\\\
\\\
\int_{x}^{\,{}^{x}\\!/_{a}}\nu_{L}(y)dy&\;\hbox{$x>0$}\end{array}\right.\quad
a=e^{-bt}.$ (14)

##  3 OU-BCTS and OU-CGMY processes

In this section we study OU-BCTS and OU-CGMY processes with finite variation
and distinguish the case where the BDLP is of infinite activity,
$0<\alpha_{p}<1$ and $0<\alpha_{n}<1$, to that of finite activity namely when
$Z(\cdot)$ is a compound Poisson process. We do not discuss the setting
$\alpha_{p}=\alpha_{n}$ because it is already covered in Sabino [45] and
corresponds Variance Gamma driven OU process, therefore of infinite activity
and finite variation.

###  3.1 Infinite activity

Apparently, the study of the transition law of a OU-BCTS process $X(\cdot)$
coincides with the study of the process $Z(\cdot)$ and in particular of the
processes $Z_{p}(t)$ and $Z_{n}(t)$ defined in (4). One of course can rely on
these last two processes to build OU-CTS processes.

Cufaro Petroni and Sabino [21] and Qu et al. [43] designed an exact
decomposition of the transition law of OU-CTS processes as the convolution of
two independent rv’s plus a degenerate term. For simplicity, we report this
result here below in addition because, such a OU process is driven by a CTS
subordinator, we consider only $Z_{p}(\cdot)$.

######  Proposition 3.1.

For $0<\alpha_{p}<1$, and at every $t>0$, the pathwise solution of an OU-CTS
equation (3) with $X(0)=X_{0},\;\bm{P}\hbox{-\emph{a.s.}}$ is in distribution
the sum of three independent rv’s

$X(t)=aX_{0}+Z_{p}(t)\,\stackrel{{\scriptstyle
d}}{{=}}\,aX_{0}+X_{1}+X_{2}\qquad\quad a=e^{-b\,t}$ (15)

where $X_{1}$ is distributed according to the law
$\mathcal{CTS}\\!\left(\alpha_{p},\frac{\beta_{p}}{a},c_{p}\,\frac{1-a^{\alpha}_{p}}{\alpha_{p}\,b}\right)$,
whereas

$X_{2}=\sum_{k=1}^{N_{a}}J_{k}$

is a compound Poisson rv where $N_{a}$ is an independent Poisson rv with
parameter

$\Lambda_{a}=\frac{c_{p}\,\beta_{p}^{\alpha_{p}}\Gamma(1-\alpha_{p})}{b\,\alpha_{p}^{2}a^{\alpha_{p}}}\,\left(1-a^{\alpha_{p}}+a^{\alpha_{p}}\log
a^{\alpha_{p}}\right)$ (16)

and $J_{k},k>0$ are iid rv’s with density

$f_{J}(x)=\frac{\alpha_{p}\,a^{\alpha_{p}}}{1-a^{\alpha}_{p}+a^{\alpha_{p}}\log
a^{\alpha_{p}}}\int_{1}^{\frac{1}{a}}\frac{x^{-\alpha_{p}}\left(\beta_{p}\,v\right)^{1-\alpha_{p}}e^{-\beta_{p}v\,x}}{\Gamma(1-\alpha_{p})}\,\frac{v^{\alpha_{p}}-1}{v}\,dv$
(17)

namely, a mixture of a gamma law and a distribution with density

$f_{V}(v)=\frac{\alpha_{p}\,a^{\alpha}_{p}}{1-a^{\alpha_{p}}+a^{\alpha_{p}}\log
a^{\alpha_{p}}}\,\frac{v^{\alpha_{p}}-1}{v}\qquad\quad 1\leq
v\leq\,^{1}/_{a}.$ (18)

The extension to the bilateral OU-BCTS process is straightforward, for
instance the simulation algorithms consist of repeating the procedure for an
OU-CTS process two times.

The main contribution of this subsection is the derivation of the lch and
hence the chf and the moment generating function of $Z(t)$ that will be
instrumental to find the risk neutral conditions for market models based on
OU-BCTS processes and to the pricing of derivative contracts using FFT
methods.

######  Proposition 3.2.

The lch $\psi_{Z_{p}}(u,t)$, $u\in\mathbb{R}$ with $0<\alpha_{p}<1$ can be
represented as:

$\psi_{Z_{p}}(u,t)=-\frac{c_{p}\,\beta_{p}^{\alpha_{p}}\,\Gamma(1-\alpha_{p})}{\alpha_{p}\,b}\left[I\left(u,\alpha_{p},\beta_{p},\frac{\beta_{p}}{a}\right)+\log
a\right],\quad a=e^{-b\,t}$ (19)

with

$\displaystyle I\left(u,\alpha,\beta_{1},\beta_{2}\right)$ $\displaystyle=$
$\displaystyle\int_{\beta_{1}}^{\beta_{2}}z^{-1-\alpha}(z-iu)^{\alpha}dz=$
(20) $\displaystyle=$
$\displaystyle-\frac{1}{\alpha}\left[\left(\frac{u}{i\,\beta_{2}}\right)^{\alpha}\,_{2}F_{1}\left(-\alpha,-\alpha,1-\alpha,-\frac{i\,\beta_{2}}{u}\right)\right.-$
$\displaystyle=$
$\displaystyle\left.\left(\frac{u}{i\,\beta_{1}}\right)^{\alpha}\,_{2}F_{1}\left(-\alpha,-\alpha,1-\alpha,-\frac{i\,\beta_{1}}{u}\right)\right]$

where ${}_{2}F_{1}(a,b,c,x)$ is the hypergeometric function, $0<\alpha<1$,
$\beta_{1}>0$ and $\beta_{2}>0$. Finally, taking $\alpha=1/2$ the
hypergeometric function in (20) can be expressed in terms of elementary
functions as

${}_{2}F_{1}\left(-\frac{1}{2},-\frac{1}{2},\frac{1}{2},-x\right)=\sqrt{x+1}-\sqrt{x}\sinh^{-1}\left(\sqrt{x}\right)$

where $\sinh^{-1}x=\log\left(\sqrt{x^{2}+1}+x\right)$.

###### Proof.

From (14) it also results

$\nu_{Z_{p}}(u,t)=\frac{c_{p}}{b\,x}\int_{x}^{\frac{x}{a}}\frac{e^{-\beta_{p}y}}{y^{\alpha_{p}+1}}\,dy\;=\;\frac{c_{p}}{b}\int_{1}^{\frac{1}{a}}\frac{e^{-\beta_{p}wx}}{x^{\alpha_{p}+1}\,w^{\alpha_{p}+1}}\,dw$

therefore because of the Lévy-Khintchin theorem

$\psi_{Z_{p}}(u,t)=\frac{1}{b}\int_{1}^{\frac{1}{a}}w^{-\alpha_{p}-1}dw\int_{0}^{\infty}c_{p}\,\left(e^{i\,u\,x}-1\right)\frac{e^{-\beta_{p}wx}}{x^{\alpha_{p}+1}}\,dx$

The second integral is the lch of a
$\mathcal{CTS}(\alpha_{p},\beta_{p},c_{p})$ law with $0<\alpha_{p}<1$
therefore from Lemma 2.5 of Küchler and Tappe [33] we have

$\displaystyle\psi_{Z_{p}}(u,t)$ $\displaystyle=$
$\displaystyle\frac{c_{p}\,\Gamma(-\alpha_{p})}{b}\int_{1}^{\frac{1}{a}}\frac{\left(\beta_{p}\,w-i\,u\right)^{\alpha_{p}}-(\beta_{p}\,w)^{\alpha_{p}}}{w^{\alpha_{p}+1}}dw$
$\displaystyle=$
$\displaystyle\frac{c_{p}\,\Gamma(-\alpha_{p})}{b}\left[\int_{1}^{\frac{1}{a}}\frac{\left(\beta_{p}\,w-i\,u\right)^{\alpha_{p}}}{w^{\alpha_{p}+1}}dw-\beta_{p}^{\alpha_{p}}\int_{a}^{\frac{1}{a}}\frac{dw}{w}\right]$
$\displaystyle=$
$\displaystyle\frac{c_{p}\,\beta_{p}^{\alpha}\,\Gamma(-\alpha_{p})}{b}\left[\int_{\beta_{p}}^{\frac{\beta_{p}}{a}}z^{-1-\alpha_{p}}(z-iu)^{\alpha_{p}}dz+\log
a\right]$

where of course $-\alpha_{p}\Gamma(-\alpha_{p})=\Gamma(1-\alpha_{p})$ and in
last step we have used the change of variables $\beta_{p}\,w=z$. In order to
conclude the proof, we first observe that under the special case
$\gamma=\alpha+1$ the derivative of the hypergeometric function is

$\frac{d}{d\,x}\,_{2}F_{1}\left(\alpha,\beta,\alpha+1,x\right)=\frac{d}{d\,x}\,_{2}F_{1}(\beta,\alpha,\alpha+1,x)=-\frac{\alpha\left((1-x)^{\beta}-\,_{2}F_{1}(\alpha,\beta,\alpha+1,x\right))}{z}$

then with some algebra we get

$\frac{d}{d\,x}\left(z^{-\alpha}\,_{2}F_{1}\left(-\alpha,-\alpha,1-\alpha,-\frac{i\,z}{u}\right)\right)=\alpha\left(\frac{i}{u}\right)^{\alpha}z^{\alpha+1}\,\left(z-i\,u\right),$

therefore we can write the integral

$\displaystyle I\left(u,\alpha,\beta_{p},\frac{\beta_{p}}{a}\right)$
$\displaystyle=$
$\displaystyle-\frac{1}{\alpha}\left[\left(\frac{u}{i\,\beta_{2}}\right)^{\alpha}\,_{2}F_{1}\left(-\alpha,-\alpha,1-\alpha,-\frac{i\,\beta_{2}}{u}\right)\right.-$
$\displaystyle=$
$\displaystyle\left.\left(\frac{u}{i\,\beta_{1}}\right)^{\alpha}\,_{2}F_{1}\left(-\alpha,-\alpha,1-\alpha,-\frac{i\,\beta_{1}}{u}\right)\right]$

as claimed. ∎

######  Remark 1.

Several transformation and recursion formulas are applicable to the
hypergeometric functions. Using 9.131.1 in Gradshteyn and Rizhik [25] we have

${}_{2}F_{1}(-\alpha,-\alpha,1-\alpha,x)=(1-x)^{\alpha+1}\,_{2}F_{1}(1,1,1-\alpha,x)$

and accordingly (20) becomes

$\displaystyle I\left(u,\alpha,\beta_{1},\beta_{2}\right)$ $\displaystyle=$
$\displaystyle-\frac{i}{\alpha\,u}\left[\beta_{2}^{-\alpha}\left(\beta_{2}-i\,u\right)^{\alpha+1}\,_{2}F_{1}\left(1,1,1-\alpha,-\frac{i\,\beta_{2}}{u}\right)\right.-$
(21)
$\displaystyle\left.\beta_{1}^{-\alpha}\left(\beta_{1}-i\,u\right)^{\alpha+1}\,_{2}F_{1}\left(1,1,1-\alpha,-\frac{i\,\beta_{1}}{u}\right)\right]$

######  Corollary 3.3.

The cgf $m_{Z_{p}}(s,t)=\ln\bm{E}\left[{e^{s\,Z_{p}(t)}}\right]$ of $Z_{p}(t)$
with $0<\alpha_{p}<1$ exists for $s<\beta_{p}$ and is:

$m_{Z_{p}}(s,t)=-\frac{c\,\beta_{p}^{\alpha_{p}}\,\Gamma(1-\alpha_{p})}{\alpha_{p}\,b}\left[\tilde{I}\left(s,\alpha_{p},\beta_{p},\frac{\beta_{p}}{a}\right)+\log
a\right]\quad s<\beta_{p}$ (22)

where

$\displaystyle\tilde{I}\left(s,\alpha,\beta_{1},\beta_{2}\right)$
$\displaystyle=$
$\displaystyle\int_{\beta_{1}}^{\beta_{2}}z^{-1-\alpha}(z-s)^{\alpha}dz=$ (23)
$\displaystyle\frac{1}{\alpha\,s}\left[\beta_{2}^{-\alpha}\left(\beta_{2}-s\right)^{\alpha+1}\,_{2}F_{1}\left(1,1,1-\alpha,\frac{\beta_{2}}{s}\right)\right.-$
$\displaystyle\left.\beta_{1}^{-\alpha}\left(\beta_{1}-s\right)^{\alpha+1}\,_{2}F_{1}\left(1,1,1-\alpha,-\frac{\beta_{1}}{s}\right)\right]$

Note that setting $\psi_{Z_{p}}(-i\,s,t)=m_{Z_{p}}(s,t)$ in (20) one may claim
that the cgf assumes complex values which is obviously wrong and it explains
why we have preferred to rely on (21) to write (23).

In virtue of (8), Proposition 3.2 and Corollary 3.3 can be easily extended to
cope with $Z(t)$ defined for OU-BCTS processes.

######  Corollary 3.4.

For a OU-BCTS process, the lch $\psi_{Z}(u,t)$, $u\in\mathbb{R}$ can be
represented as:

$\displaystyle\psi_{Z}(u,t)$ $\displaystyle=$
$\displaystyle-\frac{c_{p}\,\beta_{p}^{\alpha_{p}}\,\Gamma(1-\alpha_{p})}{\alpha_{p}\,b}\left[I\left(u,\alpha_{p},\beta_{p},\frac{\beta_{p}}{a}\right)+\log
a\right]-$ (25)
$\displaystyle\frac{c\,\beta_{n}^{\alpha_{n}}\,\Gamma(1-\alpha_{n})}{\alpha_{n}\,b}\left[I\left(-u,\alpha_{n},\beta_{n},\frac{\beta_{n}}{a}\right)+\log
a\right]$

Accordingly,

######  Corollary 3.5.

The cgf at time $t$ exists for $-\beta_{n}<s<\beta_{p}$ and is:

$\displaystyle m_{Z}(s,t)$ $\displaystyle=$
$\displaystyle-\frac{c\,\beta_{p}^{\alpha_{p}}\,\Gamma(1-\alpha_{p})}{\alpha_{p}\,b}\left[\tilde{I}\left(s,\alpha_{p},\beta_{p},\frac{\beta_{p}}{a}\right)+\log
a\right]-$ (26)
$\displaystyle\frac{c\,\beta_{n}^{\alpha_{n}}\,\Gamma(1-\alpha_{n})}{\alpha_{n}\,b}\left[\tilde{I}\left(-s,\alpha_{n},\beta_{n},\frac{\beta_{n}}{a}\right)+\log
a\right]$

######  Remark 2.

In contrast to Proposition 3.1 that is valid under the condition
$0<\alpha_{p}<1$, we will show in the next subsection that Proposition 3.2 and
consequently all corollaries are also valid for $\alpha_{p}<0$ and
$\alpha_{n}<0$ as well.

###  3.2 Finite activity

When $\alpha_{p}<0$ the BDLP of a CTS process turns out to be a compound
Poisson process, indeed the integral of the Lévy density is convergent. In
more detail we have:

$\int_{0}^{\infty}\nu_{L_{p}}(x)dx=c_{p}\int_{0}^{\infty}x^{-\alpha_{p}-1}e^{-\beta_{p}\,x}dx=c_{p}\,\Gamma(-\alpha_{p})\,\beta_{p}^{-\alpha_{p}}=\lambda_{p}$
(27)

where now $-\alpha$ is positive. It results then that

$L_{p}(t)=\sum_{k=0}^{N_{p}(t)}J_{k},\quad J_{0}=0,\bm{P}\hbox{-\emph{a.s.}},$
(28)

where $N_{p}(t)$ is a Poisson process with the intensity $\lambda_{p}$ defined
in (27) and jumps sizes $J_{k}$ independent on $N_{p}(t)$ distributed
according to a gamma law with shape parameter $-\alpha_{p}$ and rate parameter
$\beta_{p}$ indeed the pdf of each copy of $J_{k}$ is

$f_{J}(x)=\frac{\nu_{L}(x)}{\lambda}=\frac{\beta_{p}^{-\alpha_{p}}\,x^{-\alpha_{p}-1\,e^{-\beta_{p}\,x}}}{\Gamma(-\alpha_{p})},\quad
x>0.$

This last representation is consistent with Jørgensen [30] that observed that
a compound Poisson process with gamma distributed jumps follows a Tweedie
distribution that is actually a CTS law.

######  Proposition 3.6.

For $\alpha_{p}<0$, and at every $t>0$, the pathwise solution of an OU-CTS
equation (3) with $X(0)=X_{0},\;\bm{P}\hbox{-\emph{a.s.}}$ is in distribution
the sum of two independent rv’s

$X(t)=aX_{0}+Z_{p}(t)\,\stackrel{{\scriptstyle
d}}{{=}}\,aX_{0}+X_{1}\qquad\quad a=e^{-b\,t}$ (29)

where $X_{1}$ can be written as $\sum_{k=0}^{N_{p}(t)}\tilde{J}_{k}$.
$N_{p}(t)$ is a Poisson process with intensity $\lambda_{p}$ given by equation
(27) and $\tilde{J}_{k},k>0$ are iid jumps distributed according to a mixture
of gamma law and a uniform distribution with pdf

$\int_{0}^{1}\frac{\left(\beta_{p}\,e^{b\,v\,t}\right)^{-\alpha_{p}}\,x^{-\alpha_{p}-1}e^{-\beta_{p}\,e^{b\,v\,t}\,x}}{\Gamma(-\alpha_{p})}dv$
(30)

or equivalently
$\tilde{J_{k}}\sim\Gamma(-\alpha_{p},\beta_{p}\,e^{b\,U_{k}t}),k>0$ where
$U_{k}\sim\mathcal{U}(0,1)$

###### Proof.

According to the definition of the OU-CTS process for $\alpha_{p}<0$ and the
representation (28) we can write

$X(t)=a\,X_{0}+\sum_{k=0}^{N_{p}(t)}J_{k}\,e^{-b(t-\tau_{k})}$

where $\tau_{k}$ are the jump times of the Poisson process $N_{p}(t)$ with
intensity $\lambda_{p}$. On the other hand, as observed by Lawrance [35] in
the context of Poisson point processes, for every $t>0$ we have

$\sum_{k=0}^{N_{p}(t)}J_{k}e^{-b(t-\tau_{k})}\;\stackrel{{\scriptstyle
d}}{{=}}\;\sum_{k=0}^{N_{p}(t)}J_{k}e^{-b\,t\,U_{k}},\quad
U_{0}=0,\bm{P}\hbox{-\emph{a.s.}}$

irrespective of the law of $J_{k}$, where $U_{k}\sim\mathcal{U}([0,1]),k>0$
form a sequence of iid uniformly distributed rv’s in $[0,1]$, also independent
on $J_{k}$. Knowing that for any gamma distributed random variable
$Y\sim\Gamma(\alpha,\beta)$, and $A>0$,
$A\,Y\sim\Gamma\left(\alpha,\frac{\beta}{A}\right)$, it results that
$\tilde{J}_{k}\sim\Gamma(-\alpha_{p},\beta_{p}\,e^{b\,U_{k}}),k>0$ that
concludes the proof. ∎

Because with $\alpha_{p}<0$, $Z_{p}(t)$ is distributed according to a compound
Poisson rv its lch is

$\displaystyle\psi_{Z_{p}}(u,t)$ $\displaystyle=$
$\displaystyle\lambda_{p}\,t\left(\int_{0}^{1}\left(\frac{\beta_{p}\,e^{b\,v\,t}}{\beta_{p}\,e^{b\,v\,t}-i\,u}\right)^{-\alpha_{p}}dv-1\right)=$
$\displaystyle=$
$\displaystyle\frac{\lambda_{p}}{b}\left(\int_{\beta_{p}}^{\beta_{p}\,e^{b\,t}}z^{-\alpha_{p}-1}(z-i\,u)^{\alpha_{p}}dz-b\,t\right)$

where in the last step we used the change of variables
$z=\beta_{p}\,e^{b\,v\,t}$. Replacing $a=e^{-b\,t}$ and $\lambda_{p}$ in (27)
we get the same representation of $\psi_{Z_{p}}(u,t)$ as that in Proposition
3.2; of course all corollaries of subsection 3.1 follow accordingly.

##  4 Simulation Algorithms

The sequential generation of the skeleton of an OU-BCTS or an OU-CGMY process
on a time grid $t_{0},t_{1},\dots,t_{I}$ simply consists in implementing the
following recursive procedure with initial condition $X(t_{0})=x_{0}$ taking
$a_{i}=e^{-b(t_{i}-t_{i-1})},\;i=1,\dots,I$:

$X(t_{i})=a_{i}X(t_{i-1})+Z_{a_{i}},\quad\qquad i=1,\dots,I.$ (31)

Cufaro Petroni and Sabino [21] have already discussed algorithms tailored for
OU-CTS processes of infinite activity, the extension to bilateral OU-BCTS or
OU-CGMY processes is straightforward.

In this section we illustrate the simulation procedure when these processes
are of finite activity that to our knowledge has not been investigated so far.
To this end, the simulation steps to generate the skeleton of a OU-CTS process
with parameters $b,\alpha_{p},\beta_{p},c_{p}$ is summarized in Algorithm 1

Algorithm 1

1$X_{0}\leftarrow x$

2for $i=1,\dots,I$ do

3 $\Delta t_{i}=t_{i}-t_{i-1}$, $a\leftarrow e^{-b\Delta t_{i}}$

4 $n\leftarrow N\sim\mathcal{P}(\lambda_{p}\,\Delta t_{i})$, $\triangleright$
Generate an independent Poisson rv with $\lambda$ in (27)

5 $u_{m}\leftarrow U_{m}\sim\mathcal{U}(0,1),m=1,\dots,n$ $\triangleright$
Generate $n$ iid uniform rv’s

6 $\tilde{\beta}_{m}\leftarrow\beta_{p}\,e^{b\,u_{m}},m=1,\dots,n$

7
$\tilde{j_{m}}\leftarrow\tilde{J_{m}}\sim\mathcal{G}(1-\alpha,\,\tilde{\beta}_{m}),m=1,\dots,n$
$\triangleright$ Generate $n$ independent gamma rv’s with scale $-\alpha_{p}$
and random rates $\tilde{\beta}_{m}$

8 $x_{1}\leftarrow\sum_{m=1}^{n}\tilde{j_{m}}$

9 $X(t_{i})\leftarrow a\,X(t_{i-1})+x_{1}$.

10end for

We remark that when $-\alpha_{p}=1$ the BDLP of the OU-CTS process is a
compound Poisson process with exponential jumps that corresponds to a OU
process with a gamma stationary law. For this configuration it is preferable
to use the faster and more efficient algorithm detailed in Sabino and Cufaro
Petroni [20, 47].

Finally, the procedure the generate the skeleton of OU-BCTS and OU-CGMY
processes simply entails to repeat steps $5$ to $8$ two times and add their
outcome to step $9$.

###  4.1 Numerical Experiments

In this section, we will assess the performance and the effectiveness of the
algorithms for the simulation of OU-BCTS process. All the simulation
experiments in the present paper have been conducted using _Python_ with a
$64$-bit Intel Core i5-6300U CPU, 8GB. The performance of the algorithms is
ranked in terms of the percentage error relatively to the first four cumulants
denoted _err_ % and defined as

$\text{err \%}=\frac{\text{true value}-\text{estimated value}}{\text{true
value}}$

Taking advantage of (10) one can calculate the cumulants
$c_{X,k}(x_{0},t),\;k=1,2,\ldots$ of $X(t)$ for $X_{0}=x_{0}$ from the
cumulants $c_{L,k}$ of the BCTS law according to

$\displaystyle c_{X,1}(x_{0},t)$ $\displaystyle=$
$\displaystyle\bm{E}\left[{X(t)|X_{0}=x_{0}}\right]\;=\;x_{0}e^{-b\,t}+\frac{c_{L,1}}{b}\left(1-e^{-b\,t}\right),\qquad
k=1$ (32) $\displaystyle c_{X,k}(x_{0},t)$ $\displaystyle=$
$\displaystyle\frac{c_{L,k}}{k\,b}\left(1-e^{-k\,b\,t}\right),\qquad\qquad\qquad\qquad\qquad\qquad
k=2,3,\ldots$ (33)

where

$c_{L,k}=\int_{-\infty}^{+\infty}x^{k}\nu_{L}(x)\,dx=c_{p}\,\beta_{p}^{\alpha_{p}-k}\Gamma(k-\alpha_{p})+(-1)^{k}c_{n}\,\beta_{n}^{\alpha_{n}-k}\Gamma(k-\alpha_{n})$
(34)

In our numerical experiments we consider a OU-CTS process with parameters
$\left(b,\beta_{p},c_{p}\right)=\left(0.5,1.5,0.3\right)$ and a OU-CGMY
process with $\left(b,C,G,M\right)=(0.5,0.3,0.5,1.5)$ with
$\alpha_{p}=Y\in\\{-0.5,-1.5,-2.5,-3.5\\}$.

The Tables 1 and 2 compare then the true values of the first four cumulants
$c_{X,k}(0,\Delta t)$ with their corresponding estimates from $10^{6}$ Monte
Carlo (MC) simulations respectively for the OU-CTS process with $\Delta
t=1/12$ and for the OU-CGMY process with $\Delta t=1/2$, each of the two with
the aforementioned parameters. We can conclude therefrom that the proposed
Algorithm 1 and its adaptation to the bilateral case produce unbiased
cumulants that are very close to their theoretical values. For the sake of
brevity, we do not report the additional results obtained with different
parameter settings that anyhow bring us to the same findings. Overall, from
the numerical results reported in this section, it is evident that the
Algorithm 1 proposed above can achieve a very high level of accuracy as well
as a conspicuous efficiency.

| $c_{X,1}(0,\Delta t)$ | $c_{X,2}(0,\Delta t)$ | $c_{X,3}(0,\Delta t)$ | $c_{X,4}(0,\Delta t)$
---|---|---|---|---
$\alpha_{p}$ | true | MC | err % | true | MC | err % | true | MC | err % | true | MC | err %
$-0.5$ | $1.181$ | $1.199$ | $-1.5$ | $1.157$ | $1.182$ | $-2.2$ | $1.889$ | $1.910$ | $-1.1$ | $4.320$ | $4.504$ | $-4.3$
$-1.5$ | $1.181$ | $1.205$ | $-2.0$ | $1.929$ | $2.020$ | $-4.7$ | $4.409$ | $4.616$ | $-4.7$ | $12.960$ | $13.422$ | $-3.6$
$-2.5$ | $1.969$ | $1.929$ | $2.0$ | $4.500$ | $4.378$ | $2.7$ | $13.226$ | $12.636$ | $4.5$ | $47.520$ | $45.830$ | $3.6$
$-3.5$ | $4.594$ | $4.539$ | $1.2$ | $13.500$ | $13.335$ | $1.2$ | $48.496$ | $48.099$ | $0.8$ | $205.921$ | $206.598$ | $-0.3$

Table 1: Comparing the first four true cumulants with their corresponding MC-
estimated values (multiplied by $100$) obtained with $10^{6}$ simulations and
$\Delta t=1/12$, $(b,\beta_{p},c_{p})=(0.5,1.5,0.3)$.

| $c_{X,1}(0,\Delta t)$ | $c_{X,2}(0,\Delta t)$ | $c_{X,3}(0,\Delta t)$ | $c_{X,4}(0,\Delta t)$
---|---|---|---|---
$Y$ | true | MC | err % | true | MC | err % | true | MC | err % | true | MC | err %
$-0.5$ | $-0.269$ | $-0.270$ | $-0.4$ | $0.945$ | $0.953$ | $-0.9$ | $-3.883$ | $-3.936$ | $-1.4$ | $25.13$ | $25.59$ | $-1.8$
$-1.5$ | $-0.934$ | $-0.935$ | $-0.1$ | $4.533$ | $4.526$ | $0.2$ | $-27.58$ | $-27.41$ | $0.6$ | $225.1$ | $221.1$ | $1.8$
$-2.5$ | $-4.883$ | $-4.890$ | $-0.1$ | $31.29$ | $31.35$ | $-0.2$ | $-249.4$ | $-249.3$ | $0.0$ | $2473$ | $2463$ | $0.4$
$-3.5$ | $-34.68$ | $-34.65$ | $0.1$ | $280.3$ | $280.2$ | $0.0$ | $-2747$ | $-2746$ | $0.1$ | $32127$ | $31636$ | $1.5$

Table 2: Comparing the first four true cumulants with their corresponding MC-
estimated values obtained with $10^{6}$ simulations and $\Delta t=1/2$,
$(b,C,G,M)=(0.5,0.3,0.5,1.5)$.

##  5 Financial Applications

In the following subsections we illustrate the application of the results
shown in Section 3 and of the simulation algorithm of Section 4 to the pricing
of derivative contracts in energy markets using models driven by OU-BCTS and
OU-CGMY processes. Energy markets and wider commodities markets exhibit mean-
reversion, seasonality and spikes, this last feature is particularly difficult
to capture with a pure Gaussian framework and motivates the use of Lévy
process. To this end, the literature is very rich of alternatives, for
instance Cartea and Figueroa [17] assumes that the evolution of the spot
prices follows a jump-diffusion OU process, whereas Meyer-Brandis and Tankov
[29] investigate the use of generalized OU processes.

Our model is similar to that of Benth et al. [7] and Benth and Benth [10]
where instead of NIG processes, we consider BCTS or CGMY processes as BDLP’s.
Our main goal, is to give the basis for the theoretical pricing and to provide
an efficient and exact simulation procedure rather to focus on the parameter
calibration and the model selection. Indeed, such dynamics based on OU-BCTS
and OU-CGYM processes can also find application in other financial contexts.

Our financial applications consider the pricing of a strip of call options
with a FFT-based approach, the evaluation of a forward start Asian option with
MC simulations and finally the pricing of a swing option using a modified
version of version of the Least-Squares Monte Carlo (LSMC), introduced in
Longstaff-Schwartz [36], detailed in Boogert and C. de Jong [12, 13].

We assume that the spot price is driven by the following one-factor process

$S(t)=F(0,t)\,e^{h(t)+X(t)}$ (35)

where $h(t)$ is a deterministic function, $F(0,t)$ is the forward curve
derived from quoted products and $X(t)$ is a OU-BCTS process. This market can
easily be turned into a multi-factor one, for instance adding a second CTS
process obtaining a tempered stable version of the two factor Gaussian model
of Schwartz and Smith [49]. We nevertheless focus on the model (35) to better
highlight the results obtained for the OU-BCTS and OU-CGMY processes.

Using Lemma 3.1 in Hambly et al. [26], the risk-neutral conditions are met
when the deterministic function $h(t)$ is consistent with forward curve such
that

$h(t)=-m_{X}(1,t)$ (36)

where $m_{X}(s,t)$ is the cgf $m_{X}(s,t)=s\,e^{-bt}+m_{Z}(s,t)$ and
$m_{Z}(s,t)$ is given by (26), therefore

$\displaystyle h(t)$ $\displaystyle=$
$\displaystyle\frac{c\,\beta_{p}^{\alpha_{p}}\,\Gamma(1-\alpha_{p})}{\alpha_{p}\,b}\left[\tilde{I}\left(1,\alpha_{p},\beta_{p},\frac{\beta_{p}}{a}\right)+\log
a\right]+$ (37)
$\displaystyle\frac{c\,\beta_{n}^{\alpha_{n}}\,\Gamma(1-\alpha_{n})}{\alpha_{n}\,b}\left[\tilde{I}\left(-1,\alpha_{n},\beta_{n},\frac{\beta_{n}}{a}\right)+\log
a\right]$

with $\beta_{p}>1$ and $\beta_{n}>0$. When $\alpha_{p}=\alpha_{n}=1/2$, the
integrals can be written in terms of the logarithmic function as follows

$\displaystyle\tilde{I}\left(1,\frac{1}{2},\beta_{1},\beta_{2}\right)$
$\displaystyle=$
$\displaystyle\int_{\beta_{1}}^{\beta_{2}}z^{-\frac{3}{2}}(z-1)^{\frac{1}{2}}dz=$
(38) $\displaystyle=$ $\displaystyle
2\,\beta_{2}^{\alpha}\,\left(\log\left(\sqrt{\beta_{2}}+\sqrt{\beta_{2}-1}\right)-\sqrt{\frac{\beta_{2}-1}{\beta_{2}}}\right)-$
$\displaystyle
2\,\beta_{1}^{\alpha}\,\left(\log\left(\sqrt{\beta_{1}}+\sqrt{\beta_{1}-1}\right)-\sqrt{\frac{\beta_{1}-1}{\beta_{1}}}\right),$
$\displaystyle\tilde{I}\left(-1,\frac{1}{2},\beta_{1},\beta_{2}\right)$
$\displaystyle=$
$\displaystyle\int_{\beta_{1}}^{\beta_{2}}z^{-\frac{3}{2}}(z+1)^{\frac{1}{2}}dz=$
(39) $\displaystyle=$ $\displaystyle
2\,\beta_{2}^{\alpha}\left(\log\left(\sqrt{\beta_{2}}+\sqrt{\beta_{2}+1}\right)-\sqrt{\frac{\beta_{2}+1}{\beta_{2}}}\right)-$
$\displaystyle
2\,\beta_{1}^{\alpha}\left(\log\left(\sqrt{\beta_{1}}+\sqrt{\beta_{1}+1}\right)-\sqrt{\frac{\beta_{1}+1}{\beta_{1}}}\right).$

###  5.1 Call Options

We consider a daily strip of $M$ call options with maturity $T$ and strike $K$
namely, a contract with payoff

$C(K,T)=\sum_{m=1}^{M}(S(t_{m})-K)^{+}=\sum_{m=1}^{M}c_{m}(K,t_{m}),\quad
t_{1},t_{2},\dots t_{M}=T.$

Such a contract is commonly used for hedging purposes or for the parameters
calibration. It normally encompasses monthly, quarterly and yearly maturities
but is not very liquid and is generally offered by brokers.

We assume that the market model (35) is driven by a full seven-parameters OU-
BCTS process with infinite activity and finite variation. We price the strip
of calls using the FFT-based technique of Carr and Madam [16] given the chf
$\phi(u,t)$ of the of $\log\,S(t)=\log F(0,t)+h(t)+X(t)$

$\phi(u,t)=F(0,t)e^{i\,u\,h(t)}\varphi_{X}(u,t)=F(0,t)e^{i\,u\,\left(h(t)+aX(0)\right)+\psi_{Z}(u,t)},\quad
a=e^{-b\,t}$

where $h(t)$ is given by (37) and $\psi_{Z}(u,t)$ by (25). We refer the reader
to Carr and Madam [16] for the details on the method.

The calibration and the parameters estimation is not the focus of this study,
instead we rather illustrate the applicability of our theoretical results
taking parameters sets available in the literature. In this example, we take
those of Poirot and Tankov [42] (plus $b$ and $c_{n}$) and let $\alpha_{p}$
and $\alpha_{n}$ vary:
$(b,\beta_{p},\beta_{n},c_{p},c_{n})=(0.1,2.5,3.5,0.5,1)$; for simplicity we
consider a flat forward curve with $F(0,t)=20,t>0$.

$\alpha_{n}$ $\alpha_{p}$ | 0.1 | 0.3 | 0.5 | 0.7 | 0.9
---|---|---|---|---|---
$0.1$ | $3.504$ | $3.540$ | $3.609$ | $4.262$ | $5.770$
$0.3$ | $4.865$ | $4.917$ | $5.008$ | $5.205$ | $6.290$
$0.5$ | $6.690$ | $6.757$ | $6.869$ | $7.073$ | $7.560$
$0.7$ | $9.058$ | $9.136$ | $9.261$ | $9.474$ | $9.879$
$0.9$ | $12.108$ | $12.192$ | $12.322$ | $12.535$ | $12.907$
Table 3: Strip of $M=30$ daily call options calculated with FFT, $T=1/12$
$(b,\beta_{p},c_{p})=(0.1,1.5,0.3)$.

Table 3 shows the values relatively to a strip of $M=30$ daily at-the-money
call options with maturity $T=1/12$ with different pairs of
$\alpha_{p},\alpha_{n}$. We observe that fixing one of $\alpha_{p}$ or
$\alpha_{n}$, the value of the option is increasing when the other one
increases. Moreover, Figure 1(a) illustrates the variability of the option
price with respect to the strike price $K$ where the dotting lines represent
the values obtained with $N=10^{5}$ MC simulations plus and minus three times
the estimation error (the root-mean squared error divided by $\sqrt{N}$). In
addition, Figure 1(b) compares the price of at-the-money options
$C_{m}=c_{m}(K,t_{m}),m=1,\dots,M$, $K=20$ obtained with the FFT method to
those estimated once again with $N=10^{5}$ MC simulations. In these last two
examples we have selected $\alpha_{p}=\alpha_{n}=0.5$.

As far as the MC method is concerned, the simulation of the skeleton of the
process is accomplished running the procedure explained in Cufaro Petroni and
Sabino [21] based on Proposition 4 two times because of the bilateral OU-BCTS;
the acceptance rejection step to draw from the law of $V$ in (18) assumes a
piece-wise approximation of the dominating functions into $L=100$ terms.

The results calculated with the FFT-method and with the MC method are totally
consistent on the other hand, it is well-know that the FFT approach is faster.
Nevertheless, a side-product of the MC approach are percentiles or other
statistics which are widely used by practitioners for risk-management
purposes.

(a) Effect on the strike

(b) $c_{m}(K,t_{m}),m=1,\dots,30$

Figure 1: Call-options values calculated with FFT and MC with $N=10^{5}$,
$(b,\beta_{p},\beta_{n},c_{p},c_{n})=(0.1,2.5,3.5,0.5,1)$,
$\alpha_{p}=\alpha_{n}=0.5$, $K=20$

###  5.2 Asian Options

As a second financial application we consider the pricing of Asian options
with MC simulations. In contrast to the previous example, we assume that the
market dynamic is driven by a OU-CGMY process with infinite activity and
finite variation with $C=c_{p}=c_{n}$, $G=\beta_{n}$, $M=\beta_{p}$ and
$Y=\alpha_{p}=\alpha_{n}$.

MC methods are known to be sometimes slower than FFT techniques that can also
be tailored to the pricing of Asian options (see Zhang and C. Oosterlee [51]).
Nevertheless, the former approach provides a view on the distribution of the
potential cash-flows of derivative contracts giving a precious information to
risk managers or to trading units.

We recall that the payoff at maturity $T$ of an Asian option with European
style and strike price $K$ is

$A(K,T)=\left(\frac{\sum_{i=1}^{I}}{I}S(t_{i})-K\right)^{+}.$

In this second example we consider once again a flat forward curve
$F(0,t)=20,t>0$ and a different parameter set. We select $b=10$ and
$\left(C,G,M\right)=\left(2,15,5\right)$ also used in Ballotta and Kyriakou
[2] and let $Y$ vary. Figure 2 displays a sample of four trajectories with
these parameters generated using the procedure of Cufaro Petroni and Sabino
[21] as done in the case of the daily strips of call options.

In addition to this simulation procedure we consider here two approximations:
the first boils down to simply neglect $X_{2}$ in the Proposition 3.1 and
accordingly to the bilateral case (_Approximation 1_); the second – in the
same vein of Benth et al. [8] dealing with the _normal inverse Gaussian_
-driven OU processes – takes advantage of the approximation of the law of
$Z(t)$ in (3) with that of $e^{-k\,t}L(t)$ (_Approximation 2_).

It turns our that under _Approximation 1_
$X_{1}(\Delta\,t)\sim\mathcal{BCTS}\\!\left(Y,Y,\frac{M}{a},\frac{G}{a},C\,\frac{1-a^{Y}}{Y\,b},C\,\frac{1-a^{Y}}{Y\,b}\right)$
whereas under _Approximation 2_
$X_{1}(\Delta\,t)\sim\mathcal{BCTS}\left(Y,Y,\frac{M}{a},\frac{G}{a},C\,\Delta\,t,C\,\Delta\,t\right)$,
$a=e^{-b\,\Delta\,t}$ where $\Delta\,t=1/360$ because we assume daily
settlements; for simplicity we adopt the convention that there are $360$ days
per year.

In order to highlight the differences between the estimations returned by the
exact method and those with the two alternatives, we consider two Asian
options both of them with $I=90$ daily settlements, the second option however
however, is a forward start contract whose first settlement date occurs after
$30$ days. The MC option values and their relative errors are reported in
Table 4 and Table 5 with different $Y$’s and number of simulations $N$.

Irrespective to the combination of $Y$ and $N$, for the option that start
settling after one day the exact solution and _Approximation 1_ return very
close values, whereas _Approximation 1_ is slightly biased. In contrast, for
the forward start contract, although the time steps for $m>1$ coincide and are
very small, for the simple fact that the first time step is relatively high,
the estimated prices returned by the two non-exact simulation schemes are very
biased and do not offer an acceptable alternative any longer. More important,
the bias cannot be controlled increasing the number of simulations as shown in
Table 5.

The cause of this difference comes from the fact that $X_{2}$ in Proposition
3.1 can be neglected when the time step is small (accordingly for the
bilateral case). Indeed taking the Taylor expansion of the parameter
$\Lambda_{a}$ in Proposition 3.1

$\Lambda_{a}=\frac{c\Gamma(1-\alpha)b\,\beta^{\alpha}}{2}\,\Delta\,t^{2}+o\big{(}\Delta\,t^{2}\big{)}.$

As mentioned, the parameters calibration is not the focus of this study,
nevertheless these observations could lead to a convenient strategy combining
parameters estimation and exact simulation of the OU-BCTS processes. Assuming
that the data could be made available with a fine enough time-granularity
(e.g. daily $\Delta\,t=1/360$), one could base the parameters estimation on
the likelihood methods by approximating the exact transition pdf of a OU-BCTS
process with that of a BCTS law. In alternative, one could also use the
generalized method of moments to historical data taking the cumulants from the
formulas (33) and (34). Instead, to avoid being forced to always simulate the
OU-BCTS processes on a fine time-grid allowing the approximations, the
generation of the skeleton of such processes will be preferably based on the
exact method.

Figure 2: Sample trajectories of OU-CGMY processes with
$\left(b,C,G,M\right)=\left(10,2,15,5\right)$ and $Y\in\\{0.3,0.5,0.7,0.9\\}$

Table 4: Asian option. $K=20$, $T=1/4$ | Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
---|---|---|---|---|---|---
$N$ | price | error | price | error | price | error | price | error | price | error | price | error
| $Y=0.3$ | $Y=0.5$
$1000$ | $0.3293$ | $0.0224$ | $0.3771$ | $0.0275$ | $0.3671$ | $0.0253$ | $0.4439$ | $0.0278$ | $0.4461$ | $0.0300$ | $0.4865$ | $0.0326$
$10000$ | $0.3894$ | $0.0087$ | $0.3709$ | $0.0085$ | $0.3655$ | $0.0085$ | $0.4699$ | $0.0098$ | $0.4673$ | $0.0096$ | $0.4499$ | $0.0096$
$20000$ | $0.3692$ | $0.0060$ | $0.3743$ | $0.0062$ | $0.3651$ | $0.0060$ | $0.4571$ | $0.0069$ | $0.4404$ | $0.0065$ | $0.4552$ | $0.0068$
$50000$ | $0.3808$ | $0.0040$ | $0.3814$ | $0.0039$ | $0.3718$ | $0.0038$ | $0.4634$ | $0.0043$ | $0.4594$ | $0.0043$ | $0.4545$ | $0.0043$
$100000$ | $0.3792$ | $0.0028$ | $0.3715$ | $0.0027$ | $0.3697$ | $0.0027$ | $0.4627$ | $0.0031$ | $0.4559$ | $0.0031$ | $0.4467$ | $0.0030$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | price | error | price | error | price | error | price | error | price | error | price | error
| $Y=0.7$ | $Y=0.9$
$1000$ | $0.5446$ | $0.0325$ | $0.5317$ | $0.0332$ | $0.5298$ | $0.0325$ | $0.6910$ | $0.0427$ | $0.7031$ | $0.0377$ | $0.7246$ | $0.0395$
$10000$ | $0.5902$ | $0.0114$ | $0.5630$ | $0.0110$ | $0.5545$ | $0.0106$ | $0.7102$ | $0.0130$ | $0.7227$ | $0.0134$ | $0.6912$ | $0.0129$
$20000$ | $0.5851$ | $0.0080$ | $0.5631$ | $0.0075$ | $0.5544$ | $0.0077$ | $0.7147$ | $0.0091$ | $0.7066$ | $0.0091$ | $0.6973$ | $0.0090$
$50000$ | $0.5647$ | $0.0049$ | $0.5631$ | $0.0048$ | $0.5413$ | $0.0047$ | $0.7209$ | $0.0058$ | $0.6933$ | $0.0057$ | $0.6920$ | $0.0057$
$100000$ | $0.5701$ | $0.0035$ | $0.5603$ | $0.0034$ | $0.5491$ | $0.0034$ | $0.7148$ | $0.0041$ | $0.7117$ | $0.0041$ | $0.6897$ | $0.0040$
Table 5: Forward start Asian option. $K=20$, $T=1/3$ | Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
---|---|---|---|---|---|---
$N$ | price | error | price | error | price | error | price | error | price | error | price | error
| $Y=0.3$ | $Y=0.5$
$1000$ | $0.4361$ | $0.0314$ | $0.3793$ | $0.0274$ | $0.3536$ | $0.0261$ | $0.4734$ | $0.0300$ | $0.4621$ | $0.0301$ | $0.4226$ | $0.0312$
$10000$ | $0.4400$ | $0.0098$ | $0.3751$ | $0.0086$ | $0.3689$ | $0.0089$ | $0.5351$ | $0.0107$ | $0.4527$ | $0.0094$ | $0.4161$ | $0.0091$
$20000$ | $0.4522$ | $0.0071$ | $0.3884$ | $0.0062$ | $0.3775$ | $0.0063$ | $0.5296$ | $0.0075$ | $0.4602$ | $0.0067$ | $0.4294$ | $0.0066$
$50000$ | $0.4458$ | $0.0043$ | $0.3731$ | $0.0038$ | $0.3694$ | $0.0038$ | $0.5303$ | $0.0048$ | $0.4551$ | $0.0043$ | $0.4340$ | $0.0042$
$100000$ | $0.4480$ | $0.0031$ | $0.3766$ | $0.0027$ | $0.3685$ | $0.0027$ | $0.5312$ | $0.0034$ | $0.4632$ | $0.0031$ | $0.4274$ | $0.0030$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | price | error | price | error | price | error | price | error | price | error | price | error
| $Y=0.7$ | $Y=0.9$
$1000$ | $0.5575$ | $0.0347$ | $0.6160$ | $0.0376$ | $0.5185$ | $0.0329$ | $0.8454$ | $0.0477$ | $0.7726$ | $0.0429$ | $0.6411$ | $0.0399$
$10000$ | $0.6397$ | $0.0120$ | $0.5672$ | $0.0110$ | $0.5268$ | $0.0107$ | $0.7941$ | $0.0141$ | $0.7403$ | $0.0132$ | $0.6129$ | $0.0121$
$20000$ | $0.6517$ | $0.0086$ | $0.5789$ | $0.0079$ | $0.4967$ | $0.0072$ | $0.8092$ | $0.0101$ | $0.7450$ | $0.0093$ | $0.6013$ | $0.0086$
$50000$ | $0.6574$ | $0.0055$ | $0.5721$ | $0.0050$ | $0.5129$ | $0.0048$ | $0.8124$ | $0.0064$ | $0.7429$ | $0.0060$ | $0.5976$ | $0.0054$
$100000$ | $0.6500$ | $0.0039$ | $0.5719$ | $0.0035$ | $0.5094$ | $0.0033$ | $0.8142$ | $0.0046$ | $0.7399$ | $0.0042$ | $0.5970$ | $0.0038$

###  5.3 Swing Options

A swing option is a type of contract used by investors in energy markets that
lets the option holder buy a predetermined quantity of energy at a
predetermined price (strike), while retaining a certain degree of flexibility
in both the amount purchased and the price paid.

Let the maturity date $T$ be fixed and the payoff at time $t<T$ be given by
$(S(t)-K)^{+}$ where $K$ denotes the strike price, in addition we assume only
one unit of the underlying can be exercised any time period. Let $V(n,s,t)$
denote the price of such a swing option at time $t$ given the spot price $s$
which has $n$ out of $N$ exercise rights left. For $m=1,\dots,M_{S}$, the
dynamic programming principle allows us to write (see Bertsekas [11])

$V(n,s,t_{m})=\max\left\\{\begin{array}[]{ll}\bm{E}\left[{V(n,S(t_{m+1}),t_{m+1})|S(t_{m})=s}\right],\\\
\\\
\bm{E}\left[{V(n,S(t_{m+1}),t_{m+1})|S(t_{m})=s}\right]+(s-K)^{+}\end{array}\right\\},\quad
n<N$ (40)

and $V(n,s,T)=(S(T)-K)^{+}$, $n\leq N$ and $V(0,s,t)=0$. In order to solve the
recursion equation we rely on the modified version of the LSMC, introduced in
Longstaff-Schwartz [36], detailed in Boogert and C. de Jong [12, 13] where the
continuation value is approximated with a linear regression with
$m=1,\dots,M_{S}$

$\bm{E}\left[{V(n,S(t_{m+1}),t_{m+1})|S(t_{m})=s}\right]\simeq
a_{0}+a_{1}S(t_{m})+\dots,+a_{B}S^{B}(t_{m}),\quad n<N.$

In our experiments, we used simple power polynomials with $B=3$, but the
regression may be performed on a different set of basis functions as well (see
Boogert and de Jong [13] for a comparison with other basis functions).

Several other approaches have been proposed: for instance one may solve the
recursion by adapting the method of Ben-Ameur et al. [6] or might use the
quantization technique of Bardou et al. [3]. In alternative, one can also use
the tree method of Jaillet et al. [28] or the Fourier cosine expansion in
Zhang and C. Oosterlee [50] taking advantage of the explicit form of the chf
of OU-BCTS process.

In this last example we assume a OU-CGMY driven market model with $Y<0$,
namely a combination of mean-reverting compound Poisson processes with
positive and negative jumps. We consider a different set of parameters
compared to the cases illustrated so far, namely we take
$(b,C,G,M)=(25,80,10.5,15.5)$ and let $Y$ vary. The parameters are very
different than the other two examples and are chosen to mimic realistic price
path as shown in Figure 3. We also remark that, due to the fact that energy
markets are very seasonal and spikes occur in clusters due to for instance,
cold spells, one could assume that the intensity of the compound Poisson
processes is a seasonal time-dependent function. The results in Subsection 3.2
and the simulation algorithms in Section 4 can be easily adapted taking a
step-wise approximation of the intensity function.

Table 6 shows the values and MC errors relatively to the pricing of a
$120-120$ swing option with maturity $T=1$ and strike price $K=20$, namely the
holder has $N=120$ rights and must exercise all of them. We observe that the
LSMC combined with Algorithm 1 produces unbiased results for all selected
$Y$’s and apparently $2\times 10^{4}$ simulations are required to attain an
acceptable convergence. In contrast to the Asian option case, it does not make
sense to adopt the approximation of the law of $Z(t)$ in (3) with that of
$e^{-k\,t}L(t)$ (_Approximation 2_ in Subsection 5.2) because this approach
returns another compound Poisson process and therefore does not provide a
computational advantage. Overall, it is evident that our newly developed
approach can achieve high accuracy as well as efficiency.

Figure 3: Sample trajectories of OU-CGMY processes with
$\left(b,C,G,M\right)=\left(10,10,1.75,1.25\right)$ and
$Y\in\\{-0.3,-0.5,-0.7,-0.9\\}$

Table 6: 120-120 Swing option. $K=20$, $T=1$ | $Y=0.3$ | $Y=0.5$ | $Y=0.7$ | $Y=0.9$
---|---|---|---|---
$N$ | price | error | price | error | price | error | price | error
$1000$ | $101.965$ | $2.873$ | $80.500$ | $2.446$ | $66.634$ | $2.222$ | $51.117$ | $1.983$
$10000$ | $99.565$ | $0.902$ | $79.848$ | $0.778$ | $63.451$ | $0.682$ | $49.801$ | $0.599$
$20000$ | $98.328$ | $0.633$ | $79.774$ | $0.550$ | $64.322$ | $0.489$ | $49.958$ | $0.426$
$50000$ | $98.187$ | $0.399$ | $79.585$ | $0.348$ | $63.686$ | $0.307$ | $50.223$ | $0.269$
$100000$ | $98.270$ | $0.283$ | $79.284$ | $0.245$ | $63.949$ | $0.217$ | $50.487$ | $0.192$

###  5.4 Application to Forward Markets

So far, we have discussed the pricing of derivative contracts depending on the
day-ahead price which is modeled as a OU process. On the other hand, the main
point of Section 3 is the study of the additive process
$Z(t)=\int_{0}^{t}e^{-b\,(t-u}\,dL(u)$ irrespective of the construction of a
OU-BCTS or a OU-CGMY process. For instance, it is common practice to model the
forward price as a geometric Brownian motion with a time-dependent volatility
function that captures the Samuelson effect. For instance, Kiesel et al. [31]
have considered a two-factor market dynamics where one of the two factors
depends on $\int_{0}^{t}e^{-b\,(t-u)}\,dW(u)$ where $W(\cdot)$ is a standard
Wiener process.

Beyond the Gaussian world, Piccirilli et al. [41] have recently proposed a
class of models, named Non-Overlapping-Arbitrage models (NOA), with the aim or
capturing the Samuelson effect and reproducing the different levels and shapes
of the implied volatility profiles displayed by options.

In particular, they assume a stochastic evolution of a generic future price at
time $t$, maturity $T$, $t\leq T\leq T_{1}<T_{2}$, and with delivery period
$[T_{1},T_{2}]$ is described by

$\displaystyle F(t,T_{1},T_{2})$ $\displaystyle=$ $\displaystyle
F(0,T_{1},T_{2})+\int_{0}^{t}\Gamma_{1}(u,T_{1},T_{2})dL_{1}(u)+\Gamma_{2}(T_{1},T_{2})\,X_{2}(t)$
(41) $\displaystyle=$ $\displaystyle
F(0,T_{1},T_{2})+X_{1}(t,T_{1},T_{2})+\Gamma_{2}(T_{1},T_{2})\,X_{2}(t)$

where $L_{1}(\cdot)$ and $X_{2}(\cdot)$ are two independent Lévy processes.
Moreover,

$\Gamma_{1}(u,T_{1},T_{2})=\frac{\gamma_{1}}{b\,(T_{2}-T_{1})}\left(e^{-b\,(T_{1}-u)}-e^{-b\,(T_{2}-u)}\right).$
(42)
$\Gamma(T_{1},T_{2})=\frac{1}{T_{2}-T_{1}}\int_{T_{1}}^{T_{2}}\gamma(u)\,du$
(43)

are two deterministic functions that are meant to capture the Samuelson effect
in option pricing (see also Jaeck and Lautier [27]). Indeed, in the spirit of
Benth et al. [9] and Latini et al. [34], the special form of the coefficients
arises from the implicitly underlying assumption that the future can be
written as the average over an underlying artificial futures price with
instantaneous delivery.

Although Piccirilli et al. [41] illustrate the application of their model
under the assumption that $L_{1}(\cdot)$ and $X_{2}(\cdot)$ are centered NIG
processes, the setting can be modified taking two independent BCTS or CGYM
processes. Of course, such models are related to the additive process studied
in Section 3, because, after some algebra it results

$\displaystyle X_{1}(t,T_{1},T_{2})$ $\displaystyle=$
$\displaystyle\frac{\gamma_{1}}{b\,(T_{2}-T_{1})}\left(e^{-b\,(T_{1}-T)}-e^{-b\,(T_{2}-T)}\right)\int_{0}^{t}e^{-b\,(t-u)}dL_{1}(u)$
$\displaystyle=$ $\displaystyle\Gamma_{1}(T,T_{1},T_{2})Z(t),$

hence the chf and in the particular, the simulation procedure of the skeleton
of the additive process $X_{1}(\cdot,T_{1},T_{2})$ can be derived from those
of $Z(\cdot)$. It is worthwhile noticing that the Piccirilli et al. [41] found
an explicit form of the chf of $X_{1}(t,T_{1},T_{2})$ when $L_{1}(\cdot)$ is a
centered NIG process, whereas they do not provide any procedure to simulate
such a process. In alternative, taking BCTS or CGMY processes and slightly
modifying our results, one can get both the chf and the simulation procedure
giving the possibility to price other derivative contracts via Monte Carlo
simulations. We omit an explicit proof to avoid overloading the paper with
lengthy details.

##  6 Concluding Remarks

In this study we have investigated the pricing of energy derivatives in
markets driven by classical tempered stable and CGMY processes of OU type with
finite variation. To this end, we have derived the chf of the transition law
of such processes in closed form such that we can obtain the non-arbitrage
conditions and spot prices that are consistent with the forward curve. In
addition, extending the work of Cufaro Petroni and Sabino [21], we have
detailed efficient algorithms for the simulation of the skeleton of classical
tempered stable and CGMY processes of OU type with particular focus to the
case when whey coincide with compound Poisson processes. We have illustrated
the applicability of these results to the pricing of three common derivative
contracts in energy markets, namely a strip of daily call options, an Asian
option with European style and a swing option. In our numerical experiments we
have selected a one-factor model in order to better highlight the features of
our finding, nevertheless the extension to two-factor models in the same vein
of Schwartz and Smith [49] is straightforward. In the first example we have
made use of the explicit knowledge of the chf to implement the pricing with
the FFT-based technique of Carr and Madan [16] and have compared the outcomes
to those obtained via MC simulations. In the second example, we have priced
Asian options with MC simulations where we have also adopted two common
approximations techniques. These approximations provide reliable values if the
time steps of the time grid are relatively small but if one considers a
forward start contract the outcome is really biased. Although the parameter
calibration and the model selection is not the main focus of this study, these
observations give an indication of how one could conceive a simplified
procedure for the parameters estimation. In addition, we have shown that the
proposed simulation algorithm, combined with the LSMC approach of Boogert and
C. de Jong [12, 13], provides an efficient and accurate pricing of a one year
$120-120$ swing option. Furthermore, our results are not restricted to OU
processes and to the modeling of spot prices. Indeed, in the spirit of Benth
et al. [9], Latini et al. [34] and Piccirilli et al. [41] they can be adapted
to capture the Samuelson effect and to volatility smiles.

Finally, future studies could cover the extension to a multidimensional
framework for instance adopting the view of Luciano and Semeraro [38],
Ballotta and Bonfiglioli [1] or the recent approaches of Gardini et al. [22,
23] and Lu [37]. A last topic deserving further investigation is the time-
reversal simulation of the OU processes generalizing the results of Pellegrino
and Sabino [40] and Sabino [46] to the case of classical tempered stable and
CGMY processes.

## Acknowledgements

I would like to express my gratitude to Matteo Gardini and Nicola Cufaro
Petroni for their help relatively to the application of the FFT method.

## References

  * [1] L. Ballotta and E. Bonfiglioli.  Multivariate Asset Models Using Lévy Processes and Applications.  The European Journal of Finance, 13(22):1320–1350, 2013.
  * [2] L. Ballotta and I. Kyriakou.  Monte Carlo Simulation of the CGMY Process and Option Pricing.  Journal of Futures Markets, 34(12):1095–1121, 2014.
  * [3] O. Bardou, S. Bouthemy, and G. Pagés.  Optimal Quantization for the Pricing of Swing Options.  Applied Mathematical Finance, 16(2):183–217, 2009.
  * [4] O. E. Barndorff-Nielsen, J. L. Jensen, and M. Sørensen.  Some Stationary Processes in Discrete and Continuous Time.  Advances in Applied Probability, 30(4):989–1007, 1998.
  * [5] O.E. Barndorff-Nielsen and N. Shephard.  Non-Gaussian Ornstein-Uhlenbeck-based Models and some of their Uses in Financial Economics.  Journal of the Royal Statistical Society: Series B, 63(2):167–241, 2001.
  * [6] H. Ben-Ameur, M. Breton, L. Karoui, and P. L’Ecuyer.  A Dynamic Programming Approach for Pricing Options Embedded in Bonds.  Journal of Economic Dynamics and Control, 31(7):2212–2233, July 2007.
  * [7] F.E. Benth, J. Kallsen, and T. Meyer-Brandis.  A non-gaussian ornstein-uhlenbeck process for electricity spot price modeling and derivatives pricing.  Applied Mathematical Finance, 14(2):153–169, 2007.
  * [8] F.E. Benth, L. Di Persio, and S. Lavagnini.  Stochastic Modeling of Wind Derivatives in Energy Markets.  Risks, MDPI, Open Access Journal, 6(2):1–21, 2018.
  * [9] F.E. Benth, M. Piccirilli, and T. Vargiolu.  Mean-reverting Additive Energy Forward Curves in a Heath–Jarrow–Morton Framework.  Mathematics and Financial Economics, 13:543–577, 2019.
  * [10] F.E. Benth and J. Šaltyté Benth.  The Normal Inverse Gaussian Distribution and Spot Price Modelling in Energy Markets.  International Journal of Theoretical and Applied Finance, 07(02):177–192, 2004.
  * [11] D. P. Bertsekas.  Dynamic Programming and Optimal Control, Volume I.  Athena Scientific, Belmont, Mass., third edition, 2005.
  * [12] A. Boogert and C. de Jong.  Gas Storage Valuation Using a Monte Carlo Method.  Journal of Derivatives, 15:81–91, 2008.
  * [13] A. Boogert and C. de Jong.  Gas Storage Valuation using a Multifactor Price Model.  The Journal of Energy Markets, 4:29–52, 2011.
  * [14] P. Carr and J. Crosby.  A Class of Lévy Process Models with almost Exact Calibration to both Barrier and Vanilla FX Options.  Quantitative Finance, 10(10):1115–1136, 2010.
  * [15] P. Carr, H. Geman, D.B. Madan, and M. Yor.  The fine structure of asset returns: An empirical investigation.  The Journal of Business, 75(2):305–332, 2002.
  * [16] P. Carr and D.B. Madan.  Option Valuation Using the Fast Fourier Transform.  Journal of Computational Finance, 2:61–73, 1999.
  * [17] A. Cartea and M. Figueroa.  Pricing in Electricity Markets: a Mean Reverting Jump Diffusion Model with Seasonality.  Applied Mathematical Finance, No. 4, December 2005, 12(4):313–335, 2005.
  * [18] R. Cont and P. Tankov.  Financial Modelling with Jump Processes.  Chapman and Hall, London, 2004.
  * [19] N. Cufaro Petroni.  Self-decomposability and Self-similarity: a Concise Primer.  Physica A, Statistical Mechanics and its Applications, 387(7-9):1875–1894, 2008.
  * [20] N. Cufaro Petroni and P. Sabino.  Fast Pricing of Energy Derivatives with Mean-reverting Jump-diffusion Processes.  Available at: https://arxiv.org/abs/1908.03137.
  * [21] N. Cufaro Petroni and P. Sabino.  Tempered Stable Distribution and Finite Variation Ornstein-Uhlenbeck Processes.  Available at: https://arxiv.org/abs/2011.09147.
  * [22] M. Gardini, P. Sabino, and E. Sasso.  A Bivariate Normal Inverse Gaussian Process with Stochastic Delay: Efficient Simulations and Applications to Energy Markets, 2020.  Available at www.arxiv.org.
  * [23] M. Gardini, P. Sabino, and E. Sasso.  Correlating Lévy Processes with Self-decomposability: Applications to Energy Markets, 2020.  Available at www.arxiv.org.
  * [24] M. Grabchak.  Tempered Stable Distributions.  Springer International Publishing, 2016.
  * [25] I. S. Gradshteyn and I. M. Ryzhik.  Table of Integrals, Series, and Products.  Elsevier/Academic Press, Amsterdam, seventh edition, 2007.
  * [26] B. Hambly, S. Howison, and T. Kluge.  Information-Based Models for Finance and Insurance.  Quantitative Finance, 9(8):937–949, 2009.
  * [27] E. Jaeck and D. Lautier.  Volatility in Electricity Derivative Markets: The Samuelson Effect Revisited.  Energy Economics, 59:300–313, 2016.
  * [28] P. Jaillet, E.I. Ronn, and S. Tompaidis.  Valuation of Commodity-Based Swing Options.  Management Science, 50(7):909–921, 2004.
  * [29] J.Kallsen and P. Tankov.  Characterization of Dependence of Multidimensional Lévy Processes Using Lévy Copulas.  Journal of Multivariate Analysis, 97(7):1551–1572, 2006.
  * [30] B. Jørgensen.  The Theory of Dispersion Models.  Chapman & Hall, 1997.
  * [31] R. Kiesel, G. Schindlmayr, and R.H. Börger.  A Two-factor Model for the Electricity Forward Market.  Quantitative Finance, 9(3):279–287, 2009.
  * [32] I. Koponen.  Analytic Approach to the Problem of Convergence of Truncated Lévy Flights Towards the Gaussian Stochastic Process.  Phys. Rev. E, 52:1197–1199, Jul 1995.
  * [33] U. Küchler and S. Tappe.  Tempered Stable Distribution and Processes.  Stochastic Processes and their Applications, 123(12):4256 – 4293, 2013.
  * [34] L. Latini, M. Piccirilli, and T. Vargiolu.  Mean-reverting No-arbitrage Additive Models for Forward Curves in Energy Markets.  Energy Economics, 79:157–170, 2019.  Energy Markets Dynamics in a Changing Environment.
  * [35] A.J Lawrance.  Some Autoregressive Models for Point Processes.  In P. Bartfai and J. Tomko, editors, Point Proceses and Queueing Problems (Colloquia Mathematica Societatis János Bolyai 24), volume 24, pages 257–275. North Holland, Amsterdam, 1980.
  * [36] F. A. Longstaff and E.S. Schwartz.  Valuing American Options by Simulation: a Simple Least-Squares Approach.  Review of Financial Studies, 14(1):113–147, 2001.
  * [37] K.W. Lu.  Calibration for Multivariate Lévy-Driven Ornstein-Uhlenbeck Processes with Applications to Weak Subordination, 2020.  Available at www.arxiv.org.
  * [38] 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.
  * [39] D. B. Madan and E. Seneta.  The Variance Gamma (V.G.) Model for Share Market Returns.  The Journal of Business, 63(4):511–24, 1990.
  * [40] T. Pellegrino and P. Sabino.  Enhancing Least Squares Monte Carlo with Diffusion Bridges: an Application to Energy Facilities.  Quantitative Finance, 15(5):761–772, 2015.
  * [41] M. Piccirilli, M.D. Schmeck, and T. Vargiolu.  Capturing the Power Options Smile by an Additive Two-factor Model for Overlapping Futures Prices.  Energy Economics, 95:105006, 2021.
  * [42] J. Poirot and Peter P. Tankov.  Monte Carlo Option Pricing for Tempered Stable (CGMY) Processes.  Asia-Pacific Financial Markets, 13(4):327–344, 2006.
  * [43] Y. Qu, A. Dassios, and H. Zhao.  Exact Simulation of Ornstein–Uhlenbeck Tempered Stable Processes.  Journal of Applied Probability, 0(0), 2021.  Forthcoming.
  * [44] Jan Rosinski.  Tempering Stable Proceses.  Stochastic Processes and their Applications, 117(6):677 – 707, 2007\.
  * [45] P. Sabino.  Exact Simulation of Variance Gamma Related OU Proceses: Application to the Pricing of Energy Derivatives.  Applied Mathematical Finance, 27(3):207–227, 2020.
  * [46] P. Sabino.  Forward or Backward Simulations? A Comparative Study.  Quantitative Finance, 20(7):1213–1226, 2020.
  * [47] 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\.
  * [48] K. Sato.  Lévy Processes and Infinitely Divisible Distributions.  Cambridge U.P., Cambridge, 1999.
  * [49] P. Schwartz and J.E. Smith.  Short-term Variations and Long-term Dynamics in Commodity Prices.  Management Science, 46(7):893–911, 2000.
  * [50] B. Zhang and C.W. Oosterlee.  An Efficient Pricing Algorithm for Swing Options based on Fourier Cosine Expansions.  Journal of Computational Finance, 16(4):1–32, 2013.
  * [51] B. Zhang and C.W. Oosterlee.  Efficient Pricing of European-style Asian Options under Exponential Lévy Processes based on Fourier Cosine Expansions.  SIAM J. Financial Math., 4(1):399–426, 2013.