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

1. 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?

2. 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?

• You are given $C_M$ (the market price of the option) and you try various values of $\sigma$ which give corresponding $C_{BS}$ values (the result of plugging sigma into the BS formula). When the absolute value of the difference between $C_M$ and the Black Scholes price $C_{BS}$ is sufficiently small (say <0.001) you have found the correct IV. – Alex C Nov 12 '19 at 15:38
• @AlexC Part of my goal is to ultimately calculate the market value of the option, so it wouldn’t be a given in my situation. Is there a way to derive it without trial and error comparison against another source? – pez Nov 12 '19 at 15:40
• @Pez: Hi. if you don't have the market value of the option, then you can't calculate the implied volatility. you need the market value. if you don't have the market value, then you generally price the option by constructing some kind of volatility estimate and plugging it in. So, in the latter case, if two people have different volatility estimates, then they will each obtain different market values for the option. – mark leeds Nov 12 '19 at 15:57
• @markleeds thanks, so the only option is to take current market value of the option from an external source (like CBOE or a brokerage?) and reverse engineer the IV? I don’t have much experience in this field so it seems odd to me there’s no way to derive without already having the answer to compare against. Almost seems like a Catch 22! – pez Nov 12 '19 at 17:21
• Emmanuel Derman wrote a paper a number of years back that I used to do the calculation in a past position...it basically gave a step by step implementation, and I don't remember the implementation being all that challenging. I'd do some searching to see if I could find a copy it if I were you. – Chris Nov 13 '19 at 6:21

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 = 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


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.