% Variance swap calculation using the replication algorithm of
% Demeterfi, Derman, Kamal, and Zou (1999).
% Example using SPX options on 2/28/2011. Linear interpolation of implied
% volatility to create a continuum of IV along a fine grid.
% This example ignores dividends on the S&P 500.
% By Fabrice Douglas Rouah

clc; clear;

% Select the maturity (t = 1 through 9).
t = 5;

% Auxiliary Function
f = @(S, Sb, T) 2/T*((S - Sb) / Sb - log(S/Sb));

% Input the SPX spot price on 2/28/2011 and define the cutoff
S  = 1327.22;
Sb = S;

% Input the implied volatility surface for SPX on 2/28/2011
IV = [ ...
0.66249	0.46257	0.35304	0.24611	0.14728	0.10517	0.10375	0.10375	0.10375	0.10375;...
0.48812	0.35515	0.28751	0.22307	0.16360	0.12388	0.11477	0.11477	0.11477	0.11477;...
0.41881	0.31711	0.26753	0.22104	0.17847	0.14578	0.12511	0.12266	0.12266	0.12266;...
0.38830	0.30064	0.25893	0.22012	0.18468	0.15602	0.13465	0.12496	0.12491	0.12491;...
0.37170	0.29163	0.25422	0.21965	0.18821	0.16231	0.14205	0.12896	0.12559	0.12555;...
0.32688	0.27052	0.24523	0.22215	0.20128	0.18354	0.16873	0.15639	0.14591	0.13695;...
0.31740	0.27327	0.25354	0.23552	0.21915	0.20442	0.19115	0.17915	0.16815	0.15797;...
0.32155	0.28460	0.26801	0.25280	0.23892	0.22632	0.21485	0.20440	0.19474	0.18575;...
0.33407	0.30398	0.29032	0.27771	0.26612	0.25551	0.24580	0.23693	0.22870	0.22103];

% Input the strikes (ATM is column 5)
Strike = [0.5 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5].*S;

%% Input the maturities and the date for 2/28/2011
Today = '02/28/11';
Mat  = {'03/28/11'; '05/31/11'; '08/29/11'; '11/28/11'; '02/28/12'; '02/28/14'; '02/29/16'; '02/28/18'; '03/01/21'};
Time = (datenum(Mat) - datenum(Today))/365;

% Input the discount factors based on the 2/28/2011 yield curve
Dates = {...
'03/03/11'; '03/05/11'; '03/09/11'; '03/16/11'; '03/23/11'; '04/02/11'; '05/02/11'; ...
'06/02/11'; '06/16/11'; '09/15/11'; '12/21/11'; '03/21/12'; '06/21/12'; '09/20/12'; ...
'12/19/12'; '03/19/13'; '06/20/13'; '09/19/13'; '12/18/13'; '03/18/14'; '03/02/15'; ...
'03/01/16'; '03/01/17'; '03/01/18'; '03/01/19'; '03/02/20'; '03/01/21'; '03/01/22'; ...
'03/01/23'; '03/02/26'; '03/03/31'; '03/03/36'; '03/01/41'; '03/01/51'; '03/01/61'};

Years = (datenum(Dates) - datenum(Today))/365;

DF = [...
0.99998500014869800 0.99997500044912600 0.99995111308640100 0.99991222912297800 0.99987334818351100 ...
0.99979198571696100 0.99954964292588800 0.99923668663870500 0.99910794722545700 0.99818605657608200 ...
0.99697717643439500 0.99544225977217400 0.99323396957980800 0.99017081381990600 0.98625297354899800 ...
0.98154157399382700 0.97598196677327000 0.96990691843316500 0.96333781311287600 0.95622588310726600 ...
0.92563237401242000 0.88771864100926600 0.84756561049366600 0.80726790504638900 0.76781566721062900 ...
0.72959238518582100 0.69284469670556100 0.65686678080626400 0.62404851186287800 0.53360097103237400 ...
0.41623864372750800 0.32729883450264000 0.25866419647973800 0.16729143399603000 0.12188702483655400];

% Extract the interest rate for the selected maturity
T  = Time(t);
df = interp1(Years,DF,T);
r  = -log(df)/T;;

%% Select and interpolate the strikes and implied volatilies
% Select the implied volatility for OTM calls and puts
CallV = IV(t,5:10);
PutV  = IV(t,1:5);

% Select OTM strikes for calls, create fine grid, and interpolate the IV
CallK = Strike(5:10);
CallKI = [CallK(1)+0.1:0.1:CallK(end)];
CallVI = interp1(CallK,CallV,CallKI,'linear');

% Select OTM strikes for puts, create fine grid, and interpolate the IV
PutK = Strike(1:5);
PutKI = [PutK(1):0.1:PutK(end)-0.1];
PutVI = interp1(PutK,PutV,PutKI,'linear');

%% Do the required calculations on calls.  
% Rename CallK and CallV for convenience.
n = length(CallVI);
K = CallKI;
V = CallVI;

for i=1:n-1
	Temp(i) = (f(K(i+1), Sb, T) - f(K(i), Sb, T)) / (K(i+1) - K(i));
	if i==1 
		CallWeight(1) = Temp(1);
	end
	CallValue(i) = BSPrice(S, K(i), r, T, V(i), 'Call');
	if i>1
		CallWeight(i) = Temp(i) - Temp(i-1);
	end;
    CallContrib(i) = CallValue(i)*CallWeight(i);
end
Pi1 = sum(CallContrib);

%% Do the calculations on puts. Flip the Vectors for Convenience
n = length(PutVI);
K = fliplr(PutKI);
V = fliplr(PutVI);

for i=1:n-1
    Temp2(i) = (f(K(i+1), Sb, T) - f(K(i), Sb, T)) / (K(i) - K(i+1));
	if i==1
		PutWeight(1) = Temp2(1);
	end;
	PutValue(i) = BSPrice(S, K(i), r, T, V(i), 'Put');
	if i>1
		PutWeight(i) = Temp2(i) - Temp2(i-1);
	end
    PutContrib(i) = PutValue(i) * PutWeight(i);
end
Pi2 = sum(PutContrib);

% Total cost of the portfolio
Pi_CP = Pi1 + Pi2;

%% Results of the replication
% Estimate of fair variance
Kvar = 2/T*(r*T - (S/Sb*exp(r*T) - 1) - log(Sb/S)) + exp(r*T)*Pi_CP;

% Estimate of fair volatility
Kvol = sqrt(Kvar)*100;

% Estimates from Numerix
Numerix = [...
17.4346541536 19.5067089310 20.9119726457 21.3518611465 21.7475857520 ...
21.6785038366 21.4961533188 21.4257801140 20.8642345295];

disp(char('For maturity on '))
disp(Mat(t))
disp(char('Matlab fair volatility'))
disp(Kvol)
disp(char('Numerix fair volatility'))
disp(Numerix(t))

