Modelos lineales de Series de Tiempo y Predicción

Una clase importante de modelos lineales de series de tiempo es la familia de los Modelos Autorregresivos Integrados de Medias Móviles (ARIMA). Propuesto por Box y Jenkis (1976). En éstos modelos se supone que el valor actual (Yt) puede depender sólo de los valores pasados de la propia serie de tiempo (Yt-1, Yt-2, ...) o de los valores pasados de algún término de error (et-1, et-2, ...). Vamos a ver cómo modelar procesos ARIMA en R.

De acuerdo con los autores del modelo ARIMA, la construcción de éste consta de tres etapas:

1. Identificar el modelo: Determinar el orden (número de términos pasados y número de errores pasados para incluir en el análisis) usando métodos gráficos o criterios de información.
2. Estimación del modelo: Los parámetros del modelo deben ser estimados usando mínimos cuadrados o máxima verosimilitud.
3. Comprobación del modelo: Examinar el modelo para diagnosticar posibles insuficiencias ( los residuos del modelo deber ser ruido blanco o de lo contrario hay dependencia lineal).

 Identificación del modelo

Vamos a analizar los precios de cierre de Facebook. Para ello usamos la librería quantmod. Si no la tienes instalada, usa el comando install.packages("quantmod").

library(quantmod) #carga la librería
getSymbols("FB", from="2012-05-18", to="2016-07-18")

Estacionariedad en media y varianza

Luego usamos transformaciones para que el proceso sea estacionario en media y en varianza. La transformación usual para lograr estacionariedad en varianza es log Pt Para lograr estacionariedad en media, realizamos una prueba de raíz unitaria. Existen varias pruebas de raíz unitaria, vamos a elegir en esta ocasión el test ADF.

Ho: Yt tiene una raíz unitaria

      Ha: Yt no tiene una raíz unitaria.

En R desarrollamos esta prueba mediante la librería "tseries". Si no la tienes instalada debes usar el comando install.packages("tseries"). Luego hacemos la prueba para una hipótesis alternativa de que la serie tiene tendencia lineal y errores estacionarios:

library(tseries)
adf.test(log(FB[,4]), alternative="explosive")

Nos arroja un mensaje diciendo que el valor-p = 0.5567. Este valor es superior para un 5% de significancia, por lo que aceptamos la hipótesis nula de que Yt tiene raíz unitaria, por lo tanto no es estacionario en media. Ahora hacemos la prueba bajo la hipótesis alternativa de que la tendencia lineal no es significativa:

adf.test(log(FB[,4]), alternative="stationary")

El valor-p = 0.4433. Aceptamos la hipótesis nula de que Yt tiene raíz unitaria. Por lo que no es estacionario. Ahora, como en ambas situaciones ha sido rechazada (Con tendencia significativa y sin ella) Debemos diferenciar la serie una vez. Recordemos que al diferenciar, el primer elemento del nuevo vector contiene NA. Por lo que no debemos tenerlo en cuenta. Luego realizamos la prueba para una hipótesis alternativa de que la serie tiene tendencia lineal y errores estacionarios:

r =  diff(log(FB[,4]))
adf.test(r[-1], alternative="explosive")

El valor-p = 0.99. Aceptamos la hipótesis nula de que Yt tiene raíz unitaria. Por lo tanto no es estacionario. Ahora hacemos la prueba bajo la hipótesis alternativa de que la tendencia lineal no es significativa:

adf.test(r[-1], alternative="stationary")

El valor-p = 0.01 es inferior al 5% de significancia. Por lo tanto rechazamos la Hipótesis nula. Los retornos logarítmicos de Facebook siguen un proceso estacionario en media. Hasta aquí logramos obtener un proceso que inicialmente no era estacionario, pero mediante diferenciación se pudo transformar en una nueva serie estacionaria en media. Ahora podemos encontrar un modelo ARMA para los retornos de Facebook.

$${ r }_{ t }=\quad \triangle log{ P }_{ t }\quad \sim \quad ARIMA(p,0,q)\quad =\quad ARMA(p,q)$$

Trazar ACF y PACF de la serie estacionaria en media y varianza

Observamos la función de autocorrelación (ACF) y función  de autocorrelación parcial (PACF):

acf(r[-1], ylim=c(-0.1,0.1)) #muestra el eje y entre -0.1 y 0.1

pacf(r[-1], ylim=c(-0.1,0.1))

Observamos en ambos casos que el rezago de orden 5 es significativo. Estimaremos el mejor ARMA.

Estimación del Modelo

Para estimar modelos ARIMA recomiendo la librería "rugarch", ya que es la librería más completa que permite realizar modelos más generales como los ARFIMAX y modelos de la familia GARCH. Para comenzar, si no tienes la librería instalada usa el comando install.packages("rugarch"). A continuación se muestra el mejor modelo estimado.

#cargamos la librería

library(rugarch)

# especificamos el modelo, indicando que no estime parámetros innecesarios mediante el argumento fixed.pars

arma.especificacion = arfimaspec(

    mean.model = list(armaOrder = c(0,5), include.mean = TRUE),

    distribution.model = "norm", 

    fixed.pars = list(ma1=0, ma2=0, ma3=0, ma4=0))

#Estimar el modelo

arma.estimar = arfimafit(arma.especificacion, r[-1])

La salida del modelo es bastante grande. En resumen, siguiendo el criterio de información de Akaike el mejor modelo estimado es un ARMA(0,5). Donde se ha verificado que sus parámetros sean significativos y que los errores sean ruido blanco. Los modelos estimados fueron los siguientes:

Mean Model      : ARFIMA(5,0,5), Akaike       -4.4712

Mean Model      : ARFIMA(0,0,5), Akaike       -4.4719

Mean Model      : ARFIMA(5,0,0), Akaike       -4.4716

La estimación se realiza por el método de máxima verosimilitud. La parte principal del cálculo de probabilidad se realiza en código C para dar mayor velocidad.

Comprobación del modelo

Como se mencionó anteriormente, los coeficientes estimados son significativos. Podemos verificarlo en la salida al escribir arma.estimar la salida es bastante grande. En resumen el coeficiente que nos interesa es ma5 que podemos ver tiene un valor-p inferior a un 5% de significancia.

          Estimate      Std. Error    t value    Pr(>|t|)
ma5    0.065788    0.031935      2.0601   0.039393

Ahora miremos si los errores distribuyen ruido blanco.

#comprobar si los errores son ruido blanco mediante un método gráfico.

acf(residuals(arma.estimar), ylim=c(-0.1,0.1))

pacf(residuals(arma.estimar), ylim=c(-0.1,0.1))



Los coeficientes de autocorrelación y autocorrelación parcial se encuentran dentro de las bandas de aceptación. En este caso los errores distribuyen ruido blanco. 

Predicción

Para predecir los siguientes 7 días (del 19 de Julio al 25 de Julio de 2016)  usamos el siguiente comando que pertenece a  la librería rugarch:

r.pron = arfimaforecast(arma.estimar, n.ahead=7)

r.pron

Horizon: 7

Roll Steps: 0

Out of Sample: 0

0-roll forecast: 

      T+1       T+2           T+3            T+4           T+5           T+6           T+7 

0.0010356 0.0002817 0.0013327 0.0007175 0.0023858 0.0010771 0.0010771 

La predicción en T+1 está basada en los datos anteriores, mientras que las predicciones posteriores están basados en la media incondicional del modelo.  Si queremos observar el gráfico de los retornos y del pronóstico de los retornos escribimos:

plot(r)

lines(r.pron, col="red")