## European Option Pricing in Julia and MATLAB

October 7th, 2012

I felt like playing with Julia and MATLAB this Sunday morning.  I found some code that prices European Options in MATLAB using Monte Carlo simulations over at computeraidedfinance.com and thought that I’d port this over to Julia.  Here’s the original MATLAB code

function V = bench_CPU_European(numPaths)
%Simple European
steps = 250;
r = (0.05);
sigma = (0.4);
T = (1);
dt = T/(steps);
K = (100);

S = 100 * ones(numPaths,1);

for i=1:steps
rnd = randn(numPaths,1);
S = S .* exp((r-0.5*sigma.^2)*dt + sigma*sqrt(dt)*rnd);
end
V = mean( exp(-r*T)*max(K-S,0) )

I ran this a couple of times to see what results I should be getting and how long it would take for 1 million paths:

tic;bench_CPU_European(1000000);toc
V =
13.1596
Elapsed time is 6.035635 seconds.
>> tic;bench_CPU_European(1000000);toc
V =
13.1258
Elapsed time is 5.924104 seconds.
>> tic;bench_CPU_European(1000000);toc
V =
13.1479
Elapsed time is 5.936475 seconds.

The result varies because this is a stochastic process but we can see that it should be around 13.1 or so and takes around 6 seconds on my laptop. Since it’s Sunday morning, I am feeling lazy and have no intention of considering if this code is optimal or not right now. I’m just going to copy and paste it into a julia file and hack at the syntax until it becomes valid Julia code. The following seems to work

function bench_CPU_European(numPaths)

steps = 250
r = 0.05
sigma = .4;
T = 1;
dt = T/(steps)
K = 100;

S = 100 * ones(numPaths,1);

for i=1:steps
rnd = randn(numPaths,1)
S = S .* exp((r-0.5*sigma.^2)*dt + sigma*sqrt(dt)*rnd)
end
V = mean( exp(-r*T)*max(K-S,0) )
end

I ran this on Julia and got the following

julia> tic();bench_CPU_European(1000000);toc()
elapsed time: 36.259000062942505 seconds
36.259000062942505

julia> bench_CPU_European(1000000)
13.114855104505445

The Julia code appears to be valid, it gives the correct result of 13.1 ish but at 36.25 seconds is around 6 times slower than the MATLAB version.  The dog needs walking so I’m going to think about this another time but comments are welcome.

Update (9pm 7th October 2012):   I’ve just tried this Julia code on the Linux partition of the same laptop and 1 million paths took 14 seconds or so:

tic();bench_CPU_European(1000000);toc()
elapsed time: 14.146281957626343 seconds

I built this version of Julia from source and so it’s at the current bleeding edge (version 0.0.0+98589672.r65a1 Commit 65a1f3dedc (2012-10-07 06:40:18). The code is still slower than the MATLAB version but better than the older Windows build

Update: 13th October 2012

Over on the Julia mailing list, someone posted a faster version of this simulation in Julia

function bench_eu(numPaths)
steps = 250
r = 0.05
sigma = .4;
T = 1;
dt = T/(steps)
K = 100;

S = 100 * ones(numPaths,1);

t1 = (r-0.5*sigma.^2)*dt
t2 = sigma*sqrt(dt)
for i=1:steps
for j=1:numPaths
S[j] .*= exp(t1 + t2*randn())
end
end

V = mean( exp(-r*T)*max(K-S,0) )
end

On the Linux partition of my test machine, this got through 1000000 paths in 8.53 seconds, very close to the MATLAB speed:

julia> tic();bench_eu(1000000);toc()
elapsed time: 8.534484148025513 seconds

It seems that, when using Julia, one needs to unlearn everything you’ve ever learned about vectorisation in MATLAB.

Update: 28th October 2012

Members of the Julia team have been improving the performance of the randn() function used in the above code (see here and here for details).  Using the de-vectorised code above, execution time for 1 million paths in Julia is now down to 7.2 seconds on my machine on Linux.  Still slower than the MATLAB 2012a implementation but it’s getting there.  This was using Julia version  0.0.0+100403134.r0999 Commit 099936aec6 (2012-10-28 05:24:40)

tic();bench_eu(1000000);toc()
elapsed time: 7.223690032958984 seconds
7.223690032958984
• Laptop model: Dell XPS L702X
• CPU:Intel Core i7-2630QM @2Ghz software overclockable to 2.9Ghz. 4 physical cores but total 8 virtual cores due to Hyperthreading.
• RAM: 8 Gb
• OS: Windows 7 Home Premium 64 bit and Ubuntu 12.04
• MATLAB: 2012a
• Julia: Original windows version was Version 0.0.0+94063912.r17f5, Commit 17f50ea4e0 (2012-08-15 22:30:58).  Several versions used on Linux since, see text for details.

## Some Basic Stock Analysis with Mathematica 8

September 21st, 2012

One of my favourite investment news sites is The Motley Fool which frequently run articles such as 10 Shares Trading Near 52 week lows and 15 Shares Trading Near 52 week Highs.  The idea behind such filtering is to seek out shares that have done particularly badly (or well) over the last year and then subject them to further analysis in order to find opportunities.  Thanks to Mathematica’s FinancialData command, it is rather easy to generate these lists yourself whenever you like.

15 Shares Trading Near 52 Week Highs

The original article selected the 15 largest cap shares from the FTSE All Share Index that were trading within 3% of their 52 week high at the time of publication.  Let’s see how to do that using Mathematica.

The following code returns the tickers of all shares from the FTSE All Share Index that are trading within 3% of their 52 week high.

percentage = 3;
all52weekHighs =
Select[FinancialData["^FTAS", "Members"],
Abs[FinancialData[#, "FractionalChangeHigh52Week"]] < (percentage/100.) &];

The variable all52weekHighs contains a list of stock tickers (e.g. LLOY.L) that meet our criteria.  The next thing to do is to find the market cap of each one:

all52WeekHighsWithCaps =
Map[{#, FinancialData[#, "MarketCap"]} &, all52weekHighs];

This works fine for most shares. LloydsTSB for example returns {“LLOY.L”, 2.7746*10^10} at the time of writing but the MarketCap query fails for some tickers. For example, the Market Cap for HSL.L is not available and we get {“HSL.L”, Missing[“NotAvailable”]}.  Let’s discard these by insisting that we only consider stocks that have a numeric market cap.

Goodall52WeekHighsWithCaps =
Select[all52WeekHighsWithCaps, NumberQ[#[[2]]] &];

We sort the list according to MarketCap:

sorted = Sort[Goodall52WeekHighsWithCaps, #1[[2]] > #2[[2]] &];

Let’s prettify the list a little by iterating over all tickers and replacing the ticker with the associated stock name. Also, let’s divide the market cap by 1 million to make it more readable

finallist =
Map[{FinancialData[#[[1]], "Name"], #[[2]]/1000000} &, sorted];

Now, you may be wondering why I haven’t been showing you the output of these commands. This is simply because even this final list is rather large at 118 entries at the time of writing

Length[finallist]

118

The original article only considered the top 15 sorted by Market Cap so let’s show those. Market Caps are given in millions.

top15 = finallist[[1;;15]]//Grid

HSBC Holdings PLC    94159.
National Grid    24432.
Prudential PLC    21775.
Centrica PLC    17193.
Rolls Royce Group    16363.
WPP Plc    10743.
Experian PLC    10197.
Old Mutual PLC    8400.
Legal & General Group PLC    8036.
Wolseley PLC    7955.
Standard Life    6662.
J Sainsbury plc    6401.
Aggreko PLC    6391.
Land Securities Group PLC    6180.
British Land Co PLC    4859.

and we are done.

## Black-Scholes Option Pricing in MATLAB using the NAG Toolbox

July 23rd, 2012

A MATLAB user at Manchester University contacted me recently asking about Black-Scholes option pricing.  The MATLAB Financial Toolbox has a range of functions that can calculate Black-Scholes put and call option prices along with several of the sensitivities (or ‘greeks‘) such as blsprice, blsdelta and so on.

The user’s problem is that we don’t have any site-wide licenses for the Financial Toolbox.  We do, however, have a full site license for the NAG Toolbox for MATLAB which has a nice set of option pricing routines.  Even though they calculate the same things, NAG Toolbox option pricing functions look very different to the Financial Toolbox ones and so I felt that a Rosetta Stone type article might be useful.

For Black-Scholes option pricing, there are three main differences between the two systems:

1. The Financial Toolbox has separate functions for calculating the option price and each greek (e.g. blsprice, blsgamma, blsdelta etc) whereas NAG calculates the price and all greeks simultaneously with a single function call.
2. Where appropriate, The MATLAB functions calculate Put and Call values with one function call whereas with NAG you need to explicitly specify Call or Put.
3. NAG calculates more greeks than MATLAB.

The following code example pretty much says it all.  Any variable calculated with the NAG Toolbox is prefixed NAG_ whereas anything calculated with the financial toolbox is prefixed MW_.  When I developed this, I was using MATLAB 2012a with NAG Toolbox Mark 22.

%Input parameters for both NAG and MATLAB.
Price=50;
Strike=40;
Rate=0.1;
Time=0.25;
Volatility=0.3;
Yield=0;

%calculate all greeks for a put using NAG
[NAG_Put, NAG_PutDelta, NAG_Gamma, NAG_Vega, NAG_PutTheta, NAG_PutRho, NAG_PutCrho, NAG_PutVanna,...
NAG_PutCharm, NAG_PutSpeed, NAG_PutColour, NAG_PutZomma,NAG_PutVomma, ifail] =...
s30ab('p', Strike, Price, Time, Volatility, Rate, Yield);

%calculate all greeks for a Call using NAG
[NAG_Call, NAG_CallDelta, NAG_Gamma, NAG_Vega, NAG_CallTheta, NAG_CallRho, NAG_CallCrho, NAG_CallVanna,...
NAG_CallCharm, NAG_CallSpeed, NAG_CallColour,NAG_CallZomma, NAG_CallVomma, ifail] = ...
s30ab('c', Strike, Price, Time, Volatility, Rate, Yield);

%Calculate the Elasticity (Lambda)
NAG_CallLambda = Price/NAG_Call*NAG_CallDelta;
NAG_PutLambda = Price/NAG_Put*NAG_PutDelta;

%Calculate the same set of prices and greeks using the MATLAB Finance Toolbox
[MW_Call, MW_Put] = blsprice(Price, Strike, Rate, Time, Volatility, Yield);
[MW_CallDelta, MW_PutDelta] = blsdelta(Price, Strike, Rate, Time, Volatility, Yield);
MW_Gamma = blsgamma(Price, Strike, Rate, Time, Volatility, Yield);
MW_Vega = blsvega(Price, Strike, Rate, Time, Volatility, Yield);
[MW_CallTheta, MW_PutTheta] = blstheta(Price, Strike, Rate, Time,Volatility, Yield);
[MW_CallRho, MW_PutRho]= blsrho(Price, Strike, Rate, Time, Volatility,Yield);
[MW_CallLambda,MW_PutLambda]=blslambda(Price, Strike, Rate, Time, Volatility,Yield);

Note that NAG doesn’t output the elasticity (Lambda) directly but it is trivial to obtain it using values that it does output. Also note that as far as I can tell, NAG outputs more Greeks than the Financial Toolbox does.

I’m not going to show the entire output of the above program because there are a lot of numbers.  However, here are the Put values as calculated by NAG shown to 4 decimal places. I have checked and they agree with the Financial Toolbox to within numerical noise.

NAG_Put =0.1350
NAG_PutDelta =-0.0419
NAG_PutLambda =-15.5066
NAG_Gamma =0.0119
NAG_Vega =2.2361
NAG_PutTheta =-1.1187
NAG_PutRho =-0.5572
NAG_PutCrho = -0.5235
NAG_PutVanna =-0.4709
NAG_PutCharm =0.2229
NAG_PutSpeed =-0.0030
NAG_PutColour =-0.0275
NAG_PutZomma =0.0688
NAG_PutVomma =20.3560