Here is a complete solution. I put some comments for easy understanding.

***create corr error data***;

data corr_err (type=cov);
infile cards;
input _type_ $ _Name_ $ e1 e2;
cards;
cov e1 1 0.8
cov e2 0.8 1
;
run;

proc simnormal data=corr_err out=binormal numreal = 1000 seed=999087;
  var e1 e2;
run;

***check the correlation**;
proc corr data=binormal;
var e1 e2;
run;

***generate simulation data**;
data simudata;
  set binormal;
  x1=rannor(123);x2=rannor(123);
  y1=1+2*x1>e1;
  y2=1+1*x2>e2;

  if y1=1 and y2=1 then z=1;
  else z=0;
  run;

  proc freq data=simudata;
  table y1*y2 z/list;
  run;

**maximum likelihood estimation***;
proc nlmixed data=simudata(keep=x1 x2 z);
parms a1 b1 a2 b2 =0.1 rho=0.5;

bounds -0.99<=rho<=0.99;

xbeta1=a1+b1*x1;
xbeta2=a2+b2*x2;

p=PROBBNRM(xbeta1,xbeta2,rho);

if z=1 then cdf=p;
else        cdf=1-p;

ll=log(cdf);

model z ~ general(ll);

CONTRAST "True value b1=2" b1-2 ;
CONTRAST "True value b2=1" b2-1 ;
CONTRAST "rho=0.8" rho-0.8 ;

estimate "True value b2-b1=-1" b2-b1;

run;

不好意思 之前账户密码掉了 非常谢谢你

谢谢您的耐心解答，我还是有不懂的地方，比如您看我EXCEL中的实际数据，Z 代表的意思是：Zi=1 indicates firm-ihas committed fraud and has been detected, and Zi=0 if firm-i hasnot committed fraud or has committed fraud but has not been detected.

P(Zi=1)=P(Fi=1and Di=1)=P(Fi=1) P(Di=1|Fi=1)= Φ(XF,iβF,XD,iβD)

在现实生活中，wedo not directly observe the realizations of Fi and Di.What we observe is Zi= Fi* Di.

Fi= XF,iβF +μi;

Di= XD,iβD +ʋi;

Where XF,i is a rowvector with elements that explain firm-i’s incentive to commit fraud, and XD,icontains variables that explain the firm’s potential for getting caught. μiand ʋi are zero-mean disturbances with a bivariate normaldistribution.

For fraud occurrence, I transform  into a binaryvariable Fi, where Fi=1 if Fi* >0, and Fi =0 otherwise. For frauddetection (conditional on occurrence), I transform   into a binaryvariable Di, where Di=1 if Di*>0, and Di=0 otherwise.

我用variable r1 r2 r3 r4 r6 r7 r9 explain Fi, variable Br1 Br2 Br3 Br4 Er3 Er4 Er5 explian Di.

可是问题是，我们只有Z的数据，如何来确定Fi和Di呢？

如果确认了Fi和Di,就可以用proc qlim来做了。是我想的有问题吗？希望您能帮忙解答。

请问你的问题解决了吗，能稍微解释一下吗？

您好,
首先感谢bobguy提供了这么宝贵的范本让我有了方向,
只是,我尝试了放入我的variables后发现出现了些讯息
基本上,我的变数为X1~X28,
并且其中可拆成2类.
XBETA1可用X1~X15.
XBETA2目前我是放X1~X28,这部分也有尝试过只放X16~X28...
我的尝试如下:

PROC NLMIXED DATA=C1(KEEP=X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 X18 X19 X20 X21 X22 X23 X24 X25 X26 X27 X28 Y06 Y07 Y08 Y09 Y10 Y11 IND_STEEL IND_FOOD IND_PLASTIC IND_TEXTILE IND_ELECTRIC IND_ELECTRON IND_TOUR Z);

PARMS A1 B1 C1 D1 E1 F1 G1 H1 I1 J1 K1 L1 M1 N1 O1 P1 Q1 R1 S1 T1 U1 V1 W1 X_1 Y1 Z1 BB1 CC1 DD1
A2 B2 C2 D2 E2 F2 G2 H2 I2 J2 K2 L2 M2 N2 O2 P2 Q2 R2 S2 T2 U2 V2 W2 X_2 Y2 Z2 BB2 CC2=0.1 RHO=0.5;

BOUNDS -0.99<=RHO<=0.99;

XBETA1=A1+B1*X1+C1*X2+D1*X3+E1*X4+F1*X5+G1*X6+H1*X7+I1*X8+J1*X9+K1*X10+L1*X11+M1*X12+N1*X13+O1*X14+P1*X15+Q1*Y06+R1*Y07+S1*Y08+T1*Y09+U1*Y10+V1*Y11+W1*IND_STEEL+X_1*IND_FOOD+Y1*IND_PLASTIC+Z1*IND_TEXTILE+BB1*IND_ELECTRIC+CC1*IND_ELECTRON+DD1*IND_TOUR;

XBETA2=A2+B2*X16+C2*X17+D2*X18+E2*X19+F2*X20+G2*X21+H2*X22+I2*X23+J2*X24+K2*X25+L2*X26+M2*X27+N2*X28+O2*X1+P2*X2+Q2*X3+R2*X4+S2*X5+T2*X6+U2*X7+V2*X8+W2*X9+X_2*X10+Y2*X11+Z2*X12+BB2*X13+CC2*X14;

PROBIT=PROBBNRM(XBETA1,XBETA2,RHO);
IF Z=1 THEN CDF=PROBIT;
ELSE CDF=1-PROBIT;
LL=LOG(CDF);
MODEL Z~GENERAL(LL);
RUN;

以上,结果总是跑不出来,
而SAS的LOG提供的讯息是,

NOTE: A FINITE DIFFERENCE APPROXIMATION IS USED FOR THE DERIVATIVE OF THE 'PROBBNRM' FUNCTION.
NOTE: EXECUTION ERROR FOR OBSERVATION 1.
NOTE: PROCEDURE NLMIXED USED (TOTAL PROCESS TIME):
      REAL TIME           0.66 SECONDS
      CPU TIME            0.09 SECONDS


不知道是否有哪位高人可以帮忙解惑?