Si vivieramos en un mundo perfecto todas las series temporales serían estacionarias, con su media y su varianza constantes, no existiría la incertidumbre de su evolución. Pero, en el mundo real la mayoría de las series temporales NO son estacionarias. Si la diferenciamos n veces la serie original se transforma en estacionaria, y entonces ya podremos aplicar metodologías Box-Jenkins, dentro de la cuales nos centraremos en los modelos ARIMA.

¿Qué es un modelo ARIMA?

ARIMA es la abreviatura de “AutoRegressive Integrated Moving Average”, consiste en un modelo estadístico que utiliza variaciones y regresiones de datos estadísticos con el fin de encontrar patrones para una predicción hacia el futuro, es decir, las estimaciones futuras vienen explicadas por los datos del pasado y no por variables independientes.

Componentes

- AutoReressive (AR(p)) componente autorregresivo, se refiere a la utilización de valores pasados en la ecuación de regresión para la serie de \(Y_{t}\) . El parámetro \(\rho\) autoregresivo especifica el número de retardos utilizados en el modelo.

- Integrated (I(d)) representa el grado de diferenciación, es decir, el número de diferencias utilizadas para hacer que las series temporales sean estacionarias.

- Moving Average (MA(q)) representa el error del modelo como una combinación de los términos de error anteriores \(\varepsilon_{t}\). El orden \(q\) determina el número de términos a incluir en el modelo.

Los componentes de diferenciación, autorregresión y promedio móvil conforman un modelo ARIMA(p,d,q), no estacional, que se puede escribir como una ecuación lineal: \[Y_{t}=c+\varnothing_{1}y_{dt-1}+\varnothing_{p}y_{dt-p}+...+\theta_{1}e_{t-1}++\theta_{q}e_{t-q}+e_{t}\] donde \(y_{d}\) es \(Y\) diferenciado \(d\) veces y \(c\) es una constante

Pasos a seguir para el modelado

1.Realizar un análisis exploratorio de los datos.

Examinamos los datos buscando patrones o irregularidades, limpiamos valores atípicos o faltantes,descomponemos la serie para observar si existe tendencia y estacionalidad. (Es conveniente disponer de 50 o más datos).

2.Ajuste del modelo.

Para ello examinamos la estacionariedad, si lo requiere, aplicamos las transformaciones para convertir la serie a estacionaria, por último, estudiamos los diagramas de correlación (ACF y PACF) seleccionamos el orden del modelo.

3.Evaluación el modelo ARIMA (p,d,q) seleccionado.

Si los parámetros de orden del modelo y la estructura se especifican correctamente, no deberían existir autocorrelaciones significativas, por ello, verificamos los residuos y las trazas ACF/PACF.

4.Reajuste del modelo (si fuera necesario)

Seleccionamos unos nuevos parámetros y compara los errores del modelo y los criterios de ajuste como AIC o BIC.

5.Predicción

Una vez validado el modelo podremos realizar previsiones de futuros valores.

Ejemplo en R

Análisis exploratorio de los datos

El objetivo es el análisis, modelado y predicción de una base de datos que muestra las concentraciones de CO2 mensuales a lo largo de 14 años. Lo primero que haremos será leer nuestros datos, los convertiremos en una serie temporal y la representaremos gráficamente.

#Cargamos los datos co2
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)
print(co2ts)
##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 1959 333.13 332.09 331.10 329.14 327.36 327.29 328.23 329.55 330.62 331.40
## 1960 333.92 333.43 331.85 330.01 328.51 328.41 329.25 330.97 331.60 332.60
## 1961 334.68 334.17 332.96 330.80 328.98 328.57 330.20 331.58 332.67 333.17
## 1962 336.82 336.12 334.81 332.56 331.30 331.22 332.37 333.49 334.71 335.23
## 1963 337.95 338.00 336.37 334.47 332.46 332.29 333.76 334.80 336.00 336.63
## 1964 339.05 339.27 337.64 335.68 333.77 334.09 335.29 336.76 337.77 338.26
## 1965 341.47 341.31 339.41 337.74 336.07 336.07 337.22 338.38 339.32 340.41
## 1966 343.02 342.54 340.88 338.75 337.05 337.13 338.45 339.85 340.90 341.70
## 1967 344.28 343.42 342.02 339.97 337.84 338.00 339.20 340.63 341.41 342.68
## 1968 345.92 345.40 344.16 342.11 340.11 340.15 341.38 343.02 343.87 344.59
## 1969 347.38 346.78 344.96 342.71 340.86 341.13 342.84 344.32 344.88 345.62
## 1970 348.53 347.87 346.00 343.86 342.55 342.57 344.11 345.49 346.04 346.70
## 1971 349.93 349.26 347.44 345.55 344.21 343.67 345.09 346.27 347.33 347.82
## 1972 351.71 350.94 349.10 346.77 345.73                                   
##         Nov    Dec
## 1959 331.87 333.18
## 1960 333.57 334.72
## 1961 334.86 336.07
## 1962 336.54 337.79
## 1963 337.93 338.95
## 1964 340.10 340.88
## 1965 341.69 342.51
## 1966 342.70 343.65
## 1967 343.04 345.27
## 1968 345.11 347.07
## 1969 347.23 347.62
## 1970 347.38 349.38
## 1971 349.29 350.91
## 1972
#Trazamos la serie de tiempo co2ts
autoplot(co2ts, ts.colour = "blue", ts.linetype = "dashed")

Una vez trazada la serie co2ts se observa que presenta una tendencia ascendente (la serie no es estacionaria en media). Con esto podemos afirmar que la serie no es estacionaria. Por otro lado, no se aprecia gran variabiliad, esto quiere decir que la serie es estacionaria en cuanto a su varianza. Para confirmar lo dicho trazamos la gráfica de su función de autocorrelación.

autoplot(acf(co2ts, plot = FALSE))

En esta gráfica apreciamos 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 en el tiempo. También apreciamos que existe la componente estacional, con un periodo de 12, por lo tanto, analizaremos las componentes por separado, para ello descomponemos la serie.

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

Ajuste del modelo

El análisis previo revela que tenemos que eliminar la tendencia y la estacionalidad de la serie para que pueda ser estacionaria. El modelo que vamos a utilizar es un modelo ARIMA, diferenciando la serie lograremos que se convierta a estacionaria. Utilizaremos las funciones ndiffs y nsdiffs, que calculan el número de diferenciaciones regulares y estacionales, respectivamente, que se necesita llevar a cabo para que la serie sea estacionaria.

ndiffs(co2ts)
## [1] 1
nsdiffs(co2ts)
## [1] 1

Los respectivos cálculos nos muestran que la serie necesita una diferenciación regular y otra estacional.

diff.co2ts<-autoplot(diff(co2ts), ts.linetype = "dashed", ts.colour = "blue")
diff.co2ts

Como podemos apreciar ya se ha eliminado la componente de tendencia, graficamos nuevamente la función de autocorrelación.

autoplot(acf(diff(co2ts), plot = FALSE))

Podemos apreciar que la serie sigue sin ser estacionaria, en la parte estacional, se vuelve a ver la componente estacional de perido 12.

Ahora eliminaremos esta componente.

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

Una vez eliminadas la tendencia y la componente estacional, observamos que la serie se parece bastante a una serie estacionaria, puesto que parece ser constante en media y varianza, para asegurarnos aplicamos el test ADF (Dickey-Fuller).

Contrastamos la hipótesis para un \(\alpha=0.05\)

\(H_{0}:\) La serie no es estacionaria
\(H_{1}:\) La serie es estacionaria

adf<-adf.test(diff.co2ts.12)
adf$p.value
## [1] 0.02201398

El p-valor obtenido es de 0.022 < 0.05. Rechazamos H0, confirmamos que la serie ya es estacionaria.

con el fin de seleccionar los órdenes (p) y (q) del modelo ARIMA examinamos las tramas de ACF y PACF (autocorrelaciones simple y parcial).

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

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

A partir de la visualización de de los diagramas de autocorrelación simple y parcial podríamos probar modelos AR o MA 1 y 2 .

Evaluación del modelo ARIMA selecionado

Voy a comparar el AIC de una serie de modelos propuestos para seleccionar el mejor de ellos y validarlo.

#Mis modelos
arima1<- arima(co2ts,order=c(2,1,0), method = "ML")
arima1$aic
## [1] 380.3415
arima2<- arima(co2ts,order=c(2,1,1), method = "ML")
arima2$aic
## [1] 345.036
arima3<- arima(co2ts,order=c(2,1,2), method = "ML")
arima3$aic
## [1] 371.076
arima4<- arima(co2ts,order=c(1,1,2), method = "ML")
arima4$aic
## [1] 380.2622

El modelo que posee el menor AIC es el Modelo ARIMA (2,1,1).

Antes de proceder a su validación voy a utilizar la función auto.arima, esta función proporciona automáticamente el mejor Modelo ARIMA (orden óptimo (p,d,q)) según el valor AIC, AICc o BIC, considerando: estacionariedad, estacionalidad, diferencias, entre otras.

mod_arima <- auto.arima(co2ts, seasonal = FALSE)
mod_arima
## Series: co2ts 
## ARIMA(2,1,1) 
## 
## Coefficients:
##          ar1      ar2      ma1
##       1.5430  -0.8356  -0.8137
## s.e.  0.0439   0.0422   0.0375
## 
## sigma^2 estimated as 0.4842:  log likelihood=-168.52
## AIC=345.04   AICc=345.29   BIC=357.34

Confirmamos que el modelo que posee el mejor ajuste es el Modelo ARIMA (2,1,1).Ahora procedemos a validarlo y Para ello en primer lugar se grafican los correlogramas de los residuos para comprobar que son ruido blanco:

tsdisplay(residuals(mod_arima), lag.max=12, main='Residuos (2,1,1)')

El análisis de los residuos indica que existe un patrón en ACF que se repite en el retardo 3. Esto indica que nuestro modelo podría mejorar con una especificación diferente como p=3 o q=3. Podriamos repetir el proceso de ajuste que permite el componente MA (3) y examinar los diagramas de diagnóstico nuevamente

En nuestro caso, puesto que la componente estacional era de periodo vamos a ajustar a un Modelo ARIMA ESTACIONAL o SARIMA. Este modelo incorpora factores no estacionales y estacionales en un modelo multiplicativo. Una notación abreviada para el modelo es: \[ ARIMA ( p, d, q ) × ( P, D, Q ) S \]

REAJUSTE DEL MODELO

Probamos un SARIMA

Vamos a volver a comprobar la estacionariedad pero esta vez mediante el test de KPSS (Kwiatkowski-Phillips-Schmidt-Shin).

Contrastamos la hipótesis para un \(\alpha=0.05\)

\(H_{0}:\) La serie no es estacionaria
\(H_{1}:\) La serie es estacionaria

kpss<-kpss.test(diff.co2ts.12)
## Warning in kpss.test(diff.co2ts.12): p-value greater than printed p-value
kpss$p.value
## [1] 0.1

El p-valor obtenido es de 0.1 > 0.05. ** NO rechazamos H0**, con lo cual nuestra serie continuaba sin ser estacionaria. Por este motivo vamos a implementar componentes estacionales en nuesto modelo.

Para no extender mucho más el análisis, vamos a seleccionar el modelo automáticamente

mod_sarima <- auto.arima(co2ts, stepwise = TRUE, approximation = TRUE) #Selección por pasos y estimacion max. verosimilitud
mod_sarima
## Series: co2ts 
## ARIMA(1,0,1)(2,1,2)[12] with drift 
## 
## Coefficients:
##          ar1      ma1     sar1     sar2     sma1     sma2   drift
##       0.8880  -0.4252  -0.4862  -0.0029  -0.2211  -0.5112  0.1197
## s.e.  0.0529   0.0949   0.3324   0.1608   0.3334   0.3084  0.0029
## 
## sigma^2 estimated as 0.09842:  log likelihood=-34.46
## AIC=84.92   AICc=85.95   BIC=108.96
coeftest(mod_sarima)
## 
## z test of coefficients:
## 
##         Estimate Std. Error z value  Pr(>|z|)    
## ar1    0.8880411  0.0528888 16.7907 < 2.2e-16 ***
## ma1   -0.4251980  0.0948555 -4.4826 7.374e-06 ***
## sar1  -0.4862294  0.3323727 -1.4629    0.1435    
## sar2  -0.0029438  0.1607876 -0.0183    0.9854    
## sma1  -0.2210923  0.3334456 -0.6631    0.5073    
## sma2  -0.5112314  0.3084164 -1.6576    0.0974 .  
## drift  0.1196767  0.0029152 41.0524 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Una vez estimado el modelo, en este caso ARIMA(1,0,1)(2,1,2), con un AIC = 81.16 y un BIC = 96.18 se procede a validarlo. Representando los correlogramas de los residuos para comprobar que son ruido blanco:

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

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

Se puede apreciar en los correlogramas que no hay ningún retardo significativo (aparte del 0 del ACF, que es 1 por definición) que denote ningún tipo de estructura, por lo tanto podemos decir que los residuos son ruido blanco.

ggtsdiag(mod_sarima)

Los p-valores para la prueba Q de Ljung-Box están por encima de 0,05, lo que indica “no significativo”.

Realizamos el contraste de hipótesis:

\(H_{0}:\) Los datos se distribuyen de forma independiente
\(H_{1}:\) Los datos no se distribuyen de forma independiente

independencia <- Box.test(mod_sarima$residuals, type="Ljung-Box") # Test de Ljung-Box
independencia$p.value
## [1] 0.6187975

Efectivamente, los datos se distribuyen de forma independiente. P-valor = 0.619>0.05

qqnorm(mod_sarima$residuals)
qqline(mod_sarima$residuals) 

Gráficamente observamos que los datos siguen una distribución normal aunque en las colas los datos se alejan y se dispersan un poco, para una mayor seguridad comprobamos la normalidad mediante el test Shapiro Wilk:

\(H_{0}:\) Los datos se distribuyen normalmente
\(H_{1}:\) Los datos no se distribuyen normalmente

normalidad <-shapiro.test(mod_sarima$residuals)    # Test de Shapiro-Wilk
normalidad$p.value  
## [1] 0.1456028

Comprobamos que el p-valor del test es 0.146>0.05 no rechazamos nuestra hipótesis nula, los residuos siguen una distribución normal.

Por lo tanto, damos el modelo por VÁLIDO

Predicción

prediccion <- forecast(mod_sarima, h=36) #nivel confianza 95%, h = periodos
autoplot(prediccion)

Los pronósticos se muestran como una línea azul, entre los intervalos de predicción del 80% y del 95% como un área sombreada.