¿Por qué son los Procesos Gaussianos modelos estadísticos válidos para el pronóstico de series temporales?

Question

Renuncia de duplicados: Sé sobre la pregunta

Predicción de series temporales usando regresión de Procesos Gaussianos

pero este no es un duplicado, porque esa pregunta solo se preocupa por modificaciones en la función de covarianza, mientras que yo argumento que en realidad también se debe modificar el término de ruido.

Pregunta: en un modelo de series temporales, la suposición usual de regresión no paramétrica (datos iid) falla, porque los residuos están autocorrelacionados. Ahora, la Regresión de Procesos Gaussianos es una forma de regresión no paramétrica, porque el modelo subyacente es, asumiendo que tenemos una muestra iid $D=\{(t_i,y_i)\}_{i=1}^N$:

$$y_i = \mathcal{GP}(t_i)+\epsilon_i,\ i=1,\dots,N$$

donde $\mathcal{GP}(t)$ es un Proceso Gaussiano, con una función de media $mu(t)$ (usualmente asumida como 0) y una función de covarianza $k(t,t')$, mientras que $\epsilon\sim\mathcal{N}(0,\sigma)$. Posteriormente utilizamos inferencia bayesiana para calcular la distribución predictiva posterior, dada la muestra $D$.

Sin embargo, en datos de series temporales, la muestra no es iid. Por lo tanto:

1. ¿cómo justificamos el uso de tal modelo?
2. dado que la función de media para la predicción de series temporales con Procesos Gaussianos usualmente se asume como cero, cuando calculo una predicción lo suficientemente lejos en el futuro, espero que regrese a la media de los datos. Esto parece una elección particularmente pobre, porque me gustaría poder (en principio) predecir tan lejos en el futuro como quiera, y que el modelo logre acertar con la tendencia general del tiempo, con solo un aumento en la incertidumbre de la predicción (ver el caso a continuación con un modelo ARIMA):

¿cómo se maneja esto al utilizar Procesos Gaussianos para la predicción de series temporales?

Answer 1

Algunos conceptos relevantes pueden surgir en la pregunta ¿Por qué incluir latitud y longitud en un GAM corrige la autocorrelación espacial?

Si utilizas Procesamiento gaussiano en regresión entonces incluyes la tendencia en la definición del modelo $y(t) = f(t,\theta) + \epsilon(t)$ donde esos errores son $\epsilon(t) \sim \mathcal{N}(0,{\Sigma})$ con $\Sigma$ dependiendo de alguna función de la distancia entre los puntos.

En el caso de tus datos, niveles de CO2, puede ser que el componente periódico sea más sistemático que solo ruido con una correlación periódica, lo que significa que es mejor incorporarlo en el modelo

Demostración utilizando el modelo DiceKriging en R.

La primera imagen muestra un ajuste de la línea de tendencia $y(t) = \beta_0 + \beta_1 t + \beta_2 t^2 +\beta_3 \sin(2 \pi t) + \beta_4 \sin(2 \pi t)$.

El intervalo de confianza del 95% es mucho más pequeño en comparación con la imagen arima. Pero nota que el término residual también es muy pequeño y hay muchos puntos de datos. Para comparación se realizan otros tres ajustes.

• Se ajusta un modelo más simple (lineal) con menos puntos de datos. Aquí puedes ver el efecto del error en la línea de tendencia que causa que el intervalo de confianza de la predicción aumente al extrapolar más lejos (este intervalo de confianza también es tan correcto como lo sea el modelo).
• Un modelo de mínimos cuadrados ordinarios. Se puede ver que proporciona más o menos el mismo intervalo de confianza que el modelo de proceso gaussiano
• Un modelo de Kriging ordinario. Este es un proceso gaussiano sin la tendencia incluida. Los valores predichos serán iguales a la media cuando se extrapolan lejos. La estimación del error es grande porque los términos residuales (dato-media) son grandes.

library(DiceKriging)
library(datasets)

# data
y <- as.numeric(co2)
x <- c(1:length(y))/12

# design-matrix 
# the model is a linear sum of x, x^2, sin(2*pi*x), and cos(2*pi*x)
xm <- cbind(rep(1,length(x)),x, x^2, sin(2*pi*x), cos(2*pi*x))
colnames(xm) <- c("i","x","x2","sin","cos")

# fitting non-stationary Gaussian processes 
epsilon <- 10^-3
fit1 <- km(formula= ~x+x2+sin+cos, 
          design = as.data.frame(xm[,-1]), 
          response = as.data.frame(y),
          covtype="matern3_2", nugget=epsilon)

# fitting simpler model and with less data (5 years)
epsilon <- 10^-3
fit2 <- km(formula= ~x, 
           design = data.frame(x=x[120:180]), 
           response = data.frame(y=y[120:180]),
           covtype="matern3_2", nugget=epsilon)

# fitting OLS
fit3 <- lm(y~1+x+x2+sin+cos, data = as.data.frame(cbind(y,xm)))

# ordinary kriging 
epsilon <- 10^-3
fit4 <- km(formula= ~1, 
           design = data.frame(x=x), 
           response = data.frame(y=y),
           covtype="matern3_2", nugget=epsilon)

# predictions and errors
newx <- seq(0,80,1/12/4)
newxm <- cbind(rep(1,length(newx)),newx, newx^2, sin(2*pi*newx), cos(2*pi*newx))
colnames(newxm) <- c("i","x","x2","sin","cos")
# using the type="UK" 'universal kriging' in the predict function
# makes the prediction for the SE take into account the variance of model parameter estimates
newy1 <- predict(fit1, type="UK", newdata = as.data.frame(newxm[,-1]))
newy2 <- predict(fit2, type="UK", newdata = data.frame(x=newx))
newy3 <- predict(fit3, interval = "confidence", newdata = as.data.frame(x=newxm))
newy4 <- predict(fit4, type="UK", newdata = data.frame(x=newx))

# plotting
plot(1959-1/24+newx, newy1$mean,
 col = 1, type = "l",
 xlim = c(1959, 2039), ylim=c(300, 480),
 xlab = "time [years]", ylab = "atmospheric CO2 [ppm]")
polygon(c(rev(1959-1/24+newx), 1959-1/24+newx), c(rev(newy1$lower95), newy1$upper95), 
        col = rgb(0,0,0,0.3), border = NA)
points(1959-1/24+x, y, pch=21, cex=0.3, col=1, bg="white")
title("Proceso gaussiano con función polinómica + trigonométrica para la tendencia")

# plotting
plot(1959-1/24+newx, newy2$mean,
 col = 2, type = "l",
 xlim = c(1959, 2010), ylim=c(300, 380),
 xlab = "time [years]", ylab = "atmospheric CO2 [ppm]")
polygon(c(rev(1959-1/24+newx), 1959-1/24+newx), c(rev(newy2$lower95), newy2$upper95), 
        col = rgb(1,0,0,0.3), border = NA)
points(1959-1/24+x, y, pch=21, cex=0.5, col=1, bg="white")
points(1959-1/24+x[120:180], y[120:180], pch=21, cex=0.5, col=1, bg=2)
title("Proceso gaussiano con función lineal para la tendencia")

# plotting
plot(1959-1/24+newx, newy3[,1],
     col = 1, type = "l",
     xlim = c(1959, 2039), ylim=c(300, 480),
     xlab = "time [years]", ylab = "atmospheric CO2 [ppm]")
polygon(c(rev(1959-1/24+newx), 1959-1/24+newx), c(rev(newy3[,2]), newy3[,3]), 
        col = rgb(0,0,0,0.3), border = NA)
points(1959-1/24+x, y, pch=21, cex=0.3, col=1, bg="white")
title("Regresión lineal ordinaria con función polinómica + trigonométrica para la tendencia")

# plotting
plot(1959-1/24+newx, newy4$mean,
 col = 1, type = "l",
 xlim = c(1959, 2039), ylim=c(300, 480),
 xlab = "time [years]", ylab = "atmospheric CO2 [ppm]")
polygon(c(rev(1959-1/24+newx), 1959-1/24+newx), c(rev(newy4$lower95), newy4$upper95), 
        col = rgb(0,0,0,0.3), border = NA, lwd=0.01)
points(1959-1/24+x, y, pch=21, cex=0.5, col=1, bg="white")
title("kriging ordinario")

Answer 2

Una de las principales suposiciones para el GP es que los datos deben ser estacionarios. Tus datos tienen una clara tendencia, por lo tanto no son estacionarios. La forma correcta de utilizar el GP en series temporales (y en cualquier otro tipo de datos) es primero eliminar algunas tendencias obvias, luego aplicar el GP sobre el residual.