Calculates BlackScholes Implied Volatility for Full Surface at High Speed
Calculates BlackScholes Implied Volatility Surface for an Option Price Matrix.
Uses Li's Rational Function Approximator for the Initial Estimate, followed by
3rdOrder Householder's Root Finder (i.e. using vega,vomma & ultima) for greater
convergence rate and wider domainofconvergence relative to NewtonRaphson. Both
Li's Approximator and the Root Finder are calculated matrixwise (i.e.
fully vectorized) for increased efficiency.
1.12  Amended to facilitate input variables "cp", "q" & "r" as matrices


1.11  Amended Code allowing matrix inputs for "cp", "r" and "q" cp Call[+1],Put[1] [m x n],[1 x 1]


1.9  Changed comment


1.7  Refactoring of anonymous functions in Householder root solver to recalculate derivatives & for only those points which have not converged. This increases speed by 2050% (depending on particular surface). 

1.6  Added Comments 

1.1  Fixed bug in PutCall Parity line P = P + S.*exp(q.*T)  K.*exp(r.*T); % Convert Put to Call by Parity Relation 
JFz (view profile)
JFz (view profile)
Hi, Mark,
Thanks for this file. Looks excellent.
But when I test it with a simple example,
cp = 1; S = 100; K = 200; T = 1; r = 0; q = 0; p = 3;
it gives me this error:
147 B = [ones(size(X,1),1),X]\Y;
Error using \
Matrix dimensions must agree.
How to resolve this issue?
Liqun (view profile)
Good work, Mark!
To the new user, the two issues mentioned below still remain as of 12/11/2016, which means that Mark has not update the code yet. You should take care of that bu yourself after download the code.
Ravi Teja (view profile)
Henrik Dam (view profile)
I have been fiddling with the new Jaeckel, but am nowhere near your speed (I probably can not write as efficiently).
Have you seen:
http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1027282
although the idea seems to be completely different.
Finally have you considered simply using (allowing) the IV you are calibrating to as initial value out of the range of the approximation?
best regards,
Henrik
Mark Whirdy (view profile)
Hi Henrik
Li's initial approximation serves only to reduce the number of householder iterations really  which is ultimately what I'm trying to achieve as, especially where the function is used within an objection function for calibration of a vol surface.
The bivar regression for outofLidomain values was the bestbad solution I had (having scoured the literature for approximation schemes for deep OTM/ITM strikes). If you come across better please let me know. Currently, I have no better solution (other than arbitrarily choosing say 0.2 and letting householder do the legwork).
A more robust solution might involve repeating Li's research over a wider domain of moneyness to derive a rational approximation in many more terms (very time consuming), or else having a look at ... http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1027282
... or ...
http://www.pjaeckel.webspace.virginmedia.com/LetsBeRational.pdf
All the best
Mark
Henrik Dam (view profile)
Hi Mark.
Thanks for your fast response.
Ah okay. I also have other data, so that probably explains why it sometimes work. So basically the algorithm has a region where it works and you use information and regression to deduce values out of that region? I haven't read the article carefully, but is that suggested or your idea? It is very interesting. What do you intend to replace it with? (I'm just curious).
best regards,
Henrik
Mark Whirdy (view profile)
Hi Henrik
It is due to having no indomain (with respect to Li's approximator) values to use as regressors for estimating the sigma for outofdomain values.
Put simply, you need to supply some neartothemoney values to allow it calculate farfromthemoney values.
Using a multivariate regression for outofdomain values is a pretty brutal approach here anyway, and I will update the file as soon as I can.
In the interim, just replace the regression section with
sigma = 0.2*ones(size(sigma));
then you'll get
sigma =
0.188411203301236
Henrik Dam (view profile)
Hi Mark.
I have the exact same problem as Tomas I will have look my self as well to see if I can find an error. This is actual (rounded) observed data.
Call price: 2.2 *10^(6)
S0: 1100
K: 1450
T= 1/12
r = 0.01
q=0.02
Matlab's own calculator yields an implied vol of: 0.188.
Thank you for your very nice code and your time.
Mark Whirdy (view profile)
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.
Best
Mark
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.
Best,
T.
Mark Whirdy (view profile)
Hi Hans
its referenced in the file
type
"help calcBSImpVol"
http://papers.ssrn.com/sol3/papers.cfm?abstract_id=952727
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.
Mark Whirdy (view profile)
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 CarrMadan FFT pricing, shortdated deepOTM expiries can result in negative prices.
Gabriele Pompa (view profile)
Well, in that code I didn't care that much to efficiency since it belonged to a preprocessing 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)
p=0.5*(1.+erf(x./sqrt(2)));

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 ;)
Mark Whirdy (view profile)
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 "indexingreplication" which should shave a few milliseconds off this
q > q_mat = repmat(q,1,size(P,2))
q > q_mat = q(:,ones(1,h))
Gabriele Pompa (view profile)
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 148to154 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).
Gabriele Pompa (view profile)
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
Marco (view profile)
I worked with it for a couple of hours and so far so good. Fast and smooth.
Feargal Deery (view profile)
simple and clean, this code dramatically speeds up calculating option implied vols. Especially useful for large surfaces. Thanks