残差gls确实具有相同的自相关结构,但这并不意味着系数估计值及其标准误差没有适当调整。(显然也没有要求是对角线。)这是因为残差定义为。的协方差矩阵等于,则无需使用 GLS!Ωe=Y−Xβ^GLSeσ2I
简而言之,你没有做错任何事,不需要调整残差,例程都运行正常。
predict.gls在形成预测向量的标准误差时确实考虑了协方差矩阵的结构。但是,它没有方便的“预测一些观察结果”的特性predict.Arima,该特性在生成预测值时考虑了数据序列末尾的相关残差和残差的结构。 arima能够在估计中加入预测变量矩阵,如果您对提前几步进行预测感兴趣,它可能是一个更好的选择。
编辑:在 Michael Chernick (+1) 的评论提示下,我添加了一个示例,将 GLS 与 ARMAX (arima) 结果进行比较,显示系数估计、对数似然等都相同,至少到小数点后四位地点(考虑到使用了两种不同的算法,合理的准确度):
# Generating data
eta <- rnorm(5000)
for (j in 2:5000) eta[j] <- eta[j] + 0.4*eta[j-1]
e <- eta[4001:5000]
x <- rnorm(1000)
y <- x + e
> summary(gls(y~x, correlation=corARMA(p=1), method='ML'))
Generalized least squares fit by maximum likelihood
Model: y ~ x
Data: NULL
AIC BIC logLik
2833.377 2853.008 -1412.688
Correlation Structure: AR(1)
Formula: ~1
Parameter estimate(s):
Phi
0.4229375
Coefficients:
Value Std.Error t-value p-value
(Intercept) -0.0375764 0.05448021 -0.68973 0.4905
x 0.9730496 0.03011741 32.30854 0.0000
Correlation:
(Intr)
x -0.022
Standardized residuals:
Min Q1 Med Q3 Max
-2.97562731 -0.65969048 0.01350339 0.70718362 3.32913451
Residual standard error: 1.096575
Degrees of freedom: 1000 total; 998 residual
>
> arima(y, order=c(1,0,0), xreg=x)
Call:
arima(x = y, order = c(1, 0, 0), xreg = x)
Coefficients:
ar1 intercept x
0.4229 -0.0376 0.9730
s.e. 0.0287 0.0544 0.0301
sigma^2 estimated as 0.9874: log likelihood = -1412.69, aic = 2833.38
编辑:由 anand (OP) 的评论提示,这里是来自与上述相同基本数据结构的预测的比较,gls并arima删除了一些无关的输出行:
df.est <- data.frame(list(y = y[1:995], x=x[1:995]))
df.pred <- data.frame(list(y=NA, x=x[996:1000]))
model.gls <- gls(y~x, correlation=corARMA(p=1), method='ML', data=df.est)
model.armax <- arima(df.est$y, order=c(1,0,0), xreg=df.est$x)
> predict(model.gls, newdata=df.pred)
[1] -0.3451556 -1.5085599 0.8999332 0.1125310 1.0966663
> predict(model.armax, n.ahead=5, newxreg=df.pred$x)$pred
[1] -0.79666213 -1.70825775 0.81159072 0.07344052 1.07935410
正如我们所看到的,预测值是不同的,尽管随着我们向未来移动,它们正在趋同。这是因为gls在形成预测时不会将数据视为时间序列并考虑观察 995 处残差的特定值,但arima会。obs 处残差的影响。995 随着预测范围的增加而减少,从而导致预测值的收敛。
因此,对于时间序列数据的短期预测,arima会更好。