# 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?

• are BetaH, BetaM, BetaL also parameters (i.e. things to estimate)? Sep 20, 2010 at 18:31
• Kwak -- yes; I'll clarify that in the question statement.
– kyle
Sep 20, 2010 at 18:40
• I looks to me like you should only need 3 guesses, one for when each of $H'$, $M'$ and $L'$ is largest. But, I may simply misunderstand the question. Sep 20, 2010 at 20:55

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.

Your likelihood function is non-concave (i.e. the Hessian matrix of your likelihood function is not SDN). From this it follows that

1. You will only find a local maximae to your likelihood function (no garantuee of global optimality)
2. This maxima will always depend on your choice of starting point.
3. 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.

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.

• > this is not the same model (i.e. Kelly's model is non linear in the parameters). Sep 22, 2010 at 22:35
• Paul-- believe it or not, we tried this as a starting point for the optimization (taking into account the concerns you raised), but the estimator was way off from where we'd expect it to converge to. No good explanation why. We do have a reasonable workaround, which I'll be sharing with the question shortly.
– kyle
Sep 23, 2010 at 1:15
• @kwak:> what do you mean by "this"? Not the same model as 2sls or not the same model as what the GHK simulator is usually applied to, ie not a multinomial probit? I know it's not a multinomial probit, that's why I suggested a slight variant of the GHK simulator. As to whether 2sls applies, doesn't the equation for W have some error term, say, e, and doesn't the full model imply E[e|Z,X] = 0 and E[(H,M,L)|Z] not 0?
– paul
Sep 23, 2010 at 1:17
• > you state You could avoid the problem altogether by simply estimating which is certainly correct. Nonetheless, the model you propose is a different one than the one Kyle wanted to fit (it is not a reformulation of the same model). Sep 23, 2010 at 7:37