Code to Calculate Volatility Smile

Volatility smiles are created by implied volatility changing as the underlying asset moves more ITM or OTM. The more an option is ITM or OTM, the greater its implied volatility becomes. Implied volatility tends to be lowest with ATM options.

Packages used:

import pandas as pd
import numpy as np
from math import sqrt, exp, log, pi
from scipy.stats import norm
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import interpolate

Data is in the form of Days to maturity (t), Strike (K), Spot (S), and Option Price (C)

A Volatility Column is added!

Functions to calculate the values of d1 and d2 as well as the call price.

To extend to puts, one could just add a function that calculates the put price, or combine calls and puts into a single function that takes an argument specifying which type of contract one is dealing with.

def d(sigma, S, K, r, q, t):

    d1 = 1 / (sigma * sqrt(t)) * ( log(S/K) + (r - q + sigma**2/2) * t)

    d2 = d1 - sigma * sqrt(t)

    return d1, d2

def call_price(sigma, S, K, r, q, t, d1, d2):

    C = norm.cdf(d1) * S * exp(-q * t)- norm.cdf(d2) * K * exp(-r * t)

    return C

def put_price(sigma, S, K, r, q, t, d1, d2):

    P = norm.cdf(-d2) * K * exp(-(r) * t) - norm.cdf(-d1) * S * exp(-(q) * t)
    
    return P

The Code that calculates Volatility using Numerical Methods assuming the solution converges

r = 0
q = 0.02
t = 357/365
S = 2941.76
call_strikes_list = []

call_implied_vol = []

no_of_iterations = []

for i in range(len(call_df)):
    
    K = call_df['Strike'].iloc[i]
    
    C0 = call_df['Option price'].iloc[i]

    call_strikes_list.append(K)

    #  We need to provide an initial guess for the root of our function
    
    vol = 0.5

    
    epsilon = 1.0   #  Define variable to check stopping conditions
    
    abstol = 0.0001 #  Stop calculation when abs(epsilon) < this number
    
    count = 0               #  Variable to count number of iterations
    
    max_iter = 1000         #  Max number of iterations before aborting

    while epsilon > abstol:
        #  if-statement to avoid getting stuck in an infinite loop.
        if count > max_iter:
            print ('Program failed to find a root.  Exiting.')
            break

        count += 1
        
        orig_vol = vol
        
        d1, d2 = d(vol, S, K, r, q, t)
        
        function_value = call_price(vol, S, K, r, q, t, d1, d2) - C0
        
        vega = S * norm.pdf(d1) * sqrt(t)* exp(-q * t)
        
        vol = - function_value/vega + vol
        
        epsilon = abs(function_value)
    
    call_implied_vol.append(vol)
    
    no_of_iterations.append(i)

Here we get the Volatility column which we plot by the following code in our call_df

call_zip = list(zip(call_strikes_list, call_implied_vol))
plt.plot(*zip(*call_zip))
plt.show()

Similarly the code for put

r = 0
q = 0.02
t = 357/365
S = 2941.76
put_strikes_list = []

put_implied_vol = []

no_of_iterations = []

for i in range(len(put_df)):
    
    K = put_df['Strike'].iloc[i]
    
    P0 = put_df['Option price'].iloc[i]

    put_strikes_list.append(K)

    
    if  i > 80:     
        vol = 0.05
    else:
        vol = 0.5
    
    epsilon = 1.0          #  Define variable to check stopping conditions
    
    abstol = 0.0001          #  Stop calculation when abs(epsilon) < this number
    
    count = 0                  #  Variable to count number of iterations
    
    max_iter = 1000         #  Max number of iterations before aborting

    while epsilon > abstol:
        #  if-statement to avoid getting stuck in an infinite loop.
        if i > max_iter:
            print ('Program failed to find a root.  Exiting.')
            break

        count += 1
        
        orig_vol = vol
        
        d1, d2 = d(vol, S, K, r, q, t)
        
        function_value = put_price(vol, S, K, r, q, t, d1, d2) - P0
        
        vega = S * norm.pdf(d1) * sqrt(t)* exp(-q * t)
        
        vol = -function_value/vega + vol
        
        epsilon = abs(function_value)
    
    put_implied_vol.append(vol)
    
    no_of_iterations.append(i)
    
put_zip = list(zip(put_strikes_list, put_implied_vol))
plt.plot(*zip(*put_zip))
plt.show()