Calculating Implied Volatility... Fast

How to calculate IV in a computationally effective manner.

May 08, 2024

∙ Paid

Introduction

I wasn’t entirely sure this topic would be interesting to most readers—it sure wasn’t to me at first. Despite this, I still managed to encounter it in my work, and my bet is that you’ll probably find it useful, too. This is if you ever work with a large (they’re almost all massive so this is perhaps redundant) dataset of option prices and require IV values, which I think isn’t too unlikely for most researchers. The libraries available are great, but knowing how it works under the hood to fit your needs will save you a lot of sitting about waiting for calculations to finish.

It is often overlooked that there is a lot of computational work that goes into efficiently calculating implied volatilities from option prices. If you are working with orderbook data, and performing research on the data - you will 100% encounter the issue of calculating Greeks and implied volatility being a compute bottleneck. Any large-scale data analysis using Greeks will require this, and if you need to accurately figure out how your risk has shifted after every market event, even the beefiest of trading systems will become clogged without fast calculations.

Fast Greek calculations are not simply about being cutting edge on the latency front - they are often about being able to do certain data analysis tasks or have specific risk management update frequencies AT ALL.

In this article, we will look at how to calculate this efficiently using Rust and Python. We will compare multiple methods and benchmark their performance across various option prices.

Index

Introduction
Index
Newton
Really Fast Newton???
SR
Low Accuracy
Let’s Be Rational
How does Optiver etc. do it

Newton

The Newton-Raphson method is a root-finding algorithm that uses iterations to find successively better approximations to the roots (or zeroes) of a real-valued function. The formula for the Newton-Raphson method is:

\(x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \)

Where:

x_n is the current guess,
f(x_n) is the function value at 𝑥_𝑛,
𝑓′(x_n) is the derivative of the function at 𝑥_𝑛,
𝑥_𝑛+1 is the next guess.

When applied to finding the implied volatility, the function f becomes the difference between the market price of the option and the theoretical price of the option calculated using the Black-Scholes formula:

\(f(\sigma) = C_{BS}(\sigma) - C_{market} \)

Here 𝜎 is the implied volatility, 𝐶_𝐵𝑆 is the Black-Scholes price of the option, and 𝐶_𝑚𝑎𝑟𝑘𝑒𝑡 is the actual market price of the option. The derivative 𝑓′(𝜎) is the Vega of the option or the rate of change of the option price with respect to changes in volatility.

Visually:

Newton-Raphson Method— Explained and Visualised | Towards Data Science

Implementing a demonstration of this in Python:

from scipy.optimize import newton
import numpy as np
import scipy.stats as si

def bs_call(S, K, r, sigma, T):
    """Calculate Black-Scholes call option price."""
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S * si.norm.cdf(d1) - K * np.exp(-r * T) * si.norm.cdf(d2)

def bs_put(S, K, r, sigma, T):
    """Calculate Black-Scholes put option price."""
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return K * np.exp(-r * T) * si.norm.cdf(-d2) - S * si.norm.cdf(-d1)

def bs_vega(S, K, r, sigma, T):
    """Calculate Vega of an option, which is the derivative of the price with respect to sigma."""
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    return S * si.norm.pdf(d1) * np.sqrt(T)

def nr_solve(option_type, price, S, K, r, T, initial_guess = 0.2):
    """Function to estimate the implied volatility for both call and put options using Newton-Raphson method."""
    if option_type == 'call':
        price_fn = bs_call
    elif option_type == 'put':
        price_fn = bs_put
    else:
        raise ValueError("Invalid option type. Must be 'call' or 'put'.")

    obj_fn = lambda sigma: price_fn(S, K, r, sigma, T) - price
    dfn = lambda sigma: bs_vega(S, K, r, sigma, T)

    try:
        result = newton(func=obj_fn, x0=initial_guess, fprime=dfn, tol=1e-5, maxiter=100)
        return result
    except RuntimeError as e:
        print("Failed to converge:", e)
        return None

We get these performance results:

It isn’t exactly the fastest, that’s for sure, but it’s very accurate. This is our first method, so don’t be too disappointed; we’ll get Newton up to standard in no time. We also wrote it in Python… let’s try some Rust instead.

Really Fast Newton

How can we make Newton faster? Well, we already got a speed boost from using Scipy, but we can further improve it by using Rust and nrfind::find_root. Something that might look a bit like this:

The Quant Stack