function y=BS_MC_AV_call(S0,K,T,r,sigma,M)
    stacksize('max')
    tic();
    Y=rand(1,M/2,'normal');
    S=S0*exp(sigma*sqrt(T)*Y+(r-sigma^2/2)*T);
    S2=S0*exp(-sigma*sqrt(T)*Y+(r-sigma^2/2)*T);
    C=exp(-r*T)*max(S-K,0);
    C2=exp(-r*T)*max(S2-K,0);
    C=0.5*(C+C2);
    price=sum(C)/(M/2);
    VarEst=sum((C-price)^2)/(M/2-1);
    RMSE=sqrt(VarEst)/sqrt(M/2);
    CI95=[price-1.96*RMSE,price+1.96*RMSE];
    CI99=[price-2.58*RMSE,price+2.58*RMSE];
    time=toc();
    y=[price time RMSE; CI95 0; CI99 0]
endfunction