Surface - suhasghorp.com

Option Volatility and Binomial Model

August 4, 2017January 7, 2021 suhasghorp Leave a comment

In my previous post Options and Volatility Smile , we used Black-Scholes formula to derive Implied Volatility from given option strike and tenor. Most of the options traded on exchanges are American (with a few index options being European) and can be exercised at any time prior to expiration. Whether it is optimal to exercise an option early depends on whether the stock pays dividend or the level of interest rates and is a very complex subject. What I want to focus on is using Binomial model to price an American option. I summarize below Binomial model theory from the excellent book Derivative Markets 3rd Ed. by Robert McDonald

There are many flavors of the Binomial model but they all have following steps in common:

Simulate future prices of underlying stock at various points in time until expiry
Calculate the payoff of the option at expiry
Discount the payoff back to today to get the option price

1. Simulate future prices of the underlying

Assume the price of a stock is $S_0$ and over a small time period $\Delta t$, the price could either be $S_0u$ or $S_0d$ where u is a factor by which price rises and d is a factor by which price drops. The stock is assumed to follow a random walk, also assume that $p$ is the probability of the stock price rising and $(1-p)$ is the probability of it falling. There are many ways to approach the values of $u$,$d$ and $p$ and the various Binomial models differ in the ways these three parameters are calculated.

In Cox-Ross-Rubinstein (CRR) model, $u = \frac{1}{d}$ is assumed. Since we have 3 unknowns, 2 more equations are needed and those come from risk neutral pricing assumption. Over a small $\Delta t$ the expected return of the stock is

$$pu + (1-p)d = e^{r \Delta t}$$

and the expected variance of the returns is

$$pu^2 + (1-p)d^2 – (e^{r \Delta t})^2 = \sigma ^2 \Delta t$$

Solving for $u$, $d$ and $p$, we get

$$p = \frac{e^{r\Delta t} – d}{u-d}$$

$$u = e^{\sigma \sqrt{\Delta t}}$$

$$d = e^{-\sigma\sqrt{\Delta t}}$$

The CRR model generates a following tree as we simulate multi step stock price movements, this a recombining tree centered around $S_0$.

2. Calculating payoffs at expiry

In this step, we calculate the payoff at each node that corresponds to expiry.

For a put option, $payoff = max(K – S_N, 0)$

For a call option, $payoff = max(S_N – K, 0)$

where $N$ is node at expiry with a stock price $S_N$ and $K$ is the strike.

3. Discounting the payoffs

In this step, we discount the payoffs at expiry back to today using backward induction where we start at expiry node and step backwards through time calculating option value at each node of the tree.

For American put, $V_n = max(K – S_n, e^{-r \Delta t} (p V_u + (1-p) V_d))$

For American call, $V_n = max(S_n – K, e^{-r \Delta t} (p V_u + (1-p) V_d))$

$V_n$ is the option value at node n

$S_n$ is the stock price at node n

$r$ is risk free interest rate

$\Delta t$ is time step

$V_u$ is the option value from the upper node at n+1

$V_d$ is the option value from the lower node at n+1

All the variants of Binomial model, including CRR, converge to Black-Scholes in the limit $\Delta t \to 0$ but the rate of convergence is different. The variant of Binomial model that I would like to use is called Leisen-Reimer Model which converges much faster. Please see the original paper for formulas and a C++ implementation at Volopta.com which I have ported to Python in the next section.

The python code is going to look very similar to Options and Volatility Smile post except that we will swap out Black-Scholes framework with Leisen-Reimer model. We will also use the same option chain data AAPL_BBG_vols.csv

A note on Python code

I usually do not write code like below, I am purposely avoiding using any classes so that the focus remains on the objective which is to understand the technique.

import math
import numpy as np
import matplotlib as mpl
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.optimize import brentq
from scipy.interpolate import interp2d
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

def BinomialLeisenReimer(S0, K, r, T, N, sigma, div):

    is_call = False

    odd_N = N if (N % 2 == 1) else (N + 1)

    d1 = (math.log(S0 / K) + ((r - div) + (sigma ** 2) / 2.) * T) / (sigma * math.sqrt(T))

    d2 = (math.log(S0 / K) + ((r - div) - (sigma ** 2) / 2.) * T) / (sigma * math.sqrt(T))

    dt = T / odd_N
    r_n = math.exp((r - div) * dt)
    method = 2.
    Term1 = math.pow((d1 / (N + 1.0 / 3.0 - (1 - method) * 0.1 / (N + 1))), 2) * (N + 1.0 / 6.0)
    pp = 0.5 + math.copysign(1, d1) * 0.5 * math.sqrt(1 - math.exp(-Term1))

    Term2 = math.pow((d2 / (N + 1.0 / 3.0 - (1 - method) * 0.1 / (N + 1))), 2) * (N + 1.0 / 6.0)
    p = 0.5 + math.copysign(1, d2) * 0.5 * math.sqrt(1 - math.exp(-Term2))

    u = r_n * pp / p
    d = (r_n - p * u) / (1 - p)

    qu = p
    qd = 1 - p

    #Initialize price tree

    STs = [np.array([S0])]
    # Simulate the possible stock prices path
    for i in range(N):
        prev_branches = STs[-1]
        st = np.concatenate((prev_branches * u, [prev_branches[-1] * d]))
        STs.append(st)  # Add nodes at each time step

    #Initialize payoff tree
    payoffs = np.maximum(0, (STs[N]-K) if is_call else (K-STs[N]))

    def check_early_exercise(payoffs, node):
        early_ex_payoff = (STs[node] - K) if is_call else (K - STs[node])
        return np.maximum(payoffs, early_ex_payoff)

    #Begin tree traversal
    for i in reversed(range(N)):
        # The payoffs from NOT exercising the option
        payoffs = (payoffs[:-1] * qu + payoffs[1:] * qd) * df
        # Payoffs from exercising, for American options
        early_ex_payoff = (STs[i] - K) if is_call else (K - STs[i])
        payoffs = np.maximum(payoffs, early_ex_payoff)

    return payoffs[0]

headers = ['Date','Strike','Call_Bid','Call_Ask','Call','Maturity','Put_Bid','Put_Ask','Put']
dtypes = {'Date': 'str', 'Strike': 'float', 'Call_Bid': 'float', 'Call_Ask': 'float',
          'Call':'float','Maturity':'str','Put_Bid':'float','Put_Ask':'float','Put':'float'}
parse_dates = ['Date', 'Maturity']
data = pd.read_csv('AAPL_BBG_vols.csv',skiprows=1,header=None, names=headers, dtype=dtypes, parse_dates=parse_dates)
data['LR_Imp_Vol'] = 0.0
data['TTM'] = 0.0
r = 0.0248 #risk-free rate
S0 = 153.97 # spot price as of 5/11/2017
div = 0.0182
for i in data.index:
    t = data['Date'][i]
    T = data['Maturity'][i]
    K = data['Strike'][i]
    Put = data['Put'][i]
    sigma_guess = 0.2
    time_to_maturity = (T - t).days/365.0
    data.loc[i, 'TTM'] = time_to_maturity
    def func_LR(sigma):
        return BinomialLeisenReimer(S0, K, r, time_to_maturity, ((T - t).days) - 1, sigma, div) - Put

    lr_imp_vol = brentq(func_LR, 0.03, 1.0)
    data.loc[i, 'LR_Imp_Vol'] = lr_imp_vol

ttm = data['TTM'].tolist()
strikes = data['Strike'].tolist()
imp_vol = data['LR_Imp_Vol'].tolist()
f = interp2d(strikes,ttm,imp_vol, kind='cubic')

plot_strikes = np.linspace(data['Strike'].min(), data['Strike'].max(),25)
plot_ttm = np.linspace(0, data['TTM'].max(), 25)
fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(plot_strikes, plot_ttm)
Z = np.array([f(x,y) for xr, yr in zip(X, Y) for x, y in zip(xr,yr) ]).reshape(len(X), len(X[0]))
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0.1)
ax.set_xlabel('Strikes')
ax.set_ylabel('Time-to-Maturity')
ax.set_zlabel('Implied Volatility')
fig.colorbar(surf, shrink=0.5, aspect=5)

100

101

102

import math

import numpy as np

import matplotlib as mpl

import pandas as pd

import matplotlib.pyplot as plt

from scipy.integrate import quad

from scipy.optimize import brentq

from scipy.interpolate import interp2d

from mpl_toolkits.mplot3d import Axes3D

from matplotlib import cm

def BinomialLeisenReimer(S0, K, r, T, N, sigma, div):

is_call = False

odd_N = N if (N % 2 == 1) else (N + 1)

d1 = (math.log(S0 / K) + ((r - div) + (sigma ** 2) / 2.) * T) / (sigma * math.sqrt(T))

d2 = (math.log(S0 / K) + ((r - div) - (sigma ** 2) / 2.) * T) / (sigma * math.sqrt(T))

dt = T / odd_N

r_n = math.exp((r - div) * dt)

method = 2.

Term1 = math.pow((d1 / (N + 1.0 / 3.0 - (1 - method) * 0.1 / (N + 1))), 2) * (N + 1.0 / 6.0)

pp = 0.5 + math.copysign(1, d1) * 0.5 * math.sqrt(1 - math.exp(-Term1))

Term2 = math.pow((d2 / (N + 1.0 / 3.0 - (1 - method) * 0.1 / (N + 1))), 2) * (N + 1.0 / 6.0)

p = 0.5 + math.copysign(1, d2) * 0.5 * math.sqrt(1 - math.exp(-Term2))

u = r_n * pp / p

d = (r_n - p * u) / (1 - p)

qu = p

qd = 1 - p

#Initialize price tree

STs = [np.array([S0])]

# Simulate the possible stock prices path

for i in range(N):

prev_branches = STs[-1]

st = np.concatenate((prev_branches * u, [prev_branches[-1] * d]))

STs.append(st) # Add nodes at each time step

#Initialize payoff tree

payoffs = np.maximum(0, (STs[N]-K) if is_call else (K-STs[N]))

def check_early_exercise(payoffs, node):

early_ex_payoff = (STs[node] - K) if is_call else (K - STs[node])

return np.maximum(payoffs, early_ex_payoff)

#Begin tree traversal

for i in reversed(range(N)):

# The payoffs from NOT exercising the option

payoffs = (payoffs[:-1] * qu + payoffs[1:] * qd) * df

# Payoffs from exercising, for American options

early_ex_payoff = (STs[i] - K) if is_call else (K - STs[i])

payoffs = np.maximum(payoffs, early_ex_payoff)

return payoffs[0]

headers = ['Date','Strike','Call_Bid','Call_Ask','Call','Maturity','Put_Bid','Put_Ask','Put']

dtypes = {'Date': 'str', 'Strike': 'float', 'Call_Bid': 'float', 'Call_Ask': 'float',

'Call':'float','Maturity':'str','Put_Bid':'float','Put_Ask':'float','Put':'float'}

parse_dates = ['Date', 'Maturity']

data = pd.read_csv('AAPL_BBG_vols.csv',skiprows=1,header=None, names=headers, dtype=dtypes, parse_dates=parse_dates)

data['LR_Imp_Vol'] = 0.0

data['TTM'] = 0.0

r = 0.0248 #risk-free rate

S0 = 153.97 # spot price as of 5/11/2017

div = 0.0182

for i in data.index:

t = data['Date'][i]

T = data['Maturity'][i]

K = data['Strike'][i]

Put = data['Put'][i]

sigma_guess = 0.2

time_to_maturity = (T - t).days/365.0

data.loc[i, 'TTM'] = time_to_maturity

def func_LR(sigma):

return BinomialLeisenReimer(S0, K, r, time_to_maturity, ((T - t).days) - 1, sigma, div) - Put

lr_imp_vol = brentq(func_LR, 0.03, 1.0)

data.loc[i, 'LR_Imp_Vol'] = lr_imp_vol

ttm = data['TTM'].tolist()

strikes = data['Strike'].tolist()

imp_vol = data['LR_Imp_Vol'].tolist()

f = interp2d(strikes,ttm,imp_vol, kind='cubic')

plot_strikes = np.linspace(data['Strike'].min(), data['Strike'].max(),25)

plot_ttm = np.linspace(0, data['TTM'].max(), 25)

fig = plt.figure()

ax = fig.gca(projection='3d')

X, Y = np.meshgrid(plot_strikes, plot_ttm)

Z = np.array([f(x,y) for xr, yr in zip(X, Y) for x, y in zip(xr,yr) ]).reshape(len(X), len(X[0]))

surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0.1)

ax.set_xlabel('Strikes')

ax.set_ylabel('Time-to-Maturity')

ax.set_zlabel('Implied Volatility')

fig.colorbar(surf, shrink=0.5, aspect=5)

Options and Volatility Smile

July 30, 2017August 6, 2017 suhasghorp Leave a comment

An equity option represents the right to buy (“call” option) or sell (“put” option) a unit of underlying stock at a pre-specified price (strike) at a predetermined maturity date (European option) or at any time up to the predetermined date (American option).

Option writer sells an option and option holder buys an option.

For a European call option on an index with strike 8,000 and index level of 8200 at maturity, the option holder receives the difference of $200 from option writer. This is called the instrinsic value or payoff of the option from the holder’s point of view.

The payoff function for a call option is $$ h_{T}(S,K) = max[S_{T}-K, 0] \tag{Eq. 1}$$

where T = maturity date, $\ S_T $ is the index level at maturity and K is the strike price.

In-the-money: a call (put) is in-the-money when S > K (S < K)
At-the-money: call or put is at-the-money when $\ S \approx K $
Out-of-the-money: a call is out-of-the-money when S < K (S > K)

A fair present value (is different than payoff) of a European call option is given by Black-Scholes formula:

$\ C_{0}^{*} = C^{BSM}(S_{0},K,T,r,\sigma) \tag{Eq. 2}$

$\S_{0} $ current index level (spot)
K strike price of the option
T time-to-maturity of the option
r risk-free short rate
$\ \sigma $ volatility or the std dev of the index returns

$\ C^{BSM} = S_{t} . N(d_{1}) – e^{r(T-t)} . K. N(d_{2})\tag{Eq. 3} $

where

$\displaystyle N(d) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{d} e^{-\frac{1}{2}x^{2}} dx$

$\displaystyle d1 = \frac{\log\frac{S_{t}}{K} + (r + \frac{\sigma^2}{2})(T-t)}{\sigma\sqrt{T-t}}$

$\displaystyle d2 = \frac{\log\frac{S_{t}}{K} + (r – \frac{\sigma^2}{2})(T-t)}{\sigma\sqrt{T-t}}$

# BSM option valuation
import math
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.optimize import brentq
from scipy.interpolate import interp2d
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# BSM option valuation

import math

import numpy as np

import matplotlib as mpl

import matplotlib.pyplot as plt

from scipy.integrate import quad

from scipy.optimize import brentq

from scipy.interpolate import interp2d

from mpl_toolkits.mplot3d import Axes3D

from matplotlib import cm

## Helper functions ##
def dN(x):
    ''' Probability density function of standard normal random variable x.'''
    return math.exp(-0.5 * x ** 2) / math.sqrt(2 * math.pi)

def N(d):
    ''' Cumulative density function of standard normal random variable x. '''
    return quad(lambda x: dN(x), -20, d, limit=50)[0]

def d1f(St, K, t, T, r, sigma):
    ''' Black-Scholes-Merton d1 function.
        Parameters see e.g. BSM_call_value function. '''
    d1 = (math.log(St / K) + (r + 0.5 * sigma ** 2)
        * (T - t)) / (sigma * math.sqrt(T - t))
    return d1

def BSM_call_value(St, K, t, T, r, sigma):
    ''' Calculates Black-Scholes-Merton European call option value
        Parameters
    ==========
    St: float
    stock/index level at time t
    K: float
    strike price
    t: float
    valuation date
    T: float
    date of maturity/time-to-maturity if t = 0; T > t
    r: float
    constant, risk-less short rate
    sigma: float
    volatility
    Returns
    =======
    call_value: float
    European call present value at t
    '''
    d1 = d1f(St, K, t, T, r, sigma)
    d2 = d1 - sigma * math.sqrt(T - t)
    call_value = St * N(d1) - math.exp(-r * (T - t)) * K * N(d2)
    return call_value

def BSM_put_value(St, K, t, T, r, sigma):
    ''' Calculates Black-Scholes-Merton European put option value.
    Parameters
    ==========
    St: float
    stock/index level at time t
    K: float
    strike price
    t: float
    valuation date
    T: float
    date of maturity/time-to-maturity if t = 0; T > t
    r: float
    constant, risk-less short rate
    sigma: float
    volatility
    Returns
    =======
    put_value: float
    European put present value at t
    '''
    put_value = BSM_call_value(St, K, t, T, r, sigma) - St + math.exp(-r * (T - t)) * K
    return put_value

## Helper functions ##

def dN(x):

''' Probability density function of standard normal random variable x.'''

return math.exp(-0.5 * x ** 2) / math.sqrt(2 * math.pi)

def N(d):

''' Cumulative density function of standard normal random variable x. '''

return quad(lambda x: dN(x), -20, d, limit=50)[0]

def d1f(St, K, t, T, r, sigma):

''' Black-Scholes-Merton d1 function.

Parameters see e.g. BSM_call_value function. '''

d1 = (math.log(St / K) + (r + 0.5 * sigma ** 2)

* (T - t)) / (sigma * math.sqrt(T - t))

return d1

def BSM_call_value(St, K, t, T, r, sigma):

''' Calculates Black-Scholes-Merton European call option value

Parameters

==========

St: float

stock/index level at time t

K: float

strike price

t: float

valuation date

T: float

date of maturity/time-to-maturity if t = 0; T > t

r: float

constant, risk-less short rate

sigma: float

volatility

Returns

=======

call_value: float

European call present value at t

'''

d1 = d1f(St, K, t, T, r, sigma)

d2 = d1 - sigma * math.sqrt(T - t)

call_value = St * N(d1) - math.exp(-r * (T - t)) * K * N(d2)

return call_value

def BSM_put_value(St, K, t, T, r, sigma):

''' Calculates Black-Scholes-Merton European put option value.

Parameters

==========

St: float

stock/index level at time t

K: float

strike price

t: float

valuation date

T: float

date of maturity/time-to-maturity if t = 0; T > t

r: float

constant, risk-less short rate

sigma: float

volatility

Returns

=======

put_value: float

European put present value at t

'''

put_value = BSM_call_value(St, K, t, T, r, sigma) - St + math.exp(-r * (T - t)) * K

return put_value

# Test Option Payoff and Time value
K = 8000 # Strike price
T = 1.0 # time-to-maturity
r = 0.025 # constant, risk-free rate
vol = 0.2 # constant volatility

# Generate spot prices 
S = np.linspace(4000, 12000, 150)
h = np.maximum(S - K, 0) # payoff of the option
C = [BSM_call_value(Szero, K, 0, T, r, vol) for Szero in S] #BS call option values

plt.figure()
plt.plot(S, h, 'b-.', lw=2.5, label='payoff') # plot inner value at maturity
plt.plot(S, C, 'r', lw=2.5, label='present value') # plot option present value
plt.grid(True)
plt.legend(loc=0)
plt.xlabel('index level $S_0$')
plt.ylabel('present value $C(t=0)$')

# Test Option Payoff and Time value

K = 8000 # Strike price

T = 1.0 # time-to-maturity

r = 0.025 # constant, risk-free rate

vol = 0.2 # constant volatility

# Generate spot prices

S = np.linspace(4000, 12000, 150)

h = np.maximum(S - K, 0) # payoff of the option

C = [BSM_call_value(Szero, K, 0, T, r, vol) for Szero in S] #BS call option values

plt.figure()

plt.plot(S, h, 'b-.', lw=2.5, label='payoff') # plot inner value at maturity

plt.plot(S, C, 'r', lw=2.5, label='present value') # plot option present value

plt.grid(True)

plt.legend(loc=0)

plt.xlabel('index level $S_0$')

plt.ylabel('present value $C(t=0)$')

The present value of the option is always higher than the undiscounted payoff, the difference being the time value. In other words, the option’s present value is composed of payoff plus the time value. Time value indicates that there is always a chance of option going in-the-money or more in-the-money during that time.

Simulating Returns

The Geometric Brownian motion model of the BS equation is given by

$$\displaystyle dS_{t} = rS_{t}dt + \sigma S_{t} dt dZ_{t}\tag{Eq.4}$$

The discretized version is

$$\displaystyle S_{t} = S_{t – \Delta t} e^{(r – \frac{1}{2}\sigma^2) \Delta t + \sigma \sqrt{\Delta t} z_{t}}\tag{Eq.5}$$

where t $\in {(\Delta t, 2\Delta t,…..,T)}$

Using the above discretized version, we will simulate the spot prices with $S_{0}$=100, T=10, r = 0.05 and $\sigma$=0.2

import pandas as pd

# model parameters
S0 = 100.0 # initial index level
T = 10.0 # time horizon
r = 0.05 # risk-less short rate
vol = 0.2 # instantaneous volatility
# simulation parameters
np.random.seed(250000)
#generate a pd array with business dates, ignores holidays
gbm_dates = pd.DatetimeIndex(start='10-05-2007',end='10-05-2017',freq='B')
M = len(gbm_dates) # time steps
I = 1 # index level paths
dt = 1 / 252. # 252 business days a year
df = math.exp(-r * dt) # discount factor

# stock price paths
rand = np.random.standard_normal((M, I)) # random numbers
S = np.zeros_like(rand) # stock matrix
S[0] = S0 # initial values
for t in range(1, M): # stock price paths using Eq.5
    S[t] = S[t - 1] * np.exp((r - vol ** 2 / 2) * dt + vol * rand[t] * math.sqrt(dt))
#create a pd dataframe with date as index and a column named "spot"
gbm = pd.DataFrame(S[:, 0], index=gbm_dates, columns=['spot']) 
gbm['returns'] = np.log(gbm['spot'] / gbm['spot'].shift(1)) #log returns
# Realized Volatility 
gbm['realized_var'] = 252 * np.cumsum(gbm['returns'] ** 2) / np.arange(len(gbm))
gbm['realized_vol'] = np.sqrt(gbm['realized_var'])
print gbm.head()
gbm = gbm.dropna()

import pandas as pd

# model parameters

S0 = 100.0 # initial index level

T = 10.0 # time horizon

r = 0.05 # risk-less short rate

vol = 0.2 # instantaneous volatility

# simulation parameters

np.random.seed(250000)

#generate a pd array with business dates, ignores holidays

gbm_dates = pd.DatetimeIndex(start='10-05-2007',end='10-05-2017',freq='B')

M = len(gbm_dates) # time steps

I = 1 # index level paths

dt = 1 / 252. # 252 business days a year

df = math.exp(-r * dt) # discount factor

# stock price paths

rand = np.random.standard_normal((M, I)) # random numbers

S = np.zeros_like(rand) # stock matrix

S[0] = S0 # initial values

for t in range(1, M): # stock price paths using Eq.5

S[t] = S[t - 1] * np.exp((r - vol ** 2 / 2) * dt + vol * rand[t] * math.sqrt(dt))

#create a pd dataframe with date as index and a column named "spot"

gbm = pd.DataFrame(S[:, 0], index=gbm_dates, columns=['spot'])

gbm['returns'] = np.log(gbm['spot'] / gbm['spot'].shift(1)) #log returns

# Realized Volatility

gbm['realized_var'] = 252 * np.cumsum(gbm['returns'] ** 2) / np.arange(len(gbm))

gbm['realized_vol'] = np.sqrt(gbm['realized_var'])

print gbm.head()

gbm = gbm.dropna()

               spot   returns realized_var realized_vol 
2007-10-08 98.904295 -0.011018 0.030589    0.174898
2007-10-09 98.444867 -0.004656 0.018026    0.134261
2007-10-10 97.696364 -0.007632 0.016911    0.130041
2007-10-11 98.280594 0.005962  0.014922    0.122158

spot returns realized_var realized_vol

2007-10-08 98.904295 -0.011018 0.030589 0.174898

2007-10-09 98.444867 -0.004656 0.018026 0.134261

2007-10-10 97.696364 -0.007632 0.016911 0.130041

2007-10-11 98.280594 0.005962 0.014922 0.122158

plt.figure(figsize=(9, 6))
plt.subplot(211)
gbm['spot'].plot()
plt.ylabel('daily quotes')
plt.grid(True)
plt.axis('tight')
plt.subplot(212)
gbm['returns'].plot()
plt.ylabel('daily log returns')
plt.grid(True)
plt.axis('tight')

plt.figure(figsize=(9, 6))

plt.subplot(211)

gbm['spot'].plot()

plt.ylabel('daily quotes')

plt.grid(True)

plt.axis('tight')

plt.subplot(212)

gbm['returns'].plot()

plt.ylabel('daily log returns')

plt.grid(True)

plt.axis('tight')

Implied Volatility is the value of $\sigma$ that solves Eq. 2 given the option market quote $C_{0}^{*}$

Volatility surface is the plot of the implied volatilities for different option strikes and different option maturities on the same underlying (an option chain).

Vol Surfaces exhibit :
Smiles: option implied volatilities exhibit a smile form, i.e. for calls the OTM implied volatilities are higher than the ATM ones; sometimes they rise again for ITM options
term structure:smiles are more pronounced for short-term options than for longer-term options; a phenomenon sometimes called volatility term structure

To demonstrate Vol Surface, I will use an option chain on AAPL stock as of 5/11/2017. I have downloaded this data from a reputable vendor, you can find this file here AAPL_BBG_vols

headers = ['Date','Strike','Call_Bid','Call_Ask','Call','Maturity','Put_Bid','Put_Ask','Put']
dtypes = {'Date': 'str', 'Strike': 'float', 'Call_Bid': 'float', 'Call_Ask': 'float',
          'Call':'float','Maturity':'str','Put_Bid':'float','Put_Ask':'float','Put':'float'}
parse_dates = ['Date', 'Maturity']

data = pd.read_csv('AAPL_BBG_vols.csv',skiprows=1,header=None, names=headers, dtype=dtypes, parse_dates=parse_dates)
data['BS_Imp_Vol'] = 0.0
data['TTM'] = 0.0
r = 0.0248 #risk-free rate
S0 = 153.97 # spot price as of 5/11/2017
div = 0.0182
for i in data.index:
    t = data['Date'][i]
    T = data['Maturity'][i]
    K = data['Strike'][i]
    Put = data['Put'][i]
    time_to_maturity = (T - t).days/365.0
    data.loc[i, 'TTM'] = time_to_maturity
    def func_BS(sigma):
        return BSM_put_value(S0, K, 0, time_to_maturity, r, sigma) - Put

    bs_imp_vol = brentq(func_BS, 0.03, 1.0)
    data.loc[i,'BS_Imp_Vol'] = bs_imp_vol

headers = ['Date','Strike','Call_Bid','Call_Ask','Call','Maturity','Put_Bid','Put_Ask','Put']

dtypes = {'Date': 'str', 'Strike': 'float', 'Call_Bid': 'float', 'Call_Ask': 'float',

'Call':'float','Maturity':'str','Put_Bid':'float','Put_Ask':'float','Put':'float'}

parse_dates = ['Date', 'Maturity']

data = pd.read_csv('AAPL_BBG_vols.csv',skiprows=1,header=None, names=headers, dtype=dtypes, parse_dates=parse_dates)

data['BS_Imp_Vol'] = 0.0

data['TTM'] = 0.0

r = 0.0248 #risk-free rate

S0 = 153.97 # spot price as of 5/11/2017

div = 0.0182

for i in data.index:

t = data['Date'][i]

T = data['Maturity'][i]

K = data['Strike'][i]

Put = data['Put'][i]

time_to_maturity = (T - t).days/365.0

data.loc[i, 'TTM'] = time_to_maturity

def func_BS(sigma):

return BSM_put_value(S0, K, 0, time_to_maturity, r, sigma) - Put

bs_imp_vol = brentq(func_BS, 0.03, 1.0)

data.loc[i,'BS_Imp_Vol'] = bs_imp_vol

Once we have the implied volatilities, we will generate a grid of strikes and maturities and use Cubic interpolation to derive the missing implied volatilities needed for a smooth surface.

ttm = data['TTM'].tolist()
strikes = data['Strike'].tolist()
imp_vol = data['BS_Imp_Vol'].tolist()
f = interp2d(strikes,ttm,imp_vol, kind='cubic')

plot_strikes = np.linspace(data['Strike'].min(), data['Strike'].max(),25)
plot_ttm = np.linspace(0, data['TTM'].max(), 25)
fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(plot_strikes, plot_ttm)
Z = np.array([f(x,y) for xr, yr in zip(X, Y) for x, y in zip(xr,yr) ]).reshape(len(X), len(X[0]))
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0.1)
ax.set_xlabel('Strikes')
ax.set_ylabel('Time-to-Maturity')
ax.set_zlabel('Implied Volatility')
fig.colorbar(surf, shrink=0.5, aspect=5)

ttm = data['TTM'].tolist()

strikes = data['Strike'].tolist()

imp_vol = data['BS_Imp_Vol'].tolist()

f = interp2d(strikes,ttm,imp_vol, kind='cubic')

plot_strikes = np.linspace(data['Strike'].min(), data['Strike'].max(),25)

plot_ttm = np.linspace(0, data['TTM'].max(), 25)

fig = plt.figure()

ax = fig.gca(projection='3d')

X, Y = np.meshgrid(plot_strikes, plot_ttm)

Z = np.array([f(x,y) for xr, yr in zip(X, Y) for x, y in zip(xr,yr) ]).reshape(len(X), len(X[0]))

surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0.1)

ax.set_xlabel('Strikes')

ax.set_ylabel('Time-to-Maturity')

ax.set_zlabel('Implied Volatility')

fig.colorbar(surf, shrink=0.5, aspect=5)