Errores no correlacionados del modelo de mínimos cuadrados generalizados (GLS)

Question

Como institución financiera, a menudo nos encontramos con el análisis de datos de series temporales. Muchas veces acabamos haciendo una regresión utilizando variables de series temporales. Cuando esto ocurre, a menudo nos encontramos con residuos con estructura de serie temporal que viola el supuesto básico de errores independientes en la regresión 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) tiene claramente una estructura AR(1), como se desprende de ACF y PACF. Adopté dos enfoques diferentes, el primero fue, obviamente, ajustar el modelo mediante mínimos cuadrados generalizados gls() en R. Mi expectativa era que los residuos de gls(object) fueran un ruido blanco (errores independientes). Pero los residuos de gls(object) siguen teniendo la misma estructura ARIMA que en la regresión ordinaria. Lamentablemente, hay algo que no funciona en lo que estoy haciendo y que no he podido averiguar. Por lo tanto, decidí ajustar manualmente los coeficientes de regresión del modelo lineal (estimaciones OLS). Sorprendentemente parece que funciona cuando trazo los residuos de la regresión ajustada (los residuos son ruido blanco). Realmente quiero utilizar gls() en nlme para que la codificación sea mucho más sencilla y fácil. ¿Cuál sería el enfoque que debería adoptar aquí? ¿Debo usar REML? o ¿es mi expectativa de residuos no correlacionados (ruido blanco) del objeto gls() errónea?

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 algún problema con lo que estoy haciendo? ???????

Answer 1

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

En resumen, no has hecho nada mal, no hay necesidad de ajustar los residuos y las rutinas funcionan correctamente.

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

EDIT: A raíz de un comentario de Michael Chernick (+1), añado un ejemplo en el que se comparan los resultados de GLS con los de ARMAX (arima), mostrando que las estimaciones de los coeficientes, las probabilidades logarítmicas, etc. son todas iguales, al menos con cuatro decimales (un grado de precisión razonable 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

EDIT: A raíz de un comentario de anand (OP), he aquí una comparación de las predicciones de gls y arima con la misma estructura de datos básica que la anterior y eliminando algunas líneas de salida extrañas:

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 previstos son diferentes, aunque convergen a medida que nos alejamos en el futuro. Esto se debe a que gls no trata los datos como una serie temporal y tiene en cuenta el valor específico del residuo en la observación 995 a la hora de formar predicciones, sino que arima sí. El efecto del residuo en la obs. 995 disminuye a medida que aumenta el horizonte de previsión, lo que conduce a la convergencia de los valores previstos.

En consecuencia, para las predicciones a corto plazo de los datos de las series temporales, arima será mejor.

Answer 2

NOTA: A ver si ahora puedo publicar esta respuesta :)

---

En primer lugar, permítanme señalar que la condición de que la norma de Frobenius sea uno es irrelevante, ya que al multiplicar la matriz por alguna constante positiva se modifican tanto las entradas como los valores singulares por la misma constante multiplicativa. De todos modos, por la misma razón, podemos limitarnos al caso mencionado en la pregunta.

La solución se ve más fácilmente en el lenguaje de la información cuántica. Las entradas $A_{ij}$ de la matriz pueden verse como los coeficientes de un vector de estado bipartito $|\psi\rangle=\sum_{ij}A_{ij}|i\rangle|j\rangle$ , donde $\{|i\rangle\}$ y $\{|j\rangle\}$ son bases ortonormales (locales). Supongamos que el mayor $k$ entradas $|A_{ij}|^2$ corresponden a pares de índices $(i,j)\in I_k$ , donde $I_k$ es algo $k$ -subconjunto de elementos del conjunto $\{ (i,j)|i,j=1,\ldots,n \}$ . Entonces

$$ \sum_{(i,j)\in I_k}|A_{ij}|^2 = \sum_{(i,j)\in I_k}\textrm{Tr}(|i\rangle\langle i|\otimes |j\rangle\langle j| |\psi\rangle\langle\psi|) $$

Esta cantidad es inferior a

$$ \sum_{i:\exists j s.t. (i,j)\in I_k}\textrm{Tr}(|i\rangle\langle i|\otimes I |\psi\rangle\langle\psi|)=\textrm{Tr}_1(\sum_{i:\exists j s.t. (i,j)\in I_k}|i\rangle\langle i|\textrm{Tr}_2(|\psi\rangle\langle\psi|)) $$

donde $I$ denota la matriz/operador de identidad y hemos dividido la traza en la traza sobre el primer sistema y sobre el segundo sistema. La matriz $\textrm{Tr}_2(|\psi\rangle\langle\psi|)$ es semidefinido positivo, con valores propios correspondientes a los cuadrados de los valores singulares de $A$ (es el estado reducido de $|\psi\rangle$ en el lenguaje de la información cuántica), y tenemos

$$ \textrm{Tr}_1(\sum_{i:\exists j s.t. (i,j)\in I_k}|i\rangle\langle i|\textrm{Tr}_2(|\psi\rangle\langle\psi|))\leq \max_{P_k} \textrm{Tr}_1(P_k\textrm{Tr}_2(|\psi\rangle\langle\psi|))=\sum_{i=1}^k\sigma_k^2, $$

donde el máximo es sobre las proyecciones de rango $k$ . Como esto es válido para cualquier $k$ la conjetura queda demostrada.