Implied Volatility Calculations with Python

Implied volatility $\sigma_{imp}$ is the volatility value $\sigma$ that makes the Black-Scholes value of the option equal to the traded price of the option.

Recall that in the Black-Scholes model, the volatility parameter $\sigma$ is the only parameter that can't be directly observed. All other parameters can be determined through market data (in the case of the risk-free rate $r$ and dividend yield $q$ and when the option is quoted. This being the case, the volatility parameter is the result of a numerical optimization technique given the Black-Scholes model.

The Greek 'Vega'¶

Vega is the first derivative of $\sigma$ volatility and thus is an integral piece in the formulation of implied volatility. What follows is a quick derivation of Vega.

As Vega is the first derivative of volatility, its partial derivative takes the form $\frac{\partial C}{\partial \sigma}$. Therefore, we take the partial derivative of the Black-Scholes formula with respect to $\sigma$

$$vega(C) = \frac{\partial C}{\partial \sigma} = Se^{-q(T-t)} N'(d_1) \frac{\partial d_1}{\partial \sigma} - Ke^{-r(T-t)} N'(d_2) \frac{\partial d_2}{\partial \sigma}$$

We can then rearrange the formula into:

$$vega(C) = Se^{-q(T-t)} N'(d_1) \left(\frac{\partial d_1}{\partial \sigma} - \frac{\partial d_2}{\partial \sigma} \right)$$

We know $d_2 = d_1 - \sigma \sqrt{T - t}$, rearrange to get $d_1 - d_2 = \sigma \sqrt{T-t}$. Therefore,

$$\frac{\partial d_1}{\partial \sigma} - \frac{\partial d_2}{\partial \sigma} = \sqrt{T - t}$$

Can then replace $N'(d_1)$ as $N'(d_1) = \frac{1}{\sqrt{2\pi}} e^{\frac{-d_1^2}{2}}$ to get:

$$vega(C) = \frac{1}{\sqrt{2\pi}} Se^{-q(T-t)} e^{\frac{-d_1^2}{2}} \sqrt{T-t}$$

For a non-dividend paying asset, the vega derives to the following:

$$S \cdot N(d_1) \sqrt{T - t}$$

Python Implementation of Vega for Dividend Paying Assets¶

In [93]:
import numpy as np
import scipy.stats as si

In [23]:
def vega_div(S, K, T, r, q, sigma):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#sigma: volatility of underlying asset
#q: continuous dividend rate

d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

vega = 1 / np.sqrt(2 * np.pi) * S * np.exp(-q * T) * np.exp(-(d1 ** 2) * 0.5) * np.sqrt(T)

return vega


Python Implementation of Vega for Non-Dividend Paying Assets¶

In [94]:
def vega(S, K, T, r, sigma):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#sigma: volatility of underlying asset

d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

vega = S * si.norm.cdf(d1, 0.0, 1.0) * np.sqrt(T)

return vega

In [13]:
vega(50, 100, 1, 0.05, 0.25)

Out[13]:
0.35953928148824355

Formulation¶

Dividend Paying Assets¶

Consider $C$ is the quoted value of a call. We then define the value of a Black-Scholes call option as:

$$C_{value}(S, K, T, \sigma_{imp}, r, q) = C$$

With the usual definitions of the parameters.

• $S$ = spot price of option
• $K$ = strike price
• $T$ = time to maturity
• $q$ = continuous dividend rate
• $r$ = interest rate
• $\sigma_{imp}$ = implied volatility

This can also be used for the price of a put option:

$$P_{value}(S, K, T, \sigma_{imp}, r, q) = P$$

As these are functions of volatility, values of both call and put options are increasing as:

$$\frac{\partial C_{value}}{\partial \sigma} = \frac{\partial P_{value}}{\partial \sigma} = vega(C_{value}) = vega(P_{value})$$

$$= Se^{-qT} \sqrt{T} \frac{1}{\sqrt{2\pi}} e^{\frac{-d_1^2}{2}} > 0$$

For implied volatility to exist and be strictly positive, the given value of $C$ of the call option must be arbitrage-free:

$$Se^{-qT} - Ke^{-rT} \leq C < Se^{-qT}$$

Also, volatility only exists and is positive if the value of the put option meets the following inequality:

$$Ke^{-rT} - Se^{-qT} \leq P < Ke^{-rT}$$

As mentioned previously, implied volatility is the only parameter in the Black-Scholes model that isn't directly observable. For a call or put option, the maturity and strike of the option are given, and when the option is traded, the price and spot price of the underlying is known as well.

If a lognormal model for the underlying is assumed, the price of the option given by Black-Scholes for $t = 0$ is:

$$C = Se^{-qT}N(d_1) - Ke^{-rT}N(d_2)$$

$$P = Ke^{-rT}N(-d_2) - Se^{-qT}N(-d_1)$$

We can look at the call and put functions as a function of the volatility parameter $\sigma$. Finding implied volatility requires solving the nonlinear problem $f(x) = 0$ where $x = \sigma$ and:

$$f(x) = Se^{-qT}N(d_1(x)) - Ke^{-rT}N(d_2(x)) - C$$

$$f(x) = Ke^{-rT}N(-d_2(x)) - Se^{-qT}N(-d_1(x)) - P$$

For a call and put option, respectively.

Implied volatility is often higher when deep out of or in the money than at the money options. This is called the volatility smile.

To solve the function when $f(x) = 0$, Newton's method is employed. Interestingly, differentating the call function above is the same as computing the vega of the option.

$$f'(x) = \frac{1}{\sqrt{2\pi}} Se^{-qT} \sqrt{T} exp \left(-\frac{(d_1(x))^2}{2} \right)$$

$f'(x) \geq 0$ for any $x$, the function $f(x)$ is increasing and therefore has at most one solution.

Recursion for Newton's method for solving is stated as:

$$x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}$$

Since we are dealing with implied volatilities, we can write it using $\sigma$

$$\sigma^{imp}_{k + 1} = \sigma^{imp} - \frac{f(\sigma^{imp}_k)}{f'(\sigma^{imp}_k)}, \forall k \geq 0$$

With $f(x) = Se^{-qT}N(d_1(x)) - Ke^{-rT}N(d_2(x)) - C$ and $f'(x) = \frac{1}{\sqrt{2\pi}}Se^{-qT}\sqrt{T} exp \left(-\frac{(d_1(x))^2}{2} \right)$

A good first estimate of volatility is $\sigma = 0.25$

Newton's Method for Nonlinear Equations¶

Non-dividend paying asset¶

In [73]:
def newton_vol_call(S, K, T, C, r, sigma):

#S: spot price
#K: strike price
#T: time to maturity
#C: Call value
#r: interest rate
#sigma: volatility of underlying asset

d1 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

fx = S * si.norm.cdf(d1, 0.0, 1.0) - K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0) - C

vega = (1 / np.sqrt(2 * np.pi)) * S * np.sqrt(T) * np.exp(-(si.norm.cdf(d1, 0.0, 1.0) ** 2) * 0.5)

tolerance = 0.000001
x0 = sigma
xnew  = x0
xold = x0 - 1

while abs(xnew - xold) > tolerance:

xold = xnew
xnew = (xnew - fx - C) / vega

return abs(xnew)

In [74]:
newton_vol_call(25, 20, 1, 7, 0.05, 0.25)

Out[74]:
0.33647905603288009
In [75]:
def newton_vol_put(S, K, T, P, r, sigma):

d1 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

fx = K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0) - S * si.norm.cdf(-d1, 0.0, 1.0) - P

vega = (1 / np.sqrt(2 * np.pi)) * S * np.sqrt(T) * np.exp(-(si.norm.cdf(d1, 0.0, 1.0) ** 2) * 0.5)

tolerance = 0.000001
x0 = sigma
xnew  = x0
xold = x0 - 1

while abs(xnew - xold) > tolerance:

xold = xnew
xnew = (xnew - fx - P) / vega

return abs(xnew)

In [76]:
newton_vol_put(25, 20, 1, 7, 0.05, 0.25)

Out[76]:
0.35295018049449461

Dividend paying asset¶

In [89]:
def newton_vol_call_div(S, K, T, C, r, q, sigma):

d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

fx = S * np.exp(-q * T) * si.norm.cdf(d1, 0.0, 1.0) - K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0) - C

vega = (1 / np.sqrt(2 * np.pi)) * S * np.exp(-q * T) * np.sqrt(T) * np.exp((-si.norm.cdf(d1, 0.0, 1.0) ** 2) * 0.5)

tolerance = 0.000001
x0 = sigma
xnew  = x0
xold = x0 - 1

while abs(xnew - xold) > tolerance:

xold = xnew
xnew = (xnew - fx - C) / vega

return abs(xnew)

In [90]:
newton_vol_call_div(25, 20, 1, 7, 0.05, 0.10, 0.25)

Out[90]:
0.40788841179831004
In [91]:
def newton_vol_put_div(S, K, T, P, r, q, sigma):

d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

fx = K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0) - S * np.exp(-q * T) * si.norm.cdf(-d1, 0.0, 1.0) -  P

vega = (1 / np.sqrt(2 * np.pi)) * S * np.exp(-q * T) * np.sqrt(T) * np.exp((-si.norm.cdf(d1, 0.0, 1.0) ** 2) * 0.5)

tolerance = 0.000001
x0 = sigma
xnew  = x0
xold = x0 - 1

while abs(xnew - xold) > tolerance:

xold = xnew
xnew = (xnew - fx - P) / vega

return abs(xnew)

In [92]:
newton_vol_put_div(25, 20, 1, 7, 0.05, 0.10, 0.25)

Out[92]:
0.037184024697147447