Ejemplo sobre análisis y predicción de una serie temporal

Más info: https://rpubs.com/palominoM/series

Presentación del ejemplo y análisis inicial

Ejemplo: Analizamos una serie sobre las emisiones de CO2

Cargamos las librerías necesarias

library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)
library(ggfortify)
## Loading required package: ggplot2
## Registered S3 methods overwritten by 'ggfortify':
##   method                 from    
##   autoplot.Arima         forecast
##   autoplot.acf           forecast
##   autoplot.ar            forecast
##   autoplot.bats          forecast
##   autoplot.decomposed.ts forecast
##   autoplot.ets           forecast
##   autoplot.forecast      forecast
##   autoplot.stl           forecast
##   autoplot.ts            forecast
##   fitted.ar              forecast
##   fortify.ts             forecast
##   residuals.ar           forecast

Cargamos la dataset

co2<-read.csv("co2.csv", header = TRUE, sep = ";")
#Transformamos los datos en una serie temporal 
co2ts<-ts(co2$CO2, start = c(1959,1), frequency = 12)

Ahora de la libreria ggfortify dibujamos la serie temporal

autoplot(co2ts)

Es evidente que se trata de una serie temporal con tendencia alcista y que ciclos internos que deberemos analizar.

autoplot(acf(co2ts, plot = FALSE))

Se puede apreciar en esta gráfica como la serie no es estacionaria ya que el valor de la función de autocorrelación no decae de manera exponencial a medida que aumentan los rezagos (períodos) en el tiempo. También apreciamos que existe la componente estacional, con un periodo de 12 (evaluamos al analizar los máximos y el número de períodos entre máximos), por lo tanto se analizará las componente por separado para ello descomponemos la serie.

autoplot(stl(co2ts, s.window = "periodic"))

Como ya vimos en la primera parte de esta serie de tutoriales con R, la descomposición entre una serie temporal se compone de: - Tendencia (media móvil) - Estacionalidad (períodos entre máximos) - Ruido.

Análisis ARIMA

Analizamos las funciones ndiffs y nsdiffs, que calculan cada una el número de diferenciaciones regulares y estacionales respectivamente que se necesita llevar a cabo para que la serie sea estacionaria.

autoplot(diff(co2ts), linetype = "dashed", colour = "blue")

Como podemos ver ahora en la serie ya se ha eliminado la componente de tendencia, una ves eliminada esta componente volvemos a gráficar la función de autocorrelación.

autoplot(acf(diff(co2ts), plot = FALSE), colour = 'brown')

Vuelven a pasar 12 rezagos entre máximos, no es estacionaria en la parte estacional.

Podemos ilustrar los valores entre meses con box-plots o con gráficas lineales.

monthplot(diff(co2ts), col = "midnightblue")

O box-plots

boxplot(diff(co2ts)~cycle(diff(co2ts)))

Existen picos en abril (mes 4) y debemos eliminar esta componente estacional.

diff.co2ts.12<-diff(co2ts, lag = 12)
autoplot(diff.co2ts.12, colour = "darkorange4", linetype = "dashed")

Una vez eliminado tanto la componente de tendencia y la estacional, observemos que esta serie se parece bastante a una serie estacionaria, ya que parece ser constante en media y varianza, para asegurarnos aplicamos un test de estacionariedad como el test ADF (Dickey-Fuller) y el test de KPSS (Kwiatkowski-Phillips-Schmidt-Shin).

La prueba de Dickey-Fuller aumentada es un tipo de prueba estadística llamada prueba de raíz unitaria .

En teoría de probabilidad y estadística, una raíz unitaria es una característica de algunos procesos estocásticos (como los paseos aleatorios) que pueden causar problemas en la inferencia estadística que involucra modelos de series de tiempo. En términos simples, la raíz unitaria no es estacionaria pero no siempre tiene un componente de tendencia .

La prueba del ADF se realiza con los siguientes supuestos:

Hipótesis nula (HO) : la serie no es estacionaria o tiene una raíz unitaria. Hipótesis alternativa (HA) : la serie es estacionaria o no tiene raíz unitaria. Si no se rechaza la hipótesis nula, esta prueba puede proporcionar evidencia de que la serie no es estacionaria.

Condiciones para rechazar la hipótesis nula (HO)

Si el estadístico de prueba <valor crítico y el valor p <0,05, rechace la hipótesis nula (HO), es decir, la serie temporal no tiene raíz unitaria, lo que significa que es estacionaria. No tiene una estructura dependiente del tiempo. Prueba de Dickey-Fuller Antes de pasar a la prueba ADF, primero comprendamos qué es la prueba Dickey-Fuller.

Utiliza un modelo autorregresivo y optimiza un criterio de información en múltiples valores de retraso diferentes. Una prueba de Dickey-Fuller es una prueba de raíz unitaria que prueba la hipótesis nula de que α=1 en la siguiente ecuación modelo. alphaes el coeficiente del primer rezago en Y.

Hipótesis nula (H0): alfa=1

dónde,

y(t-1) = rezago 1 de la serie temporal delta Y(t-1) = primera diferencia de la serie en el momento (t-1) Fundamentalmente, tiene una hipótesis nula similar a la prueba de raíz unitaria. Es decir, el coeficiente de Y(t-1) es 1, lo que implica la presencia de una raíz unitaria. Si no se rechaza, la serie se considera no estacionaria.

La prueba de Dickey-Fuller aumentada evolucionó basándose en la ecuación anterior y es una de las formas más comunes de la prueba de raíz unitaria.

adf.test(diff.co2ts.12)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff.co2ts.12
## Dickey-Fuller = -3.7783, Lag order = 5, p-value = 0.02201
## alternative hypothesis: stationary
kpss.test(diff.co2ts.12)
## Warning in kpss.test(diff.co2ts.12): p-value greater than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  diff.co2ts.12
## KPSS Level = 0.11873, Truncation lag parameter = 4, p-value = 0.1

Resultados contradictorios, el primer test nos indica que la serie es estacionaria y la segunda lo contrario.

¿Cuándo elegir la prueba ADF o KPSS? Podría haber mucha confusión sobre cuándo se debe utilizar la prueba ADF o la prueba KPSS y qué prueba daría un resultado correcto. Una mejor solución es aplicar/ejecutar ambas pruebas y asegurarse de que la serie sea realmente estacionaria.

Los siguientes son los posibles resultados de la aplicación de ambas pruebas.

Caso 1 : Ambas pruebas concluyen que la serie dada es estacionaria – La serie es estacionaria Caso 2 : Ambas pruebas concluyen que la serie dada no es estacionaria – La serie no es estacionaria Caso 3 : ADF concluye que no es estacionario y KPSS concluye que es estacionario : la serie es de tendencia estacionaria. Para que la serie sea estrictamente estacionaria, en este caso es necesario eliminar la tendencia. Luego se verifica la estacionariedad de la serie sin tendencia. Caso 4 : ADF concluye estacionario y KPSS concluye no estacionario : la serie es estacionaria en diferencias. La diferenciación se utilizará para hacer que las series sean estacionarias. Luego se verifica la estacionariedad de la serie diferenciada.

autoplot(acf(diff.co2ts.12, plot = FALSE)) # Observar el lag order 5

autoplot(pacf(diff.co2ts.12, plot = FALSE))

Vistas la gráficas acf y pacf pueden plantear varios modelos tentativos para el análisis, como lo serían los ARIMA(0,1,2)(0,1,1) o ARIMA(1,1,0)(2,1,0) deducidos, o realizando una combinación de estos como ARIMA(1,1,2)(2,1,1)) o incluso variando algunos de los componentes como puede ser con ARIMA(1,1,1)(2,1,1)) o ARIMA(1,1,2)(1,1,1)) o ARIMA(0,1,1)(0,1,1) o ARIMA(1,1,0)(1,1,0), etc.

arima1<- Arima(co2ts, order=c(0,1,2), seasonal=list(order=c(0,1,1),period=12))
arima2<- Arima(co2ts, order=c(1,1,0), seasonal=list(order=c(2,1,0),period=12))
arima3<- Arima(co2ts, order=c(1,1,2), seasonal=list(order=c(2,1,1),period=12))
arima4<- Arima(co2ts, order=c(1,1,1), seasonal=list(order=c(2,1,1),period=12))
arima5<- Arima(co2ts, order=c(1,1,2), seasonal=list(order=c(1,1,1),period=12))
arima6<- Arima(co2ts, order=c(0,1,1), seasonal=list(order=c(0,1,1),period=12))
arima7<- Arima(co2ts, order=c(1,1,0), seasonal=list(order=c(1,1,0),period=12))

VAmos a evaluar los 7 modelos mediante el parámetro de testeo AIC. El criterio de información de Akaike (AIC) o aic bic es un estimador de la calidad relativa del modelo que tiene en cuenta su complejidad. A medida que se va aumentando el número de parámetros de entrada o de grados de un polinomio, el valor de R va a ser mejor, porque el error cuadrático medio disminuye. El criterio de información Akaike (aic metric) penaliza los modelos complejos en favor del los sencillos para evitar el sobreajuste. A menores valores, mejor es el modelo.

AIC(arima1,arima2,arima3,arima4,arima5,arima6,arima7)
##        df       AIC
## arima1  4  85.72537
## arima2  4  93.70131
## arima3  7  90.55929
## arima4  6  89.02326
## arima5  6  88.62985
## arima6  3  83.82626
## arima7  3 106.19086

El Criterio de Información Bayesiano , a menudo abreviado BIC , es una métrica que se utiliza para comparar la bondad del ajuste de diferentes modelos de regresión. En la práctica, ajustamos varios modelos de regresión al mismo conjunto de datos y elegimos el modelo con el valor BIC más bajo como el modelo que mejor se ajusta a los datos. A menores valores, mejor es el modelo.

BIC(arima1,arima2,arima3,arima4,arima5,arima6,arima7)
##        df       BIC
## arima1  4  97.71422
## arima2  4 105.69016
## arima3  7 111.53978
## arima4  6 107.00654
## arima5  6 106.61312
## arima6  3  92.81789
## arima7  3 115.18250

Los valores máximos se alcanzan en el modelo “arima7” con una combinación ARIMA(1,1,0) y los mínimos se alcanzan con arima6 ARIMA(0,1,1) y eso quiere decir que no tiene componente autoregresiva.

Vamos a graficar los correlogramas de los residuos para comprobar que son ruido blanco. Un Ruido Blanco es una serie temporal que no tiene un patrón y por lo cual no se puede predecir el futuro.

Para que la serie de tiempo sea Ruido Blanco debe tener:

Media constante.

Varianza constante.

Sin autocorrelaciones en ningún período.

El Ruido Blanco es una secuencia de datos aleatorios donde cada valor tiene un período de tiempo asociado y la serie de tiempo se tiene un comportamiento aleatorio y no se puede proyectar el futuro.

autoplot(acf(arima6$residuals, plot = FALSE))

autoplot(pacf(arima6$residuals, plot = FALSE))

Ni en la gráfica de autocorrelación (ACF) ni en la parcial (PACF) se aprecian significancia (no atraviesan la línea azul). Por lo tanto podemos decir que los residuos son ruido blanco. Ahora, vamos a pintar también una serie de gráficas sobre los residuos:

ggtsdiag(arima6)

Al llevar a cabo el modelo de pronóstico ARIMA siempre es necesario analizar el comportamiento de los residuos, en este caso nos interesa analizar si estos residuos se comportan como ruido blanco. Para ello se recomienda llevar a cabo la prueba de autocorrelación de forma conjunta de Ljung-Box (la 3a gráfica, superior). No se aprecia significancia

Box.test(arima6$residuals, type="Ljung-Box") # Lo confirmamos de manera analítica
## 
##  Box-Ljung test
## 
## data:  arima6$residuals
## X-squared = 0.1922, df = 1, p-value = 0.6611

La prueba estadística de Box-Pierce es una versión simplificada de la estadística de Ljung-Box para los cuales los estudios de simulación posteriores han demostrado un rendimiento deficiente.

Box.test(arima6$residuals) # Test de Box-Pierce
## 
##  Box-Pierce test
## 
## data:  arima6$residuals
## X-squared = 0.18866, df = 1, p-value = 0.664

Auto-Arima

Existe una función en R para calcular la mejor combinación de órdenes para ARIMA que comprobaremos:

auto.arima(co2ts, stepwise = FALSE, approximation = FALSE)
## Series: co2ts 
## ARIMA(1,0,1)(0,1,1)[12] with drift 
## 
## Coefficients:
##          ar1      ma1     sma1   drift
##       0.8999  -0.4391  -0.7686  0.1193
## s.e.  0.0479   0.0885   0.0930  0.0034
## 
## sigma^2 = 0.09903:  log likelihood = -35.58
## AIC=81.16   AICc=81.58   BIC=96.18

El AIC del modelo arima6 arrojaba un resultado de 83 (major que el calculado 81) y el BIC fue de 92 (peor aquí con 96). Comparamos la predicción de los dos modelos propuestos Un primer modelo del tipo (0,1,1)(0,1,1) que corresponde al modelo arima6, y otro calculado automáticamente (0,1,0)(0,1,1).

Por defecto, se utiliza CSS (suma condicional de cuadrados), es posible que los coeficientes autorregresivos no sean estacionarios (es decir, queden fuera de la región para procesos estacionarios).

Forzamos el método MLE (maximum likelihood estimation), esto es más lento (inapreciable en nuestro caso) pero proporciona mejores estimaciones y siempre devuelve un modelo estacionario.

arima6<- Arima(co2ts, order=c(0,1,1), seasonal=list(order=c(0,1,1),period=12), method="ML")
forecast1<-forecast(arima6, level = c(95), h = 50)
autoplot(forecast1)

arima_auto<- Arima(co2ts, order=c(1,0,1), seasonal=list(order=c(0,1,1),period=12), method="ML")
forecast2<-forecast(arima_auto, level = c(95), h = 50)
autoplot(forecast2)

AIC(arima6,arima_auto)
## Warning in AIC.default(arima6, arima_auto): models are not all fitted to the
## same number of observations
##            df      AIC
## arima6      3 83.82626
## arima_auto  4 89.22061

Pues mejor arima6

BIC(arima6,arima_auto)
## Warning in BIC.default(arima6, arima_auto): models are not all fitted to the
## same number of observations
##            df       BIC
## arima6      3  92.81789
## arima_auto  4 101.23640

Está claro que el estimador automático no ha funcionado.