Option pricing techniques\(:\) A performance-based comparative study of the randomized quasi-Monte Carlo method and Fourier cosine method
Mohammad Rafiqul Islam
Florida State Universitymislam@math.fsu.com
2023-07-06
Abstract
Pricing financial derivatives such as options with desired accuracy can be hard due to the nature of the functions and complicated integrals required by the pricing techniques. In this paper we investigate the pricing methodology of the European style options using two advanced numerical methods, namely, Quasi-Monte Carlo and Fourier Cosine (COS). For the RQMC method, we use the random-start Halton sequence. We use the Black-Scholes-Merton model to measure the pricing quality of both of the methods. For the numerical results we compute the option price of the call option and we found a few reasons to prefer the RQMC method over the COS method to approximate the European style options.
Introduction
Option pricing is still an active research area in the field of
financial mathematics. The price of an option highly rely on the
volatility of the underlying asset and the market. There are many models
in the option pricing literature but many of them do not match with the
real-world situation. The very first model goes to Louis Bachelier
(Sullivan and Weithers (1991)), who
introduced the Brownian motion and used it to model the stock price in
1900. Bachelier’s model was based on the assumption that the stock price
process follows normal distribution. Later it was discovered that stock
price process follows log-normal distribution, rather normal (Antoniou et al. (2004)). Based on this the most
famous option pricing model was given by the Black-Scholes model. Again,
due the nature of some real-world complecated options BS model is not
sufficient (Janková (2018)). Therefore,
alternative methods have been developed, including the Randomized Quasi
Monte Carlo (RQMC) method and the Fourier COS method.
The RQMC method is a Monte Carlo-based technique that produces more
precise estimates of the option value by employing a low-discrepancy
sequence of random numbers. By increasing the approximation’s accuracy,
this method, as compared to traditional Monte Carlo simulations, reduces
error. In contrast, the Fourier COS method decomposes the option value
into a sum of cosine functions using the Fourier transform. For specific
types of options, such as European options, this strategy has proven to
be particularly effective and precise.
The performance of the RQMC and Fourier COS methods for pricing various
types of options will be compared and analyzed in this term paper. We
will start by going through the Black-Scholes model and the theory
underlying option pricing. The RQMC and Fourier COS approaches will next
be thoroughly covered, along with their benefits and drawbacks. The
accuracy and effectiveness of these methods compared to the traditional
Black Scholes method for pricing options with various characteristics
will next be demonstrated numerically.
Overall, in this term paper, we provide a comprehensive analysis of the
RQMC and Fourier COS methods for option pricing, with a focus on their
performance in comparison to the Black-Scholes model. Our findings
contribute to the growing body of literature on alternative methods for
option pricing and can help inform investment decisions in real-world
scenarios.
Mathematical Models and Theoretical Framework
Let \((\Omega,\mathcal F,\mathbb
Q)\) be a probability space where the filtration \(\{\mathcal F_t\}_{t\ge 0}\) contains all
the price information of the underlying asset \(S_t\) up to and including time \(t\), \(\mathbb
Q\) is the risk-neutral probability measure. Under this
probability measure, the Black-Scholes-Merton model follows the
following SDE
\[\begin{equation}\label{eq:bs}
dS_t=rS_tdt+S_t\sigma dW_t^{\mathbb Q}
\end{equation}\] where, \(r\)
and \(\sigma\) are the drift and
diffusion terms but here in this case it is the risk-free interest rate
and standard deviation, respectively, and \(W\) is the standard Brownian motion. The
solution to (\(\ref{eq:bs}\)) can be
obtained applying Itô’s formula to the function \(Z_t=f(t,S_t)=\ln S_t\). \[\begin{equation}\label{eq:ito}
\begin{split}
dZ_t&=\frac{\partial f}{\partial t}dt+\frac{\partial f}{\partial
S_t}dS_t+\frac{1}{2}\frac{\partial^2 f}{\partial S_t^2}(dS_t)^2\\
&=\frac{1}{S_t}(rS_tdt+S_t\sigma dW_t^{\mathbb
Q})+\frac{1}{2}\left(-\frac{1}{S_t^2}\right)\sigma^2S_t^2dt\\
&=\left(r-\frac{1}{2}\sigma^2\right)dt+\sigma dW_t^{\mathbb Q}
\end{split}
\end{equation}\] For \(Z_0=\ln
S_0\) we have \[\begin{equation}\label{eq:gbm}
S_t=S_0\cdot\exp{\left(\left(r-\frac{1}{2}\sigma^2\right)t+\sigma
W_t^{\mathbb Q}\right)}
\end{equation}\] Since forward start options are path-dependent
exotic options, therefore, equation (\(\ref{eq:gbm}\)) would be very useful for
both of the methods used to price these options.
Fourier Cosine Method (COS)
Any function supported on \([0,\pi]\) has the cosine expansion in the
following form
\[\begin{equation}\label{eq:cos1}
f(\eta)=\sideset{}{'}\sum_{k=1}^{\infty}\alpha_k\cdot\cos(k\eta)
\end{equation}\] where \[\alpha_k=\frac{2}{\pi}\int_{0}^{\pi}f(\eta)\cos
(k\eta)d\eta\]
where the notation \(\sideset{}{'}\sum\) implies that the
first term in the summation is weighted by one-half (Fang and Oosterlee (2009)). If the support of
the function is different than \([0,\pi]\), say \([a,b]\) then we use the following scale
variable \[\eta=\frac{x-a}{b-a}\pi,\hspace{4mm}
x=\frac{b-a}{\pi}\eta+a\]
Then we have the following form
\[\begin{equation}\label{eq:cos2}
\begin{split}
f(x)&=\sideset{}{'}\sum_{k=1}^{\infty}\alpha_k\cdot\cos\left(k\pi\frac{x-a}{b-a}\right)\\
\alpha_k&=\frac{2}{b-a}\int_{a}^{b}f(x)\cos\left(k\pi\frac{x-a}{b-a}\right)
\end{split}
\end{equation}\] Expectations or densities are defined on
infinite support but when we express in terms of the cosine expansion we
truncate the infinite support to a finite support yet we do not loose
accuracy significantly due to the conditions for the existence of a
Fourier transform, the integrands in densities have to decay to zero at
\(\pm\infty\).
Now suppose \([a,b]\in \mathbb R\) is
an interval such that the truncated integral approximates the infinite
integral
\[\begin{equation}\label{eq:cos3}
\hat{\phi}_X(\omega)=\int_{a}^{b}e^{i\omega y}f_{X}(y)dy\approx
\int_{\mathbb R}e^{i\omega y} f_{X}(y)dy=\phi_X(\omega)
\end{equation}\] Now we use Euler formula to simplify
further
\[e^{i\omega}=\cos \omega+i\sin
\omega\]
and the real part of \(e^{i\omega}\) is
\(\text{Re}(e^{i\omega})=\cos \omega\).
Therefore, for any random variable \(X\), and a constant \(a\in \mathbb R\) we have \[\begin{align*}
\phi_X(\omega)e^{ia}&=\mathbb E\left[e^{i(\omega X+ a)}\right]\\
&=\int_{-\infty}^{\infty}e^{i(\omega y+ a)}f_{X}(y)dy\\
\implies
Re\left(\phi_X(\omega)e^{ia}\right)&=Re\left(\int_{-\infty}^{\infty}e^{i(\omega
y+ a)}f_{X}(y)dy\right)\\
&=\int_{-\infty}^{\infty}\cos (\omega y+a)f_{X}(y)dy
\end{align*}\] Now we use the Fourier argument \(\omega=\frac{k\pi}{b-a}\) and multiply the
characteristic function in (\(\ref{eq:cos3}\)) by \(e^{-i\frac{ka\pi}{b-a}}\) we obtain \[\begin{align*}
\hat{\phi}_X\left(\frac{k\pi}{b-a}\right)\cdot
e^{-i\frac{ka\pi}{b-a}}&=\int_{a}^{b}e^{i\left(y\frac{k\pi}{b-a}-\frac{ka\pi}{b-a}\right)}f_X(y)dy\\
\implies Re\left(\hat{\phi}_X\left(\frac{k\pi}{b-a}\right)\cdot
e^{-i\frac{ka\pi}{b-a}}\right)&=Re\left(\int_{a}^{b}e^{i\left(y\frac{k\pi}{b-a}-\frac{ka\pi}{b-a}\right)}f_X(y)dy\right)\\
&=\int_{a}^{b}\cos\left(y\frac{k\pi}{b-a}-\frac{ka\pi}{b-a}\right)f_X(y)dy\\
&=\int_{a}^{b}\cos\left(k\pi\frac{y-a}{b-a}\right)f_X(y)dy
\end{align*}\] Therefore, matching with equation (\(\ref{eq:cos2}\)) \[\begin{equation}\label{eq:cos4}
\alpha_k=\frac{2}{b-a}\int_{a}^{b}\cos\left(k\pi\frac{y-a}{b-a}\right)f_X(y)dy=\frac{2}{b-a}Re\left(\hat{\phi}_X\left(\frac{k\pi}{b-a}\right)\cdot
e^{-i\frac{ka\pi}{b-a}}\right)
\end{equation}\] So from equation (\(\ref{eq:cos3}\)) we have \(\alpha_k\approx F_K\) where \[\begin{equation}\label{eq:cos5}
F_k=\frac{2}{b-a}Re\left(\phi_X\left(\frac{k\pi}{b-a}\right)\cdot
e^{-i\frac{ka\pi}{b-a}}\right)
\end{equation}\] By replacing \(\alpha_k\) by \(F_k\) we have \[\begin{equation}\label{eq:cos6}
\hat{f}_X(y)\approx \sideset{}{'}\sum_{k=0}^{\infty} F_k \cos
\left(k\pi\frac{y-a}{b-a}\right)
\end{equation}\] For the truncated version of equation (\(\ref{eq:cos6}\)) we will have \[\begin{equation}\label{eq:cos7}
\hat{f}_X(y)\approx \sideset{}{'}\sum_{k=0}^{N-1} F_k \cos
\left(k\pi\frac{y-a}{b-a}\right)
\end{equation}\]
Pricing European Options Under COS Method
We denote the random variable \(X\) as the log of the asset price, i.e., \(X_t=\log S_t\) and consider two points \(t\) and \(T\) in \([t,T]\) where \(t\) is the beginning of the interval and \(T\) is the maturity of the options and set the corresponding random variable as \(X_t=\log S_t=x\) and \(X_T=\log S_T=y\). Then the risk-neutral price of the option at time \(t\) can be given by the following formula \[\begin{equation}\label{eq:op1} P(t,x)=e^{-r(T-t)}\mathbb{E}^{\mathbb R}[V(T,y)\Big| \mathcal F_t]=e^{-r(T-t)}\int_{\mathbb R}V(T,y)f_{X_T}(y)dy \end{equation}\] where \(f_{X_T}(y)\) is the probability density function of the asset price process \(X_t\) and \(X_t\) normally distributed with with parameters \(\left((r-\frac{1}{2}\sigma^2)t,\sigma^2t\right)\). Therefore, we can truncate the integral as follows \[\begin{equation}\label{eq:op2} \hat{P}(t,x)=e^{-r(T-t)}\int_{a}^{b}V(T,y)f_{X_T}(y)dy \end{equation}\] as \(f_{X_T}(y)\to 0\) when \(y\to 0\). Now we represent the integral in equation (\(\ref{eq:op2}\)) in terms of Fourier Cosine (COS) method as follow \[\begin{equation}\label{eq:op3} \begin{split} \hat{P}(t,x)&=e^{-r(T-t)}\int_{a}^{b}V(T,y)\sideset{}{'}\sum_{k=0}^{\infty}\alpha_k(x)\cos\left(k\pi\frac{y-a}{b-a}\right)dy\\ &=\frac{b-a}{2}e^{-r(T-t)}\sideset{}{'}\sum_{k=0}^{\infty}\alpha_k(x)H_k \end{split} \end{equation}\] where, \[\begin{equation}\label{eq:op4} H_k=\frac{2}{b-a}e^{-r(T-t)}\int_{a}^{b}V(T,y)\cos\left(k\pi\frac{y-a}{b-a}\right)dy \end{equation}\] In a similar fashion in equation (\(\ref{eq:cos5}\)) we can approximate \(\alpha_k\) by the Fourier coefficient \(F_k\) and we get the option price as below\[\begin{equation}\label{eq:op5} P(t,x)\approx \hat{P}(t,x)=e^{-r(T-t)}\sideset{}{'}\sum_{k=0}^{N-1}Re\left(\phi_X\left(\frac{k\pi}{b-a}\right)e^{-ik\pi\frac{a}{b-a}}\right)H_k \end{equation}\] where, \(\phi_X(x)\) is the characteristic function of \(X\). The next thing is to find a closed form expression of the payoff function \(H_k\) to approximate the price. As Fang and Oosterlee (2009) showed in their paper, we can find that form as follows: The payoff function \(V(T,y)\) can be expressed as \[V(T,y)=\begin{cases}K(e^y-1)&\text{for a call option}\\K(1-e^y)&\text{for a put option}\end{cases}\] Now we need to use the following two results to deal with \(H_k\):
So for a call option we have
\[\begin{equation}\label{eq:op8}
\begin{split}
H_k^{\text{call}}&=\frac{2}{b-a}\int_{0}^{b} K(e^y-1)\cos
\left(k\pi\frac{y-a}{b-a}\right)dy\\
&=\frac{2}{b-a}K(\chi_k(0,b)-\psi_k(0,b))
\end{split}
\end{equation}\] Similarly, for a put option we have \[\begin{equation}\label{eq:op9}
\begin{split}
H_k^{\text{put}}&=\frac{2}{b-a}\int_{0}^{b} K(1-e^y)\cos
\left(k\pi\frac{y-a}{b-a}\right)dy\\
&=\frac{2}{b-a}K(-\chi_k(0,b)+\psi_k(0,b))
\end{split}
\end{equation}\] The last thing we need to compute the price of
an European call or put option in this COS method is the choice of
integration range: on what interval we integrate so that the truncation
error is not too big and the computational cost is not too much. As
proposed by Fang and Oosterlee (2009), the
efficient interval can be given as \[\begin{equation}\label{eq:op10}
[a,b]=[(x+\zeta_1)-L\sqrt{\zeta_2+\sqrt{\zeta_4}},(x+\zeta_1)+L\sqrt{\zeta_2+\sqrt{\zeta_4}}]
\end{equation}\] where, \(L\in
[6,12]\) is a number depending on a user-defined tolerance level
and the variables \(\zeta_1,\cdots,\zeta_4\) are the specific
cumulants of the underlying stochastic process, \(x=\frac{\log S_t}{K}\). An alternative form
of the integration range was proposed by Oosterlee and Grzelak (2019) as below \[\begin{equation}\label{eq:op11}
[a,b]=[-L\sqrt{T-t},L\sqrt{T-t}]
\end{equation}\] In the comparative analysis section we apply
this COS method to compare the the pricing accuracy with the randomized
quasi-Monte Carlo method along with the error analysis.
Higher Dimensional Randomized Quasi Monte Carlo Method
Pricing of a security derivative involves the computation of the
discounted expected value of the payoff function. In continuous time
setting, it is evaluating the integrals for the expectation. Suppose we
want to get an estimate \(\theta\) for
the expected payoff of a derivative. Then we can express this as follows
\[\mathbb{E}[\text{payoff}]=\theta=\int_{I^s}f(x)dx\]
where \(I^s\) is the \(s\)- dimensional unit cube. The primary
benefit of Monte Carlo and quasi-Monte Carlo techniques is that they
allow us to assess integrals using sums and to utilize statistical
techniques to analyze the error in the context of uniformly distributed
modulo 1 sequences (Ökten and Eastman
(2004)). That is, we can express the integrals as sums of the
form
\[\int_{I^s}f(x)dx=\lim_{Nt\to
\infty}\frac{1}{N}\sum_{n=1}^{N}f(q_n)\] where, \(\{q_n\}_{n\ge1}\) is \(s-\)dimensional low discrepancy sequence.
The quasi-Monte Carlo method differs with regular Monte Carlo in the
sense that \(q_n\) is a \(s-\)dimensional random vector for the MC
method and on the contrary, it is the \(n\)th term of of a \(s-\)dimensional low-discrepency sequence
for the QMC method.
We start developing the theoretical background with the definition of a
uniformly distributed sequence of modulo \(1\)
We can extend the definition from one dimension to a finite dimension by the following theorem
The next thing is to know what we mean by a low discrepancy sequence. If we have a sequence \(\{x_n\}_{n\ge 1}\in I^s\) its extreme discrepancy is defined as\[D_N(x_n)=\sup_{[a,b)}\Big|\frac{1}{N}\sum 1_{[a,b)}(x_n)-\lambda_s([a,b)\Big|\] If the supremum is taken over all the intervals of the form \([0,a)\) then the extreme discrapncy is said to be star-discrepancy, i.e., \[D^*_N=\sup_{[0,a)}\Big|\frac{1}{N}\sum 1_{[0,a)}(x_n)-\lambda_s([0,a)\Big|\]
The following theorem relates the star-discrepancy with u.d. mod \(1\) sequence.
Now we are ready to introduce the \(s-\)dimensional randomized quasi-Monte
Carlo method to evaluate integrals by the quadrature rule \[\begin{equation}\label{eq:quad}I\approx
Q(q_u)=\frac{1}{N}\sum_{n=1}^{N}f(q_u^{(n)})
\end{equation}\] where \(u\) is
the randomization index. Our goal is to find the expected value of \(Q(q_u)\) and reduce the variance of \(Q(q_u)\) as much as possible by generating
low discrepancy sequences. By Ökten and Eastman
(2004), we want the discrepancy to be as small as possible, since
the accuracy of our quadrature rule (\(\ref{eq:quad}\)) depends on it by the
Koksma-Hlawka inequality \[\begin{equation}\label{eq:kh}
\Big|Q(q_u)-I\Big|\le V_{HK}(f)D^*_N(q_u)
\end{equation}\] where, \(V_{HK}(f)\) is the variation of \(f\) in the sense of Hardy and Krause (Niederreiter (1992)). The inequality in (\(\ref{eq:kh}\)) tells us that the more
smaller variation of \(f\) and
star-discrepancy the better accuracy.
There are a quite good number of low-discrepancy sequences in the
literature but for our forward start option we apply mainly two
low-discrepancy sequences, random-start Halton sequence and scrambled
Faure sequences.
Random-start Halton Sequence
Random-start Halton sequence is based on the van der Corput and Halton sequence. For any base \(b\), any positive integer \(n\) can be expressed as \[n=(a_k,a_{k-1},\cdots,a_1,a_0)b=a_0+a_1b+a_2b^2+\cdots+a_kb^k\] then the van der Corput sequence in base \(b\) is defined as \[\begin{equation}\label{eq:vdc} \begin{split} \phi_b(n)&=(.a_0,a_{1},\cdots,a_{k-1},a_k)b\\ &=\sum_{i=0}^{k}\frac{a_i}{b^{i+1}} \end{split} \end{equation}\] The Halton sequence is a generalization of the van der Corput sequence to the higher dimension. The random start idea was introduced by Struckmeier (1995) that relies on an alternative derivation of the Halton sequence which we call given by Lambert (1985) based on von Neuman–Kakutani transformation.If \(\mathbf b=(b_1,b_2,\cdots,b_s)\) be a vector of positive integers that are pairewise relatively prime and we have a vector \(\mathbf x \in I^s\) then the \(s-\)dimensional von Neuman-Kakutani transformation is given by \[\begin{equation}\label{eq:vnk} \mathbf T_{\mathbf b}(\mathbf x)=(T_{b_1}(x_1),T_{b_2}(x_2),\cdots,T_{b_s}(x_s)) \end{equation}\] For \(b=\frac{1}{2},\hspace{1mm} \mathbf T_{\mathbf b}(\mathbf x)\) is the Halton sequence in base \(b_1,b_2,\cdots,b_s\) if the orbit vector \(\mathbf x\) is set to \(\mathbf 0\). Let \(q_u\) be the orbit if a randomly selected “starting values” \(\mathbf u\in U(0,1)^s\) is used under \(\mathbf T_{\mathbf u}\). Then the sequence we get will be the random-start permuted Halton sequence. So we redefine the quadrature equation (\(\ref{eq:quad}\)) as follows \[\begin{equation}\label{eq:rshquad} Q(q_{\mathbf u})=\frac{1}{N}\sum_{i=0}^{N-1}f\left(\mathbf T_{\mathbf b}^i(\mathbf u)\right) \end{equation}\] By Wang and Hickernell (2000), the quadrature rule is provided by the following theorem.
The discrepancy for this \(q_{\mathbf u}^{(n)}\) sequences is given by Bouleau and Lepingle (1994) (p. 77) \[D_N(q_{\mathbf u})\le \left(1+\Pi_{i=1}^{k}(b_i-1)\frac{\log {(b_iN)}}{\log b_i}\right)\]
Comparative Analysis
In this section we provide the performance of the both method. We show our results for only for the European style call options. However, both of the methods are designed to price both put or call options. We consider a hypothetical non-dividend paying stock with the initial price \(S_0=100\), risk-free interest rate rate \(r=0.1\), volatility \(\sigma=0.25\), and time to maturity \(\tau=T-t=0.1\).
COS Method Implementation
First we set the integration interval as mentioned in equation (\(\ref{eq:op11}\)) and choose the parameter \(L=10\) as mentioned in that section. We present the results for three different strike prices \(K=80,100,\) and \(120\).
COS method prices are shown against the BS prices.
The figure above (\(\ref{fig1}\)) depicts the price computed for three different strikes \(K=80,100,\) and \(120\) where the blue solid line represents the price computed by the BS model and the red dashed line repesents the prices computed by the COS method. Next we represent the absolute error and CPU time usuage for the COS method as the measure of the performance.
| N | N16 | N32 | N64 | N128 |
|---|---|---|---|---|
| msec. | 3.66e-01 | 3.710000e-01 | 4.164400e-01 | 4.241371e-01 |
| abs. err. K=80 | 3.66e+00 | 7.441700e-01 | 1.922192e-02 | 1.305295e-07 |
| abs. err. K=100 | 3.88e+00 | 5.279164e-01 | 1.515078e-02 | 3.868614e-07 |
| abs. err. K=120 | 3.17e+00 | 8.132994e-01 | 2.135178e-02 | 3.496943e-07 |
As we can see from table (\(\ref{tab1}\)), the number of point estimation \(N\) does not increase significantly. Also the errors reduces as we increase \(N\).
RQMC Method Implementation
In this section we present a similar numerical results that we provided for the COS method. For the inverse transformations, we use the box muller transformation. Here we generate \(N=10,000\) asset price paths and estimate 40 of the option prcies for \(K=100\) and we obtain the following figure
40 Estimates of the option price for K=100 by RQMC method
The mean of the 40 estimates is very close to the true mean by Black-Scholes option price for case \(K=100\). We run few more simulations changing \(K\) values and obtained similar close approximations. However, run-time is significantly higher for the RQMC approximations.
Comparision
Now we present the combined summary of the numerical experiments we had so far. In table (\(\ref{tab2}\)) and (\(\ref{tab3}\)) we see that both of the methods works very well. However, run time is higher for the randomized quasi-Monte Carlo method, but absolute error is less than that of the COS method.
| K | BS | COS | COS.Abs.Err | COS.RT.msec |
|---|---|---|---|---|
| 80 | 20.80 | 21.54 | 7.4e-01 | 0.36600 |
| 100 | 3.66 | 3.13 | 5.3e-01 | 0.37100 |
| 120 | 0.04 | 0.85 | 8.1e-01 | 0.41644 |
| K | BS | RQMC | RQMC.Abs.Err | RQMC.RT.sec |
|---|---|---|---|---|
| 80 | 20.80 | 20.81 | 1.0e-02 | 8.32 |
| 100 | 3.66 | 3.73 | 7.0e-02 | 8.38 |
| 120 | 0.04 | 0.18 | 1.4e-01 | 8.69 |
In figure (\(\ref{fig3}\)) we see the combined option price estimations for those three strike prices. We notice that the RQMC gives more closer approximation to the BS estimate.
Price comparision of BS method, COS method, and RQMC method
Conclusion
In conclusion, we conclude that RQMC provides a better approximation for estimating the European option price. Even though it takes a bit longer to converge, still RQMC may be preferred over the COS method. Since using a large number of point estimations in the COS method provides overfitting and a too-small interval of the estimation, therefore, RQMC method may be chosen over the COS method. Moreover, the COS method requires more parameters to estimate the option price, and choosing the truncation interval parameter \(L\) may cause an increase in the error estimation.