White paper
Comparing option pricing methods in q
In this paper, we compare the use of both Monte Carlo (MC) and QuasiMonte Carlo (QMC) methods in the process of pricing European and Asian options. In doing so, we consider the use of two discretization schemes  standard discretization and Brownianbridge construction. Results produced by the different methods are compared with the deterministic BlackScholes price for each option type, using Global Sensitivity Analysis (SA). Note that the methods demonstrated below follow the work presented by S. Kucherenko et al. 2007.
S. Kucherenko et al. 2007, “The Importance of Being Global – Application of Global Sensitivity Analysis in Monte Carlo Option Pricing”, Wilmott, pp. 82–91
BlackScholes
The most common model used to calculate the price of options is BlackScholes, where the formula for each market is derived from the BlackScholes equation. In this paper, we look specifically at the BlackScholes models for European and Asian call options. The standard BlackScholes model for European options assumes a payoff based on the underlying price at exercise. The modified model for Asian options assumes a payoff based on the average underlying price over a predefined time period. In each case, the BlackScholes model produces a closedform solution with a deterministic result.
For European call options, the price of the corresponding option at time \(t\), \(P(S_{t},t)\), is given by:
Where \(T\) is the expiry, \(S_{t}\) is the price of the underlying asset at time \(t\), \(K\) is the strike price of the option, \(\sigma\) is the volatility and \(r\) is the interest rate. Note that the price is discounted by the dividends, \(q\), throughout.
For Asian call options, we implement the same formula, using an adjusted \(S_{t}\), \(\sigma^{2}\) and drift rate, \(\mu\):
Where \(n\) is the number of timesteps.
Monte Carlo and QuasiMonte Carlo simulations
Within the financial industry, there is a need to price complex financial instruments. Despite this need, there are a lack of analytical solutions to do so. MC methods are used with the financial industry to mimic the uncertainty associated with the underlying price of an instrument and to subsequently generate a value based on the possible underlying input values. One example of where MC is used in finance, is in evaluating an option on an equity. For each underlying asset, an MC simulation is used to create thousands of random price paths, with an associated payoff. The option price for each path is calculated by taking the average over the future payoffs and discounting them to the present.
These models are based on pseudorandom numbers which, despite being commonly used, exhibit very slow convergence, with a rate of \(O(1/\sqrt{N})\) where \(N\) is the number of sampled points. To improve upon these models, QMC methods have been developed which use lowdiscrepancy sequences (LDS) to produce a rate of convergence ~ \(O(1/N)\). LDS are deterministic uniformly distributed sequences which are specifically designed to place sample points as uniformly as possible. Practical studies have shown that the most effective QMC method for applicaton in financial engineering is based on Sobol' LDS.
S. Kucherenko et al. 2001, “Construction and Comparison of HighDimensional Sobol’ Generators”, Wilmott, Nov, pp. 6479
broda.co.uk
P. Jäckel 2001, Monte Carlo Methods In Finance, pp. 122.
P. Glasserman 2003, Monte Carlo Methods in Financial Engineering, Springer.
Wiener path
The starting point for asset price simulation is the construction of a Wiener path (or Brownian motion). Such paths are built from a set of independent Gaussian variates, using either standard discretization or Brownianbridge construction.
In the standard approximation, the Wiener path is found by taking the cumulative sum of the Gaussian variates.
When constructing a Brownian bridge, the last step of the Wiener path is calculated first, followed by the midstep, and then the space left between steps is bisected until all steps have been determined.
An example of building up a Brownian bridge is shown in the diagram below, where we have a total of 14 timesteps (from 1 to 14). Note that we also have an additional timestep 0
, which is assumed to have a value of 0.
The construction of a Brownian bridge over 14 steps. See Jäckel, 2001, op. cit.
Both standard discretization and Brownianbridge construction share the same variance and therefore the same resulting convergence when used with MC models. However, performance differs between the two when QMC methods are introduced, with faster convergence seen for Brownianbridge construction.
In order to showcase how performant the QMC simulation is when paired with Brownianbridge construction, we use Global SA as outlined in S. Kucherenko et al. 2007. This method allows us to estimate the contribution of individual input parameters in the final variance of the output over a number of experiments. In each experiment, we:
 Randomly generate \(n\) random numbers, either pseudorandom (MC) or Sobol’ sequence (QMC). (See also broda.co.uk)
 Convert into a normal distribution.
 Convert into a Wienerpath random walk using standard discretization or Brownianbridge construction.

Convert into an assetprice path based on parameters:
s
: Asset price at \(t=0\)v
: Volatilityr
: Interest rateq
: Dividendst
: Expiry

Convert into an option price based on the option type and strike price,
k
.
The prices produced are then averaged to find a final predicted price.
Implementation
In the following sections, we compare the methods of option pricing mentioned above. The BlackScholes price for each market is compared to an average price generated using the following combinations of simulation and discretization methods:
 Pseudorandom number generation (MC) with standard discretization.
 Sobol’ sequences (QMC) with standard discretization.
 Sobol’ sequences (QMC) with Brownianbridge construction.
The BlackScholes function for each market produces a closedform solution with a deterministic result, while the MC/QMC functions perform a number of random experiments and return an average price, based on the option type and the strike price.
Once both the BlackScholes and MC/QMC prices have been calculated for each market, the root mean square error (RMSE) is calculated between the two. This is demonstrated in the final example below, where the process is repeated for an increasing number of paths, with resulting errors compared.
The technical dependencies required for the below work are as follows:
 Option Pricing kdb+/q library
 embedPy
 Sobol’ C++ library 
SobolSeq1024
function provided in the Option Pricing kdb+/q library with max dimension of 1024.  matplotlib
Utility functions
For simplicity, utility functions are omitted from the code snippets below. These can be found within the Option Pricing library linked above.
Load scripts
As mentioned previously, the implementations of option pricing methods outlined below are based on original C++ scripts used in S. Kucherenko et al. 2007. All code is contained within the optionpricing repository:
Wrappers for the C++ pseudorandom and Sobol’ sequence number generators (see also broda.co.uk) are contained within rand.q
, along with the cumulative and inverse cumulative normal distribution functions in norm.q
.
To run the below examples, q scripts are loaded including the C++ wrappers and graphics functions used throughout.
\l code/q/rand.q
\l code/q/norm.q
\l notebooks/graphics/graphics.q
BlackScholes option pricing
The following functions provide q implementations of the BlackScholes formula for each market, outlined above. They take in parameter dictionary pd
as an argument, containing the parameters s
, v
, r
, q
, k
, and t
detailed in the previous section. Note that the BlackScholes price of an Asian call option depends on the number of timesteps n
, which must also be passed as an argument.
/ European
bsEuroCall:{[pd]
/ Calculate volatility*sqrt delta T coefficient
coeff:(v:pd`v)*sqrt t:pd`t;
/ Calculate d1
d1:(log[pd[`s]%pd`k]+t*(pd[`r]pd`q)+.5*v*v)%coeff;
/ Calculate d2
d2:d1coeff;
/ Calculate the option price  P(S,t)
(pd[`s]*exp[neg t*pd`q]*cnorm1 d1) 
pd[`k]*exp[neg t*pd`r]*cnorm1 d2 }
/ Asian
bsAsiaCall:{[n;pd]
/ Calculate adjusted drift rate
adjmu:.5*((r:pd`r).5*v2:v*v:pd`v)*n1:1+1.%n;
/ Calculate adjusted volatility squared
adjv2:(v2%3)*n1*1+.5%n;
/ Calculate adjusted price
adjS :pd[`s]*exp(t:pd`t)*(hv2:.5*adjv2)+adjmur;
/ Calculate d1
d1:(log[adjS%k:pd`k]+t*(rq:pd`q)+hv2)%rtv2:sqrt adjv2*t;
/ Calculate d2
d2:d1rtv2;
/ Calculate the option price  P(S,t)
(adjS*exp[neg q*t]*cnorm1 d1)k*exp[neg r*t]*cnorm1 d2 }
The outputs of these functions are demonstrated below for 512 timesteps.
nsteps:512 / number of timesteps
pd:`s`k`v`r`q`t!100 100 .2 .05 0 1 / parameter dictionary
/ Calculate BS price for European/Asian options
"European Black Scholes Price: ",string bseuro:bsEuroCall pd
"Asian Black Scholes Price: ",string bsasia:bsAsiaCall[nsteps]pd
European Black Scholes Price: 10.45058
Asian Black Scholes Price: 5.556009
Monte Carlo and QuasiMonte Carlo option pricing
Generate random numbers
The first stage in predicting an option price is to generate a set of random numbers using either MC or QMC methods. In the example below we generate 512 pseudorandom and Sobol’ sequence numbers, with results plotted for comparison.
Random numbers are generated using the Mersenne Twister number generator which has one parameter, the number of steps.
The Sobol’ sequence generator takes two arguments:
 the index of the point (0 < i < 2^{31}  1)
 the dimension of the Sobol’ sequence, i.e. the number of steps (0 < d < 1025).
/ Function to generate n random numbers in d dimensions
rdmngen:{[n;d](d;n)#mtrand3 d*n}
/ Function to generate n sobol numbers in d dimensions
sobngen:{[n;d]flip sobolrand[d]each 1+til n}
/ Generate n random and sobol numbers in 2D
data:(rdmngen;sobngen).\:nsteps,2
subplot[data;("Random";"Sobol");2;2#`scatter]
It is clear that the pseudorandom numbers are not evenly distributed, with points clustering together in some sections, while leaving large portions of white space in others.
In contrast, the Sobol’ sequence plot exhibits a much more even distribution, with few points clumping together.
Convert to a Gaussian distribution
The generated sequences are converted from a uniform distribution to a Gaussian distribution. Following this conversion, around 68% of the values lie within one standard deviation, while two standard deviations account for around 95% and three account for 99.7%.
Gaussian distribution
In the example below we convert the uniform generated Sobol’ sequence to a Gaussian distribution, using the inverse cumulative normal function, invcnorm
.
/ Convert sobol sequence to normal distribution
zsob:invcnorm each sob:last data
subplot[(sob;zsob);("Sobol Uniform";"Sobol Gaussian");2;2#`scatter]
The differences between the Gaussian distributions produced for random and Sobol’ sequences are best demonstrated for a small number of timesteps, e.g. 64. Below we plot the 1D Gaussian distributions for both random and Sobol’ number generation across 64 timesteps.
/ Returns 1D Gaussian distribution to plot
gausscnv:{[g;n;d]first invcnorm each$[g~`rdm;rdmngen;sobngen][n;d] }
/ Calculates Gaussian variates for 64 steps, in 2 dimensions
dist:gausscnv[;64;2]each`rdm`sob
subplot[dist;("Random";"Sobol");2;2#`hist]
As expected, the Sobol’ sequence exhibits a Gaussian curve with much better statistical properties than the randomnumber sequence.
Convert into a Wienerpath random walk
The q code to build both a Brownianbridge and Wienerpath random walk is shown below.
Brownian bridge:
/* n = number of timesteps
/* dt = length of timesteps
bbridge:{[n;dt]
/ create initial brownian bridge with indices for all timesteps
bb:first flip(n1).[util.initbb n]\(`bidx`ridx`lidx!3#n1;((n1)#0b),1b);
/ calculate weights and sigma value for each point in the path
bb:update lwt:bidxlidx,rwt:ridxbidx,sigma:ridxlidx from bb;
bb:update lwt%sigma,rwt%sigma,sigma:sqrt dt*lwt*rwt%sigma from bb;
/ create a projection for weiner path creation containing new bbridge
util.buildpath .[bb;(0;`sigma);:;sqrt n*dt] }
Wiener path:
/ Performs cumulative sum
/ or inverse cumulative normal for random/sobol numbers
/* u = gaussian variates
/* d = dictionary containing bbridge and boolean for sobol/random numbers
wpath:{[u;d]$[(::)~d`bb;sums;d`bb]invcnorm u }
An example of how the Brownian bridge is built is shown below using bbdemo
. The function outputs a table with n
timesteps.
/ Demonstrates build up of bbridge indices
bbdemo:{[n]
/ Create matrix showing steps taken
/ Step already taken = 1b, 0b otherwise
x:1b,'enlist[n#0b],
last flip(n1).[util.initbb n]\(`bidx`ridx`lidx!3#n1;((n1)#0b),1b);
/ Print "X" where 1b, showing path taken
flip(`$"i",'string til count x)!x:flip(" X")x }
An example is shown below using 8 timesteps, showing the order in which steps are added to the path. Note that i0
was added here, where we assume that it has a value equal to 0.
q)bbdemo 8
i0 i1 i2 i3 i4 i5 i6 i7 i8

X
X X
X X X
X X X X
X X X X X
X X X X X X
X X X X X X X
X X X X X X X X
X X X X X X X X X
When recording the order of steps in the path, we take note of the left and right weights and indexes, and the corresponding sigma value for each step in the sequence. This is shown for 512 timesteps and 1 unit of time, with the sigma value for each index in the Brownian bridge subsequently plotted.
q)dt:1
q)10#b:last value bbex:bbridge[nsteps;dt]
bidx ridx lidx lwt rwt sigma

511 511 511 22.62742
255 511 1 0.5 0.5 11.31371
127 255 1 0.5 0.5 8
383 511 255 0.5 0.5 8
63 127 1 0.5 0.5 5.656854
191 255 127 0.5 0.5 5.656854
319 383 255 0.5 0.5 5.656854
447 511 383 0.5 0.5 5.656854
31 63 1 0.5 0.5 4
95 127 63 0.5 0.5 4
Once the Brownian bridge has been intialized, it can be used to transform Gaussian variates into a Wienerpath random walk. Below, a Wiener path with 512 timesteps is constructed using a Sobol’ sequence (of length 512) and the Brownian bridge constructed previously. Note that the function wpath
takes two arguments:
 Sequence of generated numbers, Sobol’ or random.
 Dictionary indicating whether to use standard discretization or Brownianbridge construction, and whether to use Sobol’ sequences (
1b
) or pseudorandom numbers (0b
). If using a Brownian bridge, the initial Brownian bridge must be passed in, if not use(::)
.
q)d:`bb`sobol!(bbex;1b)
q)show w:wpath[sobolrand[nsteps;2];d]
0.4450092 0.06385387 0.1017726 1.221271 0.9099617 1.552524 0.56..
q)plt[`:title]"Wiener path random walk";
q)plt[`:plot]w;
q)plt[`:show][];
Convert into asset price path
At this point, the Wiener path is converted into an assetprice path using the methods outlined in S. Kucherenko et al. 2007, where the assetprice path is calculated as:
Where \(S_0\) and \(S(t)\) are the asset prices at time \(0\) and \(t\) respectively, \(r\) is the interest rate, \(\sigma\) is the volatility and \(W(t)\) is a Wiener path up to time \(t\).
This process can be done using the function spath
, detailed below.
/* u = guassian variates
/* n = number of timesteps
/* d = dictionary with bbridge and boolean for random/sobol
/* pd = dictionary of parameters s,v,r,t,q
spath:{[u;n;d;pd]
/ Original asset price
pd[`s]*
/ Wiener path*volatility*time for one timestep
exp(wpath[u;d]*pd[`v]*sqrt dt)+
/ Calculate sum of interest rate (discounted by dividends)
/ and half volitility squared for each timestep
(1+til n)*(pd[`r]pd[`q]+.5*v*v:pd`v)*dt:pd[`t]%n }
Here we calculate six different assetprice paths and overplot them for comparison. We start by generating the Sobol’ sequences for 8 paths with 512 timesteps, incrementing the Sobol’ index each time. Brownianbridge approximation is also used.
1"\nGenerated sequences: \n";
show u:sobolrand[nsteps;]each 2+til 8
plt[`:title]"Asset Price Path"
plt[`:plot] each spath[;nsteps;d;pd]each u
plt[`:show][]
Generated sequences
0.25 0.75 0.25 0.75 0.25 0.75 0.25 0.25 0.75 0.75 ..
0.75 0.25 0.75 0.25 0.75 0.25 0.75 0.75 0.25 0.25 ..
0.375 0.625 0.125 0.875 0.875 0.125 0.625 0.125 0.875 0.625 ..
0.875 0.125 0.625 0.375 0.375 0.625 0.125 0.625 0.375 0.125 ..
0.125 0.375 0.375 0.125 0.625 0.875 0.875 0.375 0.125 0.375 ..
0.625 0.875 0.875 0.625 0.125 0.375 0.375 0.875 0.625 0.875 ..
0.3125 0.3125 0.6875 0.5625 0.1875 0.0625 0.9375 0.5625 0.0625 0.8125..
0.8125 0.8125 0.1875 0.0625 0.6875 0.5625 0.4375 0.0625 0.5625 0.3125..
Convert into option price
Lastly, to find a single option price, an average is taken across the assetprice path for the MC/QMC method. This allows for comparison between the predicted price and the BlackScholes equivalent. The formulae for both the European and Asian options are outlined in S. Kucherenko et al. 2007.
For a European call option, the final MC/QMC price is calculated using:
Where \(C\) is the final price of the call option, \(r\) is the interest rate, \(T\) is the length of timestep, \(N\) is the finite number of simulated price paths and \(K\) is the strike price. Note that \(max(S^{(i)}_{T}K,0)\) represents the payoff for a call option.
This has been translated into the below function:
/ MC/QMC Price of European Call Option
/* m = number of paths
/* u = sequence of random numbers
/* d = dictionary with bbridge and boolean for random/sobol
/* pd = dictionary of parameters s,v,r,t,q
mcEuroCall:{[u;n;d;pd]
exp[neg pd[`r]*pd`t]*avg 0 
(last each spath[;n;d;pd]each u)pd`k }
Similarly, for Asian call options, the below is used:
Where we integrate over unit hypercube \(H^{n}\). In this case, the payoff for an geometric average Asian call option is calculated as the maximum between 0 and the geometric average of the underlying price, \(S(t)\), minus the strike price, \(K\).
The final price for an Asian call option can therefore be determined generated by using the below function:
/ MC/QMC Price of Asian Call Option
/* m = number of paths
/* n = number of timesteps
/* d = dictionary with bbridge and boolean for random/sobol
/* pd = dictionary of parameters s,v,r,t,q
mcAsiaCall:{[u;n;d;pd]
exp[neg pd[`r]*pd`t]*avg 0 
(prd each xexp[;1%n]spath[;n;d;pd]each u)pd`k }
We also need a numbergenerator function for l
trials, m
paths and n
steps which can be used with the Sobol’ or randomnumber generators.
numgen:{[ng;l;m;n]ng@''$[ng~mtrand3;(l;first m)#n;(0N;m)#1+til l*m]}
Here we demonstrate how to run these functions below for 512 timesteps, 256 paths and 5 trials. Sequences are generated for Sobol’ sequences using the above numgen function which will produce a sequence for each path and each trial.
ntrials:5
npaths:256
"\nGenerated sequences:\n"
5#u:first numgen[sobolrand nsteps;ntrials;npaths;nsteps]
"European Monte Carlo Price: ",
string mcEuroCall[u;nsteps;`bb`sobol!(bbex;1b);pd]
"Asian Monte Carlo Price: ",
string mcAsiaCall[u;nsteps;`bb`sobol!(bbex;1b);pd]
Generated sequences:
0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5..
0.25 0.75 0.25 0.75 0.25 0.75 0.25 0.25 0.75 0.75 0.25 0.2..
0.75 0.25 0.75 0.25 0.75 0.25 0.75 0.75 0.25 0.25 0.75 0.7..
0.375 0.625 0.125 0.875 0.875 0.125 0.625 0.125 0.875 0.625 0.125 0.3..
0.875 0.125 0.625 0.375 0.375 0.625 0.125 0.625 0.375 0.125 0.625 0.8..
European Monte Carlo Price: 10.28224
Asian Monte Carlo Price: 5.365942
Remembering that the BlackScholes option prices for the same number of timesteps were:
"European Black Scholes Price: ",string bseuro
"Asian Black Scholes Price: ",string bsasia
European Black Scholes Price: 10.45058
Asian Black Scholes Price: 5.556009
Example
In this section we deploy all the aforementioned techniques and compare the results.
Multiple threads
The example below can be run from the terminal across multiple threads using the following commands:
$ q s 8
q)\l op.q
q)loadfile`:init.q
q)loadfile`:code/q/run.q
where we load in the functions contained within the Option Pricing library using the first two commands and run the example by loading in run.q
.
Parameters
As shown previously, a dictionary of parameters is created which contains the initial asset price s
, volatility v
, interest rate r
, dividends q
, expiry t
and strike price k
.
q)show pd:`s`k`v`r`q`t!100 100 .2 .05 0 1
s 100
k 100
v 0.2
r 0.05
q 0
t 1
Additional parameters are also initialized for the number of paths (experiments), steps and trials.
q)l:20 / Number of trials
q)n:1024 / Number of steps
q)show m:"j"$xexp[2;3+til 8] / Number of paths
8 16 32 64 128 256 512 1024
Given that the initial Brownian bridge is the same throughout, it is also initialized and passed in as an argument.
q)10#last value bb:bbridge[n;1]
bidx ridx lidx lwt rwt sigma

1023 1023 1023 32
511 1023 1 0.5 0.5 16
255 511 1 0.5 0.5 11.31371
767 1023 511 0.5 0.5 11.31371
127 255 1 0.5 0.5 8
383 511 255 0.5 0.5 8
639 767 511 0.5 0.5 8
895 1023 767 0.5 0.5 8
63 127 1 0.5 0.5 5.656854
191 255 127 0.5 0.5 5.656854
Run experiments
The functions below calculate the RMSE between the BlackScholes and MC/QMC prices for each market and each MC/QMC technique. Note that we reset the sobolrand
index after each set of trials have been completed.
/ Run all techniques for option pricing
/* bb = initial brownian bridge
/* pd = dictionary of parameters
/* l = number of trials
/* n = number of timesteps
/* m = number of paths
runall:{[bb;pd;l;n;m]
/ Start timer for European options
st:.z.p;
/ Output column names
0N!util.rcol;
/ Run experiments for European options
e:util.run[`euro;bsEuroCall pd;bb;pd;l;n]each m;
/ Output total time taken for European
1"European: time taken = ",string[.z.pst],"\n";
/ Start timer for Asian options
st:.z.p;
/ Output column names
0N!util.rcol;
/ Run experiments for Asian options
a:util.run[`asia;bsAsiaCall[n;pd];bb;pd;l;n]each m;
/ Output total time taken for Asian
1"Asian: time taken = ",string .z.pst;
/ Return table with European and Asian prices and errors
e,a }
/ Dictionary keys
util.d:`bb`sobol!
/ Output column names
util.rcol:`mkt`npaths`rmse_bb_sobol`rmse_std_sobol`rmse_std_rdm,
`prx_bb_sobol`prx_std_sobol`prx_std_rdm`prx_bs
/ RMSE
util.rmse:{sqrt avg x*x:y}
/ Run each technique for a specific market
/* mkt = market, European/Asian
/* bs = BlackScholes price
util.run:{[mkt;bs;bb;pd;l;n;m]
/ Create project with correct MC function for each market
mc:$[mkt~`asia;mcAsiaCall;mcEuroCall][;n;;pd];
/ Generate MC option price and calculate error for bbridge and sobol
ea:util.rmse[bs]a:mc[;util.d(bb;1b)]each sob:numgen[sobolrand n;l;m;n];
/ Generate MC option price and calculate error for standard and sobol
eb:util.rmse[bs]b:mc[;util.d(::;1b)]each sob;
/ Generate MC option price and calculate error for bbridge and random
ec:util.rmse[bs]c:mc[;util.d(bb;0b)]each numgen[mtrand3;l;m;n];
/ Return dictionary of results
util.rcol!0N!(mkt;m;ea;eb;ec;last a;last b;last c;bs) }
Compare results
At this stage it is possible to plot the results obtained for the option prices, RMSE and log RMSE values.
q)r:runall[bb;pd;l;n;m]
q)select from r where mkt=`euro / European RMSEs and prices
mkt npaths rmse_bb_sobol rmse_std_sobol rmse_std_rdm prx_bb_sobol pr..
..
euro 8 2.218523 3.328543 4.537168 13.86836 4...
euro 16 1.345787 2.442911 3.505794 9.206011 12..
euro 32 0.6865024 1.623545 3.618555 9.879788 11..
euro 64 0.3774031 0.9891046 1.93851 10.90519 11..
euro 128 0.2089234 0.5977986 1.532864 10.34505 11..
euro 256 0.117329 0.3648233 0.9575077 10.52265 10..
euro 512 0.05984563 0.3127605 0.7618771 10.50504 9...
euro 1024 0.03176637 0.2853521 0.5670563 10.43112 10..
q)/ Asian RMSEs and prices
q)select from r where mkt=`asia
mkt npaths rmse_bb_sobol rmse_std_sobol rmse_std_rdm prx_bb_sobol pr..
..
asia 8 1.044296 2.126291 2.421675 6.461112 5...
asia 16 0.6879741 1.37292 1.80831 4.369775 6...
asia 32 0.3959254 0.90278 1.167337 5.445392 6...
asia 64 0.2453828 0.4006613 0.8905137 5.641087 5...
asia 128 0.1543742 0.3089822 0.5973851 5.473975 5...
asia 256 0.0771557 0.2283313 0.4139539 5.590241 5...
asia 512 0.03863931 0.1614974 0.3102061 5.576155 5...
asia 1024 0.01975347 0.166499 0.1831748 5.544304 5...
Option prices
The plot below shows the option prices produced for each number of paths, compared to the BlackScholes equivalent (blackdashed line). It is clear that the SobolBrownian bridge method converges the fastest.
q)prxerrplot[r;`prx]
RMSE
We can also plot the RMSE produced by comparing the prices for each method as they converge to the relative BlackScholes price. The expected result is again exhibited, where the SobolBrownian bridge method converges the fastest.
q)prxerrplot[r;`rmse]
Log RMSE
Lastly, we can look at the log RMSE plot as another means of comparison between the methods. Similarly, we see that the SobolBrownian bridge method (blue) exhibits superior performance.
q)prxerrplot[r;`logrsme]
Conclusion
In this paper we demonstrated that it is possible to calculate option prices using both BlackScholes and Monte Carlo/QuasiMonte Carlo methods in q. The Monte Carlo/QuasiMonte Carlo methods deployed different implementations of both Wienerpath approximation and randomnumber generation.
Looking at the results produced, it is clear that both the option price produced and the resulting RMSE/log RMSE converged fastest when compared with the BlackScholes price for the QuasiMonte Carlo approach, with Sobol’ sequence number generation and Brownianbridge construction.
Additionally, by plotting results we have shown that the q implementation replicates the original results produced in C++, presented in the paper by S. Kucherenko et al. 2007.
Author
Deanna Morgan joined First Derivatives in June 2018 as a data scientist in the Capital Markets Training Program and currently works as a machinelearning engineer in London.
Other papers by Deanna Morgan
Acknowledgements
I gratefully acknowledge Sergei Kucherenko for allowing us to create a version of the C++ Option Pricing library in q and for providing technical knowledge throughout the project. I would additionally like to acknowledge my colleagues in the KX Machine Learning team for their guidance in the technical aspects of this paper.