Introducción a las series temporales (2): modelo ARIMA con R

En este artículo vamos a ahondar más en técnicas de modelado de series temporales. En concreto veremos los modelos ARIMA.

Si eres nuevo/a en este tema, te recomiendo que visites primero este post en el que hemos introducido el concepto de serie temporal. Además, se ilustró el modelado estructural, es decir, descomponer la serie en componentes (tendencia, estacionalidad, etc)

En este artículo, veremos otras formas de modelarlas. Trabajaremos ejemplos prácticos con R. Además haremos predicciones (forecasting).

El proceso de modelado consiste en encontrar una forma matemática que explique el comportamiento de la serie para:

• Predecir futuros valores (forecasting)
• Aceptar o rechazar una teoría previa, por ejemplo, en macroeconomía.

Conceptos básicos

En primer lugar, consideramos las series temporales como procesos estocásticos. Cada observación en un proceso estocástico es una variable aleatoria y las observaciones evolucionan acorde a leyes probabilísticas.

Es decir, un proceso estocástico es simplemente una colección de variables aleatorias ordenadas en el tiempo.

Ruido blanco

Para que una señal o serie temporal sea considerada un ruido blanco (denotado con la letra epsilon), debe cumplir las siguientes características:

E(\epsilon_{t})=0

Var(\epsilon_{t}) = \sigma^2

Cov(\epsilon_{t},\epsilon_{t+k}) =0

E es la esperanza o media, que debe ser 0. Var es la varianza, que debe ser constante, ya que sigma es un valor fijo. Y cov es la covarianza entre dos muestras separas k unidades temporales.

El ruido blanco es importante, porque la parte que no se pueda modelar de una serie temporal, se considerará ruido blanco. Aunque trabajemos con series no estacionarias, se intentará buscar la forma de transformarla en estacionaria de forma que pueda ser modelada.

Estacionariedad

En primer lugar, no debemos confundir con el concepto explicado de estacionalidad, la cuál es la variación periódica y predecible de la misma con un periodo inferior o igual a un año.

Para que una serie temporal sea considerada estacionaria debe cumplir las siguientes propiedades:

 E(y_{t}) = \mu

Var(y_{t}) = V_{0}

Cov(y_{t},y_{t+k}) =\gamma_{k}

En otras palabras:

• Media constante
• Varianza constante. También conocido como homocedasticidad.
• Autocovarianza constante

Estas series en rojo, no son estacionarias. La primera, la de la izquierda, no cumple la propiedad de media constante, ya que vemos que tiene una tendencia positiva y ascendente. La serie del medio no tiene varianza constante, puesto que la amplitud de onda en el medio es superior que en los extremos. Por último, la última, no tiene autocovarianza constante, puesto que la longitud de onda es inferior en el medio que en los extremos.

A continuación mostramos una serie completamente estacionaria:

Operador de retardo o diferencias

Cuando nos enfrentamos a una serie temporal no estacionaria, es muy habitual aplicar diferencias sobre la serie original para eliminar la tendencia de la serie. De forma matemática el operador de retardos (L) se define así:

Ly_t = y_{t-1}

Es decir, al aplicar el operador L se obtiene las observaciones retardadas un período. De forma general podemos decir que:

L^d y_t = y_{t-d}

Es decir, retardamos d veces la muestra actual.

Por otra parte, podemos definir el operador diferencia como (1-L), de forma que:

y_t - y_{t-1} = (1-L)y_t

Como hemos comentado, esta técnica de aplicar diferencias se empleará cuando la serie sea no estacionaria. Además también podremos aplicar diferencias sobre el logaritmo de la serie. Veamos un ejemplo.

Análisis con un ejemplo práctico

Veamos un ejemplo de cómo transformar una serie eliminando la no estacionariedad utilizando R.

# Cargamos los datos de ventas de tractores
datos = read.csv("http://ucanalytics.com/blogs/wp-content/uploads/2015/06/Tractor-Sales.csv")
# Transformamos el dataframe, en una serie temporal con la función de R ts. Le indicamos la fecha de inicio y la periodicidad de los datos. En este caso como son datos mensuales, la periodicidad es 12.
serie <- ts(datos[, 2], start = c(2003,1), frequency = 12)

Representamos la serie anteriormente cargada para ver su aspecto. Se puede apreciar como no es estacionaria, puesto que tiene no tiene una media constante (tendencia ascendente).

Ahora vamos a diferenciar la serie para ver si conseguimos reducir la no estacionalidad. Para ello utilizaremos la función nsdiffs que nos dice el número de diferencias necesaria para hacer la serie estacionaria. También usaremos la función diff para ejecutarla con código.

ndif <- nsdiffs(serie)
seriedif <- diff(serie,lag = 1,differences = ndif)

Representamos y vemos como hemos eliminado la pendiente, pero por contra la varianza ahora no es constante, ya que aumente con el tiempo.

Aplicamos otra diferencia más sobre la serie anterior para ver si se corrige el problema anterior.

ndif2 <- nsdiffs(seriedif)
seriedif2 <- diff(seriedif,lag = 1,differences = ndif2)

Vemos como no hemos conseguido una serie estacionaria aplicando diferencias sobre la serie original. Por tanto, ahora vamos a comprobar si con la segunda técnica (logaritmo y diferencias) lo conseguimos.

logserie=log10(serie)
logserie_dif=diff(logserie,lag = 1,differences = ndif)

Graficamos la serie y…¡aleluya! Hemos conseguido una serie estacionaria.

Modelo ARIMA

Se trata de un modelo modelo autorregresivo integrado de promedio móvil. ARIMA es el acrónimo del ingés autoregressive integrated moving average. Seguramente no hayas entendido nada, así que vamos por partes.

Modelo ARMA

Se trata de un modelo mixto con componente autorregresiva y otra de media móvil.

En primer lugar el modelo autorregresivo se basa en la idea de que la observación actual puede explicarse con valores pasados. Estos modelos se denotan como AR(p) siendo p el número de muestras que explican la muestra actual. Es decir, un AR(1) quiere decir que sólo la muestra anterior y el ruido contribuye a la salida. Un AR(2) las dos anteriores y el ruido, etc. La expresión general del AR(p) es:

y_{t} = \phi_1 y_{t-1} + \phi_2 y_{t-2}+......+\phi_p y_{t-p} + \epsilon_t

\phi_1...\phi_p

En segundo lugar el modelo de medias móviles consiste en realizar una aproximación a la serie utilizando únicamente ruido blanco. La expresión general de los modelos MA(q):

y_t = \epsilon_t +\theta_{1} \epsilon_{t-1} +\theta_2 \epsilon_{t-2} +.....+ \theta_q \epsilon_{t-q}

\theta_{1}....\theta_{q}

Juntando ambos términos obtenemos los modelos ARMA (p,q), los cuáles determinan la serie en función su pasado hasta el retardo p, del ruido actual y el pasado del ruido hasta el retardo q. Hemos comentado anteriormente que lo que no consiguiésemos modelar de forma estacionaria sería ruido blanco. Por tanto, con este modelo la componente del ruido se modelará con un modelo MA(q).

Estos modelos ARMA(p,q) modelan una serie estacionaria, pero cuando queremos modelar directamente una no estacionaria utilizamos un modelo ARIMA.

ARIMA(p,d,q)

Hemos visto las componentes p y q, pero ¿Qué es la d? Se trata del número de diferencias necesarias para obtener la estacionariedad.

La I hace referencia al concepto Integrado. En inglés integración y diferenciación son conceptos similares, de ahí la sigla.

Una serie estacionaria que sigue modelo ARMA(p,q) visto anteriormente, se trata por tanto de una serie que sigue modelo ARIMA con I(0). I(d) es la forma de denotar el concepto de integración, siendo d el número de diferencias.

En otras palabras, una serie ARIMA(p,d,q) hay que diferenciarlo d veces para conseguir una serie ARMA(p,q).

SARIMA o ARIMA(p,d,q)(P,D,Q)

Hemos visto los parámetros p, d y q en minúsculas, pero ¿Que significan la P, D y Q en mayúsculas? Básicamente hacen referencia a un modelo ARIMA estacional. Por eso también se le conoce como SARIMA, la s de estacional (seasonal en inglés).

Como hemos comentado, es muy habitual trabajar con series temporales con períodos de medida inferiores a un año, en las que es frecuente encontrarse patrones estacionales o cíclicas.

Por tanto los términos del primer paréntesis (p,d,q) son para modelar la parte regular, es decir, la dependencia asociada a observaciones consecutivas, y los del segundo paréntesis (P,D,Q) la estacionalidad, que está asociada a observaciones separadas por n periodos.

Memoria larga

Las series o procesos con I(d) con d>0 siendo d cualquier número real se denominan series de memoria larga. Los procesos I(0) se denominan de memoria corta, por lo que los ARMA son de este tipo.

Esto quiere decir que el valor actual depende de toda su historia pasada aunque esté muy alejada en el tiempo.

Forecasting usando ARIMA

Una vez hemos definido de forma teórica los modelos ARIMA, vamos a proceder a hacer un ejemplo práctico de predicción o forecasting.

Para ello vamos a utilizar la función auto.arima de la librería forecast de R, ya que nos proporciona una opción rápida para seleccionar los parámetros p,d,q. Esta función evalúa todos los posibles modelos (model selection) considerando diferentes factores, como la estacionariedad, estacionalidad, diferencias, etc.

Cómo métricas o criterios de selección de modelos la función auto.arima utiliza el AIC y BIC. Sin entrar en detalle, buscará aquel modelo con menor AIC y BIC. Como diferencia recalcaremos que el criterio BIC penaliza más los modelos con un número mayor de parámetros estimados. Por tanto, será más óptimo utilizar el criterio BIC para modelos ARIMA.

En primer lugar dividimos la serie temporal en datos de entranmiento y de test. Como en una serie temporal no se pueden coger muestras de forma aleatoria, porque dependen del tiempo, dividimos «a mano» los datos.

train <- window(serie, start = c(2003,1), end = c(2011,12))
test <- window(serie, start = c(2012,1))
arima <- auto.arima(train, approximation = FALSE, trace = FALSE)
summary(arima)

El resultado de la función auto.arima lo mostramos en la imagen siguiente. Se puede apreciar como el mejor modelo es un ARIMA(1,1,0)(,1,1,0). Es decir, este modelo tiene q=0, sólo tiene componente AR(1) y se aplican primeras diferencias para obtener una serie estacionaria.

library(forecast)
result <- forecast(arima, h = 12*3)
# Ploteamos la predicción vs los datos de test o reales
plot(result[[4]], col = "red", lwd = 2, main = "Predicción vs real")
lines(serie, col = "blue", lwd = 2, lty = 2)

Con el código anterior hemos hecho una predicción de las ventas de tractores a 3 años vista con la función forecast de R. El parámetro h indica el número de períodos que queremos.

En la siguiente imagen mostramos nuestras predicciones en rojo, y en azul el valor real de la serie.

Esperemos que os haya gustado. Si necesitáis que ampliemos más alguno de los temas abordados en este post, dejarlo en los comentarios!