function xxx = lect02ivsurfcal
clear;
format compact;
% initial values
par = [0.0100   10.0000    1.0000   -0.6000   0.0100]; % init
parL= [0.0010    0.5000    0.0100   -0.9000   0.0010]; % lower bounds
parU= [Inf      Inf        Inf       0.9000   Inf   ]; % upper bounds
% data load as on 17-Sept-93
pts = load('c:/VolSurf0.dat');
K = pts(:,1);
F = pts(:,2);
C = pts(:,3);
T = pts(:,4);
IV = pts(:,5);
r  = 0.0303; % interest rate
% calibration
par
x = lsqnonlin(@sumsq,par,parL,parU,optimset('Display','iter','TolFun',0.00001),r,K,F,C,T,IV);
x
Vinit  = x(5);
Vbar   = x(1);
kappav = x(2);
thetav = Vbar*kappav;
sigmav = x(3);
rho    = x(4);
% presentation of results and plots
% create grid and theoretical IVs
[KK,TT] = meshgrid(0.8:0.005:1.2,0.05:0.05:0.80);
CC = [];
for indx=1:size(TT(:,1))
    t = TT(indx,1); % Take times
    GRD = SVprices(Vinit,r,t,Vbar,kappav,thetav,sigmav,rho);
    GK = GRD(:,1);  % Strikes
    GC = GRD(:,2);  % Option Prices
    FC = interp1(GK,GC,KK(1,:));
    IVs = impvol(KK(1,:),FC,r, t);
    CC = [CC; IVs'];
end
mesh(KK,TT,CC); hold;
plot3(K./F.*exp(r*T/365),T/365,IV,'ro');hold off;

%%%%% SUM OF SQUARED DIFFERENCE (THEORETICAL IVs - ACTUAL IVs)^2 %%%%%
function SS = sumsq(par,r,K,F,C,T,IV)
Vinit  = par(5);
Vbar   = par(1);
kappav = par(2);
thetav = Vbar*kappav;
sigmav = par(3);
rho    = par(4);
m = size(K,1); % m: rows
TI = unique(T);
SS = []; % errors
for jndx=1:size(TI)
    T0 = TI(jndx);
    t  = T0/365;
    KI = [];
    CI = [];
    for indx=1:m
        % take correct maturity times
        if (T(indx)==T0)
            KI = [KI; K(indx)/F(indx)*exp(r*t)]; % moneyness
            CI = [CI; C(indx)/F(indx)];          % option price
        end      
    end
    GRD = SVprices(Vinit,r,t,Vbar,kappav,thetav,sigmav,rho);
    GK = GRD(:,1);
    GC = GRD(:,2);
    GGK = [];
    GGC = [];
    for indx=1:size(GK)
        if (GK(indx) > 0.7 & GK(indx)<1.3)
            GGK = [GGK;GK(indx)];
            GGC = [GGC;GC(indx)];
        end
    end
    IVs = impvol(GGK',GGC',r, t);
    FV = interp1(GGK,IVs,KI);
    
    SS = [SS; FV];
end
SS = SS - IV;

%%%%% FUNCTIONS FOR CHAR. FUNCTION STOCHASTIC VOL. OPTION PRICING %%%%%
function y = SVprices(V0,r,t,Vbar,kappav,thetav,sigmav,rho)
% parameters
a = 600;            % endpoint of char fun grid (0,+a)
N =  4096;          % number of grid points
eta = a/N;          % char fun grid size
b = pi/eta;         % +/- bounds for the log-strike grid (-b,+b)
lambda = 2*pi/a;    % log-strike grid size
% dumping parameter as in Carr-Madan
aa = 1.25;
% create grids
u = (0:N-1) * eta;
x = -b + (0:N-1) * lambda;
% create char fun and invert
h = cfn(u, V0, r, t, kappav, thetav, sigmav, rho, aa);
h2 = exp(i*b*u) .* h * eta;
g   = fft(h2);
g2  = real( g .* exp(-aa*x) / pi);  % prices based on char fun
y = [exp(x)' g2'];

function y = cfn(th, V0, r, t, kappav, thetav, sigmav, rho, aa)
% The characteristic function of the SV model
% Adjusted to incorporate the Simpson weighting scheme
th1 = th - (aa+1)*i;
f1 = exp(-r*t) * cf0(V0, r, t, kappav, thetav, sigmav, rho, th1);
f2 = aa^2 + aa - (th.^2) + i*(2*aa+1)*th;
y  = f1 ./ f2;
% Create Simpson's weights
N  = size(th,2);
q1 = (-1).^(1:N);
q2 = eye(1,N);
S  = ( 3 + q1 - q2 )/3;
y  = y .* S;

function y = cf0(V0, r, t, kappav, thetav, sigmav, rho, u)
% The unadjusted char function
ksi  = sqrt( (kappav - (1 + i*u)*rho*sigmav).^2 - i*u.*(i*u + 1)*sigmav^2 );
ksi0 = sqrt( (kappav - i*u*rho*sigmav).^2 - i*u.*(i*u - 1)*sigmav^2 ); 
f21  = -2*thetav/sigmav^2* log(1 - ((ksi0 - kappav + i*u*rho*sigmav).*(1 - exp(-ksi0*t)))./(2*ksi0)); 
f22  = -(thetav/sigmav^2)*(ksi0 - kappav + i*u*rho*sigmav)*t; 
f24  = (i*u.*(i*u - 1).*(1 - exp(-ksi0*t)))./(2*ksi0 - (ksi0 - kappav + i*u*rho*sigmav).*(1 - exp(-ksi0*t)))*V0; 
  y  = exp(r*t*i*u + f21 + f22 + f24); 
   
%%%%% FUNCTIONS FOR IMPLIED VOLATILITIES %%%%%
function y = bs(x, r, sig, t)
% The simple Black-Scholes European call option Price
d1 = ( -log(x)+( r + .5*sig^2 )*t ) / sig/sqrt(t);
d2 = d1 - sig*sqrt(t);
n1 = normcdf(d1);
n2 = normcdf(d2);
y = n1 - x.*n2*exp(-r*t);

function y = dbs(sig, x, r, t, c)
y = bs(x, r, sig, t) - c;

function y = impvol(x, c, r, t)
y = [];
for indx=1:size(x,2)
    xx = x(1,indx);
    cc = c(1,indx);
    yy = fzero(@dbs,0.1,optimset('fzero'),xx,r,t,cc);
    if (yy < 0.01)
        yy = NaN;
    end
    y  = [y ; yy];
end