Implied Volatility

The implied volatility (IV) of an option contract is that value of the volatility of the underlying instrument which, when input in an option pricing model (such as Black–Scholes), will return a theoretical value equal to the current market price of that option. The VIX , in contrast, is a model-free estimate of Implied Volatility. The latter is viewed as being important because it represents a measure of risk for the underlying asset. Elevated Implied Volatility suggests that risks to underlying are also elevated. Ordinarily, to estimate implied volatility we rely upon Black-Scholes (1973). This implies that we are prepared to accept the assumptions of Black Scholes (1973). To estimate Implied Volatility we need, first and foremost, to set up a fully functional Option Pricing model. The simplest way to do this is to work with excel. A good way to way to acquaint one self with the steps involved in implementing Black Scholes is set up the estimation in excel .

Implied Volatility estimated with Goal Seek

Implied Volatility is important to data scientists because because it provides an ex ante measure of expected risk - which can be different diverge from historical time series projections gleaned from GARCH style models. Again, I use excel here to set out a simplified treatment invoking goal-seek. Please follow link.

Excel VBA code for Black Scholes

To estimate implied volatility for a number of options goal seek can become a bit clunky. To automate that estimation it is worth setting out VBA code for the Black Schole Call and couple this with a bisection function to run to invoke a trial and error iterative search for the Implied Volatility consistent with the market price of the option.

Function BSC(S, K, r, q, sigma, T)

Dim dOne, dTwo, Nd1, Nd2

dOne = (Log(S / K) + (r - q + 0.5 * sigma ^ 2) * T) / (sigma * Sqr(T))

dTwo = dOne - sigma * Sqr(T)

Nd1 = Application.NormSDist(dOne)

Nd2 = Application.NormSDist(dTwo)

BSC = Exp(-q * T) * S * Nd1 - Exp(-r * T) * K * Nd2

End Function

to then estimate implied volatility we need the Bisection Function

Function BSCImVol(S, K, r, q, T, callmktprice)

H = 5

L = 0

Do While (H - L) > 0.00000001

If BSC(S, K, r, q, (H + L) / 2, T) > callmktprice Then

H = (H + L) / 2

Else: L = (H + L) / 2

End If

Loop

BSCImVol = (H + L) / 2

End Function

Please see video below.

C++ code for Implied Volatility

To implement Implied Volatility in C++ please use the following code:

//Implied Volatility using Black Scholes with Dividend

#include <iostream>

#include <cmath>

#include <vector>

using namespace std;

// N(0,1) density

double

f (double x)

{

double pi = 4.0 * atan (1.0);

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

}

// Boole's Rule

double

Boole (double StartPoint, double EndPoint, int n)

{

vector < double >X (n + 1, 0.0);

vector < double >Y (n + 1, 0.0);

double delta_x = (EndPoint - StartPoint) / double (n);

for (int i = 0; i <= n; i++)

{

X[i] = StartPoint + i * delta_x;

Y[i] = f (X[i]);

}

double sum = 0;

for (int t = 0; t <= (n - 1) / 4; t++)

{

int ind = 4 * t;

sum +=

(1 / 45.0) * (14 * Y[ind] + 64 * Y[ind + 1] + 24 * Y[ind + 2] +

64 * Y[ind + 3] + 14 * Y[ind + 4]) * delta_x;

}

return sum;

}

// N(0,1) cdf by Boole's Rule

double

N (double x)

{

return Boole (-10.0, x, 240);

}

// Black-Scholes Call Price

double

BSPrice (double S, double K, double r, double T, double q, double v, char PutCall)

{

double d = (log (S / K) + T * (r - q + 0.5 * v * v)) / (v * sqrt (T));

double call = S *exp(-q*T)* N (d) - exp (-r * T) * K * N (d - v * sqrt (T));

if (PutCall == 'C')

return call;

else

// Put Parity

return call - S*exp (-q * T) + K * exp (-r * T);

}

// Bisection Algorithm

double BisecBSV(double S, double K, double r, double T, double q,

double a, double b, double MktPrice, char PutCall) {

const int MaxIter = 50000;

double Tol = 0.0000001;

double midP = 0.0, midCdif;

double lowCdif = MktPrice - BSPrice(S, K, r, T, q, a, PutCall);

double highCdif = MktPrice - BSPrice(S, K, r, T, q, b, PutCall);

if (lowCdif*highCdif > 0)

return -1;

else

for (int i = 0; i <= MaxIter; i++) {

midP = (a + b) / 2.0;

midCdif = MktPrice - BSPrice(S, K, r, T, q, midP, PutCall);

if (abs(midCdif)<Tol) goto LastLine;

else {

if (midCdif>0) a = midP;

else b = midP;

}

LastLine:

return midP;

}

// Main program

int main() {

double S = 100; // Spot Price

double K = 100; // Exercise

double r = 0.05; // Interest Rate

double T = 1; // Maturity in Years

double q = 0.0;

double a = 0.00000001; // Bisection algorithm starting value for lower volatility

double b = 7.0; // Bisection algorithm starting value for higher volatility

double MktPrice = 10.45;

char PutCall = 'C'; // C

double IV = BisecBSV(S, K, r, T, q, a, b, MktPrice, PutCall);

cout << " Implied Volatility = " << IV << " " << endl;

system("PAUSE");

}

Google Colaboratory

Javascript Code for Black Scholes and Implied Volatility

And for good measure I have included the following Javacsript code for Black Scholes and Implied Volatility. Please check below video for estimation.

/* Returns probability of occuring below and above target price. */

function probability(price, target, days, volatility) {

var p = price;

var q = target;

var t = days / 365;

var v = volatility;

var vt = v*Math.sqrt(t);

var lnpq = Math.log(q/p);

var d1 = lnpq / vt;

var y = Math.floor(1/(1+.2316419*Math.abs(d1))*100000)/100000;

var z = Math.floor(.3989423*Math.exp(-((d1*d1)/2))*100000)/100000;

var y5 = 1.330274*Math.pow(y,5);

var y4 = 1.821256*Math.pow(y,4);

var y3 = 1.781478*Math.pow(y,3);

var y2 = 0.356538*Math.pow(y,2);

var y1 = 0.3193815*y;

var x = 1-z*(y5-y4+y3-y2+y1);

x = Math.floor(x*100000)/100000;

if (d1<0) {x=1-x};

var pabove = Math.floor(x*1000)/10;

var pbelow = Math.floor((1-x)*1000)/10;

return [[pbelow],[pabove]]

}

// JavaScript adopted from Bernt Arne Odegaard's Financial Numerical Recipes

// http://finance.bi.no/~bernt/gcc_prog/algoritms/algoritms/algoritms.html

// by Steve Derezinski, CXWeb, Inc. http://www.cxweb.com

// This program is free software; you can redistribute it and/or

// modify it under the terms of the GNU General Public License

// as published by the Free Software Foundation.

// This program is distributed in the hope that it will be useful,

// but WITHOUT ANY WARRANTY; without even the implied warranty of

// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

// GNU General Public License for more details.

// http://www.fsf.org/copyleft/gpl.html

function ndist(z) {

return (1.0/(Math.sqrt(2*Math.PI)))*Math.exp(-0.5*z);

//?? Math.exp(-0.5*z*z)

}

function N(z) {

b1 = 0.31938153;

b2 = -0.356563782;

b3 = 1.781477937;

b4 = -1.821255978;

b5 = 1.330274429;

p = 0.2316419;

c2 = 0.3989423;

a=Math.abs(z);

if (a>6.0) {return 1.0;}

t = 1.0/(1.0+a*p);

b = c2*Math.exp((-z)*(z/2.0));

n = ((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t;

n = 1.0-b*n;

if (z < 0.0) {n = 1.0 - n;}

return n;

}

function fraction(z) {

// given a decimal number z, return a string with whole number + fractional string

// i.e. z = 4.375, return "4 3/8"

var whole = Math.floor(z);

var fract = z - whole;

var thirtytwos = Math.round(fract*32);

if (thirtytwos == 0) {return whole + " ";} //(if fraction is < 1/64)

if (thirtytwos == 32) {return whole + 1;} //(if fraction is > 63/64)

//32's non-trivial denominators: 2,4,8,16

if (thirtytwos/16 == 1) { return whole + " 1/2";}

if (thirtytwos/8 == 1) { return whole + " 1/4";}

if (thirtytwos/8 == 3) { return whole + " 3/4";}

if (thirtytwos/4 == Math.floor(thirtytwos/4)) {return whole + " " + thirtytwos/4 + "/8";}

if (thirtytwos/2 == Math.floor(thirtytwos/2)) {return whole + " " + thirtytwos/2 + "/16";}

else return whole + " " + thirtytwos + "/32";

} //end function

function black_scholes(call,S,X,r,v,t) {

// call = Boolean (to calc call, call=True, put: call=false)

// S = stock prics, X = strike price, r = no-risk interest rate

// v = volitility (1 std dev of S for (1 yr? 1 month?, you pick)

// t = time to maturity

// define some temp vars, to minimize function calls

var sqt = Math.sqrt(t);

var Nd2; //N(d2), used often

var nd1; //n(d1), also used often

var ert; //e(-rt), ditto

var delta; //The delta of the option

d1 = (Math.log(S/X) + r*t)/(v*sqt) + 0.5*(v*sqt);

d2 = d1 - (v*sqt);

if (call) {

delta = N(d1);

Nd2 = N(d2);

} else { //put

delta = -N(-d1);

Nd2 = -N(-d2);

}

ert = Math.exp(-r*t);

nd1 = ndist(d1);

gamma = nd1/(S*v*sqt);

vega = S*sqt*nd1;

theta = -(S*v*nd1)/(2*sqt) - r*X*ert*Nd2;

rho = X*t*ert*Nd2;

return ( S*delta-X*ert *Nd2);

} //end of black_scholes

function option_implied_volatility(call,S,X,r,t,o) {

// call = Boolean (to calc call, call=True, put: call=false)

// S = stock prics, X = strike price, r = no-risk interest rate

// t = time to maturity

// o = option price

// define some temp vars, to minimize function calls

sqt = Math.sqrt(t);

MAX_ITER = 100;

ACC = 0.0001;

sigma = (o/S)/(0.398*sqt);

for (i=0;i<MAX_ITER;i++) {

price = black_scholes(call,S,X,r,sigma,t);

diff = o-price;

if (Math.abs(diff) < ACC) return sigma;

d1 = (Math.log(S/X) + r*t)/(sigma*sqt) + 0.5*sigma*sqt;

vega = S*sqt*ndist(d1);

sigma = sigma+diff/vega;

}

return "Error, failed to converge";

} //end of option_implied_volatility

//function call_iv(s,x,r,t,o) { return option_implied_volatility(true,s,x,r/100,t/365,o) }

Python Code for Black Scholes Implied Volatility Calls and Puts

The code below I obtained from Kevin Mooney's github. I made a few changes for dividends. Kevin does not use a bisection technique here. Rather he uses a root finding approach. Please check out Kevin's youtube channel. When selecting parameters please ensure that values used satisfy option properties. Also some parameter values may not produce convergence - no matter how many iterations are applied.

# This Python code can be used to calculate the value for Implied

# Volatility for a European call

from math import sqrt, exp, log, pi

from scipy.stats import norm

# Function 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

# Option parameters

S = 100.0

K = 100

t = 1

r = 0.05

q = 0.05

C0 =10.45

# Tolerances

tol = 1e-3

epsilon = 1

# Variables to log and manage number of iterations

count = 0

max_iter = 1000

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

vol = 0.50

while epsilon > tol:

# Count how many iterations and make sure while loop doesn't run away

count += 1

if count >= max_iter:

print('Breaking on count')

break;

# Log the value previously calculated to computer percent change

# between iterations

orig_vol = vol

# Calculate the vale of the call price

d1, d2 = d(vol, S, K, r,q, t)

function_value = call_price(vol, S, K, r, q, t, d1, d2) - C0

# Calculate vega, the derivative of the price with respect to

# volatility

vega = S * norm.pdf(d1) * sqrt(t)* exp(-q * t)

# Update for value of the volatility

vol = -function_value / vega + vol

# Check the percent change between current and last iteration

epsilon = abs( (vol - orig_vol) / orig_vol )

# Print out the results

print('Sigma = ', vol)

print('Code took ', count, ' iterations')

The Python code below follows Kevin's code above. I just make a few changes that invoke put call parity.

# This Python code can be used to calculate the value for Implied

# Volatility for a European put

from math import sqrt, exp, log, pi

from scipy.stats import norm

# Function to calculate the values of d1 and d2 as well as the call

# price. To extend to puts, define the put from put call parity

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

# From Put call Prity

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

P = - S * exp(-q * t) + K * exp(-r * t) + call_price(sigma, S, K, r, q, t, d1, d2)

return P

# Option parameters

S = 100.0

K = 100

t = 1

r = 0.05

q = 0.0

P0 =5.57

# Tolerances

tol = 1e-3

epsilon = 1

# Variables to log and manage number of iterations

count = 0

max_iter = 1000

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

vol = 0.50

while epsilon > tol:

# Count how many iterations and make sure while loop doesn't run away

count += 1

if count >= max_iter:

print('Breaking on count')

break;

# Log the value previously calculated to computer percent change

# between iterations

orig_vol = vol

# Calculate the vale of the call price

d1, d2 = d(vol, S, K, r,q, t)

function_value = put_price(vol, S, K, r, q, t, d1, d2) - P0

# Calculate vega, the derivative of the price with respect to

# volatility

vega = S * norm.pdf(d1) * sqrt(t)* exp(-q * t)

# Update for value of the volatility

vol = -function_value / vega + vol

# Check the percent change between current and last iteration

epsilon = abs( (vol - orig_vol) / orig_vol )

# Print out the results

print('Sigma = ', vol)

print('Code took ', count, ' iterations')

The video below outlines the how to run code snippets above in Anaconda Spyder. We then compare results in Python OnlineGBD and then test results against C++ code that we know is correct. In that way we test the robustness of the Python Code.

The video below uses Amer Lu's Python code to estimate Black Scholes/Implied Volatility/Greeks in a Jupyter Notebook

Python code for Black Scholes Implied Volatility using Bisection

The code below can be used to estimate implied volatility using Bisection.

# https://github.com/YuChenAmberLu?tab=repositories

from math import log, sqrt, pi, exp

from scipy.stats import norm

## define two functions, d1 and d2 in Black-Scholes model

# change sigma*sqrt(T) to (sigma*sqrt(T))

def d1(S,K,T,r,q,sigma):

return(log(S/K)+(r - q +sigma**2/2.)*T)/(sigma*sqrt(T))

def d2(S,K,T,r,q,sigma):

return d1(S,K,T,r,q,sigma)-sigma*sqrt(T)

## define the call options price function

def bs_call(S,K,T,r,q,sigma):

return S*exp(-q*T)*norm.cdf(d1(S,K,T,r,q,sigma))-K*exp(-r*T)*norm.cdf(d2(S,K,T,r,q,sigma))

## define the put options price function

def bs_put(S,K,T,r,q,sigma):

return K*exp(-r*T)-S*exp(-q*T)+bs_call(S,K,T,r,q,sigma)

# Implied Volatility using bisection

# quantconnect.com/forum/discussion/2269/generate-volatility-surface-plot-by-interpolation/p1

def implied_vol(option_type, option_price, S, K, r, T, q):

# apply bisection method to get the implied volatility by solving the BSM function

precision = 0.00001

upper_vol = 50.0

max_vol = 50.0

min_vol = 0.0001

lower_vol = 0.0001

iteration = 0

while 1:

iteration +=1

mid_vol = (upper_vol + lower_vol)/2.0

if option_type == 'c':

price = bs_call(S,K,T,r,q,mid_vol)

lower_price = bs_call(S,K,T,r,q,lower_vol)

if (lower_price - option_price) * (price - option_price) > 0:

lower_vol = mid_vol

else:

upper_vol = mid_vol

if abs(price - option_price) < precision: break

if mid_vol > max_vol - 5 :

mid_vol = 0.0001

break

elif option_type == 'p':

price = bs_put(S,K,T,r,q,mid_vol)

upper_price = bs_put(S,K,T,r,q,upper_vol)

if (upper_price - option_price) * (price - option_price) > 0:

upper_vol = mid_vol

else:

lower_vol = mid_vol

if abs(price - option_price) < precision: break

if iteration > 100: break

return mid_vol

print (implied_vol('p', 5.57, 100, 100, 0.05, 1, 0))

Google Colaboratory