st: get expected error terms in multivariate probit model

[Zhi Su <su.zh@husky.neu.edu>]

I want to get expected error terms in multivariate probit model.
What I do is draw errors from N(0,V), where V is the variance and
covariance matrix calculated the simulated Y through predicting latent
Y = predicted XB+random-drawn errors
keep errors that satisfy simulated Y=Actual Y
But the variance and covariance matrix of simulated errors are
different from the one estimated by mvprobit.
Can anyone help me to review the codes and tell me why?

Here are my codes

*estimate the multivariate probit model using mvprobit,`basic' is a set of X
mvprobit (T1=`basic') (T2=`basic') (T3=`basic') (T4=`basic'),robust

*save the xb
mvppred xbm,xb

*start mata workplace
mata

/*number of repetions*/
N=3000

/*number of observations*/
O=677

/*I get variance and covariance matrix from MVProbit
/*I input the values from results of MVPROBIT*/

V=(1,0.6524083,0.6021572,0.7241811\0.6524083,1,0.4410809,0.421958\0.6021572,0.4410809,1,0.4515065\
0.7241811,0.421958,0.4515065,1)

/*build a column vector and fill in values of means of simulated error
terms for each observation, which is calculated by the following
process*/
M_u=J(O,1,.)

/*build a matrix and fill in values of means of simulated error terms
for each observation*/
M_m=J(O,4,.)

/* Loop through the simulations for 677 obsertions used in calculation
of mvprobit*/
for(k=1;k<=O;k++) {

/*Multiple draws from multivariate normal distribution with variance
and covariance V*/
/*the matrix is drawn from N(0,V)*/
U=(invnormal(uniform(N,cols(V)))*cholesky(V)')

/*create a copy of data for each observations (from 1 to 677) into
mata workplace*/
Y=st_data(k,("T1","T2","T3","T4"))

XBM=st_data(k,("xbm1", "xbm2", "xbm3","xbm4"))

/*generate simulated latent Y matrix*/
Y_si=XBM:+U

/* make an N by 4 matrix of 1 to prepare for a matrix of Y_si index*/
Y_in= J(N, 4, 1)

/* loop through the rows and columns of the matrix Y_si, setting the
value in a cell of Y_in*/
/* equal to 0 if Y_si is less than 0*/

for(i = 1; i <= rows(Y_si); i++) {
      for(j = 1; j<=cols(Y_si); j++) {
         	if (Y_si[i,j]<0) {
                  Y_in[i,j]=0
         	}
       }
}

/* make an N by 1 matrix of 0 to prepare for Y_in matrix selection rule*/
D= J(N,1, 0)

/*loop through each row of Y_in and compare it with actual Y*/
/*set corresponding row of D to 1 if two rows are equal*/

for(i = 1; i <= rows(Y_si); i++) {
   if  (Y_in[i,.]==Y) {
      D[i]=1
    }
}

/* only simulated error U that satisfy Y_in=actual Y are obtained,
when rows of D==1*/
U_s=select(U, D:==1)

/*column means of simulated error matrix U_s*/
M_c=mean(U_s)

/*fill in rows of M_m*/
M_m[k,.]=M_c

/*store the calculated values in the Stata variable "u*"*/
st_store((1,rows(M_m)),"u1", M_m[.,1])
st_store((1,rows(M_m)),"u2", M_m[.,2])
st_store((1,rows(M_m)),"u3", M_m[.,3])
st_store((1,rows(M_m)),"u4", M_m[.,4])

/*store the calculated values in the Stata variable "res_si"*/
st_store((1,rows(M_u)),"res_si",M_u)

/*end mata worksplace*/
end

*I standardize u1 u2 u3 u4

egen Su1=std(u1)
egen Su2=std(u2)
egen Su3=std(u3)
egen Su4=std(u4)

corrmat Su1 Su2 Su3 Su4

           Su1        Su2        Su3        Su4
Su1          1
Su2  .81469367          1
Su3  .77865427  .62243427  .99999999
Su4  .90029139  .65351567  .66685678          1

The variance and coviance matrix is different from the original one I set, why?
          u1           u2          u3     u4
u1         1
u2 0.6524083            1
u3 0.6021572    0.4410809          1
u4 0.7241811    0.421958   0.4515065       1

Are there any mistakes in the codes?
Thank you!
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/