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

Print["ECall= ",EuropeanCall[58.5,60,0.29,0.04,0.02,0.3]];
Print["EPut= ",EuropeanPut[58.5,60,0.29,0.04,0.02,0.3]];

SOR[U_,b_,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;
];Print["error=",error];
Print["ni=",ni]]

EuropeanThetaSchemeSOR[u_,h_,k_,n_,m_,theta_,omega_,eps_]:=Module[
{tau,alfa}, tau=0.0;alfa=k/(h*h);
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)*(u[i+1,j-1]-2*u[i,j-1]+u[i-1,j-1]);
U[i]=u[i,j-1]];
    g[-n]=Lbc[x[-n],tau];
    U[-n]=g[-n];
    g[n]=Rbc[x[n],tau];
    U[n]=g[n];
    b[-n+1]+=alfa*theta*g[-n];
    b[n-1]+=alfa*theta*g[n];
    SOR[U,b,n,alfa,theta,omega,eps];
    For[i=-n,i<=n,i++,
        u[i,j]=U[i];
        ]
    ]
]

(* Implicit scheme *)

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]
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;
m=5;
h=0.1;
k=T*0.5*sigma^2/m;
theta=1.0;
omega=1;
eps=0.000001;
EuropeanThetaSchemeSOR[u,h,k,n,m,theta,omega,eps]

g1=Plot[Payoff[x,m*k],{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]

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]

Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i], X, sigma, r, d,T])^2,{i,-n,n}]/(2 n +1)]

n=20;
m=10;
h=0.05;
k=T*0.5*sigma^2/m;
theta=1.0;
omega=1;
eps=0.000001;
EuropeanThetaSchemeSOR[u,h,k,n,m,theta,omega,eps]

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

n=20;
m=20;
h=0.05;
k=T*0.5*sigma^2/m;
theta=1.0;
omega=1;
eps=0.000001;
EuropeanThetaSchemeSOR[u,h,k,n,m,theta,omega,eps]

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 scheme *)

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]
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;
m=5;
h=0.1;
k=T*0.5*sigma^2/m;
theta=0.5;
omega=1.0;
eps=0.000001;
EuropeanImplicit[u,h,k,n,m,theta,omega,eps]

g1=Plot[Payoff[x,m*k],{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]

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]

Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i], X, sigma, r, d,T])^2,{i,-n,n}]/(2 n +1)]

Sqrt[Sum[(Vacd[i,m]-EuropeanCall[Ssur[i], X, sigma, r, d,T])^2,{i,-n,n}]/(2 n +1)]

n=20;
m=10;
h=0.05;
k=T*0.5*sigma^2/m;
theta=0.5;
omega=1;
eps=0.000001;
EuropeanImplicit[u,h,k,n,m,theta,omega,eps]

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