function y=BS_MC_CV_call(S0,K,T,r,sigma,M)
    stacksize('max')
    tic();
    X=rand(1,M,'normal');
    S=S0*exp(sigma*sqrt(T)*X+(r-sigma^2/2)*T);
    m=sum(S)/M;
    VarEstS=sum((S-m)^2)/(M-1);
    P=sum(max(S-K,0))/M;
    Lstar=sum((max(S-K,0)-P).*(S-m))/(M-1)/VarEstS;
    C=exp(-r*T)*(max(S-K,0)-Lstar*S);
    price=sum(C)/M+Lstar*S0;
    VarEst=sum((C+Lstar*S0-price)^2)/(M-1);
    RMSE=sqrt(VarEst)/sqrt(M);
    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