On the number of terms in the COS method for European option pricing

Abstract

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

1 Introduction

To calibrate stock price models, it is crucial to price European options quickly because stock price models are typically calibrated to given prices of liquid call and put options by minimizing the mean-square-error between model prices and given market prices. During the optimization routine, the model prices of call and put options need to be evaluated often for different model parameters.

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

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

The COS method has been extensively extended and applied, see [4, 15, 16, 18, 25,26,27, 34, 39, 45]. Other Fourier pricing techniques are discussed e.g., by [8, 28, 35, 36].

With respect to these papers we make the following three main contributions: we develop an explicit, useful and rigorous bound for N; we analyze the order of convergence of the COS method in detail; and we rigorously analyze how the Greeks of an option can be approximated by the COS method.

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

One may apply a simple triangle inequality to bound the error of the three approximations and obtain:

$$\begin{aligned} \bigg \Vert f-\mathop {{\mathop {\sum }\nolimits '}}\limits _{k=0}^{N}c_{k}e_{k}^{L}\bigg \Vert _{2}\le \big \Vert f-f_{L}\big \Vert _{2}+\Big \Vert f_{L}-\mathop {{\mathop {\sum }\nolimits '}}\limits _{k=0}^{N}a_{k}e_{k}^{L}\Big \Vert _{2}+\Big \Vert \mathop {{\mathop {\sum }\nolimits '}}\limits _{k=0}^{\infty }(a_{k}-c_{k})e_{k}^{L}\Big \Vert _{2}.\nonumber \\ \end{aligned}$$

(1.1)

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

It is well known that the series truncation error, i.e., the second term on the right-hand side of Inequality (1.1), can be bounded using integration by parts, see [7]. One contribution of this article is to use this idea in order to find an explicit and useful bound for N, provided the density of the log-returns is smooth. Our bound for N is provably large enough to ensure that the COS method converges within a predefined error tolerance. There are many financial models with smooth densities having semi-heavy tails. Examples include the Black-Scholes (BS) model, see [6], the Heston model, see [20], the Normal Inverse Gaussian (NIG) model, see [5] and the CGMY model with parameter $Y\in (0,1)$, see [2, 3, 10, 23]. The density of the log-returns in the Variance Gamma (VG) model, see [32], is not smooth for some parameters, and our methodology cannot be applied to the VG model. We also compare the solution for finding N with another solution proposed by Aimi et al. [1].

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

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

However, for fat-tailed and smooth densities, such as the density of the log-returns in the Finite Moment Log Stable (FMLS) model (see [9]), the truncation error due to (i) and (iii) becomes much more relevant compared to densities with semi-heavy tails. We show theoretically that the COS method converges at least as fast as $O(N^{-\alpha })$ for $N\rightarrow \infty$, where $\alpha >0$ is the Pareto tail index, e.g., for the FMLS model $\alpha \in (1,2)$. Empirical experiments indicate that the COS method converges for such densities as fast as $O(N^{-\alpha })$, i.e., the theoretical bound is sharp and the COS method does not converge exponentially but the order of convergence is $\alpha$.

Greeks, also known as option sensitivities, play an important role in risk management. The Greek letters Delta or Gamma respectively represent the first and second derivatives of the price of the option with respect to the current price of the underlying asset. There are formulas in the literature on how to approximate the Delta and Gamma of the option using the COS method, see [14, 25, 38]. Another contribution of this article is to provide explicit formulas for the truncation range and the number of terms for the Greeks Delta and Gamma.

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

2 Overview: the COS method for option pricing

We model the stock price over time by a semimartingale $(S_{t})_{t\ge 0}$ on a filtered probability space $(\Omega ,\mathcal {F},P,(\mathcal {F}_{t})_{t\ge 0})$. The filtration $(\mathcal {F}_{t})_{t\ge 0}$ satisfies the usual conditions and $\mathcal {F}_{0}=\{\Omega ,\emptyset \}$. We assume that there is a bank account paying continuous compound interest $r\in \mathbb {R}$ and there is a risk-neutral measure Q. All expectations are taken under Q. All densities are risk-neutral.

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

In several places, we assume that the payoff function is bounded. The prices of European call options are not bounded. If we want to approximate the price of a call option with a certain error tolerance, we need only approximate the price of a put option within that error tolerance and apply the put-call parity.

Fix some $t_{0}\in [0,T)$. The price of the European option with payoff w at time $t_{0}$ is then given by

$$\begin{aligned} e^{-r(T-t_{0})}E[w(S_{T})|\mathcal {F}_{t_{0}}]. \end{aligned}$$

(2.1)

Since we only consider European options, we will focus on the time-0 price of the option and set $t_{0}=0$ for the remainder of the article.

If we know the characteristic function $\varphi _{\log (S_{T})}$ of $\log (S_{T})$ in closed form or we are able to obtain it numerically efficiently, the COS method is able to price the European option numerically very quickly, as follows: we denote by

$$\begin{aligned} X_{T}:=\log (S_{T})-E[\log (S_{T})],\quad t\ge 0 \end{aligned}$$

the centralized log-returns. The characteristic function $\varphi$ of $X_{T}$ is then equal to

$$\begin{aligned} \varphi (u)=\varphi _{\log (S_{T})}(u)\exp (-iuE[\log (S_{T})]),\quad u\in \mathbb {R}, \end{aligned}$$

where $E[\log (S_{T})]=-i\varphi _{\log (S_{T})}^{\prime }(0).$ We assume that $X_{T}$ has a density f, but the exact structure of f need not be known. Since $E[X_T]=0$, the density of $X_{T}$ is centered around zero and it is justified to truncate the density f on a symmetric truncation range $[-L,L]$. Define

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

The time-0 price of the European option with payoff w is then given by

$$\begin{aligned} e^{-rT}E[w(S_{T})]&=e^{-rT}\int _{\mathbb {R}}w(\exp (x+E[\log (S_{T})]))f(x)dx=\int _{\mathbb {R}}v(x)f(x)dx. \end{aligned}$$

(2.2)

We need some abbreviations to discuss the COS method: suppose f is $J+1$ times continuously differentiable for $J\ge 0$. We will approximate f by cosine functions to solve the integral at the right-hand side of Eq. (2.2) numerically. We also approximate the derivatives of f by cosine functions in order to approximate the Greeks, i.e., the sensitivities of the option, numerically.

By $f^{(j)}$ we denote the $j^{th}$-derivative of f. We use the convention $f^{(0)}\equiv f$. For $L>0$, let $f_{L}^{(j)}:=1_{[-L,L]}f^{(j)}$, $j=0,\ldots ,J+1$. Suppose that $f^{(j)}$ is integrable and vanishes at $\pm \infty$. By integration by parts, the Fourier transform of $f^{(j)}$ is given by

$$\begin{aligned} u\mapsto (-iu)^{j}\varphi (u),\quad u\in \mathbb {R},\quad j=0,\ldots ,J+1. \end{aligned}$$

(2.3)

Define the basis functions

$$\begin{aligned} e_{k}^{L}(x)=1_{[-L,L]}(x)\cos \left( k\pi \frac{x+L}{2L}\right) ,\quad x\in \mathbb {R},\quad k=0,1,\ldots \end{aligned}$$

The Fourier coefficients of $f_{L}^{(j)}$ are defined by $a_{k}^{j}$ and approximated by $c_{k}^{j}$, where

$$\begin{aligned} a_{k}^{j}:=&\frac{1}{L}\int _{-L}^{L}f^{(j)}(x)e_{k}^{L}(x)dx,\\ c_{k}^{j}:=&\frac{1}{L}\int _{\mathbb {R}}f^{(j)}(x)\cos \big (k\pi \frac{x+L}{2L}\big )dx,\quad k=0,1,\ldots ,\quad j=0,\ldots ,J+1. \end{aligned}$$

We also write $a_{k}$ and $c_{k}$ instead of $a_{k}^{0}$ and $c_{k}^{0}$, respectively. Intuitively, we then have

$$\begin{aligned} f^{(j)}\approx f_{L}^{(j)}=\sum _{k=0}^{\infty }{}^{\prime }a_{k}^{j}e_{k}^{L}\approx \sum _{k=0}^{N}{}^{\prime }a_{k}^{j}e_{k}^{L}\approx \sum _{k=0}^{N}{}^{\prime }c_{k}^{j}e_{k}^{L}, \end{aligned}$$

where $\sum {}^{\prime }$ indicates that the first summand (with $k=0$) is weighted by one-half. A little analysis shows that

$$\begin{aligned} c_{k}^{j}=\frac{1}{L}\Re \bigg \{\bigg (-i\frac{k\pi }{2L}\bigg )^{j}\varphi \bigg (\frac{k\pi }{2L}\bigg )e^{i\frac{k\pi }{2}}\bigg \},\quad k=0,1,\ldots ,\quad j=0,\ldots ,J+1, \end{aligned}$$

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

$$\begin{aligned} v_{k}:=\int _{-M}^{M}v(x)e_{k}^{L}(x)dx,\quad k=0,1,\ldots \end{aligned}$$

(2.4)

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

$$\begin{aligned} \int _{\mathbb {R}}v(x)f(x)dx\approx \int _{-M}^{M}v(x)\mathop {{\mathop {\sum }\nolimits '}}\limits _{k=0}^{N}c_{k}e_{k}^{L}(x)dx=\sum _{k=0}^{N}{}^{\prime }c_{k}v_{k}. \end{aligned}$$

(2.5)

The coefficients $c_{k}$ are given in closed form when $\varphi$ is given analytically and the coefficients $v_{k}$ can also be computed explicitly in important cases, e.g., for plain vanilla European put or call options and digital options, see [14]. This makes the COS method numerically very efficient and robust.

In Lemma 2.1 we give the approximation in line (2.5) a precise meaning. To do so, we need a bound for the term

$$\begin{aligned} B_{f}(L):=\sum _{k=0}^{\infty }\frac{1}{L}\left| \int _{\mathbb {R}{\setminus }[-L,L]}f(x)\cos \left( k\pi \frac{x+L}{2L}\right) dx\right| ^{2},\quad L>0. \end{aligned}$$

(2.6)

Integrable functions f with $B_{f}(L)\rightarrow 0$, $L\rightarrow \infty$ are called COS-admissible. The class of COS-admissible densities is very large; in particular, it includes bounded densities with existing first and second moments and stable densities, see [21].

Lemma 2.1

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

$$\begin{aligned} \int _{\mathbb {R}{\setminus }[-M,M]}v(x)f(x)dx\le \frac{\varepsilon }{2}. \end{aligned}$$

(2.7)

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

$$\begin{aligned} \left\| f-f_{L}\right\| _{2}\le \frac{\varepsilon }{6\xi }\quad \text {and}\quad \sqrt{B_{f}(L)}\le \frac{\varepsilon }{6\xi }. \end{aligned}$$

(2.8)

Choose N large enough so that

$$\begin{aligned} \left\| f_{L}-\sum _{k=0}^{N}{}^{\prime }a_{k}e_{k}^{L}\right\| _{2}\le \frac{\varepsilon }{6\xi }. \end{aligned}$$

(2.9)

Then it follows that

$$\begin{aligned} \left| \int _{\mathbb {R}}v(x)f(x)dx-\sum _{k=0}^{N}{}^{\prime }c_{k}v_{k}\right| \le \varepsilon . \end{aligned}$$

Proof

[21, Cor. 8]. $\square$

Remark 2.2

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

3 On the choice of N for smooth densities

We summarize the assumptions about the density f of the log-returns in order to find explicit expressions for M, L and N. We denote by $C_{b}^{J+1}(\mathbb {R})$ the set of bounded functions from $\mathbb {R}$ to $\mathbb {R}$ which are $(J+1)$-times, continuously differentiable with bounded derivatives. By $\left\| .\right\| _{\infty }$ and $\left\| .\right\| _{2}$ we denote the supremum norm and the $L^{2}$ norm, i.e.,

$$\begin{aligned} \left\| f\right\| _{\infty }=\sup _{x\in \mathbb {R}}|f(x)|,\quad \left\| f\right\| _{2}=\sqrt{\int _{\mathbb {R}}(f(x))^{2}dx}. \end{aligned}$$

Let $\mathbb {N}_{0}=\{0\}\cup \mathbb {N}$ and $J\in \mathbb {N}_{0}$. Let $C_{1}>0$, $C_{2}>0$ and $C_{3}>0$ be suitable constants. Let $L_{0}>0$. Assume $f\in C_{b}^{J+1}(\mathbb {R})$. We say that f and its derivatives have semi-heavy tails if

$$\begin{aligned} |f^{(j)}(x)|\le C_{1}C_{2}^{j}e^{-C_{2}|x|},\quad j=0,\ldots ,J+1,\quad |x|\ge L_{0}. \end{aligned}$$

(3.1)

We say that f and its derivatives have heavy tails with index $\alpha >0$ if

$$\begin{aligned} |f^{(j)}(x)|\le C_{3}\prod _{m=1}^{j}(\alpha +m)\,|x|^{-1-\alpha -j},\quad j=0,\ldots ,J+1,\quad |x|\ge L_{0}. \end{aligned}$$

(3.2)

We suppose that f satisfies one of the following assumptions:

Assumption A1

$f\in C_{b}^{J+1}(\mathbb {R})$ and f and its derivatives have semi-heavy tails.

Assumption A2

$f\in C_{b}^{J+1}(\mathbb {R})$ and f and its derivatives have heavy tails with index $\alpha >0$.

Remark 3.1

Assume the functions $u\mapsto |u^{j}\varphi (u)|$, $j=0,\ldots ,J+1$, are integrable. By Fourier inversion we have that $f(x)=\frac{1}{2\pi }\int _{\mathbb {R}}e^{iux}\varphi (u)du$. By [17, Lemma 2.8] it follows that $f\in C_{b}^{J+1}(\mathbb {R})$.

Remark 3.2

In exceptional cases, the constants $C_{1}$ and $C_{2}$ are explicitly known; see Example 3.4. However, it should be pointed out that in Theorem 3.8 we obtain bounds for M, L and N for models with semi-heavy tails (e.g., BS, VG, Heston, NIG and CGMY) without knowing $C_{1}$ or $C_{2}$.

Remark 3.3

In Theorem 3.8 we treat models with Pareto-tails; i.e., the density for the log-returns behaves like

$$\begin{aligned} f(x)\sim |x|^{-1-\alpha },\quad x\rightarrow \pm \infty . \end{aligned}$$

(3.3)

The right-hand side of Inequality (3.2) is obtained by differentiating the right-hand side of (3.3). We assume in Theorem 3.8 that $C_{3}$ and $\alpha$ are known. The exact tail-behavior of the density of the log-returns is indeed known for the stable law, in particular for the FMLS model.

Example 3.4

Let $\sigma >0$. In the Laplace model, see [30], the centralized log-returns at maturity $T>0$ are Laplace distributed with variance $\sigma ^{2}T$. To ensure stock prices are finite, we need $\sigma \sqrt{T}<\sqrt{2}$, see [19, Example 3]. It holds that

$$\begin{aligned} |f_{\text {Lap}}^{(j)}(x)|=\frac{1}{\sqrt{2}\sigma \sqrt{T}}\left( \frac{\sqrt{2}}{\sigma \sqrt{T}}\right) ^{j}e^{-\frac{\sqrt{2}}{\sigma \sqrt{T}}|x|},\quad x\in \mathbb {R}{\setminus }\{0\},\quad j=0,1,2,\ldots \end{aligned}$$

Let $L_{0}>0$. Choose $C_{1}=\frac{1}{\sqrt{2}\sigma \sqrt{T}}$ and $C_{2}=\frac{\sqrt{2}}{\sigma \sqrt{T}}$. Then $f_{\text {Lap}}^{(j)}$ satisfies Inequality (3.1).

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

Lemma 3.5

Let $J\in \mathbb {N}_{0}$. Suppose $f\in C_{b}^{J+1}(\mathbb {R})$. It holds for $J\ge 1$ that

$$\begin{aligned} \left\| f_{L}-\sum _{k=0}^{N}{}^{\prime }a_{k}e_{k}^{L}\right\| _{2}&\le \sum _{j=1}^{J}\frac{2^{j+1}}{j\pi ^{j+1}}\frac{L^{j+\frac{1}{2}}}{N^{j}}\left( |f^{(j)}(-L)|+|f^{(j)}(L)|\right) \nonumber \\&\quad +\frac{2^{J+2}\Vert f^{(J+1)}\Vert _{\infty }}{J\pi ^{J+1}}\frac{L^{J+\frac{3}{2}}}{N^{J}} \end{aligned}$$

(3.4)

and for $J=0$ that

$$\begin{aligned} \left\| f_{L}-\sum _{k=0}^{N}{}^{\prime }a_{k}e_{k}^{L}\right\| _{2}\le \frac{4\Vert f^{(1)}\Vert _{\infty }}{\pi }\frac{L^{\frac{3}{2}}}{\sqrt{N}}. \end{aligned}$$

Proof

It holds for any $\nu >0$ by integration by parts, see [29, Eq. (1.3)], that

$$\begin{aligned} \int _{-L}^{L}f(x)e^{i\nu x}dx&= \sum _{j=0}^{J}\frac{i^{j+1}}{\nu ^{j+1}}\left( e^{-i\nu L}f^{(j)}(-L)-e^{i\nu L}f^{(j)}(L)\right) \nonumber \\&\quad +\frac{i^{J+1}}{\nu ^{J+1}}\int _{-L}^{L}f^{(J+1)}(x)e^{i\nu x}dx. \end{aligned}$$

(3.5)

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

$$\begin{aligned} |a_{k}|&= \frac{1}{L}\left| \int _{-L}^{L}f(x)\cos \left( k\pi \frac{x+L}{2L}\right) dx\right| \\&= \frac{1}{L}\left| \Re \left\{ e^{i\frac{k\pi }{2}}\int _{-L}^{L}f(x)e^{i\frac{k\pi }{2L}x}dx\right\} \right| \\&= \frac{1}{L}\bigg |\Re \bigg \{\sum _{j=0}^{J}i^{j+1}\frac{(2L)^{j+1}}{(k\pi )^{j+1}}\left( f^{(j)}(-L)-(-1)^{k}f^{(j)}(L)\right) \\&\quad +e^{i\frac{k\pi }{2}}i^{J+1}\frac{(2L)^{J+1}}{(k\pi )^{J+1}}\int _{-L}^{L}f^{(J+1)}(x)e^{i\frac{k\pi }{2L}x}dx\bigg \}\bigg |\\&\le \sum _{j=1}^{J}\frac{2^{j+1}}{\pi ^{j+1}}\frac{L^{j}}{k^{j+1}}\left( |f^{(j)}(-L)|+|f^{(j)}(L)|\right) +\frac{2^{J+2}\left\| f^{(J+1)}\right\| _{\infty }}{\pi ^{J+1}}\frac{L^{J+1}}{k^{J+1}}. \end{aligned}$$

Note that $\langle e_{k}^{L},e_{l}^{L}\rangle =L\delta _{k,l}$ for $(k,l)\ne (0,0)$ and hence

$$\begin{aligned} \left\| f_{L}-\sum _{k=0}^{N}{}^{\prime }a_{k}e_{k}^{L}\right\| _{2}=\sqrt{L\sum _{k=N+1}^{\infty }|a_{k}|^{2}}\le \sqrt{L}\sum _{k=N+1}^{\infty }|a_{k}|. \end{aligned}$$

(3.6)

By the integral test for convergence,

$$\begin{aligned} \sum _{k=N+1}^{\infty }\frac{1}{k^{j+1}}\le \int _{N}^{\infty }x^{-j-1}dx=\frac{1}{j}N^{-j},\quad j\ge 1, \end{aligned}$$

(3.7)

which implies Eq. (3.4) for $J\ge 1$. If $J=0$, apply the first Equality from (3.6) and the Inequality (3.7). $\square$

Given $L>0$, we need to find an upper bound for $\left\| f^{(j)}\right\| _{\infty }$ to estimate the series truncation error by Inequality (3.4). It follows by the inverse Fourier transform and Eq. (2.3) that

$$\begin{aligned} \Vert f^{(j)}\Vert _{\infty }\le \frac{1}{2\pi }\int _{\mathbb {R}}|u|^{j}|\varphi (u)|du,\quad j=0,1,\ldots ,J+1. \end{aligned}$$

(3.8)

Inequality (3.8) provides an explicit expression to find a bound for the term $\Vert f^{(j)}\Vert _{\infty }$, $j=0,\ldots ,J+1$ for several models.

Example 3.6

In the symmetric NIG model with parameters $\alpha >0$, $\beta =0$ and $\delta >0$, the centralized log-returns at time $T>0$ have density $f_{\text {NIG}}\in C_{b}^{\infty }(\mathbb {R})$, which can be expressed in terms of the modified Bessel-function of the third kind. The characteristic function is given by

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

see [5] and [41, Sec. 5.3.8]. We obtain, by Inequality (3.8) and using $\alpha ^{2}+u^{2}\ge u^{2}$, that

$$\begin{aligned} \Vert f_{\text {NIG}}^{(j)}\Vert _{\infty }\le \frac{\exp (T\delta \alpha )j!}{(T\delta )^{j+1}\pi },\quad j=0,1,2,\ldots \end{aligned}$$

We need Lemma 3.7 to obtain a bound for $B_{f}(L)$. We use the following abbreviation: the maximum of two real numbers x, y, is denoted by $x\vee y$.

Lemma 3.7

Assume $f\in C_{b}^{1}(\mathbb {R})$ is integrable such that $xf^{2}(x)\rightarrow 0$, $x\rightarrow \pm \infty$. Let $L>0$. If

$$\begin{aligned} f^{\prime }(x)\ge 0,\quad x\le -L\quad \text {and}\quad f^{\prime }(x)\le 0,\quad x\ge L \end{aligned}$$

(3.9)

then

$$\begin{aligned} B_{f}(L)\le \frac{1}{L}\left( \int _{\mathbb {R}{\setminus }[-L,L]}f(x)dx\right) ^{2}+\frac{8}{3}L\left( f^{2}(L)\vee f^{2}(-L)\right) . \end{aligned}$$

Proof

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

$$\begin{aligned} \int _{L}^{\infty }\left| f^{\prime }(x)\right| dx=-\int _{L}^{\infty }f^{\prime }(x)dx=f(L). \end{aligned}$$

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

$$\begin{aligned} \int _{L}^{\infty }f(x)\cos \left( k\pi \frac{x+L}{2L}\right) dx&= \underbrace{\left[ f(x)\sin \left( k\pi \frac{x+L}{2L}\right) \frac{2L}{k\pi }\right] _{L}^{\infty }}_{=0}\\&\quad -\frac{2L}{k\pi }\int _{L}^{\infty }f^{\prime }(x)\sin \left( k\pi \frac{x+L}{2L}\right) dx, \end{aligned}$$

implying

$$\begin{aligned} \frac{1}{L}\left| \int _{L}^{\infty }f(x)\cos \left( k\pi \frac{x+L}{2L}\right) dx\right| ^{2}&= \frac{4}{\pi ^{2}k^{2}}L\left| \int _{L}^{\infty }f^{\prime }(x)\sin \left( k\pi \frac{x+L}{2L}\right) dx\right| ^{2}\\&\le \frac{4}{\pi ^{2}k^{2}}Lf^{2}(L). \end{aligned}$$

Similarly,

$$\begin{aligned} \frac{1}{L}\left| \int _{-\infty }^{-L}f(x)\cos \left( k\pi \frac{x+L}{2L}\right) dx\right| ^{2}\le \frac{4}{\pi ^{2}k^{2}}Lf^{2}(-L). \end{aligned}$$

Using

$$\begin{aligned} \sum _{k=1}^{\infty }\frac{1}{k^{2}}=\frac{\pi ^{2}}{6}\quad \text {and}\quad |a+b|^{2}\le 4(a^{2}\vee b^{2}), \end{aligned}$$

we arrive at

$$\begin{aligned}&\sum _{k=0}^{\infty } \frac{1}{L}\left| \int _{\mathbb {R}{\setminus }[-L,L]}f(x)\cos \left( k\pi \frac{x+L}{2L}\right) dx\right| ^{2}\\&\quad \le \frac{1}{L}\left| \int _{\mathbb {R}{\setminus }[-L,L]}f(x)dx\right| ^{2}+4\left( \frac{4\pi ^{2}}{6\pi ^{2}}Lf^{2}(L)\vee \frac{4\pi ^{2}}{6\pi ^{2}}Lf^{2}(-L)\right) \\&\quad =\frac{1}{L}\left| \int _{\mathbb {R}{\setminus }[-L,L]}f(x)dx\right| ^{2}+\frac{8}{3}L\left( f^{2}(L)\vee f^{2}(-L)\right) . \end{aligned}$$

$\square$

In Theorem 3.8 we provide explicit formulas for M, N and L when f satisfies Assumption A1 or A2 to ensure that the COS method approximates the true price within a predefined error tolerance $\varepsilon >0$. We also include the derivatives of f in Theorem 3.8 in order to be able to approximate the sensitivities (Greeks) of the price of the option, see Sect. 5. To approximate the time-0 price of the option by Theorem 3.8, set $\ell =0$.

Theorem 3.8

(Find M, L and N) Let $v:\mathbb {R}\rightarrow \mathbb {R}$ be bounded, with $|v(x)|\le K$ for all $x\in \mathbb {R}$ and some $K>0$. Let $\varepsilon >0$ be small enough. Suppose $J\in \mathbb {N}_{0}$.

(i) Assume density f satisfies Assumption A1. For some even $n\in \mathbb {N}$ let

$$\begin{aligned} L=M=\root n \of {\frac{2K\mu _{n}}{\varepsilon }}, \end{aligned}$$

(3.10)

where $\mu _{n}$ is the $n^{th}-$moment of f, i.e., $\mu _{n}=\frac{1}{i^{n}}\left. \frac{\partial ^{n}}{\partial u^{n}}\varphi (u)\right| _{u=0}$. Let $\xi =\sqrt{2M}K$. If $J\ge 1$, let $\ell \in \{0,\ldots ,J-1\}$, $k\in \{1,\ldots ,J-\ell \}$ and

$$\begin{aligned} N&\ge \left( \frac{2^{k+2}\left\| f^{(k+1+\ell )}\right\| _{\infty }L^{k+\frac{3}{2}}}{k\pi ^{k+1}}\frac{12\xi }{\varepsilon }\right) ^{\frac{1}{k}}. \end{aligned}$$

(3.11)

If $J=0$, let $\ell =0$ and

$$\begin{aligned} N\ge \left( \frac{4\left\| f^{(1)}\right\| _{\infty }L^{\frac{3}{2}}}{\pi }\frac{6\xi }{\varepsilon }\right) ^{2}. \end{aligned}$$

(3.12)

(ii) Assume density f satisfies Assumption A2, $J\ge 1$ and f is unimodal. Let $M=\left( \frac{4C_{3}K}{\varepsilon \alpha }\right) ^{\frac{1}{\alpha }}$, $\xi =\sqrt{2M}K$ and

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

(3.13)

Let $\ell =0$ and $k\in \{1,\ldots ,J\}$ and define N as in Equation (3.11). In both cases, i) and ii), it holds that

$$\begin{aligned} \left| \int _{\mathbb {R}}v(x)f^{(\ell )}(x)dx-\sum _{k=0}^{N}{}^{\prime }c_{k}^{\ell }v_{k}\right| \le \varepsilon . \end{aligned}$$

(3.14)

Proof

We start with case (i). For $\varepsilon$ small enough, M is large enough so that Assumption A1 holds, i.e., $L_{0}\le M=L$. It follows that

$$\begin{aligned} \int _{\mathbb {R}{\setminus }[-M,M]}\big |v(x)f^{(\ell )}(x)\big |dx\le K\int _{\mathbb {R}{\setminus }[-M,M]}C_{1}C_{2}^{\ell }e^{-C_{2}|x|}(x)dx\le \frac{\varepsilon }{2}; \end{aligned}$$

(3.15)

the last inequality holds true if

$$\begin{aligned} M\ge -\frac{1}{C_{2}}\log \left( \frac{1}{KC_{1}C_{2}^{\ell -1}}\frac{\varepsilon }{4}\right) . \end{aligned}$$

(3.16)

Further,

$$\begin{aligned} \left\| f^{(\ell )}-f_{L}^{(\ell )}\right\| _{2}&\le \sqrt{2C_{1}^{2}C_{2}^{2\ell }\int _{L}^{\infty }e^{-2C_{2}x}dx}\nonumber \\&=\frac{C_{1}C_{2}^{\ell }}{\sqrt{C_{2}}}e^{-C_{2}L}\end{aligned}$$

(3.17)

$$\begin{aligned}&\le \frac{\varepsilon }{6\xi }; \end{aligned}$$

(3.18)

the last inequality holds true if

$$\begin{aligned} L\ge -\frac{1}{C_{2}}\log \left( \frac{\sqrt{C_{2}}}{C_{1}C_{2}^{\ell }}\frac{\varepsilon }{6\xi }\right) . \end{aligned}$$

(3.19)

Proposition 2 in Junike and Pankrashkin [21] shows that

$$\begin{aligned} B_{f^{(\ell )}}(L)\le \frac{2}{3}\frac{\pi ^{2}}{L^{2}}\int _{\mathbb {R}{\setminus }[-L,L]}|xf^{(\ell )}(x)|^{2}dx. \end{aligned}$$

We then have

$$\begin{aligned} \sqrt{B_{f^{(\ell )}}(L)}&\le \sqrt{\frac{4}{3}\frac{\pi ^{2}C_{1}^{2}C_{2}^{2\ell }}{L^{2}}\int _{L}^{\infty }x^{2}e^{-2C_{2}x}dx}\nonumber \\&= \frac{2\pi C_{1}C_{2}^{\ell }}{\sqrt{6C_{2}}}e^{-C_{2}L}\sqrt{1+\frac{1}{LC_{2}}+\frac{1}{2L^{2}C_{2}^{2}}}\end{aligned}$$

(3.20)

$$\begin{aligned} \le&\frac{\varepsilon }{6\xi }; \end{aligned}$$

(3.21)

the last inequality holds true if

$$\begin{aligned} L\ge -\frac{1}{C_{2}}\log \left( \left( \frac{2\pi C_{1}C_{2}^{\ell }}{\sqrt{6C_{2}}}\sqrt{1+\frac{1}{LC_{2}}+\frac{1}{2L^{2}C_{2}^{2}}}\right) ^{-1}\frac{\varepsilon }{6\xi }\right) . \end{aligned}$$

(3.22)

Assume $J\ge 1$. For $\varepsilon$ small enough, we have

$$\begin{aligned} N\ge \frac{3LC_{2}}{\pi }, \end{aligned}$$

(3.23)

because the right-hand side of Inequality (3.23) is of order $\varepsilon ^{-\frac{1}{n}}$ while the right-hand side of Inequality (3.11) is of order $\varepsilon ^{-(\frac{1}{n}+\frac{2}{kn}+\frac{1}{k})}$. By Inequality (3.23) it follows that

$$\begin{aligned}&\sum _{j=1}^{k}\frac{2^{j+1}}{j\pi ^{j+1}}\frac{L^{j+\frac{1}{2}}}{N^{j}}\left( |f^{(j+\ell )}(-L)|+|f^{(j+\ell )}(L)|\right) \\&\quad \le \sum _{j=1}^{k}\frac{2^{j+1}}{j\pi ^{j+1}}\frac{L^{j+\frac{1}{2}}}{N^{j}}\left( 2C_{1}C_{2}^{j+\ell }e^{-C_{2}L}\right) \\&\quad \le \frac{8}{\pi ^{2}}C_{1}C_{2}^{1+\ell }e^{-C_{2}L}\frac{L^{\frac{3}{2}}}{N}\sum _{j=0}^{\infty }\left( \frac{2LC_{2}}{\pi N}\right) ^{j}\\&\quad \le \frac{8}{\pi ^{2}}C_{1}C_{2}^{1+\ell }e^{-C_{2}L}\frac{\pi \sqrt{L}}{3C_{2}}\frac{1}{1-\frac{2LC_{2}}{\pi N}}\\&\quad \le \frac{8C_{1}C_{2}^{\ell }}{\pi }e^{-C_{2}L}\sqrt{L}\le \frac{\varepsilon }{12\xi }; \end{aligned}$$

the last inequality holds true if

$$\begin{aligned} L\ge -\frac{1}{C_{2}}\log \left( \frac{\pi }{8C_{1}C_{2}^{\ell }}\frac{\varepsilon }{12\xi \sqrt{L}}\right) . \end{aligned}$$

(3.24)

By Equation (3.10), M and L are of order $\varepsilon ^{-\frac{1}{n}}$ $.$ Hence, for $\varepsilon$ small enough, Inequalities (3.16), (3.19), (3.22) and (3.24) are indeed satisfied because the right-hand sides of these Inequalities are of order $\log (\varepsilon )$. By the definition of N in Inequality (3.11), we also have

$$\begin{aligned} \frac{2^{k+2}\left\| f^{(k+1+\ell )}\right\| _{\infty }}{k\pi ^{k+1}}\frac{L^{k+\frac{3}{2}}}{N^{k}}\le \frac{\varepsilon }{12\xi }. \end{aligned}$$

By Lemma 3.5 it follows that

$$\begin{aligned} \left\| f_{L}^{(\ell )}-\sum _{k=0}^{N}{}^{\prime }a_{k}^{\ell }e_{k}^{L}\right\| _{2}\le \frac{\varepsilon }{6\xi }. \end{aligned}$$

(3.25)

As $B_{f^{(\ell )}}(L)\rightarrow 0$, $L\rightarrow \infty$, $f^{(\ell )}$ is COS-admissible. Inequalities (3.15), (3.18), (3.21), (3.25) and Lemma 2.1 imply Inequality (3.14). Now assume $J=0$. By the definition of N in Inequality (3.12), Lemma 3.5 again implies Inequality (3.25). Apply Lemma 2.1 to conclude.

Next, we treat the case (ii). For $\varepsilon$ small enough, M is large enough so that Assumption A2 holds, i.e., $L_{0}\le M\le L$. The inequality

$$\begin{aligned} \int _{\mathbb {R}{\setminus }[-M,M]}\big |v(x)f(x)\big |dx\le K\int _{\mathbb {R}{\setminus }[-M,M]}f(x)dx\le \frac{2KC_{3}}{\alpha }M^{-\alpha }\le \frac{\varepsilon }{2} \end{aligned}$$

is satisfied by the definition of M. A little calculation shows that the definition of L implies

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

therefore

$$\begin{aligned} \big \Vert f-f_{L}\big \Vert _{2}&\le \sqrt{2\int _{L}^{\infty }C_{3}^{2}x^{-2-2\alpha }dx}\nonumber \\&=\frac{\sqrt{2}C_{3}}{\sqrt{1+2\alpha }}L^{-\frac{1+2\alpha }{2}}\nonumber \\&\le \frac{\varepsilon }{6\xi }. \end{aligned}$$

(3.26)

Since f is a unimodal density satisfying Assumption A2, f also satisfies the assumption of Lemma 3.7. To see this, note that the unimodality implies the assertion in line (3.9). Further,

$$\begin{aligned} |xf^{2}(x)|\le C_{3}^{2}|x|^{-1-2\alpha }\rightarrow 0,\quad x\rightarrow \pm \infty . \end{aligned}$$

Using the bound for $B_{f}(L)$ from Lemma 3.7, we obtain

$$\begin{aligned} \sqrt{B_{f}(L)}&\le \sqrt{\frac{1}{L}\left( 2\int _{L}^{\infty }C_{3}x^{-1-\alpha }dx\right) ^{2}+\frac{8}{3}LC_{3}^{2}L^{-2\alpha -2}}\nonumber \\&=\sqrt{\frac{4C_{3}^{2}}{\alpha ^{2}}L^{-2\alpha -1}+\frac{8}{3}C_{3}^{2}L^{-2\alpha -1}}\nonumber \\&=2C_{3}\sqrt{\frac{1}{\alpha ^{2}}+\frac{2}{3}}L^{-\frac{1+2\alpha }{2}}\nonumber \\&\le \frac{\varepsilon }{6\xi }, \end{aligned}$$

(3.27)

the last Inequality holds by the definition of L. For $\varepsilon$ small enough, N is large enough and we have

$$\begin{aligned} \frac{\alpha +k}{N-\frac{2(\alpha +k)}{\pi }}\le \sqrt{\frac{1}{\alpha ^{2}}+\frac{2}{3}}. \end{aligned}$$

(3.28)

Use Inequality (3.28), the definition of L and $\frac{96}{\pi ^{2}}\le 12$, to see that

$$\begin{aligned} L\ge \left( \frac{96C_{3}\xi }{\pi ^{2}\varepsilon }\frac{\alpha +k}{N-\frac{2(\alpha +k)}{\pi }}\right) ^{\frac{2}{1+2\alpha }}. \end{aligned}$$

(3.29)

Using $\prod _{m=1}^{j}(\alpha +m)\le (\alpha +k)^{j}$ for $j\le k$, it follows that

$$\begin{aligned}&\sum _{j=1}^{k}\frac{2^{j+1}}{j\pi ^{j+1}}\frac{L^{j+\frac{1}{2}}}{N^{j}}\left( |f^{(j)}(-L)|+|f^{(j)}(L)|\right) \\&\quad \le \sum _{j=1}^{k}\frac{2^{j+1}}{j\pi ^{j+1}}\frac{L^{j+\frac{1}{2}}}{N^{j}}\left( 2C_{3}(\alpha +k)^{j}L^{-1-\alpha -j}\right) \\&\quad \le \frac{8}{\pi ^{2}}C_{3}L^{-\frac{1}{2}-\alpha }\frac{\alpha +k}{N}\sum _{j=0}^{\infty }\left( \frac{2(\alpha +k)}{N\pi }\right) ^{j}\\&\quad \le \frac{8}{\pi ^{2}}C_{3}L^{-\frac{1}{2}-\alpha }\frac{\alpha +k}{N-\frac{2(\alpha +k)}{\pi }}\\&\quad \le \frac{\varepsilon }{12\xi }, \end{aligned}$$

the last inequality holds by Inequality (3.29). By the definition of N, we also have

$$\begin{aligned} \frac{2^{k+2}\left\| f^{(k+1)}\right\| _{\infty }}{k\pi ^{k+1}}\frac{L^{k+\frac{3}{2}}}{N^{k}}\le \frac{\varepsilon }{12\xi }. \end{aligned}$$

By Lemma 3.5 it follows that

$$\begin{aligned} \left\| f_{L}-\sum _{k=0}^{N}{}^{\prime }a_{k}e_{k}^{L}\right\| _{2}\le \frac{\varepsilon }{6\xi }. \end{aligned}$$

As $B_{f}(L)\rightarrow 0$, $L\rightarrow \infty$, f is COS-admissible. Apply Lemma 2.1 to conclude. $\square$

Remark 3.9

If v is a European put option with maturity T, K can be set to the strike of the option times $e^{-rT}$. The error tolerance is described by $\varepsilon$. Numerical experiments in Junike and Pankrashkin [21] suggest choosing $n\in \{4,6,8\}$ for $\mu _{n}$. According to Theorem 3.8, any $k\in \{1,\ldots ,J-\ell \}$ is allowed to define N by Inequality (3.11). In the applications, one could minimize N over all admissible k. But this could be time-consuming, and in the applications, we set k to a fixed value, e.g., for the BS model, $k=40$ is a reasonable choice, see Sect. 6.1. For other models, another choice for k might be more suitable. Bounds for $\Vert f^{(k+1)}\Vert _{\infty }$ are explicitly known for some models, e.g., the BS, NIG and FMLS models. These bounds can also be estimated numerically, e.g., for the Heston model. Section 6 contains examples indicating that the bound for N is often sharp and very useful in applications.

Remark 3.10

If a density satisfies Assumption A1, it also satisfies Assumption A2, i.e., theoretically, case (i) in Theorem 3.8 is included in case (ii). However, in (i) we do not need to know the exact tail behavior of the density, i.e., the constants $C_{1}$ and $C_{2}$ from Assumption A1, in order to estimate the truncation range because we apply Markov’s inequality to find a bound for M. This approach is not applicable to densities with heavy tails because higher moments usually do not exist. In (ii), we assume the tail behavior of the density is known precisely, i.e., we have to know $C_{3}$ and $\alpha$ in Assumption A2 to estimate M, L or N. The constants $C_{3}$ and $\alpha$ are known, for example, for the FMLS model.

4 Order of convergence

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

$$\begin{aligned} \text {err}(L,N)=\left| \int _{\mathbb {R}}v(x)f(x)dx-\sum _{k=0}^{N}{}^{\prime }c_{k}v_{k}\right| . \end{aligned}$$

Let $M=L=L(N)$. We say the error of the COS method converges with order $\rho >0$, if there is a constant $\kappa >0$ such that for all $N\in \mathbb {N}$ it holds that

$$\begin{aligned} \text {err}(L(N),N)\le \frac{\kappa }{N^{\rho }}. \end{aligned}$$

(4.1)

The error is of infinite order or converges exponentially, if Inequality (4.1) holds for any $\rho$, see [7, Sec. 2.3]. We use big O notation: for functions $g,h:\mathbb {N}\rightarrow \mathbb {R}$, we write $h(N)=O(g(N))$ as $N\rightarrow \infty$, if the absolute value of h(N) is at most a positive constant multiple of g(N) for all sufficiently large values of N.

Theorem 4.1

(Bounds for the order of convergence) Let $v:\mathbb {R}\rightarrow \mathbb {R}$ be bounded, with $|v(x)|\le K$ for all $x\in \mathbb {R}$ and some $K>0$. Assume $J\in \mathbb {N}$.

(i) Assume density f satisfies Assumption A1. Let $\beta \in \left( 0,\frac{J}{J+3}\right)$ and $\gamma >0$. If $M=L=\gamma N^{\beta }$ it holds that

$$\begin{aligned} \text {err}(L(N),N)\le O\left( N^{-J(1-\beta )+2\beta }\right) ,\quad \text {as }N\rightarrow \infty . \end{aligned}$$

(ii) Assume density f is unimodal and satisfies Assumption A2. Let $\beta \in \left( 0,\frac{J}{J+2+2\alpha }\right)$ and $\gamma >0$. If $M=L=\gamma N^{\beta }$ it holds that

$$\begin{aligned} \text {err}(L(N),N)\le O\left( N^{-\alpha \beta }\right) ,\quad \text {as }N\rightarrow \infty . \end{aligned}$$

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

Proof

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

$$\begin{aligned} \Big |\int _{\mathbb {R}}v(x)f(x)dx-\mathop {{\mathop {\sum }\nolimits '}}\limits _{k=0}^{N}c_{k}v_{k}\Big |\le&A_{0}(L)+\Vert v_{L}\Vert _{2}\left( A_{1}(L)+A_{2}(L,N)+\sqrt{B_{f}(L)}\right) , \end{aligned}$$

(4.2)

where $B_{f}(L)$ as in Eq. (2.6) and

$$\begin{aligned} A_{0}(L)&= \int _{\mathbb {R}{\setminus }[-L,L]}\big |v(x)f(x)\big |dx,\\ A_{1}(L)&= \left\| f-f_{L}\right\| _{2},\\ A_{2}(L,N)&= \big \Vert f_{L}-\mathop {{\mathop {\sum }\nolimits '}}\limits _{k=0}^{N}a_{k}e_{k}^{L}\big \Vert _{2}. \end{aligned}$$

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

$$\begin{aligned} \Vert v_{L}\Vert _{2}\le \sqrt{2L}K=\sqrt{2\gamma }KN^{\frac{\beta }{2}}. \end{aligned}$$

We now prove (i). For L large enough, we have

$$\begin{aligned} A_{0}(L)\le 2KC_{1}\int _{L}^{\infty }e^{-C_{2}x}dx=\frac{2KC_{1}}{C_{2}}e^{-C_{2}L}. \end{aligned}$$

Further, by Inequality (3.17), it holds that

$$\begin{aligned} A_{1}(L)&\le \frac{C_{1}}{\sqrt{C_{2}}}e^{-C_{2}L}. \end{aligned}$$

Assuming $L\ge 1$ and applying Inequality (3.20), it follows that

$$\begin{aligned} \sqrt{B_{f}(L)}\le&\frac{2\pi C_{1}}{\sqrt{6C_{2}}}\sqrt{1+\frac{1}{C_{2}}+\frac{1}{2C_{2}^{2}}}e^{-C_{2}L}. \end{aligned}$$

Inequality (3.4) implies

$$\begin{aligned} A_{2}(L,N)&\le \sum _{j=1}^{J}\frac{2^{j+1}}{j\pi ^{j+1}}\frac{L^{j+\frac{1}{2}}}{N^{j}}\left( 2C_{1}C_{2}^{j}e^{-C_{2}L}\right) +\frac{2^{J+2}\left\| f^{(J+1)}\right\| _{\infty }}{J\pi ^{J+1}}\frac{L^{J+\frac{3}{2}}}{N^{J}}\\&\le e^{-C_{2}L}\sqrt{L}\sum _{j=1}^{J}\frac{2^{j+2}C_{1}C_{2}^{j}\gamma ^{j}}{j\pi ^{j+1}}+\frac{2^{J+2}\left\| f^{(J+1)}\right\| _{\infty }}{J\pi ^{J+1}}\frac{L^{J+\frac{3}{2}}}{N^{J}}, \end{aligned}$$

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

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

$$\begin{aligned} \Big |\int _{\mathbb {R}}v(x)f(x)dx-\mathop {{\mathop {\sum }\nolimits '}}\limits _{k=0}^{N}c_{k}^{L}v_{k}^{M}\Big |&\le b_{1}e^{-C_{2}\gamma N^{\beta }}+b_{2}N^{\frac{\beta }{2}}\big (b_{3}N^{\frac{\beta }{2}}e^{-C_{2}\gamma N^{\beta }}+b_{4}N^{\beta (J+\frac{3}{2})-J}\big )\nonumber \\&\le O\left( N^{-J(1-\beta )+2\beta }\right) ,\quad N\rightarrow \infty . \end{aligned}$$

(4.3)

$\beta <\frac{J}{J+3}$ implies $-J(1-\beta )+2\beta <0$ and the right-hand side of (4.3) converges to zero for $N\rightarrow \infty$. We show (ii). It holds for L large enough using Inequalities (3.26) and (3.27) that

$$\begin{aligned} A_{0}(L)&\le 2KC_{3}\int _{L}^{\infty }x^{-1-\alpha }dx\le \frac{2KC_{3}}{\alpha }L^{-\alpha },\\ A_{1}(L)&\le \frac{\sqrt{2}C_{3}}{\sqrt{1+2\alpha }}L^{-\frac{1}{2}-\alpha },\\ \sqrt{B_{f}(L)}&\le 2C_{3}\sqrt{\frac{1}{\alpha ^{2}}+\frac{2}{3}}L^{-\frac{1}{2}-\alpha }. \end{aligned}$$

Inequality (3.4) and Assumption A2 imply for N large enough and for some suitable constants $a_{1},\ldots ,a_{J}>0$ that

$$\begin{aligned} A_{2}(L,N)\le&\sum _{j=1}^{J}a_{j}N^{\beta (-\frac{1}{2}-\alpha )-j}+\frac{2^{J+2}\left\| f^{(J+1)}\right\| _{\infty }}{J\pi ^{J+1}}\frac{L^{J+\frac{3}{2}}}{N^{J}}. \end{aligned}$$

By Inequality (4.2), it holds for some suitable constants $b_{1},\ldots ,b_{4}>0$ that

$$\begin{aligned}&\Big |\int _{\mathbb {R}}v(x)f(x)dx-\mathop {{\mathop {\sum }\nolimits '}}\limits _{k=0}^{N}c_{k}v_{k}\Big |\\&\quad \le b_{1}N^{-\alpha \beta }+b_{2}N^{\frac{\beta }{2}}\bigg (b_{3}N^{-\beta (\alpha +\frac{1}{2})}+\sum _{j=1}^{J}a_{j}N^{\beta (-\frac{1}{2}-\alpha )-j}+b_{4}N^{\beta (J+\frac{3}{2})-J}\bigg )\\&\quad \le b_{1}N^{-\alpha \beta }+b_{2}\bigg (b_{3}N^{-\beta \alpha }+\sum _{j=1}^{J}a_{j}N^{-\beta \alpha -j}+b_{4}N^{\beta \left( J+2\right) -J}\bigg )\\&\quad \le O\left( N^{-\alpha \beta }\right) ,\quad N\rightarrow \infty . \end{aligned}$$

The last inequality can be seen as follows: as $\beta \le \frac{J}{J+2+\alpha }$, it follows that

$$\begin{aligned} -\alpha \beta \ge \beta (J+2)-J. \end{aligned}$$

$\square$

Remark 4.2

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

$$\begin{aligned} -J(1-\beta )+2\beta \rightarrow -\infty ,\quad J\rightarrow \infty . \end{aligned}$$

Remark 4.3

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

Remark 4.4

Theorem 4.1 cannot be applied to densities that are non-differentiable, e.g., the density of the VG model is non-differentiable for some parameters. To improve the order of convergence of the COS method if the density of the log-returns is non-differentiable, Ruijter et al. [38] apply spectral filters and consider the filter-COS method, $\sum _{k=0}^{N}{}^{\prime }\hat{s}(\frac{k}{N})c_{k}v_{k}$, where $\hat{s}$ is a spectral filter, i.e., a smooth function with support $[-1,1]$ and $\hat{s}(0)=1$. For an analysis of the order of convergence and some error bounds for the filter-COS method, see [38] and references therein.

5 On the Greeks

The Greeks or sensitivities of a European option play an important role in hedging and risk management. The most important Greeks are Delta and Gamma, which are the first and second derivative of the price of a European option with respect to the current underlying price $S_{0}>0$.

Fang and Oosterlee [14, Remark 3.2] state formulas for the approximation of Delta and Gamma by the COS method. We proof these formulas and discuss how to choose M, L and N for the Greeks.

In this section we assume $S_{t}=S_{0}\bar{S}_{t}$ for some stochastic process $(\bar{S}_{t})_{t\ge 0}$, which does not depend on $S_{0}$ anymore. This assumption is a very typical one, see [31, Sec. 3.1.2]. As in Sect. 2, we consider a European option with maturity $T>0$ and payoff $w(S_{T})$ for some $w:[0,\infty )\rightarrow \mathbb {R}$. Let $\eta :=E[\log (S_{T})]-\log (S_{0})$. Then $\eta$ does not depend on $S_{0}$. Define

$$\begin{aligned} v(x,s):=e^{-rT}w(\exp (x+\log (s)+\eta )),\quad x\in \mathbb {R},\quad s>0. \end{aligned}$$

(5.1)

The time-0 price of the European option is then given by

$$\begin{aligned} e^{-rT}E[w(S_{T})]&=\int _{\mathbb {R}}v(x,S_{0})f(x)dx, \end{aligned}$$

where, as before, f is the density of the centralized log-returns. Delta and Gamma are defined by

$$\begin{aligned} \frac{\partial ^{\ell }}{\partial S_{0}^{\ell }}\int _{\mathbb {R}}v(x,S_{0})f(x)dx,\quad \ell =1,2 \end{aligned}$$

(5.2)

if the partial derivatives exist. The next lemma provides some conditions to interchange integration and differentiation in Eq. (5.2).

Lemma 5.1

Let w be bounded. Let v be defined as in Eq. (5.1). Assume $J\in \mathbb {N}_{0}$ and density f satisfies Assumption A1. It follows that

$$\begin{aligned} \frac{\partial }{\partial S_{0}}\int _{\mathbb {R}}v(x,S_{0})f(x)dx=-\frac{1}{S_{0}}\int _{\mathbb {R}}v(x,S_{0})f^{(1)}(x)dx \end{aligned}$$

and if $J\ge 1$

$$\begin{aligned} \frac{\partial ^{2}}{\partial S_{0}^{2}}\int _{\mathbb {R}}v(x,S_{0})f(x)dx=\frac{1}{S_{0}^{2}}\int _{\mathbb {R}}v(x,S_{0})\big (f^{(1)}(x)+f^{(2)}(x)\big )dx. \end{aligned}$$

Proof

Let $K>0$ such that v is bounded by K. Let $\underline{s},\overline{s}\in \mathbb {R}$ such that $0<\underline{s}<S_{0}<\overline{s}$. Let

$$\begin{aligned} h(x,s):=v(x,1)f(x-\log (s)),\quad (x,s)\in \mathbb {R}\times (\underline{s},\overline{s}). \end{aligned}$$

Then $x\mapsto h(x,s)$ is integrable for all $s\in (\underline{s},\overline{s})$ and the partial derivative

$$\begin{aligned} \frac{\partial }{\partial s}h(x,s)=-\frac{v(x,1)}{s}f^{(1)}(x-\log (s)) \end{aligned}$$

exists for all $(x,s)\in \mathbb {R}\times (\underline{s},\overline{s})$. Define $x_{0}:=L_{0}+|\log (\overline{s})|+|\log (\underline{s})|$. Then,

$$\begin{aligned} |x-\log (s)|\ge L_{0},\quad (x,s)\in \mathbb {R}{\setminus }[-x_{0},x_{0}]\times (\underline{s},\overline{s}). \end{aligned}$$

Let

$$\begin{aligned} m(x)&:={\left\{ \begin{array}{ll} K\underline{s}^{-1}\left\| f^{(1)}\right\| _{\infty } &{} ,x\in [-x_{0},x_{0}]\\ K\underline{s}^{-C_{2}-1}C_{1}C_{2}e^{C_{2}x} &{} ,x<-x_{0}\\ K\underline{s}^{-1}\overline{s}^{C_{2}}C_{1}C_{2}e^{-C_{2}x} &{} ,x>x_{0}. \end{array}\right. } \end{aligned}$$

Then $|\frac{\partial }{\partial s}h(x,s)|\le m(x)$ for all $(x,s)\in \mathbb {R}\times (\underline{s},\overline{s})$ and $x\mapsto m(x)$ is integrable. Interchanging differentiation and integration is allowed by the dominated convergence theorem, see e.g., [17, Lemma 2.8], and it follows that

$$\begin{aligned} \frac{\partial }{\partial S_{0}}\int _{\mathbb {R}}v(x,S_{0})f(x)dx&=\frac{\partial }{\partial S_{0}}\int _{\mathbb {R}}v(x,1)f(x-\log (S_{0}))dx\\&=-\frac{1}{S_{0}}\int _{\mathbb {R}}v(x,1)f^{(1)}(x-\log (S_{0}))dx\\&=-\frac{1}{S_{0}}\int _{\mathbb {R}}v(x,S_{0})f^{(1)}(x)dx. \end{aligned}$$

If $J\ge 1$, f is twice differentiable. Apply the arguments above to $f^{(1)}$ to conclude. $\square$

In Theorem 5.2 we provide explicit formulas for M, N and L when f satisfies Assumption A1 to ensure that the COS method approximates the time-0 price and the Greeks Delta and Gamma within a predefined error tolerance $\varepsilon >0$. One can use the same parameters M, N and L to obtain both the price and the Greeks. We define $v_{k}:=\int _{-M}^{M}v(x,S_{0})e_{k}^{L}(x)dx$ as in Equation (2.4).

Theorem 5.2

(M, L and N for the time-0 price, Delta and Gamma) Let w be bounded. Let v be defined as in Equation (5.1). Let $\varepsilon >0$ be small enough. Let $\gamma =\min \big \{\varepsilon ,\,\varepsilon S_{0},\frac{\varepsilon S_{0}^{2}}{2}\big \}$. Suppose $J\ge 3$. Assume density f satisfies Assumption A1. For some even $n\in \mathbb {N}$ define

$$\begin{aligned} L=M=\root n \of {\frac{2K\mu _{n}}{\gamma }}, \end{aligned}$$

(5.3)

where $\mu _{n}$ is the $n^{th}-$moment of f, i.e., $\mu _{n}=\frac{1}{i^{n}}\left. \frac{\partial ^{n}}{\partial u^{n}}\varphi (u)\right| _{u=0}$. Let $\xi =\sqrt{2M}K$, $k\in \{1,\ldots ,J-2\}$ and

$$\begin{aligned} N&\ge \left( \frac{2^{k+2}\left( \max _{\ell \in \{0,1,2\}}\left\| f^{(k+1+\ell )}\right\| _{\infty }\right) L^{k+\frac{3}{2}}}{k\pi ^{k+1}}\frac{12\xi }{\gamma }\right) ^{\frac{1}{k}}. \end{aligned}$$

(5.4)

It follows that the time-0 price and the Greeks Delta and Gamma can be approximated by the COS method, i.e.,

$$\begin{aligned} \bigg |\int _{\mathbb {R}}v(x,S_{0})f(x)dx-\sum _{k=0}^{N}{}^{\prime }c_{k}v_{k}\bigg |&\le \varepsilon ,\\ \bigg |\frac{\partial }{\partial S_{0}}\int _{\mathbb {R}}v(x,S_{0})f(x)dx-\bigg (-\frac{1}{S_{0}}\sum _{k=0}^{N}{}^{\prime }c_{k}^{1}v_{k}\bigg )\bigg |&\le \varepsilon ,\\ \bigg |\frac{\partial ^{2}}{\partial S_{0}^{2}}\int _{\mathbb {R}}v(x,S_{0})f(x)dx-\frac{1}{S_{0}^{2}}\sum _{k=0}^{N}{}^{\prime }(c_{k}^{1}+c_{k}^{2})v_{k}\bigg |&\le \varepsilon . \end{aligned}$$

Proof

Theorem 3.8 ensures that the time-0 price can be approximated by the COS method. By Lemma 5.1 and Theorem 3.8, it holds that

$$\begin{aligned}&\bigg |\frac{\partial }{\partial S_{0}}\int _{\mathbb {R}}v(x,S_{0})f(x)dx -\bigg (-\frac{1}{S_{0}}\sum _{k=0}^{N}{}^{\prime }c_{k}^{1}v_{k}\bigg )\bigg |\\&\quad =\frac{1}{S_{0}}\left| \int _{\mathbb {R}}v(x,S_{0})f^{(1)}(x)dx-\sum _{k=0}^{N}{}^{\prime }c_{k}^{1}v_{k}\right| \\&\quad \le \frac{\gamma }{S_{0}}\le \varepsilon . \end{aligned}$$

Using the triangle inequality, we can see that

$$\begin{aligned}&\bigg |\frac{\partial ^{2}}{\partial S_{0}^{2}}\int _{\mathbb {R}}v(x,S_{0})f(x)dx -\frac{1}{S_{0}^{2}}\sum _{k=0}^{N}{}^{\prime }(c_{k}^{1}+c_{k}^{2})v_{k}\bigg |\\&\quad = \frac{1}{S_{0}^{2}}\bigg |\int _{\mathbb {R}}v(x,S_{0})\big (f^{(1)}(x)+f^{(2)}(x)\big )dx-\sum _{k=0}^{N}{}^{\prime }(c_{k}^{1}+c_{k}^{2})v_{k}\bigg |\\&\quad \le \frac{1}{S_{0}^{2}}\bigg (\bigg |\int _{\mathbb {R}}v(x,S_{0})f^{(1)}(x)dx-\sum _{k=0}^{N}{}^{\prime }c_{k}^{1}v_{k}\bigg |\\&\qquad +\bigg |\int _{\mathbb {R}}v(x,S_{0})f^{(2)}(x)dx-\sum _{k=0}^{N}{}^{\prime }c_{k}^{2}v_{k}\bigg |\bigg )\\&\quad \le \frac{2\gamma }{S_{0}^{2}}\le \varepsilon . \end{aligned}$$

$\square$

6 Numerical experiments

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

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

6.1 BS model

We consider the BS model with volatility $\sigma >0$ and maturity $T>0$. The density $f_{\text {BS}}$ of the log-returns is normally distributed and belongs to the family of stable laws. An upper bound for $\Vert f_{\text {BS}}^{(k+1)}\Vert _{\infty }$ can be obtained directly from Inequality (6.3) setting $\alpha =2$ and $c=\frac{\sigma \sqrt{T}}{\sqrt{2}}$. Let $n\in \mathbb {N}$ be even. In the BS model, the $n^{th}$ moment of the centralized log-returns is given by

$$\begin{aligned} \mu _{n}=(\sigma \sqrt{T})^{n}\left( (n-1)(n-3)(n-5)\cdot \cdot \cdot 3\cdot 1\right) . \end{aligned}$$

The formula for N in Inequality (3.11) does not depend on the volatility $\sigma \sqrt{T}$ if we set $\ell =0$ and if we bound $\Vert f_{\text {BS}}^{(k+1)}\Vert _{\infty }$ using Inequality (6.3). In Fig. 1, we show the dependency of N on k. At the beginning, the number of terms N decreases sharply as k increases and stabilizes approximately for $k\ge 40$. We also see that N increases as the error tolerance $\varepsilon$ decreases or the bound K increases.

How sharp is the bound for N in Theorems 3.8 and 5.2? We make the following experiment. Consider a put option with parameters

$$\begin{aligned} \varepsilon =10^{-8},\quad \sigma =0.2,\quad T=1,\quad S_{0}=K=100,\quad r=0. \end{aligned}$$

(6.1)

We set $n=8$ and $k=40$ to obtain $M=L=6.94$ and $N=218$ by Theorem 5.2. Other values for k and N are reported in Table 1. Theorem 5.2 indicates that M,  L and N serves to approximate both the time-0 price and the Greeks Delta and Gamma using the COS method. This can be confirmed by an experiment: $N_{\min }=120$ is the minimal number of terms such that the absolute differences of the approximation by the COS method and the closed form solution by the Black-Scholes formula for time-0 price, Delta and Gamma, are less than the error tolerance. N is about twice as larger as $N_{\min }$.

How can the number of terms N be estimated if there are no closed form solutions available for the bounds of the derivatives of the density of the log-returns? We suggest solving the right-hand side of Inequality (3.8) numerically to find a bound for $\Vert f^{(k+1)}\Vert _{\infty }$. Here, we use R’s integrate function with default values. The CPU time of the COS method and the numerical integration routine to obtain a bound for $\Vert f^{(k+1)}\Vert _{\infty }$ are of similar magnitude; see Table 1.

The value of N does not change when using numerical integration to obtain a bound for $\Vert f^{(k+1)}\Vert _{\infty }$ compared to the closed form solution for the bound of $\Vert f^{(k+1)}\Vert _{\infty }$.

Table 1 Approximation of the time-0 price of a put option by the COS method: N depending on k by Inequality (3.11)

Fig. 1

N depending on k by Inequality (3.11) for various $\varepsilon$ and K. We set $n=8$ and $\ell =0$. In the BS model, N does not depend on the standard deviation of the log-returns, i.e., on the value of $\sigma \sqrt{T}$

6.1.1 Last coefficient rule-of-thumb

By Boyd [7, Sec. 2.12], the following rule is called the Last Coefficient Rule-of-Thumb: The series truncation error approximating $f_{L}$ by a finite cosine series is bounded by $\sum _{k>N}|a_{k}|$. If N is large enough, the order of magnitude of the series truncation error is expected to scale like $|a_{N}|$. [1, Sec. 3.2] propose to select the number of terms for the BS model with volatility $\sigma >0$ and maturity T by the smallest natural number $N_{F}$ satisfying the following inequality

$$\begin{aligned} |a_{N_{F}}|\approx |c_{N_{F}}|=\frac{2}{b-a}e^{-\frac{1}{2}(\frac{\sigma \pi }{b-a})^{2}TN_{F}^{2}}\le \varepsilon _{N_{F}}, \end{aligned}$$

(6.2)

where $[a,b]=[-L,L]$ is the truncation range and $\varepsilon _{N_{F}}$ is some error tolerance. The solution to find the number of terms by (6.2) works even better than Inequality (3.11) if $\varepsilon _{N_{F}}$ is small enough, see Table 2. However, the rule by Inequality (6.2) does not work if $2^{-1}(b-a)\varepsilon _{N_{F}}\ge 1$ because we would then choose $N_{F}=0$.

Table 2 Last coefficient rule-of-thumb applied to a put option with $K=S_{0}=100$, $T=1$, $r=0$ and different volatilities

6.2 Heston model

There is an integral expression for the density of the log-returns in the Heston model and the density is smooth and has semi-heavy tails, see [12]. In this section, we provide numerical experiments under two different parameter sets, M1 and M2. The parameter set M1 is taken from [14] and the parameter set M2 is borrowed from [42]. We compute the time-0 prices of put options for different strikes K and different maturities T. Reference prices are obtained by the Carr–Madan formula.

Table 3 compares N obtained by Inequality (3.11) and the minimal $N_{\min }$ such that the absolute difference of the approximation by the COS method and the reference price is less than the error tolerance $\varepsilon$. On average, N is about six times larger than $N_{\min }$. The CPU time using N instead of $N_{\min }$ increases roughly by the factor three. To obtain N, we estimate $\Vert f^{(k+1)}\Vert _{\infty }$ by solving the right-hand side of Inequality (3.8) by numeric integration using R’s function integrate with default values. The CPU times of the COS method and the numerical integration routine are of similar magnitudes.

Table 3 Number of terms N and CPU time to approximate the time-0 price of a put option by the COS method depending on the strike K, the maturity T and the parameters of the model

6.3 VG model

Using the VG model as an example, this section shows that Theorem 3.8 does not help in finding the number of terms N for the COS method if the density of the log-returns is not sufficiently smooth. The density of the VG model at time $T>0$ has semi-heavy tails, see [2, Example 7.5]. It can be expressed by means of the Whittaker function and is $(J+1)$-times continuously differentiable if $J+2<\frac{2T}{\nu }$, see [22].

Let $f_{\text {VG}}$ denote the density of the log-returns in the VG model at maturity $T>0$ with parameters $\sigma >0$, $\theta \in \mathbb {R}$ and $\nu >0$. If $T<\frac{\nu }{2}$, the density $f_{\text {VG}}$ is unbounded and Theorem 3.8 cannot be used to find N. If $T\in \left( \nu ,\frac{3\nu }{2}\right)$, the derivative of $f_{\text {VG}}$ is continuous, but the second derivative of $f_{\text {VG}}$ is unbounded, see [22, Theorem. 4.1 and Theorem. 6.1].

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

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

By Eq. (3.10), we set $M=L=0.91$ using $n=4$ moments. By numerically optimizing the derivative of the density $f_{\text {VG}}$, we obtain $\Vert f^{(1)}\Vert _{\infty }=218$, and Theorem 3.8 suggests $N\approx 4\cdot 10^{14}$.

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

6.4 FMLS model

The stable law has been used to model stock prices since Mandelbrot [33] and Fama [13]. A representation of stable densities by special functions can be found in Zolotarev [47]. The density $f_{\text {Stable}}\in C_{b}^{\infty }(\mathbb {R})$ of the family of stable distributions with stability parameter $\alpha \in (0,2]$, skewness $\beta \in [-1,1]$, scale $c>0$ and location $\theta \in \mathbb {R}$ has the characteristic function

$$\begin{aligned} u\mapsto&\exp \left( iu\theta -|uc|{}^{\alpha }(1-i\beta \text {sgn}(u)\Phi _{\alpha }(u)\right) ,\\&\Phi _{\alpha }(u)={\left\{ \begin{array}{ll} \tan \frac{\pi \alpha }{2} &{} ,\alpha \ne 1\\ -\frac{2}{\pi }\log (|u|) &{} ,\alpha =1, \end{array}\right. } \end{aligned}$$

see [46]. It follows by Inequality (3.8) that

$$\begin{aligned} \Vert f_{\text {Stable}}^{(j)}\Vert _{\infty }&\le \frac{1}{\pi \alpha c^{j+1}}\int _{0}^{\infty }t^{\frac{j+1}{\alpha }-1}e^{-t}dt\nonumber \\&\le \frac{\Gamma \left( \frac{j+1}{\alpha }\right) }{\pi \alpha c^{j+1}},\quad j=0,1,2,\ldots , \end{aligned}$$

(6.3)

where $\Gamma$ denotes the gamma function. The density of the stable law is unimodal, see [44]. Let $F_{\text {Stable}}$ be the cumulative distribution function for the stable density $f_{\text {Stable}}$. By Samorodnitsky and Taqqu [40, Property 1.2.15] it holds that

$$\begin{aligned} \lim _{x\rightarrow \infty }x^{\alpha }(1-F_{\text {Stable}}(x))&=C_{\alpha }\frac{1+\beta }{2}c^{\alpha },\\ \lim _{x\rightarrow \infty }x^{\alpha }(F_{\text {Stable}}(-x))&=C_{\alpha }\frac{1-\beta }{2}c^{\alpha }, \end{aligned}$$

where

$$\begin{aligned} C_{\alpha }={\left\{ \begin{array}{ll} \frac{1-\alpha }{\Gamma (2-\alpha )\cos \big (\frac{\pi \alpha }{2}\big )} &{},\alpha \ne 1\\ \frac{2}{\pi } &{},\alpha =1. \end{array}\right. } \end{aligned}$$

For stable densities we therefore suggest to set $C_{3}$ in Theorem 3.8 at least as large as $\alpha C_{\alpha }\frac{1+|\beta |}{2}c^{\alpha }$ to obtain M and L.

The FMLS model with parameters $\sigma >0$ and $\alpha \in (1,2)$ describes the log-returns by a stable process with maximum negative skewness. The centralized log-returns in the FMLS model at time $T>0$ are stably distributed with stability parameter $\alpha$, skewness $\beta =-1$, scale $c=\sigma T^{\frac{1}{\alpha }}$ and location $\theta =0$.

The density of the log-returns in the FMLS model has a heavy left tail with Pareto index $\alpha$, i.e., the left tail decays like $|x|^{-1-\alpha }$, $x\rightarrow -\infty$, but the right tail decays exponentially, see [9]. This makes the FMLS very attractive from a theoretical point of view: put and call option prices and all moments of the underlying stock $S_{T}$ exist. The expectation of the log-returns also exists, but the log-returns do not have finite variance. Fitting the FMLS model to real market data shows very satisfactory results. Carr and Wu [9] calibrated the FMLS model to real market data and estimated $\alpha =1.5597$ and $\sigma =0.1486$. We test the formulas for M, L and N with these values for $\alpha$ and $\sigma$ for a European call option with the following parameters:

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

The reference price is 9.7433708 by the Carr–Madan formula, and we confirm the reference price by the COS method with $M=L=10^{5}$ and $N=10^{7}$.

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

Fig. 2

Density of the log-returns in the FMLS model and asymptotic tail behavior

To apply the COS method we define by Theorem 3.8, $M=69$, $L=176$ and $N=5815$ setting $k=40$. According to Theorem 3.8, any other choice for k is possible, but another value for k does not significantly improve N.

With these parameters, the COS method prices the call option within the error tolerance in about 1.5 ms. The minimal N to stay below the error tolerance is $N_{\min }=1200$ and the CPU time using $N_{\min }$ is about 0.4 ms.

We also apply the Carr–Madan formula with the “default parameters”, see [8], which are also recommended by Madan and Schoutens [31, Sec. 3.1], i.e., 4096 terms, a damping factor equal to 1.5 and a Fourier truncation range of 1024. With these parameters, the Carr–Madan formula prices the option within the error tolerance in about 1.1 ms. We also implemented the Carr–Madan formula using the fast Fourier transform but we found no speed advantage, compare also with Crisóstomo [11].

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

6.5 Order of convergence

In this section, we empirically analyze the order of convergence of the COS method for the BS model and the FMLS model and compare the results with Theorem 4.1.

The order of convergence in (4.1) can be estimated straightforwardly by a simulation, see [24]. As

$$\begin{aligned} \log \left( \frac{\kappa }{N^{\rho }}\right) =\log (\kappa )-\rho \log (N), \end{aligned}$$

the negative slope of a straight line obtained from a log-log plot of the error $\text {err}(L(N),N)$ against N can be used as an indicator for $\rho$. As in Sect. 4, we assume $M=L$ in this section.

6.5.1 BS model

In the BS model with parameter $\sigma >0$, the density of the log-returns is arbitrarily smooth, and the tails decay even faster than exponentially. In Fig. 3 we analyze how the error of the COS method behaves in the BS model for large N and for different choices of L.

We consider a call option with parameters $\sigma =0.2,$ $S_{0}=K=100$, $r=0$ and $T=1$. We see in Fig. 3 that the COS method seems to converge exponentially at the beginning for moderate N if we choose L constant, i.e., independent of N. But for constant L, e.g., $L=4\sigma$ or $L=6\sigma$, the COS method does not converge for $N\rightarrow \infty$ but the error remains constant for a certain level of N. This can be explained by Inequality (1.1): while the second term on the right-hand side of Inequality (1.1) converges to zero for $N\rightarrow \infty$ and L fixed, the first and third terms on the right-hand side of Inequality (1.1) do not improve as N is increased but L is kept constant.

Fig. 3

Order of convergence of the COS method for a call option in the BS model with different choices for the truncation range $[-L,L]$

For large L, such as $L=20\sigma$, this effect disappears somewhat because the tails of the Gaussian distribution decay so rapidly that, up to fixed-precision arithmetic,¹ the Gaussian density can be thought of as a density with finite support. Using arbitrary-precision arithmetic instead should show that even for $L=20\sigma$, the error of the COS method does not converge to zero for $N\rightarrow \infty$.

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

6.5.2 FMLS model

We test Theorem 4.1(ii) for the density of the log-returns in the FMLS model, which belongs to the family of stable densities and has heavy tails. For the FMLS model we use the parameters $\alpha =1.5597$ and $\sigma =0.1486$: these values are taken from Carr and Wu [9], who calibrated the FMLS model to real market data. We use the same reference price for the time-0 price of a European call option with maturity $T=1$, $S_{0}=K=100$ and $r=0$ as in Sect. 6.4.

Fig. 4

Order of convergence of the COS method with different choices for the truncation range $[-L,L]$ for the FMLS model

We compute $L_{\text {optimal}}$ for a fixed $N\in \mathbb {N}$ minimizing $\text {err}(L,N)$, i.e., for each N we choose the truncation range such that the error of the COS method is minimal. In particular, for each $N\in \{2^{m},\,m=4,5,\ldots ,20\}$ we define

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

where

$$\begin{aligned} \mathbb {L}=\{\exp (0.07m),\quad m=0,1,\ldots ,200\}, \end{aligned}$$

is a sufficiently fine grid of the interval $[1,10^{6}]$.

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

If we choose $L=\frac{1}{100}N$, the empirical order of convergence is 1.57, i.e., the error behaves like $O(N^{-1.57})$ for large N, which is very close to its theoretical bound of $O(N^{-1.5597})$. The empirical order of convergence does not differ much if we choose $L_{\text {optimal}}$ instead of $L=\frac{1}{100}N$.

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

In Fig. 4 we also plot the empirical order of convergence of the Carr–Madan formula for the FMLS model with damping factor of 0.7 and a Fourier truncation range of 1200. The Figure indicates an exponential order of convergence for the Carr–Madan formula.

7 Conclusions

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

In both cases, we found explicit and useful bounds for L and N and showed in numerical experiments the usefulness of these formulas in applications to obtain the time-0 price of an option and the Greeks Delta and Gamma. The densities of the log-returns are smooth for many models in finance, such as the BS, NIG, Heston and FMLS models.

If the density is not differentiable, Theorem 3.8 cannot be used to find a bound for N. If the density is only differentiable a few times, which is the case for the VG model for some parameters and short maturities, our bound for N is too large to be useful in most practical applications.

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