## Introduction to volatility models with Matlab (Implied volatility)

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

1. Compute using (3).
2. Compute theoretical value of an option using (2).
3. Set initial values(initial guess) for secant method, that is 4. and .
5. These temporary variables will change in the secant method loop. Keep in mind that the iterator in (4) of the secant loop is (not !).
6. Compute approximation error(last product on the left side of (4)).
7. Compute new using (4).
8. Compute new using (2)
9. If error < accuracy threshold then last is our from secant method alias implied volatility. Proceed to the next .
10. 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);

% 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 ). 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.