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

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.


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 Junike and Pankrashkin [21], the choice of N is largely based on trial and error.
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 a k e L k , where a k are Fourier coefficients of f L and e L k are cosine basis functions.iii) Approximate a k by some coefficients c k which can be obtained directly from the characteristic function of f .Thus, to apply the COS method, two decisions must to be made: find a suitable truncation range [−L, L] and identify the number N of cosine functions.
One may apply a simple triangle inequality to bound the error of the three approximations and obtain: . ( 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 Boyd [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 Black and Scholes [6], the Heston model, see Heston [20], the Normal Inverse Gaussian (NIG) model, see Barndorff-Nielsen [5] and the CGMY model with parameter Y ∈ (0, 1), see Carr et al. [10], Albin and Sundén [2], Küchler and Tappe [23], Asmussen [3].The density of the log-returns in the Variance Gamma (VG) model, see Madan et al. [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 Carr and Madan [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 Carr and Wu [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 −α ) for N → ∞, where α > 0 is the Pareto tail index, e.g., for the FMLS model α ∈ (1,2).Empirical experiments indicate that the COS method converges for such densities as fast as O(N −α ), i.e., the theoretical bound is sharp and the COS method does not converge exponentially but the order of convergence is α.
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,38,25].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: Section 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 Sections 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.

Overview: the COS method for option pricing
We model the stock price over time by a semimartingale (S t ) t≥0 on a filtered probability space (Ω, F , P, (F t ) t≥0 ).The filtration (F t ) t≥0 satisfies the usual conditions and F 0 = {Ω, ∅}.We assume that there is a bank account paying continuous compound interest r ∈ 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, ∞) → R. For example, a European put option with strike K > 0 can be described by the payoff w(x) = max(K − x, 0), x ≥ 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 ∈ [0, T ).The price of the European option with payoff w at time t 0 is then given by e 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 ϕ log(ST ) 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 where E[log(S T )] = −iϕ ′ log(ST ) (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 The time-0 price of the European option with payoff w is then given by We need some abbreviations to discuss the COS method: suppose f is J + 1 times continuously differentiable for J ≥ 0. We will approximate f by cosine functions to solve the integral at the right-hand side of Equation (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.

The Fourier coefficients of f (j)
L are defined by a j k and approximated by c j k , where We also write a k and c k instead of a 0 k and c 0 k , respectively.Intuitively, we then have where ′ indicates that the first summand (with k = 0) is weighted by one-half.A little analysis shows that , k = 0, 1, ..., j = 0, ..., J + 1, i.e., the coefficients c j k can be obtained explicitly if ϕ is given in closed form.Here, ℜ(z) denotes the real part of a complex number z and i the imaginary unit.For 0 < M ≤ L define To keep the notation simple, we suppress the dependence of a j k and c j k 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 The coefficients c k are given in closed form when ϕ 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 Fang and Oosterlee [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 (2.7) Choose N large enough so that (2.9) Proof.Junike and Pankrashkin [21,Cor. 8].
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.

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 J+1 b (R) the set of bounded functions from R to R which are (J + 1)-times, continuously differentiable with bounded derivatives.By .∞ and . 2 we denote the supremum norm and the L 2 norm, i.e., We say that f and its derivatives have semi-heavy tails if We say that f and its derivatives have heavy tails with index α > 0 if We suppose that f satisfies one of the following assumptions: 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 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 α 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.
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π].Here, we explicitly need the dependence of the series-truncation error on the truncation range [−L, L], so we give the full proof.
Given L > 0, we need to find an upper bound for f (j) ∞ to estimate the series truncation error by Inequality (3.4).It follows by the inverse Fourier transform and Equation (2.3) that Inequality (3.8) provides an explicit expression to find a bound for the term f (j) ∞ , j = 0, .., J + 1 for several models.
Example 3.6.In the symmetric NIG model with parameters α > 0, β = 0 and δ > 0, the centralized log-returns at time T > 0 have density f NIG ∈ C ∞ b (R), which can be expressed in terms of the modified Bessel-function of the third kind.The characteristic function is given by [5] and Schoutens [41,Sec. 5.3.8].We obtain, by Inequality (3.8) and using 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 ∨ y.
Note that for all k ∈ N, Similarly, we arrive at 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 ε > 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 Section 5. To approximate the time-0 price of the option by Theorem 3.8, set ℓ = 0. Theorem 3.8.(Find M , L and N ).Let v : R → R be bounded, with |v(x)| ≤ K for all x ∈ R and some K > 0. Let ε > 0 be small enough.Suppose J ∈ N 0 .i) Assume density f satisfies Assumption A1.For some even n ∈ N let . (3.13) Let ℓ = 0, k ∈ {1, ..., J} and define N as in Equation (3.11).
In both cases, i) and ii), it holds that Proof.We start with case i).For ε small enough, M is large enough so that Assumption A1 holds, i.e., L 0 ≤ M = L.It follows that the last inequality holds true if Further, the last inequality holds true if Proposition 2 in Junike and Pankrashkin [21] shows that We then have the last inequality holds true if Assume J ≥ 1.For ε small enough, we have because the right-hand side of Inequality (3.23) is of order ε − 1 n while the right-hand side of Inequality (3.11) the last inequality holds true if . Next, we treat the case ii).For ε small enough, M is large enough so that Assumption A2 holds, i.e., L 0 ≤ M ≤ L. The inequality is satisfied by the definition of M .A little calculation shows that the definition of L implies 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, Using the bound for B f (L) from Lemma 3.7, we obtain the last Inequality holds by the definition of L. For ε small enough, N is large enough and we have Use Inequality (3.28), the definition of L and 96 π 2 ≤ 12, to see that . (3.29) the last inequality holds by Inequality (3.29).By the definition of N , we also have .
By Lemma 3.5 it follows that As B f (L) → 0, L → ∞, f is COS-admissible.Apply Lemma 2.1 to conclude.
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 ε.Numerical experiments in Junike and Pankrashkin [21] suggest choosing n ∈ {4, 6, 8} for µ n .According to Theorem 3.8, any k ∈ {1, ..., J − ℓ} 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 Section 6.1.For other models, another choice for k might be more suitable.Bounds for f (k+1) ∞ 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 α in Assumption A2 to estimate M , L or N .The constants C 3 and α are known, for example, for the FMLS model.

Order of Convergence
Theorem 4.1 describes the order of convergence of the COS method if we allow N → ∞ 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., Let M = L = L(N ).We say the error of the COS method converges with order ρ > 0, if there is a constant κ > 0 such that for all N ∈ N it holds that The error is of infinite order or converges exponentially, if Inequality (4.1) holds for any ρ, see Boyd (ii) Assume density f is unimodal and satisfies Assumption A2.Let β ∈ 0, J J+2+2α and γ > 0.

In both cases we have err(L(N
As in the proof of Corollary 8 in Junike and Pankrashkin [21], one can show that where B f (L) as in Equation (2.6) and 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 We now prove (i).For L large enough, we have Further, by Inequality (3.17), it holds that Assuming L ≥ 1 and applying Inequality (3.20), it follows that

Inequality (3.4) implies
where we used L N ≤ γ.Let b 1 , ..., b 4 > 0 be suitable constants.By Inequality (4.2), it follows for N large enough that < 0 and the right-hand side of (4.3) converges to zero for N → ∞.We show (ii).It holds for L large enough using Inequalities (3.26) and (3.27) that Inequality (3.4) and Assumption A2 imply for N large enough and for some suitable constants a 1 , ..., a J > 0 that By Inequality (4.2), it holds for some suitable constants b The last inequality can be seen as follows: as β ≤ J J+2+α , it follows that −αβ ≥ β(J + 2) − J.
Remark 4.2.The COS method converges exponentially, i.e., faster than O(N −ρ ) as N → ∞ for any ρ > 0 if f is smooth and has semi-heavy tails.To see this, let 0 < β < 1: then, 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, where ŝ is a spectral filter, i.e., a smooth function with support [−1, 1] and ŝ(0) = 1.For an analysis of the order of convergence and some error bounds for the filter-COS method, see Ruijter et al. [38] and references therein.

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 St for some stochastic process ( St ) t≥0 , which does not depend on S 0 anymore.This assumption is a very typical one, see Madan and Schoutens [31,Sec. 3.1.2].As in Section 2, we consider a European option with maturity T > 0 and payoff w(S T ) for some w : [0, ∞) → R. Let η := E[log(S T )] − log(S 0 ).Then η does not depend on S 0 .Define v(x, s) := e −rT w(exp(x + log(s) + η)), x ∈ R, s > 0. (5.1) The time-0 price of the European option is then given by where, as before, f is the density of the centralized log-returns.Delta and Gamma are defined by if the partial derivatives exist.The next lemma provides some conditions to interchange integration and differentiation in Equation (5.2).
Lemma 5.1.Let w be bounded.Let v be defined as in Equation (5.1).Assume J ∈ N 0 and density f satisfies Assumption A1.It follows that Then x → h(x, s) is integrable for all s ∈ (s, s) and the partial derivative for all (x, s) ∈ R × (s, s) and x → m(x) is integrable.Interchanging differentiation and integration is allowed by the dominated convergence theorem, see e.g., Grubb [17,Lemma 2.8], and it follows that If J ≥ 1, f is twice differentiable.Apply the arguments above to f (1) to conclude.
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 ε > 0. One can use the same parameters M , N and L to obtain both the price and the Greeks.We define v k := M −M v(x, S 0 )e L k (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 ε > 0 be small enough.Let γ = min ε, εS 0 , εS 2 0 2 .Suppose J ≥ 3. Assume density f satisfies Assumption A1.For some even n ∈ N define where µ n is the n th −moment of f , i.e., µ n = 1

.4)
It follows that the time-0 price and the Greeks Delta and Gamma can be approximated by the COS method, i.e., 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 Using the triangle inequality, we can see that

Numerical experiments
Some numerical experiments are compared with the Carr-Madan formula, see Carr and Madan [8], which is applicable when the characteristic function of the log-returns is given in closed form and when E[S 1+γ T ] is finite for some γ > 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 γ = 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.

BS model
We consider the BS model with volatility σ > 0 and maturity T > 0. The density f BS of the log-returns is normally distributed and belongs to the family of stable laws.An upper bound for f (k+1) BS ∞ can be obtained directly from Inequality (6.3) setting α = 2 and c = σ √ T √ 2 .Let n ∈ N be even.In the BS model, the n th moment of the centralized log-returns is given by The formula for N in Inequality (3.11) does not depend on the volatility σ √ T if we set ℓ = 0 and if we bound f (k+1) BS ∞ using Inequality (6.3).In Figure 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 ≥ 40.We also see that N increases as the error tolerance ε decreases or the bound K increases.
How sharp is the bound for N in Theorem 3.8 and Theorem 5.2?We make the following experiment.Consider a put option with parameters 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 f (k+1) ∞ .Here, we use R's integrate function with default values.The CPU time of the COS method and the numerical integration routine to obtain a bound for f (k+1) ∞ are of similar magnitude; see Table 1.
The value of N does not change when using numerical integration to obtain a bound for f (k+1) ∞ compared to the closed form solution for the bound of f (k+1) ∞ .1: Approximation of the time-0 price of a put option by the COS method: N depending on k by Inequality (3.11).CPU time is measured in milliseconds.The COS method with N min = 120 takes about 0.068 milliseconds.We set n = 8 and ℓ = 0.In the BS model, N does not depend on the standard deviation of the log-returns, i.e., on the value of σ √ T .

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 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 |.Aimi et al. [1,Sec. 3.2] propose to select the number of terms for the BS model with volatility σ > 0 and maturity T by the smallest natural number N F satisfying the following inequality where [a, b] = [−L, L] is the truncation range and ε NF is some error tolerance.The solution to find the number of terms by (6.2) works even better than Inequality (3.11) if ε NF is small enough, see Table 2.However, the rule by Inequality (6.

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 Dragulescu and Yakovenko [12].In this section, we provide numerical experiments under two different parameter sets, M1 and M2.The parameters of models M1 and M2 are taken from Fang and Oosterlee [14] and Schoutens et al. [42], respectively.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 ε.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 f (k+1) ∞ 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.the COS method depending on the strike K, the maturity T and the parameters of the model.We set ε = 10 −3 , n = 4, k = 20, S 0 = 100 and r = 0 to obtain N by Inequality (3.11).CPU time is measured in milliseconds.The CPU time to estimate f (k+1) ∞ by numeric integration is also provided in the last row.

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 Albin and Sundén [2,Example 7.5].It can be expressed by means of the Whittaker function and is (J + 1)-times continuously differentiable if J + 2 < 2T ν , see Küchler and Tappe [22].Let f VG denote the density of the log-returns in the VG model at maturity T > 0 with parameters σ > 0, θ ∈ R and ν > 0. If T < ν 2 , the density f VG is unbounded and Theorem 3.8 cannot be used to find N .If T ∈ ν, 3ν  2 , the derivative of f VG is continuous, but the second derivative of f VG is unbounded, see Küchler and Tappe [22,Thm. 4.1 and Thm. 6.1].
We can apply Theorem 3.8 with J = 0 if T ∈ ν, 3ν 2 , 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: We calculated a reference price π 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.

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 Stable ∈ C ∞ b (R) of the family of stable distributions with stability parameter α ∈ (0, 2], skewness β ∈ [−1, 1], scale c > 0 and location θ ∈ R has the characteristic function see Zolotarev [46].It follows by Inequality (3.8) that α παc j+1 , j = 0, 1, 2, ..., (6.3)where Γ denotes the gamma function.The density of the stable law is unimodal, see Yamazato [44]. where For stable densities we therefore suggest to set C 3 in Theorem 3.8 at least as large as αC α The density of the log-returns in the FMLS model has a heavy left tail with Pareto index α, i.e., the left tail decays like |x| −1−α , x → −∞, but the right tail decays exponentially, see Carr and Wu [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 α = 1.5597 and σ = 0.1486.We test the formulas for M , L and N with these values for α and σ for a European call option with the following parameters: 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 → C 3 |x| −1−α .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 Carr and Wu [9].With these parameters, the COS method prices the call option within the error tolerance in about 1.5 milliseconds.The minimal N to stay below the error tolerance is N min = 1200 and the CPU time using N min is about 0.4 milliseconds.
We also apply the Carr-Madan formula with the "default parameters", see Carr and Madan [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 milliseconds.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.

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 Leisen and Reimer [24].As log κ N ρ = log(κ) − ρ log(N ), the negative slope of a straight line obtained from a log-log plot of the error err(L(N ), N ) against N can be used as an indicator for ρ.As in Section 4, we assume M = L in this section.

BS Model
In the BS model with parameter σ > 0, the density of the log-returns is arbitrarily smooth, and the tails decay even faster than exponentially.Figure 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 σ = 0.2, S 0 = K = 100, r = 0 and T = 1.We see in Figure 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σ or L = 6σ, the COS method does not converge for N → ∞ 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 → ∞ 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.For large L, such as L = 20σ, this effect disappears somewhat because the tails of the Gaussian distribution decay so rapidly that, up to fixed-precision arithmetic1 , 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σ, the error of the COS method does not converge to zero for N → ∞.
If we choose L = L(N ) = σ √ N , Theorem 4.1(i) indicates that the error of the COS method converges exponentially to zero.This is confirmed empirically in Figure 3.Other choices for L also work well, e.g., L = σ 5 N .

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 α = 1.5597 and σ = 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 Section 6.4.We compute L optimal for a fixed N ∈ N minimizing 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 ∈ {2 m , m = 4, 5, ..., 20} we define L optimal (N ) := argmin L∈L err(L, N ), where L = {exp(0.07m),m = 0, 1, ..., 200}, is a sufficiently fine grid of the interval [1, 10 6 ].
In Figure 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 = 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 optimal instead of L = 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 Figure 4 we also plot the empirical order of convergence of the Carr-Madan formula for the FMLS model with damping factor of 0.

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 −α ), where α > 0 is the Pareto tail index.This is the case, for example, for the FMLS model where α ∈ (1, 2).Numerical experiments indicate that the bound O(N −α ) is sharp.

[ 7 ,Theorem 4 . 1 .
Sec. 2.3].We use big O notation: for functions g, h : N → R, we write h(N ) = O(g(N )) as N → ∞, if the absolute value of h(N ) is at most a positive constant multiple of g(N ) for all sufficiently large values of N .(Bounds for the order of convergence).Let v : R → R be bounded, with |v(x)| ≤ K for all x ∈ R and some K > 0. Assume J ∈ N. (i) Assume density f satisfies Assumption A1.Let β ∈ 0, J J+3 and γ > 0. If M = L = γN β it holds that err(L(N ), N ) ≤ O N −J(1−β)+2β , as N → ∞.

Remark 4 . 3 .
By Theorem 4.1, the COS method converges at least like O(N −α ) as N → ∞ if f is smooth and has heavy tails with Pareto index α > 0. Numerical experiments indicate that the COS method does not converge faster than O(N −α ) for the FMLS model, see Section 6.5.2, i.e., the bound in Theorem 4.1(ii) is sharp.

Figure 1 :
Figure 1: N depending on k by Inequality (3.11) for various ε and K.We set n = 8 and ℓ = 0.In the BS model, N does not depend on the standard deviation of the log-returns, i.e., on the value of σ √ T .

Figure 2 :
Figure 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 milliseconds.The minimal N to stay below the error tolerance is N min = 1200 and the CPU time using N min is about 0.4 milliseconds.We also apply the Carr-Madan formula with the "default parameters", see Carr and Madan[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 milliseconds.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.

Figure 3 :
Figure 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].

Figure 4 :
Figure 4: Order of convergence of the COS method with different choices for the truncation range [−L, L] for the FMLS model.
7 and a Fourier truncation range of 1200.The Figure indicates an exponential order of convergence for the Carr-Madan formula.
[21]called COS-admissible.The class of COSadmissible densities is very large; in particular, it includes bounded densities with existing first and second moments and stable densities, see Junike and Pankrashkin[21].
and f and its derivatives have semi-heavy tails.
Remark 3.2.In exceptional cases, the constants C 1 and C 2 are explicitly known; see Example 3.4.

Table 3 :
Number of terms N and CPU time to approximate the time-0 price of a put option by Let F Stable be the cumulative distribution function for the stable density f Stable .By Samorodnitsky and Taqqu [40, Property 1.2.15] it holds that lim x→∞