# How to compute the sandwich variance ML estimator in R

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.

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))


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!