我正在通过 ML 方法估计一个简单的 AR(1) 过程。我还希望计算准 MLE 标准误差，它由 Hessian 和分数的三明治形式给出（例如，请参见此处的最后一张幻灯片）

因此，我首先指定（高斯）AR（1）过程的（条件）对数似然。然后我用 R 的optim优化这个，它返回 Hessian 给我，在 MLE 估计中进行评估，我用它作为我的信息矩阵估计，以获得我的参数的标准误差。

到目前为止一切顺利（我得到的结果与 Matlab 中的统计工具箱相同）。

但是，如何继续估计 QMLE 标准误差？为此，我需要对得分函数的外积进行估计（即在 MLE 估计中评估的梯度的外积）。

在任何 R 的优化 /ML 命令中，我还没有找到任何方法来获得梯度的估计（数字）。我错过了什么吗？谢谢

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.

编辑：我也许可以尝试使用numDeriv包来获得似然函数的梯度（在每次观察时评估）。但是我被困在如何实现我的目标上，因为我不知道如何为此目的重写我的似然函数......

EDIT2：不适用

EDIT3：对不起，我的愚蠢，外积的总和当然与总和的外积不同。现在看来一致：

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

其中LLi是每个观察的似然函数：

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

这给了我标准错误：

> se2
[1] 0.41208510 0.04256279 0.10242072

哪个与 Hessian 获得的相当相同（？）：

> se
[1] 0.40621637 0.04179929 0.09874189

有什么改进建议吗？明智地编程我的方法似乎并不优雅。再次感谢。