% Calculates local volatility from implied volatility using the mapping
% described in Jim Gatheral's book.  Uses implied volatilities from
% Table III.4.4 in Carole Alexander's book "Market Risk Analysis" (Vol 3).

clc; clear;

% Input the calls, strikes, and maturities.
data = xlsread('Table_III_4_4_Carol_Alexander', 'Sheet1', 'A4:M25')';

% Extract market call prices, strikes and maturities.
calls = data(2:end,2:end);
K = data(1,2:end);
T = data(2:end,1);
PutCall = 'C';
[R C] = size(calls);

% Input the spot price, risk free rate, and dividend yield.
S = 100;
r = 0.04;
q = 0;

% Maturity in years.
Basis = 365;
T = T./Basis;

% Parameters for the bisection algorithm.
% Tolerance, Maximum iterations, and lower and upper IV values.
Tol = 1e-6;
MaxIter = 1e5;
a = 1e-5;
b = 1e2;

% Extract the implied volatility using the Bisection algorithm.
for i=1:R
	for j=1:C
		IV(i,j) = BisecBSV(PutCall,S,K(j),r,q,T(i),a,b,calls(i,j),Tol,MaxIter);
	end
end
IV

% Construct total implied variance.
for j=1:C
	w(:,j) = IV(:,j).^2 .* T;
end

% Plot the total implied variance and show arbitrage at low strikes:
% the lines cross and they should not cross.
W = w';
plot(K,W)
xlabel('Strike Price')
ylabel('Total Implied Variance')
pause

% Remove the strikes where arbitrage is evident 
plot(K(5:end)', W(5:end,:))
xlabel('Strike Price')
ylabel('Total Implied Variance')
pause
clear w

% Resize the total variance matrix, w.
w = W(5:end,:)';
[R C] = size(w);

% Forwards prices using resized strike vector.
K = K(5:end);
y = log(K./S);

% Increments for finite-difference derivatives.
dt = 10e-10;
dy = 1e-10;

% First derivatives with respect to time.
% Calculate backward, forward, and central differences, but retain central.
for j=1:C
	for i=1:R
		T0 = T(i) - dt;
		T1 = T(i);
		T2 = T(i) + dt;
		w0 = interp1(T,w(:,j),T0,'cubic');
		w1 = w(i,j);
		w2 = interp1(T,w(:,j),T2,'cubic');
		Forward   = (w2-w1)/(T2-T1);
		Backward  = (w1-w0)/(T1-T0);
		Central   = (w2-w0)/(T2-T0);
		dwdT(i,j) = Central;
	end
end

% First derivatives with respect to moneyness.
for i=1:R
	for j=1:C
		y0 = y(j) - dy;
		y1 = y(j);
		y2 = y(j) + dy;
		w0 = interp1(y,w(i,:),y0,'cubic');
		w1 = w(i,j);
		w2 = interp1(y,w(i,:),y2,'cubic');
		Forward   = (w2-w1)/(y2-y1);
		Backward  = (w1-w0)/(y1-y0);
		Central   = (w2-w0)/(y2-y0);
		dwdy(i,j) = Central;
	end
end

% Second derivative with respect to moneyness.
for i=1:R
	for j=1:C
		y0 = y(j) - dy;
		y1 = y(j);
		y2 = y(j) + dy;
		w0 = interp1(y,dwdy(i,:),y0,'cubic');
		w1 = w(i,j);
		w2 = interp1(y,dwdy(i,:),y2,'cubic');
		Forward   = (w2-w1)/(y2-y1);
		Backward  = (w1-w0)/(y1-y0);
		Central   = (w2-w0)/(y2-y0);
		dw2dy2(i,j) = Central;
	end
end

% Mapping of implied vol to local vol using the formula in Jim Gatheral's
% book.
for i=1:R
	for j=1:C
		num = dwdT(i,j);
		den = 1 - y(j)/w(i,j)*dwdy(i,j) ...
			+ 0.25*(-0.25 - 1/w(i,j) + y(j)^2/w(i,j)^2)*(dwdy(i,j))^2 ...
		    + 0.5*dw2dy2(i,j);
		LV(i,j) = sqrt(num/den);
	end
end

% Display the vendor local vol and the calculated local vol.
LV

% Plot the local vols
surf(K,T,LV);
xlabel('Strike')
ylabel('Maturity in years')
zlabel('Local Volatility')
