sum=cbnd(a,b,p)

function sum=cbnd(a,b,p)
%calculation of cumulative probability in bivariate normal distribution;
% using Z.Drezner's method;

%AA = [0.24840615, 0.39233107, 0.21141819, 0.03324666, 0.00082485334];
%BB= [0.10024215, 0.48281397, 1.0609498, 1.7797294, 2.6697604];
AA=[0.3253030 0.4211071 0.1334425 0.006374323];
BB=[0.1337764 0.6243247 1.3425378 2.2626645];
if a<=0 & b<=0 & p<=0
    inva=a/sqrt(2*(1-p^2));
    invb=b/sqrt(2*(1-p^2));
    sum=0;
    for i=1:4
        for j=1:4
            sum=sum+AA(i)*AA(j)*exp(inva*(2*BB(i)-inva)+invb*(2*BB(j)-invb)+2*p*(BB(i)-inva)*(BB(j)-invb));
        end
    end
    sum=sum*sqrt(1-p^2)/pi;
elseif a<=0 & b>=0 & p>=0
    sum=normcdf(a)-cbnd(a,-b,-p);
elseif a>=0 & b<=0 & p>=0
    sum=normcdf(b)-cbnd(-a,b,-p);
elseif a>=0 & b>=0 & p<=0
    sum=normcdf(a)+normcdf(b)-1+cbnd(-a,-b,p);
elseif a*b*p>0
    p1 = (p * a - b) * sgn(a) / sqrt(a ^ 2 - 2 * p * a * b + b^2);
    p2 = (p * b - a) * sgn(b) / sqrt(a ^ 2 - 2 * p * a * b + b^2);
    delta = (1 - sgn(a) * sgn(b)) / 4;
    sum = cbnd(a, 0, p1) + cbnd(b, 0, p2) - delta;
end

function ss=sgn(x)
if x>=0
    ss=1;
else
    ss=-1;
end
return