Implied volatility (IV) is the volatility of an asset derived from changes in value of corresponding option in such way that if we input IV into option pricing model, it will return theoretical value equal to the current option value. Contrary to historical volatility, IV is the volatility forecast for price of the underlying asset from current time to option expiration . Let’s quickly introduce the Black-Scholes option pricing model. There are multiple ways to derive BSE and you can review them here. The brief solution of BSE can be found here. So the familiar BSE form is

(1)

where

– option value,

– underlying asset spot price,

– risk-free interest,

– current time, date of expiry, time to expiry,

– volatility (in our case implied volatility).

The solution of (1) for value of an call or put option at time can be expressed as

(2)

where

is strike price of an option and is value of normal CDF at particular points . Recall that normal CDF is written

Do you see how deeply is buried? We see that and and since there doesn’t exist the closed form solution of the can’t be expressed analytically. That’s why we use numerical methods and approximation to evaluate implied volatilities. Familiar widely used numerical approaches are Newton-Raphson or secant method. Approximation of IV can be estimated by Bharadia-Christopher-Salkin(1996) model or Corado-Miller(1996) model. Among the latest and most effective estimators we can sort Hallerbach(2004), Jaeckel(2006) and Jaeckel(2013).

For purpose of this article let’s show how to implement just the most basic B-C-S model and secant method.

The B-C-S model is derived as

(3)

where

– time to option expiration,

– current option value,

– underlying asset spot price,

– option strike price( is discounted strike price).

The secant numerical method is given by

(4)

where

is theoretical option value from (2).

Algorithm itself works as follows(suppose computation for one specific ):

- Compute using (3).
- Compute theoretical value of an option using (2).
- Set initial values(initial guess) for secant method, that is
- and .
- These temporary variables will change in the secant method loop. Keep in mind that the iterator in (4) of the secant loop is (not !).
- Compute approximation error(last product on the left side of (4)).
- Compute new using (4).
- Compute new using (2)
- If error < accuracy threshold then last is our from secant method alias implied volatility. Proceed to the next .
- If error is > accuracy threshold go back to step 6 (proceed to the next ) and don’t forget you need to recompute new values. Use from last secant iteration (note that initial values from step 3. and 4. will not be used in 2nd secant iteration). On every new iteration, variables become variables and new variables are computed.

I will use historical EOD option chain data for BAC stock.

*Fig.3 Example option chain .csv file for BAC stock.*

If you are wading through my intro to volatility thoroughly you may want to download the example option chain. Option data are very expensive and tough to store and manipulate so it’s nearly impossible to find them free. So now we need to prepare data for IV calculation as follows:

- Use optionparse() function for importing the .csv option chain into matlab table.

table = optionparse('path/to/datafile','ivolatility','equity');

- Now we extract dates , adjusted stock spot prices and we choose call options with strike , time to expiry and the ‘risk-free’ interest rate be which was T-Bill rate in the early 2015. Let me point out that in practice “constant” will not be constant because we will not find such options at every time which have exactly half a year expiration. Only finite number of options is offered every . So our will be variable. At every I choose the first option with and . Following code finds forementioned options.

% find dates [~,ia,ib] = unique(table.data.date,'first','legacy'); t = table.data.date(ia); S = table.data.adjusted_stock_close_price(ia); % time to expiration in N/252 fraction, for simplicity I don't consider holidays, in practice you should Tte = yearfrac(table.data.date,table.data.expiration,13); % search options with expiration T_t and strike K T_t_const = 0.5; K = 15; % expiration filter ic_start = ia + accumarray(ib,Tte,[],@(x) find(x>=T_t_const,1)) - 1; ic_end = ia + accumarray(ib,Tte,[],@(x) find(x<=T_t_const+0.5,1,'last')) - 1; idExpiration = zeros(size(Tte)); % strike filter id = zeros(length(ic_start),1); for i=1:length(ic_start) temp = find(table.data.strike(ic_start(i):ic_end(i)) == K,1); if ~isempty(temp) id(i,1) = find(table.data.strike(ic_start(i):ic_end(i)) == K,1); else id(i,1) = NaN; end end % uncomment proper option idTK = ic_start + id - 1; % call options % idTK = ic_start + id; % put options % time to expiry at each time t T_t_var = Tte(idTK); % let the value of the option be (ask+bid)/2. V is vector with filtered % options V = (table.data.ask(idTK) + table.data.bid(idTK))/2;

- Our variables and constants are ready. I will demonstrate IV computation for .

1. Calculate initial IV estimate by B-C-S model (3).

t = 1; r = 0.001; delta = (1/2)*(S(t) - K*exp(-r*T_t_var(t))); % (3) sigma_BCS(t,1) = sqrt(2*pi-T_t_var(t)) * (V(t)-delta)/(S(t)-delta); % (3)

2. Compute theoretical value of an option using (2).

d1 = (log(S(t)/K) + (r+1/2*sigma_BCS(t,1)^2)*(T_t_var(t))) / sigma_BCS(t,1)*sqrt(T_t_var(t)); % (2) d2 = d1 - sigma_BCS(t)*T_t_var(t); BS(t) = S(t)*cdf('normal',d1,0,1) - K*exp(-r*T_t_var(t))*cdf('normal',d2,0,1); % (5) for call option

3. Set initial values(initial guess) for secant method, that is . Note that I smuggled in the initial error value and secant method accuracy.

err_thr = 1e-4; % secant method error tolerance(accuracy) err = 1; % initial error value temp_sigma(1:2,1) = [0,sigma_BCS(t,1)];

4. and .

temp_BS(1:2,1) = [0,BS(t,1)];

5. Secant method, steps described in the code.

while err >= err_thr % error update - step 6 err = (temp_BS(2,1) - V(t)) * ... ((temp_sigma(2,1) - temp_sigma(1,1)) / (temp_BS(2,1) - temp_BS(1,1))); % (7) error corresponding to last computed sigma % sigma_{i+1} update - step 7 temp_sigma(3,1) = temp_sigma(2,1) - err; % new sigma from secant iteration temp_sigma(1) = []; % old sigma no more needed in next loop % BS(sigma_{i+1}) update - step 8 d1 = (log(S(t)/K) + (r+1/2*temp_sigma(2,1)^2)*(T_t_var(t))) / temp_sigma(2,1)*sqrt(T_t_var(t)); % (2) d2 = d1 - temp_sigma(2,1)*T_t_var(t); % (2) temp_BS(3,1) = S(t)*cdf('normal',d1,0,1) - K*exp(-r*T_t_var(t))*cdf('normal',d2,0,1); % (2) for call option temp_BS(1) = []; % old BS(sigma) no more needed in next loop end

9. If error < accuracy threshold then last is our .

sigma_secant(t,1) = temp_sigma(2,1);

10. If error > accuracy threshold go to step 6 (next ).

You can download complete working BCS_secant() function here.

*Fig.4 B-C-S-secant computed daily implied volatility for Bank of America stock in 2015
*

We can slowly move to the next concept of practical importance which ensues directly from implied volatility calculation. We have computed one IV value for every and we supposed that where were fixed (“fixed” for because cannot be strictly constant). Clearly we can release this assumption and let the be variables. For every we now have not just one point but set of points which constitutes the *implied volatility surface.* If we now fix just at some time in we obtain values constituting *implied volatility smile or implied volatility skew*, if we fix just at some time we obtain *implied volatility term structure*.

Multitudes of volatility indexes (composites) are constructed using various setups in IV computation. Some of them average IV’s of stocks in indexes such a S&P500, DJIA, NASDAQ …, some of them average similar ETF’s etc. As was mentioned before such indexes measure expectation of volatility. For example have a look at the CBOE volatility indexes overview.