How is the covariance matrix calculated in the MAXLIK?

The version of the Hessian retrieved in _max_diagnostic is an intermediary calculation of the Hessian in the iterations.  It is not the final one.  The covariance matrix returned in cov is computed after convergence and will differ from the ones computed during the iterations.

If you are interested in the Hessian computed after convergence, you will find that in the global, _max_FinalHess.

The covariance matrix in the cov return argument takes into account the number of observations.

cov = {
 0.057190334  0.030671550 -0.018652539 -0.010151653,
 0.030671550  0.057763764 -0.010134082 -0.019056447,
-0.018652539 -0.010134082  0.0081342988 0.0043459482,
-0.010151653 -0.019056447  0.0043459482 0.0083873686 };

 
Here is invpd(cov)/1000 and you'll notice that it closely resembles the Hessian retrieved in _max_diagnostic.

 0.0908 -0.0441  0.2049 -0.0964
-0.0441  0.0906 -0.0964  0.2024
 0.2049 -0.0964  0.6321 -0.2985
-0.0964  0.2024 -0.2985  0.6170

Here is 1000*cov and you'll notice that it resembles the inverse of the Hessian in _max_diagnostic.

 57.1903  30.6715 -18.6525 -10.1517
 30.6715  57.7638 -10.1341 -19.0564
-18.6525 -10.1341   8.1343   4.3459
-10.1517 -19.0564   4.3459   8.3874

Thank you for the response. However, still there is some discrepancy between  the inverse of hessian and 1000*cov :

This is inverse of the Hessian

45.059278  24.224743 -15.109826 -8.5460929
24.224743  45.880519 -8.6733621 -15.237785
-15.109826 -8.6733621 7.0742131  4.0424484
-8.5460929 -15.237785 4.0424484  7.0617510

and  this is 1000*cov :

 57.1903  30.6715 -18.6525 -10.1517
 30.6715  57.7638 -10.1341 -19.0564
-18.6525 -10.1341   8.1343   4.3459
-10.1517 -19.0564   4.3459   8.3874

Is there anything else that is factored in calculating cov?

Thank you very much. Now it makes sense.