Forward Measures Resolve Bond Pricing Expectation Differences

# Deriving implied volatility programmatically


# Deriving implied volatility programmatically












I'm working on a project to calculate the value of options using `Python`. I'm using the Black-Scholes model, and I can get accurate results by plugging in a given value for implied volatility. I usually get IV either from the textbook I have, or from optionsprofitcalculator.

How can I derive IV on my own? I found two approaches so far:

- This answer which says IV can be derived by iteration using the BSM. I've read about this in other answers as well, but I don't think I'm understanding it correctly. It seems like it says to plug in guesses and adjust, until you find the IV that gives the correct value to the option. But to recognize the correct IV, I would need to already have the answer, like from the textbook or OPC, right? How would you implement this programmatically, without checking against another source?

- This question which says you can derive IV from model-free variance swaps. Would it be a better fit for my project, compared to the iteration method? If so, how would I implement it in practice, programmatically?

## Answer by amdopt (score 7, accepted)

https://quant.stackexchange.com/a/49724

Approximating implied volatility of European options can be done in a few ways--this is just one. Below is a python implementation that uses Newton Raphson. You can use the `implied_volatility` function to find the approximate implied volatility. You can then check it by plugging the output from that back into the `option_price` function.

```
import numpy as np
from scipy.stats import norm

"""
right = 'C' or 'P'
s = Spot Price
k =  Strike Price
t = Days to expiration
rfr = Risk-free Rate
sigma = volatility
div = Annual dividend rate. Defaulted to zero.
price = Known option price. Needed for implied_volatility function
"""

def d_one(s, k, t, rfr, sigma, div=0):
    """d1 calculation"""
    d_1 = (np.log(s / k) +
           (rfr - div + sigma ** 2 / 2) * t) / (sigma * np.sqrt(t))
    return d_1

def d_two(s, k, t, rfr, sigma, div=0):
    """d2 calculation"""
    d_2 = d_one(s, k, t, rfr, sigma, div) - sigma * np.sqrt(t)
    return d_2

def nd_one(right, s, k, t, rfr, sigma, div=0):
    """nd1 calculation"""
    if right == 'C':
        nd_1 = norm.cdf(d_one(s, k, t, rfr, sigma, div), 0, 1)
    elif right == 'P':
        nd_1 = norm.cdf(-d_one(s, k, t, rfr, sigma, div), 0, 1)
    return nd_1

def nd_two(right, s, k, t, rfr, sigma, div=0):
    """nd2 calculation"""
    if right == 'C':
        nd_2 = norm.cdf(d_two(s, k, t, rfr, sigma, div), 0, 1)
    elif right == 'P':
        nd_2 = norm.cdf(-d_two(s, k, t, rfr, sigma, div), 0, 1)
    return nd_2

def option_price(right, s, k, t, rfr, sigma, div=0):
    """option price"""
    right = right.upper()
    t /= 365
    if right == 'C':
        price = (s * np.exp(-div * t) *
                 nd_one(right, s, k, t, rfr, sigma, div)
                 - k * np.exp(-rfr * t) *
                 nd_two(right, s, k, t, rfr, sigma, div))
    elif right == 'P':
        price = (k * np.exp(-rfr * t) *
                 nd_two(right, s, k, t, rfr, sigma, div)
                 - s * np.exp(-div * t) *
                 nd_one(right, s, k, t, rfr, sigma, div))
    return price

def option_vega(s, k, t, rfr, sigma, div=0):
    """option vega"""
    t /= 365
    vega = (.01 * s * np.exp(-div * t) * np.sqrt(t)
            * norm.pdf(d_one(s, k, t, rfr, sigma, div)))
    return vega

def implied_volatility(right, s, k, t, rfr, price, div=0):
    """implied volatility approximation"""
    epsilon = 0.00000001
    sigma = 1.0

    def newton_raphson(right, s, k, t, rfr, sigma, price, epsilon, div=0):
        diff = np.abs(option_price(right, s, k, t, rfr, sigma, div) - price)
        while diff > epsilon:
            sigma = (sigma -
                     (option_price(right, s, k, t, rfr, sigma, div) - price) /
                     (option_vega(s, k, t, rfr, sigma, div) * 100))
            diff = np.abs(
                    option_price(right, s, k, t, rfr, sigma, div) - price)
        return sigma

    iv = newton_raphson(right, s, k, t, rfr, sigma, price, epsilon, div)
    return iv
```

## Answer by Mehdi Zare (score 2)

https://quant.stackexchange.com/a/49684

Expected volatility in the underlying price over the life of the option is a major component of the BSM option pricing model. When you calculate the volatility based on the current market price, you're figuring out what the market thinks the volatility would be, that's why it's called implied volatility.

So to answer your question, you can either assume a value for volatility and calculate the option price or use the option price to calculate implied volatility.