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

References

Sinclair, E. (2010). Option trading: pricing and volatility strategies and techniques. Hoboken, N.J: Wiley.

Natenberg, S. (2015). Option volatility and pricing: advanced trading strategies and techniques. New York: McGraw-Hill Education.

Related Posts

Categories

Recent Posts