Elaborado por Jorge Rodríguez.

Introducción

En este tutorial se explica lo relacionado al pronóstico de series de tiempo. Lo cual comprende:

Para el análisis de series de tiempo en R, existen una gran variedad de librerias. En este tutorial se emplea la libreria fpp2, la cual contiene contiene:

Esta libreria se instala y se carga mediante install.packages("fpp2")y library(fpp2) respectivamente.

Elaborando la serie de tiempo en R

Cargando librerias y datos

Recordemos que se puede leer datos directamente desde una API, SQL u otros medios. Para este tutorial empleo la lectura de datos desde Excel. Entonces para ello emplearemos las siguientes funciones:

  • La función read_excel() nos permite leer datos desde un archivo excel. Para ello es necesario instalar la libreria readxl
  • La función head() nos muestra los primeros datos de nuestra tabla de datos
library(readxl)
library(fpp2)
## Registered S3 methods overwritten by 'tibble':
##   method     from  
##   format.tbl pillar
##   print.tbl  pillar
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## -- Attaching packages ---------------------------------------------- fpp2 2.4 --
## v ggplot2   3.3.3     v fma       2.4  
## v forecast  8.14      v expsmooth 2.3
## 
dataset <- read_excel("Data serie tiempo.xlsx", sheet = 'Data1')

Creando el objeto de serie de tiempo

La serie de tiempo se debe almacenar en el objeto ts empleando la función ts(). El empleo de esta función es de la siguiente forma : serie <- ts(data, frequency= ,start=). Donde el primer agumento es un vector de datos. Luego los dos siguientes argumentos se indican de acuerdo a la siguiente tabla.

Frecuencia frequency= start=
Anual 1 2000
Trimestral 4 c(2000,2)
Mensual 12 c(2000,9)
Diario 7 o 365.25 1 o c(2000,234)
Semanal 52.18 c(2000,23)
Horario 24 o 168 o 8766 1
cada 30 min 48 o 336 o 17532 1

Por ejemplo, para los datos del archivo excel adjunto :

  • En Data1, los datos son mensuales y empieza el 2004. Entonces corresponde frequency=12 y start=2004
  • En Data2, los datos son mensuales y empieza el 2005. Entonces corresponde frequency=12 y start=2005
  • En Data3, los datos son cada 30 minutos y empieza el 2019. Entonces corresponde frequency=48 y start=1

En el siguiente fragmento de código definimos el objeto de la serie de tiempo y lo almacenaremos en data_serie.

data_serie <- ts(dataset$Serie, frequency=12, start=2004)

Nota: Los resultados de esta guía corresponden a Data1 del archivo excel. Sugiero que tambien practiques con las otras series de tiempo como Data2 o Data3.

Análisis de la serie de tiempo

Para el análisis de datos, esta libreria fpp2 nos permite emplear las siguientes funciones:

Grafico de la serie de tiempo

Para el gráfico simple de la serie de tiempo empleamos la funcion autoplot(). Donde solo se requiere ingresar el objeto de la serie de tiempo.

autoplot(data_serie)+
  labs(title = "Serie de tiempo",       
       x = "Tiempo",
       y = "Valor",
       colour = "#00a0dc")+
    theme_bw()

Descomposición de la serie de tiempo

Para la descomposición de la serie de tiempo se emplea la función decompose(). Donde se debe indicar el objeto de la serie de tiempo y el tipo de descomposición. Los tipos de descomposición que acepta esta función es additive y multiplicative.

  • Aditivo: Serie=T + S + I

  • Multiplicativo: Serie=T x S x I

    Donde:

  • T: Tendencia

  • S: Estacionalidad

  • I: Irregular o error

# Descomposición de la serie de tiempo. Se almacena en el objeto fit
fit <- decompose(data_serie, type='additive')
#fit <- decompose(data_serie, type='multiplicative')

# Para graficar esta descomposición volvemos a emplear la funcion autoplot, pero con el objeto fit
autoplot(fit)+
  labs(title = "Descomposición de la serie de tiempo",                   
       x = "Tiempo",
       y = "Valor",
       colour = "Gears")+
    theme_bw()

Grafico de la serie de tiempo con su tendencia

El siguiente fragmento de código nos permite graficar la serie de tiempo con su tendencia. Notese que emplea el objeto fit en el cual guardamos previamente los valores de la descomposición. Nótese que se emplea la funcion trendcycle() para obtener los datos de tendencia del objeto fit.

autoplot(data_serie, series="Serie tiempo") + 
    autolayer(trendcycle(fit), series="Tendencia") +
    labs(title = "Serie de tiempo",      
       x = "Tiempo",
       y = "Valor"
       ) + 
    theme_bw()
## Warning: Removed 12 row(s) containing missing values (geom_path).

Grafico de estacionalidad

Para realizar el gráfico de estacionalidad empleamos la función ggseasonplot. Donde el argumento es el objeto que contiene la serie de tiempo.

ggseasonplot(data_serie)

Pronóstico

Métodos simples

Para el pronóstico de series de tiempo mediante métodos básicos, la libreria fpp2 nos brinda las siguientes funciones:

  • naive(), metodo de naive simple
  • ses(), exponential smoothing
  • meanf(), media movil
  • snaive(), metodo naive considerando estacionalidad

El argumento a colocar en estas funciones es la serie de tiempo y el valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar. Por ejemplo si deseamos pronosticar 12 datos, se debe indicar h=12.

Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:

  • fitted(), obtiene un ajuste con la data historica
  • checkresiduals(), permite analizar los residuales
# elaborando el método
m1 <- snaive(data_serie, h=96)

# graficando el pronóstico
autoplot(m1)

# verificando el ajuste del método
autoplot(m1)+autolayer(fitted(m1), series="Ajuste")
## Warning: Removed 12 row(s) containing missing values (geom_path).

# verificando los residuales
checkresiduals(m1)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 229.01, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

Método regresión

Para el pronóstico de series de tiempo mediante regresión, la libreria fpp2 nos brinda la función tslm(). Esta función la emplearemos para crear una regresión de la serie de tiempo con los datos de la descomposición estacional y/o tendencia. Entonces:

  • Si observamos solo tendencia usaremos : tslm(data_serie ~ trend)
  • Si observamos solo estacionalidad usaremos : tslm(data_serie ~ season)
  • Si observamos tendencia y estacionalidad usaremos : tslm(data_serie ~ trend + season)

Luego con la función forecast realizamos el pronostico. El argumento a colocar en estas funcion es la regresión y el valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar.

Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:

  • fitted(), obtiene pronostico con la data historica
  • checkresiduals(), permite analizar los residuales
# elaborando la regresion
regresion <- tslm(data_serie ~ trend + season)

# elaborando el pronostico
m2 <- forecast(regresion, h=96)

# graficando el pronóstico
autoplot(m2)

# verificando el ajuste del método
autoplot(m2)+autolayer(fitted(m2), series="Ajuste")

# verificando los residuales
checkresiduals(m2)

## 
##  Ljung-Box test
## 
## data:  Residuals from Linear regression model
## Q* = 939.69, df = 11, p-value < 2.2e-16
## 
## Model df: 13.   Total lags used: 24

Método holt winters

Para el pronóstico de series de tiempo mediante holt winters, la libreria fpp2 nos brinda la función hw().

Los argumentos a colocar en esta funcion son:

  • La serie de tiempo
  • El valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar.
  • El tipo de descomposición a usar para la estacionalidad. Los tipos de descomposición que acepta esta función es additive y multiplicative.

Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:

  • fitted(), obtiene pronostico con la data historica
  • checkresiduals(), permite analizar los residuales
# elaborando el pronostico
m3 <- hw(data_serie, h=96, seasonal = 'multiplicative')

# graficando el pronóstico
autoplot(m3)

# verificando el ajuste del método
autoplot(m3)+autolayer(fitted(m3), series="Ajuste")

# verificando los residuales
checkresiduals(m3)

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt-Winters' multiplicative method
## Q* = 35.398, df = 8, p-value = 2.262e-05
## 
## Model df: 16.   Total lags used: 24

ARIMA

Para el pronóstico de series de tiempo mediante ARIMA, la libreria fpp2 nos brinda la función auto.arima().

Primero crearemos un modelo ARIMA, para ello el argumento a colocar en esta funcion es la serie de tiempo. Considerar que esta función es solo una aproximación iterativa que busca los indices de AR y MA. Pues en determinadas series de tiempo podria no encontrar los indices adecuados para un modelo ARIMA. En ese caso lo adecuado es seguir la metodología de estimación de índices ARIMA. Esta metodología no esta cubierta en esta guia.

Luego con la función forecast realizamos el pronostico. El argumento a colocar en estas funcion es el modelo ARIMA y el valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar.

Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:

  • fitted(), obtiene pronostico con la data historica
  • checkresiduals(), permite analizar los residuales
# elaborando el modelo ARIMA
modelo_arima <- auto.arima(data_serie)

# elaborando el pronostico
m4 <- forecast(modelo_arima, h=96)

# graficando el pronóstico
autoplot(m4)

# verificando el ajuste del método
autoplot(m4)+autolayer(fitted(m4), series="Ajuste")

# verificando los residuales
checkresiduals(m4)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(1,1,2)[12]
## Q* = 31.738, df = 20, p-value = 0.04618
## 
## Model df: 4.   Total lags used: 24

Red neuronal

Para el pronóstico de series de tiempo mediante una red neuronal, la libreria fpp2 nos brinda la función nnetar().

Primero crearemos un modelo de red neuronal (neural network), para ello el argumento a colocar en esta funcion es la serie de tiempo.

Luego con la función forecast realizamos el pronostico. El argumento a colocar en estas funcion es el modelo de red neuronal y el valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar.

Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:

  • fitted(), obtiene pronostico con la data historica
  • checkresiduals(), permite analizar los residuales
# elaborando el modelo de red neuronal
neural_network <- nnetar(data_serie)

# elaborando el pronostico
m5 <- forecast(neural_network, h=96)

# graficando el pronóstico
autoplot(m5)

# verificando el ajuste del método
autoplot(m5)+autolayer(fitted(m5), series="Ajuste")
## Warning: Removed 12 row(s) containing missing values (geom_path).

# verificando los residuales
checkresiduals(m5)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

Estimación de error

Para estimar los errores de pronóstico, se debe realizar con los valores ocurridos o reales. Es decir este error se mide tiempo despues de haber realizado el pronóstico.

Supongamos que los valores reales ocurridos para la Data1 son los siguientes:

Mes Valor Real
Enero 19 13487 Mayo
Febrero 19 12776 Junio
Marzo 19 13812 Julio
Abril 19 13032 Agosto

Entonces, para poder comparar nuestros datos sera necesario almacenarlo en un objeto de serie de tiempo. Entonces:

real <- c(13487, 12776, 13812, 13032, 14268, 14473, 15359, 14457 )
data_real <- ts(real, frequency=12,start=2019)

La función accuracy() determina los errores de pronostico, para ello es necesario:

Entonces lo que realizará esta función es comparar el pronóstico de los siguientes 8 datos y el valor real. Pues estamos considerando que ya pasaron 8 meses, y nos encontramos en la etapa de evaluar el error de pronóstico de nuestros modelos.

# modelo en base a métodos simples
accuracy(m1,data_real)
##                  ME     RMSE      MAE      MPE     MAPE     MASE       ACF1
## Training set 384.75 505.6560 431.1548 3.523001 3.918342 1.000000  0.4327075
## Test set     365.00 466.3137 438.0000 2.602949 3.107894 1.015877 -0.2078547
##              Theil's U
## Training set        NA
## Test set     0.5388025
# modelo en base a regresion lineal
accuracy(m2,data_real)
##                         ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -5.046141e-14 332.1779 251.4463 -0.1051868 2.220769 0.5831928
## Test set     -6.430464e+02 678.8526 643.0464 -4.6198455 4.619846 1.4914515
##                    ACF1 Theil's U
## Training set  0.7182366        NA
## Test set     -0.3644790 0.7820812
# modelo en base a holt winters
accuracy(m3,data_real)
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -18.06908 230.7608 174.6523 -0.1444046 1.559487 0.4050804
## Test set     -58.72820 271.5594 148.5207 -0.4098491 1.032395 0.3444719
##                    ACF1 Theil's U
## Training set  0.1094440        NA
## Test set     -0.2064019 0.2958201
# modelo en base a ARIMA
accuracy(m4,data_real)
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  -4.256141 210.5291 157.8358 -0.0486745 1.399429 0.3660768
## Test set     -57.295339 229.0560 153.4890 -0.4244515 1.079951 0.3559951
##                    ACF1 Theil's U
## Training set -0.0261825        NA
## Test set     -0.4049233 0.2518715
# modelo en base a red neuronal
accuracy(m5,data_real)
##                        ME     RMSE      MAE         MPE     MAPE      MASE
## Training set   0.04494342 282.4865 211.8491 -0.06337753 1.885499 0.4913528
## Test set     193.23540313 397.7726 328.8728  1.30807365 2.280096 0.7627720
##                    ACF1 Theil's U
## Training set  0.2318215        NA
## Test set     -0.0533756 0.4658539

Entonces a partir de estos errores de pronóstico podemos determinar cual es el modelo mas adecuado. En la práctica se suelen evaluar varios modelos e incluso tomar como pronóstico un valor medio del resultados de dos o más modelos.