Errores no correlacionados del modelo de mínimo cuadrado generalizado (GLS)

Como institución financiera, a menudo nos encontramos con el análisis de datos de series temporales. Muchas veces terminamos haciendo regresión usando variables de series de tiempo. A medida que esto sucede, a menudo encontramos residuos con una estructura de series de tiempo que viola la suposición básica de errores independientes en la regresión de OLS. Recientemente, estamos construyendo otro modelo en el que creo que tenemos una regresión con errores autocorrelacionados. Los residuos del modelo lineal tienen lm(object)claramente una estructura AR (1), como se desprende de ACF y PACF. Tomé dos enfoques diferentes, el primero obviamente era ajustar el modelo usando mínimos cuadrados generalizados gls()en R. Esperaba que los residuos de gls (objeto) fueran un ruido blanco (errores independientes). Pero los residuos degls(object)todavía tienen la misma estructura ARIMA que en la regresión ordinaria. Lamentablemente, hay algo mal en lo que estoy haciendo que no pude resolver. Por lo tanto, decidí ajustar manualmente los coeficientes de regresión del modelo lineal (estimaciones OLS). Sorprendentemente, eso parece estar funcionando cuando tracé los residuos de la regresión ajustada (los residuos son ruido blanco). Realmente quiero usar gls()en el nlmepaquete para que la codificación sea mucho más simple y fácil. ¿Cuál sería el enfoque que debería tomar aquí? ¿Se supone que debo usar REML? ¿o es incorrecto mi expectativa de residuos no correlacionados (ruido blanco) del objeto 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

¿Hay algo mal con lo que estoy haciendo ???????

De hecho, los residuos de glstendrán la misma estructura de autocorrelación, pero eso no significa que las estimaciones de coeficientes y sus errores estándar no se hayan ajustado adecuadamente. (Obviamente, tampoco es necesario que sea ​​diagonal.) Esto se debe a que los residuos se definen como . Si la matriz de covarianza de fuera igual a , ¡no habría necesidad de usar GLS!$\Omega$$e = Y - X\hat{\beta}^{\text{GLS}}$$e$$\sigma^2\text{I}$

En resumen, no ha hecho nada malo, no hay necesidad de ajustar los residuos, y todas las rutinas funcionan correctamente.

predict.glstoma en cuenta la estructura de la matriz de covarianza al formar errores estándar del vector de predicción. Sin embargo, no tiene la característica conveniente de "predecir algunas observaciones por delante" predict.Arima, que tiene en cuenta los residuos relevantes al final de la serie de datos y la estructura de los residuos al generar los valores pronosticados. arimatiene la capacidad de incorporar una matriz de predictores en la estimación, y si está interesado en la predicción unos pasos más adelante, puede ser una mejor opción.

EDITAR: provocado por un comentario de Michael Chernick (+1), estoy agregando un ejemplo que compara GLS con los resultados de ARMAX (arima), que muestra que las estimaciones de coeficientes, las probabilidades de registro, etc., todos salen iguales, al menos con cuatro decimales lugares (un grado razonable de precisión dado que se utilizan dos algoritmos diferentes):

# 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

EDITAR: provocado por un comentario de anand (OP), aquí hay una comparación de predicciones de glsy arimacon la misma estructura de datos básica que anteriormente y algunas líneas de salida extrañas eliminadas:

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

Como podemos ver, los valores pronosticados son diferentes, aunque convergen a medida que avanzamos hacia el futuro. Esto se debe a glsque no trata los datos como una serie temporal y no tiene en cuenta el valor específico del residuo en la observación 995 al formar predicciones, pero lo arimahace. El efecto del residuo en obs. 995 disminuye a medida que aumenta el horizonte de pronóstico, lo que conduce a la convergencia de los valores pronosticados.

En consecuencia, para predicciones a corto plazo de datos de series de tiempo, arimaserá mejor.

Desea los residuos normalizados. Ver ?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")