# Singular Fourier-Padé Series Expansion of European Option Prices

Tat Lung (Ron) Chan∗${\dagger}$
∗Corresponding author. Email<EMAIL_ADDRESS>${\dagger}$University of East
London, Water Lane, Stratford, UK, E15 4LZ

(v1.1 released May 2016)

###### Abstract

We apply a new numerical method, the singular Fourier-Padé (SFP) method
invented by Driscoll and Fornberg (2001, 2011), to price European-type options
in Lévy and affine processes. The motivation behind this application is to
reduce the inefficiency of current Fourier techniques when they are used to
approximate piecewise continuous (non-smooth) probability density functions.
When techniques such as fast Fourier transforms and Fourier series are applied
to price and hedge options with non-smooth probability density functions, they
cause the Gibbs phenomenon; accordingly, the techniques converge slowly for
density functions with jumps in value or derivatives. This seriously adversely
affects the efficiency and accuracy of these techniques. In this paper, we
derive pricing formulae and their option Greeks using the SFP method to
resolve the Gibbs phenomenon and restore the global spectral convergence rate.
Moreover, we show that our method requires a small number of terms to yield
fast error convergence, and it is able to accurately price any European-type
option deep in/out of the money and with very long/short maturities.
Furthermore, we conduct an error-bound analysis of the SFP method in option
pricing. This new method performs favourably in numerical experiments compared
with existing techniques.

###### keywords:

singular Fourier-Padé series, European-type options, Lévy processes, affine
processes

##  1 Introduction

The application of Fourier techniques to pricing and hedging financial
derivatives has flourished in computational and mathematical finance over the
last decade because most of the underlying processes driving derivatives have
a characteristic function (i.e., the Fourier transform of the probability
density function (PDF)). In general, two Fourier techniques are applied to
option pricing: fast Fourier transforms (FFTs) and Fourier series, such as
Fourier cosine (COS) series. To explain the first technique, FFTs in option
pricing, we first consider a European option pricing formula that is a risk-
neutral expectation of the discounted payoff of the option driven by an
underlying process. This formation naturally implies the integration of the
discounted payoff and the PDF of the underlying process. As the underlying
process has a characteristic function, one can apply an inverse continuous
Fourier transform to recover the PDF, discretise the continuous Fourier
integral and use an FFT algorithm that computes discrete Fourier transforms in
$O(N\log_{2}(N))$ operations to obtain the option prices (cf. Carr and Madan,
1999; Lewis, 2001; Lipton, 2002; Chourdakis, 2004; Itkin, 2005; Lord et al.,
2008). The main disadvantage of using these methods is that they require
thousands of grid points and considerable computational time to reach an
acceptable level of accuracy. Another approach is to use Fourier series to
expand a PDF in terms of a partial sum of orthogonal basis functions. The
basis function can be complex an exponential, e.g., Chan (2016), or COS, e.g.
Fang and Oosterlee (2009a), function. Each coefficient in the series is
approximated via the characteristic function of the PDF. By integrating each
basis function with the discounted payoff function of the option, we can
obtain an option price. This technique originated in Fang and Oosterlee
(2009a) to price European-style options with COS basis functions. The success
of Fourier COS series in European option pricing has lead to the rapid
expansion of pricing and hedging options with early-exercise features and
exotic options, like Asian, multi-asset or barrier options (e.g., Leentvaar
and Oosterlee, 2008; Fang and Oosterlee, 2009b; Fang and Oosterlee, 2011;
Zhang and Oosterlee, 2013). Compare with FFTs, the advantage of using Fourier
series for option pricing is that they can achieve a global spectral
(exponential) convergence rate and require fewer summation terms as long as
the governing PDF is sufficiently smooth. However, using any type of Fourier
series to represent a $\mathbf{C}^{\nu}$111A vector space in which functions
are $\nu$ continuously differentiable. piecewise continuous (non-smooth)
function222A function is called piecewise continuous on an interval if the
function is made up of a finite number of $\nu$ times differentiable
continuous pieces., e.g., a non-smooth PDF, is notoriously fraught.
Discontinuities cause the Gibbs phenomenon, which has two important
consequences for the Fourier partial sum of length $N$:

  1. 1.

failure to converge at the jump, and

  2. 2.

pointwise convergence elsewhere at the rate $O(N^{1})$.

More generally, if the function $f$ and its derivatives up to order $\nu-1$
are continuous but $f^{(\nu)}$ is discontinuous (i.e., $f$ has a jump of order
$\nu$), then the global convergence rate is $O(N^{-\nu})$. The impact of the
Gibbs phenomenon can lead to inaccurate pricing and hedging when the
approximate option prices are generated via FFT or Fourier series methods at
or around the jumps.

Ruijter et al. (2013) attempt to resolve the Gibbs phenomenon and improve the
algebraic-index convergence rate of the plain COS method by adding filters.
They call this the filter COS method. In their paper, they intensively test
six different filters, e.g., the Fejér filter and the exponential filter, and
discover that the exponential filter gives a better algebraic-index
convergence rate for pricing European options. The method is implemented to
represent the non-smooth PDF of the Variance Gamma (VG) model. Although they
provide some improvements over the current COS method, filtering only works
well away from the peak of the VG PDF; at the peak, the approximation performs
worse. The method also requires almost a thousand COS partial summation terms
to reach acceptable accuracy when the PDF is non-smooth.

According Ortiz-Gracia and Oosterlee (2013, 2016), another disadvantage of
using the COS method is that the COS exhibits periodicity in the vicinity of
the integration boundaries, and long-maturity option round-off errors may
accumulate near the domain boundaries. A series of papers of using wavelets,
e.g., B-spline and Shannon wavelets, have addressed this problem. Wavelet
methods are similar to Fourier series methods; they use wavelets to represent
a PDF and then integrate the discounted payoff function with wavelet basis
functions to reach the closed-form representation of an option pricing
formula. Compared with the COS method, the wavelet methods are flexible and
accurate for pricing options with long maturities. However, the choice of the
wavelet basis functions determines whether the method can achieve a spectral
convergence rate. In some cases, the B-spline and Haar wavelets cannot achieve
spectral convergence (cf. Table 1 in Ortiz-Gracia and Oosterlee, 2013 and
Ortiz-Gracia and Oosterlee, 2016). Although Ortiz-Gracia and Oosterlee use
Shannon wavelets (Shannon Wavelets Inverse Fourier Technique (SWIFT)) to
correct this problem and achieve spectral convergence in some cases where the
PDF is smooth (cf. Table 1 and Table 4 in Ortiz-Gracia and Oosterlee (2016)),
the accuracy of SWIFT heavily relies on a scale parameter. The scale parameter
is determined by the speed of decay of a characteristic function, which in
turn, determines the accuracy of the approximation of the corresponding PDF.
Although, in general, Ortiz-Gracia and Oosterlee (2016) consider 0 or 1 the
best choice of the scale, there is no theoretical reason for this choice, and
different scale parameters are chosen for different PDFs throughout their
paper.

Finally, in our numerical tests, the accuracy of the COS method is very
sensitive to small and large option prices. The COS method tends to be less
accurate when we measure small option prices (cf. Table 8 in Section 8).

To circumvent the aforementioned disadvantages of FFT, Fourier series and
wavelet methods, we propose a novel method, the singular Fourier-Padé (SFP)
method, which exhibits the following characteristics:

  1. 1.

global spectral convergence rate for piecewise continuous PDFs,

  2. 2.

fast error convergence with fewer partial summation terms required,

  3. 3.

accurately prices any European-type option with the features of deep in/out of
the money and very long/short maturities,

  4. 4.

consistently accurate for approximating large or small option prices
throughout, and

  5. 5.

does not require a scale parameter to adjust its accuracy.

Why do we choose the SFP method? Compared with the SFP method, the Fourier-
Padé technique is not a better choice. The Fourier-Padé technique (cf.
Chisholm and Common, 1981; Small and Charron, 1988; Geer, 1995) is a famous
technique for approximating a non-smooth function such that spectral
convergence can be achieved away from the jumps and convergence is not
globally degraded. Nevertheless, Driscoll and Fornberg (2001, 2011) note that
the fundamental limitation of the Fourier-Padé method is that the use of poles
to approximate branch cuts is inefficient, and logarithmic singularities
translated by the jumps on the unit circle are very difficult for the Padé
approximants to simulate. Accordingly, the method fails to converge at the
jump. To address this problem, Driscoll and Fornberg (2001, 2011) add proper
logarithmic branch singularity terms into the Fourier-Padé approximation
process. They call this method the singular Fourier-Padé method. In all
numerical tests (cf. Driscoll and Fornberg, 2001, 2011), the SFP can
accelerate error convergence to the true solutions for analytic periodic
functions or non-smooth functions, and the global effects of the Gibbs
phenomenon due to a jump are largely eliminated. If the function is very
difficult, when jump locations are known in advance, the method still appears
to converge globally and exponentially and offer 4–6 digits of accuracy at the
jumps. Most importantly, in addition to comparing the Fourier-Padé technique,
Driscoll and Fornberg compare the SFP method with other two famous methods of
overcoming the Gibbs phenomenon–the Gegenbauer method (Gottlieb and Shu,
1997), a method to project the Fourier partial sum onto a space spanned by a
Gegenbauer polynomial, and singularity removal (Eckhoff, 1997), a method of
singularity removal for functions with jump discontinuities. In the numerical
tests comparing these three methods (cf. Driscoll and Fornberg, 2001, 2011),
the SFP method yields better global spectral convergence and faster error
convergence than do the other methods. Based on all these benefits of the SFP
method, we chose it over the other available numerical methods.

The remainder of this paper is structured as follows. Section 1 provides an
introduction. Section 2 describes the SFP method. Section 3 introduces the
financial stochastic models we examine in this paper. Section 4 describes and
proves the formulation of the SFP option pricing/Greeks formulae for different
styles of European options. In Section 5, we describe the SFP algorithm and
the Fourier-Padé method to find the jump locations on a non-smooth function.
Section 6 describes the choice of truncated integration intervals. Section 7
provides the error analysis of the SFP method in option pricing and hedging.
Section 8 discusses, analyses and compares the numerical results of the SFP
method with those of other numerical methods. Finally, we conclude and discuss
possible future developments in Section 9.

##  2 Singular Fourier-Padé Interpretation and Correction of the Gibbs
Phenomenon

If we consider a function $f$ with a formal power series representation
$\sum_{k=0}^{\infty}b_{k}x^{k},$ and a rational function defined by
$R_{N,M}=P_{N}/Q_{M},$ where $P_{N}$ and $Q_{M}$ are the polynomials of

$\displaystyle P_{N}(x)=\sum_{n=0}^{N}p_{n}x^{n}\hbox{ and
}Q_{M}(x)=\sum_{m=0}^{M}q_{m}x^{m},$ (1)

respectively, then we say that $R_{N,M}=P_{N}/Q_{M}$ is the (linear) Padé
approximant of order (N, M) of the formal series satisfying the condition

$\displaystyle\left(\sum_{n=0}^{N}p_{n}x^{n}\right)-\left(\sum_{m=0}^{M}q_{m}x^{m}\right)\left(\sum^{M+N}_{k=0}b_{k}x^{k}\right)$
$\displaystyle=\mathcal{O}(x^{N+M+1}).$ (2)

Here, $f$ is approximated by $\sum^{M+N}_{k=0}b_{k}x^{k},$ To obtain the
approximant $R(N,M),$ we simply calculate the coefficients of polynomials
$P_{N}$ and $Q_{M}$ by solving a system of linear equations. To obtain
$\\{q_{m}\\}_{m=0}^{M},$ we first normalise $q_{0}=1$ to ensure that the
system is well-determined and has a unique solution in (2). Then, we consider
the coefficients for $x^{N+1},\ldots,x^{M+N},$ and we can yield a
Toeplitz*333A Toeplitz matrix or diagonal-constant matrix is an invertible
matrix in which each descending diagonal from left to right is constant.
linear system:

$\displaystyle\begin{bmatrix}b_{N+1}&b_{N}&b_{N-1}&\cdots&b_{N+1-M}\\\
b_{N+2}&b_{N+1}&b_{N}&\ddots&b_{N+2-M}\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\
b_{N+M}&\cdots&b_{N+2}&b_{N+1}&b_{N}\end{bmatrix}\begin{bmatrix}q_{0}\\\
q_{1}\\\ \vdots\\\ q_{M}\end{bmatrix}=0.$ (3)

Once $\\{q_{m}\\}_{m=0}^{M}$ is known, $\\{p_{n}\\}_{n=0}^{N}$ is found
through the terms of order N and less in (2). This yields
$\mathbf{p}=B\mathbf{q}$, where $b_{ij}=b_{i-j}$. For example, if $N=M,$ one
obtains

$\displaystyle\begin{bmatrix}p_{0}\\\ p_{1}\\\ \vdots\\\
p_{N}\end{bmatrix}=\begin{bmatrix}b_{0}&&&\\\ b_{1}&b_{0}&&\\\
\vdots&\ddots&\ddots&\\\
b_{N}&\cdots&b_{1}&b_{0}\end{bmatrix}\begin{bmatrix}q_{0}\\\ q_{1}\\\
\vdots\\\ q_{M}\end{bmatrix}.$ (4)

Computing Padé approximants through linear algebra, as we have done here, is
simple but not necessarily the most efficient or stable numerical method (cf.
Driscoll and Fornberg, 2001, 2011).

Now, suppose that $f$ be a piecewise analytic function defined on the interval
$[-\pi,\pi),$ with $s$ jump locations in $f$ at $t=\zeta_{s}\in[\pi,\pi)$,
$s=1,\ldots,S.$ We define a jump as an actual discontinuous point on $f$ or a
discontinuous point appearing after its derivatives. Then the complex Fourier
series (CFS) representation is given by

$\displaystyle
f(t)=\sum_{k=-\infty}^{\infty}b_{k}e^{ikt},\quad\hbox{with}\quad
b_{k}=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(t)e^{-ikt}dt.$ (5)

The transformation $z=e^{it}$, which maps the interval $[\pi,\pi)$ onto the
unit circle in the complex plane, transforms the Fourier series into the
following Laurent series in $z$, which can be split into

$\displaystyle f(z)$
$\displaystyle=\sum_{k=-\infty}^{\infty}b_{k}e^{ikt}=\sum_{k=0}^{\infty}{\vphantom{\sum}}^{\prime}b_{k}z^{k}+\sum_{k=0}^{\infty}{\vphantom{\sum}}^{\prime}b_{-k}z^{-k}$
$\displaystyle=f^{+}(z)+f^{-}(z^{-1}),$ (6)

where the prime sums indicate that the zeroth term should be halved. The
Fourier-Padé approximation of $f_{\pm}$ is comprised of polynomials

$\displaystyle
P_{N}^{\pm}(z)=Q_{M}^{\pm}(z)f^{\pm}(z)+\mathcal{O}(z^{N+M+1}),\quad
z\rightarrow 0.$ (7)

The resulting approximant is then defined as

$\displaystyle{P^{+}_{N}(z)\over Q^{+}_{M}(z)}+{P^{-}_{N}(z^{-1})\over
Q^{-}_{M}(z^{-1})}.$ (8)

However, Driscoll and Fornberg (2001, 2011) note that this approximant does
not reproduce very well at/around the jump locations of the function and make
the approximation inaccurate. Therefore, they suggest that every jump in value
of $f$ at $t=\zeta$ can be attributed to a logarithm of the form

$\displaystyle\log\left(1-{z\over e^{i\zeta}}\right).$ (9)

This logarithmic singularity in $f_{\pm}$, which is difficult for the Padé
approximant to simulate, can be exploited to enhance the approximation
process. This is the rationale behind the SFP method introduced in Driscoll
and Fornberg (2001, 2011). We modify the Fourier-Padé approximant (7) to
obtain the following condition:

$\displaystyle
P_{N}^{\pm}(z)+\mathcal{L}=Q_{M}^{\pm}(z)f^{\pm}(z)+\mathcal{O}(z^{U+1}),$
(10)

where

$\displaystyle\mathcal{L}=\sum_{s=1}^{S}L^{\pm}_{s}(z)\log\left(1-{z\over
e^{i\zeta_{s}}}\right)$ (11)

for some polynomials $L_{s},s=1,\ldots,S$, and $U$ is determined by S and the
degrees of $P_{N}$, $Q_{M}$ and $L_{s}$.

If we extend the SFP method to support any piecewise analytic real function
$f$ in a finite interval $[a,b]$ with a set of jump locations
$\\{\zeta_{s}\\}_{s=1}^{S}\in[a,b]$ appearing in $f$, the CFS representation
of the function is defined as

$\displaystyle
f(x)=\mathfrak{Re}\left[\sum_{k=-\infty}^{\infty}b_{k}e^{i\frac{2\pi}{b-a}kx}\right],\,\hbox{with}\,\,b_{k}=\frac{1}{b-a}\int_{a}^{b}f(x)e^{-i\frac{2\pi}{b-a}kx}dx.$
(12)

Here, $\mathfrak{Re}$ represents the real part of the function. As we focus on
approximating a real function, we can further obtain

$\displaystyle
f(x)=\mathfrak{Re}\left[2\sum_{k=1}^{\infty}b_{k}e^{i\frac{2\pi}{b-a}kx}+b_{0}\right].$
(13)

Based on this representation, we denote $z$ as
$\exp\left({i\frac{2\pi}{b-a}x}\right)$ and then transform $f$ into a
truncated power series of $f_{1}$ equal to

$\displaystyle\mathfrak{Re}\left[2\sum_{k=1}^{U}b_{k}z^{k}+b_{0}\right].$ (14)

The transformation $z=\exp\left({i\frac{2\pi}{b-a}x}\right)$ also suggests
that the jump location $\zeta$ translates into

$\displaystyle\log\left(1-{z\over\varepsilon}\right),\hbox{ where
}\varepsilon=e^{i\frac{2\pi}{b-a}\zeta}$ (15)

in $f_{1}.$ Finally, applying (14) and (15), we can reach the final form of
the SFP approximant given by

$\displaystyle
P_{N}^{+}(z)+\sum_{s=1}^{S}L_{N_{s}}^{+}(z)\log\left(1-z/\varepsilon_{s}\right)=\left(2\sum_{k=1}^{U}b_{k}z^{k}+b_{0}\right)Q_{M}^{+}(z)+\mathcal{O}(z^{U+1}),$
(16)

where

$\displaystyle\begin{array}[]{rclrcl}P_{N}^{+}(z)&=&\sum_{n=0}^{N}p_{n}z^{n},&Q_{M}^{+}(z)&=&\sum_{m=0}^{M}q_{m}z^{m}\neq
0,\\\
L_{N_{s}}^{+}(z)&=&\sum_{n_{s}=0}^{N_{s}}l_{n_{s}}z^{n_{s}},&s&=&1,\ldots,S,\\\
U&=&N+M+S+\sum_{s=1}^{S}N_{s}.&&&\end{array}$ (20)

Once we have determined the unknown coefficients of $\\{p_{n}\\}_{n=0}^{N},$
$\\{q_{m}\\}_{m=0}^{M}$ and $\\{l_{n_{s}}\\}_{n_{s}=0}^{N_{s}}$ (cf. Section
5), the SFP representation of $f(x)$ can be formulated as

$\displaystyle
f(x)=\mathfrak{Re}\left({P_{N}^{+}(z)+\sum_{s=1}^{S}L_{N_{s}}^{+}(z)\log\left(1-z/\varepsilon_{s}\right)\over
Q_{M}^{+}(z)}\right),\hbox{ with }z=e^{i\frac{2\pi}{b-a}x}.$ (21)

##  3 Stochastic Models of the Asset Dynamics

We assume frictionless financial equity markets and no arbitrage, and take as
given an equivalent martingale measure $\mathds{Q}$ chosen by the market. All
stochastic processes defined in the following are assumed to live on the
complete filtered probability space
$(\Omega,\mathcal{F},\\{\mathcal{F}_{t}\\}_{t\geq 0},\mathds{Q}).$ Standard
references for exponential Lévy processes can be found in Schoutens (2003) and
Cont and Tankov (2004); for affine processes, see Duffie et al. (2003).

###  3.1 Exponential Lévy Processes

The stock price process $(S_{t})_{t\geq 0}$ under $\mathds{Q}$ driven by an
exponential Lévy process can be defined:

$\displaystyle S_{T}$
$\displaystyle=S_{t}e^{(r-q)(T-t)+X_{T}-X_{t}+\omega(T-t)}$ (22)
$\displaystyle=S_{t}e^{(r-q+\omega)(T-t)+X_{T-t}},$ (23)

where $X_{T},$ $X_{t}$ and $X_{T}-X_{t}$ are all Lévy processes. As Lévy
processes have independent stationary increments, we can say that
$X_{T}-X_{t}=X_{T-t}.$ Throughout the paper, $r\geq 0$ and $q\geq 0$ denote
the constant risk-free interest rate and the constant dividend yield,
respectively; $S_{t}$ represents the known stock price at time $t$; and
$S_{T}$ represents the random stock price at time $T$. The condition that
$(S_{T}e^{-(r-q)(T-t)})_{t\geq 0}$ is a martingale will be guaranteed by an
appropriate choice of the mean-correcting compensator $\omega$ as follows:

$\displaystyle\omega={1\over T-t}\mathbb{E}\left(e^{X_{T-t}}\right),$ (24)

where $\mathbb{E}\left(e^{X_{T-t}}\right)$ is assumed to be finite for all
$0\leq t\leq T$. Given a Lévy process $(X_{T-t})_{t\geq 0}$, define the
corresponding characteristic function as follows:

$\displaystyle\varphi(u)=\mathbb{E}[e^{iu(X_{T-t})}]=e^{(T-t)\phi(u)},\quad
u\in\mathbb{R}.$ (25)

Here, $\phi(u)$ is a continuous function with the Lévy-Khintchine
representation given by

$\phi(u)={1\over 2}A^{2}u-i\gamma
u+\int_{-\infty}^{\infty}(1-e^{iu\chi}+iu\chi\mathds{1}_{|\chi|\geq
1})\nu(\mathrm{d}\chi),\quad\chi\in X_{T-t},$

with characteristic triplet $(A,\gamma,\nu).$ There is substantial
consideration of exponential Lévy processes in stock process modelling. Due to
space limitations, we cannot show all the processes in this paper; rather, we
selectively investigate the processes due to their popularity or their
numerical implementation challenge. However, this does not mean that our
method only works on the processes we selected. As the paper evolves, we can
see that as long as a process has a characteristic function, the option price
and its risk driven by the process can be properly approximated via our
method.

###  3.2 Exponential Affine Processes

An exponential affine characteristic function is given by

$\displaystyle\varphi(u)=\mathbb{E}[e^{iuX_{t}}]=\exp(C(u,t)+X_{0}D(u,t)).$
(26)

Here, the function $C(u,t)$ is fully characterised by $dC(u,t)/dt=D(u,t)$ and
the function $D(u,t)$ satisfies a Riccati equation (Duffie et al., 2003). More
precisely, we can find functions $C(u,t)$ and $D(u,t)$ with $C(u,t)=0$ and
$D(u,t)=u,$ such that

$\displaystyle M_{t}=\exp(C(u,T-t)+X_{t}D(u,T-t))$ (27)

is a martingale process. In the family of affine processes, we focus on
investigating the Heston model in this paper. The stochastic differential
equation (SDE) of the Heston model is written as

$\displaystyle dS_{t}$ $\displaystyle=(r-q)S_{t}dt+\sqrt{y_{t}}S_{t}dW_{1,t},$
(28) $\displaystyle dy_{t}$
$\displaystyle=\lambda(\bar{y}-y_{t})dt+\eta\sqrt{y_{t}}dW_{2,t},$ (29)

where $L_{t}$ and $y_{t}$ denote the stochastic log-asset price variable and
the variance of the asset price process, respectively. In this process, the
speed of mean reversion $\lambda$, the mean level of variance $\bar{y}$ and
the volatility of volatility $\eta$ are constant values greater than or equal
to zero. Additionally, the Brownian motions $W_{1,t}$ and $W_{2,t}$ are
correlated with the correlation coefficient $\rho_{s}$, and
$\omega=\mathbb{E}[e^{i(-1i)X_{t}}]$ is the mean-correcting compensator. The
model characteristic function fits in the general affine characteristic
function framework and is given by

$\displaystyle\varphi(u)$
$\displaystyle=\exp\bigg{(}iu\left((r-q)t+\omega\right)+\frac{y_{0}}{\eta^{2}}\left(\frac{1-e^{Et}}{1-Fe^{Et}}(\lambda-i\rho_{s}\eta
u-E)\right)+$
$\displaystyle\frac{\lambda\bar{y}}{\eta^{2}}\left(t(\lambda-i\rho_{s}\eta
u-E)-2\log(\frac{1-Fe^{-Et}}{1-F})\right)\bigg{)},\,\,u\in\mathbb{R},$ (30)
with $\displaystyle E=$ $\displaystyle\sqrt{\left(\lambda-i\rho_{s}\eta
u\right)+(u^{2}+iz)\eta^{2}},\,$ $\displaystyle F=$
$\displaystyle(\lambda-i\rho_{s}\eta u-E)/(\lambda-i\rho_{s}\eta u+E).$

This characteristic function is uniquely specified because we take
$\sqrt{(x+yi)}$ such that its real part is nonnegative, and we restrict the
complex logarithm to its principal branch.

In this case, as Lord and Kahl (2010) prove, the resulting characteristic
function is the correct one for all complex numbers $z$ in the analytic strip
of the characteristic function. In the SDE, we have two possible conditions
with respect to $\lambda,$ $\bar{y}$ and $\eta$:

$\displaystyle 2\lambda\bar{y}\geq\eta^{2},$ (31) $\displaystyle
2\lambda\bar{y}<\eta^{2}.$ (32)

The model satisfies the Feller property if (31) holds; otherwise, (32) holds.
If a process fulfils the property, the process never hits zero; conversely, if
it does not, the process can reach 0. Condition (32) is a very important
property for the Heston SDEs because they can only have a unique solution when
we specify a boundary condition at 0. In mathematical finance, the chosen
boundary condition is that the process remains at 0. We define this as an
absorbing boundary condition. When the process reaches 0 and is allowed to
leave 0, we call it a reflecting boundary. These two boundary conditions are
crucial for pricing early-exercise options, including American options and
barrier options.

##  4 Singular Fourier-Padé Representation of European Option Prices and their
Option Greeks

###  4.1 Option Pricing Formulae

In this section, we derive closed-form formulae for European-style options
using the SFP method. Considering the PDF $f$ of a stochastic process, the
current log-price $x:=\log S,$ the strike price of $K$ and maturity $T\geq t,$
we can express the option price $V(x,K,t)$ starting at time $t$ with its
contingent claim paying out $G(S_{T},K)$ as follows:

$\displaystyle V(x,K,t)$
$\displaystyle=e^{-r(T-t)}\mathbb{E}(G(S_{T},K)|S_{t}=e^{x})$ (33)
$\displaystyle=e^{-r(T-t)}\mathbb{E}(G(S_{t}e^{X_{T}-X_{t}},K))$ (34)
$\displaystyle=e^{-r(T-t)}\bigintssss_{-\infty}^{\infty}G(e^{x+\chi},K)f(z)dz,\quad\chi\in
X_{T}-X_{t}.$ (35)

Furthermore, if we choose an interval $[c,d]$ satisfying the condition

$\displaystyle\int_{c}^{d}f(\chi)e^{iu\chi}\rm{d}\chi\approx\mathbb{E}[e^{iu(X_{T}-X_{t})}]=\varphi(u),$
(36)

where $\varphi(u)$ is a characteristic function of $X_{T}-X_{t},$ we can
approximate the pricing formula:

$\displaystyle V(x,K,t)\approx
e^{-r(T-t)}\int_{c}^{d}G(e^{x+\chi},K)f(\chi)\mathrm{d}\chi.$ (37)

######  Theorem 1.

When a dividend-paying risky asset price process $(S_{t})_{t\geq 0}$ with a
traceable, analytical characteristic function $\varphi(\cdot)$ has a current
asset price $e^{x}=S,$ risk-free interest rate $r$ and compounded continuous
dividend $q,$ the SFP pricing formula of a European vanilla call option driven
by this process with maturity $T$ and strike price $K$ is

$\displaystyle V(x,K,t)$
$\displaystyle=e^{-r(T-t)}\mathfrak{Re}\left[{P_{N}^{+}(z)+\sum_{s=1}^{S}L^{+}_{N_{s}}(z)\log\left(1-z/\varepsilon_{s}\right)\over
Q^{+}_{M}(z)}\right],$ (38)

with

$\displaystyle\begin{array}[]{rclrcl}z&=&e^{i\frac{2\pi}{d-c}(-x+\log
K)}&\varepsilon_{s}&=&e^{i\frac{2\pi}{d-c}\zeta_{s}}\\\
P_{N}^{+}(z)&=&\sum_{n=0}^{N}p_{n}z^{n},&Q_{M}^{+}(z)&=&\sum_{m=0}^{M}q_{m}z^{m}\neq
0,\\\
L_{N_{s}}^{+}(z)&=&\sum_{n_{s}=0}^{N_{s}}l_{n_{s}}z^{n_{s}},&s&=&1,\ldots,S.\\\
\end{array}$ (42)

Here, $\zeta_{s}$ is the jump in $V(x,K,t).$

Proof: First, the payoff function of a European call option is
$G(e^{x+\chi},K)=\max\left(e^{x+\chi}-K,0\right)$ in (37). To conform to the
SFP approximation framework, we first transform the payoff:

$\displaystyle\max\left(e^{x+\chi}-K,0\right)=K\max\left(e^{x+\chi-\log
K}-1,0\right).$ (43)

Accordingly, by replacing $x+\chi-\log K$ with $y$, we have a new form of
$V(x,K,t)$ denoted as

$\displaystyle
e^{-r(T-t)}\int_{-\infty}^{\infty}K\max(e^{y}-1,0)f\left(y-x+\log
K\right)\mathrm{d}y.$ (44)

As we intend to make the SFP method more efficient, we define a truncated
computational interval $[c,d]$ (cf. Section 6), which satisfies condition
$(\ref{eqn:charAppx}),$ to replace $[-\infty,\infty].$ Then, $V(x,K,t)$ is
reformulated as

$\displaystyle e^{-r(T-t)}\int_{c}^{d}K\max(e^{y}-1,0)f\left(y-x+\log
K\right)\mathrm{d}y$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\int_{0}^{d}K(e^{y}-1)f\left(y-x+\log K\right)\mathrm{d}y.$ (45)

Using the Fourier transform shift theorem and the CFS expansion shown in (13),
we express $f\left(y-x+\log K\right)$ as

$\displaystyle\mathfrak{Re}\left[2\sum_{k=1}^{\infty}b_{k}e^{i\frac{2\pi}{b-a}ky}+b_{0}\right],$
(46)

where

$\displaystyle
b_{k}=\frac{1}{d-c}\int_{c}^{d}f(y)e^{-i\frac{2\pi}{d-c}ky}\mathrm{d}y\left(e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}\right)\quad\hbox{and}\quad
b_{0}=\frac{1}{d-c}\int_{c}^{d}f(y)\mathrm{d}y.$ (47)

We substitute (46) into (4.1) and apply Fubini’s theorem; $V(x,K,t)$ can be
computed as

$\displaystyle\quad
e^{-r(T-t)}\int_{0}^{d}K(e^{y}-1)\mathfrak{Re}\left[2\sum_{k=1}^{\infty}b_{k}e^{i\frac{2\pi}{d-c}ky}+b_{0}\right]\mathrm{d}y$
$\displaystyle=e^{-r(T-t)}\mathfrak{Re}\Bigg{[}2\sum_{k=1}^{\infty}\frac{1}{d-c}\int_{c}^{d}f(y)e^{-i\frac{2\pi}{d-c}ky}\mathrm{d}y\left(e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}\right)\int_{0}^{d}K(e^{y}-1)e^{i\frac{2\pi}{d-c}ky}\mathrm{d}y+$
$\displaystyle\quad\frac{1}{d-c}\int_{c}^{d}f(y)\mathrm{d}y\int_{0}^{d}K(e^{y}-1)\mathrm{d}y\Bigg{)}.$
(48)

In the equation above, basic calculus implies that

$\displaystyle\int_{0}^{d}K(e^{y}-1)e^{i\frac{2\pi}{d-c}ky}\mathrm{d}y$
$\displaystyle=\frac{K(d-c)}{i2\pi
k+(d-c)}\left(e^{(i\frac{2\pi}{d-c}k+1)d}-1\right)-\frac{K(d-c)}{i2\pi
k}\left(e^{(i\frac{2\pi}{d-c}k)d}-1\right),$ (49)
$\displaystyle\int_{0}^{d}K(e^{y}-1)\mathrm{d}y$
$\displaystyle=K(e^{d}-1)-Kd.$ (50)

In the meantime, because of condition (36), we can also see that

$\displaystyle\int_{c}^{d}f(y)e^{-i\frac{2\pi}{d-c}ky}\mathrm{d}y\approx\varphi\left(-\frac{2\pi}{d-c}k\right)\quad\hbox{and}\quad\int_{c}^{d}f(y)\mathrm{d}y\approx\varphi(0)=1.$
(51)

For the sake of simplicity in (4.1), we set

$\displaystyle\begin{array}[]{rclrcl}\widehat{B}_{k}&=&\frac{1}{d-c}\varphi\left(-\frac{2\pi}{d-c}k\right),&\widehat{B}_{0}&=&{1\over
d-c}\varphi(0)={1\over d-c},\\\ \widehat{G}_{k}&=&\frac{K(d-c)}{i2\pi
k+(d-c)}\left(e^{(i\frac{2\pi}{d-c}k+1)d}-1\right)-\frac{K(d-c)}{i2\pi
k}\left(e^{(i\frac{2\pi}{d-c}k)d}-1\right),&\widehat{G}_{0}&=&K(e^{d}-1)-Kd\\\
\end{array}$ (54)

to obtain a simplified form:

$\displaystyle V(x,K,t)$
$\displaystyle=e^{-r(T-t)}\mathfrak{Re}\left[2\sum_{k=1}^{\infty}\widehat{B}_{k}\widehat{G}_{k}e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}+\widehat{B}_{0}\widehat{G}_{0}\right].$ (55)

To express our final pricing formula with the SFP representation, we
approximate

$\displaystyle
2\sum_{k=1}^{\infty}\widehat{B}_{k}\widehat{G}_{k}e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}+\widehat{B}_{0}\widehat{G}_{0}$ (56)

in (55) with (16). To do so, we first replace $-x+\log K$ with a variable
$y_{1}$ in (56), and then, this allows us to set
$\exp\left({i\frac{2\pi}{d-c}y_{1}}\right)$ equal to $z.$ The transformation
$z=\exp\left({i\frac{2\pi}{d-c}y_{1}}\right)$ maps the interval $[c,d]$ onto
the unit circle in $z.$ This change also transforms the jumps $\zeta$ along
$V(x,K,t)$ into $z$ with the form of
$\varepsilon=\exp\left(i\frac{2\pi}{d-c}\zeta\right).$ Finally, expressing
(56) with a new variable of $z,$ we have

$\displaystyle
2\sum_{k=1}^{\infty}\widehat{B}_{k}\widehat{G}_{k}z^{k}+\widehat{B}_{0}\widehat{G}_{0}.$
(57)

Substituting the equation above with $f_{1}(z)$ in (16), we obtain the
approximant given by

$\displaystyle
P_{N}^{+}(z)+\sum_{s=1}^{S}L^{+}_{N_{s}}(z)\log\left(1-z/\varepsilon_{s}\right)=\left(2\sum_{k=1}^{U}\widehat{B}_{k}\widehat{G}_{k}z^{k}+\widehat{B}_{0}\widehat{G}_{0}\right)Q^{+}_{M}(z)+\mathcal{O}(z^{U+1})$
(58)
$\displaystyle\begin{array}[]{rclrcl}P_{N}^{+}(z)&=&\sum_{n=0}^{N}p_{n}z^{n},&Q_{M}^{+}(z)&=&\sum_{m=0}^{M}q_{m}z^{m}\neq
0,\\\
L_{N_{s}}^{+}(z)&=&\sum_{n_{s}=0}^{N_{s}}l_{n_{s}}z^{n_{s}},&s&=&1,\ldots,S,\\\
\varepsilon_{s}&=&e^{i\frac{2\pi}{d-c}\zeta_{s}},&U&=&N+M+\sum_{s=1}^{S}N_{s}.\end{array}$
(62)

Once we can determine the unknown coefficients of $\\{p_{n}\\}_{n=0}^{N},$
$\\{q_{m}\\}_{m=0}^{M}$ and $\\{l_{n_{s}}\\}_{n_{s}=0}^{N_{s}}$ in (58) via
the algorithm shown in Section 5 and replace

$2\sum_{k=1}^{\infty}\widehat{B}_{k}\widehat{G}_{k}e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}+\widehat{B}_{0}\widehat{G}_{0}$

with

${P_{N}^{+}(z)+\sum_{s=1}^{S}L^{+}_{N_{s}}(z)\log\left(1-z/\varepsilon_{s}\right)\over
Q^{+}_{M}(z)},\quad z=\exp\left({i\frac{2\pi}{d-c}(-x+\log K)}\right)$

in (55), we reach our first SFP pricing formula for European vanilla call
options, as in (38). Q.E.D.

We can apply the same technique to seek the SFP pricing formula of vanilla put
options. To achieve this, we first transform the put payoff function
$G(e^{x+\chi},K)=\max\left(K-e^{x+\chi},0\right)$ into

$\displaystyle-K\max\left(e^{x+\chi-\log K}-1,0\right).$ (63)

Then, we follow the proof of Theorem 1 to reach the CFS representation of
$V(x,K,t)$ given by

$\displaystyle
e^{-r(T-t)}\mathfrak{Re}\left[2\sum_{k=1}^{\infty}\widehat{B}_{k}\widehat{G}_{k}e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}+\widehat{B}_{0}\widehat{G}_{0}\right],$ (64)

with

$\displaystyle\widehat{G}_{k}$
$\displaystyle=\int_{c}^{0}-K(e^{y}-1)e^{i\frac{2\pi}{d-c}ky}\mathrm{d}y=\frac{K(d-c)}{i2\pi
k+(d-c)}\left(e^{(i\frac{2\pi}{d-c}k+1)c}-1\right)-\frac{K(d-c)}{i2\pi
k}\left(e^{(i\frac{2\pi}{d-c}k)c}-1\right),$ (65)
$\displaystyle\widehat{G}_{0}$
$\displaystyle=\int_{c}^{0}-K(e^{y}-1)\mathrm{d}y=K(e^{c}-1)-Kc.$ (66)

In the same manner shown in the proof of Theorem 1, we approximate the CFS
expansion in (64) with the SFP approximant (16). We can therefore yield the
same expression of (38) with different coefficient values of
$\\{p_{n}\\}_{n=0}^{N},$ $\\{q_{m}\\}_{m=0}^{M}$ and
$\\{l_{n_{s}}\\}_{n_{s}=0}^{N_{s}}.$

From the vanilla call and put pricing formulae, we can see that they share the
same format as the SFP approximant of (38). It is no different for other
European-type options discussed in this paper. The only difference is that the
values of $\\{p_{n}\\}_{n=0}^{N},$ $\\{q_{m}\\}_{m=0}^{M}$ and
$\\{l_{n_{s}}\\}_{n_{s}=0}^{N_{s}}$ in the SFP formula of each option are
completely different from one another. This is attributed to the jumps $\zeta$
in each option and their payoff functions, which lead to different Fourier
transformations of $\widehat{G}_{k}$ and $\widehat{G}_{0}$ in (55). For the
reader’s information, we enumerate all the options we investigate in this
paper and their corresponding payoff functions $G(e^{x+\chi},K)$ and Fourier
transformations $G_{k}$ and $G_{0}$ in Tables 1 and 2.

###  4.2 Option Greeks

Now, we turn our attention to deriving the option Greeks. Although accurately
valuing financial claims plays a key role in financial modelling, the risk
management (hedging) of these derivative instruments is equally important.
Financial institutions manage option risk when they sell options to their
clients through the analysis of the Greeks, which are defined as the
sensitivities of the option price to single-unit changes in the values of
state variables or parameters. Such sensitivities represent the different
dimensions of the risk associated with an option. In this paper, we will focus
on deriving three option Greeks—Delta, Gamma, and Vega. Delta, $\Delta,$ is
defined as the rate of change in the option value with respect to changes in
the underlying asset price; Gamma, $\Gamma$ is the rate of change of $\Delta$
with respect to changes in the underlying price; and finally, Vega is the
measurement of an option’s sensitivity to changes in the volatility of the
underlying asset price. In general, volatility measures the amount and speed
at which the price moves up and down, and it is often based on changes in the
recent, historical prices of a trading instrument. Other Greeks, such as
Theta, can be derived in a similar fashion; however, depending on the
characteristic function, the derivation expression might be rather lengthy. We
omit them here, as many terms are repeated.

Delta is the first derivative of the value of $V$ of the option with respect
to the underlying instrument price S. Hence, differentiating the CFS expansion
of $V$ (55) with respect to $S,$ we have

$\displaystyle\Delta_{t}={\partial V(x,K,t)\over\partial S}$
$\displaystyle={\partial V(x,K,t)\over\partial x}{\partial x\over\partial S}$
$\displaystyle=e^{-r(T-t)-x}\Bigg{(}\mathfrak{Re}\Bigg{[}2\sum_{k=1}^{\infty}\left(-i\frac{2\pi}{d-c}k\right)\widehat{B}_{k}\widehat{G}_{k}e^{i\frac{2\pi}{d-c}k(-x+\log
K)}\Bigg{]}\Bigg{)}.$ (67)

In a similar fashion, we can obtain $\Gamma_{t}$ by differentiating
$\Delta_{t}$ with respect to $S$ such that

$\displaystyle\Gamma_{t}={\partial^{2}V(x,K,t)\over\partial
S^{2}}={\partial\Delta_{t}\over\partial S}={\partial\Delta_{t}\over\partial
x}{\partial x\over\partial S},$ (68)

and eventually,

$\displaystyle\Gamma_{t}$
$\displaystyle=e^{-r(T-t)-2x}\Bigg{(}\mathfrak{Re}\Bigg{[}2\sum_{k=1}^{\infty}\left(i\frac{2\pi}{d-c}k\right)\left(i\frac{2\pi}{d-c}k+1\right)\widehat{B}_{k}\widehat{G}_{k}e^{i\frac{2\pi}{d-c}k(-x+\log
K)}\Bigg{]}\Bigg{)}.$ (69)

It is also easy to obtain the formula for Vega, ${\partial V\over\partial
y_{t}},$ where $y_{t}$ is the initial value of the volatility at time $t.$ For
example, for the Heston model, as $y_{0}$ is the initial value of the
volatility in (3.2), we derive Vega as follows:

$\displaystyle{\partial V(x,K,y_{0},t)\over\partial y_{0}}$
$\displaystyle=e^{-r(T-t)}\Bigg{(}\mathfrak{Re}\Bigg{[}2\sum_{k=1}^{\infty}{\partial\widehat{B}_{k}\over\partial
y_{0}}\widehat{G}_{k}e^{i\frac{2\pi}{d-c}k(-x+\log K)}\Bigg{]}\Bigg{)},$ (70)

with

$\displaystyle{\partial\widehat{B}_{k}\over\partial
y_{0}}={\partial\varphi(-\frac{2\pi}{d-c}k,y_{0})\over\partial y_{0}},$ (71)

where $\varphi$ contains the parameter $y_{0}.$

To obtain our first SFP representation of $\Delta,$ we first let
$y_{1}=-x+\log K,$ $z=\exp\left(i{2\pi\over d-c}y_{1}\right)$ and then
transform all the jumps $\zeta$ in $\Delta_{t}$ into
$\varepsilon=\exp\left(i{2\pi\over d-c}\zeta\right)$ in (4.2). Accordingly,
this transforms the CFS representation in (4.2) into the form

$\displaystyle
f_{1}(z)=2\sum_{k=1}^{U}\left(-i\frac{2\pi}{d-c}k\right)\widehat{B}_{k}\widehat{G}_{k}z^{k}.$
(72)

Based on the equation above, using (16), we can eventually obtain the SFP
approximant given by

$\displaystyle
P_{N}^{+}(z)+\sum_{s=1}^{S}L^{+}_{N_{s}}(z)\log\left(1-z/\varepsilon_{s}\right)=f_{1}(z)Q^{+}_{M}(z)+\mathcal{O}(z^{U+1}).$
(73)

Applying the approximation algorithm in Section 5 to determine the
coefficients of $P_{N}^{+},$ $Q^{+}_{M},$ and $L^{+}_{N_{s}},$ we can obtain
the SPF formula for $\Delta_{t}$ with the form

$\displaystyle
e^{-r(T-t)-x}\mathfrak{Re}\left({P_{N}^{+}(z)+\sum_{s=1}^{S}L^{+}_{N_{s}}(z)\log\left(1-z/\varepsilon_{s}\right)\over
Q^{+}_{M}(z)}\right).$ (74)

To determine the SFP approximant of $\Gamma_{t}$ and Vega, we follow the same
idea of approximating $\Delta_{t}$ but replace $f_{1}(z)$ for

$\displaystyle
2\sum_{k=1}^{U}\left(i\frac{2\pi}{d-c}k\right)\left(i\frac{2\pi}{d-c}k+1\right)\widehat{B}_{k}\widehat{G}_{k}z^{k}\quad\mbox{and}\quad
2\sum_{k=1}^{\infty}{\partial\widehat{B}_{k}\over\partial
y_{0}}\widehat{G}_{k}z^{k}.$ (75)

Table 1: Payoff functions and their transforms for a variety of financial
contingency claims

Financial Contingency Claim | Payoff Function | Transformed Payoff Function
---|---|---
| $G(S_{T},K)$ | $G(e^{x+\chi},K)$
Call | $\max(S_{T}-K,0)$ | $K\max(e^{x+\chi-\log K}-1,0)$
Put | $\max(K-S_{T},0)$ | $-K\max(e^{x+\chi-\log K}-1,0)$
Covered Call | $\min(S_{T},K)$ | $K\min(e^{x+\chi-\log K}-1,0)+K$
Cash-or-Nothing Call | $\mathds{1}_{S_{T}\geq K}$ | $\mathds{1}_{e^{x+\chi-\log K}\geq 1}$
Cash-or-Nothing Put | $\mathds{1}_{S_{T}\leq K}$ | $\mathds{1}_{e^{x+\chi-\log K}\leq 1}$
Asset-or-Nothing Call | $S_{T}\mathds{1}_{S_{T}\geq K}$ | $e^{x+\chi}\mathds{1}_{e^{x+\chi-\log K}\geq 1}$
Asset-or-Nothing Put | $S_{T}\mathds{1}_{S_{T}\leq K}$ | $e^{x+\chi}\mathds{1}_{e^{x+\chi-\log K}\leq 1}$
Asymmetric Call | $(S_{T}^{n}-K^{n})\mathds{1}_{S_{T}\geq K}$ | $K^{n}(e^{n(x+\chi-\log K)}-1)\mathds{1}_{e^{x+\chi-\log K}\geq 1}$
Asymmetric Put | $(K^{n}-S_{T}^{n})\mathds{1}_{S_{T}\leq K}$ | $-K^{n}(e^{n(x+\chi-\log K)}-1)\mathds{1}_{e^{x+\chi-\log K}\leq 1}$
Symmetric Call | $(S_{T}-K)^{n}\mathds{1}_{S_{T}\geq K}$ | $\sum\limits_{j=0}^{n}\binom{n}{j}(-1)^{(n-j)}e^{j(x+\chi)+(n-j)\log K}\mathds{1}_{e^{x+\chi-\log K}\geq 1}$
Symmetric Put | $(K-S_{T})^{n}\mathds{1}_{S_{T}\leq K}$ | $\sum\limits_{j=0}^{n}\binom{n}{j}(-1)^{(n-j)}e^{(n-j)(x+\chi)+j\log K}\mathds{1}_{e^{x+\chi-\log K}\leq 1}$

Table 2: Complex Fourier transforms for a variety of financial contingency
claims

Financial Contingency Claim | Fourier Transform | Fourier Transform
---|---|---
| $\widehat{G}_{k}$ | $\widehat{G}_{0}$
Call | $\frac{K(d-c)}{i2\pi k+(d-c)}\left(e^{(i\frac{2\pi}{d-c}k+1)d}-1\right)-\frac{K(d-c)}{i2\pi k}\left(e^{(i\frac{2\pi}{d-c}k)d}-1\right)$ | $K(e^{d}-1)-Kd$
Put | $\frac{K(d-c)}{i2\pi k+(d-c)}\left(e^{(i\frac{2\pi}{d-c}k+1)c}-1\right)-\frac{K(d-c)}{i2\pi k}\left(e^{(i\frac{2\pi}{d-c}k)c}-1\right)$ | $K(e^{c}-1)-Kc$
Covered Call | $-\frac{K(d-c)}{i2\pi k+(d-c)}\left(e^{(i\frac{2\pi}{d-c}k+1)c}-1\right)+\frac{K(d-c)}{i2\pi k}\left(e^{(i\frac{2\pi}{d-c}k)c}-1\right)$ | $-K(e^{c}-1)+Kc$
Cash-or-Nothing Call | $\frac{(d-c)}{i2\pi k}\left(e^{(i\frac{2\pi}{d-c}k)d}-1\right)$ | $d$
Cash-or-Nothing Put | $-\frac{(d-c)}{i2\pi k}\left(e^{(i\frac{2\pi}{d-c}k)c}-1\right)$ | $-c$
Asset-or-Nothing Call | $\frac{(d-c)}{i2\pi k+(d-c)}\left(e^{(i\frac{2\pi}{b-a}k+1)d}-1\right)$ | $(e^{d}-1)$
Asset-or-Nothing Put | $-\frac{(d-c)}{i2\pi k+(d-c)}\left(e^{(i\frac{2\pi}{d-c}k+1)c}-1\right)$ | $-(e^{c}-1)$
Asymmetric Call | $\frac{K^{n}(d-c)}{i2\pi k+n(d-c)}\left(e^{(i\frac{2\pi}{d-c}k+n)d}-1\right)-\frac{K^{n}(d-c)}{i2\pi k}\left(e^{(i\frac{2\pi}{d-c}k)d}-1\right)$ | ${K^{n}\over n}(e^{nd}-1)-K^{n}d$
Asymmetric Put | $\frac{K^{n}(d-c)}{i2\pi k+n(d-c)}\left(e^{(i\frac{2\pi}{d-c}k+1)c}-1\right)-\frac{K^{n}(d-c)}{i2\pi k}\left(e^{(i\frac{2\pi}{d-c}k)c}-1\right)$ | ${K^{n}\over n}(e^{nc}-1)-K^{n}c$
Symmetric Call | $\sum\limits_{j=0}^{n}\binom{n}{j}(-1)^{(n-j)}\frac{K^{n}(d-c)}{i2\pi k+j(d-c)}\left(e^{(i\frac{2\pi}{d-c}k+j)d}-1\right)$ | $\sum\limits_{j=1}^{n}\binom{n}{j}(-1)^{(n-j)}\times$
|  | ${K^{n}\over j}(e^{jd}-1)+(-1)^{n}K^{n}d$
Symmetric Put | $\sum\limits_{j=0}^{n}\binom{n}{j}(-1)^{(n-j+1)}\frac{K^{n}(d-c)}{i2\pi k+(n-j)(d-c)}\left(e^{(i\frac{2\pi}{d-c}k+n-j)c}-1\right)$ | $\sum\limits_{j=0}^{n-1}\binom{n}{j}(-1)^{(n-j+1)}\times$
|  | ${K^{n}\over n-j}(e^{(n-j)c}-1)-K^{n}c$

##  5 Singular Fourier-Padé Algorithm and Locating Singularities

The approach to computing the polynomial coefficients needed in the SFP method
is fairly straightforward. To demonstrate the algorithm, we focus on a simple
case where the option pricing and Greeks formulae are infinitely smooth apart
from jumps located at the endpoints $c$ and $d.$ As we consider
$z=\exp\left(i{2\pi\over d-c}y_{1}\right)$ in either the option pricing
formula or the Greeks formula, the jump of $c$ and $d$ in the z-plane is -1.
For sake of simplicity, we denote $f_{1}(z)$ as the CFS representation of any
European-style pricing formula or its option Greeks one. With some
superscripts dropped for clarity and knowing that $s=1,$ in (16), we have

$\displaystyle
P_{N}(z)+L_{N_{1}}(z)\log\left(1-{z\over\varepsilon_{1}}\right)=f_{1}(z)Q_{M}(z)+\mathcal{O}(z^{U+1}),$
(76)

where $N+M+N_{1}=U.$ Both $L_{N_{1}}$ and $f_{1}(z)$ have Taylor series and
CFS expansions, respectively, to determine U; therefore, their expansions are

$\displaystyle\log\left(1-{z\over\varepsilon_{s}}\right)$
$\displaystyle=\sum_{k=1}^{U}-{z^{k}\over\varepsilon_{1}^{k}}+0$ (77)
$\displaystyle f_{1}(z)$
$\displaystyle=2\sum_{k=1}^{U}\widehat{B}_{k}\widehat{G}_{k}z^{k}+\widehat{B}_{0}\widehat{G}_{0}.$
(78)

Our goal is to derive a linear system for the unknown polynomial coefficients.
Note that $Q_{M}(z)$ and $L_{N_{1}}(z)$ are determined only by terms of order
greater than $N$. Accordingly, we seek a linear solution to

$\displaystyle\begin{bmatrix}\widehat{B}\widehat{G}&-L\\\
\end{bmatrix}\begin{bmatrix}\mathbf{q}\\\ \mathbf{l}\end{bmatrix}=\mathbf{0}.$
(79)

Here, $\widehat{B}\widehat{G}$ is the $(M+N_{1}+1)\times(M+1)$ Toeplitz matrix

$\displaystyle\begin{bmatrix}\widehat{B}_{\frac{U}{2}+1}\widehat{G}_{\frac{U}{2}+1}&\widehat{B}_{\frac{U}{2}}\widehat{G}_{\frac{U}{2}}&\cdots&\widehat{B}_{1}\widehat{G}_{1}\\\
\widehat{B}_{\frac{U}{2}+2}\widehat{G}_{\frac{U}{2}+2}&\widehat{B}_{\frac{U}{2}+1}\widehat{G}_{\frac{U}{2}+1}&\ddots&\widehat{B}_{2}\widehat{G}_{2}\\\
\vdots&\vdots&\ddots&\vdots\\\
\widehat{B}_{U}\widehat{G}_{U}&\widehat{B}_{U-1}\widehat{G}_{U-1}&\cdots&\widehat{B}_{\frac{U}{2}}\widehat{G}_{\frac{U}{2}},\end{bmatrix}$
(80)

and L is the $(M+N_{1}+1)(N_{1}+1)$ matrix defined similarly using the Taylor
coefficients of log(1+z). The vectors $\mathbf{q}=\\{q_{m}\\}_{m=0}^{M}$ and
$\mathbf{l}=\\{l_{n_{1}}\\}_{{n_{1}}=0}^{N_{1}}$ hold the unknown polynomial
coefficients in order of increasing degree. As the column dimension of the
matrix in (79) is one greater than its row dimension, we can conclude that
there is one nonzero solution to (79). In many cases, this can be made into a
square system by choosing, say, $q_{0}=1$. However, if one does not want to
assume that any particular coefficient is nonzero, one can solve (79) by a
singular value decomposition (cf. Gonnet et al., 2013). Finally, the unknown
coefficients of $\mathbf{p}=\\{p_{n}\\}_{n=1}^{N}$ can be obtained by
multiplication through the following matrix system:

$\displaystyle\mathbf{p}=\begin{bmatrix}\widehat{B}_{0}\widehat{G}_{0}&&&\\\
\widehat{B}_{1}\widehat{G}_{1}&\widehat{B}_{0}\widehat{G}_{0}&&\\\
\vdots&\ddots&\ddots&\\\ \widehat{B}_{{U\over 2}}\widehat{G}_{{U\over
2}}&\cdots&\cdots&\widehat{B}_{0}\widehat{G}_{0}\end{bmatrix}\mathbf{q}-\begin{bmatrix}l_{0}&&&\\\
l_{1}&l_{0}&&\\\ \vdots&\ddots&\ddots&\\\ l_{{U\over 2}}&\cdots&\cdots&\
l_{0}\end{bmatrix}\mathbf{l}.$ (81)

If there is more than one jump location in the option pricing/Greeks curve
(76), this suggests the following modification of the equation:

$\displaystyle
P_{N}(z)+L_{N_{1}}(z)\log\left(1-{z\over\varepsilon_{1}}\right)+\ldots+L_{N_{s}}(z)\log\left(1-{z\over\varepsilon_{S}}\right)=f_{1}(z)Q_{M}(z)+\mathcal{O}(z^{U+1}).$
(82)

Accordingly, we have to modify (79) to produce a new $L$ matrix and a vector
of coefficients for each location to reflect the changes. According to
Driscoll and Fornberg (2001, 2011), there is no rigorous optimal formula for
choosing the degrees $M,$ $N,$ and $N_{1},$…,$N_{s}.$ Because the denominator
polynomial $Q_{M}$ is shared, we allow $M$ to be the largest, with the others
being equal as far as possible. For the case of just one jump location, taking
$N$ at roughly $40\%$ of the total available degrees of freedom seems to work
well. Experiments suggest that these choices can affect the observed accuracy,
occasionally by as much as an order of magnitude, but on average, there is
little variation within a broad range of choices.

The methods described above all assume that the locations of all jumps are
known in the option pricing/Greeks curve. For jumps whose locations are
unknown, we follow the suggestion of Driscoll and Fornberg (2001, 2011) to use
the Fourier-Padé algorithm to estimate their locations. As we approximate a
PDF rather than a payoff function with the CFS method (see the proof of
Theorem 1), we only consider jumps existing in the PDF. The existence of a
jump is attributed to a combination parameters and/or a very short maturity
that causes a sharp-peaked PDF curve.

Using the Fourier-Padé algorithm (cf. the Fourier-Padé approximant in Section
2) to estimate the jumps in a PDF is fairly simple. We first express the PDF
as the CFS representation:

$\displaystyle\mathfrak{Re}\left[2\sum_{k=1}^{\infty}\varphi\left(-\frac{2\pi}{d-c}k\right)e^{i\frac{2\pi}{b-a}ky}+\varphi\left(0\right)\right].$
(83)

Then, we can differentiate (83) with respect to $y$ to obtain

$\displaystyle\mathfrak{Re}\left[2\sum_{k=1}^{\infty}\left(i\frac{2\pi}{b-a}k\right)\varphi\left(-\frac{2\pi}{d-c}k\right)e^{i\frac{2\pi}{b-a}ky}\right].$
(84)

Finally, we let $z=\exp\left({i\frac{2\pi}{d-c}y}\right)$ in the two equations
above, and they are ready for the Fourier-Padé approximation. In general, when
the PDF has a jump, the sharp-peaked jump point will have an enormously large
value after differentiation. Then, based on this result, we find that point in
the original PDF. In other words, Figure 1 is a graphical illustration of the
outlooks of two PDFs (top) under the VG model (cf. Section 8.1.2) and the
first derivative (bottom) after the Fourier-Padé approximation. In the figure,
we can see that the PDFs in the top and bottom panels are smooth without any
jumps, as they do not show any point with an extremely large value. However,
in the graph on the bottom right, we can see that there is a jump that
produces a value of $15\times 10^{4}.$

Figure 1: Density functions (top) of the VG model and their first derivative
(bottom). A piecewise continuous (non-smooth) PDF with a jump after/before
differentiation is shown in the graphs on the top and bottom right. The
parameters are taken from VG–Para1.

##  6 Choice of Truncated Intervals

As we will show in Section 7, the choice of the interval $[c,d]$ is affected
by the accuracy of the SFP method. A minimum and substantial interval $[c,d]$
can capture most of the mass of a PDF such that the SFP method can in turns
produce sensible global spectral convergence rate. In this short section, we
show how to construct an interval related to the closed-form formulas of
stochastic process cumulants. The idea of using the cumulants is first
proposed by Fang and Oosterlee (2009a) to construct the definite interval
$[c,d]$ in (36). Based on their ideas, we have the following expression for
$[c,d]$:

$\displaystyle d$
$\displaystyle=\left|c_{1}+L\sqrt{c_{2}+\sqrt{c_{4}}}+\left|\log\left({S_{0}\over
K}\right)\right|\right|$ (85) $\displaystyle c$ $\displaystyle=-d,$ (86)

where $c_{1},$ $c_{2},$ and $c_{4}$ are the first, second and fourth
cumulants, respectively, of the stochastic process and $L\in[10,12].$ For
simple, less-complicated financial models, we also obtain closed-form formulas
for $c_{1},$ $c_{2},$ and $c_{4}$, which are shown in Table 13 of Appendix 10.
However, in the Heston model, we use the absolute value of $c_{2}$ and ignore
the value of $c_{4}$ due to the negative value of $c_{2}$ and the lengthy
representation of $c_{4}$ (cf. Fang and Oosterlee, 2009a). We therefore have
$c_{1}+L\sqrt{|c_{2}|}$ rather than $c_{1}+L\sqrt{c_{2}+\sqrt{c_{4}}}.$ The
truncated intervals only work for smooth PDFs without any jumps. In general,
by trial and error, if there is a jump in a PDF, we add 0.5 to (85) to allow
the better convergence of our method. Accordingly, the formulae can be
transformed as:

$\displaystyle d$
$\displaystyle=\left|c_{1}+L\sqrt{c_{2}+\sqrt{c_{4}}}+\left|\log\left({S_{0}\over
K}\right)\right|\right|+0.5$ (87) $\displaystyle c$ $\displaystyle=-d.$ (88)

######  Remark 2.

Without a doubt, in most cases, the truncated interval of (87) works for any
normal non-smooth PDF. However, when the conditions become very extreme, for
example, $T=1e-06$ in BSM–Para4 in Section 8.1.1, the PDF is extremely narrow,
thin and spiky at one point. The interval of (87) can still perform but is not
optimal. To improve it, we replace $0.5$ in the equation with $0.1.$ This is
attributed to the desire for higher accuracy and convergence with fewer terms
required over a short interval range.

##  7 Error Analysis

In this section, we demonstrate that the total error from pricing European-
style options can be made very small by choosing a suitably large interval
$[c,d].$ Furthermore, we show that global spectral convergence can be achieved
at any point along the option pricing curve even though the input PDF is a
$\mathbf{C}^{\nu}$ piecewise continuous (non-smooth) function.

In this paper, there are three types of approximation errors in any call/put
option.

  1. 1.

Integration truncation error:

$\displaystyle\epsilon_{1}:=\left|\int_{-\infty}^{\infty}G(e^{x+\chi},K)f(\chi)\mathrm{d}\chi-\int_{c}^{d}G(e^{x+\chi},K)f(\chi)\mathrm{d}\chi\right|$
(89)

  2. 2.

Error related to approximating (4.1) with (55):

$\displaystyle\epsilon_{2}$
$\displaystyle:=\Bigg{|}\mathfrak{Re}\left[2\sum_{k=1}^{\infty}\frac{1}{d-c}\int_{c}^{d}f(y)e^{-i\frac{2\pi}{d-c}ky}\mathrm{d}y\,\widehat{G}_{k}\,e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}+\frac{1}{d-c}\int_{c}^{d}f(y)\mathrm{d}y\,\widehat{G}_{0}\right]-$
$\displaystyle\quad\quad\mathfrak{Re}\left[2\sum_{k=1}^{\infty}\widehat{B}_{k}\,\widehat{G}_{k}\,e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}+\widehat{B}_{0}\,\widehat{G}_{0}\right]\Bigg{|}$
$\displaystyle=\Bigg{|}\mathfrak{Re}\Bigg{[}2\sum_{k=1}^{\infty}\left(\frac{1}{d-c}\int_{c}^{d}f(y)e^{-i\frac{2\pi}{d-c}ky}\mathrm{d}y-\widehat{B}_{k}\right)\,\widehat{G}_{k}\,e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}+$
$\displaystyle\quad\quad\left(\frac{1}{d-c}\int_{c}^{d}f(y)\mathrm{d}y-\widehat{B}_{0}\right)\,\widehat{G}_{0}\Bigg{]}\Bigg{|}$
(90)

  3. 3.

Truncated SFP series error:

$\displaystyle\epsilon_{3}$
$\displaystyle:=\left|\mathfrak{Re}\left[2\sum_{k=1}^{\infty}\widehat{B}_{k}\widehat{G}_{k}z^{k}+\widehat{B}_{0}\widehat{G}_{0}\right]-\mathfrak{Re}\left[{P_{N}^{+}(z)+\sum_{s=1}^{S}L^{+}_{N_{s}}(z)\log\left(1-z/\varepsilon_{s}\right)\over
Q^{+}_{M}(z)}\right]\right|,$ $\displaystyle z=e^{i{2\pi\over d-c}(-x+\log
K)}\hbox{ and }\varepsilon_{s}=e^{i{2\pi\over d-c}(\zeta_{s})}$ (91)

If we introduce the concept of the cumulative probability density function
(CDF) $F(\chi)$ such that $f(\chi)d\chi=dF(\chi),$ we can simplify the
integration truncation error as follows:

$\displaystyle\epsilon_{1}$
$\displaystyle=\left|\left(\int_{-\infty}^{\infty}G(e^{x+\chi},K)f(\chi)\mathrm{d}\chi-\int_{c}^{d}G(e^{x+\chi},K)f(\chi)\mathrm{d}\chi\right)\right|$
$\displaystyle=\left|\left(\int_{-\infty}^{c}G(e^{x+\chi},K)f(\chi)\mathrm{d}\chi+\int_{d}^{\infty}G(e^{x+\chi})f(\chi)\mathrm{d}\chi\right)\right|$
$\displaystyle\leq\left|\left(\int^{c}_{-\infty}{\partial
G(e^{x+\chi},K)\over\partial\chi}F(\chi)\mathrm{d}\chi\right)\right|+\left|\int_{d}^{\infty}{\partial
G(e^{x+\chi},K)\over\partial\chi}(1-F(\chi))\mathrm{d}\chi\right|$ (92)
$\displaystyle\approx 0:\quad(\hbox{if}\,\chi=c,d,-\infty,\infty).$ (93)

We can see that $\epsilon_{1}$ is bounded and approaches zero as long as
$[c,d]$ is chosen reasonably such that $1-F(d)\approx 0$ when $d<\infty$ or
$F(c)\approx 0$ when $c>-\infty.$ We are also able to adapt the same idea to
investigate the bound of $\epsilon_{2}.$ Accordingly, taking into account
$|\exp(i{2\pi k\over d-c}y)|\leq 1,$ we first investigate the error

$\underline{\epsilon}_{2}:=\left|{1\over
d-c}\int_{c}^{d}f(y)e^{-i\frac{2\pi}{d-c}ky}\mathrm{d}y-\widehat{B}_{k}\right|$

in $\epsilon_{2}.$ If we expand the equation above, we obtain

$\displaystyle\underline{\epsilon}_{2}$ $\displaystyle:=\left|{1\over
d-c}\int_{c}^{d}f(y)e^{-i\frac{2\pi}{d-c}ky}\mathrm{d}y-{1\over
d-c}\varphi\left(-i\frac{2\pi}{d-c}k\right)\right|$ (94)
$\displaystyle=\left|{1\over
d-c}\int_{-\infty}^{\infty}f(y)e^{-i\frac{2\pi}{d-c}ky}\mathrm{d}y-{1\over
d-c}\int_{c}^{d}f(y)e^{-i\frac{2\pi}{d-c}ky}\mathrm{d}y\right|$ (95)
$\displaystyle\leq\left|{1\over
d-c}\left(\int_{-\infty}^{c}f(y)\mathrm{d}y+\int_{d}^{\infty}f(y)\mathrm{d}y\right)\right|$
(96) $\displaystyle=\left|{1\over
d-c}\left(F(\infty)-F(d)+F(c)-F(-\infty)\right)\right|$ (97)
$\displaystyle\approx 0:\quad(\hbox{if}\,y=c,d,-\infty,\infty).$ (98)

Based on the result above,

$\displaystyle\epsilon_{2}$
$\displaystyle:=\Bigg{|}\mathfrak{Re}\Bigg{[}2\sum_{k=1}^{\infty}\left(\frac{1}{d-c}\int_{c}^{d}f(y)e^{-i\frac{2\pi}{d-c}ky}\mathrm{d}y-\widehat{B}_{k}\right)\,\widehat{G}_{k}\,e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}+$
$\displaystyle\quad\quad\left(\frac{1}{d-c}\int_{c}^{d}f(y)\mathrm{d}y-\widehat{B}_{0}\right)\,\widehat{G}_{0}\Bigg{]}\Bigg{|}$
$\displaystyle\leq\Bigg{|}\mathfrak{Re}\Bigg{[}2\sum_{k=1}^{\infty}\underline{\epsilon}_{2}\,\widehat{G}_{k}\,e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}+\left(\frac{1}{d-c}\int_{c}^{d}f(y)\mathrm{d}y-\frac{1}{d-c}\int_{c}^{d}f(y)\mathrm{d}y\right)\,\widehat{G}_{0}\Bigg{]}\Bigg{|}$
$\displaystyle\approx 0.$ (99)

To conclude that $\epsilon_{2}$ is approaching zero, we first note that there
is no approximate error of $G_{k}$ and $G_{0}$ because of their closed-form
expressions. Then, once $\underline{\epsilon}_{2}$ tends to zero, the first
term of the equation will also diminish to zero. The last term of the equation
tends to zero because $\widehat{B}_{0}$ equals
$\frac{1}{d-c}\int_{c}^{d}f(y)\mathrm{d}y$ (cf. (51) and (54)).

Finally, the SPF series truncation error is also bounded (cf. Driscoll and
Fornberg, 2001, 2011) and can be formulated as follows:

$\displaystyle\epsilon_{3}$
$\displaystyle:=\left|\mathfrak{Re}\left[2\sum_{k=1}^{\infty}\widehat{B}_{k}\widehat{G}_{k}e^{i\frac{2\pi}{d-c}k(-x+\log
K)}+\widehat{B}_{0}\widehat{G}_{0}\right]-\mathfrak{Re}\Bigg{[}{P_{N}^{+}(z)+\sum_{s=1}^{S}L^{+}_{N_{s}}(z)\log\left(1-z/\varepsilon_{s}\right)\over
Q^{+}_{M}(z)}\Bigg{]}\right|$
$\displaystyle=\mathcal{O}(z^{N+M+\sum_{s=1}^{N_{s}}N_{s}+1}).$ (100)

The error term of $\epsilon_{3}$ tends to zero with a global spectral rate of
$\mathcal{O}(z^{N+M+\sum_{s=1}^{N_{s}}N_{s}+1})$ even though the input PDF is
$C^{\nu}$ piecewise continuous.

Before we illustrate the total error bound when approximating any true
European-type option price $V(x,K,t)$ defined as

$\displaystyle\epsilon$
$\displaystyle:=\left|V(x,K,t)-e^{-r(T-t)}\mathfrak{Re}\left[P_{N}^{+}(z)+\sum_{s=1}^{S}L^{+}_{N_{s}}(z)\log\left(1-z/\varepsilon_{s}\right)\over
Q^{+}_{M}(z)\right]\right|,$ (101)

we first summarise the whole approximation procedure of European-type option
prices and note where $\epsilon_{1},$ $\epsilon_{2},$ and $\epsilon_{3}$ lie.
We start off by seeking a definite interval $[c,d]$ that allows us to
approximate $V(x,K,t)$ defined on $[-\infty,\infty]$ in (33) with the form

$V(x,K,t)\approx e^{-r(T-t)}\int_{c}^{d}G(e^{x+\chi},K)f(\chi)\mathrm{d}\chi.$

The interval $[c,d]$ we proposed satisfies condition (36). As a result, we
obtain our first approximation error $\epsilon_{1}.$ As $V(x,K,t)$ is now
approximated in $[c,d],$ this implies that we can construct a CFS expansion of
$V(x,K,t),$ like the one in (4.1). Then, because including a characteristic
function $\varphi(\cdot)$ in the CFS expansion allows for a more accurate
approximation, we have another CFS expansion of $V(x,K,t),$ given as in (55).
Accordingly, we have $\epsilon_{2},$ an approximation error of (4.1), being
approximated by (55). Finally, $\epsilon_{3}$ is the error of (58), which is
the formula for approximating (55) with the SFP approximant (16).

By combining the results of $\epsilon_{1},$ $\epsilon_{2}$ and $\epsilon_{3},$
we can determine the total error bound $\epsilon$; hence, we have an
inequality of

$\displaystyle\epsilon$
$\displaystyle=\left|V(x,K,t)-e^{-r(T-t)}\mathfrak{Re}\left[P_{N}^{+}(z)+\sum_{s=1}^{S}L^{+}_{N_{s}}(z)\log\left(1-z/\varepsilon_{s}\right)\over
Q^{+}_{M}(z)\right]\right|$ (102)
$\displaystyle=\left|e^{-r(T-t)}\Bigg{(}\int_{-\infty}^{\infty}G(e^{x+\chi},K)f(\chi)\mathrm{d}\chi-\mathfrak{Re}\left[P_{N}^{+}(z)+\sum_{s=1}^{S}L^{+}_{N_{s}}(z)\log\left(1-z/\varepsilon_{s}\right)\over
Q^{+}_{M}(z)\right]\Bigg{)}\right|$
$\displaystyle\leq\bigg{|}e^{-r(T-t)}\bigg{|}\Bigg{(}\left|\int_{-\infty}^{\infty}G(e^{x+\chi},K)f(\chi)\mathrm{d}\chi-\int_{c}^{d}G(e^{x+\chi},K)f(\chi)\mathrm{d}\chi\right|+$
$\displaystyle\quad\Bigg{|}\mathfrak{Re}\left[2\sum_{k=1}^{\infty}\frac{1}{d-c}\int_{c}^{d}f(y)e^{-i\frac{2\pi}{d-c}ky}\mathrm{d}y\,\widehat{G}_{k}\,e^{i\frac{2\pi}{d-c}k\left(-x+\log
K\right)}+\frac{1}{d-c}\int_{c}^{d}f(y)\mathrm{d}y\,\widehat{G}_{0}\right]-$
$\displaystyle\quad\mathfrak{Re}\left[2\sum_{k=1}^{\infty}\widehat{B}_{k}\widehat{G}_{k}e^{i\frac{2\pi}{d-c}k(-x+\log
K)}+\widehat{B}_{0}\widehat{G}_{0}\right]\Bigg{|}+$
$\displaystyle\quad\left|\mathfrak{Re}\left[2\sum_{k=1}^{\infty}\widehat{B}_{k}\widehat{G}_{k}z^{k}+\widehat{B}_{0}\widehat{G}_{0}\right]-\mathfrak{Re}\left[P_{N}^{+}(z)+\sum_{s=1}^{S}L^{+}_{N_{s}}(z)\log\left(1-z/\varepsilon_{s}\right)\over
Q^{+}_{M}(z)\right]\right|\Bigg{)}$
$\displaystyle\leq|e^{-r(T-t)}|(\epsilon_{1}+\epsilon_{2}+\epsilon_{3})|$
$\displaystyle<|e^{-r(T-t)}|(\epsilon_{1}+\epsilon_{2}+\mathcal{O}(z^{N+M+\sum^{S}_{s=1}N_{s}+1}))|$
$\displaystyle\approx 0.$ (103)

######  Remark 3.

According to Driscoll and Fornberg (2001, 2011), the rate of
$\mathcal{O}(z^{N+M+\sum_{s=1}^{N_{s}}N_{s}+1})$ is not spatially uniform
because convergence at a jump is somewhat limited by the well-known numerical
ill-conditioning of the straightforward Padé problem. Moreover, when a curve
has a discontinuous point (zeroth-order jump), it also has a convergence rate
that is the average of the one-sided limit values. Through the numerical
experiments in Driscoll and Fornberg (2001), we can see that the SFP method
can yield spectral convergence rate at jumps. This is not necessarily the case
if the jumps are very difficult to interpolate. Nevertheless, as Driscoll and
Fornberg (2011) suggest, we can still have 4–6 digits of accuracy at the jumps
when they are known in advance. When the jumps are not known in advance, we
can use the Fourier-Padé algorithm to find them (see Section 5). By
incorporating this technique, the SFP method still can perform and yield
spectral convergence or 4th-order convergence when the jumps are very
difficult to interpolate.

##  8 Numerical Results

In this section, we demonstrate the performance of the SFP method through
various numerical tests. The purpose of this section is first to test whether
the error convergence analysis presented in Section 7 is in line with the
numerical findings in this section. Second, we test the ability of the SFP
method to price any European-style option that are deep in/out of the money
and have long/short maturities. Third, we analyse whether the SFP method can
provide consistent accuracy when approximating small or large values of option
prices. Finally, in a major development of the method, we test whether the SFP
method can retain global spectral convergence even when the PDF is piecewise
continuous. A number of popular numerical methods are implemented to test the
SFP method in terms of the error convergence, convergence rate and
computational time. These methods include the COS method (a Fourier COS series
method, Fang and Oosterlee, 2009a), the filter-COS method (a COS method with
an exponential filter to resolve the Gibbs phenomenon; see Ruijter et al.,
2013), the CONV method (an FFT method, Lord et al., 2008), the Lewis-FRFT
method (a fractional FFT method, Lewis, 2001; Chourdakis, 2004), and the
B-spline, Haar wavelet and the SWIFT methods (a wavelet-based method; see
Ortiz-Gracia and Oosterlee, 2013, 2016). When we implement the CONV and Lewis-
FRFT methods, we use Simpson’s rule for the Fourier integrals to achieve
fourth-order accuracy. In the filter-COS method, we use an exponential filter
and set the accuracy parameter to $10$ as Ruijter et al. (2013) report that
this filter provides better algebraic convergence than the other options. We
also set the damping factors of the CONV and Lewis-FRFT to 0 and any value
greater than zero, respectively.

As the SFP method requests approximating jumps in logarithmic series, we
consider and apply the endpoints $c$ and $d$ as our two known jumps for all
non-smooth/smooth PDFs. Only the jump of the non-smooth PDF of the Black-
Scholes-Merton (BSM) model is known as its mean value. For the rest of the
non-smooth PDFs, we use the Fourier-Padé algorithm (cf. Section 5) to find
their locations. In all numerical experiments, we use the parameter $U$ to
denote the number of terms of the SFP method and $N$ to denote the number of
terms/grid points of the others. When we measure the approximation errors of
the numerical methods, we use absolute errors, the infinity norm errors
$R_{\infty}$ and the $L_{2}$ norm errors $R_{2}$ as the measurement units.
Moreover, to improve the accuracy of our method, we use the call-put
parity–$V_{call}(x,K,t)=V_{put}(x,K,t)+S_{0}\exp(-qT)-K\exp(-rT)$–to
approximate call prices once we have put prices ready. Finally, all the CPU
times presented (in seconds) are determined after averaging the computational
time over 120 experiments. A MacBook Pro with a 2.8 GHz Intel Core i7 CPU and
two 8 GB DDR SDRAM (cache memory) is used for all experiments. The code is
written in MATLAB R2011b. Finally, the MATLAB code of implementing the COS
method and the FFT method, such as the CONV method and the like, is retrieved
from von Sydow et al. (2015).

###  8.1 Exponential Lévy Processes

####  8.1.1 The Black-Scholes-Merton Model

The first numerical experiments are performed using the BSM model (cf. Black
and Scholes, 1973; Merton, 1973). The stock dynamics driven by the BSM model
(a geometric Brownian process) are given by

$\displaystyle S_{T}=S_{0}e^{(r-q-{1\over 2}\sigma^{2})T+\sigma W_{T}},$ (104)

where $W_{T}$ is a risk-neutral Brownian motion, and $\sigma$ is the
volatility. The characteristic function of the model is also defined as

$\displaystyle\varphi(u)=\exp\left(T\left(iu(r-q-{1\over
2}\sigma^{2})-\frac{1}{2}\sigma^{2}u^{2}\right)\right),\,\,u\in\mathbb{R}.$
(105)

The parameters for the experiments are selected from the following:

$\displaystyle\textbf{BSM--Para1}:\,S_{0}$
$\displaystyle=100,\,\sigma=0.15,\,r=0.03,\,T=1.0,\,q=0.0,$ (106)
$\displaystyle\textbf{BSM--Para2}:\,S_{0}$
$\displaystyle=100,\,\sigma=0.25,\,r=0.1,\,T=50\mbox{ or }100,\,q=0,\,K=120,$
(107) $\displaystyle\textbf{BSM--Para3}:\,K$
$\displaystyle=100,\,\sigma=0.2,\,r=0.06,\,T=1e-06,\,q=0.$ (108)

In the first numerical test (BSM–Para1), we first check for convergence
behaviour against a range of strikes $K$ from $1$ to $200$ for deep in/out-of-
the money and at-the-money vanilla put options, respectively. The parameters
are retrieved from von Sydow et al. (2015). Using the same parameters, we
check for convergence behaviour against a range of strikes $K$ from $80$ to
$120$ for cash-or-nothing put options rather than vanilla options in our
second numerical test. We compare our method with the COS method, the CONV
method and the Lewis-FRFT method in the first test but only with the COS
method in the second test. The third numerical experiment (BSM–Para2) is
devoted to comparing the performance of the COS method, the SFP method, the
SWIFT method for long maturity call options. As we sometimes encounter these
options in the insurance and pension industry, it is worth testing our method
against them. The parameters are retrieved from Ortiz-Gracia and Oosterlee
(2016) for the test. Finally, for the last numerical test (BSM–Para3), we use
the values of $\sigma,$ $r$ and $q$ from Andricopoulos et al. (2003) and check
for convergence behaviour against a range of stock prices $S_{0}$ from $80$ to
$120$ for very short maturity vanilla call options and their option Delta and
Gamma. Our method is compared with the COS, filter-COS, Lewis-FRFT and CONV
methods. Reference values for all the tests are based on the BSM analytical
formula. In each numerical test, except for the third one, we declare 250
different option prices within the range of either $K$ or $S_{0}$ to test the
efficiency of our method and the others. We represent all recovered PDFs via
the SFP method in Figure 2 for all sets of parameters. We can see that the
non-smooth recovered PDF (top right) can be obtained via BSM–Para3. We set
$L=10$ in the truncated interval proposed in Section 6. One should note that
for the truncated interval in the last numerical test BSM–Para3, we replace
$0.5$ with $0.1$ in (87) (see Remark 2), as we want our method to have higher
accuracy and a better convergence rate. It is quite rare for the maturity time
to be infinitesimal as in BSM–Para3, but to demonstrate the efficiency of our
method, it is worth performing such a test.

All the tables in this section and the others suggest that the difference in
the computational time across methods is not large. It takes less than 0.1
seconds to approximate $250$ option prices for any method when $N$ and $U$
equal 64. Without considering the accuracy of the methods, this is a quite
reasonable time frame in which to produce a substantial number of option
prices at once. However, when we consider the error convergence and
convergence rate of the methods, there are sizeable differences among them.
Figure 3 and its numerical presentation in Table 3 are dedicated to the first
numerical test (BSM–Para1). In this test, the SFP method has the smallest
$R_{\infty}$ and $R_{2}$ errors, with $R_{\infty}$ and $R_{2}$ less than
$5.8e-13$ as $U=64$ when the strike increases from 1 to 200. It also yields
global spectral convergence against which none of the other methods can
compete. This suggests that the truncated range proposed in this paper can
work over longer range than that proposed in Fang and Oosterlee (2009a). In
Table 4, we use the same parameters (BSM–Para1) to examine the ability of the
COS and SFP methods to price cash-or-nothing put options. Again, the two
methods can achieve very low convergence error and spectral convergence. Table
5 refers to the third test (BSM–Para2) and replicates Table 3 in Ortiz-Gracia
and Oosterlee (2016). In this test, the SFP method impressively provides high
accuracy over the SWIFT and COS methods with fewer terms required. For
example, when $T=50,$ we only need 32 terms to obtain the true solution with
7-digit accuracy in the SFP method, whereas the other two require more terms
to obtain 1-digit accuracy. Finally, in the last test (BSM–Para4), as the
cumulants are too small to create a meaningful truncated range for the COS and
filter-COS methods, we use (87) to allow the method to generate relevant
option prices. In Figure 4, the SFP method dominates the other four methods to
yield the global spectral convergence rate away from the jump, while the input
non-smooth PDF (see top right, Figure 2) behaves like a Dirac Delta
function444The Dirac Delta function is a generalised function or distribution
on the real number line that is zero everywhere except at zero. and has an
almost infinitely thin spike near the origin. In Table 6, when we test the
accuracy of the option prices around/at the jump—$S_{0}=99.999$—the SFP method
dominate the other methods. It has zero absolute error in the case of
$S_{0}=95$ and $6.269e-05$ error from the true solution in the case of
$S_{0}=99.999$ when $U=64.$ This is in line with the findings of Driscoll and
Fornberg (2011) as the SFP offers 4-6 digit accuracy at the jumps when the
function is very difficult. Obviously, the non-smooth PDF in the last test is
very difficult to approximate, as it behaves like a Dirac Delta. Finally,
using the same parameters of BSM–Para3, we recover the call Delta and Gamma
via the SFP approximant in Figure 5. In the graphs, the SFP provides a
solution of extremely high accuracy and global spectral convergence apart from
the jump.

Figure 2: Probability density functions are generated by three sets of parameters: BSM–Para1 (top left), BSM–Para3 (top right) and BSM–Para2 (bottom). Figure 3: Error envelopes of the CONV, Lewis-FRFT, COS and SFP methods for pricing European deep at/in/out-of-the-money put options under the BSM model with strikes increasing from 1 to 200. The number of terms/grid points increases in a sequence of 8 (red line), 16 (blue line), 32 (green line) and 64 (black line). Spectral convergence is observed in the SFP method. The parameters are taken from BSM–Para1. Table 3: Comparison of the CONV, Lewis-FRFT, COS and SFP methods in error convergence and CPU time for pricing European deep at/in/out-of-the-money put options under the BSM model. The range of $K$ is from 1 to 200. Spectral convergence is observed in the SFP method. The parameters are taken from BSM–Para1. CONV | Lewis-FRFT
---|---
N | $R_{\infty}$ | $R_{2}$ | Time | N | $R_{\infty}$ | $R_{2}$ | Time
8 | 20.49 | 180.1 | 0.00806 | 8 | 32.53 | 356 | 0.02111
16 | 6.292 | 53.54 | 0.00941 | 16 | 16.1 | 156.4 | 0.02525
32 | 1.998 | 13.83 | 0.01012 | 32 | 7.933 | 69.38 | 0.04571
64 | 0.5709 | 3.286 | 0.01641 | 64 | 4.765 | 38.81 | 0.09151
COS | SFP
N | $R_{\infty}$ | $R_{2}$ | Time | U | $R_{\infty}$ | $R_{2}$ | Time
8 | 5.046 | 10.06 | 0.00706 | 8 | 1.028e-01 | 5.893e-01 | 0.00901
16 | 5.109 | 8.447 | 0.00841 | 16 | 2.818e-04 | 1.559e-03 | 0.00951
32 | 5.111 | 8.449 | 0.01113 | 32 | 1.598e-09 | 3.732e-09 | 0.01312
64 | 5.111 | 8.449 | 0.01423 | 64 | 1.991e-13 | 5.801e-13 | 0.01731
Table 4: Comparison of the COS and SFP methods in error convergence and CPU time for pricing European cash-or-nothing put options under the BSM model. The range of $K$ is from 80 to 120. Spectral convergence is observed in the COS and SFP methods. The parameters are taken from BSM–Para1. COS | SFP
---|---
N | $R_{\infty}$ | $R_{2}$ | Time | U | $R_{\infty}$ | $R_{2}$ | Time
8 | 6.395e-02 | 7.711e-01 | 0.00716 | 8 | 3.673e-03 | 3.974e-02 | 0.00901
16 | 2.562e-03 | 2.831e-02 | 0.00831 | 16 | 3.801e-06 | 2.646e-05 | 0.00931
32 | 7.681e-08 | 8.794e-07 | 0.01123 | 32 | 5.702e-12 | 2.668e-11 | 0.01342
64 | 1.772e-15 | 4.015e-15 | 0.01523 | 64 | 1.156e-14 | 2.297e-14 | 0.01791
Table 5: Comparison of the COS, SWIFT and SFP methods in terms of absolute error for pricing a call option under the BSM model. The reference values have been computed using the Black-Scholes analytical formulae: 99.2025928525532000 ($T=50$) and 99.9945609694213000 ($T=100$). The SFP method is more accurate than the others with fewer summation terms required. The parameters are taken from BSM–Para3. Method | Error $(T=50)$ | Error $(T=100)$
---|---|---
SWIFT $(m=0)$ | $1.91e-01$ | $2.50e-05$
COS $(N=35)$ | $4.98e-01$ | $2.05e+02$
SFP $(U=32)$ | $2.653e-07$ | $7.067e-08$
SWIFT $(m=1)$ | $7.78e-09$ | $3.20e-06$
COS $(N=70)$ | $2.79e-08$ | $2.02e-05$
SFP $(U=64)$ | $2.251e-10$ | $7.037e-11$
Figure 4: Error envelopes of the CONV, Lewis-FRFT, filter-COS and SFP methods
when pricing European at/around-the-money call options under the BSM model in
an $S_{0}$ range from 80 to 120. The number of terms/grid points increases in
a sequence of 8 (red line), 16 (blue line) and 32 (green line). In the SFP
method, spectral convergence is observed away from the jump. The parameters
are taken from BSM–Para4. Table 6: Comparison of the CONV, COS, filter-COS and
SFP methods in terms of the absolute error for pricing a call option under the
BSM model. The reference values are $0$ ($S_{0}=95$) and $0.007491657716010$
($S_{0}=99.999$). The SFP method is more accurate than the others with fewer
summation terms required. The parameters are taken from BSM–Para4.

$S_{0}=95$ | CONV | filter-COS | SFP | $S_{0}=99.999$ | COS | filter-COS | SFP
---|---|---|---|---|---|---|---
$N/U$ | Error | Error | Error | $N/U$ | Error | Error | Error
8 | 2.004e+01 | 5.099 | 2.813e-04 | 8 | 2.475e-01 | 4.482e-01 | 3.849e-02
16 | 1.940e+01 | 4.987 | 7.268e-08 | 16 | 1.307e-01 | 1.944e-01 | 4.704e-03
32 | 1.311e+01 | 5.001 | 0.000 | 32 | 7.008e-02 | 1.030e-01 | 3.473e-03
64 | 6.696 | 5.000 | 0.000 | 64 | 3.960e-2 | 5.587e-02 | 6.268e-05

Figure 5: Recovered call Delta $\Delta$ (top left) and Gamma $\Gamma$ (top
right) via the SFP method under the BSM model with strikes increasing from 80
to 120 and their corresponding error envelopes: Delta $\Delta$ (bottom left)
and Gamma $\Gamma$ (bottom right). The number of terms increases in a sequence
of 8 (red line), 16 (blue line) and 32 (green line). In the SFP method,
spectral convergence is observed away from the jump. The parameters are taken
from BSM–Para4.

####  8.1.2 The Variance Gamma Model

A VG process (Madan and Seneta, 1990; Madan and Milne, 1991; Madan et al.,
1998) is an infinite activity Lévy process and is a subordinate version of
Brownian motion (cf. Cont and Tankov, 2004). The most important feature of
this model is that the Brownian motion is evaluated in random time $t^{*}$
(determined by an independent increasing Lévy process—a Gamma process) rather
than in calendar time $t$. Suppose that the VG process
$b(t^{*};\theta,\sigma)$ is defined as $\theta t^{*}+\sigma W_{t^{*}},$ where
the random time $t^{*}$ is given by a Gamma process $\rm Gamma(t;1,\upsilon)$
with a unit mean and variance $\upsilon$, $\theta$ is a drift at $t^{*}$, and
$W_{t^{*}}$ denotes a standard Brownian motion. Then, we define the stock
price dynamics driven by the VG process as follows:

$\displaystyle S_{T}=S_{0}e^{(r-q-\omega)t+\theta{\rm
Gamma(t;1,\upsilon)}+\sigma W_{{\rm Gamma(t;1,\upsilon)}}},$ (109)

and its characteristic function is given by

$\displaystyle\varphi(u)=\exp\bigg{(}iu(r-q-\omega)t\bigg{)}\left(\frac{1}{1-i\theta\upsilon
u+\frac{\sigma^{2}\upsilon}{2}u^{2}}\right)^{\frac{t}{\upsilon}},\,\,u\in\mathbb{R}.$
(110)

Here
$\omega=-{1\over\nu}\log\left(1-\theta\upsilon-\frac{\sigma^{2}\upsilon}{2}\right).$
The parameters for the experiments are selected from the followings:

$\displaystyle\textbf{VG--Para1}:\,S_{0}$
$\displaystyle=100,\,\sigma=0.12,\,\theta=-0.14,\,\nu=0.2,\,r=0.1,\,T=0.1,\,q=0,$
(111) $\displaystyle\textbf{VG--Para2}:\,K$
$\displaystyle=1,\,\sigma=0.1213,\,\theta=-0.1436,\,\nu=0.1686,\,r=0.03,\,T=1,\,q=0.01.$
(112)

The first set of parameters is chosen because relatively slow convergence was
reported for the CONV method for very short maturities in Lord et al. (2008).
The last set is from Pistorius and Stolte (2012), which originates from Madan
et al. (1998). The reference values for these tests are based on the VG
analytical formula. We present all the recovered PDFs via the SFP method in
Figure 6, and the non-smooth recovered PDF (red line) can be obtained via
VG–Para1. We set $L=10$ in the truncated intervals of (85) and (87). In first
numerical test (VG–Para1), we check for convergence behaviour against 250
different call prices within a range of $K$ from $80$ to $120$ and then
examine some individual call prices around/at the jump. Global spectral
convergence away from the jump is reported for the SFP method in Figure 7. In
Table 7, when $K=90,$ the SFP method can achieve 13 digits of accuracy with
128 terms required. This result is far better than those of the other two
methods, as they only achieve 4 or 5 digits of accuracy. In the same test,
Fang and Oosterlee (2009a) show that the COS method requires 1024 terms to
achieve an error difference of $2.52e-08$ from the true solution. Obviously,
their result indicates that the COS method cannot compete with the SFP method.
When $S_{0}=102.336,$ the call price is measured at the jump. We can see that
the convergence rate of the call price becomes algebraic. Compared with other
methods, the SFP method is still more accurate with 6 digits of accuracy and
only 128 terms required. This is in line with the finding in Driscoll and
Fornberg (2011) that the SFP method can yield 4-6 digits of accuracy at the
jump if the function is very difficult. In the last test (VG–Para2), we check
for convergence behaviour against 250 call prices within a range of $S_{0}$
from $0.5$ to $2.$ Figure 8 and its numerical representation—Table 8—suggest
that the COS and SFP method have global spectral convergence when the PDF is
smooth. However, the SFP method has the highest accuracy in terms of
$R_{\infty}$ and $R_{2}.$ Moreover, in the test, we notice that the COS method
yields less stable and accurate numerical results when it is used to
approximate small values of option prices.

Figure 6: Recovered density functions for the VG model are generated via two
sets of parameters: VG–Para1 (red line) and VG–Para2 (blue line). Figure 7:
Error envelopes of the CONV, filter-COS, COS and SFP methods for pricing
European at/around-the-money call options under the VG model. The strikes
increase from 80 to 120. The number of terms/grid points increases in a
sequence of 8 (red line), 16 (blue line), 32 (green line) and 64 (black line).
Apart from the jump, spectral convergence is observed in the SFP method. The
parameters are taken from VG–Para1. Table 7: Comparison of the COS, filter-COS
and SFP methods in terms of absolute error for pricing a vanilla call option
under the VG model. The reference values are $10.993703186728190$ ($K=90$) and
$0.689027011772653$ ($K=102.336$). The SFP method is more accurate than the
others with fewer summation terms required. The parameters are taken from
VG–Para1.

$K=90$ | COS | filter-COS | SFP | $K=102.336$ | COS | filter-COS | SFP
---|---|---|---|---|---|---|---
$N/U$ | Error | Error | Error | $N/U$ | Error | Error | Error
8 | 5.783e-01 | 3.062e-02 | 5.256e-02 | 8 | 1.326 | 2.413 | 5.506
16 | 5.843e-02 | 3.152e-01 | 1.241e-04 | 16 | 4.185e-01 | 7.921e-01 | 5.435
32 | 1.405e-02 | 4.756e-02 | 2.221e-08 | 32 | 1.153e-01 | 2.488e-01 | 2.895e-04
64 | 1.603e-03 | 7.951e-04 | 1.401e-11 | 64 | 2.976e-02 | 6.836e-02 | 2.899e-05
128 | 4.281e-04 | 8.931e-05 | 5.755e-13 | 128 | 7.595e-03 | 1.762e-02 | 1.147e-06

Figure 8: Error envelopes of the CONV, Lewis-FRFT, COS and SFP methods for pricing European at/around-the-money call options under the VG model in the range of $S_{0}$ from 0.5 to 2. The number of terms/grid points increases in a sequence of 8 (red line), 16 (blue line), 32 (green line) and 64 (black line). Spectral convergence is observed for the COS and SFP methods. The SFP method has the highest accuracy. The parameters are taken from VG–Para2. Table 8: Comparison of the CONV, Lewis-FRFT, COS and SFP methods in terms of error convergence and CPU time for pricing European at/around-the-money call options under the VG model in the range of $S_{0}$ from 0.5 to 2. Spectral convergence is observed in both the COS and SFP methods. The SFP method has the highest accuracy. The parameters are taken from VG–Para2. CONV | Lewis-FRFT
---|---
N | $R_{\infty}$ | $R_{2}$ | Time | N | $R_{\infty}$ | $R_{2}$ | Time
8 | 1.384 | 13.31 | 0.00816 | 8 | 3.147e-01 | 3.913 | 0.02111
16 | 5.933e-01 | 4.344 | 0.00941 | 16 | 1.506e-01 | 1.645 | 0.02625
32 | 1.161e-01 | 6.922e-01 | 0.01012 | 32 | 7.126e-02 | 6.727e-01 | 0.04571
64 | 7.092e-03 | 4.533e-02 | 0.01641 | 64 | 3.986e-02 | 3.245e-01 | 0.09151
COS | SFP
N | $R_{\infty}$ | $R_{2}$ | Time | N | $R_{\infty}$ | $R_{2}$ | Time
8 | 1.477e-02 | 1.292e-0 | 0.00716 | 8 | 2.174e-03 | 1.476e-02 | 0.00921
16 | 1.472e-03 | 1.131e-02 | 0.00841 | 16 | 5.234e-05 | 2.300e-04 | 0.00961
32 | 2.217e-05 | 1.596e-04 | 0.01223 | 32 | 2.409e-09 | 7.566e-09 | 0.01372
64 | 4.694e-08 | 2.857e-07 | 0.01413 | 64 | 1.541e-11 | 2.485e-11 | 0.01741

####  8.1.3 The CGMY Model

The CGMY model was developed by Carr et al. (2002) and can be seen as a
generalisation of the VG model discussed above. Carr et al. (2002) introduced
the CGMY model as a class of infinitely divisible distributions (also known as
tempered stable processes; see Cont and Tankov, 2004). The Lévy measure for
the CGMY process is given by

$\displaystyle\nu_{\rm(CGMY)}(d\chi)=\begin{cases}&C\exp(-G|\chi|)/|\chi|^{Y+1}\mathrm{d}\chi,\,\,\chi<0,\\\
&C\exp(-M|\chi|)/|\chi|^{Y+1}\mathrm{d}\chi,\,\,\chi>0,\end{cases}$ (113)

where $C>0$, $G>0$, $M>0$, and $Y<2$. The parameter $Y$ captures the fine
structure of the process. For $Y<-1$, we obtain a compound Poisson process
that has finite variation and finite activity. However, when $Y\in[0,1)$, the
process has infinite activity and finite variation, which is similar to a VG
process (we can see that when $Y=0$, this process is equivalent to a VG
process). For $Y\in[1,2)$, the process has infinite activity and infinite
variation. In this paper, we focus on a CGMY process with $Y\in(0,2)/\\{1\\}$,
so its characteristic function is defined as follows:

$\displaystyle\varphi(u)$
$\displaystyle=\exp\Bigg{(}iu(r-q+\omega)+C\Gamma(-Y)G^{Y}\left(\left(1+\frac{iu}{G}\right)^{Y}-1-\frac{iuY}{G}\right)$
$\displaystyle\,+C\Gamma(-Y)M^{Y}\left(\left(1-\frac{iu}{M}\right)^{Y}-1+\frac{iuY}{M}\right)\Bigg{)}.$
(114)

Here, $\omega=\varphi(-i).$ The parameters (CGMY–Para1) for the numerical test
are drawn from Fang and Oosterlee (2009a) and defined as follows:

$\displaystyle\textbf{CGMY--Para1}:\,$ $\displaystyle
S_{0}=100,\,C=1,\,G=5,\,M=5,\,Y=0.5\mbox{ or }1.5\mbox{ or }1.98,$ (115)
$\displaystyle r=0.1,\,T=1,\,q=0.$

We evaluate the SFP method’s convergence rate for vanilla calls and puts under
the CGMY model. Almendral and Oosterlee (2007) and Wang et al. (2007) have
reported that using the finite difference method to solve partial differential
integral equations to obtain option prices is difficult in cases when
$Y\in[1,2).$ Therefore, we evaluate the SFP method using $Y=0.5,$ $Y=1.5,$ and
$Y=1.98$ and compare the numerical results with those of the COS and CONV
methods. For the cases of $Y=0.5,$ and $Y=1.5,$ we compute the call reference
values for the numerical experiments using the COS method with $N=2^{14}.$
However, for the case of $Y=1.98,$ we use the SFP method with $U=2^{9}$ to
generate the put reference values. We use $L=10$ in the truncated interval of
(85), as all the input PDFs are smooth. In Figure 9, the recovered density
functions for the three cases are plotted. As when $Y$ tends to 2, the tails
of the PDF are fatter and heavier, and the centre of the distribution shifts.
Table 9 replicates Table 7 in Fang and Oosterlee (2009a). The table compares
our method with the CONV and COS methods, and we see that the SFP method can
achieve better accuracy than the other methods. For example, in the SFP
method, when $U$ equals 32 terms, the error difference from the true values
are $2.608e-08$ (when Y=0.5) and $5.060e-10$ (when Y=1.5). The results are far
better than those of other two methods with the same number of terms/grids.
Figure 10 is the graphical result for the test of $Y=1.98.$ In the test, we
first generate 250 different put prices in a range of strikes from 80 to 120
and check for the error convergence against the put prices in the range. In
Figure 10, we can clearly see that the SFP method obtains global spectral
convergence and fairly high accuracy.

Figure 9: Recovered density functions (top) and their zoom (bottom) for the CGMY model generated via the parameters CGMY–Para1.  Table 9: Comparison of the CONV, COS and SFP methods in terms of absolute error for pricing vanilla call options under the CGMY model with $K=100$. The reference values are $19.812948843118576$ ($Y=0.5$) and $49.790905468523860$ ($Y=1.5$). The SFP method is more accurate than others with fewer summation terms required. The parameters are taken from CGMY–Para1. $Y=0.5$ | CONV | COS | SFP | $Y=1.5$ | CONV | COS | SFP
---|---|---|---|---|---|---|---
$N/U$ | Error | Error | Error | $N/U$ | Error | Error | Error
8 | 2.320 | 4.443 | 2.984e-02 | 8 | 5.246 | 9.304e-01 | 3.020e-02
16 | 1.079 | 5.283e-01 | 3.163e-04 | 16 | 7.763e-01 | 2.863e-02 | 2.566e-05
32 | 8.172e-01 | 1.240e-02 | 2.608e-08 | 32 | 7.607e-01 | 1.240e-05 | 5.060e-10
64 | 2.089e-01 | 2.801e-05 | 7.687e-11 | 48 | 1.384 | 5.286e-12 | 8.527e-14
Figure 10: Error envelopes of the SFP methods for pricing European at/around-
the-money put options under the CGMY model in a range of $K$ from 80 to 120
when $Y=1.98.$ The number of terms increases in a sequence of 8 (red line), 16
(blue line), 32 (green line) and 64 (black line). Spectral convergence is
observed in the SFP method. The parameters are taken from CGMY–Para1.

###  8.2 Affine Processes

Finally, for affine processes, we choose the Heston model and price calls
using the following parameters:

$\displaystyle\small\textbf{Heston--Para1}:$ $\displaystyle
S=100,\,y_{0}=0.0175,\,\overline{y}=0.0398,\,\lambda=1.5768,$
$\displaystyle\eta=0.5751,\,\rho=-0.5711,\,r=0,\,T=1\mbox{ or }10\mbox{ or
}30\mbox{ or }45,\,q=0.$ (116)

Figure 11 presents the recovered density functions. It shows that $T=1$ gives
rise to a sharper peaked density than $T=10,$ as expected. In this test, we
compare the SFP method with the COS method. Most of the reference values are
generated via the COS method with $N=2^{14}.$ However, when the option value
is measured at the jump in the PDF, we generate the value using the SFP
method, as the reference value provided by other methods is not accurate. We
use $L=12$ in the truncated intervals of (85) and (87). We replicate the
numerical tests shown in Tables 4 and 5 of Fang and Oosterlee (2009a) and
display them in Table 11. In the original test, Fang and Oosterlee (2009a)
(see Appendix 11, Table 14) report that the COS method has an algebraic
convergence rate in $T=1$ and spectral convergence in $T=10$ when they
approximate the call values with a strike of 100. They also suggest that the
error convergence is very reasonable, e.g., the error difference is $3.17e-07$
with $N=192$ in the test of $T=1$ and $1.85e-10$ with $N=160$ in the test of
$T=10.$ However, this seems not to be the case when we replicate these tests.
Using Heston–Para1 and the idea of the first and second cumulants to construct
the truncated intervals (as suggested in Fang and Oosterlee, 2009a) for the
Heston model, the PDF range is fairly large at $[-11,11].$ According to Fang
and Oosterlee (2009a), this implies that more terms should be used to
compensate for the error convergence when the truncated interval is larger. It
is not clear how Fang and Oosterlee (2009a) obtain the results in Table 14
without requiring many terms, especially for the test of $T=1$.

Without making any changes, we follow the steps reported in Fang and Oosterlee
(2009a) and Ruijter et al. (2013) to compute the call prices via the COS
method and the filer-COS method, respectively, and compare their error
convergence results with those of the SFP method. In Table 10, the SFP method
provides roughly 8 digits of accuracy in the case of $T=1$ and 10 digits of
accuracy in the case of $T=10$ with only 128 terms required. This is an
immense improvement over the other methods, which provide only 1-2 digits of
accuracy with the same number of terms. Furthermore, if we compare our results
with those shown in Tables 4 and 5 of Fang and Oosterlee (2009a) (see Appendix
11, Table 14), our method yields better convergence. As the PDF is non-smooth
when $T=1,$ we compare the COS, filter-COS and SFP methods in terms of
absolute error for pricing call options around/at the jump. When $K=50,$ the
SFP method has 14 digits of accuracy, but the other two methods gain only 4
digits of accuracy with 256 terms required. When $K=105.453,$ the option price
is at the jump. The call price generated by the SFP can reach 4 digits of
accuracy. This is in line with the finding of Driscoll and Fornberg (2011).
Finally, the last test replicates Table 3 in Ortiz-Gracia and Oosterlee
(2013). In the test, the authors compare the error convergence of the Haar
wavelets method and the COS method when they price call options with long
maturities. They selected maturities $T=30$ and $T=45$, which may correspond
to pensions or mortgage contracts. Table 12 lists the error convergence
results from Table 3 in Ortiz-Gracia and Oosterlee (2013) and compares them
with those produced by the SFP method. As we see in the table, the SFP method
yields better and faster error convergence than the other methods. For
example, the absolute error of the SFP method is roughly $3.000e-06$ for both
$T=30\mbox{ and }45$ with $U=64$, but those of other two methods are range
from $2.46e-03$ to $9.68e-01$ in both cases under the same conditions.

Figure 11: Recovered density functions for the Heston model generated via the parameters Heston–Para1.  Table 10: Comparison of the CONV, COS, filter-COS and SFP methods in terms of absolute error for pricing a call option under the Heston model with $K=100$. The reference values are $5.7851554534076321$ ($T=1$) and $22.318945791154533$ ($T=10$). The SFP method is more accurate than the others with fewer summation terms required. The parameters are taken from Heston–Para1. $T=1$ | COS | filter-COS | SFP | $T=10$ | CONV | COS | SFP
---|---|---|---|---|---|---|---
$N/U$ | Error | Error | Error | $N/U$ | Error | Error | Error
8 | 19.77 | 24.44 | 2.568 | 8 | 21.65 | 23.17 | 5.557
16 | 9.148 | 13.73 | 1.111e-01 | 16 | 15.65 | 18.79 | 4.555e-01
32 | 3.096 | 5.843 | 4.725e-02 | 32 | 7.721 | 12.01 | 2.178e-04
64 | 6.381e-01 | 1.752 | 1.262e-04 | 64 | 1.878 | 4.779 | 3.231e-05
128 | 2.685e-02 | 2.68e-01 | 1.331e-08 | 128 | 7.969e-02 | 7.982e-01 | 7.529e-10
Table 11: Comparison of the COS, filter-COS and SFP methods in terms of
absolute error for pricing a call option under the Heston model with $T=1$.
The reference values are $50.070539139715081$ ($K=50$) and $3.181555642433310$
($K=105.453$). The SFP method is more accurate than others with fewer
summation terms required. The parameters are taken from Heston–Para1 .

$K=50$ | COS | filter-COS | SFP | $K=105.453$ | COS | filter-COS | SFP
---|---|---|---|---|---|---|---
$N/U$ | Error | Error | Error | $N/U$ | Error | Error | Error
8 | 5.031 | 8.546 | 3.514e-01 | 8 | 19.71 | 24.37 | 2.608
16 | 8.763e-01 | 1.202 | 1.047e-02 | 16 | 9.163 | 13.71 | 12.81
32 | 4.861e-01 | 1.197 | 6.327e-06 | 32 | 3.201 | 5.902 | 4.681
64 | 3.715e-02 | 9.661e-02 | 1.861e-08 | 64 | 7.801e-01 | 1.888 | 3.649e-02
128 | 6.003e-03 | 2.361e-02 | 2.467e-11 | 128 | 1.003e-01 | 3.941e-01 | 1.001e-03
256 | 3.922e-04 | 2.363e-04 | 8.527e-14 | 256 | 3.722e-03 | 3.801e-02 | 1.898e-04

Table 12: Comparison of the COS, Haar wavelet and SFP methods in terms of
absolute error for pricing a call option under the Heston model. The reference
value is computed via the COS method using 50000 terms. The SFP method is more
accurate than the others with fewer summation terms required. The parameters
are taken from Heston–Para1.

T=30 | T=45
---|---
Haar | COS | SFP | Haar | COS | SFP
scale | error | N | error | U | error | scale | error | N | error | U | error
3 | 2.58e+02 | 8 | 1.72e+06 | 8 | 3.783 | 3 | 3.92e+02 | 8 | 3.19e+07 | 8 | 2.584
4 | 2.72e+00 | 16 | 2.75e+05 | 16 | 2.259e-01 | 4 | 1.09e+01 | 16 | 5.45e+06 | 16 | 2.953e-01
5 | 7.94e-01 | 32 | 2.19e+03 | 32 | 3.981e-04 | 5 | 5.99e-01 | 32 | 3.10e+04 | 32 | 9.438e-04
6 | 2.46e-03 | 64 | 3.37e-01 | 64 | 1.353e-06 | 6 | 1.05e-02 | 64 | 9.68e-01 | 64 | 3.049e-06

##  9 Conclusion and discussion

In this paper, we illustrate how to use the SFP method (cf. Driscoll and
Fornberg, 2001, 2011) as an alternative approach to hedging and pricing
European-type options. We provide an error analysis to show that a good
truncated interval can yield the highest accuracy of the method and
theoretically prove that the approximate values generated by the method can
trend towards the true option prices. Through the numerical experiment
results, we first show that the SFP method has the ability to price any
European-type option with the features of deep in/out of the money and/or very
long/short maturities. Second, compared with other numerical methods of option
pricing, we show that the SFP method can has faster error convergence with
fewer partial summation terms required, and it remains consistently accurate
when approximating large and small option prices. Third, it does not require a
scale or damping factor to adjust its accuracy, e.g., as in the SWIFT, CONV
and filter-COS methods. Finally, when PDFs are smooth, the method exhibits
global spectral convergence. When PDFs are not smooth, it achieves global
spectral convergence except at the jumps. Nevertheless, the results indicate
that the SFP method offers 4–6 digits of accuracy at the jumps. This is indeed
in line with the findings of Driscoll and Fornberg (2011).

Although the theoretical analysis/numerical results presented here have
demonstrated the effectiveness of the SFP method. Further work might proceed
in three ways. First, the creation of the truncated interval (87) for a non-
smooth PDF relies on a process of trial and error. In the future, a more
theoretical, accurate algorithm for creating truncated intervals should be
used for all non-smooth PDFs. Second, whether the SFP method can be applied to
price European spread options is an interesting research question, as the SFP
worked well in with one dimension (one asset). Finally, our ultimate goal is
to extend the method to price options with early exercise and/or path-
dependant features.

## Acknowledgement

We thank Professor Toby Driscoll, Department of Mathematical Sciences,
University of Delaware and Professor Bengt Fornberg, Department of Applied
Mathematics, University of Colorado for comments that greatly improved the
manuscript.

## 10

In Table 13, we show the first $c_{1},$ second $c_{2},$ and fourth $c_{4}$
cumulants of the BSM, VG, $\rm CGMY$ and Heston models. However, as Fang and
Oosterlee (2009a) suggest, due to the lengthy representation of $c_{4},$ we
only present the first two cumulants of the Heston model. Given the
characteristic functions, the cumulants can be generally computed using

$c_{n}={1\over i^{n}}{\partial^{n}\log\varphi(z)\over\partial
z^{n}}\bigg{|}_{z=0}.$

Table 13: The first $c_{1},$ second $c_{2},$ and fourth $c_{4}$ cumulants of
various models Lévy Models
---
BSM |  $c_{1}=(r-q+\omega)t$ $c_{2}=\sigma^{2}t,$ $c_{4}=0,$ $\omega=-0.5\sigma^{2}$
VG | $c_{1}=(r-q+\theta+\omega)t$
| $c_{2}=(\sigma^{2}+\upsilon\theta^{2})t$
|
$c_{4}=3(\sigma^{4}\upsilon+2\theta^{4}\upsilon^{3}+4\sigma^{2}\theta^{2}\upsilon^{2})t$
| $\omega={1\over\upsilon}\log(1-\theta\upsilon-\sigma^{2}\upsilon/2)$
$\rm CGMY$ | $c_{1}=(r-q+\omega)t$
| $c_{2}=(C\Gamma(2-Y)(M^{Y-2}+G^{Y-2})t$
| $c_{4}=(C\Gamma(4-Y)(M^{Y-4}+G^{Y-4})t$
|
$\omega=\left(C\Gamma(-Y)G^{Y}\left(\left(1+\frac{1}{G}\right)^{Y}-1-\frac{Y}{G}\right)+C\Gamma(-Y)M^{Y}\left(\left(1-\frac{1}{M}\right)^{Y}-1+\frac{Y}{M}\right)\right)$
Affine Processes
Heston | $c_{1}=(r-q)t+(1-e^{-\lambda t}){\bar{y}-y_{0}\over 2\lambda}-0.5\bar{y}t$
| $c_{2}={1\over 8\lambda^{3}}\Big{(}\eta t\lambda\exp(-\lambda
t)(y_{0}-\bar{y})(8\lambda\rho-4\eta)+$
|  $\lambda\rho\eta(1-e^{-\lambda t})(16\bar{y}-8y_{0})+2\bar{y}\lambda
t(-4\rho\eta+\eta^{2}+4\lambda^{2})+$
|  $\eta^{2}((\bar{u}-2u_{0})e^{-2\lambda t}+\bar{y}(6e^{-\lambda
t}-7)+2y_{0})+8\lambda^{2}(y_{0}-\bar{y})(1-e^{-\lambda t})\Big{)}$

##  11 Absolute error convergence for the COS method for calls under the
Heston model

Table 14: Absolute error convergence for the COS method for calls under the Heston model with $T=1$ and $T=10.$ Reference values are 5.785155450…($T=1$) and 22.318945791…($T=10$). $T=1$ |  | $T=10$ |
---|---|---|---
$N$ | Error | $N$ | Error
64 | 4.92e-03 | 32 | 7.40e-03
96 | 2.99e-04 | 64 | 5.02e-05
128 | 1.94e-05 | 96 | 1.40e-07
160 | 2.99e-06 | 128 | 4.92e-10
192 | 3.17e-07 | 160 | 1.85e-10

## References

  * Almendral and Oosterlee (2007) Almendral, A. and Oosterlee, C.W., Accurate Evaluation of European and American Options Under the CGMY Process. SIAM J. Sci. Comput., 2007, 29, 93–117.
  * Andricopoulos et al. (2003) Andricopoulos, A.D., Widdicks, M., Duck, P.W. and Newton, D.P., Universal Option Valuation Using Quadrature Methods. J. Financ. Econ., 2003, 67, 447–471.
  * Black and Scholes (1973) Black, F. and Scholes, M., The Pricing of Options and Corporate Liabilities. J. Political Econ., 1973, 81, 637–654.
  * Carr et al. (2002) Carr, P., Geman, H., Madan, D.B. and Yor, M., The Fine Structure of Asset Returns: An Empirical Investigation. J. Bus., 2002, 75, 305–332.
  * Carr and Madan (1999) Carr, P. and Madan, D., Option Valuation Using the Fast Fourier Transform. J. Comput. Finance, 1999, 4, 61–73.
  * Chan (2016) Chan, T.L.R., An Orthogonal Series Expansions Method to Hedge and Price European-Type Options. Available online at: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2852033, 2016.
  * Chisholm and Common (1981) Chisholm, J.S.R. and Common, A.K., Generalisations of Padé Approximation for Chebyshev and Fourier Series. In E. B. Christoffel: The Influence of His Work on Mathematics and the Physical Sciences, edited by P.L. Butzer and F. Fehér, pp. 212–231, 1981 (Birkhäuser Basel: Basel).
  * Chourdakis (2004) Chourdakis, K., Option Pricing Using the Fractional FFT. J. Comput. Finance, 2004, 8, 1–18.
  * Cont and Tankov (2004) Cont, R. and Tankov, P., Financial Modelling with Jump Processes, Chapman & Hall/CRC Financial Mathematics Series, 2004 (Chapman & Hall/CRC: Boca Raton, FL).
  * Driscoll and Fornberg (2001) Driscoll, T.A. and Fornberg, B., A Padé-based Algorithm for Overcoming the Gibbs phenomenon. Numer. Algorithms, 2001, 26, 77–92.
  * Driscoll and Fornberg (2011) Driscoll, T.A. and Fornberg, B., The Gibbs Phenomenon in Various Representations and Applications, 2011 (Sampling Publishing: Potsdam, NY).
  * Duffie et al. (2003) Duffie, D., Filipovic, D. and Schachermayer, W., Affine Processes and Applications in Finance. Ann. Appl. Probab., 2003, 13, 984–1053.
  * Eckhoff (1997) Eckhoff, K.S., On a High Order Numerical Method for Solving Partial Differential Equations in Complex Geometries. J. Sci. Comput., 1997, 12, 119–138.
  * Fang and Oosterlee (2009a) Fang, F. and Oosterlee, C.W., A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions. SIAM J. Sci. Comput., 2009a, 31, 826–848.
  * Fang and Oosterlee (2009b) Fang, F. and Oosterlee, C.W., Pricing Early-exercise and Discrete Barrier Options by Fourier-Cosine Series Expansions. Numerische Mathematik, 2009b, 114, 27–62.
  * Fang and Oosterlee (2011) Fang, F. and Oosterlee, C.W., A Fourier-Based Valuation Method for Bermudan and Barrier Options Under Heston’s Model. SIAM J. Financ. Math., 2011, 2, 439–463.
  * Geer (1995) Geer, J.F., Rational Trigonometric Approximations Using Fourier Series Partial Sums. J. Sci. Comput., 1995, 10, 325–356.
  * Gonnet et al. (2013) Gonnet, P., Güttel, S. and Trefethen, L.N., Robust Padé Approximation via SVD. SIAM Rev., 2013, 55, 101–117.
  * Gottlieb and Shu (1997) Gottlieb, D. and Shu, C.W., On the Gibbs Phenomenon and Its Resolution. SIAM Rev., 1997, 39, 644–668.
  * Itkin (2005) Itkin, A., Pricing Options With VG Model Using FFT. Available online at: arXiv:physics/0503137v1, 2005.
  * Leentvaar and Oosterlee (2008) Leentvaar, C. and Oosterlee, C., Multi-Asset Option Pricing Using a Parallel Fourier-Based Technique. J. Comput. Finance, 2008, 12, 1–26.
  * Lewis (2001) Lewis, A., A Simple Option Formula for General Jump-Diffusion and other Exponential Lévy Processes. Available online at: http://optioncity.net/pubs/ExpLevy.pdf, 2001.
  * Lipton (2002) Lipton, A., The Vol Smile Problem. Risk, 2002, 15, 61–66.
  * Lord et al. (2008) Lord, R., Fang, F., Bervoets, F. and Oosterlee, C.W., A Fast and Accurate FFT-Based Method for Pricing Early-Exercise Options under Lévy Processes. SIAM J. Sci. Comput., 2008, 30, 1678–1705.
  * Lord and Kahl (2010) Lord, R. and Kahl, C., Complex Logarithms in Heston-Like Models. Math. Finance, 2010, 20, 671–694.
  * Madan et al. (1998) Madan, D.B., Carr, P. and Chang, E.C., The Variance Gamma Process and Option Pricing. Eur. Finance Rev., 1998, 2, 79–105.
  * Madan and Milne (1991) Madan, D.B. and Milne, F., Option Pricing with V. G. Martingale Components. Math. Finance, 1991, 1, 39–55.
  * Madan and Seneta (1990) Madan, D.B. and Seneta, E., The Variance Gamma (V.G.) Model for Share Market Returns. J. Bus., 1990, 63, 511–524.
  * Merton (1973) Merton, R.C., Theory of Rational Option Pricing. Bell J. Econ. Manag. Sci., 1973, 4, 141–183.
  * Ortiz-Gracia and Oosterlee (2013) Ortiz-Gracia, L. and Oosterlee, C.W., Robust Pricing of European Options with Wavelets and the Characteristic Function. SIAM J. Sci. Comput., 2013, 35, B1055–B1084.
  * Ortiz-Gracia and Oosterlee (2016) Ortiz-Gracia, L. and Oosterlee, C.W., A Highly Efficient Shannon Wavelet Inverse Fourier Technique for Pricing European Options. SIAM J. Sci. Comput., 2016, 38, B118–B143.
  * Pistorius and Stolte (2012) Pistorius, M. and Stolte, J., Fast Computation of Vanilla Prices in Time-Changed Models and Implied Volatilities Using Rational Approximations. Int. J. Theor. Appl. Finance, 2012, 15, 1250031.
  * Ruijter et al. (2013) Ruijter, M., Versteegh, M. and Oosterlee, C., On the Application of Spectral Filters in a Fourier Option Pricing Technique. Available online at: http://ta.twi.tudelft.nl/mf/users/oosterle/oosterlee/main_v6.pdf, 2013.
  * Schoutens (2003) Schoutens, W., Lévy Processes in Finance: Pricing Financial Derivatives, Wiley Series in Probability and Mathematical Statistics, 2003 (Wiley: Chichester, UK).
  * Small and Charron (1988) Small, R.D. and Charron, R.J., Continuous and Discrete Nonlinear Approximations Based on Fourier Series. IMA J. Numer. Anal., 1988, 8, 281–293.
  * von Sydow et al. (2015) von Sydow, L., Höök, L.J., Larsson, E., Lindström, E., Milovanović , S., Persson, J., Shcherbakov, V., Shpolyanskiy, Y., Sirén, S., Toivanen, J., Waldén, J., Wiktorsson, M., Levesley, J., Li, J., Oosterlee, C.W., Ruijter, M.J., Toropov, A. and Zhao, Y., BENCHOP–The BENCHmarking Project in Option Pricing. International Journal of Computer Mathematics, 2015, 92, 2361–2379.
  * Wang et al. (2007) Wang, I.R., Wan, J.W.L. and Forsyth, P.A., Robust Numerical Valuation of European and American Options Under the CGMY Process. J. Comput. Finance, 2007, 10, 31–70.
  * Zhang and Oosterlee (2013) Zhang, B. and Oosterlee, C.W., Efficient Pricing of European-Style Asian Options under Exponential Lévy Processes Based on Fourier Cosine Expansions. SIAM J. Financ. Math., 2013, 4, 399–426.