data a;
      keep y1 y2 x1 x2;
      do i = 1 to 500;
         x1 = rannor( 19283 );
         x2 = rannor( 98721 );
         u1 = rannor( 76527 );
         u2 = rannor( 65721 );
         y1l = 1 + 2 * x1 + 3 * x2 + u1;
         y2l = 3 + 4 * x1 - 2 * x2 + u1*.2 + u2;
         if ( y1l > 0 ) then y1 = 1;
         else                y1 = 0;
         if ( y2l > 0 ) then y2 = 1;
         else                y2 = 0;
         output;
      end;
   run;

   /*-- Bivariate Probit --*/
   proc qlim data=a method=qn;
      init y1.x1 2.8, y1.x2 2.1, _rho .1;
      model y1 = x1 x2;
      model y2 = x1 x2;
      endogenous y1 y2 ~ discrete;
   run;