Volatility smile modelling: two-dimensional linear regression demo

A volatility smile modeling is one of the actual problems in finance mathematics. Here is an example how to make a regression model of the volatility smile using historical data of an option.

Contents

Load the historical data

rawdata = dlmread('option_AP04.csv'); % read history file
time  = rawdata( :,1 );         % extract time values
maturity = -1/365*time;         % convert time
options= rawdata( :,2:7 );      % extract price of the option for each strike
price  = rawdata( :,8 );        % extract price of the stock
clear('rawdata', 'time');       % do not need anymore
strike = [13.50 13.00 12.50 12.00 11.50 11.00];

Calculate the impied volatility using the Black-Sholes’ model

for i=1:length(strike)
    % doc blsimpv for more information
    impvol(:,i) = blsimpv(price, strike(i), 0.075, maturity, options(:,i), 1.5, 0, [], true);
end

Plot the source data

figure; subplot(1,2,1);
% Tri = delaunay(x(:,1),x(:,2));
% trimesh(Tri,x(:,1),x(:,2),y,'EdgeColor','black');
% camlight left; lighting phong; alpha(.8);
[X,Y] = meshgrid(maturity, strike);
mesh(X,Y,impvol');
scalevec = [min(maturity) max(maturity) min(strike) max(strike) min(impvol(:)) max(impvol(:))];
axis(scalevec);
xlabel('maturity');
ylabel('strike');
zlabel('implied volatility');

Make the regression Maturity and Strike to ImpVol

% we have:  maturity [T,1]  T = 64
%           strike   [1,S]  S = 6
%           impvol   [T,S]
% let x be independent variable and y be dependeny variable

[X,Y] = meshgrid(maturity,strike);
x = [X(:),Y(:)]; % vectorize and concatenate
clear('X', 'Y');
impvolT = impvol;
y = impvolT(:);
% now we get x and y

Check it with a simple example

% mat = [1; 2; 3; 4]
% str = [22 33]
% imv= [122 133;
%      222 233;
%      322 333;
%      422 433]
% ivT = iv';
%
% [X,Y] = meshgrid(mat,str)
% x = [X(:), Y(:)]
% y = imvT(:)

Check the sample set for NaNs

idx = [ find(isnan(y)); ...
        find(isnan(sum(x,2)))]; % note that val+NaN = NaN
y(idx)   = [];
x(idx,:) = [];

Assign the model

the model is % impvol = w(1) * maturity^2 + w(2) * maturity + w(3) * strike + w(4); Making of more complex and more adequate
model is outside of this example. You can try MVR Composer to make this kind of models inductively.

% weak code
% A = [];
% for s = strike
%     S = [maturity.^2, maturity, repmat([s^2 s 1], [length(maturity) 1])];
%     A = [A; S];
% end

% better code
f = inline('[x(:,1).^2, x(:,1), x(:,2) , x(:,2).^0 ]','x');
A = f(x);
w = inv(A'*A)*A'*y;
y1 = A*w;

Plot the result

subplot(1,2,2);
Tri = delaunay(x(:,1),x(:,2));
trimesh(Tri,x(:,1),x(:,2),y1,'EdgeColor','black');
camlight left; lighting phong; alpha(.8);
xlabel('maturity');
ylabel('strike');
zlabel('implied volatility');
axis(scalevec);