Numericke riesenie Black - Scholesovej rovnice

NormalCDF[u_]:=CDF[NormalDistribution[0,1],u];(*v orednaske N(u)*)
Print[NormalCDF[u]]
Plot[NormalCDF[u],{u,-4,4}]

EuropeanCall[S_,X_,sigma_,r_,d_,tau_]:=
Module[{d1,d2},
d1=(Log[S/X]+(r-d+sigma^2/2)*tau)/(sigma*Sqrt[tau]);
d2=d1- sigma*Sqrt[tau];
Return[Exp[-d*tau]*S*NormalCDF[d1]-X*Exp[-r*tau]*NormalCDF[d2]]
];
(*Explicitny vzorec pre Europsku Put opciu so spojitou dividendou*)
EuropeanPut[S_,X_,sigma_,r_,d_,tau_]:=
EuropeanCall[S,X,sigma,r,d,tau]+X*Exp[-r tau] - S*Exp[-d tau];

X=60;
S=58.5;
tau=0.3;
sigma=0.29;
r=0.04;
d=0.02;
Print["ECall=", EuropeanCall[S,X,sigma,r,d,tau]]
Print["EPut=", EuropeanPut[S,X,sigma,r,d,tau]]

n=5;
Cm=SparseArray[{Band[{1,2}]->-a,Band[{1,1}]->1+2a,Band[{2,1}]->-a},2 n-1];
MatrixForm[Normal[Cm]]

EuropeanThetaScheme[u_,h_,k_,n_,m_,theta_]:=
Module[{tau, alfa, a, Cm, bm},
tau=0.0;
alfa=k/(h*h);
a=alfa*theta;
Cm=SparseArray[{Band[{1,2}]->-a,Band[{1,1}]->1+2a,Band[{2,1}]->-a},2 n-1];
For[i=-n,i≤n,i++,x[i]=i*h;u[i,0]=Payoff[x[i],tau]];
For[j=1,j≤m,j++,
(*Print["casovy krok cislo ",j];*)
tau=j*k;
For[i=-n+1,i<n,i++,
b[i]=u[i,j-1]+alfa*(1.0-theta)(*prava strana*)
(u[i+1,j-1]-2*u[i,j-1]+
u[i-1,j-1])];
u[-n,j]=g[-n]=Lbc[x[-n],tau];(*okrajove podmienky*)
u[n,j]=g[n]=Rbc[x[n],tau];
b[-n+1]+=a*g[-n];(*uprava pravej strany*)
b[n-1]+=a*g[n];
bm=Table[b[i],{i,-n+1,n-1}];(*tvorba pravej straby*)
U=LinearSolve[Cm,bm];
For[i=-n+1,i≤n-1,i++,
u[i,j]=U[[i+n]]]
]
]

(*Implicitna schema*)

X=60.0;
r=0.1;
sigma=0.29;
d=0.08;
T=1.0;
k1=r/(0.5*sigma^2);
k2=(r-d)/(0.5*sigma^2);
Payoff[x_,tau_]:=Exp[0.25((k2-1)^2+4*k1)*tau]*Max[Exp[0.5(k2+1)*x]-Exp[0.5(k2-1)*x],0] (*obsahuje aj tau, je to vseobecnejsie a aj pre ameriocku opciu*)
Lbc[x_,tau_]:=0;
Rbc[x_,tau_]:=Exp[0.5(k2+1)*x+0.25((k2+1)^2)*tau]-
Exp[0.5(k2-1)*x+0.25((k2-1)^2)*tau]

n=10;(*10 bodov dolava, a doprava - teda 21*)
m=5; (*pocet casovych krokov*)
h=0.1;
k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
theta=1.0;
EuropeanThetaScheme[u,h,k,n,m,theta]
g1=Plot[Payoff[x,0],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]];
g2=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];

g3=ListPlot[Table[{Ssur[i],Vacd[i,0]},{i,-n,n}],PlotStyle->RGBColor[1,0,0],Joined->True];
g4=ListPlot[Table[{Ssur[i],Vacd[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g3,g4]

g3=ListPlot[Table[{Ssur[i],Vacd[i,0]},{i,-n,n}],PlotStyle->RGBColor[1,0,0],Joined->True];
g4=ListPlot[Table[{Ssur[i],Vacd[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g3,g4]
g5=Plot[EuropeanCall[S,X,sigma,r,d,T],{S,Ssur[-n],Ssur[n]},PlotStyle->RGBColor[0,1,0]];
Show[g3,g4,g5]

g6=ListPlot[Table[{Ssur[i],Vacd[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1]];
Show[g3,g6,g5]

h=0.1;n=10;m=5;
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.05;n=20;m=10;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
theta=1.0;
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.025;n=40;m=20;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
theta=1.0;
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.0125;n=80;m=40;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
theta=1.0;
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.00625;n=160;m=80;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
theta=1.0;
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

(*vyskusat implicitnu, explicitnu a crank nicolsonivu schemu s roznymi casovymi a priestorovymu krokmi, vytvorit tabiulky, ktore dokumentuju rad konvergenicie metody*)

(*explicitna*)

n=10;(*10 bodov dolava, a doprava - teda 21*)
m=10; (*pocet casovych krokov*)
h=0.1;
k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
theta=0.;
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
g3=ListPlot[Table[{Ssur[i],Vacd[i,0]},{i,-n,n}],PlotStyle->RGBColor[1,0,0],Joined->True];
g4=ListPlot[Table[{Ssur[i],Vacd[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
g5=Plot[EuropeanCall[S,X,sigma,r,d,T],{S,Ssur[-n],Ssur[n]},PlotStyle->RGBColor[0,1,0]];
g6=ListPlot[Table[{Ssur[i],Vacd[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1]];
Show[g4,g3,g6,g5]
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

1-2k/(h*h)

theta=0;
h=0.1;n=10;m=10;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.05;n=20;m=40;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.025;n=40;m=160;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.025;n=40;m=20;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.0125;n=80;m=40;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.00625;n=160;m=80;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

(*crank-nicholson*)

n=10;(*10 bodov dolava, a doprava - teda 21*)
m=5; (*pocet casovych krokov*)
h=0.1;
k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
theta=0.5;
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
g3=ListPlot[Table[{Ssur[i],Vacd[i,0]},{i,-n,n}],PlotStyle->RGBColor[1,0,0],Joined->True];
g4=ListPlot[Table[{Ssur[i],Vacd[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
g5=Plot[EuropeanCall[S,X,sigma,r,d,T],{S,Ssur[-n],Ssur[n]},PlotStyle->RGBColor[0,1,0]];
g6=ListPlot[Table[{Ssur[i],Vacd[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1]];
Show[g4,g3,g6,g5]
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

theta=0.5;
h=0.1;n=10;m=5;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.05;n=20;m=10;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.025;n=40;m=20;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.0125;n=80;m=40;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

h=0.00625;n=160;m=80;k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

For[
h=0.1;n=10;m=5;ii=0,ii≤5,ii++,
k=T*0.5*sigma^2/m;(*dlzka cas. kroku*);
theta=1.0;
EuropeanThetaScheme[u,h,k,n,m,theta];
Table[
Table[Vacd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
implicit[ii]=Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)];

]

(*prerobit payoff, rbc, lbc pre Put opciu*)

X=60.0;
r=0.1;
sigma=0.29;
d=0.08;
T=1.0;
k1=r/(0.5*sigma^2);
k2=(r-d)/(0.5*sigma^2);
Payoff[x_,tau_]:=Exp[0.25((k2-1)^2+4*k1)*tau]*Max[0,Exp[0.5(k2-1)*x]-Exp[0.5(k2+1)*x]] (*obsahuje aj tau, je to vseobecnejsie a aj pre ameriocku opciu*)
Lbc[x_,tau_]:=-Exp[0.5(k2+1)*x+0.25((k2+1)^2)*tau]+Exp[0.5(k2-1)*x+0.25((k2-1)^2)*tau];
Rbc[x_,tau_]:=0;

n=10;(*10 bodov dolava, a doprava - teda 21*)
m=5; (*pocet casovych krokov*)
h=0.1;
k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
theta=1.0;
EuropeanThetaScheme[u,h,k,n,m,theta]
g1=Plot[Payoff[x,0],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]];
g2=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

h=0.1;n=10;m=5;
Table[
Table[Vapd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vapd[i,m]-EuropeanPut[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

g3=ListPlot[Table[{Ssur[i],Vapd[i,0]},{i,-n,n}],PlotStyle->RGBColor[1,0,0],Joined->True];
g4=ListPlot[Table[{Ssur[i],Vapd[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g3,g4]
g5=Plot[EuropeanPut[S,X,sigma,r,d,T],{S,Ssur[-n],Ssur[n]},PlotStyle->RGBColor[0,1,0]];
Show[g3,g4,g5]
g6=ListPlot[Table[{Ssur[i],Vapd[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1]];
Show[g3,g6,g5]

n=20;(*10 bodov dolava, a doprava - teda 21*)
m=10; (*pocet casovych krokov*)
h=0.05;
k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
theta=0.5;
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vapd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Sqrt[Sum[(Vapd[i,m]-EuropeanPut[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]

n=10;(*10 bodov dolava, a doprava - teda 21*)
m=20; (*pocet casovych krokov*)
h=0.1;
For[ii=1,ii≤5,ii++,
k=T*0.5*sigma^2/m;(*dlzka cas. kroku*)
theta=0.5;
EuropeanThetaScheme[u,h,k,n,m,theta]
Table[
Table[Vapd[i,j]=X*Exp[-0.5 (k2-1)*i*h-(0.25(k2-1)^2+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[i*h],{i,-n,n}];
Print[Sqrt[Sum[(Vapd[i,m]-EuropeanPut[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)]];
Implicit[ii]=Sqrt[Sum[(Vapd[i,m]-EuropeanPut[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2 n+1)];
n=n+ii*10;
m=m+ii*5;
h=h/2;
]

Table[Implicit[i],{i,1,5}]

(*pouzitie SOR algoritmu namiesto linearneho solvera*)