function y=BlackScholes_MC_CV_call(S0,K,T,r,sigma,M)
stacksize('max')
    tic();
    X=rand(1,M,'normal')
    S=S0*exp(sigma*X*sqrt(T)+(r-1/2*sigma^2)*T);
    C=max(S-K,0);
    m=1/M*sum(S);
    P=1/M*sum(C);
    VarEst=1/(M-1)*sum((S-m)^2);
    CoVarEst=1/(M-1)*sum((S-m).*(C-P));
    l=CoVarEst/VarEst;
    PCV=exp(-r*T)*(1/M)*sum(C-l*S)+l*S0;
    VarEstCV=1/(M-1)*sum((PCV-exp(-r*T)*(C-l*S)+l*S0)^2);
    CI95=[PCV-1.96*sqrt(VarEstCV)/sqrt(M),PCV+1.96*sqrt(VarEstCV)/sqrt(M)];
    CI99=[PCV-2.58*sqrt(VarEstCV)/sqrt(M),PCV+2.58*sqrt(VarEstCV)/sqrt(M)];
    time=toc();
    y=[PCV 0; CI95; CI99; time 0]
endfunction