广义最小二乘模型 (GLS) 的非相关误差

机器算法验证 r 回归 时间序列 广义最小二乘法
2022-03-20 16:49:02

作为金融机构,我们经常会遇到时间序列数据的分析。很多时候,我们最终使用时间序列变量进行回归。发生这种情况时,我们经常会遇到时间序列结构的残差,这违反了 OLS 回归中独立误差的基本假设。最近我们正在构建另一个模型,我相信我们在其中有自相关误差的回归。线性模型的残差lm(object)显然具有 AR(1) 结构,从 ACF 和 PACF 可以看出。我采用了两种不同的方法,第一种显然是使用 R 中的广义最小二乘法拟合模型gls()。我的期望是 gls(object) 的残差将是白噪声(独立误差)。但残差来自gls(object)仍然具有与普通回归相同的 ARIMA 结构。不幸的是,我正在做的事情有问题,我无法弄清楚。因此,我决定手动调整线性模型(OLS 估计)的回归系数。令人惊讶的是,当我绘制调整回归的残差(残差是白噪声)时,这似乎起作用了。我真的很想gls()nlme包中使用,这样编码会更简单、更容易。我应该在这里采取什么方法?我应该使用REML吗?还是我对 gls() 对象的非相关残差(白噪声)的期望是错误的?

gls.bk_ai <- gls(PRNP_BK_actINV ~ PRM_BK_INV_ENDING + NPRM_BK_INV_ENDING, 
                 correlation=corARMA(p=1), method='ML',  data  = fit.cap01A)

gls2.bk_ai  <- update(gls.bk_ai, correlation = corARMA(p=2))

gls3.bk_ai <- update(gls.bk_ai, correlation = corARMA(p=3))

gls0.bk_ai <- update(gls.bk_ai, correlation = NULL)

anova(gls.bk_ai, gls2.bk_ai, gls3.bk_ai, gls0.bk_ai)  
     ## looking at the AIC value, gls model with AR(1) will be the best bet

acf2(residuals(gls.bk_ai)) # residuals are not white noise

我在做什么有问题吗??????

2个回答

残差gls确实具有相同的自相关结构,但这并不意味着系数估计值及其标准误差没有适当调整。(显然也没有要求是对角线。)这是因为残差定义为的协方差矩阵等于,则无需使用 GLS!Ωe=YXβ^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) 的评论提示,这里是来自与上述相同基本数据结构的预测的比较,glsarima删除了一些无关的输出行:

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会更好。

你想要标准化的残差。?residuals.lme

#Reproducible code from ?corARMA
fm1Ovar.lme <- lme(follicles ~ sin(2*pi*Time) + cos(2*pi*Time),
                   data = Ovary, random = pdDiag(~sin(2*pi*Time)))
fm5Ovar.lme <- update(fm1Ovar.lme,
                corr = corARMA(p = 1, q = 1))

#raw residuals divided by the corresponding standard errors
acf(residuals(fm5Ovar.lme),type="partial")

#standardized residuals pre-multiplied 
#by the inverse square-root factor of the estimated error correlation matrix
acf(residuals(fm5Ovar.lme,type="normalized"),type="partial")