/*           Bivariate Normal (vanderEnde,Jun96)           */

/* THIS IS PROGRAM EndeProb.EDT:
Calculates cdfbvn with increased accuracy (with more approximating
coeffcients). Replaces my old version of bvn_pi4 & bvn_mask
with the new bvn_pi2 & bvn_mask. --Martin van der Ende
*/
/*
Note from Ron Schoenberg:
If you are you computing the logarithm of the bivarte normal cdf, it is
difficult to get the accuracy you require in the tails.  The Win95 and OS/2
versions of GAUSS, however, have new functions for the log of the normal
cdf's, lncdfn(), lncdfn2(), lnbvncdf(), and lnmvncdf(), which generate
16 places of accuracy in the tails using an algorithm developed by
Matthew Ford.  It adds a little computing time, but the added accuracy
should ultimately reduce time to convergence.
*/



/* *************************************************************
   *                                                           *
   * The procedure bvn_mask calculates the cumulative Normal   *
   * distribution function cdfbvn(h,k,rho) in a different way. *
   *                                                           *
   * The proc bvn_pi2 is a help procedure based on D.R. Divgi: *
   * "Calculation of the univariate and bivariate Normal       *
   * integral", The Annals of Statistics, 1979, 903-910.       *
   *                                                           *
   * Additions are:                                            *
   * - some simplifications,                                   *
   * - more approximating coefficients                         *
   * - an asymptotical expansion for small probabilities.      *
   *                                                           *
   * Accuracy should be over 10 digits absolute and 3-4 digits *
   * relative. Speed is similar to that of the current Gauss   *
   * cdfbvn.                                                   *
   *                                                           *
   * These procedures may be freely distributed, but not sold. *
   * I provide it hoping that it may be useful. But I disclaim *
   * any reliability for possible errors.                      *
   *                                                           *
   * Martin van der Ende                                       *
   *                                                           *
   * May 28, 1996.                                             *
   *                                                           *
   * University of Amsterdam                                   *
   * vakgroep algemene economie                                *
   * Roeterssstraat 11                                         *
   * 1018 WB Amsterdam                                         *
   * The Netherlands                                           *
   * email: mvdende@butler.fee.uva.nl                          *
   * fax: +31.20.5254254                                       *
   *                                                           *
   ************************************************************* */


proc bvn_pi2 (h,k,rho,h2,k2,rho2,rSq,rDens);
   /* bvn_pi2 returns cdfbvn(-h,-k,rho).
      cdfbvn(h,k,rho) returns normpr2(-h,-k,rho,-h2,-k2,rho2,rSq,rDens).
      Definitions:
      rho2   = sqrt(1-rho^2);
      h2     = (h-rho.*k)./rho2;
      k2     = (k-rho.*h)./rho2;
      rSq    = h^2 + k2^2;
      rDens  = exp(-rSq/2)/(2*pi);
      NOTE: h=rSina1,k2=rCosa1, k=rSina2,h2=rCosa2.
      NOTE: k2>=0 implies one can take -pi/2<=angle1<=pi/2
            h2>=0 implies one can take -pi/2<=angle2<=pi/2 */

   local vweight,angle,Shift,i,jkodd,jkeven,Uph,Upk,wx,corr0,park,tan2,mtan1;


   /* SPECIAL (SCALAR) CASES: */
/* if     maxc(maxc(abs(rho))) > 1;
          print "abs(correlation coefficient)>1"; Shift=-1;
   elseif rho==0;  Shift= Uph.*Upk;
   elseif rho==1;  Shift= Uph.*(h.>k) + Upk.*(h.<=k);
   elseif rho==-1; Shift= Uph+Upk-1; Shift= Shift.*(Shift.>0);
   else; */
     /* GENERAL CASES: */


/* vweight=zeros(21,1);
vweight[1]=  1.253314137222678;     vweight[2]=  -0.999999996322367;
vweight[3]=  0.6266570144788616;    vweight[4]=  -0.333332913298024;
vweight[5]=  0.1566622679316691;    vweight[6]=  -6.666028411513941E-2;
vweight[7]=  2.609623581290622E-2;  vweight[8]=  -9.499531778787081E-3;
vweight[9]=  3.232824755470388E-3;  vweight[10]= -1.027330989192324E-3;
vweight[11]= 3.020025482658464E-4;  vweight[12]= -8.068099400127647E-5;
vweight[13]= 1.912225364901334E-5;  vweight[14]= -3.912070127260670E-6;
vweight[15]= 6.711727217209611E-7;  vweight[16]= -9.365778519087949E-8;
vweight[17]= 1.026548544344371E-8;  vweight[18]= -8.449983408562188E-10;
vweight[19]= 4.887156068638151E-11; vweight[20]=-1.7641960413010695E-12;
vweight[21]= 2.982043919845545E-14; */

 vweight=zeros(31,1);
vweight[1]=  1.253314137315500;     vweight[2]=  -0.9999999999999754;
vweight[3]=  0.6266570686571231;    vweight[4]=  -0.3333333333248598;
vweight[5]=  0.1566642670935284;    vweight[6]=  -6.666666626444445E-2;
vweight[7]=  2.611070955316157E-2;  vweight[8]=  -9.523804503777532E-3;
vweight[9]=  3.263827020322102E-3;  vweight[10]= -1.058178793302071E-3;
vweight[11]= 3.263502098286853E-4;  vweight[12]= -9.615835274820773E-5;
vweight[13]= 2.715573937406461E-5;  vweight[14]= -7.363023032007412E-6;
vweight[15]= 1.915812162821294E-6;  vweight[16]= -4.766013394595678E-7;
vweight[17]= 1.125080846091322E-7;  vweight[18]= -2.491722870193576E-8;
vweight[19]= 5.103898425245527E-9;  vweight[20]= -9.517178862491669E-10;
vweight[21]= 1.589460355702903E-10; vweight[22]= -2.339412414181092E-11;
vweight[23]= 2.985763266418183E-12; vweight[24]= -3.249122467023523E-13;
vweight[25]= 2.958566939223091E-14; vweight[26]= -2.204022324089611E-15;
vweight[27]= 1.304589139468861E-16; vweight[28]= -5.887558091126596E-18;
vweight[29]= 1.899221719494747E-19; vweight[30]= -3.893896136973921E-21;
vweight[31]= 3.807675366775568E-23;


/* vweight=zeros(41,1);
 vweight[1]= 1.253314137315500;     vweight[2]=  -1;
 vweight[3]= 0.6266570686577501;    vweight[4]=  -0.3333333333333332;
 vweight[5]= 0.1566642671644364;    vweight[6]=  -6.666666666665691E-2;
 vweight[7]= 2.611071119401255E-2;  vweight[8]=  -9.523809523527702E-3;
 vweight[9]= 3.263838898234508E-3;  vweight[10]= -1.058201055227967E-3;
 vweight[11]=3.263838829052882E-4;  vweight[12]= -9.620008249423181E-5;
 vweight[13]=2.719863508376179E-5;  vweight[14]= -7.399976375870648E-6;
 vweight[15]=1.942724550073706E-6;  vweight[16]= -4.932964393369089E-7;
 vweight[17]=1.213895545910440E-7;  vweight[18]= -2.899419757469565E-8;
 vweight[19]=6.728519264552478E-9;  vweight[20]= -1.517131942551575E-9;
 vweight[21]=3.319108847287039E-10; vweight[22]= -7.022141565916787E-11;
 vweight[23]=1.428684402852183E-11; vweight[24]= -2.773496141377873E-12;
 vweight[25]=5.088537684255269E-13; vweight[26]= -8.730508302470239E-14;
 vweight[27]=1.385595302709042E-14; vweight[28]= -2.012439906068720E-15;
 vweight[29]=2.647206396468263E-16; vweight[30]= -3.122000822598359E-17;
 vweight[31]=3.267760772243843E-18; vweight[32]= -3.003526377626578E-19;
 vweight[33]=2.396165106310595E-20; vweight[34]= -1.637042915419164E-21;
 vweight[35]=9.422353450786448E-23; vweight[36]= -4.474458136625314E-24;
 vweight[37]=1.704342515274355E-25; vweight[38]= -4.999746081902462E-27;
 vweight[39]=1.059282556071995E-28; vweight[40]= -1.441337407752791E-30;
 vweight[41]=9.450854298284086E-33; */



     /* ANGLE REDUCTION TO (-PI/2, PI/2]:
        L(h,k;rho)
      = cdfnc(h) - L(h,-k;-rho)                  if ((k2.<0) .and (h2.>=0))
      = cdfnc(k) - L(-h,k;-rho)                  if ((h2.<0) .and (k2.>=0))
      = -1 + cdfnc(k) + cdfnc(h) + L(-h,-k; rho) if ((h2.<0) .and (k2.<0)). */

     /* APPLY ANGLE REDUCTION TO (-PI/2, PI/2] IF: */
     jkodd  = h2.<0;
     jkeven = k2.<0;

     /* ADD 'Shift' AFTER APPLYING ANGLE REDUCTION */
     Uph    = cdfnc(h);
     Shift  = Uph.*jkeven;
     Upk    = cdfnc(k);
     Shift  = Shift + Upk.*jkodd;
     Shift  = Shift - (jkodd .and jkeven);

     /* GET COSINUS AND MINUS SINUS OF REDUCED ANGLES */
     h      = h.*(1 - 2*jkeven);
     k      = k.*(1 - 2*jkodd);
     h2     = ABS(h2);
     k2     = ABS(k2);

     /* ADJUST CORRELATION */
     corr0  = ABS(rho).<1E-5;

     /* GET REDUCED ANGLE */
     angle  = -atan(rho2./(rho+corr0));
/*   angle  = -atan2(rho2,rho); */
     angle  = angle + pi*( ((h.>0) .and (k.>0) .and (angle.<0)) -
                           ((h.<0) .and (k.<0) .and (angle.>0)) );
/*   NECESSARY IF YOU REDUCED ANGLE1, ANGLE2 TO ONLY PI/2 IN MODULUS */


     /* INITIALIZE THE POLYNOMIAL APPROXIMATION */
     jkeven = angle;                                    /* now jkeven=j0 */
     jkodd  = k+h;                                     /* now jkodd=j1 */

     tan2   = h./(k2+1E-300);                  /* = tan(angle2) */
     mtan1  = k./(h2+1E-300);                  /* = -tan(angle1) */
     rho=0;rho2=0; /* clear */


     /* GET THE POLYNOMIAL APPROXIMATION */
     wx  = jkeven - vweight[1]*jkodd;
     i=3;
     do until i>31;
        h = h.*k2; k=k.*h2;
        jkeven = (h+k + (i-2)*rSq.*jkeven)/(i-1);
        h = h.*k2; k=k.*h2;
        jkodd  = (h+k + (i-1)*rSq.*jkodd)/i;
        wx  = wx-vweight[i]*jkodd-vweight[i-1]*jkeven;

        i=i+2;
     endo;

     h=0; k=0; h2=0; k2=0; /* clear */

     /* APPLY AN ASYMPTOTIC EXPANSION IF NEEDED:
          rDens.*(a2-a1 - int_a1^a2 1 - 1/(Rcos a)^2 + 3/(Rcos a)^4 -
                                       15/(Rcos a)^6 da
        = rDens.*((tan a)./R^2 - (3./R^4).*[(1/3)*(tan a)^3 + tan(a)]_a1^a2 +
                  (15./R^6).*[(1/5)*(tan a)^5 + (2/3)*(tan a)^3 + tan(a)]_a1^a2
        = rDens.*((tan a2 - tan a1).*(1./R^2 -3./R^4 + 15./R^6)  -
                  (tan^3 a2 - tan^3 a1).*(1./R^4 - 10./R^6) +
                  (tan^5 a2 - tan^5 a1).*(3./R^6) */

     jkodd = tan2^2; jkeven = mtan1^2;
     rSq = rSq+(rSq.==0);

/*   ONE-TERM ASYMPTOTIC EXPANSION
     angle = (tan2+mtan1)./rSq; */

/*   THREE-TERM ASYMPTOTIC EXPANSION
     park  = rSq^2 - 3*rSq + 15;
     park  = (tan2+mtan1).*park;
     tan2  = tan2.*jkodd; mtan1=mtan1.*jkeven;
     park  = park - (tan2+mtan1).*(rSq-10);
     tan2  = tan2.*jkodd; mtan1=mtan1.*jkeven;
     park  = park + (tan2+mtan1)*3;
     park  = park./rSq^3; */


/*   FIVE-TERM ASYMPTOTIC EXPANSION
     park  = 945 - rSq.*(105-rSq.*(15-rSq.*(3-rSq)));
     park  = (tan2+mtan1).*park;
     tan2  = tan2.*jkodd; mtan1=mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(1260 - rSq.*(105-rSq.*(10-rSq)));
     tan2  = tan2.*jkodd; mtan1=mtan1.*jkeven;
     park  = park + 3*(tan2+mtan1).*(378-rSq.*(21-rSq));
     tan2  = tan2.*jkodd; mtan1=mtan1.*jkeven;
     park  = park + 15*(tan2+mtan1).*(36-rSq);
     tan2  = tan2.*jkodd; mtan1=mtan1.*jkeven;
     park  = park + (tan2+mtan1)*105;
     park  = park./rSq^5; */


/*   SEVEN-TERM ASYMPTOTIC EXPANSION
     park  = 135135 -
             rSq.*(10395-rSq.*(945-rSq.*(105-rSq.*(15-rSq.*(3-rSq)))));
     park  = (tan2+mtan1).*park;
     tan2  = tan2.*jkodd; mtan1=mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(270270 -
             rSq.*(17325-rSq.*(1260-rSq.*(105-rSq.*(10-rSq)))));
     tan2  = 3*tan2.*jkodd; mtan1=3*mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(135135 -
             rSq.*(6930-rSq.*(378-rSq.*(21-rSq))));
     tan2  = 5*tan2.*jkodd; mtan1=5*mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(25740 -
             rSq.*(990-rSq.*(36-rSq)));
     tan2  = 7*tan2.*jkodd; mtan1=7*mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(2145-rSq.*(55-rSq));
     tan2  = 9*tan2.*jkodd; mtan1=9*mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(rSq-78);
     tan2  = 11*tan2.*jkodd; mtan1=11*mtan1.*jkeven;
     park  = (park+tan2+mtan1)./rSq^7; */

/*   NINE-TERM ASYMPTOTIC EXPANSION */
     park  = 34459425-
             rSq.*(20270-rSq.*(135135-rSq.*(10395-rSq.*(945-rSq.*
                  (105-rSq.*(15-rSq.*(3-rSq)))))));
     park  = (tan2+mtan1).*park;
     tan2  = tan2.*jkodd; mtan1=mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(91891800 -
             rSq.*(4729725-rSq.*(270270-rSq.*(17325-rSq.*(1260-rSq.*
                  (105-rSq.*(10-rSq)))))));
     tan2  = 3*tan2.*jkodd; mtan1=3*mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(64324260 -
             rSq.*(2837835-rSq.*(135135-rSq.*(6930-rSq.*(378-rSq.*
                  (21-rSq))))));
     tan2  = 5*tan2.*jkodd; mtan1=5*mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(18378360 -
             rSq.*(675675-rSq.*(25740-rSq.*(990-rSq.*(36-rSq)))));
     tan2  = 7*tan2.*jkodd; mtan1=7*mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(2552550 -
             rSq.*(75075-rSq.*(2145-rSq.*(55-rSq))));
     tan2  = 9*tan2.*jkodd; mtan1=9*mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(185640 -
             rSq.*(4095-rSq.*(78-rSq)));
     tan2  = 11*tan2.*jkodd; mtan1=11*mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(7140 -
             rSq.*(105-rSq));
     tan2  = 13*tan2.*jkodd; mtan1=13*mtan1.*jkeven;
     park  = park + (tan2+mtan1).*(136-rSq);
     tan2  = 15*tan2.*jkodd; mtan1=15*mtan1.*jkeven;
     park  = (park+tan2+mtan1)./rSq^9;

     wx=wx - (wx-park).*(rSq.>70);

     /* GET THE RETURN VALUE FOR NONZERO CORRELATION */
     wx=Shift + rDens.*wx;

     /* REPLACE THE RETURN VALUE FOR ZERO CORRELATION */
     wx=wx + (Uph.*Upk-wx).*corr0;

/* endif; */
   retp(wx);
endp;



proc bvn_mask(h,k,rho);
/* calculates cdfbvn in a different way */
local h2,k2,rho2,rSq,rDens;

  rho2=sqrt(1-rho^2)+1E-300;
  h2=(h-rho.*k)./rho2;      /* Pr(x=h conditional on y=k)=Pr(h2) */
  k2=(k-rho.*h)./rho2;      /* Pr(y=k conditional on x=h)=Pr(k2) */
  rSq=h^2 + k2^2;           /* square of distance of wedge vortex to 0 */
  rDens=exp(-rSq/2)/(2*pi); /* rho2 times pdfn(h,k,rho) */

  retp(bvn_pi2(-h,-k,rho,-h2,-k2,rho2,rSq,rDens));

endp;

************************************************************************
Economics       is entirely responsible for the information
provided above. Please direct any comments to Information Provider at:
          E-Mail Address: econ@american.edu
          Phone Number:   (202) 885-3770
