Finance
M.P. Scholl — Wednesday 03 February 2016

The COS Method¶

In this post we implement the method proposed in

A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions by F. Fang and C. W. Oosterlee.

SIAM Journal on Scientific Computing archive

Volume 31 Issue 2, November 2008

Pages 826-848

Density Function¶

We want to derive the characteristic function φ(u) of the standard normal distribution $\mathbb{N}(0, 1)$ by solving:

$$ \phi(u) = \int_{-\infty}^{+\infty} e^{iux} f(x) dx $$

For the normal distr

$$ \mathcal{F}(u) = \int_{-\infty}^{+\infty} \frac{1}{\sqrt{2\pi}}e^{\frac{x^2}{2}}e^{iux}dx $$

Separating the integral

$$ \mathcal{F}(u) = \int_{0}^{+\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}e^{iux}dx + \int_{-\infty}^{0} \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}e^{iux}dx $$

We negate the right side so the domains of the integral match

$$ \mathcal{F}(u) = \int_{0}^{+\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}e^{iux}dx + \int_{0}^{+\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}e^{-iux}dx $$

Eulers identity

$$ \mathcal{F}(u) = \int_{0}^{+\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} [\cos(ux) + i \sin(ux)] dx + \int_{0}^{+\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} [\cos(-ux) + i \sin(-ux)]dx $$

Since $\cos$ is even: $\cos(-ux)=\cos(ux)$, and $\sin$ is odd: $\sin(-ux)=-\sin(ux)$. Simplifying:

\begin{equation} \mathcal{F}(u) = 2\int_{0}^{+\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} \cos(ux)dx \end{equation}

Now $\mathcal{F}^\prime(u)$

$$ \mathcal{F}^\prime(u) = -2\int_{0}^{+\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} \sin(ux)dx $$$$ \mathcal{F}^\prime(u) = \frac{2}{\sqrt{2\pi}}\int_{0}^{+\infty} \sin(ux) d[e^{-\frac{x^2}{2}}] $$

Integrating by parts $\int u dv = uv - \int v du$.

$$ \mathcal{F}^\prime(u) = \frac{2}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\sin(ux)|^{+\infty}_0 - \frac{2}{\sqrt{2\pi}}\int_{0}^{+\infty} e^{-\frac{x^2}{2}} u\cos(ux)dx $$$$ \mathcal{F}^\prime(u) = -u 2\int_{0}^{+\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} u\cos(ux)dx $$

Using 1 we see that

\begin{equation} \mathcal{F}^\prime(u) = -u\mathcal{F}(u) \end{equation}

as required.

Equation \ref{ode} has a simple solution:

$$\phi(u) = e^{-\frac{1}{2}u^2}$$

Now the coefficients $F_k$

\begin{equation} F_k = \frac{2}{b-a}\int_a^b f(x)cos(k\pi \frac{x-a}{b-a}) \end{equation}

This leads us to the approximation (taking the finite interval a,b and taking the real part as stipulated in the assignment)

$$ F_k = \frac{2}{b-a} Re[\phi(\frac{k\pi}{b-a})e^{-i\frac{ka\pi}{b-a}}] $$

For our final function this is:

$$ f(x) \approx \frac{1}{2}F_0 + \sum_{k=1}^N F_k cos(k\pi\frac{x-a}{b-a}) $$

In the paper the notation $\sum^\prime$ is used to denote a sum where the first element is counted half. This works since at $k=0$ the $\cos$ term is 1.

\begin{equation} f(x) \approx \sum^N_{k=0}\prime F_k cos(k\pi\frac{x-a}{b-a}) \end{equation}

Convergence¶

Since we are dealing with a sum we will investigate the elements in the sum.

Lets investigate the $k$-th element in equation \ref{normalsum}. We quickly find a constant bound on the right hand term (is this the right notation?)

$$|\cos(k\pi\frac{x-a}{b-a})|=1$$

in essence $\mathcal{O}(1)$

Now the bound on the $k$-th term entirely dependent on $F_k$. Herein we have two terms dependent on $k$:

$$ \phi(\frac{k\pi}{b-a}) = \mathcal{O}(e^{-k^2}) $$

and by Eulers identity and taking the real part we also quickly find that

$$ e^{-i\frac{ka\pi}{b-a}} = \mathcal{O}(1) $$

Presuming that the method converges, we now say that increasing $N$ only adds new elements of the orders described above, and the second (non-constant) term is exponential, we can say the method has exponential convergence.

Option Price¶

$$ f(x)\approx e^{-rT}\frac{b-a}{2}\sum_{k=0}^N F_k(x)G_k $$$$ V(x,t) = e^{-rT} \int_\mathcal{R}g(y)f(y|x)dy $$

We truncate the range to a,b:

\begin{equation} v(x,t) = e^{-rT} \int_a^b g(y)f(y|x)dy \end{equation}

Firstly we find in the paper the risk neutral valuation

$$ v(x,t_0)= = e^{-r\Delta t}\mathbb{E}^\mathbb{Q}[v(y,T)|x] = e^{-rT} \int_\mathcal{R}v(y,T)f(y|x)dy $$

Only the characteristic function of $f$ is known, so we replace it by the cosine expansion

$$ f(y|x)=\sum_{k=0}^{+\infty}\prime F_k(x) \cos(k\pi\frac{x-a}{b-a}) $$

Taking equation \ref{coeff} as an example, we can express $A_k$ as

$$ F_k(x)= \frac{2}{b-a} \int_a^b f(y|x)\cos(k\pi\frac{x-a}{b-a})dy $$

We insert this into equation \ref{eqv} obtaining:

$$ v(x,t)=e^{-rT} \int_a^b v(y,T) \sum_{k=0}^NF_k(x) \cos(k\pi\frac{x-a}{b-a}) $$

Let $$ G_k = \frac{2}{b-a}\int_a^b v(y,T) \cos(k\pi\frac{x-a}{b-a})dy $$

We arrive at the required equation:

$$ f(x)\approx e^{-rT}\frac{b-a}{2}\sum_{k=0}^N F_k(x)G_k $$

In [1]:

import math
import scipy.stats as stats
import numpy as np

import matplotlib.pyplot as pyplot
pyplot.style.use('ggplot')
%matplotlib inline
from pylab import rcParams
rcParams['figure.figsize'] = 16.5, 10.5

In [2]:

def Phi(u):
    return float(np.exp(-0.5 * u * u))

def F_k(k,a,b):
    x = (k * a * np.pi) / (b - a)
    phi = Phi((k * np.pi) / (b -a))
    exponent = np.exp(complex(0, -x))
    return (2. / (b - a)) * (phi * exponent).real

def f(x, N, a, b):
    accumulator = 0.5 * F_k(0, a, b)
    for k in np.arange(1,N):
        accumulator += F_k(k, a, b) * np.cos((k * np.pi) * ((x - a) / (b - a)))
    return accumulator

def Normal(x):
    return 1. / np.sqrt(2 * np.pi) * np.exp(-0.5 * x * x)

In [3]:

fig = pyplot.figure()
fig.set_size_inches(16.5,10.5)
pw = [2**i for i in range(1,7)]
xs = np.arange(-5,5,0.05)
legends = []
for N in pw:
    a = -10
    b = 10
    ys = [f(x,N,a,b) for x in xs]
    pyplot.plot(xs, ys)
    legends.append("N="+str(N))
    
yt =[Normal(x) for x in xs]
legends.append("Reference")
pyplot.xlabel("x")
pyplot.ylabel("f(x)")
pyplot.plot(xs, yt)
pyplot.legend(legends)

pyplot.show()

Notation¶

$$ x = \log(\frac{S_0}{K}) $$$$ u = \frac{k\pi}{b-a} $$

Chi¶

$$ \chi_k = \frac{\cos(k\pi\frac{d-a}{b-a})e^d-\cos(k\pi\frac{c-a}{b-a})e^c + \frac{k\pi}{b-a}[\sin(k\pi\frac{d-a}{b-a})e^d-\sin(k\pi\frac{c-a}{b-a})e^c]}{1+(\frac{k\pi}{b-a})^2} $$

In [4]:

def Chi(k, a, b, c, d):
    u = k * (np.pi / (b - a))
    expd = np.exp(d)
    expc = np.exp(c)
    cma = (c - a ) * u
    dma = (d - a ) * u
    return (np.cos(dma) * expd 
            - np.cos(cma) * expc 
            + u * np.sin(dma) * expd 
            - u * np.sin(cma) * expc) / (1 + u * u)

Psi¶

$$ \psi_k = \left\{ \begin{array}{lr} d - c & : 0 = k\\ [\sin(k\pi\frac{d-a}{b-a}) - \sin(k\pi\frac{c-a}{b-a})]\frac{b-a}{k\pi} & : 0 < k \end{array} \right. $$

In [5]:

def Psi(k, a, b, c, d):
    if 0 == k:
        return d - c
    u = k * (np.pi / (b - a))
    cma = (c - a ) * u
    dma = (d - a ) * u
    return (np.sin(dma) - np.sin(cma)) / u

$$ \phi(u) = e^{iu(r-\frac{1}{2}\sigma^2)\tau - \frac{1}{2}\sigma^2\tau u^2} $$

In [6]:

def Phi(u, r, tau, sigma):
    return np.exp(complex(-0.5 * sigma * sigma * tau * u * u, u * (r - 0.5 * sigma * sigma) * tau))

Truncation range¶

N.b. that $a$ is less than $b$, the minus sign applies to the $a$ in the minusplus sign $\mp$

First¶

$$ a,b = \log(\frac{S_0}{K}) + r\tau \mp L\sqrt{\tau\sigma^2} $$

Second¶

$$ a,b = \log(\frac{S_0}{K}) + (r-\frac{1}{2}\sigma^2)\tau \mp L\sqrt{\tau\sigma^2} $$

In [7]:

# Truncation range
def Range(S0, K, r, tau, sigma, L=12):
    c1 = sigma * sigma * tau
    x = np.log(S0 / K) + r * tau
    return (x - L * np.sqrt(c1), x + L * np.sqrt(c1))

# Truncation range using the cumulants from Table 11 in the paper
def RangePaper(S0, K, r, tau, sigma, L=12):
    c1 = (r - sigma * sigma * 0.5 ) * tau
    c2 = sigma * sigma * tau
    x = np.log(S0 / K) + c1
    return x - L * np.sqrt(c2), x + L * np.sqrt(c2)

Series¶

$$ F_k = Re[\phi(u)e^{iu(\log(\frac{S_0}{K})-a)}] $$$$ V_k = \alpha \frac{2}{b-a} K (\chi - \psi) $$$$ \sum_{k=0}^{N-1}\prime F_k V_k $$

In [8]:

class Side:
    CALL = 0
    PUT = 1

class Exercise:
    EUROPEAN = 0
    AMERICAN = 1
    
class Option:
   'Option'
   def __init__(self, side, exercise):
      self.side = side
      self.exercise = exercise
        
class Model:
    class Price:
        def __init__(self, value, delta):
            self.value = value
            self.delta = delta
        def __str__(self):
            return "V=" + str(self.value) + ", Δ=" + str(self.delta)

    def __init__(self, call = None, put = None):
        self.call = call
        self.put = put
    def __str__(self):
        return "Call:\r\n" + str(self.call) + ":\r\nPut:\r\n" + str(self.put)

In [9]:

def FourierCosineSeries(S0, K, r, tau, sigma, side, N=64, L=12, exercise=True):
 
    (a, b)  = Range(S0, K, r, tau, sigma, L) 
    
    if not exercise:
        (a, b)  = RangePaper(S0, K, r, tau, sigma, L) 
    
    (c, d) = (0, b)  if (Side.CALL == side) else (a, 0)
    
    accumulator = 0.
    for k in range(0, N):
        u = k * (np.pi / (b - a))
        exponential = np.exp(complex(0, u * (np.log(S0 / K) - a)))
        F_k = (Phi(u, r, tau, sigma) * exponential).real
        V_k = (1 if Side.CALL == side else -1) * 2 / (b - a) * K * (Chi(k, a, b, c, d) - Psi(k, a, b, c, d))
        accumulator += (0.5 if 0 == k else 1) * F_k * V_k
        
    return np.exp(-r * tau) * accumulator

In [10]:

FourierCosineSeries(100., 100., 0.03, 1.0, 0.40, Side.CALL, N=64, L=12)

Out[10]:

17.138735220514118

In [11]:

def error_absolute(x, x0):
    return np.abs(x-x0)
    
def error_relative(x, x0):
    return np.abs(x-x0)/x

In [12]:

def Analytical(side, S, K, r, T, sigma):
    d1 = (math.log(S / K) + (r+ sigma * sigma/ 2.0) * T) / (sigma * math.sqrt(T))
    d2 = d1 - sigma * math.sqrt(T)
    if Side.CALL == side:
        return (S * stats.norm.cdf(d1) - K * math.exp(-r * T) * stats.norm.cdf(d2), math.exp(0. * T) * stats.norm.cdf(d1))
    else:
        return (K * math.exp(-r * T) * stats.norm.cdf(-d2) - S * stats.norm.cdf(-d1), - math.exp(0. * T) * stats.norm.cdf(-d1))

In [13]:

xs = range(10, 151)
pyplot.plot(xs, list(map(lambda x: Analytical(Side.CALL, float(x), 100., 0.03, 1.0, 0.40)[0], xs)), label="Analytical")
#pyplot.plot(xs,list(map(lambda x: FourierCosineSeries(float(x), 100., 0.03, 1.0, 0.40, Side.CALL, N=15), xs)), label="F(N=15)")
#pyplot.plot(xs,list(map(lambda x: FourierCosineSeries(float(x), 100., 0.03, 1.0, 0.40, Side.CALL, N=20), xs)), label="F(N=20)")
pyplot.plot(xs,list(map(lambda x: FourierCosineSeries(float(x), 100., 0.03, 1.0, 0.40, Side.CALL, N=16), xs)), label="F(N=16)")
pyplot.xlabel("S0")
pyplot.ylabel("V")
pyplot.legend()

Out[13]:

In [14]:

xs = range(1,192)
a = Analytical(Side.CALL, 100., 110., 0.04, 1.0, 0.30)[0]
pyplot.yscale('log')



pyplot.plot(xs, 
            list(map(lambda x:
                error_relative(a,FourierCosineSeries(100.,110.,0.04,1.0,0.30,Side.CALL,N=x,L=12,exercise=True)),xs)),
            label="First truncation range")
pyplot.plot(xs, 
            list(map(lambda x:
                error_relative(a,FourierCosineSeries(100.,110.,0.04,1.0,0.30,Side.CALL,N=x,L=12,exercise=False)),xs)),
            label="Second truncation range")
pyplot.xlabel("N")
pyplot.ylabel("Relative error")
pyplot.legend()

Out[14]:

In [15]:

xs = range(1,64)
a = Analytical(Side.CALL, 100., 100., 0.03, 1.0, 0.40)[0]
#pyplot.yscale('log')

pyplot.plot(xs, 
            list(map(lambda x: (a), xs)), label="Analytical result")
pyplot.plot(xs, 
            list(map(lambda x: 
                (FourierCosineSeries(100., 100., 0.03, 1.0, 0.40, Side.CALL, N=x, L=12, exercise=True)), xs)),
            label="First truncation range")
pyplot.plot(xs, 
            list(map(lambda x: 
                (FourierCosineSeries(100., 100., 0.03, 1.0, 0.40, Side.CALL, N=x, L=12, exercise=False)), xs)), 
            label="Second truncation range")
pyplot.xlabel("N")
pyplot.ylabel("Option value")
pyplot.legend()

Out[15]:

Fourier Series Expansion for Option Pricing