Thursday, July 13, 2017

Calculating Implied Volatility from an Option Price

Nowadays, with the ubiquity of computers and information, many traders wish to know how the sausage is made and work out for themselves some of the numbers appearing on the screen of their trading platforms.  This is often a simple matter of looking up the formulas involved and plugging in the numbers.

When calculating the numbers pertaining to options, however, we run into an issue with implied volatility.  The Black-Scholes model tells us what an option should be worth given its strike price, the risk-free interest rate, the remaining time until expiration, the stock's price, and the implied volatility.  For example, the price of a call, $C$, is given by,

$$C = \Phi(d_1) S - \Phi(d_2) K e^{-r t},$$

where $S$ is the price of the stock, $K$ is the strike price, $r$ is the annualized risk-free rate, and $t$ is the remaining time to expiration expressed in years.  $\Phi$ is the normal cumulative distribution function, and $d_1$ and $d_2$ are given by,
$$d_1 = \frac{1}{\sigma \sqrt{t}} \left[ \ln\left(\frac{S}{K}\right) + \left(r + \frac{\sigma.^2}{2}\right) t\right],$$
and
$$d_2 = d_1 - \sigma \sqrt{t},$$
respectively.  In the above expressions for $d_1$ and $d_2$,  $\sigma$ is the implied volatility.

If we wish to calculate $\sigma$, we run into an issue.  All the variables in the above equations are known, and we can get the call price directly from the option chain, but we are incapable of isolating $\sigma$ algebraically.  Instead, we will have to turn to numerical methods to calculate the implied volatility.

Recall from high school algebra that if $y = f(x)$, the value of $x$ for which $f(x) = 0$ is called the root of the function $f$.  We will make use of a root finding algorithm to find our volatility, $\sigma$.  Our function will be the theoretical call price from the Black-Scholes model minus the known option price.  We will insist that equals zero and find the value of $\sigma$ that makes it so.  Explicitly, we want,
$$C(\sigma) - C_0 = 0,$$
where $C_0$ is the call price from the option chain.

For the sake if this article, we will assume the underlying stock is trading for 100, the strike price is 105, there are 30 days until the contract expires, and the risk-free rate is 1%.  We will plot $C(\sigma) - C_0$ for values of implied volatility from zero to one below just to get a sense of how this function behaves.
Figure 1.  We wish to find the value of $\sigma$ where $C(\sigma) - C_0$ is zero.

We wish to find the value of $\sigma$ where the function $C(\sigma) - C_0 = 0$.  From Fig 1., we can see this happens around $\sigma = 0.38$.  We will use an iterative method for finding an approximation to the root developed by Isaac Newton and Joseph Raphson.  The method requires an initial guess to start the process.  We will denote our successive approximations of the root as $\sigma_i, i = 0,1,2,...$ and will start with a guess of $\sigma_0 = 0.5$.

The Newton-Raphson method calculates the slope of the tangent line evaluated at the first iteration point.  It then uses the point where that tangent line crosses zero as the starting point for the next iteration.  This is shown graphically in Fig. 2.

Figure 2.  Plot of the line tangent to our function at $\sigma = 0.5$.  The intercept of this line with the x-axis is used as the starting point for the next iteration.  We can see after only one iteration it is relatively close to the actual root

The process then repeats until our test value for $\sigma$ is as close to zero as we like.  Fig. 3 shows the second iteration.  Note that after only two iterations, we are very close to the root.
Figure 3.  Tangent line for the second iteration

To implement this technique, we need to be able to calculate the derivative of the function of interest.  One of the nice features of the Black-Scholes model is that it has a closed form that can be easily differentiated.  Traders refer to the derivative of the option price with respect to volatility as vega and denote it with the Greek letter $\nu$.  Vega is given by,
$$\nu = \frac{\partial C}{\partial\sigma} = S \phi(d_1)\sqrt{t},$$
where $\phi$ is the normal probability density function.  With the above equations, we have enough information to implement a program to calculate the implied volatility of an option.  We will use Python for this exercise because it is a popular, freely available programming language that has a fairly extensive math and statistics libraries.

We will make use of the scypi.stats library as well as specific functions in the math package.  After importing these, we need to write a couple of helper functions to implement the Black-Scholes model for call option prices.  The first function calculates the values of $d_1$ and $d_2$.  The second function actually calculates the theoretical price of the call.

from scipy.stats import norm
from math import sqrt, exp, log, pi

def d(sigma, S, K, r, t):
    d1 = 1 / (sigma * sqrt(t)) * ( log(S/K) + (r + sigma**2/2) * t)
    d2 = d1 - sigma * sqrt(t)
    return d1, d2

def call_price(sigma, S, K, r, t, d1, d2):
    C = norm.cdf(d1) * S - norm.cdf(d2) * K * exp(-r * t)
    return C

Next we enter in the known option information as well as our initial guess for $\sigma$.

#  S  = Stock price
#  K  = strike
#  C  = price of call as predicted by Black-Scholes model
#  r  = risk-free interest rate
#  t  = time to expiration expressed in years
#  C0 = price of call option from option chain

S = 100.0
K = 105.0
r = 0.01
t = 30.0/365
C0 = 2.30

#  We need a starting guess for the implied volatility.  We chose 0.5
#  arbitrarily.
vol = 0.5

Next we define some variables needed for bookkeeping.  These will include a variable to count the number of iterations along with a maximum number of iterations.  We will use these to make sure our program doesn't get stuck in an infinite loop.  We will also set a tolerance where we will terminate the iterations when our test root gives a value that is close enough to zero to satisfy our needs.

epsilon = 1.0          #  Define variable to check stopping conditions
abstol = 1e-4          #  Stop calculation when abs(epsilon) < this number

i = 0                  #  Variable to count number of iterations
max_iter = 1e3         #  Max number of iterations before aborting

Lastly, we will do the iteration and print out the results.

while epsilon > abstol:
    #  if-statement to avoid getting stuck in an infinite loop.
    if i > max_iter:
        print 'Program failed to find a root.  Exiting.'
        break

    i = i + 1
    orig = vol
    d1, d2 = d(vol, S, K, r, t)
    function_value = call_price(vol, S, K, r, t, d1, d2) - C0
    vega = S * norm.pdf(d1) * sqrt(t)
    vol = -function_value/vega + vol
    epsilon = abs(function_value)

print 'Implied volatility = ',  vol
print 'Code required', i, 'iterations.'

Running this produces the result,

Implied volatility =  0.368856324914
Code required 3 iterations.

We see the program converges quickly taking only three iterations to find a value within our tolerance.  The above technique can be used for puts by substituting in the Black-Scholes formula for put prices, of course.   So that's it.  Pretty simple.

Newton's method works well when the function and its derivative are well-behaved.  In our case, the function only has one root, so we don't have to worry about the possibility of having several mathematically allowable answers and setting up our code and initial guess to converge to the proper root or deciding which of several results is the right answer to our problem.

Here is the above code in its entirety.


#!/usr/bin/python

from scipy.stats import norm
from math import sqrt, exp, log, pi

def d(sigma, S, K, r, t):
    d1 = 1 / (sigma * sqrt(t)) * ( log(S/K) + (r + sigma**2/2) * t)
    d2 = d1 - sigma * sqrt(t)
    return d1, d2

def call_price(sigma, S, K, r, t, d1, d2):
    C = norm.cdf(d1) * S - norm.cdf(d2) * K * exp(-r * t)
    return C

#  S  = spot
#  K  = strike
#  C  = price of call as predicted by Black-Scholes model
#  r  = risk-free interest rate
#  t  = time to expiration expressed in years
#  C0 = price of call option from option chain

S = 100.0
K = 105.0
r = 0.01
t = 30.0/365
C0 = 2.30

#  We need a starting guess for the implied volatility.  We chose 0.5
#  arbitrarily.
vol = 0.5

epsilon = 1.0  #  Define variable to check stopping conditions
abstol = 1e-4  #  Stop calculation when abs(epsilon) < this number

i = 0   #  Variable to count number of iterations
max_iter = 1e3  #  Max number of iterations before aborting

while epsilon > abstol:
    #  if-statement to avoid getting stuck in an infinite loop.
    if i > max_iter:
        break

    i = i + 1
    orig = vol
    d1, d2 = d(vol, S, K, r, t)
    function_value = call_price(vol, S, K, r, t, d1, d2) - C0
    vega = S * norm.pdf(d1) * sqrt(t)
    vol = -function_value/vega + vol
    epsilon = abs(function_value)

print 'Implied volatility = ',  vol
print 'Code required', i, 'iterations.'

For further reading see:

Wikipedia article on the Black-Scholes model
Wikipedia article on Newton's method

Updated 11/12/2018:  Fixed text which said the strike price was 155 with a stock price of 150.  The code uses a strike of 105 and a price of 100.  Sorry for not catching this sooner.

15 comments:

  1. Great work Kevin.I have tried searching the program and kind find anything better.Thank you

    ReplyDelete
  2. Thanks for your work - great job.

    Allow me one remark on the example. You're not using the same values for S & K in the code as in the text. This might be a little bit confusing.

    ReplyDelete
  3. How will we calculate implied volatality for PUT options??

    ReplyDelete
    Replies
    1. The first equation in this article is for a call. Just substitute the equation for a put option in its place. The wikipedia article on Black-Scholes gives the put formula along with the call. I just chose to do one type of option here for simplicity.

      Delete
  4. Hi Kevin, Thank you for your work.
    how can we calibrate the SVI parameters P={a,b,p,m,s} to best fit the volatility?
    which can be accomplished as a simple least square optimization.

    Thank you in advance.

    ReplyDelete
    Replies
    1. You need a function that returns the volatility as a function of the parameters you want to fit to. Something like sigma = F(a, b, p, m, s). Then it's just a matter of doing a least squares minimization. You'll need actual volatility points though. You'd minimize sum( [F - actual vol]^2). Scipy has a curve_fit function that makes this fairly easy.

      Delete
  5. Sir is there any way I can get in touch with you. I am trying to calculate implied volatility of puts in excel but not able to do it because I am not programmer by profession hence needed your help

    ReplyDelete
  6. this is the best article on this topic (where i was stuck for a pretty good time) i could find on the internet! adding to this was your youtube video! Thank you for the help

    ReplyDelete
  7. Thank you very much Kevin ji for your kind work

    ReplyDelete