We consider the following GARCH(1, 1) model:
$y_t = h_t \epsilon_t$ where $(\epsilon_t)_{t \in \{1, \dots, n\}}$ are i.i.d. random variables with mean 0 and standard deviation 1.
$h_t = \omega + \alpha \epsilon_{t-1}^2 + \beta h_{t-1}^2$
A standard way of estimating the parameters of such a model is to use the quasi maximum likelihood method, i.e. maximising the log-likelihood of the model under the assumption that the $(\epsilon_t)_{t \in \{1, \dots, n\}}$ are gaussian, that is:
$\displaystyle\mathscr{l}(\omega, \alpha, \beta, y_1, \dots, y_n) = -\frac{1}{2}\log(2\pi) - \frac{1}{2}\sum_{t=1}^n\left(\log(h_t) + \frac{y_t^2}{h_t}\right)$
The estimator is shown to have good properties even if this assumption does no hold. The asymptotic distribution of that estimator is:
$\sqrt{n}\left(\hat{\theta}_n - \theta_0\right) \rightarrow \mathcal{N}\left(0, A(\theta_0)^{-1} B(\theta_0) A(\theta_0)^{-1}\right)$
where $\hat{\theta}_n$ is the vector of estimated parameters, $\theta_0$ the true vector of parameters, $B(\theta_0)$ the outer product of the gradient and $A(\theta_0)$ the Hessian matrix, both evaluated at the true value of the parameters.
We can get a good approximation of the Hessian asymptotically by replacing the true value of the parameter by its estimate. But when it comes to the gradient, $B(\hat{\theta}_n)$ is by definition equal to zero since $\hat{\theta}_n$ maximises $\mathscr{l}$. So it cannot be a good approximation of $B(\theta_0)$.
That leads me to my question: how do you estimate the covariance matrix of the asymptotic distribution of the parameters?