La estimación con el MLE y devolver la puntuación/gradiente (QMLE)

Question

Estoy en la estimación de una simple AR(1) el proceso por el ML enfoque. Yo también deseo para el cálculo de la Cuasi MLE de los errores estándar, el cual es dado por el sándwich forma de Hesse y la Puntuación (véase, por ejemplo, la última diapositiva aquí)

Así que, yo simplemente especificando el (condicional) registro de probabilidad (gauss) AR(1) el proceso. Entonces me optimizar este con R optim, que devuelve el estado de Hesse para mí, evaluados en el MLE-estima, que utilizo como mi información de la matriz de estimación, para obtener los errores estándar de los parámetros.

Tan lejos y tan bien (me da los mismos resultados que con las estadísticas de la caja de herramientas de Matlab).

Pero, ¿cómo debo proceder para estimar la QMLE errores estándar? Para que necesito la estimación del exterior producto de la puntuación de la función (es decir, el exterior producto de la gradiente evaluado en el MLE de las estimaciones).

No he encontrado ninguna manera de obtener una estimación (numéricamente) para el degradado en cualquiera de R optimización /ML comandos. Me estoy perdiendo algo?. Gracias

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.

EDITAR: Yo tal vez podría probar y utilizar el numDeriv paquete para obtener el gradiente de la función de probabilidad (evaluada en cada una de las observaciones). Pero estoy atascado en exactamente cómo lograr mi objetivo, como no sé cómo reescribir mi probabilidad de la función para ese propósito...

EDIT2: NA

EDIT3: lo Siento por mi estupidez, la suma de exterior de productos es, por supuesto, no es el mismo que el exterior producto de las sumas. Parece coherente ahora:

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

Donde LLi es la Probabilidad de la función de cada observación:

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

Lo que me da errores estándar:

> se2
[1] 0.41208510 0.04256279 0.10242072

Que es razonablemente idénticos(?) a los obtenidos por el estado de Hesse:

> se
[1] 0.40621637 0.04179929 0.09874189

Alguna sugerencia para mejoras? Programación de mi enfoque, no parece que elegante. Gracias de nuevo.

Answer 1

El numDeriv paquete hecho puede ser utilizado para el cálculo del gradiente y el de hesse (si es necesario). En ambos casos, el argumento y de la log-verosimilitud se pasa a través del mecanismo de puntos, usando un argumento con el nombre adecuado. Para un vector de valores de la función, el jacobian función en el mismo paquete puede ser utilizado de manera similar.

También podría considerar la computación analíticos derivados en lugar de valores numéricos.

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

Podemos calcular así el Jacobiano de una función que devuelve un vector de longitud $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)