How to compute the sandwich variance ML estimator in R

Question

I'm currently estimating a DCC-type model by maximum likelihood. Im using the command solnp and it return an object where I can compute the Hessian H evaluated at the optimal values.

If I want to find the covariance matrix of the estimators, following Davidson and MacKinnon (2004), I can compute $$ \hat{Var} ( \hat{\theta} ) = -H^{-1} ( \hat{\theta} ) $$ However, there is a robust estimation for distribution miss-especification: the sandwich estimator or the quasi-maximum likelihood estimator (QMLE). It is as follows: $$ \hat{Var} ( \hat{\theta} ) = H^{-1} ( \hat{\theta} ) G^{T} ( \hat{\theta} ) G ( \hat{\theta} ) H^{-1} ( \hat{\theta} ) $$ where G is the gradient.

As solnp does not compute the gradient, I am using the command gHgenb to compute it at the optimal values of the parameters.

The problem is that when I compute the non-robust covariance matrix I get the following results:

                     Mean           SE          t        Pvalue
[Joint]dcca1 0.0087397376 9.267334e-03  0.9430692  3.461028e-01
[Joint]dccb1 0.9653333580 1.255250e-02 76.9036529 1.723509e-278
             0.0001513922 4.115874e-05  3.6782522  2.604074e-04

And when I compute the sandwich estimators I get the following results:

                     Mean           SE          t Pvalue
[Joint]dcca1 0.0087397376 6.622702e-09  1319663.4      0
[Joint]dccb1 0.9653333580 6.706874e-08 14393193.4      0
             0.0001513922 1.538495e-10   984028.1      0

where the standard errors are extremely low.

The questions is: what should I do? Should I multiplicate any of the matrix by the number of observations (which one?)?, as Bollerslev and Wooldridge (1992) state that the non-robust matrix is: $$ \hat{Var} ( \hat{\theta} ) = -H^{-1} ( \hat{\theta} ) / T $$

Thank you so much.

References:

Bollerslev, T., & Wooldridge, J. M. (1992). Quasi-maximum likelihood estimation and inference in dynamic models with time-varying covariances. Econometric reviews, 11(2), 143-172.

Davidson, R., & MacKinnon, J. G. (2004). Econometric theory and methods (Vol. 5). New York: Oxford University Press.

Answer 1

I just found the answer to my problem by comparing the manual method versus the QMLE robust standard errors for a GARCH using the package fGarch. The solution is as follows: Let hessian be the hessian computed by the optimization at the optimal values. Now, I computed a $N\times K$ matrix (where $N$ is the number of observations and $K$ the number of parameters). Using the command grad, from the package numDeriv, I computed the gradient for every observation in the optimal parameters and with those values I filled the matrix $N\times K$, which I called hessian. Then, the Information matrix and the QMLE robust standard errors are computed as follows:

meat<-t(gradient)%*%gradient
Info<-solve(hessian)%*%meat%*%solve(hessian)
SE_robust <- sqrt(diag(Info))

Answer 2

The sandwich matrix is $-H^{-1}VH^{-1}$ for $V$ an estimate of the expectation of $G(\theta_0)^TG(\theta_0)$, where $G$ is a row vector (following the poster's notation). Using $G(\hat{\theta})^TG(\hat{\theta})$ to approximate this expectation will not provide a good answer, because $G(\hat{\theta}) = 0$ by definition of $\hat{\theta}$!

To provide slightly more detail: if $\hat{\theta}$ solves $\sum_{i=1}^{n}\psi_i(\theta) = 0$ where $\psi_i$ are independent functions of the response $Y_i$, and $\theta_0$ solves $\Psi(\theta) = 0$ where $\Psi = E(\psi_1)$, then $\sqrt{n}(\hat{\theta} - \theta_0)$ is asymptotically $N(0,V)$ where $V = H^{-1}SH^{-1}$. These matrices are $H = E\nabla_\theta\psi_1(\theta_0)$ and $S = E\psi_1(\theta_0)\psi_1(\theta_0)^T$. We approximate $H$ by $\tilde{H} = \sum_{i=1}^{n}\psi^{\prime}_{i}(\hat{\theta})$ and $V$ by $\tilde{V} = \sum_{i=1}^{n}\psi_i(\hat{\theta})\psi_i(\hat{\theta})^T$. This last step is the part that's different from the question: $\sum_{i=1}^{n}\psi_i(\hat{\theta}) = 0$ so $\left(\sum_{i=1}^{n}\psi_i(\hat{\theta})\right)\left(\sum_{i=1}^{n}\psi_i(\hat{\theta})\right)^T = 0$, but $\sum_{i=1}^{n}\psi_i(\hat{\theta})\psi_i(\hat{\theta})^T \neq 0$.

Here is a useful reference: https://www.stat.berkeley.edu/~census/mlesan.pdf

My explanation is a poor recreation of material from van der Vaart's Asymptotic Statistics, Section 5.3, page 51/52; that's where to look for the good stuff!