When an error in a selection equation is correlated with regression equations, then these correlation measures have to be corrected in the likelihood function, otherwise the results will be biased.

Here are a couple of procedures one can use.
1) nlmixed
2) qlim

here is a link for the further reference ( first 2 pages)
http://siteresources.worldbank.org/DEC/Resources/movestay.pdf

**************************************************;
%let n=500;

data Sigma (type=cov);
infile cards;
input _type_ $ _Name_ $ errz err0 err1;
cards;
cov errz 1 0.8 0.5
cov err0 0.8 1 0.3
cov err1 0.5 0.3 1
;
run;

proc simnormal data=sigma out=err_mat numreal=&n seed=134;
var errz err0 err1;
run;


data sim;
  set err_mat;
   x0=rannor(123);
   x1=rannor(123);
   z=rannor(123);
   s=1+2*z>errz;
   if s=0 then y=2+2*x0+1.41*err0;
   else        y=1+3*x1+0.2*err1;
   if s=0 then y0=y;
   else y1=y;

run;

proc nlmixed data=sim;
parms a1 a2 b1 b2 z1 z2=1 s0=1 s1=0.5 r0 r1=0.01;
/*s0=1; s1=0.2;*/
err0=(y-a1-a2*x0);
err1=(y-b1-b2*x1);

pdf0=pdf('normal',err0,0,s0);
pdf1=pdf('normal',err1,0,s1);

if s=0 then do;
    zbeta=(z1+z2*z+r0*(err0/s0))/sqrt(1-r0**2);
    cdf=cdf('normal',zbeta);
    prob=cdf;
    lik=(1-prob)*pdf0;
end;
else   do;
    zbeta=(z1+z2*z+r1*(err1/s1))/sqrt(1-r1**2);
    cdf=cdf('normal',zbeta);
    prob=cdf;
    lik=(prob)*pdf1;
end;
ll=log(lik);
model s ~ general(ll);
run;

proc qlim data=sim;
   model s = z/discrete;
   model y0 = x0 / select(s=0);
   model y1 = x1 / select(s=1);
run;