I am doing a Monte Carlo simulation involving CML (constrained maximum likeliohood) procedure. In most cases, the codes are running well, however in certain DGP the inverse of a matrix, when calculating the log-likelihood, does not exist (non positive definiteness). I have set up a flag in the main body of the codes such that if there is any error in any particular Monte Carlo sample , that sample is dropped and new sample is drawn again. Now I set up trapchk and the trap in the following way
oldval = trapchk(1); trap 1,1; rhoi=invpd(rho); trap oldval,1; if scalerr(rhoi); print "wrong rho" rho; go to fap; endif;
Where Fap is the flag I set earlier. But it does not work. I tried with _fap, no luck. Can anybody help please?
Here is a simple working example that tries to convert the inverse of a non-positive definite matrix with the GAUSS command invpd. (NOTE that you can use the GAUSS command inv to compute the inverse of a non-positive definite matrix. Though, I assume in the context of your program this matrix must be positive definite).
//Create label mylabel: //non-positive definite matrix x = reshape(seqa(1, 1, 16), 4, 4); //store previous trap state old_trap = trapchk(1); trap 1; x_inv = invpd(x); trap old_trap; if scalinfnanmiss(x_inv); print "matrix not positive definite"; goto mylabel; endif;
It might be more clear to use a do while loop in this case, however. For example:
//Flag to track if matrix is positive definite matrix_ok = 1; do while (matrix_ok == 1); //non-positive definite matrix x = reshape(seqa(1, 1, 16), 4, 4); //store previous trap state old_trap = trapchk(1); trap 1; x_inv = invpd(x); trap old_trap; if scalinfnanmiss(x_inv); print "matrix not positive definite"; //Set flag to indicate matrix not pos. def. matrix_ok = 0; endif; endo;
Many thanks for this. But, suppose your non-positive definiteness problem appears within the user defined procedure (such as log-likelihood function) and you want to stop the program when the problem is identified by "trap" and resume again where it stopped.
To be more specific, suppose you are doing 2000 Monte Carlo simulations. You have identified the problem at 1401 replication (in your procedure). Then you want to drop this particular 1401-th sample, but keeping the previous 1400 replications results and drawing a new sample for 1401-th replication. To do this, I think I need to set a flag at the beginning of the simulation.
How can I use this flag within my procedure? or can I do that?
Many thanks in advance!
Are you using a loop for your simulation? Then use break. Doing it they way I have presented below you will get 2000 good simulations. The ones with singular rhos won't get counted.
simul = 1; do until simul > 2000; oldval = trapchk(1); trap 1,1; rhoi=invpd(rho); trap oldval,1; if scalerr(rhoi); print "wrong rho" rho; break; endif; . . . // cmlmt run . . simul = simul + 1; endo;