clear all; clf;
global Nx h x r dt S0 t T T2 NK_call xx_call yy_call Xq Yq A_call
% KOSPI200 Call Option Price on 30 December 2021
market_call=[6.80 4.03 2.13 1.03 0.47; 10.25 7.50 5.13 3.45 2.20
12.80 10.50 8.05 6.65 4.15];
T(1)=datenum(2022,01,13)-datenum(2021,12,30);
T(2)=datenum(2022,02,10)-datenum(2021,12,30);
T(3)=datenum(2022,03,10)-datenum(2021,12,30);
LT=length(T); dt=1/365; t=linspace(0,dt*T(LT),T(LT)+1); S0=394.19; r=0.02;
Lx=3*S0; Nx=round(Lx); x=linspace(0,Lx,Nx); h=x(2)-x(1);
[Xq,Yq] = meshgrid(x,t); T2(1)=0;
for i=2:LT-1
T2(i)=(T(i-1)+T(i))/2;
end
T2(LT)=T(LT); ST_call=[390:5:410]; NK_call=length(ST_call);
for k=1:LT
for i=1:NK_call
A_call(NK_call*(k-1)+i,1)=ST_call(i);
A_call(NK_call*(k-1)+i,2)=T(k);
exv_call(NK_call*(k-1)+i)=market_call(k,i);
end
end
xdata_call=A_call(:,1); CV_call=exv_call'; % market_call Price
vol0_call=ones(8,1)*0.3; lb_call=0*vol0_call; ub_call=0*vol0_call+1;
vol0_call(end-2:end)=1; ub_call(end-2:end)=3;
xx_call=[0 0 0 Lx Lx Lx vol0_call(end-2)*S0 vol0_call(end-1)*S0 ...
vol0_call(end)*S0]'; yy_call=[dt*T2 dt*T2 dt*T2]';

figure(1); clf; plot(xx_call,yy_call,'o'); grid on; hold on
w_call=[vol0_call(1)*ones(3,1); vol0_call(2)*ones(3,1); vol0_call(3:5)];
F_call = scatteredInterpolant(xx_call,yy_call,w_call);
for i=1:size(Xq,2)
for j=1:size(Xq,1)
vf_call(i,j)=F_call(x(i),t(j));
end
end
max_it=5; options = optimset('MaxIter',max_it);
vol_call = lsqcurvefit(@BScallT, vol0_call, xdata_call, CV_call, lb_call, ...
ub_call, options);
w_call=[vol_call(1)*ones(3,1); vol_call(2)*ones(3,1); vol_call(3:5)];
xx_call=[0 0 0 x(end) x(end) x(end) vol_call(end-2)*S0 ...
vol_call(end-1)*S0 vol_call(end)*S0]';
F_call = scatteredInterpolant(xx_call,yy_call,w_call);
for i=1:size(Xq,2)
for j=1:size(Xq,1)
vf_call(i,j)=F_call(x(i),t(j));
end
end
plot(xx_call,yy_call,'b*'); grid on; hold on
xx_call1=xx_call; fs = 17;
figure(3); clf; hold on; grid on; set(gca,'fontsize',15)
mesh(x,t,vf_call'); axis([0 1200 0 0.2 0.1 0.25]); view([-5 45])
text('Interpreter','latex','String','$S$','position',[1080 -0.03 0.1],'fontsize',fs)
text('Interpreter','latex','String','$t$','position',[-108 0.0785 0.168],'fontsize',fs)
lw=1; fs=18; ms=10; ph=3; height=15;
figure(7); clf; hold on; box on; grid on
V_call=BScallT(vol_call, xdata_call);
plot(A_call(:,1),CV_call,'ko','linewidth',lw,'markersize',ms);
plot(A_call(:,1),V_call,'k+','linewidth',lw,'markersize',ms);
axis([ST_call(1)-.5 ST_call(end)+.5 -0.05 height]); set(gca,'fontsize',fs)
xticks(ST_call)
legend('Interpreter','latex','string',{'Market price','Numerical price', ...
'Numerical Call Put'},'location','northeast')
text('Interpreter','latex','string','$K$','fontsize',fs, ...
'position',[407.7 -1.11265 0])
text('Interpreter','latex','string','Price','fontsize',fs, ...
'position',[386.5 13.77 0])

function V=BScallT(vol_call, xdata_call)
global Nx h x t r dt S0 xx_call yy_call A_call;

w_call=[vol_call(1)*ones(3,1); vol_call(2)*ones(3,1); vol_call(3:5)];
xx_call=[0 0 0 x(end) x(end) x(end) vol_call(end-2)*S0 ...
vol_call(end-1)*S0 vol_call(end)*S0]';
F_call = scatteredInterpolant(xx_call,yy_call,w_call);
for i=1:Nx
for j=1:length(t)
vf_call(i,j)=F_call(x(i),t(j));
end
end
for j=1:length(A_call(:,1))
Nt=A_call(j,2); u(:,1)=max(x-A_call(j,1),0);
for n=1:Nt
for i=2:Nx
a(i-1)=r*x(i)/(2*h)-(vf_call(i,Nt-n+1)*x(i))^2/(2*h^2);
d(i-1)=1/dt+(vf_call(i,Nt-n+1)*x(i))^2/h^2+r;
c(i-1)=-r*x(i)/(2*h)-(vf_call(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