For information of trivariate distribution to bivariate distribution see,

http://www.sci.wsu.edu/math/faculty/genz/papers/bvnn/node3.html


proc iml;

  start fun(t) global (b2, b3, rho12, rho13, rho23);
      x1=(b2 - rho12#t)/sqrt((1-rho12#rho12));
x2=(b3 - rho13#t)/sqrt((1-rho13#rho13));
      rho=(rho23 - rho12#rho13)/sqrt((1-rho12#rho12)#(1-rho13#rho13));
   f=probbnrm(x1,x2,rho);
      v=exp(-0.5 # (t # t) ) # f;
      return(v);
  finish;

  b1=0;
  b2=0;
  b3=0;
  rho12=0.99992;
  rho13=0.64627;
  rho23=0.63975;
  *print b1 b2 b3 rho12 rho13 rho23;
  interval= .M || b1;
  call quad(z,"fun",interval);
  prob2=z/(sqrt(2#constant('pi')));
  prob_theory=0.125 +(0.25/constant('pi'))*(arsin(rho12) + arsin(rho13) + arsin(rho23) );

  
  print prob2[format=E21.14] prob_theory[format=E21.14];
  quit;