# Standard errors for Composite Marginal Likelihood

I am estimating a multivariate ordered probit model using a composite marginal likelihood (CML) approach. In other words, I replace the full likelihood function by a surrogate likelihood constructed from pair-wise bivariate ordered probit models. The code is written in MATLAB and I use the solver fmincon to maximize my likelihood under constraints.
The code works perfectly and I recover the true parameters in my empirical exercise but I am struggling now to compute the standard errors.

My CML estimator is asymptotically normally distributed. The asymptotic variance-covariance matrix for the estimator is given by the inverse of the Godambe sandwich information matrix (Godambe, 1960) defined by:

\begin{align*} & V_{CML}(\theta)= [H(\theta)]^{-1} J(\theta) [H(\theta)]^{-1} \quad \text{with} \\ & H(\theta) = E \Big[ -\frac{\partial^2 \: \: log \: L_{CML}(\theta) }{\partial \theta \: \: \partial \theta'} \Big] \\ & J(\theta) = E \Big[ \left( \frac{\partial \: \: log \: L_{CML}(\theta) }{\partial \theta } \right)^{'} \left( \frac{\partial \: \: log \: L_{CML}(\theta) }{\partial \theta } \right) \Big] \end{align*}

I tried to use H as the Hessian at the estimated parameter values and J as the Gradient vector. Both are given by fmincon. When I then take the square root of the diagonal elements, I find very low standard errors which does not make sense to me.
Note that I am simulating the data.

Am I doing something wrong?

Thank you

Alexis