Implied Volatility in Matlab

Question

I'm trying to calculate the implied volatility using the Black-Scholes formula in Matlab (2012b), but somehow have problems with some strike prices. For instance blsimpv(1558,1440,0.0024,(1/12),116.4) will return NaN. I thought it probably would be some problem with the function and therefore searched the internet for some other matlab scripts and customized it to my custom needs, but unfortunately I'm still not able to return a valid implied volatility.

function sigma=impvol(C,S,K,r,T)

   %F=S*exp((r).*T);
   %G=C.*exp(r.*T);

   %alpha=log(F./K)./sqrt(T);
   %beta=0.5*sqrt(T);

   %a=beta.*(F+K);
   %b=sqrt(2*pi)*(0.5*(F-K)-G);
   %c=alpha.*(F-K);

   %disc=max(0,b.^2-4*a.*c);
   %sigma0=(-b+sqrt(disc))./(2*a);

   i=-1000;
   while i<=5000
      sigma0=i/1000;
      sigma=NewtonMethod(sigma0);
      if sigma<=10 && sigma>=-10
          fprintf('This is sigma %f',sigma)
      end
      i=i+1;
    end
end

function s1=NewtonMethod(s0)

    s1=s0;
    count=0;
    f=@(x) call(S,K,r,x,T)-C;
    fprime=@(x) call_vega(S,K,r,x,T);

    max_count=1e4;

    while max(abs(f(s1)))>1e-7 && count<max_count
        count=count+1;
        s0=s1;
        s1=s0-f(s0)./fprime(s0);

    end

end

end

function d=d1(S,K,r,sigma,T)
   d=(log(S./K)+(r+sigma.^2*0.5).*(T))./(sigma.*sqrt(T));
end

function d=d2(S,K,r,sigma,T)
   d=(log(S./K)+(r-sigma.^2*0.5).*(T))./(sigma.*sqrt(T));
end

function p=Phi(x)
  p=0.5*(1.+erf(x/sqrt(2)));
end

function p=PhiPrime(x)
   p=exp(-0.5*x.^2)/sqrt(2*pi);
end

function c=call(S,K,r,sigma,T)
  c=S.*Phi(d1(S,K,r,sigma,T))-K.*exp(-r.*(T)).*Phi(d2(S,K,r,sigma,T));
end

function v=call_vega(S,K,r,sigma,T)
   v=S.*PhiPrime(d1(S,K,r,sigma,T)).*sqrt(T);
end

Running impvol(116.4,1558,1440,0.0024,(1/12)) will however unfortunately return the value 'Inf'. There somehow is a problem with the Newton-Rhapson method not converging but I am kind of clueless how to solve this. Does anyone know how to solve this problem or know some other way how to calculate the implied volatility?

Thank u in advance already for your help! Kind regards,

Henk

Answer 1

I would definitely suggest this code: Fast Matrixwise Black-Scholes Implied Volatility It is able to compute the entire surface in one shot and - my experience - I found it much more reliable than blsimpv() or impvol() which are other functions implemented in matlab.

Answer 2

Newton-Rhapson method does not work well for implied volatility. You should use the bisection method (not sure how it is used in Matlab). It is described in http://en.wikipedia.org/wiki/Bisection_method. For completeness, it works this way:

1) Pick an arbitrary high (impossible) volatility like high=200%/year.

2) Pick lowest possible volatility (low=0%).

2a) Calculate option premium for 0% volatility, if actual premium is lower than that, it means negative volatility (which is "impossible").

3) While implied volatility is not found:

3.1) If "high" and "low" are very near (e.g. equal up to 5th decimal), either one is your implied volatility. If not...

3.2) Calculate average between "high" and "low". avg=(high+low)/2

3.3) Calculate option premium for avg volatility.

3.4) If actual premium is higher then p(avg), make min=avg, because implied volatility must lie between avg and max.

3.4a) If actual premium is lower than p(avg), make max=avg, because implied vol must lie between min and avg.

The main disvantage of bisect is that you have to pick a maximum value, so your function won't find implied volatilities bigger than that. But something like 200%/year should be high enough for real-world usage.

I use yet another method more like Newton's method, hence not limited to a range, since vega is a derivative, but with a "linearization" fix to avoid hunting and failure due to small vegas:

def implied_volatility(type, premium, S, K, r, s_dummy, t):
    if S <= 0.000001 or K <= 0.000001 or t <= 0.000001 or premium <= 0.000001:
        return 0.0

    s = 0.35

    for cycle in range(0, 120):
        ext_premium = type(S, K, r, s, t)
        if abs(premium - ext_premium) < 0.005:
            return s
        ext_vega = type.vega(S, K, r, s, t)
        # print S, K, r, s, t, premium, ext_premium, ext_vega
        if ext_vega < 0.0000001:
            # Avoids zero division if stuck 
            ext_vega = 0.0000001

        s_new = s - (ext_premium - premium) / ext_vega
        if abs(s_new - s) > 1.00:
            # estimated s is too different from previous;
            # it is better to go linearly, since
            # vega is too small to give a good guess
            if s_new > s:
                s += 1.0
            else:
                s -= 1.0
        else:
            s = s_new

        if s < 0.0:
            # No volatility < 0%
            s = 0.0001
        if s > 99.99:
            # No point calculating volatilities > 9999%/year
            return 100.0

    return 0.0

Still, I think that bisect is your best bet.

Answer 3

I created a simple function that conducts a sort of trial and error calculation if the output from blsimpv is NaN. This slows down the computation time significantly for me but it always gives me a desirable result.

The function is shown to be used below

BSIVC(t,i)= blsimpv(S(t,i),K,r,tau(t),HestonCiter(t,i))
 if isnan(BSIVC(t,i));
  BSIVC(t,i)= secondIVcalc(HestonCiter(t,i),S(t,i),K,r,q,tau(t))
 end

The function itself is described below:

function IV= secondIVcalc(HestonC,S,K,r,q,T)
lowCdif = 1;
a=0;
while lowCdif>0.0001
a= a+0.00001
lowCdif = HestonC - BSCprice(S,K,r,q,a,T);
end
IV= a;
end

Please note that BSCprice is not an in-built function in matlab.

Just to make the code clearer- BSCprice is of the format BSCprice(Underlying Asset Price, Strike Price, interest rate, dividend yield, implied vol, time to maturity).