How can I work around "lumpiness" in simulated maximum likelihood estimation? I am attempting to estimate a model of the following form:
W = alphaH * H + alphaM * M + alphaL * L + X * beta

where H, M, L are indicators for a discrete choice variable, and beta is something like 35-dimensional. Because we believe our data/model has endogeneity issues, we have expanded the model to
W  = alphaH * H' + alphaM * M' + alphaL * L' + X * beta
H  = Z * betaH
M  = Z * betaM
L  = Z * betaL
H' = 1( H = max(H,M,L) )
M' = 1( M = max(H,M,L) )
L' = 1( L = max(H,M,L) )

where Z are instruments, and betaH, betaM, betaL are parameters to be estimated. This "subregression" corresponds to a latent utility-based choice model.
We have been able to estimate the second-stage model (estimates of H, M, L, implying H', M', L') in Stata using the mvprobit command, but can't figure out how to estimate the entire model in one fell swoop. To work around this, we wrote some code in MATLAB to estimate the model using simulated maximum likelihood, but MATLAB is choking on local minima (maxima in this problem, but MATLAB will only minimize the negative...), of which there are plenty.
We have attempted to work around this by starting from a few dozen initial conditions, none of which usually converges to the "right" answer; I say this with near certainty since we have been testing the code piecewise and have confirmation (on randomized test data) that if the optimization starts near the "correct" values (in test), it converges to reasonable values, otherwise it gets nowhere close (although the resultant outcome has a far lower overall likelihood).
Are there any tricks -- MATLAB, Stata, or otherwise -- to work around this problem? Is this an inherent issue with simulation versus closed-form analysis?
Thanks for your help.
 A: It sounds like you need to use a more robust optimization algorithm that can handle local minima. Particle swarm methods work quite well in this case. Or you could try other evolutionary optimization methods or simulated annealing.
A: Your likelihood function is non-concave (i.e. the Hessian matrix of your likelihood function is not SDN). From this it follows that


*

*You will only find a local maximae to your likelihood function (no garantuee of global optimality)

*This maxima will always depend on your choice of starting point.

*Your maximization procedure will always be an iterative one.


Without directly solving the issues aboves, one way to handle them would be thru Monte carlo optimization (a review is here) which is basically a recasting of Rob's suggestion inside the frame of statistics.
A: You could avoid the problem altogether by simply estimating
W = alphaH * H + alphaM * M + alphaL * L + X * beta
using 2sls. The fact that H,M, and L are discrete doesn't violate any of the assumptions of 2sls. Of course, using maximum likelihood will produce more efficient estimates, but it relies on more assumptions. If nothing else, the 2sls estimates should provide good starting values for you maximization algorithm. 
For maximizing the likelihood,you should try changing your simulation method to make the likelihood function smooth. I think a slight variant of the Geweke-Hajivassiliou-Keane (often just GHK) simulator would work. 
