Wednesday, February 4, 2009

Replicate VIX (CBOE Volatility Index) --MATLAB Implementation

This post is just to demonstrate how to replicate the Calculations behind the CBOE Volatility Index, Commonly called VIX. It is also commonly thought of investor gauge of fear. One can read more about at www.cboe.com

The Implementation below follows the methodology as illustrated in the following white paper:

http://www.cboe.com/micro/vix/vixwhite.pdf

M FILE is here:

http://docs.google.com/Doc?id=ddb2j6dw_12fjk57bfx




CODE Published Here:

%% ReplicateVIX
% Near-Term
% Strike Call Put
Data{1} = [775 125.48 0.11;...
800 100.79 0.41;...
825 76.70 1.30;...
850 54.01 3.60;...
875 34.05 8.64;...
900 18.41 17.98;...
925 8.07 32.63;...
950 2.68 52.33;...
975 0.62 75.16;...
1000 0.09 99.61;...
1025 0.01 124.52];
% Next Term
% Strike Call Put
Data{2} = [775 128.78 2.72;...
800 105.85 4.76;...
825 84.14 8.01;...
850 64.13 12.97;...
875 46.38 20.18;...
900 31.40 30.17;...
925 19.57 43.31;...
950 11.00 59.70;...
975 5.43 79.10;...
1000 2.28 100.91;...
1025 0.78 124.38];

%Time_To_Maturity
TM = [16;44];
%Risk_Free_Rate
Rf = 1.1625/100; %Per Annum
% Difference between Calls and Puts (Absolute Value)
DF{1} = abs(Data{1}(:,2) - Data{1}(:,3));
DF{2} = abs(Data{2}(:,2) - Data{2}(:,3));
% Current Time
CT = '08:30:00';
% FInd Hour, Minute, Second from the time using datevec function
[Year, Month, Day, Hour, Minute, Second] = datevec(CT);
%In Years
%1440 is the number of minutes in a day and 510 is the number of minutes to 8:30 AM which
%is the time the option expires on its expiration date
NumYears(1) = [1440 - (Hour * 60 + Minute + Second/60) + 510]/(1440 * 365) + [(TM(1) - 2)/365];
NumYears(2) = [1440 - (Hour * 60 + Minute + Second/60) + 510]/(1440 * 365) + [(TM(2) - 2)/365];
% In days
NumDays = NumYears .* 365;
% Find the minimum of the difference in Call and Put
% Prices and Get the corresponding Strike Price.
ATM(1,:) = Data{1}(find(DF{1}==min(DF{1})),:);
ATM(2,:) = Data{2}(find(DF{2}==min(DF{2})),:);
% Calculate Forward Price Level and Referential Strike
% Application of PUT CALL Parity
Level = ATM(:,1) + exp(Rf*NumYears(:)) .* (ATM(:,2) - ATM(:,3));
%Reference Strike
for i = 1:2
Strike = ATM(i,1);
if(ATM(i,2)>=ATM(i,3))
Ref_Strike(i)=ATM(i,1);
else
Ref_Strike(i) = Data{i}(find(Data{i}(:,1) < ATM(i,1),1,'last'),1);
end
% Differences of Strikes
Temp = diff(Data{i}(:,1));
Delta_Strike{i} = [Temp(1);Temp];
%{
? If the strike is above the “reference strike” , use the call price
? If the strike is below the “reference strike” , use the put price
? If the strike equals the “reference strike” , use the average of the call
and put prices
%}
cpval= zeros(size(Data{i},1),1);
cid = find(Data{i}(:,1) > Ref_Strike(i));
cpval(cid) = Data{i}(cid,2);
pid = find(Data{i}(:,1) < Ref_Strike(i));
cpval(pid) = Data{i}(pid,3);
Aid = find(Data{i}(:,1) == Ref_Strike(i));
cpval(Aid) = (Data{i}(Aid,2) + Data{i}(Aid,3))/2;

% Now do the math as given in the paper vixwhite.pdf

vix{i} = Delta_Strike{i} * exp(Rf*NumYears(i)) .* cpval ./(Data{i}(:,1).^2);
Var(i) = (2/NumYears(i)) * sum(vix{i}) - ((Level(i)/Ref_Strike(i) - 1).^2)/NumYears(i);

% Center the data to 30 days
if(i==1)
Term(i) = NumYears(i) * Var(i) * ((NumDays(i+1)-30)/(NumDays(i+1)-NumDays(i)));
elseif(i==2)
Term(i) = NumYears(i) * Var(i) * ((-NumDays(i-1)+30)/(NumDays(i)-NumDays(i-1)));
end


end

% Final Vix Calculation
VIX = sqrt(sum(Term) * 365/30) * 100

7 comments:

A said...

Hi,

First, I wanted to let you know that I've been using the Matlab2IB software and found it helpful with my trading.

I went through the same exercise above to recreate the VIX formula as well as the "business" VIX - see vix term structure on the cfe. I've confirmed that the results of your code match my code for the VIX. However, I'm having some difficulty recreating the business vix. This should be trivial but it seems there is something I'm missing. I was hoping you could take a look at your code and post a business vix version. I believe it's simply a matter of changing the calendar time to business time in a few places and changing the last line to:

BUSINESS_VIX = sqrt(sum((Level(1) * Var(2))/NumYears(1))) * 100

(for the nearest point of the vix term structure)

Or I can post my code but it may be easier for all to follow if we continue with your example. (The business vix is of particular interest because the cfe provides historical data for the vix term structure using the business vix, rather than the calendar vix, convention).

Many thanks.

TradingwithMatlab said...

Hello,

Thanks. I will take a look and will post the code again.

Thanks,

Anonymous said...

Hi,
Your code has some bugs. More specifically:

Var(i) = (2/NumYears(i)) * sum(vix{i}) - ((Level(i)/Ref_Strike(i) - 1).^2)/NumYears(i);

Term(i) = NumYears(i) * Var(i) * ((NumDays(i+1)-30)/(NumDays(i+1)-NumDays(i)));

TradingwithMatlab said...

I dont see any bug in the code. I tested it and got the same numbers as in the document. Can you give me the error that you got using the code??

Anonymous said...

Hi,
please check the algorithm provided in the VIX White Paper. The volatility (variance) formula according to the Paper is:

σ^2 = 2/Τ*Sum(ΔΚ/Κ^2)*exp(RT)*Q(K) - (F/K0 - 1)^2/T

But your code is:

σ^2 = 2/F*Sum(ΔΚ/Κ^2)*exp(RT)*Q(K) - (F/K0 - 1)^2/F, where F = Level(i) in your code.

However, as you said, the result is the same in both cases. Homework: Can you see why?

TradingwithMatlab said...

Yes,

Because when you multiply Var(i) and Level(i) in the term calculation, the numerator and denominator cancel each other out.

So, Level(i) or NumYears(i) will yield the same answer. I will try to edit my post later.

Thanks

Anonymous said...

Hi I am having some problems with the mat lab formula. Do you take the average of the call and puts for Data 1 and Data?