Metodologia de Analisis de Series de Tiempo
Forecast SARIMA.
Elian Gonzalez Mares
18 de abril de 2024
Email
: ElianElian@ciencias.unam.mx🦕
RPubs : Mis publicaciones
Facultad : Facultad de
Ciencias
Descripcion.
En este documento se presenta la metodología para realizar un análisis de una serie de tiempo usando la herramienta R, teniendo en cuenta los pasos a realizar y algunos conceptos limitados a dos modelos; ARIMA y SARIMA. Nuestro principal interes sera saber identificar el modelo de una serie de tiempo y ponerlo a prueba, despues, para tener fundamentado realizar una prediccion con el modelo.
Teoria.
Propiedad Estacionaria (Debilmente Estacionaria).
Definicion.
Una serie de tiempo \(X_t\) se considera débilmente estacionaria o simplemente estacionaria si cumple con las siguientes propiedades:
La media \(E(X_t)\) es constante e igual para todos los períodos, es decir, \(E(X_t) = E(X_{t+k}) = \mu\).
La varianza \(\text{Var}(X_t)\) es constante e igual para todos los períodos, es decir, \(\text{Var}(X_t) = \text{Var}(X_{t+k}) = \sigma^2\).
La covarianza \(\text{Cov}(X_t,X_{t+k})\) entre dos períodos \(t\) y \(t+k\) depende únicamente de la distancia o rezago \(k\) entre ellos, es decir, \(\text{Cov}(X_t,X_{t+k}) = \gamma_k\).
Cuando una serie de tiempo es estacionaria, sus propiedades estadísticas no cambian con el tiempo, lo que simplifica el proceso de modelado y predicción. Los modelos ARIMA, asumen estacionariedad para su aplicación.
Estacionalidad.
Definicion.
Se dice de una serie de tiempo que es estacional cuando está influenciada por factores estacionales de periodo fijo en todo el tiempo; como el dia, mes, trimestre.
La estacionalidad es la repetición de determinadas variaciones en alguna variable cada cierto período, normalmente igual o menor a un año. En períodos más amplios se suele hablar de ciclos, aunque las variaciones cíclicas no son tan frecuentes como las estacionales. Por ejemplo, cada 2 años se celebra la bienal de Venecia, cada 4 años el mundial de fútbol y cada 76 años se observa el cometa Halley. - BBVA ¿Qué es y cómo nos afecta la estacionalidad?
ARIMA.
Definicion.
ARIMA, denominado como Auto-Regresivo Integrado de Media Movil por sus siglas y denotado como ARIMA(p,d,q), es la combinación de los modelos autorregresivo y el de medias móviles, con la particularidad de que debe incluir un proceso de restablecimiento (el cual se denomina integración I) de inestabilidad original presente en una serie de tiempo.
\[\begin{equation*} x_t' = \phi_1 x'_{t-1} + \ldots + \phi_p x'_{t-p} + Z_t - \theta Z_{t-1} - \ldots - \theta_q Z_{t-q} \end{equation*}\]
\[\begin{equation*} AR(p)I(d) = MA(q) \end{equation*}\]
\(B^d X_t = X_{t-d}\) es el operador de retrasi;
\((1 - B)^d X_t = \Delta^d X_t\) es el operador diferencia, denota la serie temporal diferenciada \(d\) veces;
y la \(X'_t\) denota la serie de tiempo inducida a la estabilidad o estacionariedad.
- (AR) Parametro p. Autroregresion p
Auto Regressive Son modelos donde el valor de una variable en un periodo se relaciona con los valores de periodos previos. Se refiere a los lags p de las diferencias de los valores entre las series. Esta parte del modelo captura la relación lineal entre una observación y un número lineal de observaciones retrasadas (lagged observations). La notación AR(p) representa esta parte, donde “p” es el número de retrasos incluidos en el modelo.
Basicamente este parametro indica que tantos retrasos (lags p) tomamos en cuenta determinan el valor actual, es decir, p retrasos, p datos previos se relacionan con el valor actual.
- (I) Parametro d. Diferenciacion d (Random Walk, etc. suma acumulativa de ruidos blancos (0,0,0))
Se refiere al orden del número diferencial que se usan para hacer las series de tiempo estacionarias.
Esta parte del modelo indica el número de diferenciaciones necesarias para volver estacionaria a la serie temporal, es decir, hacer que la media y la varianza sean constantes en el tiempo. La notación I(d) representa esta parte, donde “d” es el número de diferenciaciones.
- (MA) Parametro q. Error de Regresion q
Consideran la posibilidad de una relación entre la variable y los residuales de periodos previos. Se refiere a los lags q de los errores. Tiene sentido el nombre de media movil ya que
Esta parte del modelo representa la relación lineal entre una observación y los errores residuales de predicciones pasadas en un número lineal de períodos de tiempo. La notación MA(q) representa esta parte, donde “q” es el número de términos de errores residuales incluidos en el modelo.
Basicamente es lo analogo al AR(p), solo que en este caso considera a los q errores o residuales (lags) de la regresion para definir una relacion en el periodo actual con los previos, es decir, q periodos previos de residuos tienen relacion con nuestro periodo actual.
Los datos deben ser estacionarios, o si no, deben sufrir una transformacion (Box-Cox, log, raiz, 1/x, etc.): las propiedades de la serie no dependen del momento en que se capturan. Un proceso estacionario tiene una media y una variación que no cambian con el tiempo y no tiene una tendencia (media o pendiente cambiante).
Los datos deben ser univariados: ARIMA trabaja con una sola variable. La regresión automática tiene que ver con la regresión con los valores pasados.
Linealidad: Se asume que la relación entre las observaciones pasadas y presentes es lineal. Esto significa que las relaciones entre las variables en el modelo pueden ser representadas por una combinación lineal de los parámetros del modelo.
Independencia de los Errores: Se asume que los errores de predicción son independientes entre sí y no muestran autocorrelación significativa. En otras palabras, los residuos del modelo ARIMA deben comportarse como ruido blanco.
Normalidad de los Errores: Se asume que los errores de predicción siguen una distribución normal con una media de cero (Ruido blanco con media cero). Aunque este supuesto puede ser relajado en algunos casos, es útil para interpretar los intervalos de confianza y realizar pruebas estadísticas.
Homocedasticidad: Los errores de predicción deben tener una varianza constante a lo largo del tiempo. Esto significa que la dispersión de los errores no debe cambiar a medida que avanza la serie temporal.
- Rojas Jimenez, K. 2022. Ciencia de Datos para Ciencias Naturales fragmentos del texto usados.
- Ptolomeo UNAM, Modelacion ARIMA fragmentos del texto usados.
SARIMA.
Definicion.
El modelo SARIMA (Seasonal AutoRegressive Integrated Moving Average) es una extensión del modelo ARIMA que también tiene en cuenta la estacionalidad (Seasonal). La notación general para un modelo SARIMA es SARIMA(p, d, q) (P, D, Q) (s);
\(p\): Orden del componente autoregresivo no estacional (AR).
\(d\): Orden de diferenciación no estacional, es decir, antes de ser estacional que orden de diferenciacion tenia, analogamente los demas.
\(q\): Orden del componente de media móvil no estacional (MA).
\(P\): Orden del componente autoregresivo estacional.
\(D\): Orden de diferenciación estacional.
\(Q\): Orden del componente de media móvil estacional.
\(s\): Longitud del período estacional.
\[ (1 - \phi_1 B - \phi_2 B^2 - \ldots - \phi_p B^p)(1 - B)^d(1 - \Phi_1 B^s - \Phi_2 B^{2s} - \ldots - \Phi_P B^{Ps})(1 - B^s)^D X_t = c + (1 + \theta_1 B + \theta_2 B^2 + \ldots + \theta_q B^q)(1 + \Theta_1 B^s + \Theta_2 B^{2s} + \ldots + \Theta_Q B^{Qs})Z_t \]
\(X_t\) es la serie temporal en el tiempo “t”.
\(B\) es el operador de retraso.
\(\phi_1, \phi_2, \ldots, \phi_p\) son los parámetros de autoregresión.
\(\Phi_1, \Phi_2, \ldots, \Phi_P\) son los parámetros de autoregresión estacionales.
\(\theta_1, \theta_2, \ldots, \theta_q\) son los parámetros de la parte de media móvil.
\(\Theta_1, \Theta_2, \ldots, \Theta_Q\) son los parámetros de la parte de media móvil estacional.
\(c\) es una constante.
\(Z_t\) es un ruido blanco. \end{itemize}
Analogo a los supuestos del ARIMA.
Estacionariedad.
Linealidad.
Independencia de los errores.
Normalidad de los errores.
Homocedasticidad.
Estacionalidad: Verificacion de los ciclos, y que estos sean estacionales.
Estacionalidad estacionaria: El componente estacional del modelo SARIMA debe ser estacionario. Esto implica que la estacionalidad debe ser modelada adecuadamente y no debe haber patrones sistemáticos en los residuos estacionales.
Practica.
Metodologia.
La metodologia usada en este tipo de analisis, es la siguiente mostrada en un diagrama:
Primero se identifica el modelo analizando la serie de tiempo a secas (lo que creemos que es, viendo sus ACF y PACF, etc).
Despues estimas los parametros en consecuencia de tu analisis, antes de este paso, es recomendable primero verificar que en efecto tu serie de tiempo cumpla con los supuestos de tu modelo, asi que en este paso tambien veras que tu serie de tiempo cumpla con ser estacionaria, homocedastica, univariada, etc. En este paso tambien se pueden probar diversos modelos, modelos con outliers (valores atipicos) y modelos que no tengan, esto para poder ver de manera mas representativa sus diferencias; y para su descarte o aceptacion, se puede ver un analisis adicional opcional ademas de los supuestos; AIC y BIC. *Deben de ser los minimos de todos ambos*.
Verificas los supuestos del modelo a traves de Pruebas Parametricas, verificas los mismos supuestos que se pretendian ver en la serie de tiempo a traves de un Analisis Residual, pero si es que no se cumplen, podemos optar por realizar tranformaciones o cambios a nuestro modelo o serie de tiempo.
Al final de hacer todo lo anterior, despues se usa el modelo para predecir, esto junto a un analisis de lo que representan estas predicciones en el contexto dado (buenas predicciones para el tema). Si aun con esto las Predicciones no son adecuadas, eso quiere decir que debemos cambiar nuestro enfoque de modelo, y probar otras alternativas (Modelos de Suavizado Exponencial SE, GARCH, Redes Neuronales RNN, INAR, Espacio de Estados SSM, Change Point.)
Una especial mencion al analisis de los autocorrelogramas o correlogramas o grafico de correlacion de datos.
Si en una FACP las primeras p correlaciones son significativas (distintas de cero o muy por arriba de cero), muy probablemente hablemos de un AR(p) y en su ACF decae a cero.
Si en una ACF las primeras q correlaciones son significativas (distintas de cero o muy por arriba de cero), quizas se trate de un MA(q), y PACF decae a cero.
Si en la ACF se de una convergencia a cero como se dice en la grafica, quizas se trate de un ARMA(p,q)
Juntando estos datos, con el factor de diferenciacion I(d) para volver estacionaria a tu serie de tiempo, y el factor S de estacionalidad, se obtiene el SARIMA por completo.
Supuestos.
Antes de modelar, solo tomaremos en cuenta el modelo SARIMA, primero debemos ver que la serie de tiempo original cumple con los supuestos SARIMA, es decir, cumple la Homocedasticidad, cumple la Estacionariedad, cumple con tener Ciclos Estacionales (Estacionalidad), y cumple con que su ruido blanco es gaussiano de media cero y sus errores son independientes (Normalidad e Independencia de los errores). Basicamente el supuesto de univariado, linealidad y multicolinealidad son menos explorados o se dan por hecho (triviales?).
Una vez la serie de tiempo cumple los supuestos en teoria, podemos modelarla estimando los parametros que deseemos (analizando que significa cada uno), y a este modelo SARIMA, le haremos pruebas parametricas, en general pruebas, para ver que tan robusto o aceptable es.
Pruebas de los Supuestos.
Ejemplo. Se cargaran datos, se aplicara un modelo y se le haran Pruebas, se analizara sus ACF y PACF y se hara un Analisis de Residuos.
ts_70 <- tsdl[[70]]
SARIMA <- auto.arima(ts_70, seasonal = TRUE) #, trace = TRUE si se quieren ver todos los modelos. Analisis de Residuos.
Primero hay que decir que se le puede ver a los residuos como los errores de la regresion, como los residuales precisamente o como el ruido blanco \(Z_t\) de la serie de tiempo/proceso \(X_t\).
ACF y PACF.
Primero vamos a interpretar los ACF y PACF.
El modelo anterior se dice que es un SARIMA (1,0,2)(2,1,0)[12], analizando sus ACF y PACF iniciales veremos de donde sale esto, vemos que el p no estacional nos indica que ciertamente en la PACF segun las bandas de confianza, el parametro p de lag igual a 1 es el significativo, q es algo mas incierto pero el ACF indica que ciertamente los primeros dos lags son considerablemente mas significativos que los demas dadas las bandas de confianza. Ahora para el analisis de las ACF y PACF residuales (las que nos interesan aqui), esperamos que la mayoria de las correlaciones se encuentren dentro de la banda de confianza para brindar un modelo fuerte, si hay diferencias muy grandes entre las correlaciones, podria indicar que el modelo no es lo suficientemente adecuado; Si los residuos parecen ser ruido blanco (es decir, no hay autocorrelación significativa). En este caso, es lo suficientemente aceptable.
acf(ts_70, main = "ACF de ts_70")
pacf(ts_70, main = "PACF de ts_70")
acf(residuals(SARIMA), main = "ACF residuales SARIMA")
pacf(residuals(SARIMA), main = "PACF residuales SARIMA")
Pruebas/Test de Supuestos.
Independencia.
Verificar si los residuos son mutuamente independientes: Utilizamos el estadístico Q de Ljung-Box para realizar una prueba conjunta sobre la significancia de las correlaciones de los residuos hasta un cierto lag. Si el valor p de la prueba es mayor que un nivel de significancia predefinido, entonces no hay evidencia suficiente para rechazar la hipótesis nula de que los residuos son independientes.
residuos <- residuals(SARIMA)
Box.test(residuos,lag = 20, type = "Ljung-Box")##
## Box-Ljung test
##
## data: residuos
## X-squared = 18.132, df = 20, p-value = 0.5787Aprobado.
No hay suficiente evidencia para rechazar la hipótesis nula de independencia de los residuos. Por lo tanto, según los resultados de la prueba de Box-Ljung, los residuos del modelo SARIMA parecen ser independientes en los primeros 20 rezagos/lags. Si NO cumple el supuesto, entonces, modificar parametros. Como persiste dependencia en los residuales, quiere decir que falto modelar dependencia en el modelo. Tambien se puede intentar agregando mas lags/rezagos, asi saldria mayor el p-value, pero ver hasta que punte es permisible.
Homocedasticidad.
Verificar si la varianza de los residuos es constante o no constante (heterocedasticidad): Se puede graficar los residuos contra el tiempo y observar si la dispersión de los residuos parece constante a lo largo del tiempo. O hacer pruebas, como la de Breusch-Pagan.
n <- length(residuos)
regresor <- 1:n
lmtest::bptest(residuos~regresor)##
## studentized Breusch-Pagan test
##
## data: residuos ~ regresor
## BP = 1.1217, df = 1, p-value = 0.2895Aprobado.
No hay suficiente evidencia para rechazar la hipótesis nula de homocedasticidad. Esto sugiere que la varianza de los residuos es constante a lo largo del tiempo, osease, sugiere que si se cumple homocedasticidad y rechazamos heterocedasticidad. Si NO se cumple entonces realizar una transformación Box-Cox, sobre la serie original.
Estacionariedad.
Verificar si los residuos cumplen la propiedad estacionaria: que la varianza sea constante, que la media constante, y que la correlacion solo dependa de la distancia del rezago/lag: Se considerara que si se cumple este supuesto entonces no hay presencia de tendencia/media/pendiente (es constante o se suele decir que es cero como la media del ruido blanco igual a cero).
Se pueden verificar que la media es constante a lo largo del tiempo, mediante el calculo de la media, la homocedasticidad se comprueba anteriormente y se puede probar la correlacion y la independencia de los errores para este supuesto, mediante los ACF y PACF; o se pueden hacer pruebas como Dickey-Fuller aumentada (ADF) o la prueba de Kwiatkowski-Phillips-Schmidt-Shin (KPSS). Estas pruebas evalúan si una serie temporal tiene raíces unitarias (ADF) o si la serie es estacionaria en tendencia (KPSS).
tseries::adf.test(residuos)##
## Augmented Dickey-Fuller Test
##
## data: residuos
## Dickey-Fuller = -6.5685, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationarytseries::kpss.test(residuos)##
## KPSS Test for Level Stationarity
##
## data: residuos
## KPSS Level = 0.056614, Truncation lag parameter = 5, p-value = 0.1Aprobado.
Prueba de Dickey-Fuller Aumentada (ADF): segun el warning: p-value smaller than printed p-value; el valor-p obtenido es menor que 0.01. Esto indica que hay suficiente evidencia para rechazar la hipótesis nula de presencia de una raíz unitaria en los residuos, lo que sugiere que los residuos son estacionarios.
Prueba KPSS para Estacionariedad en Nivel (Level Stationarity): - segun el warning: p-value greater than printed p-value; el valor-p obtenido es mayor que 0.1. No hay suficiente evidencia para rechazar la hipótesis nula de estacionariedad en tendencia, lo que sugiere que los residuos son estacionarios.
Si NO se cumplen estos supuestos, se pueden hacer transformaciones a la serie original, haciendola estacionaria desde un principio, etc, aunque esto simplifica los modelos.
Normalidad.
Verificar si los residuos tienen distribución normal (gaussiana): Se puede graficar un qq-plot de los residuos con una línea diagonal para verificar si los residuos siguen una distribución normal visualmente, deben estar pegados a la linea, graficar su histograma y darse idea visual de sus frecuencias y si su densidad añadida al histograma parece ser normal; o se pueden hacer pruebas de normalidad de bondad de ajuste como shapiro de normalidad que tiene buena potencia o Anderson-Darling o Lilliefors.
shapiro.test(residuos)##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.82639, p-value < 2.2e-16lillie.test(residuos)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuos
## D = 0.14942, p-value < 2.2e-16ad.test(residuos)##
## Anderson-Darling normality test
##
## data: residuos
## A = 19.77, p-value < 2.2e-16Reprobado.
En todas las pruebas, el valor-p obtenido es significativamente menor que cualquier nivel de significancia comúnmente utilizado (por ejemplo, 0.05). Por lo tanto, hay suficiente evidencia para rechazar la hipótesis nula de normalidad en todos los casos. Esto sugiere que los residuos no siguen una distribución normal y que hay desviaciones significativas de la normalidad.
En este caso, NO pasa la normalidad, asi que podemos hacer intervalos de predicción distintos a los de la distribución normal, los intervalos de predicción basados en la distribución normal pueden no ser apropiados en este caso. Algunas opciones para construir intervalos de predicción en presencia de residuos no normales son Bootstrap como metodo de remuestreo, o transformaciones para mejorar los residuos; por ejemplo podriamos ajustar transformaciones antes de modelar el SARIMA o en los residuos transformados con la funcion boxcox de la libreria MASS.
AIC y BIC.
Extra. Criterio de Información de Akaike (AIC) y el Criterio de Información Bayesiano (BIC): Trataremos de buscar el modelo que tenga los mas bajos de estos en conjunto, estos criterios, el primero penaliza menos la complejidad y el segundo lo hace mas.
Forecast.
Digamos que ya tenemos un modelo que cumple con todo, ahora sigue hacer predicciones, es decir usar el modelo, hay diferentes formas de predecir, pero en este caso usaremos las siguientes paqueterias.
pronostico<-forecast(SARIMA, h = 12)
plot(pronostico, main = "Pronóstico SARIMA", xlab = "Tiempo", ylab = "Valores")
legend("topright", legend = c("Pronóstico", "Area de confianza"), col = c("blue", "grey"), lty = 1:2)