clear all;
global Nx h x r dt S0 t T T2 NK xx yy Xq Yq A;
S0 = 100; r = 0.01; Nx = 301; Lx = 3*S0; x = linspace(0,Lx,Nx); h = x(2)-x(1);
dt = 1.0/360.0; T = [90 180 270 360]; LT = length(T);
t = linspace(0,dt*T(LT),T(LT)+1); [Xq,Yq] = meshgrid(x,t); T2(1) = 0;
for ti = 2:LT-1
T2(ti) = (T(ti-1)+T(ti))/2;
end
T2(LT) = T(LT);
for vj = 1:T(LT)+1
for vi = 1:Nx
exvol(vi,vj) = 0.00001*(x(vi)-S0)^2+0.1*cos(pi*t(vj))-0.2*t(vj)+0.4;
xxx(Nx*(vj-1)+vi) = x(vi); yyy(Nx*(vj-1)+vi)=t(vj);
exv(Nx*(vj-1)+vi) = exvol(vi,vj);
end
end
ST = [95:2.5:105]; NK = length(ST);
for k = 1:LT
A(1+NK*(k-1):k*NK,1) = ST; A(1+NK*(k-1):k*NK,2) = T(k);
end
xdata = A(:,1); xx = xxx; yy = yyy; CV = BScallT(exv,xdata);
tmp = [x(1) 0.5*S0 S0 1.5*S0 x(end)]; xx = []; yy = [];
for j = 1:LT
xx = [xx tmp]; yy(1+length(tmp)*(j-1):j*length(tmp)) = dt*T2(j);
end
vol0 = 0.5*ones(length(xx),1); lb = 0*vol0; ub = 0*vol0+1.0;
options = optimset('MaxIter',5);
vol = lsqcurvefit(@BScallT, vol0, xdata, CV, lb, ub, options);
U = griddata(xx,yy,vol,Xq,Yq); V = BScallT(vol, CV);
figure(1); clf;
mesh([x(1:10:end) x(end)],[t(1:10:end) t(end)],exvol([1:10:end end],[1:10:end end])')
axis([x(1) x(end) t(1) t(end) 0 1]); view(-15,45); grid on; box on
set(gca,'fontsize',18); title('Given local volatility surface')
xlabel('$S$','Interpreter','latex'); ylabel('$t$','Interpreter','latex');
zlabel('$\sigma$','Interpreter','latex','rotation',0)
figure(2); clf; hold on; box on; grid on
plot(A(:,1),CV,'ko','linewidth',1,'markersize',10);
plot(A(:,1),V,'kx','linewidth',1,'markersize',10);
axis([ST(1)-.5 ST(end)+.5 min([min(CV) min(V)])-.5 max([max(CV) max(V)]+.5)]);
set(gca,'fontsize',18); legend('Interpreter','latex','string',...
{'Given $\sigma(S,t)$','Reconstructed $\sigma(S,t)$'})
xlabel('$K$','Interpreter','latex'); ylabel('Price','Interpreter','latex');
figure(3); clf;
mesh([x(1:10:end) x(end)],[t(1:10:end) t(end)],U([1:10:end end],[1:10:end end]))
axis([x(1) x(end) t(1) t(end) 0 1]); view(-15,45); grid on; box on
set(gca,'fontsize',18); title('Reconstructed local volatility surface')
xlabel('$S$','Interpreter','latex'); ylabel('$t$','Interpreter','latex');
zlabel('$\sigma$','Interpreter','latex','rotation',0)
figure(4); clf; hold on; view(-15,45); grid on; box on
mesh([x(1:10:end) x(end)],[t(1:10:end) t(end)],U([1:10:end end],[1:10:end end]))
mesh([x(1:10:end) x(end)],[t(1:10:end) t(end)],exvol([1:10:end end],[1:10:end end])')
axis([x(1) x(end) t(1) t(end) 0 1]); set(gca,'fontsize',18)
title('Overlapped local volatility surfaces')
xlabel('$S$','Interpreter','latex'); ylabel('$t$','Interpreter','latex');
zlabel('$\sigma$','Interpreter','latex','rotation',0)
function V = BScallT(vol, CV)
global Nx h x r dt S0 xx yy Xq Yq A;
vf = griddata(xx,yy,vol,Xq,Yq)';
for j = 1:length(A(:,1))
Nt = A(j,2); u(:,1) = max(x-A(j,1),0);
for n = 1:Nt
for i = 2:Nx
a(i-1) = r*x(i)/(2*h)-(vf(i,Nt-n+1)*x(i))^2/(2*h^2);
d(i-1) = 1/dt+(vf(i,Nt-n+1)*x(i))^2/h^2+r;
c(i-1) = -r*x(i)/(2*h)-(vf(i,Nt-n+1)*x(i))^2/(2*h^2);
end
a(Nx-1) = a(Nx-1)-c(Nx-1); d(Nx-1) = d(Nx-1)+2*c(Nx-1);
b = u(2:Nx,n)/dt; u(2:Nx,n+1) = thomas(a,d,c,b);
end
V(j,1) = interp1(x,u(:, Nt+1),S0);
end
end
function x = thomas(alpha,beta,gamma,f)
n = length(f);
for i = 2:n
mult = alpha(i)/beta(i-1); beta(i) = beta(i)-mult*gamma(i-1); f(i) = f(i)-mult*f(i-1);
end
x(n) = f(n)/beta(n);
for i = n-1:-1:1
x(i) = (f(i)-gamma(i)*x(i+1))/beta(i);
end
end