How to get Implied Volatility?

In this article, I will introduce what is implied volatility and several methods to find it. Here are the points I will try to cover:

  • What is Implied Volatility? 
  • Dichotomy Method
  • Newton Raphson Method
  • Example in Python with a set of option prices
  • Models
  • Conclusion

Implied Volatility

Historical volatility and implied volatility, what is the difference?
Historical volatility is merely the standard deviation of a financial instrument’s returns over a given period. It is measured on past data. Implied volatility is the volatility that makes the theoretical price of your option equal to the market price. This implied volatility changes with the current price of the option and then reflects the market estimation for the future fluctuations of the underlying.

Implied volatility isn’t the same for every quoted strike. If you plot these implied volatilities against strikes you will get the volatility smile. If you plot implied volatilities against strike for a set of maturity, you will get the surface of volatility.

An example to see what it looks like:

Volatility Surface of option on Arcelor Mittal – {OVDV} on Bloomberg

How can I get this volatility smile with a set of option price?

Dichotomy Method

The dichotomy method (a.k.a bisection method or binary search method) is a numerical method (simple, robust but slow) used to find the solution of an equation with a single unknown. It takes in input a continuous function in which the root has to be found, a set [a,b] that contains the root and the desired accuracy.

Let’s find the implied volatility for an European Vanilla Option on an equity. Thus, we want to find \sigma_{imp} so that f(\sigma_{imp}) = P_{BS}(\sigma_{imp})-P_{market} = 0.

Algorithms steps

Binary Search Method – Image from wikipedia

  • Step 1 – Compute M = \frac{a+b}{2}
  • Step 2 – Compute Y = f(a)*f(M)
  • Step 3 – If Y<0 then b=M, otherwise a=M.
  • Step 4 – Repeat steps 1 to 3 until either the accuracy or the maximum number of iteration is reached.

Newton-Raphson Method

This method is less time-consuming but isn’t derivative free as the binary search method we have seen previously. With this algorithm, you need the derivative of the option price with respect to the volatility. It works well for options where the Black & Scholes model has a closed-form solution for both the price and vega (harder for American and exotic options). It takes the function f, its derivative, an initial guess of the root x_0 and the desired accuracy.

The Taylor Series of a function f at x = x_0 to the first order is:

f_{x_0}(x) \approx f(x_0) + f'(x_0)*(x-x_0)

We are looking for x so that f(x) = 0. Then we want:

0 \approx f(x_0) + f'(x_0)*(x-x_0)

x = x_0 - \frac{f(x_0)}{f'(x_0)}

As before, we want to find the implied volatility so that f(\sigma_{imp}) = 0. The schema can be written as follow :

\sigma_{n+1} = \sigma_{n} - \frac{f(\sigma_n)}{f'(\sigma_n)}

with:

f'(x) = \nu = S*\sqrt{T}*n(d_1)

Algorithm steps

Newton Raphson Algorithm – By Ralph Pfeifer

  • Step 1 – Compute x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}
  • Step 2 – Repeat step 1 until the precision either the precision or the maximum number of iterations is reached.

This method sometimes fails to converge to the root, as we will see in the following implementation.

Example of implementation

I took a dataset of option prices on Apple’s stock maturing in a month. I implemented both the Dichotomy and Newton-Raphson methods in Python, applied them to the dataset and compared the results. The code is here:

import pandas as pd
import plotly as plt
import plotly.graph_objs as go
import mpmath as math
from scipy.stats import norm
import datetime as dt

class BlackScholes():
    @staticmethod
    def __d1(S,R,K,T,Vol):
        return float((math.ln(S/K)+(R+(Vol**2)/2)*T)/(Vol*math.sqrt(T)))

    @staticmethod
    def __d2(S,R,K,T,Vol):
        return float((math.ln(S/K)+(R-(Vol**2)/2)*T)/(Vol*math.sqrt(T)))

    @staticmethod
    def call_price(S,R,K,T,Vol):
        return S*norm.cdf(BlackScholes.__d1(S,R,K,T,Vol))-K*math.exp(-R*T)*norm.cdf(BlackScholes.__d2(S,R,K,T,Vol))

    @staticmethod
    def put_price(S,R,K,T,Vol):
        return -S*norm.cdf(-BlackScholes.__d1(S,R,K,T,Vol))+K*math.exp(-R*T)*norm.cdf(-BlackScholes.__d2(S,R,K,T,Vol))

    @staticmethod
    def vega(S,R,K,T,Vol):
        return float(S*math.sqrt(T)*norm.cdf(BlackScholes.__d1(S,R,K,T,Vol)))

class Implied_Vol():
    @staticmethod
    def get_vol(data,S,T,R,nb_steps_max=5000,precision=0.0001):
        data["IV_bisection"]=0
        data["IV_newton"]=0
        for i in range(data.shape[0]):
            data.loc[i,"IV_bisection"]=Implied_Vol.__bisection_method(data.loc[i,"Last_Price"],
                                                                       S,
                                                                       data.loc[i,"Strike"],
                                                                       R,
                                                                       T,
                                                                       data.loc[i,"Type"],
                                                                       nb_steps_max,
                                                                       precision)

            data.loc[i,"IV_newton"]=Implied_Vol.__newton_raphson_method(data.loc[i,"Last_Price"],
                                                                       S,
                                                                       data.loc[i,"Strike"],
                                                                       R,
                                                                       T,
                                                                       data.loc[i,"Type"],
                                                                       nb_steps_max,
                                                                       precision)
        return data

    @staticmethod
    def __bisection_method(Market_Price,S,K,R,T,type,nb_steps_max,precision):
        #lower bound and upper bound defintion :
        a = 0.001
        b = 1
        m= (b+a)/2
        count_iter = 1

        if type=="Call":
            while count_iter<=nb_steps_max and math.fabs(b-a)>precision:
                #If f(a)*f(m)<0:
                if (BlackScholes.call_price(S,R,K,T,a)-Market_Price)*(BlackScholes.call_price(S,R,K,T,m)-Market_Price)<0:
                    b=m
                else:
                    a=m
                m = (b + a) / 2
                count_iter+=1
        else:
            while (count_iter <= nb_steps_max and math.fabs(b - a) > precision):
                # If f(a)*f(m)<0:
                if (BlackScholes.put_price(S, R, K, T, a) - Market_Price) * (BlackScholes.put_price(S, R, K, T, m) - Market_Price) < 0:
                    b = m
                else:
                    a = m
                m = (b + a) / 2
                count_iter += 1
        return m

    @staticmethod
    def __newton_raphson_method(Market_Price,S,K,R,T,type,nb_steps_max,precision):
        x0 = 0.2
        x = 0
        count_iter = 1
        if type=="Call":
            while count_iter<=nb_steps_max and math.fabs(x-x0)>precision:
                x = x0
                x0 = x0 - (BlackScholes.call_price(S,R,K,T,x0)-Market_Price)/BlackScholes.vega(S,R,K,T,x0)
                count_iter += 1
        else:
            while count_iter<=nb_steps_max and math.fabs(x-x0)>precision:
                x = x0
                x0 = x0 - (BlackScholes.put_price(S,R,K,T,x0)-Market_Price)/BlackScholes.vega(S,R,K,T,x0)
                count_iter += 1
        return float(x0)

if __name__ == "__main__":

    interest_rate = 0.01
    spot = 169.37
    maturity = dt.datetime.strptime("12/01/2018", "%d/%m/%Y")
    maturity_in_years = (maturity - dt.datetime.now()).days / 365

    #Importing the data
    data = pd.read_csv("AAPL_Options_Quotes.csv", delimiter=';')

    #Computing the vol
    data = Implied_Vol.get_vol(data, spot, maturity_in_years, interest_rate, nb_steps_max=100, precision=0.00001)

    #Sorting the data by Strike
    data = data.sort_values("Strike")

    scatter = go.Scatter(
        x=data["Strike"],
        y=data["IV_bisection"],
        name="Bisection Method"
    )

    scatter2 = go.Scatter(
        x=data["Strike"],
        y=data["IV_newton"],
        name="Newton-Raphson Method"
    )

    Layout = go.Layout(title="AAPL Volatility Smile 12/01/2018 / My Financial Markets", xaxis=dict(title="Strike"),
                       yaxis=dict(title="Implied Volatility"))

    #Plotting our implied volatility smiles
    plt.offline.plot({'data': [scatter, scatter2], 'layout': Layout})

Here is the final result:

AAPL volatility smile computed with two methods

We can see that the Newton-Raphson method fails to compute the two last implied volatility for OTM calls. It’s a case a non-convergence, the algorithm stops itself as the number of iterations reaches its maximum. The following chart shows a case where the algorithm doesn’t converge, this is helpful to visualize how it’s possible.

A case of non-convergence for the Newton-Raphson algorithm – From Math.StackExchange

Models

Black & Scholes model works with a flat volatility and therefore doesn’t reflect the reality of the market. There are two mains ways to model the volatility smile: local volatility and stochastic volatility.

In local volatility models (such as Dupire formula), volatility is a function of the time and the current value of the underlying.

In stochastic volatility models (such as Heston model or SABR), volatility is considered as a random process. Heston model considers that it follows a mean-reverting process.

We will review these models in an upcoming article!

Conclusion

We have seen two algorithms with different performance and requirements that can be used to find implied volatility from option prices. Newton Raphson algorithm needs a closed-form for the derivative but it can be approximated by a finite-difference method (this method is also known as the Secant Method).

The volatility smile coming from Apple’s options has the typical shape of equity derivatives smile (skewed). As we have seen in this article, OTM puts have higher implied volatilities than OTM calls. It comes from a fear of a crash by investors. Investors are globally long in stocks and therefore buy OTM puts as a protection (higher demand makes higher price and therefore higher implied volatilities).

Thank you!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

w

Connecting to %s