# How to handle the COS method for option pricing

Gero Junike Carl von Ossietzky Universität, Institut für Mathematik, 26129
Oldenburg, Germany, E-mail<EMAIL_ADDRESS>
###### Abstract

The Fourier Cosine Expansion (COS) method is used to price European options
numerically in a very efficient way. To apply the COS method, one has to
specify two parameters: a truncation range for the density of the log-returns
and a number of terms $N$ to approximate the truncated density by a cosine
series. How to choose the truncation range is already known. Here, we are able
to find an explicit and useful bound for $N$ as well. We further show that the
COS method has an exponential order of convergence when the density is smooth
and decays exponentially. However, when the density is smooth and has heavy
tails, as in the Finite Moment Log Stable model, the COS method does not have
exponential order of convergence. Numerical experiments confirm the
theoretical results.

##  1 Introduction

For model calibration, it is of utmost importance to price European options
very quickly because stock price models are typically calibrated to given
prices of liquid call and put options by minimizing the mean-square-error
between model prices and given market prices. During the optimization routine,
the model prices of call and put options need to be evaluated very often for
different model parameters, see Fang and Oosterlee (2009a); Junike and
Pankrashkin (2022).

To compute the price of a European option, one must solve an integral
involving the product of the density of the log-returns at maturity and the
payoff function. However, for many financial models, the density $f$ of the
log-returns is unknown. Fortunately, the characteristic function of the log-
returns is often given in closed-form and can be used to obtain the density.

In their seminal paper, Fang and Oosterlee (2009a) proposed the COS method,
which is a very efficient way to approximate the density and to compute option
prices. The COS method requires two parameters: a truncation range for the
density of the log-returns and a number of terms $N$ to approximate the
truncated density by a cosine series. While it is known how to choose the
truncation range, see Junike and Pankrashkin (2022), the choice of $N$ is
largely based on trial and error.

The COS method has been extensively extended and applied, see Fang and
Oosterlee (2009b, 2011); Grzelak and Oosterlee (2011); Ruijter and Oosterlee
(2012); Zhang and Oosterlee (2013); Leitao et al. (2018); Liu et al. (2019a,
b); Oosterlee and Grzelak (2019); Bardgett et al. (2019). Other Fourier
pricing techniques are discussed e.g. by Carr and Madan (1999); Lord et al.
(2008); Ortiz-Gracia and Oosterlee (2013, 2016).

With respect to these papers we make two main contributions: a) we develop an
explicit, useful and rigorous bound for $N$ and b) we analyze the order of
convergence of the COS method in detail.

Fang and Oosterlee (2009a) proposed to approximate the (unknown) density in
three steps: i) Truncate the density $f$, i.e. approximate $f$ by a function
$f_{L}$ with finite support on some (sufficiently large) interval $[-L,L]$.
ii) Approximate $f_{L}$ by a Fourier cosine expansion $\sum a_{k}e_{k}(x)$,
where $a_{k}$ are Fourier coefficients of $f_{L}$ and $e_{k}$ are cosine basis
functions. iii) Approximate $a_{k}$ by some coefficients $c_{k}$ which can be
obtained directly from the characteristic function of $f$. Thus, to apply the
COS method, two decisions must to be made: find a suitable truncation range
$[-L,L]$ and identify the number $N$ of cosine functions.

Junike and Pankrashkin (2022) applied a simple triangle equation to bound the
error of the three approximations and obtained:

$\bigg{\|}f-\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N}c_{k}e_{k}\bigg{\|}_{2}\leq\big{\|}f-f_{L}\big{\|}_{2}+\Big{\|}f_{L}-\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N}a_{k}e_{k}\Big{\|}_{2}+\Big{\|}\sideset{}{{}^{\prime}}{\sum}_{k=0}^{\infty}(a_{k}-c_{k})e_{k}\Big{\|}_{2}.$
(1)

The first, second and third terms at the right-hand-side of Inequality (1)
correspond to approximations due to i), ii), and iii), respectively.

Junike and Pankrashkin (2022) developed explicit expressions for bounds for
the first and third terms on the right-hand-side of Inequality (1), which
depend only on the choice of $L$, not on the choice of $N$. When $f$ has semi-
heavy tails, i.e. the tails decay exponentially, Junike and Pankrashkin (2022)
applied Markov’s inequality to obtain an explicit and useful expression for
$L$.

Recently, Aimi et al. (2023) proposed some hints on how to choose $N$. But, to
the best of our knowledge, a rigorous bound on $N$ to ensure a certain
precision of the COS method is not yet known. On the other hand, it is well
known that the series truncation error, i.e. the second term on the right-
hand-side of Inequality (1), can be bounded using integration by parts, see
Boyd (2001). In this article, we use this idea to find an explicit and useful
bound for $N$, which involves the bound of the $j^{th}$-derivative of the
density $f$.

For some models such as the Black-Scholes model, see Black and Scholes (1973),
the Normal Inverse Gaussian (NIG) model, see Barndorff-Nielsen (1997), or the
Finite Moment Log Stable (FMLS) model, see Carr and Wu (2003), we find sharp
bounds for all derivatives of the density of $f$ in closed-form. In general,
those bounds can also be estimated numerically. We obtain an explicit and
useful bound for $N$, which is provably large enough to ensure that the COS
method converges within a predefined error tolerance.

Fang and Oosterlee (2009a) also analyzed the order of convergence of the COS
method focusing on the second term on the right-hand-side of Inequality (1).
They concluded that with a properly chosen truncation range, the overall error
converges exponentially for smooth density functions and compares favorably to
the Carr-Madan formula, see Carr and Madan (1999).

In this article, we also consider the errors introduced by the truncation
range, i.e., errors due to i), ii), and iii) and establish upper bounds for
the order of convergence of the COS method. We confirm, both theoretically and
empirically, that the COS method indeed converges exponentially for smooth
density functions if in addition, the tails of the density decay at least
exponentially.

However, for fat-tailed and smooth densities such as the Cauchy density or the
density of the log-returns in the FMLS model, the truncation error due to i)
and iii) becomes much more relevant compared to densities with semi-heavy
tails. We show theoretically that the COS method converges at least as fast as
$O(N^{-\alpha})$ for $N\to\infty$, where $\alpha>0$ is the Pareto tail index,
e.g. for the Cauchy density, $\alpha=1$ and for the FMLS model
$\alpha\in(1,2)$. Empirical experiments indicate that the COS method converges
as fast as $O(N^{-\alpha})$, i.e. the theoretical bound is sharp.

This article is structured as follows: Section 2 gives an overview of the
technical details of the COS method. Section 3 gives explicit formulas for the
truncation range $L$ and the number of terms $N$. Section 4 analyzes the order
of convergence of the COS method. In Sections 3 and 4, we distinguish between
models with semi-heavy tails and models with heavy tails. Section 5 contains
numerical experiments which confirm the theoretical results. Section 6
concludes.

##  2 Overview: the COS method for option pricing.

We model the stock price over time by a stochastic process $(S_{t})_{t\geq 0}$
on a filtered probability space. We assume that there is a bank account paying
continuous compound interest $r\in\mathbb{R}$. The market is assumed to be
arbitrage free and we denote by $Q$ the risk-neutral measure. All expectations
are taken under $Q$.

There is a European option with maturity $T>0$ and payoff $g(S_{T})$ at $T$,
where $g:[0,\infty)\to\mathbb{R}$. For example, a European put option with
strike $K>0$ can be described by the payoff $g(x)=\max(K-x,0)$, $x\geq 0$.

In several places, we assume that the payoff function is bounded. The prices
of European call options are not bounded, but can be obtained from the prices
of put options by the put-call parity.

Provided that the characteristic function $\varphi_{\log(S_{T})}$ of
$\log(S_{T})$ is given in closed-form, the COS method is able to price the
European option numerically very quickly as follows: we denote by

$X:=\log(S_{T})-E[\log(S_{T})],$

the _centralized log-returns_. The expectation can be computed by

$E[\log(S_{T})]=-i\varphi_{\log(S_{T})}^{\prime}(0).$

The characteristic function $\varphi$ of $X$ is then given by

$\varphi(u)=\varphi_{\log(S_{T})}(u)\exp(-iuE[\log(S_{T})]),\quad
u\in\mathbb{R}.$

We assume that $X$ has a density $f$, but the exact structure of $f$ need not
be known. Since the density of $X$ is centered around zero, it is justified to
truncate the density $f$ on a symmetric interval $[-L,L]$, see Junike and
Pankrashkin (2022). Define

$v(x):=e^{-rT}g(\exp(x+E[\log(S_{T})]),\quad x\in\mathbb{R}.$

The price of the European option is then given by

$\displaystyle e^{-rT}E[g(S_{T})]$
$\displaystyle=\int_{\mathbb{R}}v(x)f(x)dx.$ (2)

We need some abbreviations to discuss the COS method for solving the integral
at the right-hand-side of Equation (2) numerically: for $L>0$ let
$f_{L}:=1_{[-L,L]}f$. Define the basis functions

$e_{k}(x)=\cos\left(k\pi\frac{x+L}{2L}\right),\quad k=0,1,...$

and the Fourier coefficients of $f_{L}$ are given by $a_{k}$ and approximated
by $c_{k}$, where

$a_{k}:=\frac{1}{L}\int_{-L}^{L}f(x)e_{k}(x)dx\quad\text{and}\quad
c_{k}:=\frac{1}{L}\int_{\mathbb{R}}f(x)e_{k}(x)dx,\quad k=0,1,...$

Observe

$c_{k}=\frac{1}{L}\Re\bigg{\\{}\varphi\bigg{(}\frac{k\pi}{2L}\bigg{)}e^{i\frac{k\pi}{2}}\bigg{\\}},\quad
k=0,1,...,$

i.e. the coefficients $c_{k}$ can be obtained explicitly if $\varphi$ is given
in closed-form. Here, $\Re(z)$ denotes the real part of a complex number $z$.
For $0<M\leq L$ define

$v_{k}:=\int_{-M}^{M}v(x)e_{k}(x)dx,\quad k=0,1,...$ (3)

To keep the notation simple, we suppress the dependence of $e_{k}$ and $c_{k}$
on $L$ and the dependence of $v_{k}$ on $M$. The COS method states that the
price of the European option can be approximated by

$\int_{\mathbb{R}}v(x)f(x)dx\approx\int_{-M}^{M}v(x)\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N}c_{k}e_{k}(x)dx=\sum_{k=0}^{N}{}^{\prime}c_{k}v_{k},$
(4)

where $\sum{}^{\prime}$ indicates that the first summand (with $k=0$) is
weighted by one-half. The coefficients $c_{k}$ are given in closed-form when
$\varphi$ is given analytically and the coefficients $v_{k}$ can also be
computed explicitly in important cases, e.g. for plain vanilla European put or
call options and digital options, see for instance appendix in Junike and
Pankrashkin (2022). This makes the COS method numerically very efficient and
robust.

To give the approximation in line (4) a precise meaning, we need a bound for
the term

$B(L):=\sum_{k=0}^{\infty}\frac{1}{L}\left|\int_{\mathbb{R}\setminus[-L,L]}f(x)\cos\left(k\pi\frac{x+L}{2L}\right)dx\right|^{2},\quad
L>0.$ (5)

Junike and Pankrashkin (2022) called integrable functions $f$ with $B(L)\to
0$, $L\to\infty$ _COS-admissible._ The term $\sqrt{B(L)}$ is equal to third
term at the right-hand-side of Inequality (1).

The class of COS-admissible densities is very large, in particular it includes
bounded densities with existing first and second moments. The next Lemma 2.1
is taken from Junike and Pankrashkin (2022). Theorem 3.6 and Theorem 3.10
apply Lemma 2.1 and provide explicit formulas for the parameters $M$, $L$ and
$N$ of the COS method for densities with semi-heavy tails and heavy tails,
respectively.

######  Lemma 2.1.

Let $f$ be integrable and COS-admissible. Let $v:\mathbb{R}\to\mathbb{R}$ be
bounded, with $|v(x)|\leq K$ for all $x\in\mathbb{R}$ and some $K>0$. Let
$\varepsilon>0$. Let $M>0$ so that

$\int_{\mathbb{R}\setminus[-M,M]}v(x)f(x)dx\leq\frac{\varepsilon}{2}.$ (6)

Define $\xi=\sqrt{2M}K$. Let $L\geq M$ so that

$\left\|f-f_{L}\right\|_{2}\leq\frac{\varepsilon}{6\xi},\quad\text{and}\quad\sqrt{B(L)}\leq\frac{\varepsilon}{6\xi}.$
(7)

Choose $N_{L}$ large enough so that

$\left\|f_{L}-\sum_{k=0}^{N}{}^{\prime}a_{k}e_{k}\right\|_{2}\leq\frac{\varepsilon}{6\xi}.$
(8)

Then it holds for all $N\geq N_{L}$ that

$\left|\int_{\mathbb{R}}v(x)f(x)dx-\sum_{k=0}^{N}{}^{\prime}c_{k}v_{k}\right|\leq\varepsilon.$

###### Proof.

Junike and Pankrashkin (2022, Cor. 8). ∎

Often, it is fine to choose $M=L$, e.g. when applying the COS method to
densities with semi-heavy tails. However, if the density has heavy tails, it
is usually numerically more efficient to choose $L$ and $M$ differently.

For some models, e.g. the FMLS model, the exact tail behavior of the density
of the log-returns is known and the inequalities in (6) and (7) can be solved
explicitly for $M$ and $L$. For most other models, e.g. the Heston model, we
only know that the tails of the density of the log-returns decay
exponentially, but we do not know the exact behavior of the tails. For such
densities, Junike and Pankrashkin (2022) proposed to use Markov’s inequality
to obtain explicit bounds for $M$ and $L$.

The main requirement for applying the COS method is the availability of the
characteristic function in closed-form. Many models satisfy this requirement
and we list some examples:

BS

the Black-Scholes model with volatility $\sigma$, see Black and Scholes
(1973),

Heston

the Heston model with speed of mean-reversion $\kappa>0$, level of mean-
reversion $\eta>0$, vol of var $\theta>0$, initial vol $v_{0}>0$ and
correlation $\rho\in[-1,1]$, see Heston (1993),

VG

the variance gamma model with parameters $\sigma>0$, $\theta\in\mathbb{R}$ and
$\nu>0$, see Madan et al. (1998),

NIG

the normal inverse Gaussian model with parameters $\alpha>0$,
$\beta\in(-\alpha,\alpha)$ and $\delta>0$, see Barndorff-Nielsen (1997) and
Schoutens (2003, Sec. 5.3.8).

CGMY

the CGMY model with parameters $C>0$, $G>0$, $M>0$, $Y<2$, see Carr et al.
(2002),

FMLS

the Finite Moment Log Stable model with parameters $\sigma>0$ and
$\alpha\in(1,2)$, see Carr and Wu (2003).

To use the COS method, the following parameters must be specified: $M>0$ and
$L>0$ to define the truncation range and $N\in\mathbb{N}$ to set the number of
terms. $M$ and $L$ can be estimated by making assumptions on the tail-behavior
of the density. $N$ depends on the smoothness of the density.

The densities of the models listed above have the following properties
regarding tail-behavior and smoothness: The density of the log-returns in the
BS, NIG, Heston models and the CGMY model with parameter $Y\in(0,1)$, is
infinitely many times differentiable and has semi-heavy tails, i.e. the tails
decay exponentially, see Küchler and Tappe (2013); Asmussen (2022); Albin and
Sundén (2009); Schoutens (2003); Dragulescu and Yakovenko (2002).

The density of the VG model at time $T>0$ has semi-heavy tails, see Albin and
Sundén (2009, Example 7.5) and is $(J+1)$-times continuously differentiable if
$J+2<\frac{2T}{\nu}$, see Küchler and Tappe (2008, Thm. 4.1).

The stable law has been used to model stock prices since Mandelbrot (1963) and
Fama (1965). The density of the FMLS model also belongs to the family of
stable densities with maximum negative skewness. It is infinitely many times
differentiable and has a heavy left tail with Pareto index $\alpha$, i.e. the
left tail decays like $|x|^{-1-\alpha}$, $x\to-\infty$, but the right tail
decays exponentially, see Carr and Wu (2003). This makes the FMLS very
attractive from a theoretical point of view: put and call option prices and
all moments of the underlying stock price $S_{T}$ exist. The expectation of
the log-returns also exists, but the log-returns do not have finite variance.
Fitting the FMLS model to real market data shows very satisfactory results,
see Carr and Wu (2003).

##  3 How to find $M$, $L$ and $N$

We summarize the assumptions about the density $f$ of the log-returns in order
to find explicit expressions for $L$ and $N$.

Let $\mathbb{N}_{0}=\\{0\\}\cup\mathbb{N}$ and $J\in\mathbb{N}_{0}$. We denote
by $C_{b}^{J+1}(\mathbb{R})$ the $(J+1)$-times, continuously differentiable
functions from $\mathbb{R}$ to $\mathbb{R}$ which are bounded with bounded
derivatives.

By $f^{(j)}$ we denote the $j^{th}$-derivative of $f$. We use the convention
$f^{(0)}\equiv f$. We denote by $H_{j}>0$ a bound for $f^{(j)}$, i.e.
$|f^{(j)}(x)|\leq H_{j}$, $x\in\mathbb{R}$.

Let $C_{1},C_{2},C_{3}>0$ be suitable constants. Let $L_{0}>0$ be large
enough. Suppose that $f$ satisfies one of the following assumptions:

Assumption A:

$f\in C_{b}^{J+1}(\mathbb{R})$ is integrable and $f$ has semi-heavy tails,
i.e.

$|f^{(j)}(x)|\leq C_{1}e^{-C_{2}|x|},\quad j=0,...,J+1,\quad|x|\geq L_{0}.$
(9)

Assumption B:

$f\in C_{b}^{J+1}(\mathbb{R})$ is integrable and $f$ has heavy tails with
index $\alpha>0$, i.e.

$|f^{(j)}(x)|\leq C_{3}|x|^{-1-\alpha-j},\quad j=0,...,J+1,\quad|x|\geq
L_{0}.$ (10)

We use the following abbreviations: the maximum of two real numbers $x,y$, is
denoted by $x\lor y$. By $\Gamma$ we denote the gamma function

$\thinspace\Gamma(z)=\int_{0}^{\infty}u^{z-1}e^{-u}du,\quad\Re(z)>0.$

The following lemma allows to bound the series-truncation error, which depends
only on the choice of $N$. It is known in a similar form in the literature,
see e.g. Theorem 1.39 in Plonka et al. (2018), Theorem 4.2 in Wright et al.
(2015) or Theorem 6 in Boyd (2001) and can be proved by integration by parts.
It is usually stated for functions with range $[-1,1]$ or $[0,2\pi]$. Here, we
explicitly need the dependence of the series-truncation error on the
truncation range $[-L,L]$, so we give the full proof.

######  Lemma 3.1.

Let $J\in\mathbb{N}_{0}$. If $f\in C_{b}^{J+1}(\mathbb{R})$ it holds for
$J\geq 1$

$\left\|f_{L}-\sum_{k=0}^{N}{}^{\prime}a_{k}e_{k}\right\|_{2}\leq\sum_{j=1}^{J}\frac{2^{j+1}}{j\pi^{j+1}}\frac{L^{j}}{N^{j}}\left(|f^{(j)}(-L)|+|f^{(j)}(L)|\right)+\frac{2^{J+2}H_{J+1}}{J\pi^{J+1}}\frac{L^{J+1}}{N^{J}}$
(11)

and for $J=0$

$\left\|f_{L}-\sum_{k=0}^{N}{}^{\prime}a_{k}e_{k}\right\|_{2}\leq\frac{4H_{1}}{\pi}\frac{L}{\sqrt{N}}.$

###### Proof.

It holds for any $\nu>0$ by integration by parts, see Lyness (1971, Eq. 1.3),

$\displaystyle\int_{-L}^{L}f(x)e^{i\nu x}dx=$
$\displaystyle\sum_{j=0}^{J}\frac{i^{j+1}}{\nu^{j+1}}\left(e^{-i\nu
L}f^{(j)}(-L)-e^{i\nu L}f^{(j)}(L)\right)$
$\displaystyle+\frac{i^{J+1}}{\nu^{J+1}}\int_{-L}^{L}f^{(J+1)}(x)e^{i\nu
x}dx.$ (12)

For $k\in\mathbb{N}$, we apply Equation (12) for $\nu:=\frac{k\pi}{2L}$. Then
it follows that

$\displaystyle|a_{k}|=$
$\displaystyle\frac{1}{L}\left|\int_{-L}^{L}f(x)\cos\left(k\pi\frac{x+L}{2L}\right)dx\right|$
$\displaystyle=$
$\displaystyle\frac{1}{L}\left|\Re\left\\{e^{i\frac{k\pi}{2}}\int_{-L}^{L}f(x)e^{i\frac{k\pi}{2L}x}dx\right\\}\right|$
$\displaystyle=$
$\displaystyle\frac{1}{L}\bigg{|}\Re\bigg{\\{}\sum_{j=0}^{J}i^{j+1}\frac{(2L)^{j+1}}{(k\pi)^{j+1}}\left(f^{(j)}(-L)-(-1)^{k}f^{(j)}(L)\right)$
$\displaystyle+e^{i\frac{k\pi}{2}}i^{J+1}\frac{(2L)^{J+1}}{(k\pi)^{J+1}}\int_{-L}^{L}f^{(J+1)}(x)e^{i\frac{k\pi}{2L}x}dx\bigg{\\}}\bigg{|}$
$\displaystyle\leq$
$\displaystyle\sum_{j=1}^{J}\frac{2^{j+1}}{\pi^{j+1}}\frac{L^{j}}{k^{j+1}}\left(|f^{(j)}(-L)|+|f^{(j)}(L)|\right)+\frac{2^{J+2}H_{J+1}}{\pi^{J+1}}\frac{L^{J+1}}{k^{J+1}}.$

Observe:

$\left\|f_{L}-\sum_{k=0}^{N}{}^{\prime}a_{k}e_{k}\right\|_{2}=\sqrt{\sum_{k=N+1}^{\infty}|a_{k}|^{2}}\leq\sum_{k=N+1}^{\infty}|a_{k}|.$
(13)

By the integral test for convergence

$\sum_{k=N+1}^{\infty}\frac{1}{k^{j+1}}\leq\int_{N}^{\infty}x^{-j-1}dx=\frac{1}{j}N^{-j},\quad
j\geq 1,$ (14)

which implies Eq. (11) for $J\geq 1$. If $J=0$, apply the first Equality from
(13) and the Inequality (14). ∎

Given $L>0$, we need an estimate for $H_{J+1}$ to find an upper bound for the
series truncation error by Inequality (11). $L$ can be found by applying
Markov’s inequality, see Junike and Pankrashkin (2022) or Theorem 3.6 below.
Next, we provide ideas and examples for estimating $H_{J+1}$. If $f\in
C_{b}^{J+1}(\mathbb{R})$ is a density with characteristic function $\varphi$,
then it holds that

$H_{j}:=\sup_{x\in\mathbb{R}}|f^{(j)}(x)|\leq\frac{1}{2\pi}\int_{\mathbb{R}}|u|^{j}|\varphi(u)|du,\quad
j=0,1,...,J+1.$ (15)

The Inequality (15) provides explicit expressions for the term $H_{j}$,
$j=0,..,J+1$ for several models as the next examples show. The right-hand-side
of Inequality (15) can also be solved numerically, this is done in Section
5.1.

######  Example 3.2.

The density $f_{\text{Stable}}\in C_{b}^{\infty}(\mathbb{R})$ of the family of
stable distributions with stability parameter $\alpha\in(0,2]$, skewness
$\beta\in[-1,1]$, scale $c>0$ and location $\theta\in\mathbb{R}$ has
characteristic function

$\displaystyle u\mapsto$
$\displaystyle\exp\left(iu\theta-|uc|{}^{\alpha}(1-i\beta\text{sgn}(u)\Phi_{\alpha}(u)\right),$
$\displaystyle\Phi_{\alpha}(u)=\begin{cases}\tan\frac{\pi\alpha}{2}&,\alpha\neq
1\\\ -\frac{2}{\pi}\log(|u|)&,\alpha=1.\end{cases}$

It follows by Inequality (15) that

$\displaystyle\sup_{x\in\mathbb{R}}|f_{\text{Stable}}^{(j)}(x)|$
$\displaystyle\leq\frac{1}{\pi\alpha
c^{j+1}}\int_{0}^{\infty}t^{\frac{j+1}{\alpha}-1}e^{-t}dt$
$\displaystyle\leq\frac{\Gamma\left(\frac{j+1}{\alpha}\right)}{\pi\alpha
c^{j+1}},\quad j=0,1,2,...$ (16)

######  Example 3.3.

In the BS model with volatility $\sigma>0$, the centralized log-returns at
time $T>0$ have Gaussian density $f_{\text{BS}}$ with mean zero and variance
$\sigma^{2}T$. Since the normal distribution is a member of the stable
distributions, we obtain with the help of Equation (16) that

$\displaystyle\sup_{x\in\mathbb{R}}|f_{\text{BS}}^{(j)}(x)|$
$\displaystyle\leq\begin{cases}j!\left(\sigma^{j+1}T^{\frac{j+1}{2}}\sqrt{2\pi}2^{\frac{j}{2}}\left(\frac{j}{2}\right)!\right)^{-1}&,j=0,2,4,6...\\\
2^{\frac{j-1}{2}}\left(\frac{j-1}{2}\right)!\left(\sigma^{j+1}T^{\frac{j+1}{2}}\pi\right)^{-1}&,j=1,3,5,...\end{cases}$

######  Example 3.4.

In the FMLS model with parameters $\alpha\in(1,2)$ and $\sigma>0$, the
centralized log-returns at time $T>0$ are stably distributed with stability
parameter $\alpha$, skewness $\beta=-1$, scale $c=\sigma
T^{\frac{1}{\alpha}}$, and location $\theta=0$, see Carr et al. (2002). An
upper bound for the density $f_{\text{FMLS}}$ of the log-returns and its
derivatives can be obtained directly from Inequality (16).

######  Example 3.5.

In the symmetric NIG model with parameters $\alpha>0$, $\beta=0$ and
$\delta>0$, the centralized log-returns at time $T>0$ have density
$f_{\text{NIG}}\in C_{b}^{\infty}(\mathbb{R})$, whose characteristic function
is given by

$u\mapsto\exp\left(-\delta T\sqrt{\alpha^{2}+u^{2}}+\delta
T\alpha\right),\quad u\in\mathbb{R},\quad\alpha>0,\quad\delta>0,$

see Barndorff-Nielsen (1997) and Schoutens (2003, Sec. 5.3.8). We obtain by
Equation (15) and using $\alpha^{2}+u^{2}\geq u^{2}$ that

$\displaystyle\sup_{x\in\mathbb{R}}|f_{\text{NIG}}^{(j)}(x)|$
$\displaystyle\leq\frac{e^{T\delta\alpha}j!}{\pi(T\delta)^{j+1}},\quad
j=0,1,2,...$

In Theorem 3.6 we provide explicit formulas for $N$ and $L$ when $f$ satisfies
Assumption A, to ensure that the COS method approximates the true price within
a predefined error tolerance $\varepsilon>0$. Theorem 3.6 extends Junike and
Pankrashkin (2022, Cor. 9, Rem 11).

######  Theorem 3.6.

($M$, $L$ and $N$ for densities with semi-heavy tails). Let
$v:\mathbb{R}\to\mathbb{R}$ be bounded, with $|v(x)|\leq K$ for all
$x\in\mathbb{R}$ and some $K>0$. Let $\varepsilon>0$ be small enough. Suppose
$J\in\mathbb{N}_{0}$. Let $f$ be a density satisfying Assumption A. For some
even $n\in\mathbb{N}$ let

$L=M=\sqrt[n]{\frac{2K\mu_{n}}{\varepsilon}},$ (17)

where $\mu_{n}$ is the $n^{th}-$moment of $f$, i.e.

$\mu_{n}=\frac{1}{i^{n}}\left.\frac{\partial^{n}}{\partial
u^{n}}\varphi(u)\right|_{u=0}.$ (18)

Let $\xi=\sqrt{2M}K$. If $J\geq 1$ let $j\in\\{1,...,J\\}$ and

$\displaystyle N$
$\displaystyle\geq\frac{4L}{\pi}\lor\left(\frac{2^{j+2}H_{j+1}L^{j+1}}{j\pi^{j+1}}\frac{12\xi}{\varepsilon}\right)^{\frac{1}{j}}.$
(19)

If $J=0$ let

$N\geq\left(\frac{4H_{1}L}{\pi}\frac{6\xi}{\varepsilon}\right)^{2}.$ (20)

In both cases it holds that

$\left|\int_{\mathbb{R}}v(x)f(x)dx-\sum_{k=0}^{N}{}^{\prime}c_{k}v_{k}\right|\leq\varepsilon.$
(21)

###### Proof.

The $n^{th}-$moment of $f$ can be obtained from the $n^{th}$ derivative of the
characteristic function of $f$ by Equation (18). It holds by Markov’s
inequality and the definition of $M$ that

$\int_{\mathbb{R}\setminus[-M,M]}\big{|}v(x)f(x)\big{|}dx\leq
K\int_{\mathbb{R}\setminus[-M,M]}f(x)dx\leq
K\frac{\mu_{n}}{M^{n}}\leq\frac{\varepsilon}{2}.$ (22)

For $\varepsilon$ small enough, $M$ is large enough so that Assumption A
holds, i.e. $L_{0}\leq M\leq L$. Further,

$\displaystyle\left\|f-f_{L}\right\|_{2}$
$\displaystyle\leq\sqrt{2C_{1}^{2}\int_{L}^{\infty}e^{-2C_{2}x}dx}$
$\displaystyle=\frac{C_{1}}{\sqrt{C_{2}}}e^{-C_{2}L}$ (23)
$\displaystyle\leq\frac{\varepsilon}{6\xi},$

the last inequality holds true if

$L\geq-\frac{1}{C_{2}}\log\left(\frac{\sqrt{C_{2}}}{C_{1}}\frac{\varepsilon}{6\xi}\right).$
(24)

Proposition 2 in Junike and Pankrashkin (2022) shows that

$B(L)\leq\frac{2}{3}\frac{\pi^{2}}{L^{2}}\int_{\mathbb{R}\setminus[-L,L]}|xf(x)|^{2}dx.$

We then have

$\displaystyle\sqrt{B(L)}\leq$
$\displaystyle\sqrt{\frac{4}{3}\frac{\pi^{2}C_{1}^{2}}{L^{2}}\int_{L}^{\infty}x^{2}e^{-2C_{2}x}dx}$
$\displaystyle=$ $\displaystyle\frac{2\pi
C_{1}}{\sqrt{6C_{2}}}e^{-C_{2}L}\sqrt{1+\frac{1}{LC_{2}}+\frac{1}{2L^{2}C_{2}^{2}}}$
(25) $\displaystyle\leq$ $\displaystyle\frac{2\pi
C_{1}}{\sqrt{6C_{2}}}\sqrt{1+\frac{1}{MC_{2}}+\frac{1}{2M^{2}C_{2}^{2}}}e^{-C_{2}L}$
$\displaystyle\leq$ $\displaystyle\frac{\varepsilon}{6\xi},$

the last inequality holds true if

$L\geq-\frac{1}{C_{2}}\log\left(\left(\frac{2\pi
C_{1}}{\sqrt{6C_{2}}}\sqrt{1+\frac{1}{MC_{2}}+\frac{1}{2M^{2}C_{2}^{2}}}\right)^{-1}\frac{\varepsilon}{6\xi}\right).$
(26)

Assume $J\geq 1$. By Equation (19) we have $\frac{\pi}{4}\geq\frac{L}{N}$ and
it follows that

$\displaystyle\sum_{j=1}^{J}\frac{2^{j+1}}{j\pi^{j+1}}\frac{L^{j}}{N^{j}}\left(|f^{(j)}(-L)|+|f^{(j)}(L)|\right)$
$\displaystyle\leq$
$\displaystyle\frac{8}{\pi^{2}}C_{1}e^{-C_{2}L}\frac{L}{N}\sum_{k=0}^{\infty}\left(\frac{2L}{\pi
N}\right)^{k}$ $\displaystyle\leq$
$\displaystyle\frac{8}{\pi^{2}}C_{1}e^{-C_{2}L}\frac{\pi}{4}\frac{1}{1-\frac{2L}{\pi
N}}$ $\displaystyle\leq$
$\displaystyle\frac{4C_{1}}{\pi}e^{-C_{2}L}\leq\frac{\varepsilon}{12\xi},$

the last inequality holds true if

$L\geq-\frac{1}{C_{2}}\log\left(\frac{\pi}{4C_{1}}\frac{\varepsilon}{12\xi}\right).$
(27)

By Equation (17), $L$ is of order $\varepsilon^{-\frac{1}{n}}$. Hence, for
$\varepsilon$ small enough, Inequalities (24), (26) and (27) are indeed
satisfied because the right-hand-sides of those Inequalities are of order
$\log(\varepsilon)$. By the definition of $N$ in Equation (19), we also have

$\frac{2^{j+2}H_{j+1}}{j\pi^{j+1}}\frac{L^{j+1}}{N^{j}}\leq\frac{\varepsilon}{12\xi}.$

By Lemma 3.1 it follows that that

$\left\|f_{L}-\sum_{k=0}^{N}{}^{\prime}a_{k}e_{k}\right\|_{2}\leq\frac{\varepsilon}{6\xi}.$
(28)

As $B(L)\to 0$, $L\to\infty$, $f$ is COS-admissible. Apply Lemma 2.1 arrive at
Equation (21). Now assume $J=0$. By the definition of $N$ in Equation (20),
Lemma 3.1 implies again Inequality (28). Apply as before Lemma 2.1 to
conclude. ∎

######  Remark 3.7.

If $v$ is a European put option, $K$ can be set to the strike of the option.
$\varepsilon$ denotes some error tolerance. Numerical experiments in Junike
and Pankrashkin (2022) suggest to choose $n\in\\{4,6,8\\}$ for $\mu_{n}$.
According to Theorem 3.6, any $j\in\\{1,...,J\\}$ is allowed to define $N$ by
Equation (19). In the applications, one could minimize $N$ over all admissible
$j$. But this could be time-consuming and in Sections (5.1) and (5.3) we set
$j$ to a fixed value. The bounds $H_{j}$ are explicitly known for some models,
see Examples 3.5, 3.3 and 3.4 and can also be obtained numerically, see
Section 5.1. Section 5 contain examples indicating that the bound in Equation
(19) is often sharp and very useful in applications.

To extend Theorem 3.6 to the case that $f$ satisfies Assumption B, i.e. $f$
has heavy tails, we need a bound for the term $B(L)$, see Equation (5).
Proposition 2 in Junike and Pankrashkin (2022) provides a bound for $B(L)$,
but we cannot apply it to all densities of the stable distributions because
the assumptions in that Proposition are not satisfied for all stability
parameters $\alpha\in(0,2]$ of the stable law. Therefore, we need the
following lemma:

######  Lemma 3.8.

Let $f\in C_{b}^{1}(\mathbb{R})$ be integrable and square-integrable such that

$xf^{2}(x)\to 0,\quad x\to\pm\infty.$

Let $L>0$. If

$f^{\prime}(x)\geq 0,\quad x\leq-L\quad\text{and}\quad f^{\prime}(x)\leq
0,\quad x\geq L$ (29)

then

$B(L)\leq\frac{1}{L}\left(\int_{\mathbb{R}\setminus[-L,L]}f(x)dx\right)^{2}+\frac{8}{3}L\left(f^{2}(L)\lor
f^{2}(-L)\right).$

###### Proof.

It holds $f(x)\to 0$, $x\to\pm\infty$ and

$\int_{L}^{\infty}\left|f^{\prime}(x)\right|dx=-\int_{L}^{\infty}f^{\prime}(x)dx=f(L).$

Note that for all $k\in\mathbb{N}$,

$\displaystyle\int_{L}^{\infty}f(x)\cos\left(k\pi\frac{x+L}{2L}\right)dx=$
$\displaystyle\underbrace{\left[f(x)\sin\left(k\pi\frac{x+L}{2L}\right)\frac{2L}{k\pi}\right]_{L}^{\infty}}_{=0}$
$\displaystyle-\frac{2L}{k\pi}\int_{L}^{\infty}f^{\prime}(x)\sin\left(k\pi\frac{x+L}{2L}\right)dx,$

implying

$\displaystyle\frac{1}{L}\left|\int_{L}^{\infty}f(x)\cos\left(k\pi\frac{x+L}{2L}\right)dx\right|^{2}=$
$\displaystyle\frac{4}{\pi^{2}k^{2}}L\left|\int_{L}^{\infty}f^{\prime}(x)\sin\left(k\pi\frac{x+L}{2L}\right)dx\right|^{2}$
$\displaystyle\leq$ $\displaystyle\frac{4}{\pi^{2}k^{2}}Lf^{2}(L).$

Similarly,

$\frac{1}{L}\left|\int_{-\infty}^{-L}f(x)\cos\left(k\pi\frac{x+L}{2L}\right)dx\right|^{2}\leq\frac{4}{\pi^{2}k^{2}}Lf^{2}(-L).$

Using

$\sum_{k=1}^{\infty}\frac{1}{k^{2}}=\frac{\pi^{2}}{6}\quad\text{and}\quad|a+b|^{2}\leq
4(a^{2}\lor b^{2}),$

we arrive at

$\displaystyle\sum_{k=0}^{\infty}$
$\displaystyle\frac{1}{L}\left|\int_{\mathbb{R}\setminus[-L,L]}f(x)\cos\left(k\pi\frac{x+L}{2L}\right)dx\right|^{2}$
$\displaystyle\leq\frac{1}{L}\left|\int_{\mathbb{R}\setminus[-L,L]}f(x)dx\right|^{2}+4\left(\frac{4\pi^{2}}{6\pi^{2}}Lf^{2}(L)\lor\frac{4\pi^{2}}{6\pi^{2}}Lf^{2}(-L)\right)$
$\displaystyle=\frac{1}{L}\left|\int_{\mathbb{R}\setminus[-L,L]}f(x)dx\right|^{2}+\frac{8}{3}L\left(f^{2}(L)\lor
f^{2}(-L)\right).$

∎

######  Example 3.9.

Let $f$ be a unimodal density satisfying Assumption B. Then $f$ also satisfies
the assumption of Lemma 3.8. To see this, note that $f$ is bounded and
therefore square-integrable. The unimodality implies Equation (29). Further,

$|xf^{2}(x)|\leq C_{3}^{2}|x|^{-1-2\alpha}\to 0,\quad x\to\pm\infty.$

In Theorem 3.10 we provide explicit formulas for $M$, $N$ and $L$ if $f$
satisfies Assumption B and is unimodal, to ensure that the COS method
approximates the true price within a predefined error tolerance
$\varepsilon>0$. In particular, the density of the stable law is unimodal, see
Yamazato (1978).

######  Theorem 3.10.

($M$, $L$ and $N$ for densities with heavy tails). Let
$v:\mathbb{R}\to\mathbb{R}$ be bounded, with $|v(x)|\leq K$ for all
$x\in\mathbb{R}$ and some $K>0$. Let $\varepsilon>0$ be small enough. Suppose
$J\in\mathbb{N}_{0}$. Let $f$ be a unimodal density satisfying Assumption B.
Let

$\displaystyle M$
$\displaystyle=\left(\frac{4C_{3}K}{\varepsilon\alpha}\right)^{\frac{1}{\alpha}},$
(30)

$\xi:=\sqrt{2M}K$ and

$L=M\lor\left(12C_{3}\sqrt{\frac{1}{\alpha^{2}}+\frac{2}{3}}\frac{\xi}{\varepsilon}\right)^{\frac{2}{1+2\alpha}}.$
(31)

Let $j\in\\{1,...,J\\}$ and

$\displaystyle
N\geq\frac{4L}{\pi}\lor\left(\frac{2^{j+2}H_{j+1}L^{j+1}}{j\pi^{j+1}}\frac{12\xi}{\varepsilon}\right)^{\frac{1}{j}}.$
(32)

Then

$\left|\int_{\mathbb{R}}v(x)f(x)dx-\sum_{k=0}^{N}{}^{\prime}c_{k}v_{k}\right|\leq\varepsilon.$
(33)

###### Proof.

$M$ is of order $\varepsilon^{-\frac{1}{\alpha}}$. For $\varepsilon$ small
enough, $M$ is large enough so that Assumption B holds, i.e. $L_{0}\leq M\leq
L$. The inequality

$\int_{\mathbb{R}\setminus[-M,M]}\big{|}v(x)f(x)\big{|}dx\leq
K\int_{\mathbb{R}\setminus[-M,M]}f(x)dx\leq\frac{2KC_{3}}{\alpha}M^{-\alpha}\leq\frac{\varepsilon}{2}$

is satisfied by the definition of $M$.  A little calculation shows

$L\geq\left(\frac{\sqrt{2}6C_{3}\xi}{\sqrt{1+2\alpha}\varepsilon}\right)^{\frac{2}{1+2\alpha}},$

which implies

$\displaystyle\big{\|}f-f_{L}\big{\|}_{2}$
$\displaystyle\leq\sqrt{2\int_{L}^{\infty}C_{3}^{2}x^{-2-2\alpha}dx}$
$\displaystyle=\frac{\sqrt{2}C_{3}}{\sqrt{1+2\alpha}}L^{-\frac{1+2\alpha}{2}}$
$\displaystyle\leq\frac{\varepsilon}{6\xi}.$

The density $f$ satisfies the condition of Lemma 3.8 by Example 3.9. Using the
bound for $B(L)$ from Lemma 3.8, we obtain

$\displaystyle\sqrt{B(L)}$
$\displaystyle\leq\sqrt{\frac{1}{L}\left(2\int_{L}^{\infty}C_{3}x^{-1-\alpha}dx\right)^{2}+\frac{8}{3}LC_{3}^{2}L^{-2\alpha-2}}$
$\displaystyle=\sqrt{\frac{4C_{3}^{2}}{\alpha^{2}}L^{-2\alpha-1}+\frac{8}{3}C_{3}^{2}L^{-2\alpha-1}}$
$\displaystyle=2C_{3}\sqrt{\frac{1}{\alpha^{2}}+\frac{2}{3}}L^{-\frac{1+2\alpha}{2}}$
(34) $\displaystyle\leq\frac{\varepsilon}{6\xi},$

the last Inequality holds by the definition of $L$. As in the proof of Theorem
3.6 one can show that by the definition of $N$

$\left\|f_{L}-\sum_{k=0}^{N}{}^{\prime}a_{k}e_{k}\right\|_{2}\leq\frac{\varepsilon}{6\xi}.$

As $B(L)\to 0$, $L\to\infty$, $f$ is COS-admissible. Apply Lemma 2.1 to
conclude. ∎

######  Example 3.11.

Recall that the log-returns of the FMLS model are stably distributed. Let
$F_{\text{Stable}}$ be the cumulative distribution function for the stable
density $f_{\text{Stable}}$ with stability parameter $\alpha\in(0,2)$,
skewness $\beta\in[-1,1]$, scale $c>0$ and location $\theta\in\mathbb{R}$. By
Samorodnitsky and Taqqu (2017, Property 1.2.15) it holds that

$\displaystyle\lim_{x\to\infty}x^{\alpha}(1-F_{\text{Stable}}(x))$
$\displaystyle=C_{\alpha}\frac{1+\beta}{2}c^{\alpha},$
$\displaystyle\lim_{x\to\infty}x^{\alpha}(F_{\text{Stable}}(-x))$
$\displaystyle=C_{\alpha}\frac{1-\beta}{2}c^{\alpha},$

where

$C_{\alpha}=\begin{cases}\frac{1-\alpha}{\Gamma(2-\alpha)\cos\big{(}\frac{\pi\alpha}{2}\big{)}}&,\alpha\neq
1\\\ \frac{2}{\pi}&,\alpha=1.\end{cases}$

For stable densities we therefore suggest to set $C_{3}$ in Equation (17) at
least as large as $\alpha C_{\alpha}\frac{1+|\beta|}{2}c^{\alpha}$ to obtain
$L$. The number of terms $N$ can be found using Equation (32) and the bounds
$H_{j}$ are explicitly given in Example 3.2.

##  4 Order of Convergence

The main Theorem 4.1 in this section describes the order of convergence of the
COS method if we allow $N\to\infty$ and choose $M$ and $L$ depending on $N$.
We describe only the asymptotic behavior of the COS method and we assume $M=L$
in this section to keep the notation simple. We establish a bound of the order
of convergence of the error of the COS method with parameters $L$ and $N$,
i.e.,

$\text{err}(L,N)=\left|\int_{\mathbb{R}}v(x)f(x)dx-\sum_{k=0}^{N}{}^{\prime}c_{k}v_{k}\right|.$
(35)

We use big $O$ notation: for functions $g,h:\mathbb{N}\to\mathbb{R}$, we write
$h(N)=O(g(N))$ as $N\to\infty$, if the absolute value of $h(x)$ is at most a
positive constant multiple of $g(x)$ for all sufficiently large values of $N$.

######  Theorem 4.1.

Let $v:\mathbb{R}\to\mathbb{R}$ be bounded, with $|v(x)|\leq K$ for all
$x\in\mathbb{R}$ and some $K>0$. Assume $J\in\mathbb{N}$.
(i) Let $f$ be a density satisfying Assumption A. Let
$\beta\in\left(0,\frac{2J}{2J+3}\right)$ and $\gamma>0$. If $L=\gamma
N^{\beta}$ it holds that

$\text{err}(L(N),N)\leq
O\left(N^{-J(1-\beta)+\frac{3}{2}\beta}\right),\quad\text{as }N\to\infty.$

(ii) Let $f$ be a unimodal density satisfying Assumption B. Let
$\beta\in\left(0,\frac{2J}{2J+3+2\alpha}\right)$ and $\gamma>0$. If $L=\gamma
N^{\beta}$ it holds that

$\text{err}(L(N),N)\leq O\left(N^{-\alpha\beta}\right),\quad\text{as
}N\to\infty.$

In both cases we have $\text{err}(L(N),N)\to 0$, $N\to\infty.$

###### Proof.

Let $v_{L}:=1_{[-L,L]}v$. As in the proof of Corollary 8 in Junike and
Pankrashkin (2022), one can show that

$\displaystyle\Big{|}\int_{\mathbb{R}}v(x)f(x)dx-\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N}c_{k}v_{k}\Big{|}\leq$
$\displaystyle A_{0}(L)+\|v_{L}\|_{2}(A_{1}(L)+A_{2}(L,N)+\sqrt{B(L)}),$ (36)

where $B(L)$ as in Equation (5) and

$\displaystyle A_{0}(L)=$
$\displaystyle\int_{\mathbb{R}\setminus[-L,L]}\big{|}v(x)f(x)\big{|}dx$
$\displaystyle A_{1}(L)=$ $\displaystyle\left\|f-f_{L}\right\|_{2}$
$\displaystyle A_{2}(L,N)=$
$\displaystyle\big{\|}f_{L}-\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N}a_{k}e_{k}\big{\|}_{2}.$

We will state upper bounds for $A_{0}$, $A_{1}$ and $B$ depending on the tail
behaviour of $f$, i.e. for the different cases i) and ii) in Theorem 4.1. An
upper bound for $A_{2}$ can be obtained from Lemma 3.1. Note that

$\|v_{L}\|_{2}\leq\sqrt{2L}K=\sqrt{2\gamma}KN^{\frac{\beta}{2}}.$

We now proof (i). We have for $L$ large enough

$A_{0}(L)\leq
2KC_{1}\int_{L}^{\infty}e^{-C_{2}x}dx=\frac{2KC_{1}}{C_{2}}e^{-C_{2}L}.$

Further, by Inequality (23), it holds that

$\displaystyle A_{1}(L)$
$\displaystyle\leq\frac{C_{1}}{\sqrt{C_{2}}}e^{-C_{2}L}.$

Assuming $L\geq 1$ and applying Inequality (25), it follows that

$\displaystyle\sqrt{B(L)}\leq$ $\displaystyle\frac{2\pi
C_{1}}{\sqrt{6C_{2}}}\sqrt{1+\frac{1}{C_{2}}+\frac{1}{2C_{2}^{2}}}e^{-C_{2}L}.$

Inequality (11) implies

$\displaystyle A_{2}(L,N)$
$\displaystyle\leq\sum_{j=1}^{J}\frac{2^{j+1}}{j\pi^{j+1}}\frac{L^{j}}{N^{j}}\left(2C_{1}e^{-C_{2}L}\right)+\frac{2^{J+2}H_{J+1}}{J\pi^{J+1}}\frac{L^{J+1}}{N^{J}}$
$\displaystyle\leq
e^{-C_{2}L}\sum_{j=1}^{J}\frac{2^{j+2}C_{1}\gamma^{j}}{j\pi^{j+1}}+\frac{2^{J+2}H_{J+1}}{J\pi^{J+1}}\frac{L^{J+1}}{N^{J}},$

where we used $\frac{L}{N}\leq\gamma$.

Let $b_{1},...,b_{4}>0$ be suitable constants. By Inequality (36), it follows
for $N$ large enough

$\displaystyle\Big{|}\int_{\mathbb{R}}v(x)f(x)dx-\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N}c_{k}^{L}v_{k}^{M}\Big{|}\leq$
$\displaystyle b_{1}e^{-C_{2}\gamma
N^{\beta}}+b_{2}N^{\frac{\beta}{2}}\big{(}b_{3}e^{-C_{2}\gamma
N^{\beta}}+b_{4}N^{\beta(J+1)-J}\big{)}$ $\displaystyle\leq$ $\displaystyle
O\left(N^{-J(1-\beta)+\frac{3}{2}\beta}\right),\quad N\to\infty.$ (37)

$\beta<\frac{2J}{2J+3}$ implies $-J(1-\beta)+\frac{3}{2}\beta<0$ and the
right-hand-side of (37) converges to zero for $N\to\infty$.

We show (ii). It holds for $L$ large enough that

$\displaystyle A_{0}(L)$ $\displaystyle\leq
2KC_{3}\int_{L}^{\infty}x^{-1-\alpha}dx\leq\frac{2KC_{3}}{\alpha}L^{-\alpha},$
$\displaystyle A_{1}(L)$
$\displaystyle\leq\sqrt{2C_{3}^{2}\int_{L}^{\infty}x^{-2(1+\alpha)}dx}=\sqrt{\frac{2C_{3}^{2}}{1+2\alpha}}L^{-\frac{1}{2}-\alpha}$

and using Inequality (34)

$\displaystyle\sqrt{B(L)}$ $\displaystyle\leq
2C_{3}\sqrt{\frac{1}{\alpha^{2}}+\frac{2}{3}}L^{-\frac{1}{2}-\alpha}.$

Inequality (11) implies

$\displaystyle A_{2}(L,N)$
$\displaystyle\leq\sum_{j=1}^{J}\frac{2^{j+2}C_{3}\gamma^{-1-\alpha}}{j\pi^{j+1}}N^{\beta(-1-\alpha)-j}+\frac{2^{J+2}H_{J+1}}{J\pi^{J+1}}\frac{L^{J+1}}{N^{J}}.$

Hence, by Inequality (36), it follows for $N$ large enough and some some
suitable constants $b_{1},...,b_{4}>0$ and $a_{1},...,a_{J}>0$ that

$\displaystyle\Big{|}\int_{\mathbb{R}}v(x)f(x)dx-\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N}c_{k}v_{k}\Big{|}$
$\displaystyle\leq
b_{1}N^{-\alpha\beta}+b_{2}N^{\frac{\beta}{2}}\bigg{(}b_{3}N^{-\beta(\alpha+\frac{1}{2})}+\sum_{j=1}^{J}a_{j}N^{\beta(-1-\alpha)-j}+b_{4}N^{\beta(J+1)-J}\bigg{)}$
$\displaystyle\leq
b_{1}N^{-\alpha\beta}+b_{2}\bigg{(}b_{3}N^{-\beta\alpha}+\sum_{j=1}^{J}a_{j}N^{-\beta(\alpha+\frac{1}{2})-j}+b_{4}N^{\beta\left(J+\frac{3}{2}\right)-J}\bigg{)}$
$\displaystyle\leq O\left(N^{-\alpha\beta}\right),\quad N\to\infty.$

The last inequality can be seen as follows: as $0<\beta$, we have

$-\beta\alpha\geq-\beta\big{(}\alpha+\frac{1}{2}\big{)}-j,\quad j=1,...,J.$

Further, as $\beta\leq\frac{J}{J+\frac{3}{2}+\alpha}$, it follows

$-\alpha\beta\geq\beta(J+\frac{3}{2})-J.$

∎

######  Remark 4.2.

The COS method converges exponentially, i.e. faster than $O(N^{-q})$ for any
$q>0$, if $f$ is smooth and has semi-heavy tails. To see this, let $0<\beta<1$
then

$-J(1-\beta)+\frac{3}{2}\beta\to-\infty,\quad J\to\infty.$

######  Remark 4.3.

By Theorem 3.10, the COS method converges at least like $N^{-\alpha}$ if $f$
is smooth and has heavy tails with Pareto index $\alpha>0$. Numerical
experiments indicate that the COS method does not converge faster than
$O(N^{-\alpha})$ for the FMLS model, see Section 5.4, i.e. the bound in
Theorem 3.10(ii) is sharp.

##  5 Numerical experiments

All numerical experiments are compared with the Carr-Madan formula, see Carr
and Madan (1999), which is applicable when the characteristic function of the
log-returns is given in closed-form and when $E[S_{T}^{1+\gamma}]$ is finite
for some $\gamma>0$, called _damping factor_. Unless otherwise stated, we use
the Carr-Madan formula with $N=2^{17}$ terms, we set the damping factor equal
to $\gamma=0.1$, and we use a Fourier truncation range of $1200$ to compute
the reference prices. We implemented the Carr-Madan formula using Simson’s
rule without applying the fast Fourier transform.

All numerical experiments are performed on a modern laptop (Intel i7-10750H)
using the software R and vectorized code without parallelization.

###  5.1 BS model

How sharp is the bound for $N$ in Theorem 3.6? For the BS model, we compute
$L$ by Equation (17) and $N$ by Equation (19) and compare $N$ with $N_{\min}$,
which is the minimal number of terms such that the absolute difference of the
price of the COS method and the Black-Scholes price are less than the error
tolerance. Table 1 shows that $N$ is only about $40\%$ larger than $N_{\min}$.

How can the number of terms $N$ be estimated if there are no closed-form
solutions available for the bounds of the derivatives of the density of the
log-returns? We suggest solving the right-hand-side of Inequality (15)
numerically to obtain $H_{j+1}$. Here, we use R’s _integrate_ function with
default values. The CPU time of the COS method and the numerical integration
routine to obtain $H_{j+1}$ are of similar magnitude, see Table 1.

The value for $N$ does not change when using numerical integration to obtain
$H_{j+1}$ compared to the closed-form solution for $H_{j+1}$ from Example 3.3.
In this example, $j=40$ works well. We also we use $j=40$ in Section 5.3.

$j$ | 10 | 20 | 30 | 40 | 50 | 60 | 70
---|---|---|---|---|---|---|---
$N$ | 897 | 271 | 200 | 179 | 172 | 170 | 171
CPU time COS method | 0.18 | 0.09 | 0.09 | 0.07 | 0.08 | 0.07 | 0.08
CPU time numeric integration for $H_{j+1}$ | 0.15 | 0.09 | 0.09 | 0.10 | 0.09 | 0.10 | 0.09
Table 1: $L$ and $N$ in the BS model for a put option with $T=1$,
$\sigma=0.2$, $r=0$, $K=S_{0}=100$ and $\varepsilon=10^{-8}$ compared to the
minimal number of terms $N_{\min}=120$ of the COS method to stay below the
error tolerance. We used $n=8$ moments to obtain $L$. CPU time is measured in
milliseconds. The COS method with $N_{\min}=120$ takes about 0.068
milliseconds.

###  5.2 VG model

Using the VG model as an example, this section shows that Theorem 3.6 does not
help in finding the number of terms $N$ for the COS method if the density of
the log-returns is not sufficiently smooth.

Let $f_{\text{VG}}$ denote the density of the log-returns in the VG model at
maturity $T>0$ with parameters $\sigma>0$, $\theta=0$ and $\nu>0$. If
$T<\frac{\nu}{2}$, the density $f_{\text{VG}}$ is unbounded and Theorem 3.6
cannot be used to find $N$. If $T\in\left(\nu,\frac{3\nu}{2}\right)$, the
derivative of $f_{\text{VG}}$ is continuous, but the second derivative of
$f_{\text{VG}}$ is unbounded, see Küchler and Tappe (2008, Thm. 4.1 and Thm.
6.1).

We can apply Theorem 3.6 with $J=0$ if $T\in\left(\nu,\frac{3\nu}{2}\right)$
but the value for $N$ is somewhat useless from a practical point of view as
the following experiment shows: Consider a European call option and the
following parameters:

$\varepsilon=0.01,\quad\sigma=0.1,\quad\nu=0.2,\quad T=0.25,\quad
S_{0}=K=100,\quad r=0.$

By Equation (17), we set $L=0.91$ using $n=4$ moments. By numerically
optimizing the derivative of the density $f_{\text{VG}}$, we obtain
$H_{1}=218$ and Theorem 3.6 suggests $N\approx 8\cdot 10^{12}$.

We calulated a reference price $\pi_{CM}=1.809833$ using the Carr-Madan
formula. Using the COS method, already $N=50$ is sufficient to approximate the
reference price within the error tolerance.

###  5.3 FMLS model

In this section, we consider the FMLS model with parameters $\alpha\in(1,2)$
and $\sigma>0$ at maturity $T>0$. Recall that the centralized log-returns at
time $T>0$ have heavy tails and are stably distributed with stability
parameter $\alpha$, skewness $\beta=-1$, scale $c=\sigma T^{\frac{1}{\alpha}}$
and location $\theta=0$.

Carr and Wu (2003) calibrated the FMLS model to real market data and estimated
$\alpha=1.5597$ and $\sigma=0.1486$. We test the formulas for $L$ and $N$ for
the FMLS model with these parameters for a European call option with the
following parameters:

$\varepsilon=10^{-2},\quad T=1,\quad S_{0}=K=100,\quad r=0.$

To apply the COS method we define $M$ and $L$ as in Equations (30) and (31).
Motivated by Example 3.11, we suggest to set the variable $C_{3}$ appearing in
Equation (30) to

$C_{3}:=\alpha\frac{1-\alpha}{\Gamma(2-\alpha)\cos\big{(}\frac{\pi\alpha}{2}\big{)}}\sigma^{\alpha}T.$

Figure 1 shows the density of the log-returns in the FMLS model and asymptotic
tail behavior, i.e. the function $x\mapsto C_{3}|x|^{-1-\alpha}$. The left
tail does indeed decay like Pareto tails, the asymptotic tail behavior is very
close to the density. The right tail decays faster, it decays exponentially,
see Carr and Wu (2003).

Figure 1: Density of the log-returns in the FMLS model and asymptotic tail
behavior.

The reference price is 9.7433708 by the Carr-Madan formula and we confirm the
reference price by the COS method with $L=10^{5}$ and $N=10^{7}$. Using the
explicit formula for $H_{j+1}$ from Example 3.4 and Equation (32) find
$N=5451$, where we set $j=40$. According to Theorem 3.10, any other choice for
$j$ is possible but another value for $j$ does not significantly improve $N$.

Applying Equations (30) and (31), we get $M=69$ and $L=176$. With these
parameters, the COS method prices the call option within the error tolerance
in about 1.5 milliseconds. The minimal $N$ to stay below the error tolerance
is $N_{\min}=1200$ and the CPU time using $N_{\min}$ is about 0.4
milliseconds.

We also apply the Carr-Madan formula with the “default parameters”, see Carr
and Madan (1999), which are also recommended by Madan and Schoutens (2016,
Sec. 3.1), i.e. $4096$ terms, a damping factor equal to $1.5$ and a Fourier
truncation range of $1024$. With these default parameters, the Carr-Madan
formula prices the option within the error tolerance in about 1.1
milliseconds.

The value for the CPU time is about three times higher when using Theorem 3.10
to get $N$ compared to the optimized value for $N$, which is acceptable from
our point of view. The computational time of the COS method using Theorem 3.10
and the Carr-Madan formula with standard parameters are of similar magnitude.

###  5.4 Order of Convergence

In this section, we empirically analyze the order of convergence of the COS
method for the BS model with parameter $\sigma>0$ and the FMLS model with
parameters $\alpha\in(1,2)$ and $\sigma>0$ and compare the results with
Theorem 4.1. We also consider the Cauchy density.

In the BS model, the density of the log-returns is arbitrarily smooth and the
tails decay even faster than exponentially. In Figure 2 we analyze how the
error of the COS method behaves in the BS model for large $N$ and different
choices of $L$.

In particular, we see that the COS method seems to converge exponentially at
the beginning for moderate $N$ if we choose $L$ constant, i.e. independent of
$N$. But for constant $L$, e.g. $L=4\sigma$ or $L=6\sigma$, the COS method
does not converge for $N\to\infty$ but the error remains constant for a
certain level of $N$. This can be explained by Inequality (1). While the
second term on the right-hand-side of Inequality (1) converges to zero for
$N\to\infty$ and $L$ fixed, the first and third terms on the right-hand-side
of Inequality (1) do not improve as $N$ is increased but $L$ is kept constant.

For large $L$, such as $L=20\sigma$, this effect disappears somewhat because
the tails of the Gaussian distribution decay so rapidly that, up to fixed-
precision arithmetic111The software package R operates with a precision of 53
bits conforming to the IEC 60559 standard, see R Core Team (2021). R: A
language and environment for statistical computing. R Foundation for
Statistical Computing, Vienna, Austria. URL https://www.R-project.org/., the
Gaussian density can be thought of as a density with finite support. Using
arbitrary-precision arithmetic instead should show that even for $L=20\sigma$,
the error of the COS method does not converge to zero for $N\to\infty$.

If we choose $L=L(N)=\sigma\sqrt{N}$, Theorem 4.1(i) indicates that the error
of the COS method converges exponentially to zero. This is confirmed
empirically in Figure 2. Other choices for $L$ also work well, e.g.
$L=\frac{\sigma}{5}N$.

Figure 2: Order of Convergence of the COS method for a call option in the BS
model with $\sigma=0.2$, $T=1$, $K=S_{0}=100$, $r=0$ and different choices for
the truncation range $[-L,L]$.

We test Theorem 4.1(ii) for the density of the Cauchy distribution and the
density of the log-returns in the FMLS model. Both densities belong to the
family of stable densities and have heavy tails. In the case of the Cauchy
density, we set

$f_{\text{Cauchy}}(x)=\frac{1}{\pi(1+x^{2})}\quad\text{and}\quad
v(x)=1_{(-\infty,1.23]}(x),\quad x\in\mathbb{R}.$

The integral $\int_{\mathbb{R}}f_{\text{Cauchy}}(x)v(x)dx$ is equal to the
value of the cumulative distribution function of $f_{\text{Cauchy}}$ at $1.23$
and is given in closed-form. The Pareto index of $f_{\text{Cauchy}}$ is
$\alpha=1$.

For the FMLS model we use the parameters $\alpha=1.5597$ and $\sigma=0.1486$,
these values are taken from Carr and Wu (2003), who calibrated the FMLS model
to real market data. We then compute the price of a European call option with
maturity $T=1$, $S_{0}=K=100$ and $r=0$ using the Carr-Madan formula.

For the Cauchy and the FMLS models, we compute $L_{\text{optimal}}$ for a
fixed $N\in\mathbb{N}$ minimizing $\text{err}(L,N)$, i.e. for each $N$ we
choose the truncation range such that the error of the COS method is minimal,
compare with Equation (35). We then compare the empirical order of convergence
of $\text{err}(L_{\text{optimal}}(N),N)$ with the theoretical bound
$O(N^{-\alpha})$ from Theorem 4.1(ii). In particular, for each
$N\in\\{2^{i},\,i=4,5,...,20\\}$ we define

$L_{\text{optimal}}(N):=\text{argmin}_{L\in\mathbb{L}}\text{err}(L,N),$

where

$\mathbb{L}=\\{\exp(0.07i),\quad i=0,1,...,200\\},$

is a sufficiently fine grid of the interval $[1,10^{6}]$. Figure 3 shows
$L_{\text{optimal}}$ over $N$ for the Cauchy and FMLS model. We see that the
relationship between $\log(N)$ and $\log(L_{\text{optimal}})$ is approximately
linear with slope $0.95$ for the Cauchy model and slope $0.86$ for the FMLS
model.

In Figures 4 and 5 we compute the order of convergence of the COS method for
different strategies to define $L$. If $L$ is constant, the COS method does
not converge at all.

For the Cauchy model, the order of convergence can be estimated from Figure 4
to be $O(N^{-0.92})$, if we choose $L=\frac{1}{10}N$. The order of convergence
using $L_{\text{optimal}}$ is of similar magnitude. This order of convergence
is very close to the theoretical bound of $O(N^{-1})$ predicted by Theorem
4.1(ii).

For the FMLS model the empirical order of convergence is $O(N^{-1.57})$ if we
choose $L=\frac{1}{100}N$, which is also very close to its theoretical bound
of $O(N^{-1.5597})$ and does not differ much if we use $L_{\text{optimal}}$
instead.

In particular the numerical experiments indicate that the order of convergence
is not exponential for heavy tail models even though the corresponding
densities are arbitrarily smooth.

In Figure 5 we also plot the empirical order of convergence of the Carr-Madan
formula for the FMLS model with damping factor of $0.1$ and a Fourier
truncation range of $1200$. The Figure indicates an exponential order of
convergence for the Carr-Madan formula.

Figure 3: $L_{\text{optimal}}$ for Cauchy and FMLS. The horizontal lines
corresponds to the grid $\mathbb{L}$ where we searched for the optimal
truncation range. Figure 4: Order of convergence of the COS method for the
Cauchy model with different choices for the truncation range $[-L,L]$. Figure
5: Order of convergence of the COS method for a call option in the FMLS model
with $\alpha=1.5597$, $\sigma=0.1486$, $T=1$, $K=S_{0}=100$, $r=0$ and
different choices for the truncation range $[-L,L]$. The reference price is
9.743370825229…

##  6 Conclusions

In this research we analyzed the COS method, which is used for efficient
option pricing. The COS method requires two parameters: the truncation range
$[-L,L]$ to truncate the density of the log-returns and the number of terms
$N$ to approximate the truncated density by cosine functions. We considered
stock price models where the density of the log-returns is smooth and has
either semi-heavy tails, i.e. the tails decay exponentially, or heavy tails,
i.e. Pareto tails.

In both cases, we found explicit and useful bounds for $L$ and $N$ and showed
in numerical experiments the usefulness of these formulas in applications. The
densities of the log-returns are smooth for many models in finance, such as
the BS, NIG, Heston or FMLS models.

If the density is not differentiable, neither Theorem 3.6 nor 3.10 can be used
to find a bound for $N$. If the density is only differentiable a few times,
which is the case for the VG model for some parameters and short maturities,
our bound for $N$ is too large to be useful in most practical applications. A
solution might be to use the NIG model instead of the VG model. Both are
flexible, three-parameter Lévy models.

We further analyzed the order of convergence of the COS method and observed
both theoretically and empirically that the models enjoy exponential
convergence when the densities of the log-returns are smooth and have semi-
heavy tails. However, when the density of the log-returns is smooth and has
_heavy_ _tails_ , the order of convergence of the COS method can be bounded by
$O(N^{-\alpha})$, where $\alpha>0$ is the Pareto tail index. This is the case,
for example, for the FMLS model where $\alpha\in(1,2)$. Numerical experiments
indicate that the bound $O(N^{-\alpha})$ is sharp.

## References

  * Aimi et al. (2023) Alessandra Aimi, Chiara Guardasoni, Luis Ortiz-Gracia, and Simona Sanfelici.  Fast Barrier Option Pricing by the COS BEM Method in Heston Model (with Matlab Code).  _Computational Methods in Applied Mathematics_ , (0), 2023\.
  * Albin and Sundén (2009) J. M. P. Albin and Mattias Sundén.  On the asymptotic behaviour of Lévy processes, Part I: Subexponential and exponential processes.  _Stochastic processes and their applications_ , 119(1):281–304, 2009.
  * Asmussen (2022) Søren Asmussen.  On the role of skewness and kurtosis in tempered stable (CGMY) Lévy models in finance.  _Finance and Stochastics_ , 26(3):383–416, 2022\.
  * Bardgett et al. (2019) Chris Bardgett, Elise Gourier, and Markus Leippold.  Inferring volatility dynamics and risk premia from the S&P 500 and VIX markets.  _Journal of Financial Economics_ , 131(3):593–618, 2019.
  * Barndorff-Nielsen (1997) Ole E Barndorff-Nielsen.  Normal inverse Gaussian distributions and stochastic volatility modelling.  _Scandinavian Journal of Statistics_ , 24(1):1–13, 1997.
  * Black and Scholes (1973) Fischer Black and Myron Scholes.  The pricing of options and corporate liabilities.  _Journal of Political Economy_ , 81(3):637–654, 1973.
  * Boyd (2001) John P Boyd.  _Chebyshev and Fourier spectral methods_.  Courier Corporation, 2001.
  * Carr and Madan (1999) Peter Carr and Dilip Madan.  Option valuation using the fast Fourier transform.  _Journal of Computational Finance_ , 2(4):61–73, 1999.
  * Carr and Wu (2003) Peter Carr and Liuren Wu.  The finite moment log stable process and option pricing.  _The Journal of Finance_ , 58(2):753–777, 2003\.
  * Carr et al. (2002) Peter Carr, Hélyette Geman, Dilip Madan, and Marc Yor.  The fine structure of asset returns: An empirical investigation.  _The Journal of Business_ , 75(2):305–332, 2002\.
  * Dragulescu and Yakovenko (2002) Adrian A Dragulescu and Victor M Yakovenko.  Probability distribution of returns in the Heston model with stochastic volatility.  _Quantitative Finance_ , 2(6):443, 2002.
  * Fama (1965) Eugene F Fama.  The behavior of stock-market prices.  _The Journal of Business_ , 38(1):34–105, 1965\.
  * Fang and Oosterlee (2009a) Fang Fang and Cornelis W Oosterlee.  A novel pricing method for European options based on Fourier-cosine series expansions.  _SIAM Journal on Scientific Computing_ , 31(2):826–848, 2009a.
  * Fang and Oosterlee (2009b) Fang Fang and Cornelis W Oosterlee.  Pricing early-exercise and discrete barrier options by Fourier-cosine series expansions.  _Numerische Mathematik_ , 114(1):27, 2009b.
  * Fang and Oosterlee (2011) Fang Fang and Cornelis W Oosterlee.  A Fourier-based valuation method for Bermudan and barrier options under Heston’s model.  _SIAM Journal on Financial Mathematics_ , 2(1):439–463, 2011.
  * Grzelak and Oosterlee (2011) Lech A Grzelak and Cornelis W Oosterlee.  On the Heston model with stochastic interest rates.  _SIAM Journal on Financial Mathematics_ , 2(1):255–286, 2011.
  * Heston (1993) Steven L Heston.  A closed-form solution for options with stochastic volatility with applications to bond and currency options.  _The Review of Financial Studies_ , 6(2):327–343, 1993.
  * Junike and Pankrashkin (2022) Gero Junike and Konstantin Pankrashkin.  Precise option pricing by the cos method–how to choose the truncation range.  _Applied Mathematics and Computation_ , 421:126935, 2022\.
  * Küchler and Tappe (2008) Uwe Küchler and Stefan Tappe.  On the shapes of bilateral Gamma densities.  _Statistics & Probability Letters_, 78(15):2478–2484, 2008.
  * Küchler and Tappe (2013) Uwe Küchler and Stefan Tappe.  Tempered stable distributions and processes.  _Stochastic Processes and their Applications_ , 123(12):4256–4293, 2013.
  * Leitao et al. (2018) Álvaro Leitao, Cornelis W Oosterlee, Luis Ortiz-Gracia, and Sander M Bohte.  On the data-driven COS method.  _Applied Mathematics and Computation_ , 317:68–84, 2018\.
  * Liu et al. (2019a) Shuaiqiang Liu, Anastasia Borovykh, Lech A Grzelak, and Cornelis W Oosterlee.  A neural network-based framework for financial model calibration.  _Journal of Mathematics in Industry_ , 9:1–28, 2019a.
  * Liu et al. (2019b) Shuaiqiang Liu, Cornelis W Oosterlee, and Sander M Bohte.  Pricing options and computing implied volatilities using neural networks.  _Risks_ , 7(1):16, 2019b.
  * Lord et al. (2008) Roger Lord, Fang Fang, Frank Bervoets, and Cornelis W Oosterlee.  A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes.  _SIAM Journal on Scientific Computing_ , 30(4):1678–1705, 2008.
  * Lyness (1971) J. N. Lyness.  Adjusted forms of the Fourier coefficient asymptotic expansion and applications in numerical quadrature.  _Mathematics of Computation_ , 25(113):87–104, 1971.
  * Madan and Schoutens (2016) Dilip Madan and Wim Schoutens.  _Applied Conic Finance_.  Cambridge University Press, 2016.
  * Madan et al. (1998) Dilip Madan, Peter P Carr, and Eric C Chang.  The variance gamma process and option pricing.  _Review of Finance_ , 2(1):79–105, 1998.
  * Mandelbrot (1963) Benoit B Mandelbrot.  The variation of certain speculative prices.  _The Journal of Business_ , 38(1):394–419, 1963\.
  * Merton (1976) Robert C Merton.  Option pricing when underlying stock returns are discontinuous.  _Journal of Financial Economics_ , 3(1-2):125–144, 1976.
  * Oosterlee and Grzelak (2019) Cornelis W Oosterlee and Lech A Grzelak.  _Mathematical modeling and computation in finance: with exercises and Python and MATLAB computer codes_.  World Scientific, 2019.
  * Ortiz-Gracia and Oosterlee (2013) Luis Ortiz-Gracia and Cornelis W Oosterlee.  Robust pricing of European options with wavelets and the characteristic function.  _SIAM Journal on Scientific Computing_ , 35(5):B1055–B1084, 2013.
  * Ortiz-Gracia and Oosterlee (2016) Luis Ortiz-Gracia and Cornelis W Oosterlee.  A highly efficient Shannon wavelet inverse Fourier technique for pricing European options.  _SIAM Journal on Scientific Computing_ , 38(1):B118–B143, 2016.
  * Plonka et al. (2018) Gerlind Plonka, Daniel Potts, Gabriele Steidl, and Manfred Tasche.  _Numerical Fourier Analysis_.  Springer, 2018.
  * Ruijter and Oosterlee (2012) Marjon J Ruijter and Cornelis W Oosterlee.  Two-dimensional Fourier cosine series expansion method for pricing financial options.  _SIAM Journal on Scientific Computing_ , 34(5):B642–B671, 2012.
  * Samorodnitsky and Taqqu (2017) Gennady Samorodnitsky and Murad S Taqqu.  _Stable non-Gaussian random processes: stochastic models with infinite variance: stochastic modeling_.  Routledge, 2017.
  * Schoutens (2003) Wim Schoutens.  _Lévy Processes in Finance: Pricing Financial Derivatives_.  Wiley Online Library, 2003.
  * Wright et al. (2015) Grady B Wright, Mohsin Javed, Hadrien Montanelli, and Lloyd N Trefethen.  Extension of Chebfun to periodic functions.  _SIAM Journal on Scientific Computing_ , 37(5):C554–C573, 2015.
  * Yamazato (1978) Makoto Yamazato.  Unimodality of infinitely divisible distribution functions of class L.  _The Annals of Probability_ , (6):523–531, 1978\.
  * Zhang and Oosterlee (2013) Bowen Zhang and Cornelis W Oosterlee.  Efficient pricing of European-style Asian options under exponential Lévy processes based on Fourier cosine expansions.  _SIAM Journal on Financial Mathematics_ , 4(1):399–426, 2013.