NormalCDF[u_]:=CDF[NormalDistribution[0,1],u];(*v prednaske 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]]
];

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]]

SOR[U_,h_,n_,alfa_,theta_,omega_,eps_]:=Module[{error,y,ni,a},
{
ni=0;
a=alfa*theta;
error=1.0;
While[error>eps*eps,
error=0.0;
ni++;
i=-n+1;
y=(b[i]+a*U[i+1])/(1.0+2.0*a);
y=U[i]+omega*(y-U[i]);
error+=(U[i]-y)*(U[i]-y);
U[i]=y;
For[i=-n+2,i<n-1,i++,
y=(b[i]+a*(U[i-1]+U[i+1]))/(1.0+2.0*a)];
y=U[i]+omega*(y-U[i]);
error+=(U[i]-y)*(U[i]-y);
U[i]=y];
i=n-1;
y=(b[i]+a*U[i-1])/(1.0+2.0*a);
y=U[i]+omega*(y-U[i]);
error+=(U[i]-y)*(U[i]-y);
U[i]=y;
}
]

EuropeanThetaSchemeSOR[u_,h_,k_,n_,m_,theta_,omega_,eps_]:=
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[i]=u[i,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 - pridane dirichlet podm na lavo a napravo*)
b[n-1]+=a*g[n];
bm=Table[b[i],{i,-n+1,n-1}];(*tvorba pravej straby*)
U=SOR[U,b,n,alfa,theta,omega,eps];
For[i=-n+1,i≤n-1,i++,
u[i,j]=U[[i+n]]]
]
]

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] (*urcuje o aku opciu ide
*)(*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; (*dokopy 21 bodov - 10 nalavo 10 napravo od 0ly*)
m=20;  (*casovy krok*)
h=0.1;
k=T*0.5*sigma^2/m;
theta=1.0;   (*urcuje o aky typ schemy ide - implicit explicit , cranc-nicosonova...*)
EuropeanThetaSchemeSOR[u,h,k,n,m,theta,1.6,0.0001]
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]
(*spetna transformacia*)
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[g3,g6,g5]
Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]