I am estimating a simple AR(1) process by the ML approach. I also wish to compute the Quasi MLE standard errors, which is given by the sandwich form of the Hessian and the Score (see for example the last slide here)

So, I start by just specifying the (conditional) log likelihood for the (gaussian) AR(1) process. Then I optimise this with R's optim, which returns the Hessian to me, evaluated at the MLE-estimates, which I use as my information matrix estimate, to get the standard errors of my parameters.

So far so good (I get the same results as with the stats toolbox in Matlab).

But, how do I proceed to estimate the QMLE standard errors? For that I need the estimate of the outer product of the score function (i.e. the outer product of the gradient evaluated at the MLE estimates).

I have not found any way to get an estimate (numerically) for the gradient in any of R's optimization /ML commands. Am I missing something?. Thank you

data = read.table("Data/AR.txt", header=FALSE)
y = as.vector(data$V1) # A simple vector of observations: n1, n2, ... , nT

#Conditional LH
loglik = function(theta, y) {
  T = length(y)
  L = sum (dnorm(y[2:T], 
       mean = theta[1] + theta[2]*y[1:T-1], 
       sd = theta[3], log = TRUE))
  return (-L)

start=c(2.5, 0.6, 3)
b = optim(start, loglik, y=y, hessian=TRUE)
I = solve(b$hessian)
se = sqrt(diag(I)) # All good. The same MLE estimates and SE's as in Matlab.

EDIT: I could perhaps try and use the numDeriv package to get the gradient of the likelihood function (evaluated at every observation). But I am stuck on exactly how to accomplish my goal, as I don't know how to rewrite my likelihood function for that purpose...


EDIT3: Sorry for my stupidity, the sum of outer products is of course not the same as the outer product of the sums. It seems consistent now:

sum = numeric(3)
for (t in 2:length(y)) {
  g = grad(LLi, x=b$par, y=y, t=t)
  sum = sum + outer(g,g)

I2 = solve(sum)
se2 = sqrt(diag(I2))

Where LLi is the Likelihood function for each observation:

LLi = function(theta, y, t) {
  L = dnorm(y[t], 
                 mean = theta[1] + theta[2]*y[t-1], 
                 sd = theta[3], log = TRUE)
  return (L)

Which gives me standard errors:

> se2
[1] 0.41208510 0.04256279 0.10242072

Which is reasonably identical(?) to those obtained by the Hessian:

> se
[1] 0.40621637 0.04179929 0.09874189

Any suggestions for improvements? Programming wise my approach doesnt seem that elegant. Thanks again.

  • $\begingroup$ Maybe you could use simulated data rather that file data to make your example reproducible? E.g. y <- arima.sim(n = 100, list(ar = 0.9)) $\endgroup$
    – Yves
    Apr 29, 2015 at 9:50
  • $\begingroup$ I think that you could clean up your question by removing edit #2. It could be simpler and more efficient to use jacobian on a, say, 'logDens' function (with no 't' argument) returning the vector of the $T-1$ log-density values.The jacobian will be a $(T-1) \times 3$ matrix on which an 'apply' can be used. The analytical values of the derivatives could also be computed. $\endgroup$
    – Yves
    Apr 29, 2015 at 12:43
  • $\begingroup$ I really don't want to use analytical derivatives since I want to be able to use this approach on more difficult likelihoods later (I just went trough this simple example to make sure I got it working). I don't understand how I could use the Jacobian, could you please give an example? Thanks $\endgroup$
    – luffe
    Apr 29, 2015 at 12:49
  • $\begingroup$ I have made a 'logDens' function and gotten the Jacobian (T-1)x3 matrix, but what 'apply' function to I use on this? $\endgroup$
    – luffe
    Apr 29, 2015 at 13:14
  • 1
    $\begingroup$ I added the Jacobian example. Note that 1:T-1 is not the same thing as the 1:(T-1) that you expect, so parentheses are needed here. An common trap with R indices. $\endgroup$
    – Yves
    Apr 29, 2015 at 17:37

1 Answer 1


The numDeriv package can indeed be used to compute the gradient and the hessian (if needed). In both cases the argument y of the log-likelihood is passed through the dots mechanism, using an argument with the suitable name. For a vector-valued function, the jacobian function of the same package can be used similarly.

You could also consider computing analytical derivatives rather than numerical ones.

H <- hessian(func = loglik, x = b$par, y = y)
g <- grad(func = loglik, x = b$par, y = y)

We can compute as well the Jacobian of a function returning a vector of length $T-1$.

mlogDens <- function(theta, y) {
  T <- length(y)
  -dnorm(y[2:T], mean = theta[1] + theta[2] * y[1:(T-1)], 
                 sd = theta[3], log = TRUE)
## a matrix with dim (T-1, 3)
G <- jacobian(func = mlogDens, x = b$par, y = y)
## a matrix with dim (9, T-1)
GG <- apply(G, MARGIN = 1, FUN = tcrossprod)
## a vector with length 9 representing a symmetric mat.
GGsum0 <- apply(GG, MARGIN = 1, FUN = sum)
## a symmetric mat.
GGsum <- matrix(GGsum0, nrow = 3, ncol = 3)
  • 2
    $\begingroup$ Instead of doing all those apply's, could I just not do GGsum = t(G) %*% Gor is that not preferred? Thanks again, $\endgroup$
    – luffe
    Apr 30, 2015 at 10:40
  • 1
    $\begingroup$ Oh yes, this is simpler and may even be better. $\endgroup$
    – Yves
    Apr 30, 2015 at 11:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.