GHK simulation for multivariate probit model

[Shoummo Sen Gupta]

Dear all,
I am using the following program for fitting a trivariate probit model using maximum likelihood method (ml maximize).
The detailed description of the program is given in https://journals.sagepub.com/doi/pdf...867X0600600202.

program define myll
args lnf x1 x2 x3 c21 c31 c32
tempvar sp2 sp3 k1 k2 k3
quietly {
gen double ‘k1’ = 2*$ML_y1 - 1
gen double ‘k2’ = 2*$ML_y2 - 1
gen double ‘k3’ = 2*$ML_y3 - 1
tempname cf21 cf22 cf31 cf32 cf33 C1 C2
su ‘c21’, meanonly
scalar ‘cf21’ = r(mean)
su ‘c31’, meanonly
scalar ‘cf31’ = r(mean)
su ‘c32’, meanonly
scalar ‘cf32’ = r(mean)
scalar ‘cf22’ = sqrt(1 - ‘cf21’^2)
scalar ‘cf33’ = sqrt(1 - ‘cf31’^2 - ‘cf32’^2)
mat ‘C1’ = (1, 0 , 0 \ ‘cf21’, ‘cf22’, 0 \ ‘cf31’ , ‘cf32’ , ‘cf33’)
mat ‘C2’ = (1, 0 \ ‘cf21’, ‘cf22’)
egen ‘sp3’ = mvnp(‘x1’ ‘x2’ ‘x3’) if
$ML_y1==1, ///
chol(‘C1’) dr($dr) ml prefix(z) signs(‘k1’ ‘k2’ ‘k3’)
egen ‘sp2’ = mvnp(‘x1’ ‘x2’ ) if
$ML_y1==0, ///
chol(‘C2’) dr($dr) ml prefix(z) signs(‘k1’ ‘k2’)
replace ‘lnf’= cond($ML_y1, ln(‘sp3’), ln(‘sp2’), .)
}
end

My question is that whenever I am executing the simulation using ml maximize, I am getting the following error.

ml maximize][Initial: Log likelihood = -57130.947

Rescale: Log likelihood = -57130.947
Rescale eq: Log likelihood = -57074.341
Iteration 0: Log likelihood = -57074.341 (not concave)
Iteration 1: Log likelihood = -56074.987
Iteration 2: Log likelihood = -55879.375 (not concave)
Iteration 3: Log likelihood = -55877.227 (not concave)
Iteration 4: Log likelihood = -55876.993 (not concave)
Iteration 5: Log likelihood = -55876.9 (not concave)
Iteration 6: Log likelihood = -55876.855 (not concave)
Iteration 7: Log likelihood = -55876.462 (not concave)
Iteration 8: Log likelihood = -55874.438 (not concave)
Iteration 9: Log likelihood = -55873.21 (not concave)
Iteration 10: Log likelihood = -55871.72 (not concave)
Iteration 11: Log likelihood = -55870.165 (not concave)
Iteration 12: Log likelihood = -55868.719 (not concave)
Iteration 13: Log likelihood = -55868.593 (not concave)
Iteration 14: Log likelihood = -55866.979 (not concave)
Iteration 15: Log likelihood = -55866.221 (not concave)
Iteration 16: Log likelihood = -55865.626 (not concave)
Iteration 17: Log likelihood = -55865.286 (not concave)
Iteration 18: Log likelihood = -55865.044 (not concave)
Iteration 19: Log likelihood = -55864.622 (not concave)
Iteration 20: Log likelihood = -55863.761 (not concave)
Iteration 21: Log likelihood = -55863.667 (not concave)
Iteration 22: Log likelihood = -55863.462 (not concave)
Iteration 23: Log likelihood = -55863.237 (not concave)
Iteration 24: Log likelihood = -55863.055 (not concave)
Iteration 25: Log likelihood = -55862.994 (not concave)
Iteration 26: Log likelihood = -55862.949 (not concave)
Iteration 27: Log likelihood = -55862.879 (not concave)
Iteration 28: Log likelihood = -55862.851 (not concave)
cannot compute an improvement -- discontinuous region encountered
r(430);

Can you please point out the source of this error?

Thank you.