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:

*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 so that .

*Algorithms steps*

- Step 1 – Compute
- Step 2 – Compute
- 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 and the desired accuracy.

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

We are looking for so that . Then we want:

As before, we want to find the implied volatility so that . The schema can be written as follow :

with:

*Algorithm steps*

- Step 1 – Compute
- 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:

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.

*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).