function xxx = lect03ivrnd
clear;
format compact;
pts = load('c:/VolSurf0.dat');
K = pts(:,1);
F = pts(:,2);
C = pts(:,3);
T = pts(:,4)/252;
IV = pts(:,5);
r = 0.03;
VM = mean(IV)
% fit volatility surface
MN = mness(K, F, T);
pars = [ VM, -1., 3., 0.];
opts = optimset('Display','iter','TolFun',0.00001);
x = lsqnonlin(@iverrors,pars,[], [], opts, MN, T, IV);
x
figure(1);
plot3(MN,T,IV,'o'); hold on;
[DM,TM] = meshgrid(-.2:.02:.6,0:.02:1.2);
ZM = ivfitted(x, DM, TM);
mesh(DM, TM, ZM); hold off;
% fit spline to the implied volatilities
T0 = 50/252; % maturity
F0 = exp(r*T0);
% create strikes
K0 = 0.8:.001:1.2;
h  = 1e-6; % the derivative spacing
KH = K0 + h;
KL = K0 - h;
% create moneyness
MN0 = mness(K0, F0, T0);
MNH = mness(KH, F0, T0);
MNL = mness(KL, F0, T0);
% create fitted volatility curve
V0 = ivfitted(x, MN0, T0); % 
VH = ivfitted(x, MNH, T0);
VL = ivfitted(x, MNL, T0);
% create call prices
C0 = exp(r*T0) * bs(K0, F0, V0, r, T0);
CH = exp(r*T0) * bs(KH, F0, VH, r, T0);
CL = exp(r*T0) * bs(KL, F0, VL, r, T0);
% Breeden and Litzenberg approximation
PDF = (CH + CL - 2 * C0) / (h^2);
% Log-normal for ATM vol
LN = lognpdf(K0, r*T0, ivfitted(x, 0, T0)*sqrt(T0));
figure(3);
plot(K0, LN,'g');
hold on;
plot(K0, PDF,'b');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y = ivfitted(pars, MN, T)
y = pars(1) + pars(4) * T + pars(3) * exp( pars(2) * T) .* MN;

function y = iverrors(pars, MN, T, IV)
y = ivfitted(pars, MN, T) - IV;

function y = bs(K, F, sig, r, T)
% The simple Black-Scholes European call option Price
d1 = ( log(F./K)+.5*sig.*sig.*T ) ./ sig ./ sqrt(T);
d2 = d1 - sig.*sqrt(T);
n1 = normcdf(d1);
n2 = normcdf(d2);
y = ( F.*n1 - K.*n2 ) .* exp(-r.*T);

function y = bsdelta(K, F, sig, T)
% The simple Black-Scholes European call option Price
d1 = ( log(F./K)+.5*sig.*sig.*T ) ./ sig ./ sqrt(T);
y = normcdf(d1);

function y = bsidelta(D, F, sig, T)
d1 = norminv(D);
y = F .* exp( - d1.*sig.*sqrt(T) + .5*sig.*sig.*T );

function y = mness(K, F, T)
% The moneyness function
y = log(F./K) ./ sqrt(T);

function y = imness(MN, F, T)
% The inverse moneyness function
y = F .* exp( - MN .* sqrt(T) );