# Normal Tempered Stable Processes and the Pricing of Energy Derivatives

Piergiacomo Sabino<EMAIL_ADDRESS>
The 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 is modeled with a normal tempered stable driven Ornstein-
Uhlenbeck process. Such processes are the generalization of normal inverse
Gaussian processes that are widely used in energy finance applications. We
first specify their statistical properties calculating their characteristic
function 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 without relying on numerical approximations or on numerical
integration. Moreover, we conceive an efficient algorithm for the exact
generation of the trajectories which gives the possibility to implement Monte
Carlo simulations without approximations or bias. We illustrate the
applicability of the theoretical findings and the simulation algorithms in the
context of the pricing of 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; Normal Tempered Stable
processes; Simulations; Energy Markets; Derivative Pricing

##  1 Introduction

Energy and commodity markets exhibit mean-reversion, seasonality and sudden
spikes; mean-reversion in particular, cannot be captured by ordinary Lévy
processes. To this end, the modeling based on non-Gaussian Ornstein-Uhlenbeck
(OU) processes has received considerable attention in the recent literature in
an attempt to accommodate features such as jumps, heavy tails which are well
evident in real data.

A common approach to describe the evolution of day-ahead (spot) prices assumes
the following multi-factor dynamics

$S(t)=F(0,t)\,\exp\left(h(t)+X(t)\right)=F(0,t)\,\exp\left(h(t)+\sum_{n}^{N}X_{n}(t)+\sum_{m}^{L}Y_{m}(t)\right).$

$F(0,t)$ represents the forward curve derived from quoted products and
reflects the seasonality, $X_{n}(\cdot),n=1,\dots$ are independent mean-
reverting Lévy-driven OU processes and $Y_{m}(\cdot),m=1,\dots$ are
independent plain Lévy processes. Due to the Lemma 3.1 in Hambly et al. [24],
the risk-neutral conditions are met when the deterministic function $h(t)$ is
consistent with the forward curve such that

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

where $m_{X}(s,t)$ is the logarithm of the moment generating function of
$X(t)$. Standard examples in the Gaussian framework are the one factor model
($N=1$, $M=0$) of Lucia and Schwarz [34] and the two factor version ($M=N=1$)
of Schwartz and Smith [46]. The extension to OU-based models driven by a
normal inverse Gaussian (NIG) process or by a variance gamma process can be
found in Benth et al [6], Cummins et al. [15, 16] and Sabino [42]. Finally, a
third type of market models is based on mean-reverting jump-diffusion OU
processes as originally suggested in Cartea and Figueroa [11] and further
analyzed for instance in Kjaer [30], Humbly et al [25] and recently in Sabino
and Cufaro Petroni [44].

On the other hand, in order to select a realistic and viable model for
derivative pricing, one has to first answer the following questions.

  * •

Do we know the characteristic function and the moment generating function of
$X(t)$ to determine the non-arbitrage conditions.

  * •

Do we know anything about the statistics of the process $X(\cdot)$? How viable
is the parameters estimation?

  * •

Do we know how to generate the trajectories of $X(\cdot)$ in order to
implement Monte Carlo simulations?

The availability of simulation techniques of easy implementation is important
for analysis, validation and estimation purposes. Indeed, direct likelihood
analysis is often impracticable, whereas Monte Carlo (MC) based techniques and
generalized method of moments (GMM) approaches can be a viable route to
estimate the model parameters.

In this paper we study normal tempered stable (NTS) processes, which
generalize NIG processes, and focus on their OU counterpart. Benth and Benth
[9], Benth et al. [6, 7] illustrate the suitability of NIG-driven OU models in
various contexts, namely in gas, oil, power markets and also in the pricing of
wind derivatives. On the other hand, these applications rely on approximated
solutions to the answers mentioned above that may produce biased results, for
instance in the pricing of forward start contracts as we shall show.

We address each of the aforementioned questions for models that are based on
symmetric normal tempered stable driven OU processes. The first contribution
then is the derivation of the characteristic function of their transition law
in closed form. This result is instrumental to determine their statistical
properties and compared to the current state of affairs (see once again Benth
and Benth [9] and Benth et al. [6]) we can derive non-arbitrage conditions
without resorting to a numerical approximation or a numerical integration.

The second contribution consists in defining an efficient algorithm for the
exact simulation of the trajectories of symmetric normal tempered stable (and
therefore NIG) driven OU processes which is particularly suitable for forward
start contracts where the standard Euler scheme would return biased results.
Indeed, even though it is common practice to rely on such an approximated
scheme, we show it can lead to biased results. We also propose an alternative
approximation scheme that apparently outperforms the Euler approximation
having the same computational effort.

We illustrate the applicability of our 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 the first example
we consider a one-factor model purely driven by a single mean-reverting OU-NTS
process basically coinciding with that in Benth and Benth [9] and Benth et al.
[6]. We make use of the explicit knowledge of the characteristic function to
implement the pricing with the FFT-based technique of Carr and Madan [10], and
then compare the outcomes to those obtained via MC simulations. We observe
that such a model is not really suitable for the pricing of contracts having
long maturities which in contrast require a second NTS factor. The resulting
dynamics is a NTS version of the model of de Jong and Schneider [17] in the
Gaussian framework.

In the second example, we calibrate the two factor model taking the German
day-ahead and month-ahead NCG prices and evaluate a forward start Asian
option. Based on our theoretical results, we suggest to estimate the
parameters using the alternative approximation, whereas to employ the exact
simulation scheme for the pricing of energy derivatives especially if they are
forward start contracts. In addition, we have shown that the proposed
simulation algorithm, combined with the Least Squares Monte Carlo approach of
Hambly et al [25], provides an efficient and accurate pricing of a one year
$120-120$ swing option.

Finally, our results are not only restricted to OU-processes and spot models.
Indeed, we shall detail how they can be easily extended to forward markets in
order to capture the Samuelson effect and different implied volatility
profiles displayed by options futures in the spirit of the works of Benth et
al. [8], Latini et al. [31] and Piccirilli et al. [39]. The paper is organized
as follows: Section 2 introduces the notation and presents the mathematical
preliminaries. In Section 3 we focus on the symmetric normal tempered stable
processes of OU type, we derive the characteristic function of their
transition law, and present the algorithms for their exact simulation. We also
carry out numerical experiments demonstrating their efficiency. The financial
application of these results is illustrated in Section 4 in the context of the
pricing of energy derivative contracts, namely daily strips of call options,
Asian options with European exercise and swing options written on the day-
ahead spot price. Furthermore, we detail how our findings are not restricted
to the modeling of the spot dynamics and can be easily extended to future
markets. Finally Section 5 concludes the paper with an overview of future
investigations and possible further applications.

##  2 Notations and preliminary remarks

In this section we present the concepts of normal variance-mean mixtures and
of non-Gaussian Ornstein-Uhlenbeck processes that will be instrumental for the
modeling of energy markets. We also introduce the notation and the shortcuts
that will be used throughout the paper.

###  2.1 Notation

The function $\Gamma(x)$ represents the Euler gamma function, in addition, we
write $\mathcal{N}(\mu,\sigma^{2})$ to denote the Gaussian distribution with
mean $\mu$ and variance $\sigma^{2}$. 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.2 Normal variance-mean mixtures

Consider a standard normal rv $X\sim\mathcal{N}(0,1)$ and an independent rv
$V$ with pdf $f_{V}(x)$ defined on the positive real axis $\mathbb{R}^{+}$,
then a rv $Y$ is said to be distributed according to a _normal variance-mean
mixture_ with mixture pdf $f_{V}(x)$ if it has the form

$Y=\mu+\theta\,V+\sigma\sqrt{V}\,X$ (1)

where $\theta$, $\mu$ and $\sigma>0$ are real numbers. The conditional
distribution of $Y$ given $V$ is thus a Gaussian distribution
$\mathcal{N}(\mu+\theta V,\sigma^{2}\,V)$, whereas the unconditional pdf and
chf denoted $f_{Y}(x)$ and $\varphi_{Y}(u)$ respectively, are

$f_{Y}(x)=\int_{0}^{\infty}\frac{1}{\sqrt{2\,\pi\,\sigma^{2}\,v}}\exp\left(-\frac{\left(x-\mu-\theta\,v\right)^{2}}{2\,\sigma^{2}\,v}\right)dv$
(2)
$\varphi_{Y}(u)=e^{i\,u\,\mu}\varphi_{V}\left(\theta\,u+\frac{i\,\sigma^{2}\,u^{2}}{2}\right)$
(3)

where $\varphi_{V}(u)$ denotes the chf of $V$.

Assuming for simplicity $\mu=0$, the law of $Y$ coincides with the
distribution of the position of a subordinated Brownian motion (BM) where the
law of the subordinator at a fixed time is given by $f_{V}(x)$. We recall that
a subordinator $L(\cdot)$ is a pure jump Lévy process with non-decreasing
trajectories and Lévy measure $\nu_{L}$ such that

$\int_{0}^{\infty}(x\wedge 1)\nu_{L}(dx)<\infty,\quad\nu_{L}((-\infty,0))=0.$

Denoting $W(\cdot)$ the standard Wiener process, a subordinated BM

$Y(t)=\theta\,L(t)+\sigma\,W(L(t))$ (4)

is thus a new Lévy process with time-change given by $L(\cdot)$.

Furthermore, a Lévy process $L(\cdot)$ is said to be a classic tempered stable
(dubbed TS) subordinator if its Lévy measure has density

$\ell_{L}(x)=c\,\frac{e^{-\beta\,x}}{x^{1+\alpha}}\mathds{1}_{x\geq 0}$ (5)

where $c,\beta$ are all positive numbers, $0<\alpha<1$. Hereafter
$\mathcal{TS}(\alpha,\beta,c\,t)$ denotes the law of such a subordinator at
time $t$. This class is quite flexible and contains the inverse Gaussian (IG)
distribution that corresponds to $\alpha=1/2$ and the gamma distribution for
the limiting case $\alpha=0$. For sake of completeness, we remark that we are
referring to classical TS processes because alternative tempering functions
rather than the exponential function can be used in (5) (see Rosiński [40]).
An overview of such processes, named general tempered stable processes, can be
found in Grabchak [22].

By time-changing a BM with drift $\theta$ and volatility $\sigma$ with a TS
subordinator, we obtained the so-called _normal (classical) tempered stable_
process, hereafter dubbed NTS; taking $\alpha=1/2$ we have the well-known NIG
process. On the other hand, due to the scaling properties of the TS processes,
it is sufficient to consider subordinators with $\bm{E}\left[{L(t)}\right]=t$
in which case form a one-parameter family. Accordingly, it is convenient to
represent their Lévy density as follows (see Cont and Tankov [12] equation
4.19)

$\ell_{L}(x)=\frac{1}{\Gamma(1-\alpha)}\left(\frac{1-\alpha}{\nu}\right)^{1-\alpha}\,\frac{e^{-(1-\alpha)/\nu}}{x^{1+\alpha}}\mathds{1}_{x\geq
0}$ (6)

where $\nu=\bm{V}\left[{L(1)}\right]$, $\beta=\frac{(1-\alpha)}{\nu}$ and
$c=\frac{1}{\Gamma(1-\alpha)}\left(\frac{1-\alpha}{\nu}\right)^{1-\alpha}$.
Furthermore, denoting $\varphi_{Y}(u,t)$ the chf of $Y(t)$ and
$\psi_{Y}(u)=t^{-1}\log(\varphi_{Y}(u,t))$ the characteristic exponent (che)
of $Y(\cdot)$, based on (3) it results

$\psi_{Y}(u)=\frac{1-\alpha}{\alpha\,\nu}\left[1-\left(1+\frac{\nu(u^{2}\sigma^{2}/2-i\theta\,u)}{1-\alpha}\right)^{\alpha}\right].$
(7)

###  2.3 OU-TS processes

Consider 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$
(8)

namely,

$\displaystyle X(t)$ $\displaystyle=$ $\displaystyle
X_{0}\,e^{-bt}+Z_{L}(t),\qquad\quad Z_{L}(t)=\int_{0}^{t}e^{-b\,(t-s)}dL(s)$
(9)

where $L(\cdot)$ is a Lévy process, named _Background Driving Lévy Process_
(BDLP) accordingly. Following the convention in Barndorff-Nielsen and Shephard
[4], if $L(\cdot)$ is a TS subordinator $X(\cdot)$ is called OU-TS 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. [3], the solution process (9) 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 chf of the rv $Z(t)$ in (9) 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 that a law with chf $\eta(u)$ is said to be sd (see Sato [45] or
Cufaro Petroni [13]) 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).$ (10)

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}$ (11)

Hereafter the rv $Z_{a}$ will be called the _$a$ -remainder_ of $X$ and in
general has an id (see Sato [45]). This last statement apparently means that
the law of $Z_{L}(t)$ in the solution (9) 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 (9) that the lch $\psi_{X}(u,t)$
of the time homogeneous transition law with a degenerate initial condition
$X(0)=x_{0},\;\bm{P}\hbox{-\emph{a.s.}}$ is

$\psi_{X}(u,t|x_{0})=iux_{0}e^{-bt}+\psi_{Z_{L}}(u,t).$ (12)

and can also be written in terms of the corresponding che $\psi_{L}(u)$ in the
form

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

##  3 Symmetric OU-NTS processes and their exact simulation

In this section we focus on OU processes whose BDLP $Y(\cdot)$ is a symmetric
NTS process (hereafter dubbed SNTS) with $\theta=0$ and with $L(\cdot)$ a TS
subordinator. Although the law of the SNTS process at a fix time $t$ boils
down into a three-parameters distribution, in the following we will keep the
notation $\mathcal{SNTS}(\sigma,\alpha,\beta,c\,t)$ with four parameters under
the assumption that

$\beta=\frac{1-\alpha}{\nu},\quad\quad
c=\frac{1}{\Gamma(1-\alpha)}\left(\frac{1-\alpha}{\nu}\right)^{1-\alpha}=\frac{\beta^{1-\alpha}}{\Gamma(1-\alpha)}.$

In virtue of (9) we write

$\displaystyle N(t)$ $\displaystyle=$ $\displaystyle N(0)\,e^{-bt}+Z_{Y}(t),$
(14) $\displaystyle Z_{Y}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}e^{-b\,(t-s)}dY(s)=\sigma\int_{0}^{t}e^{-b\,(t-s)}dW(L(s))$
(15)

where $N(\cdot)$ denotes a OU-SNTS process relatively to which the following
proposition holds.

######  Proposition 3.1.

Denoting $a=e^{-b\,t}$ and $\omega=a^{2}$, the pathwise solution (14) with
$N(0)=N_{0},\bm{P}\hbox{-\emph{a.s.}}$, is in distribution the sum of three
independent rv’s

$N(t)=a\,N_{0}+Z_{Y}(t)\stackrel{{\scriptstyle d}}{{=}}a\,N_{0}+N_{1}+N_{2}$
(16)

where

$N_{1}\stackrel{{\scriptstyle d}}{{=}}\sigma\sqrt{M_{1}}\,X_{1}$ (17)

with $X_{1}$, $M_{1}$ independent, $X_{1}\sim\mathcal{N}(0,1)$ and
$M_{1}\sim\mathcal{TS}\\!\left(\alpha,\frac{\beta}{\omega},c\,\frac{1-\omega^{\alpha}}{2\,\alpha\,b}\right)$,
therefore $N_{1}$ is distributed according to a
$\mathcal{SNTS}\left(\sigma,\alpha,\frac{\beta}{\omega},c\,\frac{1-\omega^{\alpha}}{2\,\alpha\,b}\right)$
law. Moreover,

$N_{2}\stackrel{{\scriptstyle d}}{{=}}\sigma\sqrt{M_{2}}\,X_{2}$ (18)

with $X_{2}$, $M_{2}$ independent and independent of $X_{1}$ and $M_{1}$. In
its turn,

$M_{2}=\sum_{k=0}^{P_{\omega}}J_{k},\quad J_{0}=0,\bm{P}\hbox{-\emph{a.s.}}$

is a compound Poisson rv where $P_{\omega}$ is an independent Poisson rv with
parameter

$\Lambda_{\omega}=\frac{c\,\beta^{\alpha}\Gamma(1-\alpha)}{2\,b\,\alpha^{2}\omega^{\alpha}}\,\left(1-\omega^{\alpha}+\omega^{\alpha}\log\omega^{\alpha}\right)$
(19)

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

$f_{J}(x)=\frac{\alpha\,\omega^{\alpha}}{1-\omega^{\alpha}+\omega^{\alpha}\log\omega^{\alpha}}\int_{1}^{\frac{1}{\omega}}\frac{x^{-\alpha}\left(\beta\,v\right)^{1-\alpha}e^{-\beta
v\,x}}{\Gamma(1-\alpha)}\,\frac{v^{\alpha}-1}{v}\,dv$ (20)

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

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

###### Proof.

Based on the definition of the chf $\varphi_{Z_{Y}}(u,t)$ of $Z_{Y}(t)$ and on
the fact that $L(\cdot)$ and $W(\cdot)$ are independent, it results

$\displaystyle\varphi_{Z_{Y}}(u,t)$ $\displaystyle=$
$\displaystyle\bm{E}\left[{exp\left(i\,u\,\sigma^{2}\int_{0}^{t}e^{-\,b\,(t-s)}dW(L(s))\right)}\right]=$
(22)
$\displaystyle\bm{E}\left[{\bm{E}\left[{exp\left(i\,u\,\sigma^{2}\int_{0}^{t}e^{-\,b\,(t-s)}dW(L(s))\right)}\right]\Bigg{|}L(t)}\right]=$
$\displaystyle\bm{E}\left[{exp\left(-\frac{\sigma^{2}\,u^{2}}{2}\int_{0}^{t}e^{-2\,b\,(t-s)}dL(s)\right)}\right]=\zeta\left(\frac{i\,\sigma^{2}\,u^{2}}{2},t\right),$

where $\zeta(u,t)$ represents the chf of
$\tilde{Z}_{L}(t)=\int_{0}^{t}e^{-2\,b\,(t-s)}L(s)$, therefore $Z_{Y}(t)$ is
distributed according to a normal variance-mean mixture distribution whose
mixture law is that with chf $\zeta(u,t)$. On the other hand, from Proposition
5.1 in Cufaro Petroni and Sabino [14], $\tilde{Z}_{L}(t)$ can be written as
the sum of two independent rv’s

$\tilde{Z}_{L}(t)\stackrel{{\scriptstyle
d}}{{=}}M_{1}+M_{2},\quad\quad\zeta(u,t)=\zeta_{1}(u)\,\zeta_{2}(u)$

where $\zeta_{1}(u)$ is the chf of
$M_{1}\sim\mathcal{TS}\\!\left(\alpha,\frac{\beta}{\omega},c\,\frac{1-\omega^{\alpha}}{2\,\alpha\,b}\right)$,
whereas $\zeta_{2}(u)$ is the chf of $M_{2}$ in the hypothesis. Consequently
$\varphi_{Z_{Y}}(u,t)$ has the form

$\varphi_{Z_{Y}}(u,t)=\zeta_{1}\left(\frac{i\,\sigma^{2}\,u^{2}}{2}\right)\zeta_{2}\left(\frac{i\,\sigma^{2}\,u^{2}}{2}\right)$

hence

$Z_{Y}(t)\stackrel{{\scriptstyle d}}{{=}}N_{1}+N_{2}$

where $N_{1}$ and $N_{2}$ are two independent rv’s each distributed according
to a normal variance mixture law whose mixture laws are the distributions of
$M_{1}$ and $M_{2}$, respectively and that concludes the proof. ∎

######  Remark 1.

Because of (22), it also holds

$Z_{Y}(t)\stackrel{{\scriptstyle d}}{{=}}X\,\sqrt{M_{1}+M_{2}}$ (23)

where $X\sim\mathcal{N}(0,1)$ and it turns out to be a more convenient and a
computationally efficient way to simulate the skeleton of $Y(\cdot)$.

######  Proposition 3.2.

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

$\psi_{Z_{Y}}(u,t)=-\frac{c\,\beta^{\alpha}\,\Gamma(1-\alpha)}{2\,\alpha\,b}\left[J\left(u,\sigma,\alpha,\beta,\frac{\beta}{\omega}\right)+\log\omega\right],\quad\omega=e^{-2\,b\,t}$
(24)

with

$\displaystyle J\left(u,\sigma,\alpha,\beta_{1},\beta_{2}\right)$
$\displaystyle=$
$\displaystyle\int_{\beta_{1}}^{\beta_{2}}z^{-1-\alpha}\left(z+\frac{\sigma^{2}}{2}\,u\right)^{\alpha}dz=$
(25) $\displaystyle=$
$\displaystyle-\frac{1}{\alpha}\left[\left(\frac{\sigma^{2}\,u^{2}}{2\,\beta_{2}}\right)^{\alpha}\,_{2}F_{1}\left(-\alpha,-\alpha,1-\alpha,-\frac{2\,\beta_{2}}{\sigma^{2}\,u^{2}}\right)\right.-$
$\displaystyle\left.\left(\frac{\sigma^{2}\,u^{2}}{2\,\beta_{1}}\right)^{\alpha}\,_{2}F_{1}\left(-\alpha,-\alpha,1-\alpha,-\frac{2\,\beta_{1}}{\sigma^{2}\,u^{2}}\right)\right]$

where ${}_{2}F_{1}(a,b,c,x)$ is the hypergeometric function, $0<\alpha<1$,
$\beta_{1}>0$ and $\beta_{2}>\beta_{1}$. Finally, taking $\alpha=1/2$ the
hypergeometric function in (25) 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.

Because of Proposition 3.2 in Sabino [41], the lch $\rho(u,t)$ of
$\tilde{Z}_{L}(t)=\int_{0}^{t}e^{-2\,b\,(t-s)}L(s)$ is

$\rho(u,t)=-\frac{c\,\beta^{\alpha}\,\Gamma(1-\alpha)}{2\,\alpha\,b}\left[I\left(u,\alpha,\beta,\frac{\beta}{\omega}\right)+\log\omega\right],\quad\omega=e^{-2\,b\,t}$
(26)

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=$
(27) $\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].$

On the other hand, we know from (22) that
$\varphi_{Z_{Y}}(u,t)=\zeta\left(\frac{i\,\sigma^{2}\,u^{2}}{2},t\right)$
therefore, simply adapting this form to the lch $\rho(u,t)$ concludes the
proof. ∎

######  Remark 2.

Using the transformation 9.131.1 in Gradshteyn and Rizhik [23] we have

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

and accordingly (25) 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.-$
(28)
$\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].$

We make use of this different representation in the following corollary.

######  Corollary 3.3.

The cgf $m_{Z_{Y}}(s,t)=\ln\bm{E}\left[{e^{s\,Z_{Y}(t)}}\right]$ of $Z_{Y}(t)$
with $0<\alpha<1$ exists for
$-\frac{\sqrt{2\,\beta}}{\sigma}<s<\frac{\sqrt{2\,\beta}}{\sigma}$ and is:

$m_{Z_{Y}}(s,t)=-\frac{c\,\beta^{\alpha}\,\Gamma(1-\alpha)}{2\,\alpha\,b}\left[\tilde{J}\left(s,\sigma,\alpha,\beta,\frac{\beta}{\omega}\right)+\log\omega\right]$
(29)

where

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

###  3.1 Simulation Algorithms

Proposition 3.1 provides the theoretical basis for the sequential generation
of the skeleton of $N(\cdot)$ on a time grid $t_{0},t_{1},\dots t_{M}$, that,
assuming at each step $a_{m}=e^{-b(t_{m}-t_{m-1})},m=0,\dots,M$, can be
accomplished by the following recursion with initial condition
$N(t_{0})=n_{0}$:

$N(t_{m})=a_{m}N(t_{m-1})+Z_{Y}(t_{m-1}),\quad\qquad m=1,\dots,M,$ (31)

the implementation steps are also detailed in Algorithm 1.

Algorithm 1

1$N(t_{0})\leftarrow n_{0}$.

2for $m=1,\dots,M$ do

3 $a_{m}=e^{-b(t_{m}-t_{m-1})}$

4 $\omega_{m}=a_{m}^{2}=e^{-2\,b(t_{m}-t_{m-1})}$

5 $x\leftarrow X\sim\mathcal{N}(0,1)$.

6 $m_{1}\leftarrow
M_{1}\sim\mathcal{TS}\\!\left(\alpha,\frac{\beta}{\omega_{m}},c\,\frac{1-\omega_{m}{\alpha}}{2\,\alpha\,b}\right)$

7 $p\leftarrow P_{\omega_{m}}\sim\mathcal{P}(\Lambda_{\omega_{m}})$,
$\triangleright$ Generate an independent Poisson rv with
$\Lambda_{\omega_{m}}$ in (19)

8 $v_{i}\leftarrow V_{i},\quad i=1,\dots,p$$\triangleright$ Generate
independent rv’s with pdf (21)

9 $\tilde{\beta}_{i}\leftarrow\beta\,v_{i},\quad i=1,\dots,p$

10 $j_{i}\leftarrow J_{i}\sim\mathcal{G}(1-\alpha,\,\tilde{\beta}_{i}),\quad
i=1,\dots,n$ $\triangleright$ Generate $p$ independent gamma rv’s all with the
same scale $1-\alpha$ and random rates

11 $m_{2}\leftarrow\sum_{i=1}^{p}j_{i}$

12 $z_{m}\leftarrow x\,\sqrt{m_{1}+m_{2}}$

13 $N(t_{m})\leftarrow a_{i}\,N(t_{m-1})+z_{m}$.

14end for

The generation of the rv’s in Algorithm 1 is standard except that for the TS
distributed rv $M_{1}$ and for $V$. To this end, the sampling from a TS law
has been widely studied by several authors (see for instance Devroye [19] and
Hofert [26]), of course, taking $\alpha=0.5$ the TS law coincides with an IG
law hence, one can rely on the many-two-one transformation method of Michael
et al. [37]. With regards to the simulation of $V$, an efficient algorithm is
detailed in Cufaro Petroni and Sabino [14] which is based on the
decomposition-rejection method illustrated in Devroye [18] page 67.
Furthermore, we remark that using a Taylor expansion for $\Lambda_{\omega}$ we
find

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

and therefore the compound Poisson $M_{2}$ has a relevant impact only when
$\Delta t$ is not too small. This observation gives then some hints on
possible approximations. We will compare the performance of our approach to
two alternatives: the first simply ignores $M_{2}$ (Approximation 1), the
second, in the same vein of Benth et al. [7], relies on the Euler scheme
applied to the stochastic differential equation of an OU-SNTS process, namely
it approximates the law of $Z_{Y}(\Delta t)$ with that of $L_{Y}(\Delta
t)\sim\mathcal{SNTS}\left(\sigma,\alpha,\beta,c\,\Delta t\right)$
(Approximation 2) hence, in both cases $Z_{Y}(\Delta t)$ is approximated by a
SNTS distributed rv.

###  3.2 Numerical Experiments

In this section, we analyze the performance and the effectiveness of the
algorithms for the simulation of OU-SNTS processes. All the numerical studies
in this paper have been conducted using _Python_ with a $64$-bit Intel Core
i5-6300U CPU, 8GB. The performance of the algorithms is measured in terms of
the percentage error relatively to the second and fourth cumulants denoted
_err_ % and defined as

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

As a consequence of Equation (13) one can calculate the cumulants
$c_{N,k}(n_{0},t),\;k=1,2,\ldots$ of $N(t)$ for $N_{0}=n_{0}$ from the
cumulants $c_{L_{Y},k}$ of the symmetric NTS law according to

$\displaystyle c_{N,1}(n_{0},t)$ $\displaystyle=$
$\displaystyle\bm{E}\left[{N(t)|N_{0}=n_{0}}\right]\;=\;n_{0}e^{-b\,t}+\frac{c_{L_{Y},1}}{b}\left(1-e^{-b\,t}\right),\qquad
k=1$ (32) $\displaystyle c_{N,k}(n_{0},t)$ $\displaystyle=$
$\displaystyle\frac{c_{L_{Y},k}}{k\,b}\left(1-e^{-k\,b\,t}\right)\qquad\qquad\qquad\qquad\qquad
k=2,3,\ldots$ (33)

On the other hand, $c_{L_{Y},k}=0$ when $k$ is an odd integer and setting
$k=2\,n,n>0$ $c_{L_{Y},2\,n}$ can be computed from the Lévy density as follows

$\displaystyle c_{L_{Y},2\,n}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}x^{2\,n}\ell_{L_{Y}}(x)\,dx=$ (34)
$\displaystyle
C(\alpha,\nu,\sigma)\int_{-\infty}^{+\infty}x^{2\,n}|x|^{-\alpha-1/2}K_{\alpha+1/2}\left(A(\alpha,\nu,\sigma)\,|x|\right)\,dx=$
$\displaystyle
C(\alpha,\nu,\sigma)2^{2\,n-\alpha-1/2}\,A(\alpha,\nu,\sigma)^{-2\,n-\alpha-1/2}\,\Gamma\left(\frac{2\,n+1}{2}\right)\,\Gamma\left(n-\alpha\right),$

where $K_{\eta}(x)$ represents the modified Bessel function of the second
type,

$\displaystyle C(\alpha,\nu,\sigma)$ $\displaystyle=$
$\displaystyle\frac{2^{\alpha/2+5/4}}{\sqrt{2\,\pi}\,\Gamma(1-\alpha)}\sigma^{\alpha-1/2}\beta^{5/4-\alpha/2},$
$\displaystyle A(\alpha,\nu,\sigma)$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2\,\beta}}{\sigma}$

and in the last step we have used 6.561.16 at page 676 in Gradshteyn and
Rizhik [23]. Of course, with $\alpha=1/2$ we retrieve the cumulants of a
symmetric NIG law shown in Table 4.5 in Cont and Tankov [12].

Figure 1: Sample trajectories of a OU-SNTS process with
$\left(b,\sigma,\nu\right)=\left(5,0.3,2.5\right)$.

Tables 1 and 2 show the comparison between the true values of
$c_{N,2}(n_{0},\Delta t)$ and $c_{N,4}(n_{0},\Delta t)$ and their
corresponding MC estimates. We have considered $n_{0}=0$,
$\left(b,\sigma,\nu\right)=\left(5,0.3,2.5\right)$ with $\Delta t=1/365$ and
$\Delta t=1/12$, respectively, and let $\alpha$ vary; similar conclusions have
been derived with different parameter sets but are not reported for brevity.
Moreover, Figure 1 displays some path trajectories changing the parameter
$\alpha$.

Irrespective to the values of $\Delta t$ and $\alpha$, one can observe that
the convergence is relatively slow and a large number of simulation is
required to get reliable estimates. Indeed, with the exact scheme $10^{5}$
simulations are needed to achieve roughly five percent of percentage error
both for the second and forth cumulant.

Besides this fact, if $\Delta t$ equals $1/365$ (one day in financial
applications) the approximated schemes and the exact scheme return similar
estimates with little bias, although Approximation 1 seems better than
Approximation 2. It is well known that the quality of the MC estimation
deteriorates with the increase of the order of the cumulants, indeed in our
numerical example we need $10^{6}$ scenarios to get a small percentage error.

In contrast, when $\Delta t=1/12$ \- e.g. one month - one can observe in Table
2 that the estimates of the approximated schemes are totally biased and no
matter how high is the number of simulations, the MC estimates do not converge
to the true values. On the other hand, the exact scheme performs well and
exhibits a similar behavior to the one with a small $\Delta t$.

These observations give a few hints on how to conceive the parameter
calibration for a OU-SNTS process and provide guidance for the path
generation. For instance, one can take a time series with relatively small
frequency, e.g. daily, and estimate the parameters assuming that the
transition density is that of SNTS according to Approximation 1. In particular
for $\alpha=1/2$, the law is a symmetric NIG and some software packages based
on the maximum-likelihood method are available to accomplish the parameter
estimation222See for instance, the package _GeneralizedHyperbolic_ in R.. On
the other hand, if the skeleton of the process is to be generated on a not so
fine time grid, the exact scheme is the only valid alternative. More
important, the pricing of forward start options, for instance those contracts
whose first settlement date is in one month, do not require the simulation of
the process before the first settlement date, leading to a remarkable
computational advantage even though the exact scheme relies on an additional
step compared to the two approximated procedure. We will comment more of these
facts in the following sections dedicated to the financial applications.

$\alpha=0.1$
---
| $c_{N,2}(n_{0}=0,\Delta\,t)=2.466\times 10^{-4}$ | $c_{N,2}c_{N,n}(n_{0}=0,\Delta\,t)=1.641\times 10^{-4}$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err %
$1000$ | $2.972$ | $-20.5\%$ | $3.524$ | $-42.9\%$ | $1.344$ | $45.5\%$ | $0.785$ | $52.2\%$ | $0.914$ | $44.3\%$ | $0.935$ | $43.0\%$
$10000$ | $2.043$ | $17.2\%$ | $3.268$ | $-32.5\%$ | $3.313$ | $-34.3\%$ | $0.796$ | $51.5\%$ | $0.888$ | $45.9\%$ | $0.918$ | $44.1\%$
$20000$ | $2.233$ | $9.4\%$ | $2.048$ | $17.0\%$ | $1.813$ | $26.5\%$ | $1.284$ | $21.8\%$ | $1.908$ | $-16.3\%$ | $2.076$ | $-26.5\%$
$50000$ | $2.577$ | $-4.5\%$ | $2.109$ | $14.5\%$ | $2.788$ | $-13.1\%$ | $1.384$ | $15.7\%$ | $1.464$ | $10.8\%$ | $2.032$ | $-23.8\%$
$100000$ | $2.503$ | $-1.5\%$ | $2.320$ | $5.9\%$ | $2.580$ | $-4.6\%$ | $1.464$ | $10.8\%$ | $1.494$ | $9.0\%$ | $1.821$ | $-11.0\%$
$1000000$ | $2.483$ | $-0.7\%$ | $2.515$ | $-2.0\%$ | $2.418$ | $2.0\%$ | $1.684$ | $-2.6\%$ | $1.550$ | $5.5\%$ | $1.532$ | $6.6\%$
$\alpha=0.3$
| $c_{N,2}(n_{0}=0,\Delta\,t)=2.466\times 10^{-4}$ | $c_{N,2}c_{N,n}(n_{0}=0,\Delta\,t)=1.641\times 10^{-4}$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err %
$1000$ | $0.196$ | $92.1\%$ | $1.009$ | $59.1\%$ | $1.249$ | $49.4\%$ | $0.665$ | $59.4\%$ | $0.187$ | $88.6\%$ | $0.890$ | $45.7\%$
$10000$ | $2.915$ | $-18.2\%$ | $3.299$ | $-33.8\%$ | $3.299$ | $-33.8\%$ | $0.957$ | $41.7\%$ | $2.261$ | $-37.8\%$ | $0.957$ | $41.7\%$
$20000$ | $2.717$ | $-10.2\%$ | $2.006$ | $18.6\%$ | $3.120$ | $-26.5\%$ | $1.969$ | $-20.0\%$ | $1.980$ | $-20.7\%$ | $2.030$ | $-23.7\%$
$50000$ | $2.293$ | $7.0\%$ | $2.171$ | $12.0\%$ | $2.096$ | $15.0\%$ | $1.380$ | $15.9\%$ | $1.361$ | $17.1\%$ | $1.892$ | $-15.3\%$
$100000$ | $2.341$ | $5.1\%$ | $2.619$ | $-6.2\%$ | $2.239$ | $9.2\%$ | $1.519$ | $7.4\%$ | $1.871$ | $-14.0\%$ | $1.815$ | $-10.6\%$
$1000000$ | $2.424$ | $1.7\%$ | $2.513$ | $-1.9\%$ | $2.384$ | $3.3\%$ | $1.603$ | $2.3\%$ | $1.607$ | $2.1\%$ | $1.515$ | $7.7\%$
$\alpha=0.5$
| $c_{N,2}(n_{0}=0,\Delta\,t)=2.466\times 10^{-4}$ | $c_{N,2}c_{N,n}(n_{0}=0,\Delta\,t)=1.641\times 10^{-4}$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err %
$1000$ | $0.514$ | $79.1\%$ | $4.215$ | $-70.9\%$ | $1.437$ | $41.7\%$ | $0.617$ | $62.4\%$ | $0.110$ | $93.3\%$ | $2.991$ | $-82.3\%$
$10000$ | $1.552$ | $37.1\%$ | $1.275$ | $48.3\%$ | $3.332$ | $-35.1\%$ | $2.265$ | $-38.0\%$ | $2.381$ | $-45.1\%$ | $2.389$ | $-45.6\%$
$20000$ | $2.185$ | $11.4\%$ | $3.117$ | $-26.4\%$ | $3.032$ | $-22.9\%$ | $1.265$ | $22.9\%$ | $1.216$ | $25.9\%$ | $2.030$ | $-23.7\%$
$50000$ | $2.293$ | $7.0\%$ | $2.171$ | $12.0\%$ | $2.171$ | $12.0\%$ | $1.915$ | $-16.7\%$ | $1.892$ | $-15.3\%$ | $1.890$ | $-15.2\%$
$100000$ | $2.341$ | $5.1\%$ | $2.247$ | $8.9\%$ | $2.247$ | $8.9\%$ | $1.761$ | $-7.3\%$ | $1.471$ | $10.3\%$ | $1.510$ | $8.0\%$
$1000000$ | $2.513$ | $-1.9\%$ | $2.531$ | $-2.6\%$ | $2.350$ | $4.7\%$ | $1.698$ | $-3.5\%$ | $1.699$ | $-3.5\%$ | $1.516$ | $7.6\%$
$\alpha=0.3$
| $c_{N,2}(n_{0}=0,\Delta\,t)=2.466\times 10^{-4}$ | $c_{N,2}c_{N,n}(n_{0}=0,\Delta\,t)=1.641\times 10^{-4}$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err %
$1000$ | $0.923$ | $62.6\%$ | $1.264$ | $48.7\%$ | $1.276$ | $48.3\%$ | $0.010$ | $99.4\%$ | $0.532$ | $67.6\%$ | $0.377$ | $77.0\%$
$10000$ | $1.557$ | $36.9\%$ | $1.910$ | $22.5\%$ | $1.792$ | $27.3\%$ | $0.816$ | $50.3\%$ | $2.389$ | $-45.6\%$ | $2.168$ | $-32.1\%$
$20000$ | $2.185$ | $11.4\%$ | $1.881$ | $23.7\%$ | $1.850$ | $25.0\%$ | $1.620$ | $1.3\%$ | $1.222$ | $25.5\%$ | $1.317$ | $19.8\%$
$50000$ | $2.342$ | $5.0\%$ | $2.006$ | $18.6\%$ | $2.085$ | $15.5\%$ | $1.333$ | $18.8\%$ | $1.322$ | $19.4\%$ | $1.377$ | $16.1\%$
$100000$ | $2.559$ | $-3.8\%$ | $2.227$ | $9.7\%$ | $2.714$ | $-10.1\%$ | $1.560$ | $5.0\%$ | $1.652$ | $-0.7\%$ | $1.749$ | $-6.6\%$
$1000000$ | $2.440$ | $1.0\%$ | $2.407$ | $2.4\%$ | $2.375$ | $3.7\%$ | $1.639$ | $0.1\%$ | $1.654$ | $-0.8\%$ | $1.627$ | $0.9\%$
$\alpha=0.3$
| $c_{N,2}(n_{0}=0,\Delta\,t)=2.466\times 10^{-4}$ | $c_{N,2}c_{N,n}(n_{0}=0,\Delta\,t)=1.641\times 10^{-4}$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err %
$1000$ | $1.276$ | $48.3\%$ | $1.412$ | $42.7\%$ | $1.451$ | $41.2\%$ | $0.902$ | $45.1\%$ | $0.909$ | $44.6\%$ | $0.190$ | $88.4\%$
$10000$ | $2.006$ | $18.6\%$ | $2.170$ | $12.0\%$ | $1.805$ | $26.8\%$ | $2.190$ | $-33.5\%$ | $0.506$ | $69.2\%$ | $2.381$ | $-45.1\%$
$20000$ | $2.233$ | $9.4\%$ | $2.089$ | $15.3\%$ | $1.905$ | $22.8\%$ | $2.061$ | $-25.6\%$ | $2.244$ | $-36.8\%$ | $1.835$ | $-11.8\%$
$50000$ | $2.288$ | $7.2\%$ | $2.694$ | $-9.2\%$ | $2.066$ | $16.2\%$ | $1.835$ | $-11.8\%$ | $1.844$ | $-12.4\%$ | $1.279$ | $22.1\%$
$100000$ | $2.350$ | $4.7\%$ | $2.398$ | $2.8\%$ | $2.216$ | $10.1\%$ | $1.730$ | $-5.4\%$ | $1.744$ | $-6.3\%$ | $1.491$ | $9.2\%$
$1000000$ | $2.502$ | $-1.5\%$ | $2.519$ | $-2.1\%$ | $2.566$ | $-4.0\%$ | $1.699$ | $-3.5\%$ | $1.652$ | $-0.6\%$ | $1.601$ | $2.4\%$

Table 1: MC-estimated cumulants multiplied by $10000$ and comparison to their
true values. $\left(b,\sigma,\nu\right)=\left(5,0.3,2.5\right)$, $\Delta
t=1/365$.

$\alpha=0.1$
---
| $c_{N,2}(n_{0}=0,\Delta\,t)=5.089\times 10^{-3}$ | $c_{N,2}c_{N,n}(n_{0}=0,\Delta\,t)=2.461\times 10^{-3}$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err %
$1000$ | $6.154$ | $-20.9\%$ | $2.350$ | $53.8\%$ | $2.702$ | $46.9\%$ | $2.951$ | $-19.9\%$ | $0.216$ | $91.2\%$ | $0.370$ | $85.0\%$
$10000$ | $5.156$ | $-1.3\%$ | $3.659$ | $28.1\%$ | $2.780$ | $45.4\%$ | $2.832$ | $-15.1\%$ | $1.018$ | $58.7\%$ | $0.501$ | $79.6\%$
$20000$ | $5.067$ | $0.4\%$ | $3.491$ | $31.4\%$ | $2.845$ | $44.1\%$ | $2.743$ | $-11.5\%$ | $0.939$ | $61.8\%$ | $0.564$ | $77.1\%$
$50000$ | $5.073$ | $0.3\%$ | $3.464$ | $31.9\%$ | $3.158$ | $37.9\%$ | $2.616$ | $-6.3\%$ | $0.937$ | $61.9\%$ | $0.842$ | $65.8\%$
$100000$ | $5.080$ | $0.2\%$ | $3.544$ | $30.4\%$ | $3.179$ | $37.5\%$ | $2.520$ | $-2.4\%$ | $1.026$ | $58.3\%$ | $0.873$ | $64.5\%$
$1000000$ | $5.083$ | $0.1\%$ | $3.419$ | $32.8\%$ | $3.266$ | $35.8\%$ | $2.426$ | $1.4\%$ | $0.962$ | $60.9\%$ | $1.010$ | $59.0\%$
$\alpha=0.3$
| $c_{N,2}(n_{0}=0,\Delta\,t)=5.089\times 10^{-3}$ | $c_{N,2}c_{N,n}(n_{0}=0,\Delta\,t)=2.461\times 10^{-3}$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err %
$1000$ | $5.723$ | $-12.5\%$ | $4.456$ | $12.4\%$ | $4.472$ | $12.1\%$ | $2.183$ | $11.3\%$ | $1.240$ | $49.6\%$ | $1.874$ | $23.9\%$
$10000$ | $5.827$ | $-14.5\%$ | $3.703$ | $27.2\%$ | $3.169$ | $37.7\%$ | $2.742$ | $-11.4\%$ | $0.826$ | $66.4\%$ | $0.659$ | $73.2\%$
$20000$ | $5.480$ | $-7.7\%$ | $3.824$ | $24.9\%$ | $3.275$ | $35.6\%$ | $2.249$ | $8.6\%$ | $0.953$ | $61.3\%$ | $1.171$ | $52.4\%$
$50000$ | $5.348$ | $-5.1\%$ | $3.624$ | $28.8\%$ | $3.274$ | $35.7\%$ | $2.590$ | $-5.2\%$ | $0.943$ | $61.7\%$ | $0.928$ | $62.3\%$
$100000$ | $5.050$ | $0.8\%$ | $3.668$ | $27.9\%$ | $3.495$ | $31.3\%$ | $2.529$ | $-2.8\%$ | $1.083$ | $56.0\%$ | $1.046$ | $57.5\%$
$1000000$ | $5.076$ | $0.2\%$ | $3.686$ | $27.6\%$ | $3.247$ | $36.2\%$ | $2.450$ | $0.5\%$ | $1.071$ | $56.5\%$ | $0.945$ | $61.6\%$
$\alpha=0.5$
| $c_{N,2}(n_{0}=0,\Delta\,t)=5.089\times 10^{-3}$ | $c_{N,2}c_{N,n}(n_{0}=0,\Delta\,t)=2.461\times 10^{-3}$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err %
$1000$ | $5.710$ | $-12.2\%$ | $5.651$ | $-11.0\%$ | $2.320$ | $54.4\%$ | $2.742$ | $-11.4\%$ | $6.369$ | $-158.8\%$ | $0.177$ | $92.8\%$
$10000$ | $5.633$ | $-10.7\%$ | $4.585$ | $9.9\%$ | $3.445$ | $32.3\%$ | $3.073$ | $-24.8\%$ | $1.954$ | $20.6\%$ | $1.388$ | $43.6\%$
$20000$ | $5.430$ | $-6.7\%$ | $4.393$ | $13.7\%$ | $3.635$ | $28.6\%$ | $2.678$ | $-8.8\%$ | $1.478$ | $39.9\%$ | $1.353$ | $45.0\%$
$50000$ | $5.260$ | $-3.4\%$ | $4.075$ | $19.9\%$ | $3.350$ | $34.2\%$ | $2.616$ | $-6.3\%$ | $1.128$ | $54.2\%$ | $1.038$ | $57.8\%$
$100000$ | $5.050$ | $0.8\%$ | $4.034$ | $20.7\%$ | $3.400$ | $33.2\%$ | $2.367$ | $3.8\%$ | $1.055$ | $57.1\%$ | $1.125$ | $54.3\%$
$1000000$ | $5.104$ | $-0.3\%$ | $4.043$ | $20.6\%$ | $3.245$ | $36.2\%$ | $2.467$ | $-0.3\%$ | $1.124$ | $54.3\%$ | $0.979$ | $60.2\%$
$\alpha=0.3$
| $c_{N,2}(n_{0}=0,\Delta\,t)=5.089\times 10^{-3}$ | $c_{N,2}c_{N,n}(n_{0}=0,\Delta\,t)=2.461\times 10^{-3}$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err %
$1000$ | $5.545$ | $-9.0\%$ | $4.008$ | $21.2\%$ | $3.118$ | $38.7\%$ | $2.008$ | $18.4\%$ | $0.268$ | $89.1\%$ | $0.290$ | $88.2\%$
$10000$ | $4.702$ | $7.6\%$ | $4.576$ | $10.1\%$ | $3.322$ | $34.7\%$ | $2.742$ | $-11.4\%$ | $0.930$ | $62.2\%$ | $0.776$ | $68.4\%$
$20000$ | $5.421$ | $-6.5\%$ | $4.471$ | $12.1\%$ | $3.255$ | $36.0\%$ | $2.195$ | $10.8\%$ | $0.797$ | $67.6\%$ | $0.703$ | $71.5\%$
$50000$ | $5.329$ | $-4.7\%$ | $4.472$ | $12.1\%$ | $3.179$ | $37.5\%$ | $2.671$ | $-8.5\%$ | $1.431$ | $41.9\%$ | $0.663$ | $73.1\%$
$100000$ | $5.050$ | $0.8\%$ | $4.477$ | $12.0\%$ | $3.208$ | $37.0\%$ | $2.513$ | $-2.1\%$ | $1.270$ | $48.4\%$ | $0.875$ | $64.4\%$
$1000000$ | $5.099$ | $-0.2\%$ | $4.370$ | $14.1\%$ | $3.266$ | $35.8\%$ | $2.485$ | $-1.0\%$ | $1.144$ | $53.5\%$ | $1.010$ | $59.0\%$
$\alpha=0.3$
| $c_{N,2}(n_{0}=0,\Delta\,t)=5.089\times 10^{-3}$ | $c_{N,2}c_{N,n}(n_{0}=0,\Delta\,t)=2.461\times 10^{-3}$
| Exact | Approximation 1 | Approximation 2 | Exact | Approximation 1 | Approximation 2
$N$ | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err % | MC |  err %
$1000$ | $4.659$ | $8.4\%$ | $3.891$ | $23.6\%$ | $2.061$ | $59.5\%$ | $2.942$ | $-19.5\%$ | $0.182$ | $92.6\%$ | $0.002$ | $99.9\%$
$10000$ | $4.975$ | $2.2\%$ | $4.191$ | $17.6\%$ | $2.909$ | $42.8\%$ | $2.833$ | $-15.1\%$ | $0.222$ | $91.0\%$ | $0.400$ | $83.8\%$
$20000$ | $4.848$ | $4.7\%$ | $4.651$ | $8.6\%$ | $3.014$ | $40.8\%$ | $2.678$ | $-8.8\%$ | $0.876$ | $64.4\%$ | $0.385$ | $84.4\%$
$50000$ | $4.970$ | $2.3\%$ | $4.773$ | $6.2\%$ | $3.102$ | $39.0\%$ | $2.616$ | $-6.3\%$ | $1.208$ | $50.9\%$ | $0.431$ | $82.5\%$
$100000$ | $5.022$ | $1.3\%$ | $4.718$ | $7.3\%$ | $3.259$ | $36.0\%$ | $2.367$ | $3.8\%$ | $1.013$ | $58.8\%$ | $0.706$ | $71.3\%$
$1000000$ | $5.101$ | $-0.2\%$ | $4.843$ | $4.8\%$ | $3.282$ | $35.5\%$ | $2.419$ | $1.7\%$ | $1.271$ | $48.3\%$ | $1.108$ | $55.0\%$

Table 2: MC-estimated cumulants multiplied by $1000$ and comparison to their
true values. $\left(b,\sigma,\nu\right)=\left(5,0.3,2.5\right)$, $\Delta
t=1/12$.

##  4 Financial Applications

The day-ahead (also called spot) price of power or gas and in general of
commodities exhibit mean-reversion, seasonality and spikes, this last feature
is particularly difficult to be modeled in a pure Gaussian world. Different
approaches have been investigated in order to somehow extend the classical
Gaussian framework introduced in Lucia and Schwarz [34] and Schwartz and Smith
[46]. Among others, Cartea and Figueroa [11], Meyer-Brandis and Tankov [36]
and Sabino and Cufaro Petroni [44] have studied mean-reverting jump-diffusions
to model sudden spikes, whereas Cummins et al. [15, 16], and Sabino [41, 42]
have considered variance gamma and in general tempered stable and CGMY
processes to price power or gas derivative contracts.

In this section we consider a setting similar to Benth and Benth [9] and Benth
et al. [6] where the evolution of the spot dynamics is driven by a OU-SNTS
process and eventually by a second factor. Based on the results presented in
Section 3, compared to Benth and Benth [9] and Benth et al. [6] we provide
non-arbitrage conditions with no numerical approximations and then we
illustrate the effectiveness our findings in the pricing energy derivatives,
namely a strip of call options using the FFT-based technique of Carr and Madan
[10], a forward start Asian option with MC simulations and a swing option with
the Least-Squares Monte Carlo (LSMC) presented in Hambly et al. [24]. Finally,
we describe how these findings can be adapted to forward markets as well in
order to capture the Samuelson effect in spirit of the works of Benth et al.
[8], Latini et al. [31] or Piccirilli et al. [39].

###  4.1 Call Options

Consider a spot price dynamics driven by the following one factor process

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

where $h(t)$ is a deterministic function, $F(0,t)$ is the forward curve
derived from quoted products and reflects the seasonality, whereas $N(t)$ is
the OU-SNTS process in (14); as already mentioned, this market model basically
coincides with that discussed in Benth and Benth [9] and Benth et al. [6] when
$\alpha=1/2$.

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

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

where $m_{N}(s,t)$ is the cgf of $N(t)$,
$m_{N}(s,t)=s\,N(0)e^{-bt}+m_{Z_{Y}}(s,t)$ and $m_{Z_{Y}}(s,t)$ is given by
(29) therefore, assuming for simplicity $N(0)=0$, we have

$\displaystyle h(t)$ $\displaystyle=$
$\displaystyle\frac{c\,\beta^{\alpha}\,\Gamma(1-\alpha)}{\alpha\,b}\left[\tilde{J}\left(1,\sigma,\alpha,\beta,\frac{\beta}{\omega}\right)+\log\omega\right]$
(37)

with $\frac{\sqrt{2\,\beta}}{\sigma}>1$ or
$\sqrt{\frac{2\,(1-\alpha)}{\sigma^{2}\,\nu}}>1$. Taking $\alpha=1/2$, the
integrals can be written in terms of the logarithmic function, indeed

$\displaystyle\tilde{J}\left(1,\frac{1}{2},\beta_{1},\beta_{2}\right)$
$\displaystyle=$
$\displaystyle\int_{\beta_{1}}^{\beta_{2}}z^{-\frac{3}{2}}\left(z-\frac{\sigma^{2}}{2}\right)^{\frac{1}{2}}dz=\int_{\frac{2\,\beta_{1}}{\sigma^{2}}}^{\frac{2\,\beta_{2}}{\sigma^{2}}}y^{-\frac{3}{2}}\left(y-1\right)^{\frac{1}{2}}dy$

where we used the change of variable $z=\frac{\sigma^{2}\,y}{2}$, in addition
using equation (35) in Sabino [41] we have

$\displaystyle\tilde{J}\left(1,\frac{1}{2},\beta_{1},\beta_{2}\right)$
$\displaystyle=$ $\displaystyle
2\,\left(\log\left(\sqrt{\frac{2\,\beta_{2}}{\sigma^{2}}}+\sqrt{\frac{2\,\beta_{2}}{\sigma^{2}}-1}\right)-\sqrt{\frac{2\,\beta_{2}-\sigma^{2}}{2\,\beta_{2}}}\right)-$
(38) $\displaystyle
2\,\left(\log\left(\sqrt{\frac{2\,\beta_{1}}{\sigma^{2}}}+\sqrt{\frac{2\,\beta_{1}}{\sigma^{2}}-1}\right)-\sqrt{\frac{2\,\beta_{1}-\sigma^{2}}{2\,\beta_{1}}}\right).$

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 normally encompasses monthly, quarterly and yearly maturities,
but is not very liquid and is generally offered by brokers.

As a consequence of Proposition 3.2, we can make use of the explicit form of
the chf of the of $\log\,S(t)=\log F(0,t)+h(t)+X(t)$ to compute the price of
such a strip of calls using the FFT-based technique of Carr and Madan [10].

For this specific example we consider a set with realistic parameters,
$(b,\nu,\sigma)=(10,0.7,0.2)$, similar to that of Cummins et al. [15] and
Gardini et al. [20] and let $\alpha$ vary. In addition, in order to better
highlight the dependency on the maturity and $\alpha$, we take a flat forward
curve $F(0,t)=20$, indeed the seasonality changes the moneyness of the strip
of call options and could shadow some features of the price model; therefore
we also assume an at-the-money strike $K=20$.

Table 3 shows the numerical results relatively to different $\alpha$’s and
different maturities spanning from one month to one year. Table 3 also
compares the results obtained with the FFT method and with the MC method with
$10^{5}$ simulations. As far as this last method is concerned, we also report
the standard errors defined as the sample standard deviation divided by the
square root of the relative number of simulations.

As expected, we observe that the FFT and the MC approaches return comparable
results, of course the FFT method is known to be faster. Moreover, as also
mentioned in Subsection 3.1, except for $\alpha=0.5$, which corresponds to an
IG law, the simulation relies on the acceptance rejection method and is thus
more computationally expensive. Moreover, the call options prices exhibit a
non-negligible monotonic and increasing dependence on the stability parameter
$\alpha$. Taking the values for $\alpha=0.5$ (OU-NIG) as the benchmark, given
the used set of remaining parameters, the difference is of circa plus or minus
fifteen percent. On the other hand, Figure 2(b) shows that the price
difference is non-negligible for almost at-the-money options, whereas it
disappears for deep out and in-the-money options. Besides these observations,
an important fact is visible in Figure 2(a) which shows how each term of the
strip depends on the time-to-maturity for an at the money-option: apparently,
the time-to-maturity does not to play a relevant role for values above one
month. The cause of this evident drawback is the fact that we are using a one-
factor model and we remark that we can identify this issue because we have
switched the seasonality off. Indeed, even though we are using a lower value
for the mean reversion level $b$ than those presented in other studies (see
for instance once again Benth and Benth [9], Benth et al. [6]), such a value
is large enough to “tail off” the long term variability, hence a second factor
is required for the pricing of contracts with long maturities. In the
following section, we will illustrate the calibration and application of a
two-factor model to the pricing of a forward start Asian option and a swing
option. Of course, the FFT and MC methods discussed in this section can be
easily adapted to a two-factor price dynamics as well.

| $\alpha=0.1$ | $\alpha=0.3$ | $\alpha=0.5$ | $\alpha=0.7$ | $\alpha=0.9$
---|---|---|---|---|---
| FFT | MC | FFT | MC | FFT | MC | FFT | MC | FFT | MC
$T=1/12$ | $3.3259$ | $3.281\pm 0.0335$ | $3.8392$ | $3.7348\pm 0.034$ | $4.5342$ | $4.3131\pm 0.0334$ | $5.3861$ | $5.1903\pm 0.0324$ | $6.5152$ | $6.31\pm 0.0319$
$T=3/12$ | $16.481$ | $16.294\pm 0.104$ | $18.017$ | $18.241\pm 0.106$ | $19.905$ | $19.569\pm 0.105$ | $22.184$ | $22.115\pm 0.102$ | $25.308$ | $25.328\pm 0.098$
$T=1/2$ | $38.078$ | $37.5626\pm 0.1738$ | $40.881$ | $40.405\pm 0.179$ | $44.305$ | $43.853\pm 0.177$ | $48.515$ | $48.691\pm 0.173$ | $54.445$ | $54.868\pm 0.166$
$T=8/12$ | $59.749$ | $59.069\pm 0.2242$ | $63.799$ | $64.553\pm 0.231$ | $68.745$ | $68.227\pm 0.229$ | $74.878$ | $75.248\pm 0.224$ | $83.606$ | $84.387\pm 0.216$
$T=1$ | $81.421$ | $80.9306\pm 0.2655$ | $86.716$ | $86.762\pm 0.274$ | $93.186$ | $92.575\pm 0.271$ | $101.24$ | $101.79\pm 0.27$ | $112.77$ | $112.94\pm 0.26$

Table 3: Strip of at-the money call option prices obtained with the FFT and
the MC method with $10^{5}$ simulations under the assumption that
$\left(b,\sigma,\nu\right)=\left(10,0.2,0.7\right)$.

(a) $c_{m}(K,t_{m}),m=1,\dots,360$

(b) Effect on the strike

###  4.2 A Two-factor Model

The previously defined one-factor model (35) can easily be turned into a two-
factor one, adding a second process, for instance a NIG, a Variance Gamma or
any another NTS process, obtaining a normal tempered stable version of the two
factors Gaussian model of Schwartz and Smith [46]. For the reasons illustrated
in the previous section, such a model is of course more realistic, the two
factors are often associated by practitioners to a short and a long term
variability. We consider the following evolution of the spot prices

$S(t)=F(0,t)\,e^{h_{1}(t)+h_{2}(t)+N_{1}(t)+N_{2}(t)}$ (39)

where $h_{1}(t)$, $N_{1}(\cdot)$ and $F(0,t)$ play the same role of as the
terms in (35), whereas $N_{2}(\cdot)$ is an independent NTS process of the
type of Equation (4), and $h_{2}(t)$, similarly to $h_{1}(t)$, is a
deterministic function related to the che of $N_{2}(\cdot)$ given by (7). We
denote $(\alpha_{1},\sigma_{1},\nu_{1},\theta_{1},b)$ and
$(\alpha_{2},\sigma_{2},\nu_{2},\theta_{2})$ the parameters of the two
processes respectively. Further, we adopt the point of view of de Jong and
Schneider [17] and assume that the process $N_{2}(\cdot)$ is related to the
month-ahead forward price $F_{MA}(t)$ at time $t$ obtained rolling the front-
month forward once expired. Writing then

$F_{MA}(t)=F_{MA}(0)e^{h_{2}(t)+N_{2}(t)}$ (40)

has the advantage to relate the process $N_{2}(\cdot)$ to observable prices
and at the same time to simplify the parameter estimation of (39).

The unknown parameters $(\sigma_{2},\nu_{2},\theta_{2})$ can be estimated with
the generalized method of moments or with simulation based approaches because,
except for $\alpha=0.5$ (the NIG process), the transition density of
$N_{2}(\cdot)$ is quite complicated and one cannot always rely on the maximum-
likelihood method. In the following, we simplify the framework and take
$\alpha_{1}=\alpha_{2}=0.5$, namely, we assume that $N_{1}(\cdot)$ is a OU-NIG
process, whereas $N_{2}(\cdot)$ is a plain NIG process. Such a simplification
allows us to accomplish the parameter estimation with the maximum likelihood
method implemented by the package _GeneralizedHyperbolic_ in R. Nevertheless,
in the following section, we will illustrate the sensitivity of the Asian
option value to different parameters $\alpha$. We consider the historical data
of the day-ahead and month-ahead forward prices of the German gas market NCG
in the four years time window from January 1st 2016 to December 31st 2019 for
a total of $1461$ values.

After calibrating the parameters of the process $N_{2}(\cdot)$ we proceed with
the estimation of those of the process $N_{1}(\cdot)$, referred as a short
term factor, removing the seasonality component from the day-ahead time-series
at first, using a decomposition with a double-cosine plus a linear trend (see
Figure 3(b)333We use the package _STL_ in R). Taking then the logarithm of
each value of the deseasonalized time-series of the day-ahead and of the
month-ahead NCG prices denoted $s_{k}$ and $f_{k}$, $k=1,\dots,1461$,
respectively, we have

$s_{k+1}-f_{k+1}=(s_{k}-f_{k})\,e^{-b\,\Delta t}+\epsilon_{k+1}$ (41)

where $\Delta t=1/365$ (one day) and $\epsilon_{k},k=1,\dots$ is distributed
according to the law of $Z_{Y}(\Delta t)$ of (15) when $Y(\cdot)$ is a NIG
process. From the sample statistics of the residual it turns out that their
distribution is symmetric, thus can rely on the modeling via a symmetric OU-
NIG process and we can neglect the parameter $\theta_{1}$. We remark that due
to Proposition 3.1, $\epsilon_{k},k=1,\dots$ is not a NIG distributed random
sequence but rather the convolution of a NIG law plus a compound Poisson law.
On the other hand, because of the discussion in Section 3.2, it is reasonable
to neglect the latter contribution if $\Delta t$ is small, as it is in our
case, so that one can accomplish the parameter estimation with the maximum-
likelihood method. However, as observed once again in Section 3.2,
Approximation 1 is more reliable than the plain Euler scheme, and therefore,
we adopt such a choice. For sake of completeness, the time-series of the
month-ahead forward contract is available on trading days only, whereas that
of the day-ahead is available every day. In order to overcome this
inconsistency, we create fictitious values for the month-ahead time-series on
non-trading days by a random sampling from a NIG law with the estimated
parameters $(\hat{\sigma}_{2},\hat{\nu}_{2},\hat{\theta}_{2})$.

Table 4 reports the results of the estimation procedure, whereas Figures 3(a)
and 3(c) show that the NIG distribution fit the data of both the month-ahead
and the short term component of the day-ahead series quite well. Nevertheless,
the approximation that we have adopted in the calibration of the parameters of
the latter factor is not always reliable in the context of the pricing of
energy derivatives especially if one deals with forward start contracts. To
this end, in order to avoid to introduce an intrinsic bias, one must implement
the simulation of the skeleton of the OU-NIG component (and OU-NTS) with the
exact simulation scheme or at least select such a procedure for those time
intervals that are larger than one day. In the following application, we adopt
such a strategy for the pricing of a forward start Asian option and a swing
option.

$\hat{\sigma_{1}}$ | $\hat{\nu_{1}}$ | $\hat{b}$ | $\hat{\sigma_{2}}$ | $\hat{\nu_{2}}$ | $\hat{\theta_{2}}$
---|---|---|---|---|---
$0.2835$ | $0.0804$ | $39.86$ | $0.3142$ | $0.1023$ | -0.019
Table 4: Estimated Parameters

(a) Month-ahead

(b) Detrend

(c) Short-term

###  4.3 Asian Options

As a second financial application we consider the pricing of Asian options
with MC simulations in a two-factor market model using the calibration of
Section 4.2 and let also $\alpha$ vary. Of course, the parameter estimation
carried out in Section 4.2 assumed that the two factors were both NIG
distributed with $\alpha=0.5$, one might change $\alpha$ for each factor. In
addition, one might argue that for a general NTS distributed factor the
remaining parameters should also be different that those with $\alpha=0.5$.
Nevertheless, we only let the stability parameter change for both the two
factors and study its impact on the option price.

It is worthwhile noting that knowing the chf’s of the two factors of the
market model, one could also price Asian options with FFT methods (see Zhang
and C. Oosterlee [48]) that are faster than MC methods. On the other hand, 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)^{+}.$

Our numerical example considers a forward start contract with $I=90$ that
starts settling at $T_{1}=1/4$ with maturity $T=1/2$ and $K=11.5$.

Since the forward start feature and the observations of Section3.2, we rely on
the exact simulation scheme only for the OU-SNTS component because, even
though the presence of the second factor could shadow the impact, one single
large time point could create a bias in the option price.

Figure 5 shows the prices and the errors defined as the sample standard
deviation divided by the number of scenarios $N$ changing the parameter
$\alpha$. We can conclude that the MC procedure is convergent but that at
least $N=5\times 10^{4}$ scenarios are required to attain a reliable price and
that $\alpha$ has a non-negligible impact. Indeed, the price difference with
two consecutive $\alpha$’s in Table 5 is of circa ten percent that is a large
number for an Asian option that settles over three months only.

| $\alpha=0.1$ | $\alpha=0.3$ | $\alpha=0.5$ | $\alpha=0.7$ | $\alpha=0.9$
---|---|---|---|---|---
$N$ | price | error | price | error | price | error | price | error | price | error
$1000$ | $0.830$ | $0.044$ | $0.834$ | $0.044$ | $0.847$ | $0.045$ | $0.822$ | $0.044$ | $0.855$ | $0.047$
$10000$ | $0.800$ | $0.015$ | $0.801$ | $0.013$ | $0.808$ | $0.014$ | $0.792$ | $0.014$ | $0.810$ | $0.014$
$20000$ | $0.7812$ | $0.0100$ | $0.7894$ | $0.0010$ | $0.7989$ | $0.0099$ | $0.8036$ | $0.0100$ | $0.8154$ | $0.0099$
$50000$ | $0.7719$ | $0.0063$ | $0.7814$ | $0.0063$ | $0.7906$ | $0.0062$ | $0.8011$ | $0.0063$ | $0.8096$ | $0.0062$
$100000$ | $0.7703$ | $0.0043$ | $0.7806$ | $0.0043$ | $0.7909$ | $0.0044$ | $0.8008$ | $0.0044$ | $0.8114$ | $0.0044$

Table 5: Forward start Asian option.

###  4.4 Swing Options

As a last application we consider the pricing of swing options that are
typical components of gas contracts, which offer the opportunity to vary the
contracted volume under a number of restrictions. We consider the day-ahead
market model illustrated in Section 4.2 and its relative estimated parameters,
however, in contrast to the previous example, we take NIG-based factors as per
the original calibration only.

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, we also assume that only
one unit of the underlying can be exercised any time period. Let $V(n,s,t)$
denote the price of 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

$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$ (42)

and $V(n,s,T)=(S(T)-K)^{+}$, $n\leq N$ and $V(0,s,t)=0$. We solve the
recursion equation with the modified version of the LSMC, introduced in
Longstaff-Schwartz [32], detailed in Hambly et al [25] where the continuation
value is approximated with a linear regression using power polynomials 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.$

We used $B=3$, but the regression may be performed on a different set of basis
functions as well.

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

Table 6 shows the values and MC errors relatively to the pricing of a
$120-120$ forward start swing option with maturity $T=1+1/3$ and strike price
$K=11.5$, namely the holder has $N=120$ rights and must exercise all of them.
We observe that the LSMC combined with Algorithm 1 for the OU-NIG factor - and
a standard procedure for the generation of the NIG process- produces unbiased
results and apparently $2\times 10^{4}$ simulations are required to attain an
acceptable convergence. Overall, it is evident that our newly developed
approach can achieve high accuracy as well as efficiency.

$N$ | $1000$ | $10000$ | $20000$ | $50000$ | $100000$
---|---|---|---|---|---
price | $451.1$ | $454.37$ | $457.90$ | $459.34$ | $459.69$
error | $15.7$ | $5.20$ | $3.69$ | $2.33$ | $1.66$
Table 6: $120-120$ swing option.

###  4.5 Application to Forward Markets

The spectrum of financial applications of the results and the simulation
procedures illustrated in Section 3 is not restricted to the modeling of day-
ahead commodity prices and to examples involving mean-reverting OU processes.
For instance, it is common practice to model the evolution $f(t,T)$ of a
forward contract at time $t$ with delivery at $T$ as a geometric BM with a
time-dependent volatility function in order captures the Samuelson effect. In
particular, Kiesel et al. [29] have considered the following two-factor market
dynamics

$\displaystyle f(t,T)$ $\displaystyle=$ $\displaystyle
f(0,T)\,\exp\left\\{h_{1}(t,T)+h_{2}(t,T)+\sigma_{1}\int_{0}^{t}e^{-b\,(T-u)}\,dW_{1}(u)+\sigma_{2}\,W_{2}(t)\right\\}$
$\displaystyle=$ $\displaystyle
f(0,T)\,\exp\left\\{h_{1}(t,T)+h_{2}(t,T)+X_{1}(t,T)+X_{2}(t)\right\\}$

where $W_{1}(\cdot)$ and $W_{2}(\cdot)$ are correlated or independent Brownian
motions and $h_{1}(t,T)$ and $h_{2}(t,T)$ are two time-dependent functions
that ensure that $f(\cdot,T)$ is a martingale. From such a model one can
obtain day-ahead price evolution letting $t=T$. On the other hand, leaving the
Gaussian framework by replacing the two Wiener processes by two general Lévy
processes $L_{1}(\cdot)$ and $L_{2}(\cdot)$, one can rewrite the first factor
as
$X_{1}(t,T)=e^{-b\,T}\int_{0}^{t}e^{-b\,(t-u)}\,dL_{1}(u)=e^{-b\,T}Z_{1}(t)$.
Of course, when $L_{1}(\cdot)$ is a SNTS process, $Z_{1}(\cdot)$ becomes the
additive process studied in Section 3, and therefore we can easily adapt all
the results obtained so far.

The contract above is a fictitious contract often used to build the stochastic
evolution of a traded contract $F(t,T_{1},T_{2})$, sometimes called swap in
energy markets, at time $t$, maturity $T$, $t\leq T\leq T_{1}<T_{2}$, and with
delivery period $[T_{1},T_{2}]$. In alternative, it is also common practice to
directly model $F(t,T_{1},T_{2})$ as

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

where $m_{1}(t,T_{1},T_{2})$ and $m_{2}(t,T_{1},T_{2})$ are two deterministic
functions that ensure that $F(\cdot,T_{1},T_{2})$ is a martingale.

Instead of an exponential setting, recently Piccirilli et al. [39] 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 NIG setting where the

$\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)=$
(43) $\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 centered NIG
processes. In this specific case

$\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).$
(44)
$\Gamma(T_{1},T_{2})=\frac{1}{T_{2}-T_{1}}\int_{T_{1}}^{T_{2}}\gamma(u)\,du$
(45)

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. [8] and Latini et al. [31], the special form of the coefficients
arises from the implicitly underlying assumption that the swap can be written
as the average over an underlying artificial futures price with instantaneous
delivery.

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)$ under the assumption that $L_{1}(\cdot)$ is symmetric. It is
worthwhile noticing that the Piccirilli et al. found an explicit form of the
chf of $X_{1}(t,T_{1},T_{2})$ when $L_{1}(\cdot)$ is a centered NIG process
that can be extended to symmetric NTS processes based on Proposition 3.2, we
omit an explicit proof to avoid overloading the paper with lengthy details.

We remark then that our simulation procedure enlarges the applicability of NOA
models and other models based on the additive process $Z_{1}(\cdot)$ driven by
NIG and in general NTS processes. Indeed, trading units rely on quoted options
as a starting point to calibrate the parameters of their models to then price
other derivative contracts or to implement trading strategies. Some of them
can be evaluated without MC methods however, simulation approaches are
employed by risk managers to derive information about the distribution of the
traded portfolios and a view on their risk profile. To this end, it is
fundamental to rely on exact and non-biased simulation schemes with
sufficiently fast computational speed that can be run with different parameter
settings and market conditions.

##  5 Concluding Remarks

In this study we have investigated the pricing of energy derivatives in
markets driven by NTS processes that generalize NIG processes which are
commonly used in the modeling of energy markets. To this end, the first
contribution of our study is the full description of NTS processes of OU type,
along with the derivation of characteristic function of their transition law
in closed form. This result is instrumental to determine their statistical
properties and compared to the current state of affairs (see for instance
Benth and Benth [9] and Benth et al. [6]) we derive non-arbitrage conditions
for markets driven by such processes without resorting to numerical
approximations or numerical integrations.

The second contribution consists in defining an efficient algorithm for the
exact simulation of the trajectories of symmetric OU-NTS and OU-NIG processes
that is particularly suitable for forward start contracts, where the standard
Euler scheme would return biased results. Indeed, even though it is common
practice to rely on such an approximated scheme, we have shown that its
performance quickly deteriorates. We also propose an alternative approximation
scheme that apparently outperforms the Euler approximation having the same
computational effort.

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 the first
example we have considered a one-factor model purely driven by a single mean-
reverting OU-NTS process. We have made use of the explicit knowledge of the
chf to implement the pricing with the FFT-based technique of Carr and Madan
[10], and then have compared the outcomes to those obtained via MC
simulations. We have noticed that such a model, although performing reasonably
well in a wide range of applications (see once again Benth and Benth [9] and
Benth et al. [6]), is not suitable for the evaluation of contracts with long
maturities which drove us to add a second factor like in the setting proposed
by de Jong and Schneider [17] in the Gaussian framework.

In the second example, we have calibrated the two factor model taking the
German day-ahead NCG prices and have evaluated a forward start Asian option.
Based on our theoretical results, we have conceived a joint calibration-
simulation strategy: under the assumption that one of the factors is a
symmetric OU-NIG process we have suggested to estimate its parameters using
the proposed approximation, whereas we have proposed to use the exact
simulation scheme for the pricing of energy derivatives especially if they are
forward start contracts. We also have studied the impact of moving from a NIG
model to a general NTS model and have found that this has a non-negligible
impact on the Asian option price.

In addition, we have shown that the proposed exact simulation algorithm,
combined with the LSMC approach of Hambly et al [25], provides an efficient
and accurate pricing of an 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. [8], Latini et al. [31] and Piccirilli
et al. [39] 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 [35],
Ballotta and Bonfiglioli [1] or the recent approaches of Gardini et al. [20,
21] and Lu [33].

We remark, that all the algorithms that we have discussed are based on the
sequential-forward generation of processes. On the other hand, as shown in
Pellegrino and Sabino [38] and Sabino [43], the backward simulation is more
suitable to the LSMC method and provides a computational advantage in the
pricing of swings and storages. Therefore, a last topic deserving further
investigation is the study of the simulation of OU-NTS processes backward in
time.

## Acknowledgements

I would like to express my gratitude to Prof. Simone Boehrer for her critical
and constructive comments which helped to improve this article.

## 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] O. Bardou, S. Bouthemy, and G. Pagés.  Optimal Quantization for the Pricing of Swing Options.  Applied Mathematical Finance, 16(2):183–217, 2009.
  * [3] 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.
  * [4] 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.
  * [5] 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.
  * [6] 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.
  * [7] 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.
  * [8] 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.
  * [9] 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.
  * [10] P. Carr and D.B. Madan.  Option Valuation Using the Fast Fourier Transform.  Journal of Computational Finance, 2:61–73, 1999.
  * [11] 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.
  * [12] R. Cont and P. Tankov.  Financial Modelling with Jump Processes.  Chapman and Hall, London, 2004.
  * [13] N. Cufaro Petroni.  Self-decomposability and Self-similarity: a Concise Primer.  Physica A, Statistical Mechanics and its Applications, 387(7-9):1875–1894, 2008.
  * [14] N. Cufaro Petroni and P. Sabino.  Tempered Stable Distribution and Finite Variation Ornstein-Uhlenbeck Processes.  Available at: https://arxiv.org/abs/2011.09147.
  * [15] M. Cummins, G. Kiely, and B. Murphy.  Gas Storage Valuation under Lévy Processes using Fast Fourier Transform.  Journal of Energy Markets, 4:43–86, 2017.
  * [16] M. Cummins, G. Kiely, and B. Murphy.  Gas Storage Valuation under Multifactor Lévy Processes.  Journal of Banking and Finance, 95:167–184, 2018.
  * [17] C. de Jong and S. Schneider.  Cointegration between Gas and Power Spot Prices.  The Journal of Energy Markets, 2(3):27–46, 2009.
  * [18] L. Devroye.  Non-Uniform Random Variate Generation.  Springer-Verlag, New York, 1986.
  * [19] L. Devroye.  Random Variate Generation for Exponential and Polynomially Tilted Stable Distributions.  ACM Transactions on Modeling and Computer Simulation, 19(4), 2009\.  Article No. 18.
  * [20] 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.
  * [21] M. Gardini, P. Sabino, and E. Sasso.  Correlating Lévy Processes with Self-decomposability: Applications to Energy Markets, 2020.  Available at www.arxiv.org.
  * [22] M. Grabchak.  Tempered Stable Distributions.  Springer International Publishing, 2016.
  * [23] I. S. Gradshteyn and I. M. Ryzhik.  Table of Integrals, Series, and Products.  Elsevier/Academic Press, Amsterdam, seventh edition, 2007.
  * [24] B. Hambly, S. Howison, and T. Kluge.  Information-Based Models for Finance and Insurance.  Quantitative Finance, 9(8):937–949, 2009.
  * [25] B. Hambly, S. Howison, and T. Kluge.  Modelling Spikes and Pricing Swing Options in Electricity Markets.  Quantitative Finance, 9(8):937–949, 2009.
  * [26] M. Hofert.  Sampling Exponentially Tilted Stable Distributions.  ACM Transactions on Modeling and Computer Simulation, 22(1), 2012\.
  * [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] 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.
  * [30] M. Kjaer.  Pricing of Swing Options in a Mean Reverting Model with Jumps.  Applied Mathematical Finance, 15(5-6):479–502, 2008.
  * [31] 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.
  * [32] 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.
  * [33] K.W. Lu.  Calibration for Multivariate Lévy-Driven Ornstein-Uhlenbeck Processes with Applications to Weak Subordination, 2020.  Available at www.arxiv.org.
  * [34] J.J. Lucia and E.S. Schwartz.  Electricity Prices and Power Derivatives: Evidence from the Nordic Power Exchange.  Review of Derivatives Research, 5(1):5–50, Jan 2002.
  * [35] 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.
  * [36] T. Meyer-Brandis and P. Tankov.  Multi-factor Jump-diffusion Models of Electricity Prices.  International Journal of Theoretical and Applied Finance, 11(5):503–528, 2008.
  * [37] J. R. Michael, W. R. Schucany, and R. W. Haas.  Generating Random Variates Using Transformations with Multiple Roots.  The American Statistician, 30(2):88–90, 1976.
  * [38] 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.
  * [39] 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.
  * [40] Jan Rosinski.  Tempering Stable Proceses.  Stochastic Processes and their Applications, 117(6):677 – 707, 2007\.
  * [41] P. Sabino.  Pricing Energy Derivatives in Markets Driven by Tempered Stable and CGMY Processes of Ornstein-Uhlenbeck Type.  available at www.arxiv.org.
  * [42] 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.
  * [43] P. Sabino.  Forward or Backward Simulations? A Comparative Study.  Quantitative Finance, 20(7):1213–1226, 2020.
  * [44] P. Sabino and N. Cufaro Petroni.  Fast Pricing of Energy Derivatives with Mean-Reverting Jump-diffusion Processes.  Applied Mathematical Finance, 0(0):1–22, 2021.
  * [45] K. Sato.  Lévy Processes and Infinitely Divisible Distributions.  Cambridge U.P., Cambridge, 1999.
  * [46] P. Schwartz and J.E. Smith.  Short-term Variations and Long-term Dynamics in Commodity Prices.  Management Science, 46(7):893–911, 2000.
  * [47] 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.
  * [48] 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.