st: simluated errors in multivariate probit model

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

*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 has values of 1 on the leading diagonal and correlations ρjk = ρkj
as off-diagonal elements*/

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 varriance is different from the original variance and covariance
matrix "V" 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!



-- 
Zhi Su
348 Holmes Hall
Northeastern University
360 Huntington Avenue
Boston, MA 02115
Office:1-617-373-2316
email:su.zh@husky.neu.edu

*
*   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/