function [y]=N(x)
  [y,Q]=cdfnor("PQ",x,0,1);
endfunction

function [x] = CALL(u,S_0,K,sigma,r,T)
  G   = cdfnor("X",zeros(u),ones(u),u,ones(u)-u);
  S_T = S_0 * exp((r - sigma^2/2)*T + sqrt(T) * sigma * G);
  x=exp(-r*T) * max(S_T-K,0);
endfunction

function [res] = BS_Call(S_0,K,sigma,r,T)
  d1=(log(S_0/K)+(r+sigma^2/2)*T)/(sigma*sqrt(T));
  d2=d1-sigma*sqrt(T);
  res=S_0*N(d1)-K*exp(-r*T)*N(d2);
endfunction

S_0=100;K=100;
r=0.05;sigma=0.3;T=1;
CALL_VDC=CALL(VDC,S_0,K,sigma,r,T);
CALL_RAND=CALL(Rand,S_0,K,sigma,r,T);
res=BS_Call(S_0,K,sigma,r,T);
for i=1:100
    conv_vdc(i) = abs(mean(CALL_VDC(1:i*100))-res);
    conv_ran(i)	= abs(mean(CALL_RAND(1:i*100))-res);
end
plot2d(conv_ran);
plot2d([conv_ran,conv_vdc]);