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
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}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.
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 $$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
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)
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()
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)
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
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
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)
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)
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
FourierCosineSeries(100., 100., 0.03, 1.0, 0.40, Side.CALL, N=64, L=12)
def error_absolute(x, x0):
return np.abs(x-x0)
def error_relative(x, x0):
return np.abs(x-x0)/x
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))
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()
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()
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()