Draft:The COS method

The COS method is a numerical technique in computational finance to efficiently price European plain vanilla put and call options provided the characteristic function of the logarithmic returns is given in closed-form. The main idea of the COS method is to approximate the (usually unknown) probability density function of the logarithmic returns by a cosine series expansion. The method has been developed in 2008 by Fang Fang and Cornelis W. Oosterlee. The COS method converges for a large family of densities. The COS method converges exponentially if the density of the logarithmic returns is smooth and has semi-heavy tails.

Introduction

Consider a financial market with a bank-account paying a risk-free rate ${\textstyle r}$ and stock with deterministic price ${\textstyle S_{0}>0}$ today and random price ${\textstyle S_{T}}$ at some future date ${\textstyle T>0}$. Let ${\textstyle \mu }$ be the expectation of ${\textstyle \log(S_{T})}$ under the risk-neutral measure and assume the characteristic function ${\textstyle \varphi }$ of the centralized log-returns ${\textstyle X:=\log(S_{T})-\mu }$ is given in closed-form. The characteristic function ${\textstyle \varphi }$ is given explicitly for many models and processes, e.g., the Black-Scholes model, the variance gamma process, the CGMY model or the Heston model. The price of a European put option with maturity ${\textstyle T>0}$ and strike ${\textstyle K>0}$ is given by the fundamental theorem of asset pricing by

$$\int _{\mathbb {R} }e^{-rT}\max \left(K-e^{x+\mu },0\right)f(x)dx,$$

where ${\textstyle f}$ is the probability density function of ${\textstyle X}$. The price of a call option can be obtained by the put-call parity. Very often, ${\textstyle f}$ is not given explicitly and the COS method is used to approximate ${\textstyle f}$ and the price of the option.

The method

For some ${\textstyle L>0}$, the density ${\textstyle f}$ is truncated and the truncated density is approximated by a cosine series expansion:

$$f(x)\approx {\frac {c_{0}}{2}}+\sum _{k=1}^{N}c_{k}\cos \left(k\pi {\frac {x+L}{2L}}\right),\quad x\in [-L,L],$$

where for ${\textstyle k=0,1,...,N}$ the coefficients ${\textstyle c_{k}}$ are defined by

$${\begin{aligned}c_{k}:={\frac {1}{L}}\int _{\mathbb {R} }f(x)\cos {\big (}k\pi {\frac {x+L}{2L}}{\big )}dx=&{\frac {1}{L}}\Re {\bigg \{}\varphi {\bigg (}{\frac {k\pi }{2L}}{\bigg )}e^{i{\frac {k\pi }{2}}}{\bigg \}}.\end{aligned}}$$

${\textstyle \Re (z)}$ denotes the real part of a complex number ${\textstyle z}$ and ${\textstyle i={\sqrt {-1}}}$ is the imaginary unit. The price of a European put option can be approximated by

$$\int _{\mathbb {R} }e^{-rT}\max \left(K-e^{x+\mu },0\right)f(x)dx\approx {\frac {c_{0}v_{0}}{2}}+\sum _{k=1}^{N}c_{k}v_{k},$$

where

$${\begin{aligned}v_{k}&:=\int _{-L}^{L}e^{-rT}\max \left(K-e^{x+\mu },0\right)\cos {\big (}k\pi {\frac {x+L}{2L}}{\big )}dx.\end{aligned}}$$

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

Let ${\textstyle d:=\min \left(\log \left(K\right)-\mu ,L\right)}$. For a European put option it holds that ${\textstyle v_{k}=0}$ if ${\textstyle d\leq -L}$ and otherwise

$$v_{k}=e^{-rT}\left(K\Psi _{0}(k)-e^{\mu }\Psi _{1}(k)\right)$$

where

$$\Psi _{0}(k)={\begin{cases}d+L&,k=0\\{\frac {2L}{k\pi }}\sin {\big (}k\pi {\dfrac {d+L}{2L}}{\big )}&,k>0\end{cases}}$$

and

$${\begin{aligned}\Psi _{1}(k)&={\frac {e^{d}\left({\frac {k\pi }{2L}}\sin \left(k\pi {\frac {d+L}{2L}}\right)+\cos \left(k\pi {\frac {d+L}{2L}}\right)\right)-e^{-L}}{1+\left({\frac {k\pi }{2L}}\right)^{2}}},\quad k\geq 0.\end{aligned}}$$

To price a call option it is numerically more stable to price a put option instead and use the put-call parity.

Parameters

To apply the COS method, one has to specify the truncation range ${\textstyle [-L,L]}$ and the number of terms ${\textstyle N}$. Let ${\textstyle \varepsilon }$ be the error tolerance between the true price and its approximation by the COS method. If ${\textstyle \varepsilon }$ is small enough and ${\textstyle f}$ has semi-heavy tails, the truncation range of a put option can be chosen using Markov’s inequality by

$$L=\left({\frac {2Ke^{-rT}\mu _{n}}{\varepsilon }}\right)^{\frac {1}{n}},$$

where ${\textstyle n\in \mathbb {N} }$ is even and ${\textstyle \mu _{n}}$ is the ${\textstyle n}$-th moment of ${\textstyle X}$, which can be obtained using a computer algebra system and the following relation of ${\textstyle \mu _{n}}$ and the characteristic function:

$$\mu _{n}={\frac {1}{i^{n}}}\left.{\frac {\partial ^{n}}{\partial u^{n}}}\varphi (u)\right|_{u=0}.$$

Often, ${\textstyle n\in \{4,6,8\}}$ is a reasonable choice. In particular, the fourth moment can also be obtained by the Kurtosis and variance of ${\textstyle X}$ by ${\textstyle \mu _{4}={\text{Kurt}}(X)\left({\text{Var}}(X)\right)^{2}}$.

If ${\textstyle f}$ is also a ${\textstyle s+1\in \mathbb {N} }$ times differentiable function with bounded derivatives, the number of terms can then be chosen using integration by parts by

$$N\geq \left({\frac {2^{s+{\frac {5}{2}}}\left\Vert f^{(s+1)}\right\Vert _{\infty }L^{s+2}}{s\pi ^{s+1}}}{\frac {12Ke^{-rT}}{\varepsilon }}\right)^{\frac {1}{s}}.$$

${\textstyle \|.\|_{\infty }}$ denotes the uniform norm. There is the following useful inequality:

$$\|f^{(s+1)}\|_{\infty }\leq {\frac {1}{2\pi }}\int _{\mathbb {R} }|u|^{s+1}|\varphi (u)|du.$$

The last integral can be solved numerically with standard techniques and it is also given explicitly in some cases.

Black-Scholes model

In the Black-Scholes model with volatility ${\textstyle \sigma >0}$, the stock price is modelled at maturity ${\textstyle T>0}$ by

$$S_{T}=S_{0}e^{rT-{\frac {1}{2}}\sigma ^{2}T+\sigma {\sqrt {T}}Z},$$

where ${\textstyle Z}$ is a standard normal random variable. The expectation of ${\textstyle \log(S_{T})}$ under the risk-neutral measure is given by

$$\mu =\log(S_{0})+rT-{\frac {1}{2}}\sigma ^{2}T$$

and the characteristic function of ${\textstyle X=\log(S_{T})-\mu }$ is given by

$$\varphi (u)=e^{-{\frac {\sigma ^{2}T}{2}}u^{2}},\quad u\in \mathbb {R} ,$$

since ${\textstyle X}$ follows a normal distribution. The ${\textstyle n}$-th moment of ${\textstyle X}$ for even ${\textstyle n}$ in the Black-Scholes model is given by

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

The derivative ${\textstyle f^{(s+1)}}$ of the probability density function of ${\textstyle X}$ is bounded by

$$\|f^{(s+1)}\|_{\infty }\leq {\frac {\Gamma \left({\frac {s}{2}}+1\right)2^{\frac {s}{2}}}{\pi (\sigma {\sqrt {T}})^{s+2}}},$$

where ${\textstyle \Gamma }$ is the Gamma function.

Simple implementation

In the following code, the COS method is implemented in R for the Black-Scholes model to price European call option and the Heston model to price a put option.

########### Black-Scholes model ###########
#Characteristic function of the centralized log-returns in the BS model. params is volatility
phi = function(u, mat, params, S0, r){
  return(exp(-params^2 * mat * u^2 / 2))
}

#mu is equal to E[log(S_T)]
mu = function(mat, params, S0, r){
    return(log(S0) + r * mat - 1/2 * params^2 * mat)
}   

#cosine coefficients of the density.
ck = function(L, mat, N, params, S0, r){
  k = 0:N
  return(1 / L * Re(phi(k * pi / (2 * L), mat, params, S0, r) *  exp(1i * k * pi/2)))
}

#cosine coefficients of a put option, see Appendix Junike and Pankrashkin (2022).
vk = function(K, L, mat, N, params, S0, r){
  mymu = mu(mat, params, S0, r)   #mu = E[log(S_T)]
  d = min(log(K) - mymu, L)
  if(d <= -L)
    return(rep(0, N + 1)) #Return zero vector
  k = 0:N 
  psi0 = 2 * L / (k * pi) * (sin(k * pi * (d + L) / (2 * L)))
  psi0[1] = d + L
  tmp1 = k * pi / (2 * L) * sin( k * pi * (d + L) / (2 * L))
  tmp2 = cos(k * pi * (d + L) / (2 * L))
  tmp3 = 1 + (k * pi / (2 * L))^2
  psi1 = (exp(d) * (tmp1 + tmp2) - exp(-L)) / tmp3
  return(exp(-r * mat) * (K * psi0 - exp(mymu) * psi1))
}

#approximation of put option by COS method
put_COS = function(K, L, mat, N, params, S0, r){
  tmp = ck(L, mat, N, params, S0, r) * vk(K, L, mat, N, params, S0, r)
  tmp[1] = 0.5 * tmp[1] #First term is weighted by 1/2
  return(sum(tmp))
}

#approximation of call option by COS method using put-call parity
call_COS = function(K, L, mat, N, params, S0, r){
  return(put_COS(K, L, mat, N, params, S0, r) + S0 - K * exp(-r * mat))
}

#Price a call option in BS model by the COS method
eps = 10^-4 #error tolerance
K = 90 #strike
S0 = 100 #current stock price
r = 0.1 #interest rates
params = 0.2 #volatility
mat = 0.7 #maturity
n = 8 #number of moments to get truncation range
mu_n = (params * sqrt(mat))^n * prod(seq(1, n, by = 2)) #n-th moment of log-returns
L = (2 * K * exp(-r * mat) * mu_n / eps)^(1 / n) #Truncation range, Junike (2024, Eq. (3.10))
s = 40 #number of derivatives to determine the number of terms
boundDeriv = gamma(s / 2 + 1) * 2^(s / 2) / (pi * (params * sqrt(mat))^(s + 2))
tmp = 2^(s + 5 / 2) * boundDeriv * L^(s + 2) * 12 * K * exp(-r * mat)
N = ceiling((tmp / (s * pi^(s + 1) * eps))^(1 / s)) #Number of terms, Junike (2024, Sec. 6.1)
call_COS(K, L, mat, N, params, S0, r) #The price of call option is 17.24655.

########### Heston model ###########
#Use functions vk, put_COS, call_COS from the Black-Scholes model above.

#characteristic function of log-returns in the Heston with parameters params.
#The characteristic function is taken from Schoutens et. al (2004).
psiLogST_Heston = function(u, mat, params, S0, r){
  kappa = params[1] #speed of mean reversion
  theta = params[2] #level of mean reversion
  xi = params[3] #vol of vol
  rho = params[4] #correlation vol stock
  v0 = params[5] #initial vol
  d = sqrt((rho * xi * u * 1i - kappa)^2 - xi^2 * (-1i * u - u^2))
  mytmp = kappa - rho * xi * u * 1i
  g = (mytmp - d) / (mytmp + d)
  expdmat = exp(-d * mat)
  tmp0 = 1i * u * (log(S0) + r * mat)
  tmp1 = (mytmp - d) * mat - 2 * log((1 - g * expdmat) / (1 - g))
  tmp2 = theta * kappa * xi^(-2) * tmp1
  tmp3 = v0 * xi^(-2) * (mytmp - d) * (1 - expdmat) / (1 - g * expdmat)
  exp(tmp0 + tmp2 + tmp3)
}

library(Deriv) #There are much faster alternatives like SageMath.
psiLogST_Heston1=Deriv(psiLogST_Heston, "u") 

#mu is equal to E[log(S_T)]
mu = function(mat, params, S0, r){
    Re(-1i * psiLogST_Heston1(0, mat, params, S0, r))
}
#Characteristic function of centralized log-returns in the Heston model.
phi = function(u, mat, params, S0, r){
    psiLogST_Heston(u, mat, params, S0, r) * exp(-1i * u * mu(mat, params, S0, r))
}

#Derivatives of the characteristic function of the centralized log-returns in the Heston model.
phi1 = Deriv(phi, "u")
phi2 = Deriv(phi1, "u")
phi3 = Deriv(phi2, "u") #Takes very long but has to be done only once.
phi4 = Deriv(phi3, "u") #Takes very long but has to be done only once.
save(phi4, file = "phi4.RData") #save for later use. Load with load("phi4.RData").

#Price a put option in Heston model by the COS method.
eps = 10^-6 #error tolerance
K = 90 #strike
S0 = 100 #current stock price
r = 0.1 #interest rates
params = c(0.6067, 0.0707, 0.2928, -0.7571, 0.0654)
mat = 0.7 #maturity
mu_n = abs(phi4(0, mat, params, S0, r)) #4-th moment of log-returns. 
L = (2 * K * exp(-r * mat) * mu_n / eps)^(1 / 4) #Truncation range, Junike (2024, Eq. (3.10)). It is 36.99.
s = 20 #number of derivatives to determine the number of terms
integrand = function(u){1 / (2 * pi) * abs(u)^(s + 1) * abs(phi(u, mat, params, S0, r))} 
boundDeriv = integrate(integrand, -Inf, Inf)$value
tmp = 2^(s + 5 / 2) * boundDeriv * L^(s + 2) * 12 * K * exp(-r * mat)
N = ceiling((tmp / (s * pi^(s + 1) * eps))^(1 / s)) #Number of terms, Junike (2024, Sec. 6.1). It is 4418.   
put_COS(K, L, mat, N, params, S0, r) #The price of put option is 2.773954.

See also

• Carr–Madan formula