function y= Corr_AV_MC_call(S0,K,T,r,sigma,M)
    X=rand(1,M,'normal');
    S1=S0*exp(sigma*X*sqrt(T)+(r-1/2*sigma^2)*T);
    S2=S0*exp(sigma*(-X)*sqrt(T)+(r-1/2*sigma^2)*T);
    C1=exp(-r*T)*max(S1-K,0);
    C2=exp(-r*T)*max(S2-K,0);
    P1=1/M*sum(C1);
    P2=1/M*sum(C2);
    VarEst1=1/(M-1)*sum((C1-P1)^2);
    VarEst2=1/(M-1)*sum((C2-P2)^2);
    CoVarEst=1/(M-1)*sum((C1-P1).*(C2-P2));
    CorrEst=CoVarEst/sqrt(VarEst1*VarEst2);
    y=CorrEst
endfunction
