Tat Lung (Ron) Chan

Hedging and Pricing European-type, Early-Exercise and Discrete Barrier Options
using an Algorithm for the Convolution of Legendre Series Tat Lung (Ron) Chan
School of Business, University of East London, Water Lane, Stratford, UK, E15
4LZ<EMAIL_ADDRESS>Nicholas Hale Department of Mathematical Sciences,
Stellenbosch University, Stellenbosch, South Africa<EMAIL_ADDRESS>
This paper applies an algorithm for the convolution of compactly supported
Legendre series (the CONLeg method) (cf. Hale and Townsend 2014a), to
pricing/hedging European-type, early-exercise and discrete-monitored barrier
options under a Lévy process. The paper employs Chebfun (cf. Trefethen et al.
2014) in computational finance and provides a quadrature-free approach by
applying the Chebyshev series in financial modelling. A significant advantage
of using the CONLeg method is to formulate option pricing and option Greek
curves rather than individual prices/values. Moreover, the CONLeg method can
yield high accuracy in option pricing and hedging when the risk-free smooth
probability density function (PDF) is smooth/non-smooth. Finally, we show that
our method can accurately price/hedge options deep in/out of the money and
with very long/short maturities. Compared with existing techniques, the CONLeg
method performs either favourably or comparably in numerical experiments.

Convolution, Legendre series, European options, early-exercise options,
discrete-monitored barrier options, Lévy process This paper was first
submitted on the 12th November, 2018.

##  1 Introduction

Applying robust numerical techniques in option pricing/hedging and model
calibration provides interesting research questions in financial markets. The
techniques must be not only highly accurate, but also efficient.

Suppose we consider the well-known European vanilla option pricing formula
driven by a stochastic stock price process $(S_{t})_{t\geq 0}:$

$\displaystyle\small
V(x,K,t)=e^{-r(T-t)}\mathds{E}(U(S_{T},K)|S_{t}=e^{x})=e^{-r(T-t)}\int_{-\infty}^{+\infty}U(e^{x+\chi},K)g(\chi)\mathrm{d}\chi,$
(1)

where $V$ denotes the option value at an initial date of $t$ and with a strike
price of $K,$ $U$ stands for a payoff function at the maturity of $T,$
$\mathds{E}$ is the expectation operator under the risk-neutral measure, $x$
and $\chi$ are the log-price and state variable respectively, $S_{T}$ can be
decomposed into $S_{t}=e^{x}$ and $e^{\chi},$ $g$ is the probability density
function of the process, and finally $r$ is a risk-neutral interest rate.

In the aforementioned formula, $V$ can be seen as a convolution integral, more
precisely a cross-correlation integral, and the fast Fourier transform method
(FFT method), a numerical integration-based method, (e.g., Carr and Madan
1999, Lewis 2001, Lipton 2002, Chourdakis 2004, Jackson et al. 2008, Lord et
al. 2008), is a popular method for pricing European vanilla options as well as
more exotic options, such as the American option, under Lévy processes. This
is because the characteristic function of the underlying dynamics can be
easily transformed into a risk-free probability density function (PDF) via the
FFT method. The seminal work of these papers leads to the extension of
combining the FFT with other transformation methods, e.g., the Hilbert
transform or Gaussian transform, in pricing exotic options under the (time-
changed) Lévy process or stochastic volatility models (e.g., Broadie and
Yamamoto 2003, 2005, Feng and Linetsky 2008, Cai and Kou 2011, Wong and Guan
2011, Zeng and Kwok 2014). Within the same numerical integration-based
framework, numerical quadrature methods, such as Gaussian quadrature, are
proposed mainly by Andricopoulos et al. (2003), O’Sullivan (2005),
Andricopoulos et al. (2007), Chen et al. (2014) and Su et al. (2017). The
authors of these papers abbreviate their quadrature techniques as QUAD
methods. The original QUAD method was introduced in Andricopoulos et al.
(2003, 2007) and requires the transition density to be known in closed form,
which is the case in, e.g., the Black-Scholes model and Merton’s jump-
diffusion model. This requirement is relaxed in Chen et al. (2014) where the
QUAD-FFT or QUAD-CONV method is proposed. The main idea is that the PDF can be
recovered by inverting the characteristic function via the FFT method. This
helps open up the QUAD method to a much wider range of models. The latest
development of the method (Su et al. 2017) is to improve the calculation
speed, precomputing and caching PDFs and then applying the extrapolation and
smoothing techniques to the method. Under the general framework of the QUAD
method, only Delta, an option Greek, is formulated via the first-order finite
difference method (FD) (Andricopoulos et al. 2003). However, the accuracy of
the first-order FD is debatable in the context. Furthermore, other kinds of
the Greeks, e.g. Gamma or Theta, are not mentioned or developed yet in this
literature. In addition to the QUAD methods, in recent years, Pachón (2018)
has introduced the CHEB method, Clenshaw-Curtis quadrature based on an
expansion of the integrand in terms of Chebyshev polynomials, for approximate
European options with arbitrary payoffs. The method is one of the natural
applications of Chebfun (Trefethen et al. 2014), an open-source software
system for numerical computing with functions.

Beyond the FFT method and the QUAD method, Oosterlee and his collaborators
have attracted considerable attention (Leentvaar and Oosterlee 2008, Fang and
Oosterlee 2009a, b, 2011, Zhang and Oosterlee 2013, Ruijter et al. 2015). In
their work, they adopt the Fourier cosine series (COS) to price options or
derivatives that have different contingency claims and are characterised by
path dependence and/or early-exercise features. The implementation of these
methods is relatively simple but elegant and is capable of pricing options
under different stochastic processes as long as their characteristic function
exists. The main achievement of these methods is that they can, in many cases
(such as European options), maintain an exponential convergence rate when
pricing options. Moreover, these methods are also able to accurately price
options under infinite variation processes. The COS method requires an
adequate computational domain a priori, and, due to the recursion in time,
errors caused by an inadequate domain propagate, resulting in incorrect option
prices in the COS method. Moreover, based on the framework of the COS method,
Oosterlee and his collaborators further apply the wavelet method called the
Shannon wavelet inverse Fourier technique (SWIFT) method to price European,
spread, path-dependent and discrete barrier options under exponential Lévy
dynamics (e.g., Ortiz-Gracia and Oosterlee 2013, 2016, Colldeforns-Papiol et
al. 2017). The SWIFT method is used to circumvent the ineffectiveness of using
the COS method to price early-exercise options. Nevertheless, comparing with
the COS method, a theoretical proof is still lacking to show exponential
convergence for the SWIFT method since the wavelet scale $m,$ a parameter to
adjust the accuracy of the method, is still chosen heuristically for achieving
exponential convergence. Furthermore, the literature on SWIFT methods does not
consider computation of Greeks.

In this manuscript we propose the CONLeg method, which uses Chebyshev and
Legendre series to approximate a convolution, to improve the ineffectiveness
of using the aforementioned methods. First, as mentioned in Fang and Oosterlee
(2009a, b) and Chan (2018), the FFT-driven methods, like the CONV method (Lord
et al. 2008), are computationally expensive to price option prices or to
approximate the PDFs because a relatively large number of Fourier terms is
required to obtain sufficient accuracy (see Table 5 in Section 6). Second, we
question whether any kind of quadrature method is an effective method to
approximate option prices/hedging values when option/hedging pricing formulae
are treated as convolution integrals. Hale and Townsend (Section 2 and Section
6, 2014b) suggest that since a convolution integral appears as a trapezoid
(cf. Fig 1), more quadrature weights and abscissae are required to approximate
the areas towards the left and right vertices of the trapezoid. Third, we want
to provide a method that can work with any stochastic process with or without
a closed-form PDF. For example, the COS and SWIFT methods only work very well
when the process has a characteristic function and without applying the FFT,
the QUAD method only works well with a closed-form PDF. Fourth, comparing
against the CONV, COS, QUAD and SWIFT methods, the CONLeg method provides an
option pricing/hedging curve rather than an individual point value. Fifth,
unlike the COS method, we can require a proper computational domain a priori
in the CONLeg method when we apply it to price early-exercise options. Hence,
the CONLeg method will not produce incorrect option prices when we calculate
the option prices recursively backwards in time, as we allow a sufficient
computational domain for calculation. Sixth, through our numerical
experiments, we show that the CONLeg method can achieve high accuracy in
option pricing and hedging when the risk-free smooth probability density
function (PDF) is non-smooth, and can be an effective method to price/hedge
options for long/short maturities. Seventh, to fit in the Chebfun framework,
we lay out a closed-form transformation of the complex Fourier expression of a
smooth PDF into Chebyshev series when the closed-form PDF is not available
(cf. Appendix 8). In the similar fashion, we also provide a solution, the
Fourier–Padé method (cf. Appendix 9), for approximating non-smooth PDFs.
Eighth, through this paper, we are keen to promote the convenience of Chebfun
(Trefethen et al. 2014) in financial modelling. Chebfun is a robust, open
source MATLAB pack for computing with functions to 15 digits of accuracy. It
contains several state of the art algorithms for Chebyshev and other
orthogonal polynomials. Finally, the current CONLeg method is different to the
previous literature prompted Chebyshev series and interpolants in option
pricing/hedging (cf. Gaß et al. 2018, Pachón 2018). As we have mentioned in
the previous point, comparing with the CHEB method (Pachón 2018), the CONLeg
method is quadrature-free and not limited in pricing/hedging European-type
options. Also, unlike Gaß et al. (2018)’s tensorized Chebyshev interpolation
to computing Parametric Option Prices (POP), our method does not approximate
option pricing curves which must be first precomputed by any numerical method,
such as the Monte Carlo and the FFT methods, for some fixed parameter
configurations, and then compute other option prices for arbitrary parameter
constellations. On the contrary, for a set of fixed parameters, the CONLeg
method directly generates an option pricing/hedging curve via approximating a
convolution integral without applying other numerical methods for computing
the integral first.

The remainder of this paper is structured as follows. Section 2 describes the
algorithm for the convolution of Legendre series and how Legendre series can
be efficiently computed using the relationship with Chebyshev series and the
FFT. Section 3 introduces the financial stochastic models we examine in this
paper. Section 4 describes the formulation of the CONLeg option pricing/Greek
formulae for different styles of European options as well as Bermuda,
American, and discrete monitored barrier options. Section 5 describes the
choice of truncated integration intervals. Section 6 discusses, analyses, and
compares the numerical results of the CONLeg method with those of other
numerical methods discussed above. Finally, we conclude and discuss possible
future developments in Section 7.

##  2 Convolution of Legendre Series

Convolution is a fundamental operation that arises in many fields,
particularly in financial derivatives research (cf. Carr and Madan 1999, Lewis
2001, Lipton 2002, Jackson et al. 2008, Lord et al. 2008), econometrics
(Bondarenko 2003, Liu et al. 2016) and statistics (Hogg et al. 2004). Given
two integrable functions, $f$ and $g$, their convolution is a third function,
$h$, defined formally by the integral

$\displaystyle h(x)=(f\ast
g)(x)=\int_{-\infty}^{+\infty}f(y)g(x-y)\mathrm{d}y.$ (2)

In general, if both $f$ and $g$ are analytic (smooth) and periodic functions,
the FFT method, which utilises the convolution theorem and the fast Fourier
transform (FFT), is the best choice for approximating $h$ as the FFT
approximations of $f$ and $g$ do not suffer the Gibbs phenomenon.

If we now consider $f$ and $g$ $:[c,d]\rightarrow\mathds{R}$ as two compactly
supported outside of $[c,d],$ then the convolution aforementioned $h=f\ast g$
via the integral is given by

$\displaystyle h(x)=(f\ast
g)(x)=\int_{\min(d,x-c)}^{\max(c,x-d)}f(y)g(x-y)\mathrm{d}y,\quad
x\in[2c,2d],$ (3)

and $h(x)=0$ for $x\notin[2c,2d]$ 111Since convolution is a commutative
operation, we consider that only $f$ and $g$ are in the same intervals..
Without losing any generality, we can visualise each value of $x$ by the
diagram in Figure 1, and we split $h$ into the two pieces suggested by the
diagram:

$\displaystyle
h(x)=\begin{cases}h^{L}(x)=\bigintsss_{c}^{x-c}f(y)g(y-x)\mathrm{d}y\quad\quad&x\in[2c,c+d],\\\
&\\\
h^{R}(x)=\bigintsss_{x-d}^{d}f(y)g(y-x)\mathrm{d}y\quad\quad&x\in[c+d,2d].\end{cases}$
(4)

Here, we denote, once and for all, $L$ and $R$ as the left and right hand side
of the convolution, respectively.

$c+d$$2c$$2d$$c$$d$$y$$x$ Figure 1: The convolution domain for two Legendre
series on $[c,d].$

Based in (3), either $f$ or $g$ can be a non-periodic continuous function, the
FFT approximations of $f$ and $g$ suffer the Gibbs phenomenon such that there
is a permanent oscillatory overshoot in the neighbourhoods of the endpoints
$c$ and $d.$ Accordingly, to avoid the Gibbs phenomenon, we adopt the
algorithm proposed by Hale and Townsend (2014a) to approximate $h.$ The
crucial idea of the algorithm is to approximate $f$ and $g$ with finite
Legendre series and then convolve the approximations using the convolution
theorem for them. The result is a piecewise polynomial representation that can
be evaluated at any $x$ in the domain of $h$ to yield an approximation to
$h(x)$. If the polynomials used to approximate $f$ and $g$ have degree at most
$N$, their algorithm produces an approximation to $h$ in $\mathcal{O}(N^{2})$
operations. We summarise the approach below.

To illustrate Hale and Townsends’ algorithm of a convolution of two Legendre
series, we first define the Legendre series on $[-1,1]$ and then generalise to
intervals $[c,d].$ Legendre polynomials, invented by Adrien-Marie Legendre,
are the polynomial solutions $P_{n}(x)$ to Legendre’s differential equation

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}x}\left[\left(1-x^{2}\right)\frac{\mathrm{d}P_{n}(x)}{\mathrm{d}x}\right]+n(n+1),\quad
x\in[-1,1],$ (5)

with $P_{0}(x)=0,$ $P_{1}(x)=0$ and integer parameter $n\geq 0.$ $P_{n}(x)$
forms a polynomial sequence of orthogonal polynomials of degree $n$ and it can
be expressed through Rodrigues’ formula:

$\displaystyle
P_{n}(x)=\frac{1}{2^{n}n!}\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}\left(x^{2}-1\right)^{n}.$
(6)

Suppose we have $f_{M}$ and $g_{N}$ as two finite Legendre series on $[-1,1]$
of degrees $M$ and $N$ with coefficients of $\alpha_{0},\ldots,\alpha_{M}$ and
$\beta_{0},\ldots\beta_{N},$ then we may write $f_{M}$ and $g_{N}$ as

$\displaystyle f_{M}(x)=\sum\limits_{m=0}^{M}\alpha_{m}P_{m}(x),\quad
g_{N}(x)=\sum\limits_{n=0}^{N}\beta_{n}P_{N}(x).$ (7)

With $[c,d]$ equal to $[-1,1],$ the convolution in $(\ref{eqn:conv1})$ becomes

$\displaystyle h(x)=(f_{M}\ast
g_{N})(x)=\int_{\max(-1,x-1)}^{\min(1,x+1)}f_{M}(y)g_{N}(x-y)\mathrm{d}y,\quad
x\in[-2,2].$ (8)

From (4), $h$, consists of two pieces, $h^{L}$ on the left with $[-2,0]$ and
$h^{R}$ on the right with $[0,2]$, each of degree $N+M+1.$ We construct
$h^{L}$ and $h^{R}$ by computing their Legendre coefficients. We focus on the
computation of $h^{L}$ since that of $h^{R}$ is similar.

Denote by $\\{\gamma_{k}^{L}\\}_{k}^{M+N+1}$ the vector of the Legendre
coefficients of $h^{L},$ such that

$\displaystyle
h^{L}(x)=\int_{-1}^{x+1}f_{M}(y)g_{N}(x-y)\mathrm{d}y=\sum\limits_{k=0}^{M+N+1}\gamma_{k}^{L}P_{k}(x+1),\quad
x\in[-2,0].$ (9)

By the orthogonality of Legendre polynomials and the orthonormalisation
constant $(k+1/2)^{-1/2}$ for $P_{k}(x)$ for $k=0,\ldots,M+N+1,$ we have for
$P_{k}(x)$ for $k=0,...,M+N+1,$ that

$\displaystyle\gamma_{k}^{L}$ $\displaystyle={2k+1\over
2}\int_{-2}^{0}P_{k}(x+1)\int_{-2}^{0}P_{k}(x+1)\int_{-1}^{x+1}f_{M}(y)g_{N}(x-y)\mathrm{d}y\mathrm{d}x$
(10) $\displaystyle=\sum_{n=0}^{N}\beta_{n}\underbrace{\left[{2k+1\over
2}\sum_{m=0}^{M}\alpha_{m}\int_{-2}^{0}P_{k}(x+1)\int_{-1}^{x+1}P_{m}(y)P_{n}(x-y)\mathrm{d}y\mathrm{d}x\right]}_{=B_{k,n}^{L}}.$
(11)

Hale and Townsend (2014a) prove that the above relation can be expressed in
matrix form as $\underline{\gamma}=B^{L}\underline{\beta}.$ Importantly (Hale
and Townsend 2014a, Theorem 4.1), there is a three-term recurrence relation of
$B_{k,n}^{L}$ such that

$\displaystyle B_{k,n}^{L}$ $\displaystyle=-{2n+1\over
2k+3}B_{k,n+1}^{L}+{2n+1\over 2k-1}B_{k,n+1}^{L}+B_{k,n-1}^{L},\quad n,k\geq
1,$ (12) $\displaystyle B_{k,1}^{L}$
$\displaystyle=\begin{cases}B_{k-1,0}^{L}/(2k-1)-B_{k,0}^{L}-B_{k+1,0}^{L}/(2k+3),\quad
k\neq 0,\\\ -B_{1,0}^{L}/3,\quad k=0\end{cases},$ (13) $\displaystyle
B_{k,0}^{L}$ $\displaystyle=\begin{cases}{\alpha_{k-1}\over
2k-1}-{\alpha_{k+1}\over 2k+3},\quad k\neq 0,\\\ \alpha_{0}-\alpha_{1}/3\quad
k=0\end{cases},$ (14)

where $0\leq k\leq M+N+1$ and $0\leq n\leq N.$ both $B^{L}_{:,0}$ and
$B^{L}_{:,1}$ can be computed in $\mathcal{O}(M+N)$ operations, and the whole
$(M+N)\times N$ matrix, $B^{L}$, in $\mathcal{O}((M+N)N)$ operations. The
matrix-vector product $B^{L}\underline{\beta}$ can be computed with the same
cost and, accordingly, the coefficients $\gamma^{L}$ of $h^{L}$ in
$\mathcal{O}((M+N)N)$ operations. The coefficients $\gamma^{R}$ of $h^{R}$ can
be computed from $B^{R}\underline{\beta}$, for which a nearly identical
recurrence relation can be derived left since the computation for $h^{R}$ is
similar.

Now, we direct our attention to the convolution of two finite Legendre series
defined on the same interval $[c,d].$ We can define the composition of
$P_{k}\circ\psi_{[c,d]},$ where $\psi_{[c,d]}(x)=(2x-(d+c))/(d-c)$ is the
linear mapping from $[c,d]$ to $[-1,1].$ Apart from this, the Legendre series
of $f_{M}$ and $g_{N}$ on $[c,d],$ respectively, are formulated by

$\displaystyle
f_{M}(x)=\sum\limits_{m=0}^{M}\alpha_{m}P_{m}\circ\psi_{[c,d]}(x),\quad
g_{N}(x)=\sum\limits_{n=0}^{N}\alpha_{n}P_{n}\circ\psi_{[c,d]}(x)$ (15)

(Hale and Townsend 2014a, Lemma 4.2). The convolution of $(f\ast g)(x)$ of two
continuous functions of $f$ and $g$ defined on $[c,d]$ can be computed as

$\displaystyle(f\ast
g)(x)=\int_{\min(d,x-c)}^{\max(c,x-d)}f(y)g(x-y)\mathrm{d}y={d-c\over
2}\left((f\circ\psi^{-1}_{[c,d]})\ast(g\circ\psi^{-1}_{[c,d]})\right)(y),$
(16)

where $x\in[2c,2d]$ and $y=2\psi_{[2c,2d]}(x)\in[-2,2].$

The approach above is not restricted to $f$ and $g$ on the same intervals. In
contrast, the algorithm allows $f$ and $g$ which are polynomials of finite
degree with zero support outside of $[a,b]$ and $[c,d],$ respectively. In this
case (3) becomes

$\displaystyle h(x)=(f\ast
g)(x)=\int_{\min(b,x-c)}^{\max(a,x-d)}f(y)g(x-y)\mathrm{d}y,\quad
x\in[a+c,b+d].$ (17)

Depending on the difference between the length of $[a,b]$ and $[c,d],$ these
are three different ideas to compute (17) on these subintervals of
$[a+c,b+d]$. For the full detail of computing $h$ on these subintervals, we
refer interested readers to Hale and Townsend (2014a, Section 5). The
algorithm of convolution of Legendre series supported on same/general
intervals is fully implemented as conv in Chebfun222The algorithm in ${\tt
conv}$ actually computes the convolution between two Chebyshev series by using
fast Chebyshev-Legendre transform (Townsend et al. 2018) implemented in
cheb2leg..

##  3 Lévy Processes

In this section, we briefly introduce the important properties of one-
dimensional Lévy processes. Standard references for the stochastic processes
can be found in Schoutens (2003) and Cont and Tankov (2004). Since markets are
frictionless and have no arbitrage, we assume that an equivalent martingale
measure (EMM) $\mathds{Q}$ is chosen by the market. Moreover, there is a
complete filtered probability space
$\left(\Omega,\mathcal{F},\\{\mathcal{F}\\}_{t\geq 0},\mathds{Q}\right)$ on
which all processes are assumed to live.

###  3.1 Lévy processes

With $r\geq 0$ and $q\geq 0$ as the constant risk-free interest rate and the
constant dividend yield, respectively, we describe a stock process
$(S_{t})_{t\geq 0}$ driven by an exponential Lévy process $(X_{t})_{t\geq 0}$
such that

$S_{t}=S_{0}e^{X_{t}},$

where $X_{0}=0$ and $X_{t}$ has infinitely divisible marginal distributions.
Given a random variable $X_{t},$ we can define it as the corresponding
characteristic function as follows:

$\displaystyle\varphi(u)$
$\displaystyle=\mathbb{E}[e^{iuX_{t}}]=e^{t\phi(u)},\quad u\in\mathbb{R}.$
(18)

If we define a truncation function $h(\chi)=\chi\mathds{1}_{|\chi|\leq 1}$
which is a measurable function such that for every $u\in\mathbb{R},$
$\int|1-e^{iu\chi}+iuh(\chi)|\nu(\mathrm{d}\chi)<\infty,$ the characteristic
function of $X_{t}$ can be described by the Lévy–Khinchine representation such
that

$\displaystyle\phi(u)=iu(r-q+\omega)t-{1\over
2}\sigma^{2}u^{2}+\int_{-\infty}^{+\infty}\left(e^{iu\chi}-1-iuh(\chi)\right)\nu(\mathrm{d}\chi),\quad\chi\in
X_{T-t},$ (19)

Here, $\sigma^{2}\geq 0$ and $\nu$ are Lévy measures on $[-\infty,\infty]$
which do not depend on the choice of $h$ (but note that $r-q+\omega$ depends
on the choice of it). The condition that $(S_{t}e^{-(r-q)t})_{t\geq 0}$ is a
martingale will be guaranteed as long as an appropriate choice of the mean-
correcting compensator $\omega$ is calculated as follows:

$\displaystyle\omega$ $\displaystyle={1\over t}\log\varphi(-i)-(r-q).$ (20)

There are a substantial number of Lévy process examples in financial
modelling. In this paper, we focus on geometric Brownian motion (GBM),
variance gamma (VG), normal inverse Gaussian (NIG) and the
Carr–Geman–Madan–Yor process (CGMY). Their characteristic functions
$\varphi(u)$ can be defined as follows:

$\displaystyle\varphi_{\rm
GBM}(u):=\exp\left(t\left(iu(r-q+\omega)-\frac{1}{2}\sigma^{2}u^{2}\right)\right),$
(21) $\displaystyle\varphi_{\rm
NIG}(u):=\exp\bigg{(}t\bigg{(}iu(r-q+\omega)-\frac{1}{2}\sigma^{2}u^{2}+\delta\left(\sqrt{\alpha^{2}-\beta^{2}}-\sqrt{\alpha^{2}-(\beta+iu)^{2}}\right)\bigg{)}\bigg{)},$
(22) $\displaystyle\varphi_{\rm
VG}(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}},$ (23)
$\displaystyle\varphi_{\rm
CGMY}(u):=\exp\Bigg{(}t\Bigg{(}iu(r-q+\omega)+C\Gamma(-Y)G^{Y}\left(\left(1+\frac{iz}{G}\right)^{Y}-1-\frac{izY}{G}\right)$
$\displaystyle\quad+C\Gamma(-Y)M^{Y}\left(\left(1-\frac{iz}{M}\right)^{Y}-1+\frac{izY}{M}\right)\Bigg{)}\Bigg{)}.$
(24)

###### Remark 3.1

Most Lévy processes do not have a closed-form representation of their risk-
natural probability density function (PDF) $g$; however, some of the
processes, e.g. GBM, NIG and VG, have a closed-formed PDF given by:

$\displaystyle g(x)$ $\displaystyle={1\over\sqrt{2\pi t}\sigma}e^{-{1\over
2}{(x-(r-q+\omega_{GBM}))^{2}\over\sigma^{2}t}},$ (25) $\displaystyle g(x)$
$\displaystyle=\frac{\alpha\delta
tK_{1}\left(\alpha\sqrt{\delta^{2}t^{2}+(x-(r-q+\omega_{NIG})t)^{2}}\right)}{\pi\sqrt{\delta^{2}t^{2}+(x-(r-q+\omega_{NIG})t)^{2}}}\;e^{\delta
t\sqrt{\alpha^{2}-\beta^{2}}+\beta(x-(r-q+\omega_{NIG})t)},$ (26)
$\displaystyle g(x)$ $\displaystyle={2e^{\theta
Z(x)/\sigma^{2}}\over\nu^{t\over\nu}\sqrt{2\pi}\sigma\Gamma(t/\nu)}\left({Z(x)^{2}\over
2\sigma^{2}/\nu+\theta.^{2})}\right)^{{t\over 2\nu}-{1\over 4}}\times
K_{{t\over\nu}-{1\over
2}}\left({1\over\sigma^{2}}\sqrt{Z(x)^{2}\left({2\sigma^{2}\over\nu}+\theta^{2}\right)}\right)$
(27)

respectively. Here, $\omega_{GBM}=-1/2\sigma^{2},$
$\omega_{NIG}=\delta\left(\sqrt{\alpha^{2}-(\beta+1)^{2}}-\sqrt{\alpha^{2}-\beta^{2}}\right),$
$K(\cdot)$ is a modified Bessel function of a second kind, $\Gamma(\cdot)$ is
a gamma function, $Z(x)=x-(r-q+\omega_{VG})t$, and
$\omega_{VG}=-1/\nu\log\left(1-\theta\nu-{\sigma^{2}\nu\over 2}\right).$

##  4 Pricing and Hedging Option via Convolution of Legendre Series

In this section, we apply the algorithm for the convolution of Legendre series
to formulate option pricing/hedging formulas.

###  4.1 Pricing Formulae for European Type Options

Given the current log price $x:=\log S,$ the strike price of $K$ and maturity
$T\geq t,$ and the probability density function (PDF) $g$ of a stochastic
process, we can express the option price $V(x,K,t)$ starting at time $t$ with
its contingent claim paying out $U(S_{T},K)$ as follows:

$\displaystyle V(x,K,t)$
$\displaystyle=e^{-r(T-t)}\mathbb{E}(U(S_{T},K)|S_{t}=e^{x})$
$\displaystyle=e^{-r(T-t)}\mathbb{E}(U(S_{t}e^{X_{T}-X_{t}},K))$
$\displaystyle=e^{-r(T-t)}\int_{-\infty}^{+\infty}U(e^{x+\chi-\log
K},K)g(\chi)d\chi,\quad\chi\in X_{T}-X_{t}=X_{T-t}.$ (28)

By replacing $x+\chi-\log K$ with $y$, we have

$\displaystyle V(x,K,t)$
$\displaystyle=e^{-r(T-t)}\int_{-\infty}^{+\infty}U(e^{y},K)g\left(y-x+\log
K\right)\mathrm{d}y$
$\displaystyle=e^{-r(T-t)}K\int_{-\infty}^{+\infty}f(y)g^{R}\left(\tilde{x}-y\right)\mathrm{d}y,$
(29)

where $\tilde{x}=x-\log K,$ $Kf(y):=U(e^{y},K)$ is the pay-off in log-price
coordinates and $g^{R}(\tilde{x}):=g(-\tilde{x})$ is the reflected PDF
function.

If $f(y)$ is a piecewise continuous function, as is standard for most options,
e.g., vanilla put and call, then applying the FFT method for a convolution
function333The expression of
$\int\limits_{-\infty}^{+\infty}U(e^{y},K)g\left(y-x+\log K\right)\mathrm{d}y$
in (4.1) is indeed a cross-correlation integral; however, since we introduce
the idea of the reflected function $g^{R}(\tilde{x}):=g(-\tilde{x})$, we can
turn (4.1) into a convolution integral instead.,
$(f*g^{R})(\tilde{x}):=\int_{-\infty}^{+\infty}f(y)g^{R}(\tilde{x}-y)\mathrm{d}y,$
in $[-\infty,\infty]$ will cause Gibbs phenomenon and as a result, will affect
the accuracy of approximating $V$. To avoid Gibbs phenomenon and allow a good
approximation of $V,$ we replace $[-\infty,\infty]$ with an interval $[c,d]$
and employ the CONLeg method. The choice of $[c,d]$ satisfies the condition of

$\displaystyle\int_{c}^{d}g(\chi)e^{iu\chi}\rm{d}\chi\approx\int_{-\infty}^{+\infty}g(\chi)e^{iu\chi}\rm{d}\chi=\mathds{E}[e^{iu(X_{T}-X_{t})}]:=\varphi(u),$
(30)

where $\varphi(u)$ is a characteristic function of $X_{T}-X_{t}.$ Then we can
approximate the pricing formula $V$ on $[c,d],$ i.e.,

$\displaystyle V(x,K,t)$ $\displaystyle\approx
e^{-r(T-t)}K\int_{c}^{d}f(y)g^{R}(\tilde{x}-y)\mathrm{d}y.$
$\displaystyle=e^{-r(T-t)}K\int_{\min(d,\tilde{x}-c)}^{\max(c,\tilde{x}-d)}f(y)g^{R}(\tilde{x}-y)\mathrm{d}y,\quad\tilde{x}\in[2c,2d]$
$\displaystyle=e^{-r(T-t)}Kh(\tilde{x}).$ (31)

where, $h(\tilde{x})$ has compact support outside of $[2c,2d].$ The final form
of (4.1) is ready for approximating via the CONLeg method.

Most of the closed-form expressions of $g$ do not exist in the stochastic
processes. If they do not, we adopt the ideas proposed in Chan (2016, 2018) to
express $g$ in a complex Fourier series (CFS) representation such that $g$ is
approximated by:

$\displaystyle
g_{N}(y):=\mathfrak{Re}\left[2\sum_{k=1}^{N}b_{k}e^{i\frac{2\pi}{d-c}ky}+b_{0}\right],$
(32)

where $i$ is a complex number and $\mathfrak{Re}$ is the real part of a
complex number, and given the condition of (30),

$\displaystyle
b_{k}=\int_{c}^{d}g(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
b_{0}=\int_{c}^{d}g(y)\mathrm{d}y\approx\varphi(0)=1.$ (33)

For the expression of $g^{R},$ we simply put a negative sign in the basis
function $e^{i{2\pi\over d-c}ky},$ i.e.,

$\displaystyle g^{R}(y)\approx
g_{N}^{R}(y):=\mathfrak{Re}\left[2\sum_{k=1}^{N}b_{k}e^{-i\frac{2\pi}{d-c}ky}+b_{0}\right].$
(34)

If $g^{R}$ is smooth throughout on $[c,d],$ we can either directly approximate
the CFS representation with a Chebyshev series using Chebfun (cf. Trefethen et
al. 2014) or transform it into a Chebyshev series using the techniques shown
in Appendix 8. If $g^{R}$ is a piecewise continuous function containing a
singularity444We refer a singularity as a point at which is not defined, or a
point which fails to be well-behaved after differentiability, then we use the
Fourier–Padé ideas to locate the singularity in $g^{R}$ and form accurate
approximation. The details can be found in the Appendix 9.

Knowing singularities $\tilde{x}_{1}\ldots\tilde{x}_{\mathcal{K}+1}$ in
$[c,d],$ we divide $f$ and $g^{R}$ into a set of piecewise continuous
functions and then approximate them with Chebyshev series using chebfun (cf.
Trefethen et al. 2014, Chapter 1.4). Accordingly, we have a set of
$\mathcal{K}$ polynomials $f_{M}$ and $g_{N},$ each of degree at most $M$ and
$N$ on the subintervals $[\tilde{x}_{k},\tilde{x}_{k+1}],$ i.e.,

$\displaystyle
f_{M}=\sum\limits_{k=1}^{\mathcal{K}}f_{k,M}\mathds{1}_{[\tilde{x}_{k},\tilde{x}_{k+1}]},\quad
g^{R}_{N}=\sum\limits_{k=1}^{\mathcal{K}}g^{R}_{k,N}\mathds{1}_{[\tilde{x}_{k},\tilde{x}_{k+1}]}.$
(35)

Here, $\mathds{1}_{[\tilde{x}_{k},\tilde{x}_{k+1}]}$ as the indicator function
in the interval $[\tilde{x}_{k},\tilde{x}_{k+1}]$ and

$\displaystyle
f_{k,M}(\tilde{x})=\sum\limits_{m=0}^{M}\alpha_{k,m}^{cheb}T_{m}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x}),\quad
g^{R}_{k,N}(x)=\sum\limits_{n=0}^{N}\beta_{k,n}^{cheb}T_{n}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x}).$
(36)

Then, using the techniques implemented in cheb2leg (cf. Townsend et al. 2018),
we transform both $f_{N}$ and $g^{R}_{M}$ into Legendre series defined by

$\displaystyle
f_{k,M}(\tilde{x})=\sum\limits_{m=0}^{M}\alpha_{k,m}^{leg}P_{m}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x}),\quad
g^{R}_{k,N}(\tilde{x})=\sum\limits_{n=0}^{N}\beta_{k,n}^{leg}P_{n}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}.$
(37)

We use the algorithm for the convolution of Legendre series described in
Section 2 to approximate their convolution
$h(\tilde{x})=(f_{M}*g_{N})(\tilde{x})$ on the subintervals
$[\tilde{x}_{k},\tilde{x}_{k+1}].$ Finally, transforming $h(\tilde{x})$ back
into Chebyshev series using leg2cheb, $V$ can be approximated by

$\displaystyle e^{-r(T-t)}Kh(\tilde{x})$
$\displaystyle=e^{-r(T-t)}K\left(f_{M}*g_{N}\right)(\tilde{x})$ (38)
$\displaystyle=e^{-r(T-t)}K\sum\limits_{k=1}^{M_{k}}V_{N_{k}}\mathds{1}_{[\tilde{x}_{k},\tilde{x}_{k+1}]}$
(39) where, $\displaystyle V_{N_{k}}$
$\displaystyle=\sum\limits_{k=1}^{N_{k}}\gamma_{k}T_{k}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x})$
(40)

Using (38), we can generate a set of option prices with a value of $K$ and a
range of $S_{t}.$ However, in the financial markets, option price quotes
always appear with a value of $S_{t}$ and a range of $K.$ To fit in this
financial phenomenon, we modify (38) using the fact of
$K=Se^{-\tilde{x}}=e^{x-\tilde{x}}$ so that we obtain the new pricing formula
of

$\displaystyle V(x,K,t)=e^{-r(T-t)+x-\tilde{x}}h(\tilde{x}).$ (41)

###### Remark 4.1

In Chebfun, there is a built-in algorithm to detect singularities
automatically in a piecewise continuous function (cf. Pachón et al. 2010).
However, since PDFs can be easily expressed in the complex Fourier series and
then extend to the Fourier–Padé series (cf. Chan 2016, 2018), we instead use
Fourier–Padé ideas to locate singularities or approximate PDFs.

####  4.1.1 European Vanilla Call Options as illustration

We now consider pricing a European vanilla call, which can be exercised only
at its maturity, defined in (4.1) with a payoff function of

$\displaystyle U(S_{T},K)=\max\left(S_{T}-K\right).$ (42)

We first transform the payoff into

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

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

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

where $\tilde{x}=x-\log K$ and $g^{R}(\tilde{x}):=g(-\tilde{x})$ is a
reflecting function. To make the CONLeg more efficient, we define a truncated
computational interval $[c,d]$ (cf. Section 5), which satisfies condition
$(\ref{eqn:charAppx}),$ to replace $[-\infty,\infty].$ Then, $V(x,K,t)$ is
reformulated as

$\displaystyle V(x,K,t)$ $\displaystyle\approx
e^{-r(T-t)}K\int_{c}^{d}\max(e^{y}-1,0)g^{R}\left(\tilde{x}-y\right)\mathrm{d}y$
$\displaystyle=e^{-r(T-t)}K\int_{\min(d,x-c)}^{\max(c,x-d)}f(y)g^{R}(\tilde{x}-y)\mathrm{d}y,\quad\tilde{x}\in[2c,2d]$
$\displaystyle=e^{-r(T-t)}Kh(\tilde{x}),$ (45)

where $f(y):=\max(e^{y}-1,0).$ To easily digest how the CONLeg method
approximates $h(\tilde{x})$, we assume that $g^{R}$ is a piecewise smooth
function containing only one jump $y=0$ appearing in $f,$ the payoff function,
on $[c,d]$. We can use chebfun to approximate $f$ and $g^{R}$ on
$[c,0,d]$555One should note that when $y\leq 0,$ $\max(e^{y}-1,0)=0$. and
$[c,d]$ respectively. Using the techniques described in Section 4.1, the
European call option pricing formula is given by:

$\displaystyle V(x,K,t)$ $\displaystyle\approx
e^{-r(T-t)}K\left(f_{M}*g_{N}\right)(\tilde{x})$
$\displaystyle=e^{-r(T-t)}K\sum\limits_{k=1}^{4}V_{N_{k}}\mathds{1}_{[\tilde{x}_{k},\tilde{x}_{k+1}]}\quad\tilde{x}\in[2c,2d],$
(46)

where,
$V_{N_{k}}=\sum\limits_{k=1}^{N_{k}}\gamma_{k}T_{k}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x}),$
and
$\tilde{x}_{1}=2c,\tilde{x}_{2}=c,\tilde{x}_{3}=0,\tilde{x}_{4}=d,\tilde{x}_{5}=2d.$

Moreover, if we only focus on $[c,d],$ we have

$\displaystyle
V(x,K,t)=e^{-r(T-t)}K\left(V_{N_{1}}\mathds{1}_{[c,0]}+V_{N_{2}}\mathds{1}_{[0,d]}\right)\quad\tilde{x}\in[c,d],$
(47)

where $\tilde{x}_{1}=c,\tilde{x}_{2}=0,\tilde{x}_{3}=d.$ See Section 6.1 for
an example.

The CONLeg method is not limited from pricing European vanilla call option
(4.1.1); it can be readily extended into a put option or other options with
different pay-off structures, e.g., Cash-or-Nothing options. In Table 1, we
list all financial contingency claims we consider in this paper and both their
payoff functions and transformed payoff functions.

Table 1: Payoff functions and their transforms for a variety of financial contingency claims. $\mathds{1}$ represents an indicator function and $n<\infty$ is any positive integer. The singularity always exits at $y=0$ in the transformed payoff function when $y=x+\chi-\log K.$ Financial Contingency Claim | Payoff Function | Transformed Payoff Function
---|---|---
| $U(S_{T},K)$ | $U(e^{x+\chi-\log K},K)$
Call | $\max(S_{T}-K,0)$ | $K\max(e^{x+\chi-\log K}-1,0)$
Put | $\max(K-S_{T},0)$ | $K\max(1-e^{x+\chi-\log K},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}(1-e^{n(x+\chi-\log K)})\mathds{1}_{e^{x+\chi-\log K}\leq 1}$

###  4.2 Pricing Formulae for Bermuda Options

Consider now $\log S_{t}:=x_{t}$ driven by a Lévy process and a Bermudan
option with strike $K$ and maturity $T$ that can be exercised only on a given
number of exercise dates $t=t_{0}<t_{1}\leq t_{2}\leq\ldots t_{l}\leq
t_{l+1}\leq\ldots\leq t_{L}=T.$ We can write the Bermudan pricing formula for
such an option as

$\displaystyle V(x_{t_{l}},K,t_{l})$
$\displaystyle=\begin{cases}U(e^{x_{t_{l}}},K,t_{l})&l=L,\,t_{L}=T\\\
\max\left(C(x_{t_{l}},K,t_{l}),U(e^{x_{t_{l}}},K,t_{l})\right)&l=1,2,3,\dots,L-1\\\
C(x_{t_{l}},K,t_{l})&l=0\end{cases},$ (48)

where, $U(e^{x_{t_{l}}},K,t_{l})$ is the payoff function at $t_{l}.$ That is,
if the payoff function is a call, then $U(e^{x_{t_{l}}},K,t_{l})$ is
transformed into $\max\left(e^{x_{t_{l}}}-K,0\right).$ In (48),
$C(x_{t_{l}},K,t_{l})$ at each $t_{j}$ can be defined as

$\displaystyle C(x_{t_{j}},K,t_{j})$
$\displaystyle=e^{-r(t_{j+1}-t_{j})}\mathbb{E}\left(V(x_{t_{j+1}},K,t_{j+1})|x_{t_{j}}\right).$
(49)

To apply the CONLeg method to approximate $C(x_{t_{l}},K,t_{l}),$ we first
understand that

$\displaystyle
S_{t_{l}}=e^{-r(t_{l+1}-t_{l})}\mathbb{E}\left(S_{t_{l+1}}|S_{t_{l}}=e^{x_{t_{l}}}\right)=e^{-r(t_{l+1}-t_{l})}\mathbb{E}\left(e^{x_{t_{l}}+X_{t_{l+1}}-X_{t_{l}}}\right)=e^{x_{t_{l}}}$
(50)

is a martingale process. We also denote $\tilde{x}_{t_{l}}$ as $x_{t_{l}}-\log
K$ and follow Section 4.1 to approximate $C(x_{t_{l}},K,t_{l})$ as European
option prices at $t_{l}$. Then we can transform $C(x_{t_{l}},K,t_{l})$ into

$\displaystyle=e^{-r(t_{l+1}-t_{l})}\mathbb{E}\left(V(x_{t_{l+1}},K,t_{l+1})|x_{t_{l}}\right)$
$\displaystyle=e^{-r(t_{l+1}-t_{l})}\int^{+\infty}_{-\infty}V(x_{t_{l}}+\chi-\log
K,t_{l+1})g(\chi)\mathrm{d}\chi,\quad\chi\in X_{t_{l+1}}-X_{t_{l}}$
$\displaystyle=e^{-r(t_{l+1}-t_{l})}K\int_{\min(d,x-c)}^{\max(c,x-d)}f(y)g^{R}(\tilde{x}_{t_{l}}-y)\mathrm{d}y$
$\displaystyle=e^{-r(t_{l+1}-t_{l})}Kh(\tilde{x}_{t_{l}}).$ (51)

Since we approximate $h(\tilde{x}_{t_{l}})$ with the CONLeg method and use the
Chebyshev series to present the Bermuda option prices $C(x_{t_{l}},K,t_{l}),$
accordingly, we can further modify (48) with a new form of

$\displaystyle V(x_{t_{l}},K,t_{l})$
$\displaystyle=\begin{cases}K\tilde{f}(\tilde{x}_{t_{l}})&l=L,\,t_{L}=T\\\
K\max\left(e^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}}),\tilde{f}(\tilde{x}_{t_{l}})\right)&l=1,2,3,\dots,L-1\\\
Ke^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}})&l=0\end{cases},$ (52)

where $K\tilde{f}(\tilde{x}_{t_{l}}):=U(e^{x_{t_{l}}},K,t_{l}).$ Since the no-
arbitrage assumption leads to the requirement that $\partial V/\partial x$ is
continuous and $V(x_{t_{l}},K,t_{l})=U(e^{x_{t_{l}}},K,t_{l})$ at the early
exercise curve, we must determine exercise point $x^{*}_{t_{l}}$ appearing in
$V(x_{t_{l}},K,t_{l})=U(e^{x_{t_{l}}},K,t_{l}).$ One way to do this is to use
the Newton method proposed in Fang and Oosterlee (2009b) to find $x_{t_{l}}.$
However, since $V$ and $U$ are represented by piecewise smooth polynomials
(chebfun), we can apply a built-in function roots in Chebfun to efficiently
find these zeros. To do so, we first approximate
$\tilde{f}(\tilde{x}_{t_{l}})$ as a Chebyshev series, then apply the roots
function to find $\tilde{x}^{*}_{t_{l}}$ in the following equality:

$e^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}})-\tilde{f}(\tilde{x}_{t_{l}})=0.$

Once we have $\tilde{x}^{*}_{t_{l}}$ and use it as a break point, we
approximate

$\displaystyle\max\left(e^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}}),f(\tilde{x}_{t_{l}})\right)$
(53)

with two different Chebyshev series. Considering $\tilde{x}^{*}_{t_{l}}$ as a
singularity and combining other singularities,
$\tilde{x}_{t_{l},1}\ldots\tilde{x}_{t_{l},\mathcal{K}+1}\in\tilde{x}_{t_{l}},$
e.g., singularities in a non-smooth PDF and/or in a payoff function, in
$V(x_{t_{l}},K,t_{l}),$ we approximate $V(x_{t_{l}},K,t_{l}),$ with a set of
Chebyshev series given by

$\displaystyle
K\sum\limits_{k=0}^{M_{k}}V_{N_{k}}\mathds{1}_{[\tilde{x}_{t_{l},k},\tilde{x}_{t_{l},k+1}]}\hbox{
and
}V_{N_{k}}=\sum\limits_{k=1}^{N_{k}}\gamma_{k}T_{k}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x}).$
(54)

Finally, summarising the methods above, we present the pseudo-code of our
algorithm computing Bermudan option prices in Algorithm 1. A numerical example
is also presented in Section 6.2.

Result: Bermuda option price $V(x_{t},K,t)$ at time t

initialisation;

discretise $[t,T]$ into timesteps $t=t_{0},t_{1},\ldots,t_{l},\ldots,t_{L}=T$;

$t_{l}=t_{L-1}$;

while  _$t_{l}\neq t$_ do

compute $C(x_{t_{l}},K,t_{l})$ using the CONLeg method;

$C(x_{t_{l}},K,t_{l})=e^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}})$ in (52);

find $\tilde{x}^{*}_{t_{l}}$ in
$e^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}})-\tilde{f}(\tilde{x}_{t_{l}})=0$;

compute
$\max\left(e^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}}),\tilde{f}(\tilde{x}_{t_{l}})\right)$
with two Chebyshev series (54);

$V(x_{t_{l}},K,t_{l})=K\sum\limits_{k=0}^{M_{k}}V_{N_{k}}\mathds{1}_{[\tilde{x}_{t_{l},k},\tilde{x}_{t_{l},k+1}]}$;

next $t_{l}$;

end while

return $V(x_{t},K,t)$ equal to $e^{-r(t_{1}-t)}Kh(\tilde{x}_{t}),$ where
$t_{0}=t$;

Algorithm 1 Algorithm for computing Bermudan option price $V(x_{t},K,t)$ at
time $t$ based on (48).

###  4.3 Pricing Formulae for American Options

There are two basic approaches to evaluating American options based on our
method for Bermudan options. As suggested in Fang and Oosterlee (2009b), one
simple approach is to approximate an American option by a Bermudan option with
many exercise opportunities. In other words, increase the number of exercise
opportunities $L$ to a very large value. An alternative approach is to use
Richardson extrapolation on a series of Bermudan options with an increasing
number of $L$ (cf. Geske and Johnson 1984, Chang et al. 2007). We adapt the
latter approach (which is also implemented in Fang and Oosterlee (2009b)) to
price the American option here. Therefore, implementing the 4-point Richardson
extrapolation scheme (cf. Fang and Oosterlee 2009b), we have the American
option price given by

$\displaystyle V_{Amer}(L)={1\over
21}\left(64V(2^{L+3})-56V(2^{L+2})+14V(2^{L+1})-V(2^{L})\right),$ (55)

where $V_{Amer}(L)$ denotes the approximated value of the American option and
$V(\cdot)$ is the pricing formulae for Bermudan options in (52).

###  4.4 Pricing Formulae for Discretely Monitored Barrier Options

A barrier option is an early-exercise option whose payoff depends on the stock
price crossing a pre-set barrier level during the option’s lifetime. We call
the option an up-and-out, knock-out, or down-and-out option when the option’s
existence fades out after crossing the barrier level. Like European vanilla
options, these options can all be written as either put or call contracts that
have a pre-determined strike price on an expiration date. In this paper, we
only investigate two basic types of barrier options.

  1. 1.

Down-and-out barrier (DO) option: A down-and-out barrier option is an option
that can be exercised at a pre-set strike price on an expiration date as long
as the stock price that drives the option does not go below a pre-set barrier
level during the option’s lifetime. As an illustration, if the stock price
falls below the barrier, the option is “knocked-out” and immediately carries
no value.

  2. 2.

Up-and-out barrier (UO) option: Similar to a down-and-out barrier option, an
up-and-out barrier option will be knocked out when the stock price rises above
the barrier level during the option’s lifetime. Once it is knocked out, the
option cannot be exercised at a predetermined strike price on an expiration
date.

The structure of discretely monitored barrier options is the same as that of
Bermudan options. Instead of having a pre-set exercise date and an early-
exercise point like Bermudan options, barrier options have a pre-set monitored
date and a barrier level. In the case of Bermudan options, when the stock
price goes across the early exercise point, a payoff occurs, and the option
expires immediately. In the same manner, a barrier option knocks out
immediately when the barrier level is crossed. The barrier level acts exactly
the same as the exercise point in Bermudan options. However, in the case of a
barrier option without a rebate, no payoff occurs when the barrier level is
reached; otherwise, a rebate occurs when a barrier option is knocked out.

We use a rebate DO option to illustrate the CONLeg method to approximate
discretely monitored barrier option prices. Suppose that we have a rebate DO
option driven by $S_{t}$ with a barrier $B$, a rebate $R_{b},$ a strike $K$
and a series of monitoring dates $L$: $t=t_{0}<\ldots<t_{l}<\ldots<t_{L}=T;$
the option formulae can be described as

$\displaystyle V(x_{t_{l}},K,t_{l})$
$\displaystyle=\begin{cases}U(e^{x_{t_{l}}},K,t_{l})\mathds{1}_{\log
B>x_{t_{l}}}+R_{b}\mathds{1}_{\log B\leq x_{t_{l}}}&l=L,\,t_{L}=T\\\
C(x_{t_{l}},K,t_{l})\mathds{1}_{\log
B>x_{t_{l}}}+e^{-r(T-t_{l})}R_{b}\mathds{1}_{\log B\leq
x_{t_{l}}}&l=1,\dots,L-1\\\ C(x_{t_{l}},K,t_{l})&l=0\end{cases},$ (56)

where $\mathds{1}$ is an indicator function and $U(e^{x_{t_{l}}},K,t_{l})$ is
again either a call or put payoff. We follow the ideas of (4.2) and (52) in
Section 4.2 to approximate $C(x_{t_{l}},K,t_{l})$ such that

$\displaystyle
e^{-r(t_{l+1}-t_{l})}\mathbb{E}\left(V(x_{t_{l+1}},K,t_{l+1})|x_{t_{l}}\right)=e^{-r(t_{l+1}-t_{l})}Kh(\tilde{x}_{t_{l}}).$
(57)

After we apply the CONLeg method, (56) can be transformed into

$\displaystyle
V(x_{t_{l}},K,t_{l})=\begin{cases}K\left(\tilde{f}(\tilde{x}_{t_{l}})\mathds{1}_{\log(B/K)>\tilde{x}_{t_{l}}}+{R_{b}\over
K}\mathds{1}_{\log(B/K)\leq x_{t_{l}}}\right)&l=L,\,t_{L}=T\\\
K\left(e^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}})\mathds{1}_{\log({B\over
K})>\tilde{x}_{t_{l}}}+e^{-r(T-t_{l})}{R_{b}\over K}\mathds{1}_{\log({B\over
K})\leq x_{t_{l}}}\right)&l=1,\dots,L-1\\\
Ke^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}})&l=0\end{cases}.$ (58)

In (58), since there is a jump at $\log(B/K),$ the barrier, at $t_{l}$, we use
$\log(B/K)$ as a break point and approximate

$\displaystyle
e^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}})\mathds{1}_{\log({B\over
K})>\tilde{x}_{t_{l}}}+e^{-r(T-t_{l})}{R_{b}\over K}\mathds{1}_{\log({B\over
K})\leq x_{t_{l}}}$ (59)

with two Chebyshev series. Moreover, combining other singularities
$\tilde{x}_{t_{l},1}\ldots\tilde{x}_{t_{l},\mathcal{K}+1}\in\tilde{x}_{t_{l}}$
in $V(x_{t_{l}},K,t_{l})$, we can formulate $V(x_{t_{l}},K,t_{l})$ with a set
of Chebyshev series given in (54). Finally, the pseudo-code of our algorithm
calculating discretely monitored DO barrier option prices can be found in
Algorithm 2.

Result: discretely monitored barrier option price $V(x_{t},K,t)$ at time $t$

initialisation;

discretise $[t,T]$ into timesteps $t=t_{0},t_{1},\ldots,t_{l},\ldots,t_{L}=T$;

$t_{l}=t_{L-1}$;

while  _$t_{l}\neq t$_ do

compute $C(x_{t_{l}},K,t_{l})$ using the CONLeg method;

$C(x_{t_{l}},K,t_{l})=e^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}})$ in (58);

compute $\begin{cases}e^{-r(t_{l+1}-t_{l})}h(\tilde{x}_{t_{l}})&\mbox{if
}\log(B/K)>\tilde{x}_{t_{l}}\\\ e^{-r(T-t_{l})}R_{b}/K&\mbox{if
}\log(B/K)\leq\tilde{x}_{t_{l}}\end{cases}$ with two Chebyshev series (54);

$V(x_{t_{l}},K,t_{l})=K\sum\limits_{k=0}^{M_{k}}V_{N_{k}}\mathds{1}_{[\tilde{x}_{t_{l},k},\tilde{x}_{t_{l},k+1}]}$;

next $t_{l}$;

end while

return $V(x_{t},K,t)$ equal to $e^{-r(t_{1}-t)}Kh(\tilde{x}_{t}),$ where
$t_{0}=t$;

Algorithm 2 Algorithm for computing discretely monitored barrier option price
$V(x_{t},K,t)$ at time $t$ based on (56).

For the UO barrier options, we can use (58) and Algorithm 2 to compute their
prices, but we consider the condition of the option knocked out when the stock
price rises above $B,$ i.e.,

$\displaystyle
V(x_{t_{l}},K,t_{l})=\begin{cases}U(e^{x_{t_{l}}},K,t_{l})\mathds{1}_{\log
B<x_{t_{l}}}+R_{b}\mathds{1}_{\log B\geq\tilde{x}_{t_{l}}}&l=L,\,t_{L}=T\\\
C(x_{t_{l}},K,t_{l})\mathds{1}_{\log
B<x_{t_{l}}}+e^{-r(T-t_{l})}R_{b}\mathds{1}_{\log
B\geq\tilde{x}_{t_{l}}}&l=1,\dots,L-1\\\ C(x_{t_{l}},K,t_{l})&l=0\end{cases}.$
(60)

###  4.5 Hedging Formulae and Choice of Truncated Intervals

We now turn our attention to deriving the option Greek values. In particular,
we focus on deriving three option Greek values—Delta ($\Delta$), Gamma
($\Gamma$), and Vega. Delta is defined as the rate of change in the option
value with respect to changes in the underlying asset price; 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 is often based on
changes in the recent, historical prices of a trading instrument. Other Greek
values, such as Theta, can be derived in a similar fashion; however, depending
on the characteristic function, the derivation expressions might be rather
lengthy. We omit them here, as many terms are repeated.

As mentioned above, Delta is the first derivative of the value $V$ of the
option with respect to the underlying instrument price $S$. Hence,
differentiating the convolution form of $V$ in European options (38), Bermuda
options (54), American options (55) and barrier options (54) with respect to
$S,$ we have

$\displaystyle\Delta_{t}$ $\displaystyle={\partial V(x,K,t)\over\partial
S}=e^{-r(T-t)}K{\partial
h(\tilde{x})\over\partial\tilde{x}}{\partial\tilde{x}\over\partial x}{\partial
x\over\partial S},\quad\tilde{x}=x-\log K.$ (61)

Since $\partial x/\partial\tilde{x}=1$ and $\partial x/\partial S=\exp(-x),$
$\Delta_{t}$ simply becomes

$\displaystyle\quad e^{-r(T-t)-x}K{\partial
h(\tilde{x})\over\partial\tilde{x}}$
$\displaystyle=e^{-r(T-t)-x}K\sum\limits_{k=1}^{M_{k}}{\partial
V_{N_{k}}\over\partial\tilde{x}}\mathds{1}_{[\tilde{x}_{k},\tilde{x}_{k+1}]},$
(62) with $\displaystyle{\partial V_{N_{k}}\over\partial\tilde{x}}$
$\displaystyle=\sum\limits_{k=1}^{N_{k}}\gamma_{k}{\partial
T_{k}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x})\over\partial\tilde{x}}.$
(63)

To express the first derivative of the Chebyshev series in (63), we adopt the
fact of

${d\over dx}T_{n}(x)={n\over 2}{T_{n-1}(x)-T_{n+1}(x)\over 1-x^{2}}$

(cf. Mason and Handscomb 2002, (2.4.5)) and
$\mathrm{d}\psi_{[c,d]}(\tilde{x})=(2/(d-c))\mathrm{d}\tilde{x},$ such that we
have

$\displaystyle{\partial
T_{k}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x})\over\partial\tilde{x}}={2k\over
d-c}{T_{k+1}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x})-T_{k-1}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x})\over
1-\tilde{x}^{2}}.$ (64)

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

$\displaystyle\Gamma_{t}$ $\displaystyle={\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}$
$\displaystyle=e^{-r(T-t)-2x}K\left({\partial^{2}h(\tilde{x})\over\partial^{2}\tilde{x}}-{\partial
h(\tilde{x})\over\partial\tilde{x}}\right)$
$\displaystyle=e^{-r(T-t)-2x}K\Big{(}\sum\limits_{k=1}^{M_{k}}{\partial^{2}V_{N_{k}}\over\partial^{2}\tilde{x}}\mathds{1}_{[\tilde{x}_{k},\tilde{x}_{k+1}]}-\sum\limits_{k=1}^{M_{k}}{\partial
V_{N_{k}}\over\partial\tilde{x}}\mathds{1}_{[\tilde{x}_{k},\tilde{x}_{k+1}]}\Big{)}.$
(65)

To find $\partial^{2}V_{N_{k}}/\partial^{2}\tilde{x},$ we may use

${d^{2}\over dx^{2}}T_{n}(x)={n\over
4}{(n+1)T_{n-2}(x)-2nT_{n}(x)+(n-1)T_{n+2}(x)\over(1-x^{2})^{2}}$

(cf. Mason and Handscomb 2002, Problem 2.5.17), and thus find

$\displaystyle{\partial^{2}V_{N_{k}}\over\partial^{2}\tilde{x}}=\sum\limits_{k=1}^{N_{k}}\gamma_{k}{\partial^{2}T_{k}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x})\over\partial^{2}\tilde{x}}$
(66)

and

$\displaystyle{\partial^{2}T_{k}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x})\over\partial^{2}\tilde{x}}$
$\displaystyle={k\over(d-c)(1-\tilde{x}^{2})^{2}}\Bigg{(}{(k+1)T_{k-2}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x})}$
$\displaystyle\quad-{2kT_{k}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x})+(k-1)T_{k+2}\circ\psi_{[\tilde{x}_{k},\tilde{x}_{k+1}]}(\tilde{x})}\Bigg{)}.$
(67)

Likewise, we can obtain the formula for Vega, ${\partial
V\over\partial\sigma_{t}},$ where $\sigma_{t}$ is the initial value of the
volatility at time $t.$ For example, for the GBM model with $\sigma_{t}$ as
the initial value of the volatility, we derive Vega as follows:

$\displaystyle{\partial V(x,K,\sigma_{t},t)\over\partial\sigma_{t}}$
$\displaystyle=e^{-r(T-t)}K\int_{\min(d,\tilde{x}-c)}^{\max(c,\tilde{x}-d)}f(y){\partial
g^{R}(\tilde{x}-y)\over\partial\sigma_{t}}\mathrm{d}y,\quad\tilde{x}\in[2c,2d].$
(68)

After we differentiate $V$ with respect to $\sigma_{t}$ to obtain (68), we can
approximate (68) with the CONLeg method.

If the closed-formed PDF $g$ of the stochastic process does not exist, as we
mentioned before, we express $g^{R}$ with the CFS expression such that

$\displaystyle{\partial
g^{R}(\tilde{x}-y)\over\partial\sigma_{t}}=\mathfrak{Re}\left[2\sum_{k=1}^{N}{\partial
b_{k}\over\partial\sigma_{t}}e^{-i\frac{2\pi}{d-c}k(\tilde{x}-y)}\right],\hbox{
and }{\partial
b_{k}\over\partial\sigma_{t}}={\partial\varphi(-\frac{2\pi}{d-c}k,\sigma_{t})\over\partial\sigma_{t}},$
(69)

where $\varphi$ contains the parameter $\sigma_{t}.$

##  5 Choice of Truncated Intervals

In this section, we adopt the ideas of Fang and Oosterlee (2009a) and Chan
(2018) to choose the interval $[c,d]$. The choice of the interval $[c,d]$
plays the crucial role in the accuracy of the CONLeg method. If the choice of
$[c,d]$ is too big, the CONLeg method can perform inefficiently. On the
contrast, if $[c,d]$ is too small, the CONLeg method can produce inaccurate
option prices/hedging values. Accordingly, a minimum and substantial interval
$[c,d]$ can be chosen to capture most of the mass of a PDF such that our
algorithm can in turn yield the highest accuracy. 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 (30). Based on their ideas, we have the following expression for
$[c,d]$:

$\displaystyle d$
$\displaystyle=\left|c_{1}+L_{n}\sqrt{c_{2}+\sqrt{c_{4}}}\right|$
$\displaystyle c$ $\displaystyle=-d,$ (70)

where $c_{1},$ $c_{2},$ and $c_{4}$ are the first, second and fourth
cumulants, respectively, of the stochastic process and $L_{n}\in[8,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 2.

Table 2: The first, second, and fourth cumulants of various models. Lévy models | cumulants
---|---
BS |  $c_{1}=(r-q+\omega)t$ $c_{2}=\sigma^{2}t,$ $c_{4}=0,$ $\omega=-0.5\sigma^{2}$
NIG | $c_{1}=(r-q+\omega)t+\delta t\beta/\sqrt{\alpha^{2}-\beta^{2}}$
| $c_{2}=\delta t\alpha^{2}(\alpha^{2}-\beta^{2})^{-3/2}$
| $c_{4}=\delta
t\alpha^{2}(\alpha^{2}+4\beta^{2})^{-3/2}(\alpha^{2}-\beta^{2})^{-7/2}$
|
$\omega=-0.5\sigma^{2}-\delta(\sqrt{\alpha^{2}-\beta^{2}}-\sqrt{\alpha^{2}-(\beta+1)^{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/\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)$

In general, the truncated intervals in (5) work for smooth/non-smooth PDFs
with/without singularities. However, as the cumulants in each process in Table
2 contain $t$ , we notice that if $t$ is too small, the cumulants shrink and
make $[c,d]$ too small for the CONLeg to produce accurate option
prices/hedging values. To solve the problem heuristically, as long as $t$ is
less than or equal to $0.2,$ we should add extra 0.5 to $d$ in (5) in order to
produce a substantial interval $[c,d]$ for the CONLeg method.

##  6 Numerical Results

The main purpose of this section is to test the accuracy and efficiency of the
CONLeg method through various numerical tests. This involves the ability of
the method to price any options that are deep in/out of the money and have
long/short maturities. Most importantly, we show that the algorithm exhibits
good accuracy even when the PDF is smooth/non-smooth. A number of popular
numerical methods are implemented to compare the algorithm in terms of the
error convergence 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. 2015), the CONV method (an FFT method, Lord et al. 2008), the
Lewis-FRFT (a fractional FFT method, Lewis 2001, Chourdakis 2004), the QUAD-
CONV (a combination of the quadrature and CON methods; see O’Sullivan 2005,
Chen et al. 2014), and the SWIFT methods (a wavelet-based method; see Ortiz-
Gracia and Oosterlee 2013, Maree 2015, Ortiz-Gracia and Oosterlee 2016, Maree
et al. 2017). 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. (2015) report that this filter provides
better algebraic convergence than the other options. We also set the damping
factors of the CONV to 0 for pricing European options. 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. Finally, the code is written in MATLAB and also the codes 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).

Since we use the built-in Chebfun (Trefethen et al. 2014) commands, mainly
chebfun, conv, diff, roots and simplify, to price and hedge the options in
this paper, Chebfun (http://www.chebfun.org/download/) is required for our
option pricing algorithm. For a demonstration, we give the MATLAB code for
computing European options in Appendix 10. As Chebfun makes use of adaptive
procedures that aim to find the right number of points automatically so as to
represent each function to roughly machine precision, that is, about 15 digits
of relative accuracy, we allow Chebfun to approximate $f$ and $g^{R}$
automatically. As Chebfun is implemented in MATLAB, we should point out that
since Chebfun calls functions, nested functions and sub-functions
considerably, the overheads occur and accordingly, this slows the CONLeg
method’s computation speed down. Hence, as an introductory of the CONLeg
method, we weigh more on the convenience of the method over its efficiency at
this stage of development.

In all numerical experiments, the range of option prices we measure is based
on the input range of $S$ or $K.$ To allow our method to have accurate
approximation of deep in/out-of-money option prices, the range of $S$ or $K$
lie in the intervals $[K-20,K+20]$ or $[S-20,S+20]$ respectively. Moreover, as
we look at the intervals of $[K-20,K+20]$ or $[S-20,S+20],$ we set $[c,d]$
rather than $[2c,2d]$ for approximate option prices and their option Greek
values 666To achieve this, we can set a flag–“same”–in conv. See Appendix 10
for details..

We use $N$ to define the number of terms/grid points of the methods we compare
against the CONLeg method. When measuring approximation errors of the
numerical methods, we use absolute errors, the $L_{2}$ norm errors $R_{2}$ and
the infinity norm errors $R_{\infty}$ as the measurement units. Moreover, to
improve the accuracy of our method in pricing/hedging European type options,
we use the well-known call-put parity relationship,

$\displaystyle V_{call}(x,K,t)=V_{put}(x,K,t)+Se^{-q(T-t)}-Ke^{-r(T-t)},$ (71)

to approximate call prices once we have put prices ready.

###  6.1 European Type Options

We consider three different test cases based on the following PDFs and other
parameters:

$\displaystyle\textbf{GBM1}:\,S$
$\displaystyle=100,\,K=80-120,\,\sigma=0.15,\,T=1.0,\,r=0.03,\,q=0.01.$ (72)
$\displaystyle\textbf{GBM2}:\,S$
$\displaystyle=100,\,K=80-120,\,\sigma=0.25,\,T=50\hbox{ or
}100,\,r=0.1,\,q=0.$ (73) $\displaystyle\textbf{VG1}:\,S$
$\displaystyle=100,K=80-90,\,\sigma=0.12,\,\theta=-0.14,\,\nu=0.2,\,T=0.1,\,r=0.1,\,q=0.$
(74)

In all three numerical tests, the reference values for the GBM process and the
VG process are generated via MATLAB Financial Toolbox™–blsprice, blsdelta and
blsgamma–and the Singularity Fourier–Padé (SFP) method (cf. Chan 2018)
respectively. We also set $L_{n}=10$ in (5) when we create a computational
interval $[c,d]$ for pricing/hedging European options under the two processes.

In the first numerical test (GBM1)–Table 3, we first check for convergence
behaviour against a range of strikes $K$ from $80$ to $120$ for deep in/out-
of-the money and at-the-money vanilla put options. Apart from $q=0.01,$ the
parameters are retrieved from von Sydow et al. (2015). We declare 1000
different option prices within the range of either $K.$ In this test, even
without applying the put-call parity (71), the CONLeg can achieve very high
accuracy (around $R_{\infty}=10^{-14}$) when it is applied to approximate
option prices and its Delta $\Delta$ and Gamma $\Gamma.$ Moreover, since our
method aims to model option price/Greek curves rather than their values, our
method consumes only less than 0.1 seconds to formulate the curves in the
test. Using around 0.1 seconds to produce option price and Greek curves is a
quite reasonable computational cost for a method to meet the financial
standards. The second numerical experiment (GBM2) is devoted to the
performance of the CONLeg method for long maturity call options, which are
often encountered in the insurance and pension industry. The parameters are
retrieved from Ortiz-Gracia and Oosterlee (2016) for the test. Table 4 refers
to the second test (GBM2) and replicates Table 3 in Ortiz-Gracia and Oosterlee
(2016). In this test, with the help of (71), the CONLeg method impressively
provides high accuracy when we declare 1000 different option prices within the
range of $K$ from 80 to 120. The last set of parameters (VG1) is chosen in the
last numerical test (Table 5) because relatively slow convergence was reported
for the CONV method for very short maturities in Lord et al. (2008). This is
attributed to the PDF of the process being sharp-peaked with a difficult
logarithmic singularity $x_{sing}$ (see Fig. 2). Before we approximate the
closed-form VG PDF defined in (27) with Chebyshev series, we apply the
Fourier–Padé method (cf. Appendix 9) to locate the singularity in the PDF. We
use $x_{sing}$ as a breakpoint and approximate the PDF in two different
regions $[c,x_{sing}]$ and $[x_{sing},d]$ to resolve Gibbs phenomena. Since we
use the Fourier–Padé method to approximate the PDF, the number of Fourier
terms in the Fouier–Padé approximation is set to be 1024. In Table 5, again
with the aid of (71), within similar CPU time and 30 options measured within
the range of 80 and 90, the CONLeg method yields almost same accuracy to the
COS and filter-COS methods, but relatively higher accuracy than the Lewis-
FRFT, QUAD-CONV and CONV methods.

###### Remark 6.1

Fang and Oosterlee (2009a) suggest that the VG process with VG1 gives rise to
a probability density function that is not in $C^{\infty}(\mathds{R})$, and
thus, option pricing under VG with these parameter sets exhibits only an
algebraic convergence. Nevertheless, Chan (2018) has recently proposed the SFP
method to circumvent the problem and to approximate the VG PDF with the input
parameters of VG1. Through this method, we can regain global spectral
convergence away from singularities. For more details, we refer the readers to
Chan (2018).

Table 3: Measuring the CONLeg method in error convergence and CPU time for pricing European at/around-the-money put prices, Delta $\Delta_{t}$ and Gamma $\Gamma_{t}$ under the BSM model with parameters taken from GBM1. 1000 put prices are computed in a range of $K$ from 80 to 120. Price |  Delta $\Delta_{t}$ |  Gamma $\Gamma_{t}$ |
---|---|---|---
$R_{\infty}$ | $R_{2}$ | $R_{\infty}$ | $R_{2}$ | $R_{\infty}$ | $R_{2}$ | Time
5.329e-14 | 3.921e-13 | 5.645e-14 | 1.846e-13 | 8.538e-13 | 1.396e-12 | 9.11e-02
Table 4: Measuring the CONLeg method in error convergence and CPU time for pricing European at/around-the-money call prices, Delta $\Delta_{t}$ and Gamma $\Gamma_{t}$ under the BSM model with parameters taken from GBM2. 1000 call prices are computed in a range of $K$ from 80 to 120. | Price |  Delta $\Delta_{t}$ |  Gamma $\Gamma_{t}$ |
---|---|---|---|---
| $R_{\infty}$ | $R_{2}$ | $R_{\infty}$ | $R_{2}$ | $R_{\infty}$ | $R_{2}$ | Time
$T=50$ | 2.842e-14 | 2.821e-13 | 2.220e-16 | 1.005e-15 | 7.170e-17 | 5.007e-16 | 8.34e-02
| $R_{\infty}$ | $R_{2}$ | $R_{\infty}$ | $R_{2}$ | $R_{\infty}$ | $R_{2}$ | Time
$T=100$ | 2.842e-14 | 2.874e-13 | 3.841e-16 | 2.123e-16 | 6.126e-19 | 3.776e-18 | 9.12e-02
Table 5: Comparison of the Lewis-FRFT, CONV, QUAD–CONV, COS, filter-COS and CONLeg methods for pricing vanilla call options under the VG model with parameters taken from VG1. 30 call prices are computed in a range of $K$ from 80 to 90. | Lewis-FRFT | CONV | QUAD–CONV
---|---|---|---
$N$ | $R_{\infty}$ | $R_{2}$ | Time | $R_{\infty}$ | $R_{2}$ | Time | $R_{\infty}$ | $R_{2}$ | Time
1024 | 9.921e-03 | 1.121e-02 | 0.092 | 1.217e-04 | 8.817e-03 | 0.110 | 2.411e-04 | 7.827e-03 | 0.121
| COS | filter-COS | CONLeg
$N$ | $R_{\infty}$ | $R_{2}$ | Time | $R_{\infty}$ | $R_{2}$ | Time | $R_{\infty}$ | $R_{2}$ | Time
1024 | 2.271e-06 | 7.579e-05 | 0.089 | 6.596e-07 | 3.596e-06 | 0.0912 | 5.596e-07 | 4.511e-06 | 0.102

###  6.2 American, Barrier and Bermuda Options

In this section we again focus on three test cases based on the following PDFs
and other parameters:

$\displaystyle\textbf{NIG1}:\,S$
$\displaystyle=100,\,K=80-120,\,\alpha=15,\,\beta=-5,\,\delta=0.5,\,T=1,\,r=0.05,\,q=0.02.$
(75) $\displaystyle\textbf{CGMY1}:\,S$
$\displaystyle=0.5-1.5,\,K=1,\,C=1,\,G=5,\,M=5,\,Y=0.5,\,T=1,\,r=0.1,\,q=0.0.$
(76) $\displaystyle\textbf{CGMY2}:\,S$
$\displaystyle=90-100,\,K=100,\,C=4,\,G=50,\,M=60,\,Y=0.7,\,T=1,\,r=0.05,\,q=0.02.$
(77)

In this section, as suggested by Fang and Oosterlee (2009b), we first set
$L_{n}=8$ in (5) to create a computational interval $[c,d]$ for approximating
the PDFs of the NIG and CGMY processes. In our first test we price Bermudan
put options with 10 exercise dates with parameters of NIG1. In this test
(Table 6), a total of 170 option prices are generated in the COS and CONV
methods from a range of $K$ between 80 to 120. The CPU times are reported in
seconds, and all reference values are obtained by the CONV method with
$N=2^{20}.$ Fang and Oosterlee (2009b) report that the NIG PDF (26) is more
peaked at the mean with the parameters of NIG1. We then use the mean as a
breakpoint to approximate the PDF in two regions $[c,x_{sing}]$ and
$[x_{sing},d].$ From the result of Table 6, the accuracy of the COS method
stays the same when $N$ increases from 256 to 512. Within the similar CPU
time, we can see that the CONLeg has better accuracy than the COS method.

In the next two tests, the CGMY processes with parameters of CGMY1 and CGMY2
do not have a closed-form PDF and have singularities. Accordingly, again using
the techniques we mentioned above, we use the Fourier–Padé method to
approximate the CGMY PDFs and locate their singularities. The number of
Fourier terms in the Fouier–Padé approximation is set to be 1024 in the tests.

The prices of American options can be obtained using (55), a 4-point
Richardson extrapolation method on the prices of a few Bermudan options with
small $L.$ We compare the CONV, COS and CONLeg methods for pricing American
option price in Table 7. All reference values are obtained by the Fourier time
stepping method (Jackson et al. 2008) with the number of the Fourier terms
equal to $2^{16}.$ Within the similar CPU time, the CONLeg method achieves the
same accuracy with the COS method and is marginally better than the CONV
method.

For the last two tables, all reference values are obtained by the CONV method
with $N=2^{20}.$ In Table 8, we consider monthly monitored $(L=12)$ up-and-out
call and put options, (UO Call) and (UO Put) under the CGMY process with the
parameters of CGMY2. The barrier level $H$ is set to be 120 for the up-and-out
options. As we can see from the table, the CONLeg method is very comparable to
the SWIFT method in terms of accuracy. Finally, in Table 9, we focus on the
NIG process with parameters of NIG1, except S=90–110 and K=100. Considering
monthly monitored $(L=252)$ down-and-out call (DO Call) and put (DO Put)
options, we set $H$ equal to 80 for the options. Since $L=252$ is very large,
it causes the time interval difference between $t_{l+1}$ and $t_{l}$ to be
very small. Accordingly, the NIG PDF is very peaked and the COS method
requires a larger number of $N$ equal to 8192 to compensate a better
convergence in Table 9. Comparing the CONLeg and COS methods in the similar
CPU time, both methods can achieve similar accuracy.

Table 6: Comparison of the COS and CONLeg methods for pricing a Bermuda ($L=10$) put option under the NIG model with parameters taken from NIG1. 170 option prices are generated for both methods in a range of $K$ from 80 to 120 and $S$ equal to 100. | COS | CONLeg
---|---|---
$N$ | $R_{\infty}$ | $R_{2}$ | Time | $R_{\infty}$ | $R_{2}$ | Time
256 | 2.411e-07 | 5.411e-06 | 2.89 | 1.812e-13 | 8.576e-12 | 2.93
$N$ | $R_{\infty}$ | $R_{2}$ | Time | $R_{\infty}$ | $R_{2}$ | Time
512 | 2.411e-07 | 5.411e-06 | 3.323 | 1.812e-13 | 8.576e-12 | 2.93
Table 7: Comparison of the CONV, COS and CONLeg methods for pricing an American put option under the CGMY model with parameters taken from CGMY1. 58 and 88 option prices are computed for the CONV method and the COS method respectively in a range of $S$ from 0.5 to 1.5 and $K$ equal to 1.  $L$ in Eq. (55) |  | CONV | CONLeg
---|---|---|---
$N$ | $R_{\infty}$ | $R_{2}$ | Time (sec.) | $R_{\infty}$ | $R_{2}$ | Time (sec.)
2 | 32768 | 1.052e-04 | 8.051e-03 | 19.291 | 7.123e-05 | 6.113e-04 | 19.311
$L$ in Eq. (55) |  | COS | CONLeg
$N$ | $R_{\infty}$ | $R_{2}$ | Time (sec.) | $R_{\infty}$ | $R_{2}$ | Time (sec.)
2 | 32768 | 9.234e-05 | 4.012e-04 | 19.298 | 7.123e-05 | 6.113e-04 | 19.311
Table 8: Comparison of the SWIFT and CONLeg methods for pricing daily-monitored ($L=12$) UO call and UO put under the CGMY model with parameters taken from CGMY2. 725 option prices are computed in a range of $S$ from 90 to 110 and $K$ equal to 100. The barrier level $H$ is equal to 120. | SWIFT | CONLeg
---|---|---
| $scale$ | $R_{\infty}$ | $R_{2}$ | Time (sec.) | $R_{\infty}$ | $R_{2}$ | Time (sec.)
UO Call | 6 | 3.181e-10 | 7.641e-09 | 7.100 | 5.186e-09 | 6.886e-08 | 7.201
UO Put | 6 | 1.901e-11 | 8.621e-10 | 6.901 | 4.330e-10 | 6.131e-09 | 6.891
Table 9: Comparison of the COS and CONLeg methods for pricing daily-monitored ($L=252$) DO Call and DO Put under the NIG model with parameters taken from NIG1. 38 option prices are computed in the range of $S$ from 90 to 110 and $K$ equal to 100. The barrier level $H$ is equal to 80. | COS | CONLeg
---|---|---
| $N$ | $R_{\infty}$ | $R_{2}$ | Time (sec.) | $R_{\infty}$ | $R_{2}$ | Time (sec.)
DO Call | 8192 | 6.701e-08 | 4.641e-07 | 94.231 | 4.167e-09 | 8.651e-08 | 93.812
DO Put | 8192 | 7.232e-09 | 3.231e-08 | 92.532 | 1.619e-09 | 6.886e-08 | 92.634

##  7 Conclusions

In this paper, we have proposed an algorithm for hedging and pricing various
options based on approximating probability density functions by polynomials
(in particular, Legendre series) and computing the prices by performing
suitable convolutions with the given payoff function. We call this method
CONLeg. The main advantages of the CONLeg method are its ability to return the
price and Greeks as a function defined on a prescribed interval rather than
just point values, its ability to approximate different types options under a
process with/without a closed-from PDF, and the simplicity and convenience of
its implementation using Chebfun. We presented a proof of concept
implementation written in Chebfun, which demonstrates positive results in the
accuracy of the method. Since Chebfun uses adaptive approximations, it is
difficult at this time to discuss rates of convergence of the CONLeg method or
to compare it’s efficiency with respect to other methods in the literature
(which work with a fixed discretisation size), but initial results are
promising. For a fair comparison, one would need to implement a non-adaptive
version of CONLeg, which would be time-consuming, and not typically desirable
in a practical setting.

Our ultimate goal is to extend the method to price options with path-dependant
features under the (time-changed) Lévy process or stochastic volatility with
and without singularities. Research in this direction is already underway and
will be presented in a forthcoming manuscript.

##  8 A Closed-form Transformation of Complex Fourier Series into Chebyshev
Series

To transform a complex Fourier series (CFS) into a Chebyshev series, we apply
the result of Townsend (2014, Lemma A.3), i.e.

$\displaystyle\int_{-1}^{1}{\exp(ixy\pi)T_{q}(x)\over\sqrt{1-x^{2}}}\mathrm{d}x=\pi
i^{q}J_{q}(y\pi).$ (78)

Here, $J_{q}$ is the Bessel functions of the first kind with parameter $q$ and
$T_{q}(x)$ is the Chebyshev polynomial of degree $q.$ We first note that since
$g,$ a PDF, is a real function, the CFS representation of $g$ with a form of

$\displaystyle g(x)$ $\displaystyle\approx
g_{N}(x)=\mathfrak{Re}\left[2\sum_{k=1}^{N}\varphi\left(-\frac{2\pi}{d-c}k\right)e^{i\frac{2\pi}{d-c}kx}+\varphi(0)\right].$
(79)

can be interchangable into

$\displaystyle\mathfrak{Re}\left[\sum_{k=-N}^{N}\varphi\left(-\frac{2\pi}{d-c}k\right)e^{i{2\pi\over
d-c}kx}\right].$ (80)

Using (78) and (80), we may define a Chebyshev series on the same interval
$[c,d]$ such that

$\mathfrak{Re}\left[\sum_{k=-N}^{N}\varphi\left(-\frac{2\pi}{d-c}k\right)e^{i\frac{2\pi}{d-c}kx}\right]=\sum_{n=0}^{N}\alpha_{n}^{cheb}T_{n}\circ\psi_{[c,d]}(x),$

where,

$\displaystyle\alpha^{cheb}_{0}$ $\displaystyle=$
$\displaystyle\mathfrak{Re}\left[\sum_{\begin{subarray}{c}k=-N\\\
k\not=0\end{subarray}}^{N}{\Gamma\left({1\over
2}\right)\over\sqrt{\pi}}\varphi\left(-\frac{2\pi}{d-c}k\right)e^{i\frac{(d+c)}{d-c}k}J_{0}\left(-k\pi\right)+\varphi(0)\right],\,n=0,$
$\displaystyle\alpha^{cheb}_{n}$ $\displaystyle=$
$\displaystyle\mathfrak{Re}\left[2\sum_{\begin{subarray}{c}k=-N\\\
k\not=0\end{subarray}}^{N}\varphi\left(-\frac{2\pi}{d-c}k\right)e^{i\frac{(d+c)}{d-c}k}i^{n}J_{n}\left(-k\pi\right)\right],\,n>0.$

##  9 Locating Singularities in Probability Density Functions

Many PDFs of interest are not smooth but piecewise smooth. For example, see
Figure 2. If the locations of all singularities are not known in advance, we
can use Fourier–Padé ideas (cf. Driscoll and Fornberg 2011, Chan 2018) to
estimate the locations of singularities well enough to allow good
reconstruction nearly everywhere in the interval $[c,d]$.

The Fourier–Padé algorithm proposed in this paper is very simple to implement.
If we consider a function $g$ with a power series representation such that

$g(x)=\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},$ (82)

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

Here, $g$ 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 (83). Then, we consider
the coefficients for $x^{N+1},\ldots,x^{M+N},$ and we can yield a
Toeplitz*777A 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.$ (84)

Once $\\{q_{m}\\}_{m=0}^{M}$ is known, $\\{p_{n}\\}_{n=0}^{N}$ is found
through the terms of order N and less in (83). This yields
$\underline{p}=B\underline{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}.$ (85)

Now, assuming $g$ is a PDF, to find singularities in $g$ and to express $g$ in
a Fourier–Padé series, we first express $g$ with 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}{d-c}kx}+\varphi\left(0\right)\right].$
(86)

Then, we can differentiate (86) with respect to $x$ to obtain

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

Finally, we let $z=\exp\left({i\frac{2\pi}{d-c}x}\right)$ in the two equations
above, and they are ready for the Fourier–Padé approximation. In general, when
the PDF has a singularity, the sharp-peaked singularity point will have an
enormously large value after differentiation. In other words, Figure 2 is a
graphical illustration of the outlooks of the PDF (left) and the first
derivative (right) of the VG model after the Fourier–Padé approximation. In
the figure, we can see that the non-smooth PDF with a jump can produces a
value of $10\times 10^{11}$ at the singularity point after the first
derivative. See padeapprox in Chebfun.

Figure 2: Density functions (left) and its first derivative (right) of the VG
model with parameters are taken from VG1.

##  10 MATLAB Code

⬇

% COMPUTE GBM EUROPEAN OPTION PRICE AND GREEKS USING CONLEG METHOD.

% Tat Lung (Ron) Chan & Nick Hale, Janauary 2019

K = 100; % Strike price

window = [80, 120]; % Current stock price window

x = log(linspace(window(1), window(2), 5)/K); % Example stock prices

r = 0.035; t = 0.5; q = 0; sigma = 0.2; % Input parameters

call = true; % (call = false ==> put)

mu = (r-q-0.5*sigma^2)*t;

L = 10;

c1 = mu; c2 = sigma^2*t; c4 = 0; % Cumulants (see section 5)

d = abs(c1+L*sqrt(c2 + sqrt(c4))); % x-variable domain (see section 5)

c = -d; dom = [c, 0, d];

% Payoff function:

if ( call )

payoff = @(x) K*max(exp(x)-1, 0); % Call payoff

else

payoff = @(x) K*max(1-exp(x), 0); % Put payoff

end

payoff = chebfun(payoff, dom); % Chebfun representation:

% Guassian Density Function:

mypdf = @(u) 1./sqrt(2*pi*sigma^2*t)*exp(-0.5*(-u-mu)^2./(sigma^2*t));

mypdf = chebfun(mypdf, dom); % Chebfun representation:

% Continuous stock price:

stock = chebfun(@(x) K*exp(x), dom);

% Compute European option price via convolution:

price = exp(-r*t)*conv(payoff, mypdf, 'same');

price = simplify(price);

% Compute Greeks:

dP = diff(price); % 1st derivative

delta = dP./stock;

dP2T = diff(price, 2); % 2nd derivative

gamma = (dP2T-dP)./(stock.^2);

%Display figures:

subplot(1,3,1), plot(stock, payoff, stock, price, 'LineWidth', 2);
title('Price')

legend('Payoff', 'Price', 'Location', 'NW'), xlabel('S'),

axis([window, 0, price(x(end))])

subplot(1,3,2), plot(stock, delta, 'LineWidth', 2); title('Delta');

xlabel('S'), xlim(window)

subplot(1,3,3), plot(stock, gamma, 'LineWidth', 2); title('Gamma')

xlabel('S'), xlim(window), shg

% Display table:

Results = [K+0*x ; stock(x) ; price(x) ; delta(x) ; gamma(x)];

disp(table(Results, 'RowNames', {'strike', 'stock', 'price', 'delta',
'gamma'}))

###  10.1 Output

Figure 3: Figure produce by code above.

⬇

Results

____________________________________________________________

strike 100 100 100 100 100

stock 80 90 100 110 120

price 0.4069 2.1637 6.4983 13.5068 22.2810

delta 0.0833 0.2910 0.5771 0.8074 0.9311

gamma 0.0135 0.0269 0.0277 0.0176 0.0078

## References

  * Andricopoulos et al. (2003) Andricopoulos AD, Widdicks M, Duck PW, Newton DP (2003) Universal option valuation using quadrature methods. _J. Financ. Econ._ 67(3):447–471.
  * Andricopoulos et al. (2007) Andricopoulos AD, Widdicks M, Newton DP, Duck PW (2007) Extending quadrature methods to value multi-asset and complex path dependent options. _J. Financ. Econ._ 83(2):471–499.
  * Bondarenko (2003) Bondarenko O (2003) Estimation of risk-neutral densities using positive convolution approximation. _J. Econom._ 116(1):85–112.
  * Broadie and Yamamoto (2003) Broadie M, Yamamoto Y (2003) Application of the fast Gauss transform to option pricing. _Manage. Sci._ 49(8):1071–1088.
  * Broadie and Yamamoto (2005) Broadie M, Yamamoto Y (2005) A double-exponential fast Gauss transform algorithm for pricing discrete path-dependent options. _Oper. Res._ 53(5):764–779.
  * Cai and Kou (2011) Cai N, Kou SG (2011) Option pricing under a mixed-exponential jump diffusion model. _Manage. Sci._ 57(11):2067–2081.
  * Carr and Madan (1999) Carr P, Madan D (1999) Option valuation using the fast Fourier transform. _J. Comput. Finance_ 4(2):61–73.
  * Chan (2016) Chan TLR (2016) 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.
  * Chan (2018) Chan TLR (2018) Singular Fourier-Páde series expansion of European option prices. _Quant. Finance_ 18(7):1149–1171.
  * Chang et al. (2007) Chang CC, Chung SL, Stapleton RC (2007) Richardson extrapolation techniques for the pricing of American-style options. _J. Futures Markets_ 27(8):791–817, ISSN 1096-9934.
  * Chen et al. (2014) Chen D, Härkönen HJ, Newton DP (2014) Advancing the universality of quadrature methods to any underlying process for option pricing. _J. Financ. Econ._ 114(3):600–612.
  * Chourdakis (2004) Chourdakis K (2004) Option pricing using the fractional FFT. _J. Comput. Finance_ 8(2):1–18.
  * Colldeforns-Papiol et al. (2017) Colldeforns-Papiol G, Ortiz-Gracia L, Oosterlee C (2017) Two-dimensional Shannon wavelet inverse Fourier technique for pricing European options. _Appl. Numer. Math._ 117:115–138, ISSN 0168-9274.
  * Cont and Tankov (2004) Cont R, Tankov P (2004) _Financial Modelling with Jump Processes_. Chapman & Hall/CRC Financial Mathematics Series (Boca Raton, FL: Chapman & Hall/CRC).
  * Driscoll and Fornberg (2011) Driscoll TA, Fornberg B (2011) _The Gibbs Phenomenon in Various Representations and Applications_ (Potsdam, NY: Sampling Publishing).
  * Fang and Oosterlee (2009a) Fang F, Oosterlee CW (2009a) A novel pricing method for European options based on Fourier-cosine series expansions. _SIAM J. Sci. Comput._ 31(2):826–848.
  * Fang and Oosterlee (2009b) Fang F, Oosterlee CW (2009b) Pricing early-exercise and discrete barrier options by Fourier-cosine series expansions. _Numer. Math._ 114(1):27–62.
  * Fang and Oosterlee (2011) Fang F, Oosterlee CW (2011) A Fourier-based valuation method for Bermudan and barrier options under Heston’s model. _SIAM J. Financ. Math._ 2(1):439–463.
  * Feng and Linetsky (2008) Feng L, Linetsky V (2008) Pricing discretely monitored barrier options and defaultable bonds in Lévy process models: A fast Hilbert transform approach. _Math. Financ._ 18(3):337–384.
  * Gaß et al. (2018) Gaß M, Glau K, Mahlstedt M, Mair M (2018) Chebyshev interpolation for parametric option pricing. _Finance Stoch._ 22(3):701–731.
  * Geske and Johnson (1984) Geske R, Johnson HE (1984) The American put option valued analytically. _J. Finance_ 39(5):1511–1524.
  * Hale and Townsend (2014a) Hale N, Townsend A (2014a) An algorithm for the convolution of Legendre series. _SIAM J. Sci. Comput._ 36(3):A1207–A1220.
  * Hale and Townsend (2014b) Hale N, Townsend A (2014b) A fast, simple, and stable Chebyshev–Legendre transform using an asymptotic formula. _SIAM J. Sci. Comput._ 36(1):A148–A167.
  * Hogg et al. (2004) Hogg RV, McKean JW, Craig AT (2004) _Introduction to Mathematical Statistics_ (New Jersey: Prentice Hall), 6 edition.
  * Jackson et al. (2008) Jackson KR, Jaimungal S, Surkov V (2008) Fourier space time-stepping for option pricing with Lévy models. _J. Comput. Finance_ 12(2):1–29.
  * Leentvaar and Oosterlee (2008) Leentvaar C, Oosterlee C (2008) Multi-asset option pricing using a parallel Fourier-based technique. _J. Comput. Finance_ 12(1):1–26.
  * Lewis (2001) Lewis A (2001) A simple option formula for general jump-diffusion and other exponential Lévy processes. Available online at: http://optioncity.net/pubs/ExpLevy.pdf.
  * Lipton (2002) Lipton A (2002) The vol smile problem. _Risk_ 15:61–66.
  * Liu et al. (2016) Liu X, Xiao H, Chen R (2016) Convolutional autoregressive models for functional time series. _J. Econom._ 194(2):263–282.
  * Lord et al. (2008) Lord R, Fang F, Bervoets F, Oosterlee CW (2008) A fast and accurate FFT-based method for pricing early-exercise options under lévy processes. _SIAM J. Sci. Comput._ 30(4):1678–1705.
  * Maree (2015) Maree SC (2015) _Numerical pricing of Bermudan options using Shannon wavelet expansions_. Master’s thesis, Delft Institute of Applied Mathematics, Delft University of Technology, Delft, The Netherlands.
  * Maree et al. (2017) Maree SC, Ortiz-Gracia L, Oosterlee CW (2017) Pricing early-exercise and discrete barrier options by Shannon wavelet expansions. _Numer. Math._ 136(4):1035–1070, ISSN 0029-599X.
  * Mason and Handscomb (2002) Mason JC, Handscomb D (2002) _Chebyshev Polynomials_ (Florida, FL, USA: CRC Press).
  * Ortiz-Gracia and Oosterlee (2013) Ortiz-Gracia L, Oosterlee CW (2013) Robust pricing of European options with wavelets and the characteristic function. _SIAM J. Sci. Comput._ 35(5):B1055–B1084.
  * Ortiz-Gracia and Oosterlee (2016) Ortiz-Gracia L, Oosterlee CW (2016) A highly efficient Shannon wavelet inverse Fourier technique for pricing European options. _SIAM J. Sci. Comput._ 38(1):B118–B143.
  * O’Sullivan (2005) O’Sullivan C (2005) Path dependant option pricing under Lévy processes. _EFA 2005 Moscow Meetings Paper_ , 24.
  * Pachón (2018) Pachón R (2018) Numerical pricing of European options with arbitrary payoffs. _Int. J. Financ. Eng._ 5(2):1850015.
  * Pachón et al. (2010) Pachón R, Platte RB, Trefethen LN (2010) Piecewise-smooth chebfuns. _IMA J. Numer. Anal._ 30(4):898–916.
  * Ruijter et al. (2015) Ruijter M, Versteegh M, Oosterlee C (2015) On the application of spectral filters in a Fourier option pricing technique. _J. Comput. Finance_ 19(1):75–106.
  * Schoutens (2003) Schoutens W (2003) _Lévy Processes in Finance: Pricing Financial Derivatives_. Wiley Series in Probability and Mathematical Statistics (Chichester, UK: Wiley).
  * Su et al. (2017) Su H, Chen D, Newton DP (2017) Option pricing via QUAD: From Black-Scholes-Merton to Heston with jumps. _J. Deriv._ 24(3):9–27.
  * Townsend (2014) Townsend A (2014) _Computing with functions in two dimensions_. Thesis, University of Oxford, Oxford, UK, see http://math.mit.edu/~ajt/papers/thesis.pdf.
  * Townsend et al. (2018) Townsend A, Webb M, Olver S (2018) Fast algorithms for Toeplitz and Hankel matrices. _Mathematics of Computation_ 87(1):1913–1934.
  * Trefethen et al. (2014) Trefethen LN, Driscoll TA, Hale N (2014) _Chebfun Guide_ (Oxford: Pafnuty Publications), see http://www.chebfun.org/.
  * von Sydow et al. (2015) von Sydow L, Höök LJ, 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 CW, Ruijter MJ, Toropov A, Zhao Y (2015) BENCHOP–the BENCHmarking project in option pricing. _Int. J. Comput. Math._ 92(12):2361–2379.
  * Wong and Guan (2011) Wong HY, Guan P (2011) An FFT-network for Lévy option pricing. _J. Bank. Finance_ 35(4):988–999.
  * Zeng and Kwok (2014) Zeng P, Kwok YK (2014) Pricing barrier and Bermudan style options under time-changed Lévy processes: Fast Hilbert transform approach. _SIAM J. Sci. Comput._ 36(3):B450–B485.
  * Zhang and Oosterlee (2013) Zhang B, Oosterlee CW (2013) Efficient pricing of European-style Asian options under exponential Lévy processes based on Fourier cosine expansions. _SIAM J. Financ. Math._ 4(1):399–426.