5.0 | 3 ratings Rate this file 58 Downloads (last 30 days) File Size: 8.02 KB File ID: #41473 Version: 1.12
image thumbnail

Fast Matrixwise Black-Scholes Implied Volatility



23 Apr 2013 (Updated )

Calculates Black-Scholes Implied Volatility for Full Surface at High Speed

| Watch this File

File Information

Calculates Black-Scholes Implied Volatility Surface for an Option Price Matrix.
Uses Li's Rational Function Approximator for the Initial Estimate, followed by
3rd-Order Householder's Root Finder (i.e. using vega,vomma & ultima) for greater
convergence rate and wider domain-of-convergence relative to Newton-Raphson. Both
Li's Approximator and the Root Finder are calculated matrix-wise (i.e.
fully vectorized) for increased efficiency.


This file inspired Merton Jump Diffusion Option Price (Matrixwise).

Required Products MATLAB
MATLAB release MATLAB 8.0 (R2012b)
MATLAB Search Path
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (11)
11 May 2015 Mark Whirdy

Hi Tomas

I'd be happy to help but will need to see the values of the input variables to be able to replicate this here.

You can msg me by clicking on my profile, or post here if you like.


Comment only
10 May 2015 Tomas

Tomas (view profile)

Hi Mark,

The function returns this error to me:

Error using \
Matrix dimensions must agree.

Error in calcBSImpVol (line 147)
B = [ones(size(X,1),1),X]\Y;

Any idea what to do about it? It seems that X and Y are empty (as it does not read the definitions). But could not find out what to do about it.


Comment only
02 Apr 2015 Mark Whirdy

Hi Hans

its referenced in the file

"help calcBSImpVol"


Comment only
30 Mar 2015 Hans

Hans (view profile)

Hi, Mark:

Would you be so kind as to provide the detailed reference to Li's 2006 paper on rational function approximation? Thank you very much.

Comment only
02 Oct 2014 Mark Whirdy

Hi Gabriele
1) I haven't come across this actually
2) Your explanation sounds correct
3) I'd be interested in seeing the values for S,K,T,q & P which cause this?

Perhaps you can track down the particular combination which causes it to occur? Use the "Stop on Warning" option (if u have 2012b)

One total guess would be that if you are using Carr-Madan FFT pricing, short-dated deep-OTM expiries can result in negative prices.

Comment only
02 Oct 2014 Gabriele Pompa

Well, in that code I didn't care that much to efficiency since it belonged to a pre-processing step o my data.
But thank for that index trick! I'll use it.

However i'm writing here to ask if you ever experienced a crash in running the code (of course extreme situations). However I did:

Error using erf
Input must be real and full.

Error in calcBSImpVol>fcnN (line 191)


Well, after some backward induction, I came to the conclusion that the problem is in the implementation of the symmetry property of the D domain (eq. (13) and (14) of Li 2006):

v2 = fcnv(p,m,n,i,j,-x,max(exp(x).*c + 1 -exp(x),0)); % Reflection for D+ Domain (x>1)

The fact is that, due to some approximation error, the reflected normalized price c(-x,v) mat become negative, therefore you'll get a complex result in evaluating the sqrt(c) in the rational function "fcnv" definition (eq (19) Li's 2006). Therefore these spurious complex numbers propagate to "v", and then to "sigma" and therefore to both "d1fcn" and "d2fcn", which enter as argument of "fcnN" and finally to the error function "erf".

I fixed this issue naively with a max(.,0), substituting as follows v2:

v2_fixed = fcnv(p,m,n,i,j,-x,max(exp(x).*c + 1 -exp(x),0)); % Reflection for D+ Domain (x>1)

but I was wondering I you can comment:

1) have ever experienced this situation?
2) do you think my explanation is correct?
3) do you think the way I've fixed the problem can be improved (at least in efficiency)

Grazie ;-)

Comment only
30 Sep 2014 Mark Whirdy

Yes, you're right Gabriele, & thanks for your comments. I've made this & other modifications since actually.

You should, I find, avoid using the repmat function and use "indexing-replication" which should shave a few milliseconds off this

q ---> q_mat = repmat(q,1,size(P,2))

q ---> q_mat = q(:,ones(1,h))

Comment only
30 Sep 2014 Gabriele Pompa

Dear Mark, your code was very helpful for me, so thanks a lot!

One little (very little!) possible extension.
With a small modification of your code it is possible to deal also with vectorial (r,q) pairs (vectors of dim. m = number of maturities),
which is a rather realistic case in some pricing problems.

To do so:


1) reshape your [mx1] (r,q) vectors as [mxn] (r_mat, q_mat) matrices (P, is [mxn] matrix of Option Prices):

q ---> q_mat = repmat(q,1,size(P,2))
r ---> r_mat = repmat(r,1,size(P,2))

2) in calcBSImpVol, replace all (r,q) occurrences in lines 148-to-154 with ( r(C), q(C) )

3) use calcBSImpVol with (r_mat, q_mat) as input.


You can easily check it repeating constants (r,q) pairs as to form [mxn] (r_mat, q_mat) matrices:

q_mat = repmat(repmat(q,[size(P,1),1]),1,size(P,2));
r_mat = repmat(repmat(r,[size(P,1),1]),1,size(P,2));

and check whether calcBSImpVol(..., r, q) and calcBSImpVol(..., r_mat, q_mat) return the same sigma matrix of implied volatilities, as they should.

In playing this way, the execution time increased of 10^-2 sec only (well, "only" is my personal point of view).

Comment only
30 Sep 2014 Gabriele Pompa

Very fast and reliable: 10^-1 sec on a 2009 laptop with a 10 maturities x 118 strikes surface with possibly NaN values in the price matrix P (corresponding to non quoted options)

Very accurate: greater domain of convergence than both blsimpv and impvol

04 Oct 2013 Marco

Marco (view profile)

I worked with it for a couple of hours and so far so good. Fast and smooth.

24 Apr 2013 Feargal Deery

simple and clean, this code dramatically speeds up calculating option implied vols. Especially useful for large surfaces. Thanks

24 Apr 2013 1.1

Fixed bug in Put-Call Parity line

P = P + S.*exp(-q.*T) - K.*exp(-r.*T); % Convert Put to Call by Parity Relation

25 Apr 2013 1.6

Added Comments

01 May 2013 1.7

Re-factoring of anonymous functions in Householder root solver to re-calculate derivatives & for only those points which have not converged. This increases speed by 20-50% (depending on particular surface).

27 Jun 2013 1.9

Changed comment
"S Underlying Price [m x n]"
"S Underlying Price [1 x 1]"
.. as was misleading
This version of the file handles S as a scalar only.

28 Apr 2015 1.11

Amended Code allowing matrix inputs for "cp", "r" and "q"

cp Call[+1],Put[-1] [m x n],[1 x 1]
r Continuous Risk-Free Rate [m x n],[1 x 1]
q Continuous Div Yield [m x n],[1 x 1]

28 Apr 2015 1.12

Amended to facilitate input variables "cp", "q" & "r" as matrices
cp Call[+1],Put[-1] [m x n],[1 x 1]
r Continuous Risk-Free Rate [m x n],[1 x 1]
q Continuous Div Yield [m x n],[1 x 1]

Contact us