In order to understand what the implied volatility of a call option is and then calculate it, consider the Black-Scholes formula for pricing a European call option below. In this article we will focus on the mathematics behind solving for the implied volatility of a call option. If you are more interested in simply getting the number in a fast manner, this article which uses Scipy's minimize module may be more suitable.
\(Call = S_{0}N(d1) - N(d2)Ke^{-rT}\)
Where:
\( \\d1 = \frac{ln(\frac{S}{K}) + (r + \frac{\sigma^2}{2})T}{\sigma\sqrt{T}} \\ \\d2 =d1 - \sigma\sqrt{T} \)
\(S \) : current asset price
\(K\): strike price of the option
\(r\): risk free rate
\(T\) : time until option expiration
\(\sigma\): annualized volatility of the asset's returns
\(N(x) \): is the cumulative distribution function for a standard normal distribution
import numpy as np
from scipy.stats import norm
N_prime = norm.pdf
N = norm.cdf
def black_scholes_call(S, K, T, r, sigma):
'''
:param S: Asset price
:param K: Strike price
:param T: Time to maturity
:param r: risk-free rate (treasury bills)
:param sigma: volatility
:return: call price
'''
###standard black-scholes formula
d1 = (np.log(S / K) + (r + sigma ** 2 / 2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
call = S * N(d1) - N(d2)* K * np.exp(-r * T)
return call
Consider the following example to demonstrate implied volatility. You observe a stock in the market with a current price of $100, there is an option for sale, for the right to buy the stock in exactly 1 year from now for $115 , the option costs $18 and the risk free rate is 5%. Notice we know everything except the volatility(\(\sigma\)) from the equation above. Observe the example in equation form:
\(18 = 100N( \frac{ln(\frac{100}{115}) + (0.05 + \frac{\sigma^2}{2})1}{\sigma\sqrt{1}}) - N( \frac{ln(\frac{100}{115}) + (0.05 + \frac{\sigma^2}{2})1}{\sigma\sqrt{1}}-\sigma\sqrt{1})115e^{-0.05*1}\)
The \(\sigma\) parameter from above is the volatility at which the Black-Scholes formula would return a value of $18. So essentially once we solve for \(\sigma\) in the equation above we have the implied volatility of the option price. Since the formula above cannot be solved explicitly we must resort to iterative measures.
Brute Force
Let's use a brute force method below to demonstrate what iterative methods actually are. The numpy array of volatility candidates contains approximately 40,000 possible solutions to the equation described above, below we simply iterate through them all in order to find that one that minimizes the absolute difference between the observed price and the Black-Scholes price, or in math terms we are finding the root of the function.
volatility_candidates = np.arange(0.01,4,0.0001)
price_differences = np.zeros_like(volatility_candidates)
observed_price = 18
S = 100
K = 115
r = 0.05
T = 1
for i in range(len(volatility_candidates)):
candidate = volatility_candidates[i]
price_differences[i] = observed_price - black_scholes_call(S, K , T, r, candidate)
idx = np.argmin(abs(price_differences))
implied_volatility = volatility_candidates[idx]
print('Implied volatility for option is:', implied_volatility)
Implied volatility for option is: 0.5427999999999968
You can verify the iterative algorithm worked by plugging the implied volatility number back into the Black-Scholes formula.
S = 100
K = 115
r = 0.05
T = 1
sigma = implied_volatility
price = black_scholes_call(S, K, T, r, sigma)
print(price)
17.9983104365487
Notice that the value we get above is slightly different from the $18 we expected. This is inherent feature of iterative methods, think of the 0.0001 in volatility_candidates = np.arange(0.01,4,0.0001) as a tolerance for error , changing this value up or down will increase or decrease the error respectively.
Whilst the method described above works pretty well, you may have noticed that it takes quite a while to finish. Imagine having say 1000 options for which you wanted to calculate the implied volatility. It took roughly 5 seconds for the above script to finish, multiplying that by 1000 equates to well over an hour of waiting.
Also notice the search area we designated for our brute force method was bounded between 1% and 400%, there is no reason why an option wouldn't be priced higher than that which would lead trouble!
Newton Raphson Algorithm
The Newton Raphson method is a widely used algorithm for calculating the implied volatility of an option. The steps for the generic algorithm are as follows:
1) Define our function as \(f(x)\) for which we want to solve \(f(x) = 0 \)
2) Choose an initial guess
3) Iterate as follows: \(x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\)
4) Break if: \(|f(x_n)| < \epsilon\) , here \(\epsilon\) is the tolerance for error described in the brute force method
Applying the steps above to our problem of solving implied volatility.
1) \(f(\sigma) = V_{BS_\sigma} - V_{market}\)
2) Initial guess \(\sigma_0 = 0.30\) (Chosen for convenience)
3) Iterate as follows \(\sigma_{n+1} = \sigma_n - \frac{V_{BS_\sigma} - V_{market}}{\frac{\partial V_{BS_\sigma}}{\partial\sigma}}\)
4) if \(|V_{BS_\sigma} - V_{market} | < \epsilon\) , return \(\sigma_n\)
In order to code this in Python we first need to make a function that computes the partial derivative of the value of a call option with respect to volatility. This partial derivative is commonly referred to as vega.
\(vega = S\sqrt{T} N'(d1)\)
d1 in vega is the same d1 = \(\frac{ln(\frac{S}{K}) + (r + \frac{\sigma^2}{2})T}{\sigma\sqrt{T}}\) that we calculated in the Black-Scholes formula and \(N'\) is the probability density function for a standard normal.
The code snippet below implements vega:
N_prime = norm.pdf
def vega(S, K, T, r, sigma):
'''
:param S: Asset price
:param K: Strike price
:param T: Time to Maturity
:param r: risk-free rate (treasury bills)
:param sigma: volatility
:return: partial derivative w.r.t volatility
'''
### calculating d1 from black scholes
d1 = (np.log(S / K) + (r + sigma ** 2 / 2) * T) / (sigma * np.sqrt(T))
vega = S * np.sqrt(T) * N_prime(d1)
return vega
Now that we have a function to calculate vega let's create a function to implement Newton's method on our example security and compare it to the value we got using the brute force approach.
def implied_volatility_call(C, S, K, T, r, tol=0.0001,
max_iterations=100):
'''
:param C: Observed call price
:param S: Asset price
:param K: Strike Price
:param T: Time to Maturity
:param r: riskfree rate
:param tol: error tolerance in result
:param max_iterations: max iterations to update vol
:return: implied volatility in percent
'''
### assigning initial volatility estimate for input in Newton_rap procedure
sigma = 0.3
for i in range(max_iterations):
### calculate difference between blackscholes price and market price with
### iteratively updated volality estimate
diff = black_scholes_call(S, K, T, r, sigma) - C
###break if difference is less than specified tolerance level
if abs(diff) < tol:
print(f'found on {i}th iteration')
print(f'difference is equal to {diff}')
break
### use newton rapshon to update the estimate
sigma = sigma - diff / vega(S, K, T, r, sigma)
return sigma
Recall the example security we used for the brute force method:
You observe a stock in the market with a current price of $100, there is an option for sale for the right to buy the stock in exactly 1 year from now for $115 , the option costs $18 and the risk free rate rate is 5%.
observed_price = 18
S = 100
K = 115
T = 1
r = 0.05
imp_vol = implied_volatility_call(observed_price, S, K, T, r)
print('Implied volatility using Newton Rapshon is: ',imp_vol)
Which returns:
found on 2th iteration
difference is equal to -7.274655111189077e-06
Implied volatility using Newton Rapshon is: 0.5428424065162359
Compare this number to the one we calculated using the brute force approach, they are effectively identical. You can verify this number is correct by plugging the imp_vol back into the Black-Scholes formula which returns approximately 18.
It took approximately 0.002 seconds when I ran the snippet above, comparing this to 5 seconds it took when using the brute force approach. Consider again calculating the implied volatility for 1000 options using both approaches; using Newton Raphson we would be finished in 2 seconds in comparison to well over an hour for brute force approach. A significant increase in speed!
There is a useful Discussion on Stack Exchange where a contributor gives a link to a paper by Brenner and Subrahmanyam (1988) which provides a method for calculating the initial guess which is as follows:
\(\sigma \approx \sqrt{\cfrac{2\pi}{T}} . \cfrac{C}{S}\)
You may want to try replacing the initial estimate of 0.30 with the formula above.
Wikipedia also has a great page on Newton Raphson with much more detail if you are interested in learning more.
Full script for calculating Newton Raphson
import numpy as np
from scipy.stats import norm
N_prime = norm.pdf
N = norm.cdf
def black_scholes_call(S, K, T, r, sigma):
'''
:param S: Asset price
:param K: Strike price
:param T: Time to maturity
:param r: risk-free rate (treasury bills)
:param sigma: volatility
:return: call price
'''
###standard black-scholes formula
d1 = (np.log(S / K) + (r + sigma ** 2 / 2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
call = S * N(d1) - N(d2)* K * np.exp(-r * T)
return call
def vega(S, K, T, r, sigma):
'''
:param S: Asset price
:param K: Strike price
:param T: Time to Maturity
:param r: risk-free rate (treasury bills)
:param sigma: volatility
:return: partial derivative w.r.t volatility
'''
### calculating d1 from black scholes
d1 = (np.log(S / K) + (r + sigma ** 2 / 2) * T) / (sigma * np.sqrt(T))
#see hull derivatives chapter on greeks for reference
vega = S * N_prime(d1) * np.sqrt(T)
return vega
def implied_volatility_call(C, S, K, T, r, tol=0.0001,
max_iterations=100):
'''
:param C: Observed call price
:param S: Asset price
:param K: Strike Price
:param T: Time to Maturity
:param r: riskfree rate
:param tol: error tolerance in result
:param max_iterations: max iterations to update vol
:return: implied volatility in percent
'''
### assigning initial volatility estimate for input in Newton_rap procedure
sigma = 0.3
for i in range(max_iterations):
### calculate difference between blackscholes price and market price with
### iteratively updated volality estimate
diff = black_scholes_call(S, K, T, r, sigma) - C
###break if difference is less than specified tolerance level
if abs(diff) < tol:
print(f'found on {i}th iteration')
print(f'difference is equal to {diff}')
break
### use newton rapshon to update the estimate
sigma = sigma - diff / vega(S, K, T, r, sigma)
return sigma