Singular Fourier-Pad\'e Series Expansion of European Option Prices

TL;DR

This paper introduces the singular Fourier-Padé (SFP) method to efficiently and accurately price European options under Lévy and affine processes, overcoming Gibbs phenomenon issues in traditional Fourier techniques.

Contribution

The paper applies the SFP method to option pricing, providing a new approach that reduces errors caused by non-smooth densities and improves convergence speed.

Findings

SFP method effectively resolves Gibbs phenomenon in option pricing.

Requires fewer terms for fast error convergence.

Accurately prices deep in/out of the money options with various maturities.

Abstract

We apply a new numerical method, the singular Fourier-Pad\'e (SFP) method invented by Driscoll and Fornberg (2001, 2011), to price European-type options in L\'evy 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.

Tables

Table 1. 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	${\mathrm{𝟙}}_{S_T\ge K}$	${\mathrm{𝟙}}_{e^{x+\chi -\log K}\ge 1}$
Cash-or-Nothing Put	${\mathrm{𝟙}}_{S_T\le K}$	${\mathrm{𝟙}}_{e^{x+\chi -\log K}\le 1}$
Asset-or-Nothing Call	$S_T{\mathrm{𝟙}}_{S_T\ge K}$	$e^{x+\chi }{\mathrm{𝟙}}_{e^{x+\chi -\log K}\ge 1}$
Asset-or-Nothing Put	$S_T{\mathrm{𝟙}}_{S_T\le K}$	$e^{x+\chi }{\mathrm{𝟙}}_{e^{x+\chi -\log K}\le 1}$
Asymmetric Call	$\left(S_T^n-K^n\right){\mathrm{𝟙}}_{S_T\ge K}$	$K^n\left(e^{n\left(x+\chi -\log K\right)}-1\right){\mathrm{𝟙}}_{e^{x+\chi -\log K}\ge 1}$
Asymmetric Put	$\left(K^n-S_T^n\right){\mathrm{𝟙}}_{S_T\le K}$	$-K^n\left(e^{n\left(x+\chi -\log K\right)}-1\right){\mathrm{𝟙}}_{e^{x+\chi -\log K}\le 1}$
Symmetric Call	${\left(S_T-K\right)}^n{\mathrm{𝟙}}_{S_T\ge K}$	$\sum _{j=0}^n\left(\frac{n}{j}\right){\left(-1\right)}^{\left(n-j\right)}{e}^{j\left(x+\chi \right)+\left(n-j\right)\log K}{\mathrm{𝟙}}_{e^{x+\chi -\log K}\ge 1}$
Symmetric Put	${\left(K-S_T\right)}^n{\mathrm{𝟙}}_{S_T\le K}$	$\sum _{j=0}^n\left(\frac{n}{j}\right){\left(-1\right)}^{\left(n-j\right)}{e}^{\left(n-j\right)\left(x+\chi \right)+j\log K}{\mathrm{𝟙}}_{e^{x+\chi -\log K}\le 1}$

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

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

Table 3. 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 𝐾 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. 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 𝐾 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. 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 𝑇 50 T=50 ) and 99.9945609694213000 ( T = 100 𝑇 100 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 $\left(T=50\right)$	Error $\left(T=100\right)$
SWIFT $\left(m=0\right)$	$1.91e-01$	$2.50e-05$
COS $\left(N=35\right)$	$4.98e-01$	$2.05e+02$
SFP $\left(U=32\right)$	$2.653e-07$	$7.067e-08$
SWIFT $\left(m=1\right)$	$7.78e-09$	$3.20e-06$
COS $\left(N=70\right)$	$2.79e-08$	$2.02e-05$
SFP $\left(U=64\right)$	$2.251e-10$	$7.037e-11$

Table 6. 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 0 ( S 0 = 95 subscript 𝑆 0 95 S_{0}=95 ) and 0.007491657716010 0.007491657716010 0.007491657716010 ( S 0 = 99.999 subscript 𝑆 0 99.999 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

Table 7. 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 10.993703186728190 10.993703186728190 ( K = 90 𝐾 90 K=90 ) and 0.689027011772653 0.689027011772653 0.689027011772653 ( K = 102.336 𝐾 102.336 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

Table 8. 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 subscript 𝑆 0 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

Table 9. 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 𝐾 100 K=100 . The reference values are 19.812948843118576 19.812948843118576 19.812948843118576 ( Y = 0.5 𝑌 0.5 Y=0.5 ) and 49.790905468523860 49.790905468523860 49.790905468523860 ( Y = 1.5 𝑌 1.5 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

Table 10. 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 𝐾 100 K=100 . The reference values are 5.7851554534076321 5.7851554534076321 5.7851554534076321 ( T = 1 𝑇 1 T=1 ) and 22.318945791154533 22.318945791154533 22.318945791154533 ( T = 10 𝑇 10 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. 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 𝑇 1 T=1 . The reference values are 50.070539139715081 50.070539139715081 50.070539139715081 ( K = 50 𝐾 50 K=50 ) and 3.181555642433310 3.181555642433310 3.181555642433310 ( K = 105.453 𝐾 105.453 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. 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

Table 13. Table 13: The first c 1 , subscript 𝑐 1 c_{1}, second c 2 , subscript 𝑐 2 c_{2}, and fourth c 4 subscript 𝑐 4 c_{4} cumulants of various models

Lévy Models
BSM	$c_1=\left(r-q+\omega \right)t$ $c_2={\sigma }^{2}t,$ $c_4=0,$ $\omega =-0.5{\sigma }^2$
VG	$c_1=\left(r-q+\theta +\omega \right)t$
	$c_2=\left({\sigma }^2+\upsilon {\theta }^2\right)t$
	$c_4=3\left({\sigma }^4\upsilon +2{\theta }^4{\upsilon }^3+4{\sigma }^2{\theta }^2{\upsilon }^2\right)t$
	$\omega =\frac{1}{\upsilon }\log \left(1-\theta \upsilon -{\sigma }^2\upsilon /2\right)$
$\mathrm{CGMY}$	$c_1=\left(r-q+\omega \right)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 \left(-Y\right)G^Y\left({\left(1+\frac{1}{G}\right)}^Y-1-\frac{Y}{G}\right)+C\Gamma \left(-Y\right)M^Y\left({\left(1-\frac{1}{M}\right)}^Y-1+\frac{Y}{M}\right)\right)$
Affine Processes
Heston	$c_1=\left(r-q\right)t+\left(1-e^{-\lambda t}\right)\frac{\overline{y}-y_0}{2\lambda }-0.5\overline{y}t$
	$c_2=\frac{1}{8{\lambda }^3}(\eta t\lambda \exp (-\lambda t)(y_0-\overline{y})(8\lambda \rho -4\eta )+$
	$\lambda \rho \eta \left(1-e^{-\lambda t}\right)\left(16\overline{y}-8y_0\right)+2\overline{y}\lambda t\left(-4\rho \eta +{\eta }^2+4{\lambda }^2\right)+$
	${\eta }^2((\overline{u}-2u_0)e^{-2\lambda t}+\overline{y}(6e^{-\lambda t}-7)+2y_0)+8{\lambda }^2(y_0-\overline{y})(1-e^{-\lambda t}))$

Table 14. Table 14: Absolute error convergence for the COS method for calls under the Heston model with T = 1 𝑇 1 T=1 and T = 10 . 𝑇 10 T=10. Reference values are 5.785155450…( T = 1 𝑇 1 T=1 ) and 22.318945791…( T = 10 𝑇 10 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

Equations

$$P_N(x)=\sum _{n=0}^N{p}_n{x}^{n}and{Q}_M(x)=\sum _{m=0}^M{q}_m{x}^m,$$

$$\left(\sum _{n=0}^N{p}_n{x}^n\right)-\left(\sum _{m=0}^M{q}_m{x}^m\right)\left(\sum _{k=0}^{M+N}{b}_k{x}^k\right)$$

$$\left[\begin{array}{ccccc}{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}\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ b_{N+M}& \cdots & b_{N+2}& b_{N+1}& b_N\end{array}\right]\left[\begin{array}{c}{q}_0\\ q_1\\ ⋮\\ q_M\end{array}\right]=0.$$

$$\left[\begin{array}{c}{p}_0\\ p_1\\ ⋮\\ p_N\end{array}\right]=\left[\begin{array}{cccc}{b}_0& & & \\ b_1& b_0& & \\ ⋮& \ddots & \ddots & \\ b_N& \cdots & b_1& b_0\end{array}\right]\left[\begin{array}{c}{q}_0\\ q_1\\ ⋮\\ q_M\end{array}\right].$$

$$f(t)=\sum _{k=-\infty }^{\infty }{b}_k{e}^{ikt},with{b}_k=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }f(t)e^{-ikt}dt.$$

$$f(z)$$

$$=f^{+}(z)+f^{-}(z^{-1}),$$

$$P_N^{\pm }(z)=Q_M^{\pm }(z)f^{\pm }(z)+\mathcal{O}(z^{N+M+1}),z\to 0.$$

$$\frac{P_N^{+}(z)}{Q_M^{+}(z)}+\frac{P_N^{-}(z^{-1})}{Q_M^{-}(z^{-1})}.$$

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

$$P_N^{\pm }(z)+\mathcal{L}=Q_M^{\pm }(z)f^{\pm }(z)+\mathcal{O}(z^{U+1}),$$

$$\mathcal{L}=\sum _{s=1}^S{L}_s^{\pm }(z)\log \left(1-\frac{z}{e^{i{\zeta }_s}}\right)$$

$$f(x)=\mathfrak{Re}\left[\sum _{k=-\infty }^{\infty }{b}_k{e}^{i\frac{2\pi }{b-a}kx}\right],with{b}_k=\frac{1}{b-a}{\int }_a^{b}f(x)e^{-i\frac{2\pi }{b-a}kx}dx.$$

$$f(x)=\mathfrak{Re}\left[2\sum _{k=1}^{\infty }{b}_k{e}^{i\frac{2\pi }{b-a}kx}+b_0\right].$$

$$\mathfrak{Re}\left[2\sum _{k=1}^U{b}_k{z}^k+b_0\right].$$

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

$$P_N^{+}(z)+\sum _{s=1}^S{L}_{N_s}^{+}(z)\log \left(1-z/{\epsilon }_s\right)=\left(2\sum _{k=1}^U{b}_k{z}^k+b_0\right)Q_M^{+}(z)+\mathcal{O}(z^{U+1}),$$

$$\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}$$

$$f(x)=\mathfrak{Re}\left(\frac{P_N^{+}(z)+{\sum }_{s=1}^S{L}_{N_s}^{+}(z)\log \left(1-z/{\epsilon }_s\right)}{Q_M^{+}(z)}\right),withz=e^{i\frac{2\pi }{b-a}x}.$$

$$S_T$$

$$=S_t{e}^{(r-q+\omega )(T-t)+X_{T-t}},$$

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

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

$$\varphi (u)=\frac{1}{2}{A}^{2}u-i\gamma u+{\int }_{-\infty }^{\infty }(1-e^{iu\chi }+iu\chi \text{mathds}{1}_{|\chi |\ge 1})\nu (\mathrm{d}\chi ),\chi \in X_{T-t},$$

$$\phi (u)=\mathbb{E}[e^{iu{X}_t}]=\exp (C(u,t)+X_{0}D(u,t)).$$

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

$$d{S}_t$$

$$d{y}_t$$

$$\phi (u)$$

$$\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},$$

$$E=$$

$$F=$$

$$2\lambda \bar{y}\ge {\eta }^2,$$

Full text

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$:

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:

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

The resulting approximant is then defined as

[TABLE]

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

[TABLE]

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:

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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$:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

In the equation above, basic calculus implies that

[TABLE]

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

[TABLE]

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

[TABLE]

to obtain a simplified form:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

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

[TABLE]

and eventually,

[TABLE]

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:

[TABLE]

with

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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}.$

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]$:

[TABLE]

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:

[TABLE]

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.

Integration truncation error:

[TABLE] 2. 2.

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

[TABLE] 3. 3.

Truncated SFP series error:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

Based on the result above,

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

The parameters for the experiments are selected from the following:

[TABLE]

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.

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:

[TABLE]

and its characteristic function is given by

[TABLE]

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:

[TABLE]

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.

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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.

8.2 Affine Processes

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

[TABLE]

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.

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.

\appendices

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

[TABLE]

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