data a;
   keep y1 y2 x1 x2;
   do i = 1 to 500;
      x1 = rannor( 19283 );
      x2 = rannor( 19283 );
      u1 = rannor( 19283 );
      u2 = rannor( 19283 );
      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;