Pronóstico de Clientes en el Aeropuerto Internacional de la Ciudad de México

Proyecto: Pronóstico de una Serie de Tiempo

 

Presentación

 

En este trabajo se desarrollará un proyecto de analítica relacionado con el modelado de una serie de tiempo. Este tipo de proyectos es muy común en prácticamente todas las empresas, pues en su Plan de Negocios anual, es donde se especifican los objetivos y los escenarios tanto operativos como financieros y éstos se construyen a partir de pronósticos.

El trabajo consistirá en pronosticar la cantidad mensual de clientes, para el año 2019, del Aeropuerto Internacional de la Ciudad de México (AICM), en otras palabras, se hará una predicción de los pasajeros que abordaron vuelos nacionales más vuelos internacionales, pero solo las salidas. El rango histórico de datos abarca desde el año 2008 a 2019 y aunque se dispone de información más actual, por cuestiones relacionadas con la pandemia de Coivd-19, se decidió hacer el ponóstico solo para 2019, ya que esto permitirá medir la precisión del modelo sin el impacto de dicha pandemia. Posteriormente, ya con el modelo final, se hará la predicción para los meses de enero y febrero de 2020, ya que en México el confinamiento empezó a finales de marzo.

Etapas

 

Como este trabajo consistará en la creación de un modelo de aprendizaje supervisado para una serie de tiempo la etapa del análisis exploratorio es un poco más densa que los análisis exploratorios descriptivos comunes. Esta etapa se constituye de las siguientes fases:

• Análisis Exploratorio
  • Análisis visual
  • Detección de outliers
  • Correlogramas
  • Descomposición de la serie
  • Análisis de estacionariedad
  • Detección de estacionalidad

Set up

 

Iniciamos configurando las opciones generales que vamos a requerir para el desarrollo de este proyecto.

knitr::opts_chunk$set(echo = TRUE,
                      warning = FALSE,
                      message = FALSE,
                      warning = FALSE,
                      fig.align = "center",
                      out.width = "90%",
                      out.height = "70%")

paquetes <- c('tidyverse',      # manipulación datos
              'lubridate',      # manejar fechas
              'knitr',          # manejo de tablas
              'kableExtra',     # manejo de tablas
              'plotly',         # gráficos interactivos
              'gridExtra',      # layout de graficos ggplot
              'nortest',        # pruebas estadísticas
              "summarytools",   # estadísticas univariantes
              "TSstudio",       # time series
              "forecast",       # time series
              "zoo",            # time series
              "anomalize",      # anomalias
              "tseries",        # serie de tiempo
              "FitAR",          # arima
              "ggfortify",      # diagnosis residuales
              "lmtest",         # test coeficientes arima
              "strucchange",    # gráfico CUSUM test
              'prophet',        # modelo prophet
              'tidymodels')     # manejo de modelos

instalados <- paquetes %in% installed.packages()

if(sum(instalados == FALSE) > 0) {
  install.packages(paquetes[!instalados])
}

lapply(paquetes, require, character.only = TRUE)

 

Carga de datos

 

La información para este análisis es pública y los datos así como su descripción se encuentran en la página del AICM.

Es importante aclarar que los datos originales se encuentran dentro de informes en formato pdf, por lo que, para extraerlos se recurrió a la técnica conocida como Reconocimiento Óptico de Caracteres (OCR). Para el periodo 2012-2019, se utilizó R y para el resto del periodo 2008-2011 se utilizó Python, ya que ofreció los mejores resultados para esos años. Los archivos con todo el código, de R y Python, se pueden encontrar en mi repositorio de github. Habiendo hecho esta aclaración, a continuación, ya solo se importa el archivo final:

salidas <- read.csv("Salidas_Nacionales.csv")

head(salidas)

glimpse(salidas)

## Rows: 144
## Columns: 2
## $ Fecha   <chr> "01/01/2008", "01/02/2008", "01/03/2008", "01/04/2008", "01/05~
## $ Salidas <int> 1034682, 957075, 1132697, 1055504, 1087835, 1089415, 1257439, ~

En este primer resumen podemos ver que la variable Fecha está como character, por lo que, es necesario convertirla al formato correcto:

salidas$Fecha <- dmy(salidas$Fecha)
glimpse(salidas)

## Rows: 144
## Columns: 2
## $ Fecha   <date> 2008-01-01, 2008-02-01, 2008-03-01, 2008-04-01, 2008-05-01, 2~
## $ Salidas <int> 1034682, 957075, 1132697, 1055504, 1087835, 1089415, 1257439, ~

Ya tenemos la variable Fecha en el formato correcto, sin embargo, la mayoría de las funciones relacionadas con el análisis de una serie temporal requieren que los datos estén en formato de serie de tiempo (ts), por lo que, haremos la conversion:

salidas.ts <- ts(salidas[,2], start = c(2008, 1), end = c(2019, 12), frequency = 12)

ts_info(salidas.ts)

##  The salidas.ts series is a ts object with 1 variable and 144 observations
##  Frequency: 12 
##  Start time: 2008 1 
##  End time: 2019 12

Ahora si, ya tenemos nuestra serie de tiempo la cual inicia en enero de 2008 y termina en diciembre de 2019. Pasemos ahora al análisis exploratorio.

 

Análisis exploratorio

Análisis visual

Empecemos por analizar visualmente la serie para conocer qué forma tiene así como sus principales medidas estadísticas:

salidas %>% 
  ggplot(aes(Fecha, Salidas)) +
  geom_line() +
  scale_x_date(date_breaks = "6 month", date_labels = "%Y-%b") +
  theme(axis.text.x = element_text(angle = 90, size = 8)) +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "AICM: Salidas Nacionales e Internacionales 2008-2019")

Lo que se puede apreciar es que hasta finales de 2010 la tendencia es constante y hasta un poco decreciente, sin embargo, desde inicios de 2011 la tendencia es creciente y aditiva. Esta serie claramente muestra un cambio de tendencia, lo que muy probablemente ocasione que no se pueda cumplir el supuesto de estacionariedad incluso si aplicamos diferencias. Conviene de una vez comprobar ese cambio o quiebre, para lo cual haremos lo siguiente:

1. Crear el gráfico CUSUM, que nos mostrará visualmente si la serie presenta un quiebre y en qué parte.
2. Crear una variable dicotómica que tomará valores de 0 hasta enero de 2011 y posterior a esa fecha tomara valores de 1.
3. Hacer una regresión lineal cuyas variables independientes serán: la variable dicotómica creada, el tiempo (que será un consecutivo) y la interacción entre el tiempo y la variable dicotómica.

t <- c(1:length(salidas.ts))

p.salidas <- efp(salidas.ts ~ t, type = "OLS-CUSUM")

plot(p.salidas)

La gráfica nos indica que la serie sí presenta un quiebre en su tendencia al parecer desde 2009 hasta inicios de 2011. Ahora, lo comprobaremos formalmente con una regresión.

t.salidas <- salidas %>% 
  mutate(D = if_else(Fecha < "2011-02-01", 0, 1))

reg <- lm(Salidas ~ D + t + (D*t), data = t.salidas)

summary(reg)

## 
## Call:
## lm(formula = Salidas ~ D + t + (D * t), data = t.salidas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -338916  -55131   -2537   68120  306220 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1092503      40545  26.946  < 2e-16 ***
## D            -465389      54447  -8.548 1.94e-14 ***
## t              -3135       1860  -1.685   0.0942 .  
## D:t            13865       1898   7.304 1.93e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 120800 on 140 degrees of freedom
## Multiple R-squared:  0.9103, Adjusted R-squared:  0.9084 
## F-statistic: 473.7 on 3 and 140 DF,  p-value: < 2.2e-16

Los coeficientes de la variable dicotómica y de su interacción con el tiempo salieron significativos. Ahora, el coeficeinte de la variable dicotómica D, es el intercepto diferencial y el coeficiente de la interacción entre el tiempo y la variable dicotómica es el coeficiente de la pendiente diferencial. Lo anterior, significa que los dos periodos son distintos y las causas de la diferencia son tanto el intercepto como la pendiente.

Por lo anterior, la serie que vamos a modelar presenta cambio de tendencia a partir de 2011, lo cual se explica en parte debido a que en 2008-2009 se presentó la crisis financiera y, además, en 2009, en México, inició la crisis sanitaria de la Influenza H1N1, también, a inicios de 2011 empezaron los problemas financieros de la empresa Mexianana de Aviación, la cual terminó por cerrar sus operaciones a finales de agosto de ese mismo año. Todos esos acontecimientos afectaron la economía lo que explica la tendencia a la baja en el periodo 2008-2011.

Además, se puede sospechar que hay presencia de estacionalidad por el patrón que sigue la serie, también, hay que considerar que la data representa a las personas que viajan por avión, por lo que, el sentido común nos indicaría que en periodos vacacionales la cantidad de personas que salen de viaje aumentará. Más adelante lo tendremos que comprobar formalmente.

Datos atípicos (outliers)

 

Algo que afecta mucho a los modelos de series de tiempo son los llamados outliers o datos atípicos, ya que al ser observaciones demasiado grandes o pequeñas afectan la capacidad predictiva de los modelos. Por lo tanto, es de vital importancia identificarlos y tratarlos. Un primer análisis visual se hace con apoyo del boxplot:

salidas %>% 
  ggplot(aes(x = factor(1), y = Salidas)) +
  geom_boxplot(outlier.colour = "red", fill = "grey", alpha = 0.3) +
  stat_summary(fun = mean, 
               geom = "point",
               shape = 20, 
               size = 4, 
               color = "red",
               fill = "red") +
  scale_y_continuous(labels = scales::comma) +
  labs(
    title = "Boxplot Salidas",
    x = "Salidas"
  )

En el boxplot se aprecia que la serie no presenta observaciones atípicas, además, se observa que la media (punto rojo) y la mediana no están muy separadas lo que nos indica que no hay asimetria fuerte. Complementemos la gráfica con algunos estadísiticos básicos:

descr(salidas$Salidas,
      justify = "c",
      headings = T) %>% 
  kbl(align = "c", digits = 2 ,format.args = list(big.mark = ",")) %>%
  kable_paper("hover", full_width = F)

	Salidas
Mean	1,456,935.94
Std.Dev	399,210.23
Min	700,290.00
Q1	1,104,285.00
Median	1,373,431.00
Q3	1,788,912.00
Max	2,339,595.00
MAD	451,024.71
IQR	679,670.00
CV	0.27
Skewness	0.34
SE.Skewness	0.20
Kurtosis	-1.08
N.Valid	144.00
Pct.Valid	100.00

Pasemos a otro método más robusto para detectar si realmente en la serie hay observaciones atípicas. Para esto nos apoyaremos del paquete anomalize que nos permitirá calibrar el intervalo de confianza y el porcentaje máximo de posibles valores atípicos a mostrar, entre otras cosas. Algo importante a mencionar es que este paquete requiere que el input esté en formato de tibble o tibble_time, por lo que, no usaremos el objeto tipo serie de tiempo creado al inicio del análisis.

salidas.tbl <- salidas

salidas.tbl$Fecha <- as.yearmon(salidas.tbl$Fecha)
salidas.tbl$Salidas <- as.numeric(salidas.tbl$Salidas)

salidas.tbl <- as_tibble(salidas.tbl)

prep_tbl_time(salidas.tbl, message = TRUE)

salidas.tbl %>% 
  time_decompose(Salidas, method = "stl") %>% 
  anomalize(remainder, method = "gesd", max_anoms = 0.2, alpha = 0.03) %>% 
  time_recompose() %>% 
  plot_anomalies(time_recomposed = TRUE, 
                 alpha_dots = 0.4,
                 alpha_ribbon = 0.1,
                 fill_ribbon = "grey80") + 
  scale_y_continuous(labels = scales::comma) + 
  labs(title = "Detectando Anomalías", y = "Salidas")

En la gráfica anterior podemos ver que se detectaron algunas observaciones atípicas, sin embargo, al parecer solo uno de ellos está muy alejado del patrón de la serie. Tendremos que tratarlo para que no nos cause problemas más adelante.

Para identificar qué observaciones son atípicas usaremos la función tsoutliers del paquete forecast, ya que permite hacer la identificación y, además, nos propondrá los valores a sustituir.

tsoutliers(salidas.ts) %>% 
  as.data.frame() %>%
  rename(Observacion = index, Reemplazo = replacements) %>% 
  kbl(align = "c",digits = 0 ,format.args = list(big.mark = ",")) %>%
  kable_paper("hover", full_width = F)

Observacion	Reemplazo
14	878,178
17	1,039,968

Vemos que la función identifica las observaciones en las posiciones 14 y 17 como atípicas, y nos muestra los valores por los cuales haría la sustitución (que utiliza el método de interpolación lineal). El trabajo lo hace la función tsclean que no solo sustituye los valores atípicos sino que también identifica valores faltantes (que no es este el caso).

Veamos gráficamente cuáles serán los puntos que se sustituirán:

# nuevos puntos
autoplot(tsclean(ts(salidas$Salidas)), series = "clean", colour = 'red', lwd = 0.9) +
  autolayer(ts(salidas$Salidas), series = "original", color = 'gray', lwd = 1) +
  geom_point(data = tsoutliers(salidas.ts) %>% as.data.frame(), 
             aes(x = index, y = replacements), col = 'blue') +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Sustitución Outliers: nuevos puntos", x = "Tiempo", y = "Salidas")

Los puntos en color azul serán las nuevas observaciones que reemplazarán a los outliers que se identificaron. Se hace la sustitución:

salidas.ts <- tsclean(salidas.ts)

Se hace el reemplazo en la base original para volver a ver la gráfica de las anomalías:

# se actualizan los datos en el dataframe original
salidas$Salidas[14] <- 878178.1
salidas$Salidas[17] <- 1039967.7

# se repite el proceso para ver la grafica

salidas.tbl <- salidas

salidas.tbl$Fecha <- as.yearmon(salidas.tbl$Fecha)
salidas.tbl$Salidas <- as.numeric(salidas.tbl$Salidas)

salidas.tbl <- as_tibble(salidas.tbl)

#prep_tbl_time(salidas.tbl, message = TRUE)

salidas.tbl %>% 
  time_decompose(Salidas, method = "stl") %>% 
  anomalize(remainder, method = "gesd", max_anoms = 0.2, alpha = 0.03) %>% 
  time_recompose() %>% 
  plot_anomalies(time_recomposed = TRUE, 
                 alpha_dots = 0.4,
                 alpha_ribbon = 0.1,
                 fill_ribbon = "grey80") +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Detectando Anomalías", y = "Salidas")

Se siguen detectando datos atípicos, sin embargo, no están muy alejados del intervalo de confianza, por lo que, nos quedaremos con esta serie.

Algo importante a mencionar es que se debe tener bien claro el objetivo del análisis que se esté buscando, ya que en algunos casos podría ser incorrecto sustituir las observaciones atípicas.

Correlogramas

 

Para todos los modelos de series de tiempo es de suma importancia revisar las funciones de autocorrelación tanto simple como parcial con el propósito de empezar a conocer cómo está la estructura de los datos. Los correlogramas en palabras simples nos muestran la memoria de la serie, o sea, qué tanto depende el dato actual de lo que sucedió en el pasado. veamos los correlogramas de la serie con 48 rezagos, que son un tercio de la longitud de la serie (lo recomendado):

ts_cor(salidas.ts, lag.max = 48)

 

Al observar el correlograma ACF podemos sospechar que la serie no es estacionaria, o sea, la media no es constante a lo largo del tiempo, porque los valores de la función de autocorrelación no decaen rápidamente conforme aumentan los rezagos en el tiempo. Además, al parecer también hay presencia de estacionalidad, ya que en PACF hay coeficientes significativos, por lo que, tenemos que descomponer la serie para ver cada componente por separado.

Descomposición de la serie

 

Vamos a generar las gráficas de los componentes de esta serie: tendencia, estacionalidad y residuo.

ts_decompose(salidas.ts, type = "additive")

Lo que se puede apreciar en las gráficas es que la tendencia es aditiva, ya que la serie va en aumento con el tiempo y las oscilaciones no se van haciendo más grandes sino que se mantienen del mismo tamaño. Lo anterior, tenemos que confirmarlo de manera un poco más formal.

Hay presencia de estacionalidad muy marcada, ya que el mismo patrón se repite año tras año. También, hay que confirmarlo e identificar en qué meses se presenta.

Vamos a identificar formalmente el tipo de tendencia que tiene la serie, para lo cual se necesita hacer lo siguiente:

• calcular el logaritmo de la mediana para cada periodo.
• calcular el logaritmo de la diferencia entre el percentil 80 y el percentil 20 para cada periodo.
• comparar ambas series en un gráfico de dispersión.
• calcular el coeficiente de correlación.
• probar la significancia estadistica de dicho coeficiente.

tabla <- salidas %>% 
  mutate(
    Año = year(Fecha),
    Mes = month(Fecha)
  ) %>% 
  select(-Fecha) %>% 
  group_by(Año) %>% 
  summarise(ln_Mediana = log(median(Salidas)),
         ln_Rango = log((quantile(Salidas, 0.8) - quantile(Salidas, 0.2))))
  
tabla %>% 
  ggplot(aes(ln_Mediana, ln_Rango)) +
  geom_point() +
  labs(title = "Niveles vs Dispersión",
       x = "Nivel",
       y = "Dispersión")

cor(tabla$ln_Mediana, tabla$ln_Rango)

## [1] 0.8175311

cor.test(tabla$ln_Mediana, tabla$ln_Rango)

## 
##  Pearson's product-moment correlation
## 
## data:  tabla$ln_Mediana and tabla$ln_Rango
## t = 4.4892, df = 10, p-value = 0.001162
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4589703 0.9470796
## sample estimates:
##       cor 
## 0.8175311

Al ver la gráfica creada podemos ver que existe cierta dependencia y tendencia positiva entre el nivel y la dispersión, lo que nos sugiere que la tendencia es de tipo multiplicativa, sin embargo, al volver a la gráfica en niveles se observa que en los primeros años la serie es constante y a partir de 2011 se presenta tendencia creciente. Por lo anterior, los primeros años están jalando o disfrasando el verdadero tipo de tendencia de la serie. Además, hay que recordar que ya comprobamos que la serie presenta un quiebre de tendencia, lo cual afecta esta gráfica, ya que las primeras observaciones que corresponden al periodo 2008-2011 la relación es positiva y para el periodo 2012-2019 parece que no hay relación.

Análisis de estacionariedad

 

En esta parte, vamos a analizar formalmente la estacionariedad de la serie. En series de tiempo, es necesario que la serie a modelar sea estable, o sea, tenga una media y varianza constantes, sobre todo si el modelo a usar para realizar el pronóstico es un arima. Antes del análisis hay que tener presente que la serie al presentar quiebre de tendencia afectará cualquier prueba de estacionariedad que se le haga, por lo tanto, es de esperar que, aunque se apliquen diferencias, puede que no pase la prueba elegida.

Volvamos a mirar los correlogramas:

ggAcf(salidas.ts, lag.max = 48)

ggPacf(salidas.ts, lag.max = 48)

En el ACF vemos que los coeficientes decaen lentamente, lo cual es un indicativo de que la serie no es estacionaria, sin embargo, para confirmarlo formalmente haremos el test de raíz unitaria de Augmented Dickey-Fuller:

adf.test(salidas.ts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  salidas.ts
## Dickey-Fuller = -3.984, Lag order = 5, p-value = 0.01203
## alternative hypothesis: stationary

El p-valor de la prueba (0.01203) es menor al nivel de significancia de 0.05, por lo que, la prueba nos indica que la serie es estacionaria. Contrario a lo que nos mostraba el correlograma, o sea, al no tener raíz unitaria la serie no presenta un patrón sistemático impredecible. Sin embargo, vamos a determinar con apoyo de un gráfico y una regresión si debemos hacer alguna transformación a los datos, ya que la gráfica en niveles de la serie sugiere estabilizar la serie. Será la gráfica conocida como rango/media, que consiste en calcular la media de cada periodo (año) y la varianza. Posteriormente, se hace una regresión con ambas variables y si es significativo el coeficiente asociado con la media entonces hay que aplicar alguna transformación a la serie:

tabla2 <- salidas %>% 
  mutate(
    Año = year(Fecha),
    Mes = month(Fecha)
  ) %>% 
  select(-Fecha) %>% 
  group_by(Año) %>% 
  summarise(ln_Media = log(mean(Salidas)),
         ln_varianza = log(var(Salidas)))

ggplot(tabla2, aes(x = ln_Media, y = ln_varianza)) +
  geom_point() +
  labs(title = "Gráfico rango-media")

coeftest(lm(ln_varianza ~ ln_Media, data = tabla2))

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  5.99895    5.25816  1.1409 0.280504   
## ln_Media     1.24018    0.37127  3.3404 0.007486 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En el gráfico se aprecia que los primeros años están jalando la tendencia de la nube de puntos, además, el resultado de la regresión nos indica que debemos hacer alguna transformación a la serie, ya que el coeficiente asociado con la media resultó significativo.

Detección de estacionalidad

 

Ahora, pasemos al análisis de la estacionalidad empezando por el análisis visual:

salidas %>% 
  mutate(
    Año = year(Fecha),
    Mes = month(Fecha)
  ) %>% 
  select(-Fecha) %>% 
  group_by(Año, Mes) %>% 
  ggplot(aes(x = factor(Mes), y = Salidas)) +
  geom_boxplot() +
  scale_y_continuous(labels = scales::comma) +
  stat_summary(fun = mean, 
               geom = "point",
               shape = 20, 
               size = 4, 
               color = "red",
               fill = "red") +
  labs(title = "Boxplot Meses",
       x = "Meses")

ggsubseriesplot(salidas.ts) +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Salidas Promedio", x = "Meses", y = "Salidas")

Al ver el boxplot y la gráfica de promedios mensuales parece que si hay diferencias entre los meses en cada año y, por lo tanto, presencia de un patrón estacional. Sigamos con el análisis de la serie:

ggseasonplot(salidas.ts) + 
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Salidas por Mes", y = "Salidas", x = "Mes")

Las gráficas vistas hasta aquí nos indican que es en el mes de julio donde la afluencia de pasajeros es mayor, lo cual coincide con el periodo vacacional de verano. También, se aprecia que el segundo mes más alto de salidas es diciembre como era de esperarse por la temporada navideña.

Sigamos viendo algunas vistas más para confirmar el patrón estacional existente:

ts_seasonal(salidas.ts, type = "all", title = "Gráficos con la Tendencia")

El patrón estacional es el mismo a lo largo del tiempo: alta demanda en el mes de julio y diciembre y baja demanda en febrero y septiembre.

Para tener una visión más clara del patrón estacional vamos a quitar el efecto de la tendencia, ya que la serie crece año tras año desde el año 2011:

ts_seasonal(salidas.ts - decompose(salidas.ts)$trend, 
            type = "all", 
            title = "Gráfico sin Tendencia")

Sin el efecto de la tendencia se aprecia más claramente el patrón estacional. Podemos identificar que julio, agosto y diciembre son los meses de mayores salidas o vuelos, por el contrario, enero, febrero y septiembre los meses de menor demanda o salidas de pasajeros.

Otra vista diferente de la serie es mediante un heatmap:

ts_heatmap(salidas.ts, title = "Salidas")

Se aprecia que la demanda de pasajeros que sale de viaje ya sea al interior del país o al extranjero, desde la Ciudad de México, ha estado creciendo de manera sostenida desde los últimos 8 años. Además, se confirma que julio es el mes con mayor demanda seguido por el mes de diciembre.

Ahora, el siguiente paso es identificar el nivel de correlación entre la serie y sus rezagos:

ts_cor(salidas.ts, lag.max = 48)

Como era de esperarse sí hay una correlación entre la serie y sus rezagos estacionales. El retraso estacional (o retraso 12) de la serie tiene una fuerte relación lineal con la serie. Lo anterior, lo podemos ver más claramente al graficar estos componentes:

ts_lags(salidas.ts, lags = c(12, 24, 36, 48))

Solo nos queda identificar formalmente la estacionalidad. Para ello vamos a crear variables dummies, una por cada mes para, posteriormente, hacer una regresión sin incluir la constante para evitar la colinealidad. Para que el efecto estacional este presente los parámetros estimados del modelo deben ser estadísticamente distintos de cero, individualmente.

tabla3 <- salidas %>% 
  mutate(
    Mes = month(Fecha)
  ) %>% 
  mutate(
    Enero = ifelse(Mes == 1, 1,0),
    Febrero = ifelse(Mes == 2, 1,0),
    Marzo = ifelse(Mes == 3, 1,0),
    Abril = ifelse(Mes == 4, 1,0),
    Mayo = ifelse(Mes == 5, 1,0),
    Junio = ifelse(Mes == 6, 1,0),
    Julio = ifelse(Mes == 7, 1,0),
    Agosto = ifelse(Mes == 8, 1,0),
    Septiembre = ifelse(Mes == 9, 1,0),
    Octubre = ifelse(Mes == 10, 1,0),
    Noviembre = ifelse(Mes == 11, 1,0),
    Diciembre = ifelse(Mes == 12, 1,0)
  ) %>% 
  dplyr::select(-Fecha, -Mes)

modelo <- lm(Salidas ~ 0 +., data = tabla3)
summary(modelo)

## 
## Call:
## lm(formula = Salidas ~ 0 + ., data = tabla3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -550852 -351625 -101963  345302  726060 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## Enero       1354725     113425   11.94   <2e-16 ***
## Febrero     1206017     113425   10.63   <2e-16 ***
## Marzo       1454283     113425   12.82   <2e-16 ***
## Abril       1411900     113425   12.45   <2e-16 ***
## Mayo        1464396     113425   12.91   <2e-16 ***
## Junio       1443255     113425   12.72   <2e-16 ***
## Julio       1696737     113425   14.96   <2e-16 ***
## Agosto      1564176     113425   13.79   <2e-16 ***
## Septiembre  1355666     113425   11.95   <2e-16 ***
## Octubre     1465781     113425   12.92   <2e-16 ***
## Noviembre   1479518     113425   13.04   <2e-16 ***
## Diciembre   1605587     113425   14.16   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 392900 on 132 degrees of freedom
## Multiple R-squared:  0.938,  Adjusted R-squared:  0.9324 
## F-statistic: 166.5 on 12 and 132 DF,  p-value: < 2.2e-16

Los coeficientes resultaron significativos, por lo que, confirmamos la presencia de estacionalidad en la serie.

Ya para terminar esta sección, calcularemos los índices de variación estacional para darnos una idea del impacto en cada mes. Para lo cual, haremos lo siguiente:

• se descompone la serie.
• se divide la serie original entre su tendencia.
• se promedian los periodos.
• se ajustan los indices para que sumen el total de periodos (vía ponderaciones).

dc <- decompose(salidas.ts, type = "additive")

salidas %>% 
  mutate(
    serie = dc$x,
    tendencia = dc$trend,
    Año = year(Fecha),
    Mes = month(Fecha),
    v1 = serie/tendencia
  ) %>% 
  group_by(Mes) %>% 
  summarise(ie = mean(v1, na.rm = T)) %>% 
  kbl(align = "c", digits = 4) %>% kable_paper("hover", full_width = F)

Mes	ie
1	0.9537
2	0.8408
3	1.0149
4	0.9802
5	1.0114
6	0.9887
7	1.1712
8	1.0690
9	0.9147
10	0.9878
11	0.9898
12	1.0758

Los resultados nos indican que en el mes de julio la demanda es mayor al promedio anual un 17%, aproximadamente, por el contrario, en el mes de febrero la demanda de salidas es 16% menor al promedio anual, en otras palabras, en los meses de febrero, en promedio, solo se llega al 84% del promedio anual.

Transformaciones

 

Vimos que la serie en niveles tiene tendencia creciente y de tipo aditiva, porque los picos en cada año no aumentan. Por lo anterior, se tiene que aplicar alguna transformación y la más común son las diferencias.

Hay que identificar cuántas diferenciaciones aplicar a la serie. Se puede hacer aplicando 1 diferencia y revisar el correlograma hasta que veamos la serie sin tendencia, sin embargo, optaremos por un camino más directo con apoyo de la función ndiffs del paquete forecast que nos indicará cuántas diferencias aplicar a la serie:

Número de diferenciaciones a la parte regular de la serie: 1

Ya que confirmamos que hay estacionalidad en la serie tenemos que determinar cuántas diferencias aplicar. Nos apoyaremos de la función nsdiffs del paquete forecast que nos indica el orden de diferenciación (esto tambien se puede hacer manualmente diferenciando la serie en los periodos estacionales y revisando el correlograma):

Número de diferenciaciones a la parte regular de la serie: 1

Por lo tanto, hay que diferenciar 1 vez en la parte regular y 1 vez a la parte estacional de la serie en niveles. Esto se especificará en el algoritmo del modelo.

Conjunto de entrenamiento y prueba

 

Antes de generar el primer modelo, es necesario dividir la serie en una parte para entrenar el modelo y otra parte para validarlo. Cuando se trabaja con series de tiempo dicha partición no se hace de manera aleatoria como para cualquier otro modelo de machine learning, ya que las series de tiempo van ordenadas de manera cronológica, por lo tanto, el conjunto de test siempre será la parte final de la serie. Consideraremos el año 2019 para probar el modelo y el resto para entrenar el modelo.

salidas.ts.train <- window(salidas.ts, start = c(2008,1), end = c(2018,12))
salidas.ts.test <- window(salidas.ts, start = c(2019,1), end = c(2019,12))

Creación del modelo base

 

Cuando el objetivo de un proyecto sea el pronosticar una serie de tiempo, es altamente recomendable hacer un primer modelo que será la base a partir del cual se tratará de mejorarlo. Hay muchas métricas para evaluar este tipo de modelos, la mayoría relacionadas con los errores, o sea, la diferencia entre el valor real y el valor pronosticado. Ahora, lo que se debe evaluar es el dataset de test, o sea, la parte del conjunto de datos que el modelo no ha visto. En el presente análisis usaremos el Mean Absolute Percentage Error (MAPE), ya que es más fácil de interpretar.

Nuestro primer modelo será el modelo básico de suavizado exponencial ets, que está formado por: el componente de error, el componente de tendencia y el componente estacional.

ETS: ExponenTial Smoothing

fit.ets <- ets(salidas.ts.train,
               model = "ZZZ",
               damped = T,
               additive.only = T,
               lambda = "auto"
               )


fit.ets

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = salidas.ts.train, model = "ZZZ", damped = T, additive.only = T,  
## 
##  Call:
##      lambda = "auto") 
## 
##   Box-Cox transformation: lambda= 0.3839 
## 
##   Smoothing parameters:
##     alpha = 0.6174 
##     beta  = 0.0331 
##     gamma = 1e-04 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 546.3433 
##     b = -1.119 
##     s = 18.0384 -1.8706 -2.3477 -19.5167 16.0081 37.7301
##            -2.2611 2.3417 -4.531 4.3537 -37.3082 -10.6369
## 
##   sigma:  6.3844
## 
##      AIC     AICc      BIC 
## 1151.750 1157.803 1203.641

El resultado del modelo: (A,Ad,A), nos indica que modeló los tres componentes de forma aditiva. Veamos ahora los residuales, para identificar si existe autocorrelación:

checkresiduals(fit.ets, lag = 26)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 45.171, df = 9, p-value = 8.579e-07
## 
## Model df: 17.   Total lags used: 26

El p-valor de la prueba nos indica que los errores no son independientes, por lo que, el pasado está influyendo en lo que sucede en el presente. Habrá que corregir esto más adelante, por lo pronto, se sigue anializando este primer modelo.

qqnorm(fit.ets$residuals)
qqline(fit.ets$residuals)

jarque.bera.test(fit.ets$residuals)

## 
##  Jarque Bera Test
## 
## data:  fit.ets$residuals
## X-squared = 25.011, df = 2, p-value = 3.706e-06

El análisis de la normalidad muestra que los errores no se distribuyen de forma normal. Sin embargo, como es el modelo base sigamos adelante con las predicciones para los siguientes 12 meses.

fc.ets <- forecast(fit.ets, h = 12)


fc.ets %>% as.data.frame() %>% 
  kbl(align = "c", format.args = list(big.mark = ",")) %>% 
  kable_paper("hover", full_width = F)

	Point Forecast	Lo 80	Hi 80	Lo 95	Hi 95
Jan 2019	1,936,463	1,875,897	1,998,220	1,844,313	2,031,397
Feb 2019	1,749,042	1,681,361	1,818,376	1,646,196	1,855,754
Mar 2019	2,062,943	1,976,726	2,151,438	1,932,001	2,199,215
Apr 2019	2,000,518	1,905,382	2,098,526	1,856,172	2,151,582
May 2019	2,059,264	1,951,980	2,170,106	1,896,613	2,230,238
Jun 2019	2,029,460	1,913,046	2,150,139	1,853,127	2,215,770
Jul 2019	2,357,858	2,219,125	2,501,811	2,147,772	2,580,152
Aug 2019	2,184,675	2,042,108	2,333,217	1,969,024	2,414,299
Sep 2019	1,915,273	1,774,587	2,062,630	1,702,772	2,143,376
Oct 2019	2,050,692	1,894,167	2,214,944	1,814,388	2,305,069
Nov 2019	2,059,622	1,893,005	2,234,983	1,808,285	2,331,409
Dec 2019	2,223,165	2,038,292	2,418,025	1,944,401	2,525,288

forecast::accuracy(fc.ets, salidas.ts.test)

                ME     RMSE      MAE      MPE     MAPE      MASE

Training set 4523.606 34262.94 25741.90 0.319050 1.981202 0.2285163 Test set 45827.168 61529.56 48871.01 2.174499 2.304601 0.4338384 ACF1 Theil’s U Training set -0.04758451 NA Test set 0.34334424 0.3529856

autoplot(fc.ets) + scale_y_continuous(labels = scales::comma) +
  labs(title = "Pronóstico Modelo Base: ets(A,Ad,A)", y = "Salidas")

Este primer modelo presenta un MAPE de 2.3% para el conjunto de test, no es muy elevado, sin embargo, no pasó los supuestos de independencia de los errores ni normalidad. Veamos cómo se ve la gráfica de las predicciones con la serie original.

test_forecast(forecast.obj = fc.ets, 
              actual = salidas.ts, 
              test = salidas.ts.test) %>% 
  layout(legend = list(x = 0.1, y = 0.95))

Ajustemos el modelo ets para mejorarlo diferenciando 1 vez la serie:

salidas.ts.diff <- diff(salidas.ts)

salidas.ts.train.diff <- window(salidas.ts.diff, start = c(2008,2), end = c(2018,12))
salidas.ts.test.diff <- window(salidas.ts.diff, start = c(2019,1), end = c(2019,12))

fit.ets.diff <- ets(salidas.ts.train.diff,
               model = "ZZZ",
               damped = T
               )

checkresiduals(fit.ets.diff)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 33.525, df = 7, p-value = 2.113e-05
## 
## Model df: 17.   Total lags used: 24

qqnorm(fit.ets.diff$residuals)
qqline(fit.ets.diff$residuals)

jarque.bera.test(fit.ets.diff$residuals)

## 
##  Jarque Bera Test
## 
## data:  fit.ets.diff$residuals
## X-squared = 0.86664, df = 2, p-value = 0.6484

En el correlograma de los residuales podemos ver que al inicio solo un coeficiente sale de las bandas, sin embargo, el test no lo pasa. En cuanto a la prueba de normalidad vemos que ya la pasa al tener un p-valor mayor al 0.05. Nos quedamos con este modelo, por lo que, ahora se pronosticarán los próximos 12 meses.

fc.ets.diff <- forecast(fit.ets.diff, h = 12)

last_salidas <- salidas.ts.train[132]
pred <- fc.ets.diff$mean

de_diff_salidas_fore <- cumsum(pred) + last_salidas

ddiff_sfore_ts <- ts(de_diff_salidas_fore, start = 133)

ddiff_sfore_ts <- ts(ddiff_sfore_ts, start = c(2019, 1), frequency = 12)

ddiff_sfore_ts %>% as.data.frame() %>% 
  rename(Pronostico = x) %>%  mutate(Fecha = salidas$Fecha[133:144]) %>% 
  select(Fecha, Pronostico) %>% 
  kbl(align = "c", format.args = list(big.mark = ",")) %>% 
  kable_paper("hover", full_width = F)

Fecha	Pronostico
2019-01-01	1,966,643
2019-02-01	1,819,524
2019-03-01	2,072,125
2019-04-01	2,031,114
2019-05-01	2,084,730
2019-06-01	2,069,841
2019-07-01	2,335,097
2019-08-01	2,208,419
2019-09-01	2,010,270
2019-10-01	2,119,293
2019-11-01	2,136,397
2019-12-01	2,269,743

forecast::accuracy(ddiff_sfore_ts, salidas.ts.test)

##                ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set 4642.229 45350.73 32506.31 0.1398487 1.588626 0.5558924 0.2636444

ts.plot(salidas.ts, ddiff_sfore_ts,
        lty = c(1,3), col = c("steelblue", "red"), 
        main = "Pronóstico modelo en 1ras Diferencias: ets(A,Ad,A)", 
        ylab = "Salidas")

El modelo con diferencias mejoró, ya que se reduce el error de test de 2.3 a 1.58.

Ahora, optimizaremos los parámetros.

Este algoritmo tiene tres parámetros que están relacionados con el suavisado (alpha), la tendencia (beta) y la estacionalidad (gamma). Por tanto, vamos a optimizarlos para encontrar el valor que nos arroje el modelo con el menor error y una vez encontrado el valor del componente volveremos a reajustar el modelo.

Alpha: componente de suavisado

# identificar el alpha óptimo
alpha <- seq(.01, .99, by = .001)

MAPE <- NA

for(i in seq_along(alpha)) {
  fit <- ets(salidas.ts.train, 
             alpha = alpha[i],
             model = "ZZZ",
             damped = T,
             additive.only = T,
             lambda = "auto")
  
  fc.fit <- forecast(fit, h = 12)

  MAPE[i] <- forecast::accuracy(fc.fit, salidas.ts.test)[2,5]
}

# convertir a data frame e identificar el valor del alpha mínimo
alpha.fit <- tibble(alpha, MAPE)
alpha.min <- filter(alpha.fit, MAPE == min(MAPE))
alpha.min

# gráfico MAPE vs. alpha
ggplot(alpha.fit, aes(alpha, MAPE)) +
  geom_line() +
  geom_point(data = alpha.min, aes(alpha, MAPE), size = 2, color = "blue") +
  labs(title = "Mean Absolute Percentage Error")

Beta: componente de tendencia

# identificar el beta óptimo
beta <- seq(.01, .99, by = .001)

MAPE <- NA

for(i in seq_along(beta)) {
  fit <- ets(salidas.ts.train, 
             beta = beta[i],
             model = "ZZZ",
             damped = T,
             additive.only = T,
             lambda = "auto")
  
  fc.fit <- forecast(fit, h = 12)
  
  MAPE[i] <- forecast::accuracy(fc.fit, salidas.ts.test)[2,5]
}

# convertir a data frame e identificar el valor del beta mínimo
beta.fit <- tibble(beta, MAPE)
beta.min <- filter(beta.fit, MAPE == min(MAPE))
beta.min

# gráfico MAPE vs. beta
ggplot(beta.fit, aes(beta, MAPE)) +
  geom_line() +
  geom_point(data = beta.min, aes(beta, MAPE), size = 2, color = "blue") +
  labs(title = "Mean Absolute Percentage Error")

Gamma: componente de estacionalidad

# identificar el gamma óptimo
gamma <- seq(.01, .99, by = .001)

MAPE <- NA

for(i in seq_along(gamma)) {
  fit <- ets(salidas.ts.train, 
             gamma = gamma[i],
             model = "ZZZ",
             damped = T,
             additive.only = T,
             lambda = "auto")
  
  fc.fit <- forecast(fit, h = 12)
  
  MAPE[i] <- forecast::accuracy(fc.fit, salidas.ts.test)[2,5]
}

# convertir a data frame e identificar el valor del gamma mínimo
gamma.fit <- tibble(gamma, MAPE)
gamma.min <- filter(gamma.fit, MAPE == min(MAPE))
gamma.min

# gráfico MAPE vs. gamma
ggplot(gamma.fit, aes(gamma, MAPE)) +
  geom_line() +
  geom_point(data = gamma.min, aes(gamma, MAPE), size = 2, color = "blue") +
  labs(title = "Mean Absolute Percentage Error")

Al optimizar los tres componentes el que generó el menor error fue alpha con 1.53. Se reajusta del modelo con ese parámetro:

# reajuste del modelo con alpha optimo
fit.ets <- ets(salidas.ts.train,
               model = "ZZZ",
               damped = T,
               additive.only = T,
               lambda = "auto",
               alpha = alpha.min$alpha
               )

fit.ets

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = salidas.ts.train, model = "ZZZ", damped = T, alpha = alpha.min$alpha,  
## 
##  Call:
##      additive.only = T, lambda = "auto") 
## 
##   Box-Cox transformation: lambda= 0.3839 
## 
##   Smoothing parameters:
##     alpha = 0.249 
##     beta  = 0.0294 
##     gamma = 1e-04 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 546.569 
##     b = -0.9689 
##     s = 17.7407 -1.8181 -2.1958 -19.3507 15.9644 37.4089
##            -2.1358 2.2996 -4.6524 4.2632 -36.9145 -10.6095
## 
##   sigma:  7.1591
## 
##      AIC     AICc      BIC 
## 1179.983 1185.352 1228.991

checkresiduals(fit.ets, lag = 26)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 109.41, df = 9, p-value < 2.2e-16
## 
## Model df: 17.   Total lags used: 26

qqnorm(fit.ets$residuals)
qqline(fit.ets$residuals)

jarque.bera.test(fit.ets$residuals)

## 
##  Jarque Bera Test
## 
## data:  fit.ets$residuals
## X-squared = 11.782, df = 2, p-value = 0.002764

fc.ets <- forecast(fit.ets, h = 12)

forecast::accuracy(fc.ets, salidas.ts.test)

##                     ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set  6165.468 37909.77 30182.91 0.4347513 2.339992 0.2679402 0.3831377
## Test set     20412.619 45755.96 32768.79 0.9654602 1.537096 0.2908955 0.3492182
##              Theil's U
## Training set        NA
## Test set     0.2639074

test_forecast(forecast.obj = fc.ets, 
              actual = salidas.ts, 
              test = salidas.ts.test) %>% 
  layout(legend = list(x = 0.1, y = 0.95))

autoplot(fc.ets) + 
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Pronóstico ETS Salidas 2019", y = "Salidas")

ARIMA: AutoRegressive Integrated Moving Average

 

Ahora se optará por un modelo arima para ver si podemos mejorar el error de test y, además, que sus residuales pasen las pruebas de normalidad y autocorrelación. Se comenzará por la generación automática del modelo con apoyo de la función auto.arima del paquete forecast:

fitarima1 <- auto.arima(salidas.ts.train)
fitarima1

## Series: salidas.ts.train 
## ARIMA(0,1,1)(0,1,2)[12] 
## 
## Coefficients:
##           ma1     sma1    sma2
##       -0.3665  -0.8231  0.2908
## s.e.   0.0890   0.1126  0.1504
## 
## sigma^2 estimated as 1.744e+09:  log likelihood=-1437.71
## AIC=2883.42   AICc=2883.78   BIC=2894.54

coeftest(fitarima1)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1  -0.366467   0.089045 -4.1155 3.863e-05 ***
## sma1 -0.823140   0.112562 -7.3128 2.617e-13 ***
## sma2  0.290750   0.150392  1.9333    0.0532 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fitarima1)

## Series: salidas.ts.train 
## ARIMA(0,1,1)(0,1,2)[12] 
## 
## Coefficients:
##           ma1     sma1    sma2
##       -0.3665  -0.8231  0.2908
## s.e.   0.0890   0.1126  0.1504
## 
## sigma^2 estimated as 1.744e+09:  log likelihood=-1437.71
## AIC=2883.42   AICc=2883.78   BIC=2894.54
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 2775.622 39145.38 28595.76 0.1475052 2.184726 0.2538507 0.01535602

En este primer modelo los tres coeficientes son significativos. Sin embargo, hay que hacer la diagnosis a los residuales:

Autocorrelación:

checkresiduals(fitarima1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,2)[12]
## Q* = 23.542, df = 21, p-value = 0.3158
## 
## Model df: 3.   Total lags used: 24

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

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

ggtsdiag(fitarima1)

La gráfica de los residuales no muestra tendencia lo cual es bueno y en los correlogramas solo un coeficiente sale de los límites. El p-value de la prueba Ljung-Box es mayor al 0.05, lo cual nos indica que los residuales no están correlacionados. Además, los p-valores para la prueba Q de Ljung-Box están por encima de 0,05, lo que indica “no significativo”.

Normalidad:

qqnorm(fitarima1$residuals)
qqline(fitarima1$residuals)

jarque.bera.test(fitarima1$residuals)

## 
##  Jarque Bera Test
## 
## data:  fitarima1$residuals
## X-squared = 3.6914, df = 2, p-value = 0.1579

El p-value es mayor al 0.05, por lo tanto, podemos afirmar que los residuales de este modelo siguen una distribución normal.

Este modelo pasa las pruebas de autocorrelación y normalidad, por lo que, podemos darlo como válido y es apto para hacer pronósticos:

Predicción y Evaluación:

fc.fitarima1 <- forecast(fitarima1, h = 12)
forecast::accuracy(fc.fitarima1, salidas.ts.test)

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  2775.622 39145.38 28595.76 0.1475052 2.184726 0.2538507
## Test set     16397.681 38113.62 23194.48 0.7660115 1.089580 0.2059023
##                    ACF1 Theil's U
## Training set 0.01535602        NA
## Test set     0.47835744 0.2214405

El error de test es de 1.09, que es menor al mejor modelo ets que se generó anteriormente. Veamos los resultados gráficamente:

test_forecast(forecast.obj = fc.fitarima1, actual = salidas.ts, test = salidas.ts.test) %>% 
  layout(legend = list(x = 0.1, y = 0.95))

autoplot(fc.fitarima1) + 
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Pronóstico auto.arima1 Salidas 2019", y = "Salidas")

Este será nuestro modelo arima base, ya que pasa las pruebas de autocorrelación y normalidad de sus residuales.

Veamos si podemos mejorarlo modificando algunos de sus parámetros:

fitarima2 <- auto.arima(salidas.ts.train, 
                        trace = T,
                        stepwise = T,
                        approximation = T,
                        test = "adf",
                        allowdrift = T,
                        optim.method ="BFGS"
                        )

## 
##  Fitting models using approximations to speed things up...
## 
##  ARIMA(2,1,2)(1,1,1)[12]                    : Inf
##  ARIMA(0,1,0)(0,1,0)[12]                    : 2661.918
##  ARIMA(1,1,0)(1,1,0)[12]                    : 2613.67
##  ARIMA(0,1,1)(0,1,1)[12]                    : 2603.277
##  ARIMA(0,1,1)(0,1,0)[12]                    : 2644.167
##  ARIMA(0,1,1)(1,1,1)[12]                    : 2603.142
##  ARIMA(0,1,1)(1,1,0)[12]                    : 2611.132
##  ARIMA(0,1,1)(2,1,1)[12]                    : 2577.624
##  ARIMA(0,1,1)(2,1,0)[12]                    : 2576.204
##  ARIMA(0,1,0)(2,1,0)[12]                    : 2586.307
##  ARIMA(1,1,1)(2,1,0)[12]                    : Inf
##  ARIMA(0,1,2)(2,1,0)[12]                    : 2578.272
##  ARIMA(1,1,0)(2,1,0)[12]                    : 2576.742
##  ARIMA(1,1,2)(2,1,0)[12]                    : Inf
## 
##  Now re-fitting the best model(s) without approximations...
## 
##  ARIMA(0,1,1)(2,1,0)[12]                    : 2875.321
## 
##  Best model: ARIMA(0,1,1)(2,1,0)[12]

coeftest(fitarima2)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1  -0.356835   0.089151 -4.0026 6.266e-05 ***
## sar1 -0.724047   0.083730 -8.6474 < 2.2e-16 ***
## sar2 -0.526351   0.087282 -6.0305 1.635e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fitarima2)

## Series: salidas.ts.train 
## ARIMA(0,1,1)(2,1,0)[12] 
## 
## Coefficients:
##           ma1     sar1     sar2
##       -0.3568  -0.7240  -0.5264
## s.e.   0.0892   0.0837   0.0873
## 
## sigma^2 estimated as 1.591e+09:  log likelihood=-1433.48
## AIC=2874.97   AICc=2875.32   BIC=2886.09
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 3023.755 37392.03 27177.97 0.1593727 2.080588 0.2412646 0.01732966

En esta corrida el mejor modelo fue un ARIMA(0,1,1)(2,1,0), el cual tiene todos los coeficientes significativos. Analicemos sus residuales.

Autocorrelación:

checkresiduals(fitarima2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(2,1,0)[12]
## Q* = 15.748, df = 21, p-value = 0.7836
## 
## Model df: 3.   Total lags used: 24

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

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

ggtsdiag(fitarima2)

El p-value de la prueba Ljung-Box es mayor al 0.05, por lo tanto, los residuales no están correlacionados, además, en su gráfica no se aprecia ningún patrón. Todos los coeficiente están dentro de las bandas en los correlogramas. Los p-valores para la prueba Q de Ljung-Box están por encima de 0,05, lo que indica “no significativo”.

Normalidad:

qqnorm(fitarima2$residuals)
qqline(fitarima2$residuals)

jarque.bera.test(fitarima2$residuals)

## 
##  Jarque Bera Test
## 
## data:  fitarima2$residuals
## X-squared = 5.3535, df = 2, p-value = 0.06879

Los residuales no presentan autocorrelación y se distribuyen normalmente. por tanto, pasamos a hacer las predicciones para conocer el error de test:

Predicción y Evaluación:

fc.fitarima2 <- forecast(fitarima2, h = 12)
forecast::accuracy(fc.fitarima2, salidas.ts.test)

##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 3023.755 37392.03 27177.97 0.1593727 2.080588 0.2412646 0.01732966
## Test set     7730.809 45468.44 31736.19 0.3480207 1.495838 0.2817289 0.51810977
##              Theil's U
## Training set        NA
## Test set     0.2636668

Arroja un error de 1.49, ligeramente mayor que el modelo anterior (1.08)

veamos las gráficas:

test_forecast(forecast.obj = fc.fitarima2, actual = salidas.ts, test = salidas.ts.test) %>% 
  layout(legend = list(x = 0.1, y = 0.95))

autoplot(fc.fitarima2) + 
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Pronóstico auto.arima2 Salidas 2019", y = "Salidas")

Hagamos un último modelo arima, pero ahora usaremos la función arima, para construirlo por partes:

Primero recordemos la forma de los correlogramas, puesto que nos ayudarán a definir la estructura de modelo:

tsdisplay(salidas.ts.train)

Empecemos por aplicar una diferencia y ver de nuevo los correlogramas:

tsdisplay(diff(salidas.ts.train))

Ahora con 12 rezagos para ver la estructura estacional de la serie:

tsdisplay(diff(salidas.ts.train, lag = 12))

Con la información proporcionada por los correlogramas se creará un modelo arima que para la parte regular sera (2,1,2) y para la parte estacional (0,1,1):

fitarima3 <-  arima(log(salidas.ts.train), 
                  order = c(2,1,2), 
                  seasonal = list(order = c(0,1,1),
                                  period = 12),
                  method = "ML")


coeftest(fitarima3)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1   0.631598   0.082481   7.6575 1.896e-14 ***
## ar2  -0.661005   0.075036  -8.8092 < 2.2e-16 ***
## ma1  -1.068908   0.040476 -26.4085 < 2.2e-16 ***
## ma2   0.999969   0.049311  20.2789 < 2.2e-16 ***
## sma1 -0.844861   0.122300  -6.9081 4.912e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fitarima3)

## 
## Call:
## arima(x = log(salidas.ts.train), order = c(2, 1, 2), seasonal = list(order = c(0, 
##     1, 1), period = 12), method = "ML")
## 
## Coefficients:
##          ar1     ar2      ma1     ma2     sma1
##       0.6316  -0.661  -1.0689  1.0000  -0.8449
## s.e.  0.0825   0.075   0.0405  0.0493   0.1223
## 
## sigma^2 estimated as 0.0008767:  log likelihood = 238.08,  aic = -464.15
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.002089955 0.02844674 0.02053924 0.01471124 0.1458633 0.2158159
##                     ACF1
## Training set 0.003126655

Todos los coeficiente de este modelo son significativos. Ahora se analizan los residuales.

Autocorrelación:

checkresiduals(fitarima3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,2)(0,1,1)[12]
## Q* = 28.18, df = 19, p-value = 0.08003
## 
## Model df: 5.   Total lags used: 24

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

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

ggtsdiag(fitarima3)

Este modelo pasa la prueba de Ljung-Box al tener un p-valor mayor al 0.05, además, la gráfica de los residuales no muestran ningún patrón o tendencia. En los correlogramas solo un coeficiente sale de la banda superior. Los p-valores para la prueba Q de Ljung-Box están por encima de 0,05, lo que indica “no significativo”.

Normalidad:

qqnorm(fitarima3$residuals)
qqline(fitarima3$residuals)

jarque.bera.test(fitarima3$residuals)

## 
##  Jarque Bera Test
## 
## data:  fitarima3$residuals
## X-squared = 58.79, df = 2, p-value = 1.713e-13

Los residuales no presentan autocorrelación pero no se distribuyen normalmente. Pasamos a hacer las predicciones para conocer el error de test:

Predicción y Evaluación:

fc.fitarima3 <- forecast(fitarima3, h = 12)

# se invierten los logaritmos
fc.fitarima3$mean <- exp(fc.fitarima3$mean)
fc.fitarima3$fitted <- exp(fc.fitarima3$fitted)
fc.fitarima3$x <- exp(fc.fitarima3$x)

forecast::accuracy(fc.fitarima3, salidas.ts.test)

##                     ME     RMSE      MAE      MPE     MAPE      MASE
## Training set  2282.009 38588.06 28080.81 0.168258 2.057594 0.2492793
## Test set     41325.189 76980.19 59010.25 2.012655 2.784381 0.5238465
##                     ACF1 Theil's U
## Training set -0.04429377        NA
## Test set      0.35523522 0.4472975

Arroja un error de 2.78, mayor que el mejor modelo arima (1.08)

veamos las gráficas:

test_forecast(forecast.obj = fc.fitarima3, actual = salidas.ts, test = salidas.ts.test) %>% 
  layout(legend = list(x = 0.1, y = 0.95))

fc.fitarima3$lower[,1:2] <- exp(fc.fitarima3$lower[,1:2])
fc.fitarima3$upper[,1:2] <- exp(fc.fitarima3$upper[,1:2])

autoplot(fc.fitarima3) + 
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Pronóstico arima3 Salidas 2019", y = "Salidas")

Tenemos tres modelos arima los cuales todos pasaron la prueba de autocorrelación en los residuales y solo el modelo arima 3 no pasa la prueba de normalidad. El modelo que presenta el menor error es el 1 con un MAPE de 1.08. Descartaremos el modelo 3 ya que sus residuales no se distribuyen normalmente.

De estos dos modelos elegiremos el modelo 1, porque generó el menor error de test, además, ambos modelos tienen 3 coeficientes, por lo que, no podemos basarnos en el principio de parsimonia (el más simple).

TSLM: Time Series Linear Model

Se generará el modelo tslm que está diseñado para ajustar modelos lineales con datos de series de tiempo y en el cual se pueden incluir variables predictoras, en particular la tendencia y la estacionalidad.

salidas.lm  <- tslm(salidas.ts.train ~  trend + season, lambda = "auto")

summary(salidas.lm)

## 
## Call:
## tslm(formula = salidas.ts.train ~ trend + season, lambda = "auto")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.666  -9.218  -2.374   7.687  47.286 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 489.00480    5.40367  90.495  < 2e-16 ***
## trend         1.37152    0.03729  36.782  < 2e-16 ***
## season2     -26.73830    6.93206  -3.857 0.000187 ***
## season3      13.96798    6.93236   2.015 0.046170 *  
## season4       4.40207    6.93286   0.635 0.526675    
## season5      10.97165    6.93356   1.582 0.116213    
## season6       6.50020    6.93446   0.937 0.350464    
## season7      45.70716    6.93556   6.590 1.26e-09 ***
## season8      23.90810    6.93687   3.447 0.000785 ***
## season9     -11.49589    6.93837  -1.657 0.100183    
## season10      5.29325    6.94007   0.763 0.447147    
## season11      5.73127    6.94198   0.826 0.410685    
## season12     24.57878    6.94408   3.540 0.000573 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.26 on 119 degrees of freedom
## Multiple R-squared:  0.9285, Adjusted R-squared:  0.9213 
## F-statistic: 128.8 on 12 and 119 DF,  p-value: < 2.2e-16

salidas.fcast.lm <- forecast(salidas.lm, h = 12)

test_forecast(forecast.obj = salidas.fcast.lm, actual = salidas.ts, test = salidas.ts.test) %>% 
  layout(legend = list(x = 0.1, y = 0.95))

salidas.fcast.lm$x <- salidas.ts.train
forecast::accuracy(salidas.fcast.lm, salidas.ts.test)

##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  2087.944 82526.70 67371.24 -0.3447535 5.501881 0.5980688
## Test set     42027.132 59494.84 45237.35  2.0010454 2.138258 0.4015816
##                   ACF1 Theil's U
## Training set 0.8867416        NA
## Test set     0.3842807 0.3403756

autoplot(salidas.fcast.lm) + labs(title = "Pronóstico Modelo tslm1 Salidas 2019")

Este primer modelo generó un MAPE de 2.13. Haremos un segundo modelo pero ahora incluyendo la estacionalidad con variables dicotómicas:

# Usando variables dummies para modelar la estacionalidad
meses <- seasonaldummy(salidas.ts.train)

salidas.lm2  <- tslm(salidas.ts.train ~ meses + trend + season)

tsdisplay(residuals(salidas.lm2))

salidas.fcast2 <- forecast(salidas.lm2,
    data.frame(meses = I(seasonaldummy(salidas.ts.train, 12))))

test_forecast(forecast.obj = salidas.fcast2, actual = salidas.ts, test = salidas.ts.test) %>% 
  layout(legend = list(x = 0.1, y = 0.95))

autoplot(salidas.fcast2) + labs(title = "Pronóstico modelo tslm2 Salidas 2019")

forecast::accuracy(salidas.fcast2, salidas.ts.test)

##                         ME     RMSE      MAE        MPE     MAPE     MASE
## Training set -2.468953e-11 99422.98 81301.20 -0.4033314 6.569292 0.721728
## Test set      8.681211e+04 98467.53 88712.75  4.0379108 4.146165 0.787522
##                   ACF1 Theil's U
## Training set 0.8921495        NA
## Test set     0.5237033 0.5636143

El segundo modelo tslm no mejoró el error de test, ya que este tipo de algoritmo no es adecuado para el tipo de datos que tenemos. Por tanto, probaremos un último modelo para ver si podemos mejorar el error del modelo arima1, que es hasta el momento el que ha generado el menor error de test.

PROPHET

Será el modelo Prophet,creado por Facebook para hacer sus pronósticos comerciales, a ver qué resultados nos genera.

Antes de iniciar con el modelo, es preciso aclarar que este modelo necesita que las columnas donde están las fechas y la variable a pronosticar tengan determinados nombres, por lo que, antes que nada, se procede a renombrar dichas columnas:

# Se verifica el tipo Fecha para el tiempo y se renombran las columnas: "ds" y "y"
df <- salidas

cleaned_df <- df %>% 
  rename(ds = Fecha,
         y = Salidas) %>% 
  mutate(ds = as.Date(ds, "%d/%m/%Y"))

cleaned_df %>% head()

Se establece la partición de los datos en train y test:

correct_split <- initial_time_split(cleaned_df %>% arrange(ds), prop = 1-(12/144))
correct_split

## <Analysis/Assess/Total>
## <132/12/144>

cleaned_train <- cleaned_df %>%
  filter(ds <"2019-01-01")

cleaned_test <- cleaned_df %>%
  filter(ds >="2019-01-01")

# se visualiza la división
bind_rows(
  training(correct_split) %>% mutate(type = "train"),
  testing(correct_split) %>% mutate(type = "test")
) %>% 
  ggplot(aes(x = ds, y= y, color = type, group = NA)) + 
  geom_line() + 
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Split", y = "Salidas")

Se crea el modelo y se hace la predicción:

m <- prophet(df = cleaned_train %>% arrange(ds), yearly.seasonality = T)

future <- make_future_dataframe(m, periods = 12, freq = "month", include_history = TRUE)

forecast <- predict(m, future)

Se grafica la serie con el pronostico:

plot(m, forecast) + 
  add_changepoints_to_plot(m) +
  labs(title = "Modelo Prophet: Pronóstico a 12 meses")

Veamos otra vista de las predicciones, en cuántos meses quedó por debajo y por arriba del valor real:

# se valida el set de test
s <- forecast %>% 
  inner_join(cleaned_test, by = "ds") %>% 
  mutate(
    dif = (y-yhat), 
    esPositivo = dif >= 0)

s %>% 
  ggplot(aes(x = ds, y = dif, fill = esPositivo)) + 
  geom_col(position = "identity") +
  scale_fill_manual(values = c("firebrick4", "dodgerblue4")) +
  scale_y_continuous(labels = scales::comma) +
  theme_minimal() + 
  labs(title = "Diferencia Entre Salidas Esperadas y Actuales",
          subtitle = "Las Salidas pronosticadas quedaron por debajo durante la mayor parte del año", y = "Diferencia", x = "Año")

En el siguiente gráfico de dispersión, se compara el ajuste de las observaciones pronosticadas contra las observaciones reales. Se aprecia como, al igual que en los demás modelos, el pronóstico del mes de febrero quedó muy por debajo del valor real. En los demás meses los pronósticos se ajustaron relativamente bien.

qplot(s$y, s$yhat, log = "xy") +
  geom_abline() +
  scale_y_continuous(labels = scales::comma) +
  scale_x_continuous(labels = scales::comma) +
  labs(title = "Comparando las Observaciones Esperadas vs Reales", 
       y = "Salidas Esperadas",
       x = "Salidas Reales") +
  theme_minimal()

Se revisan los errores:

mape <- function(y, yhat) {
  return(mean(abs(y - yhat)/ y)*100)
}

eval <- cleaned_test %>% 
  left_join(forecast, by = "ds") %>% 
  select(ds, y, yhat, yhat_upper, yhat_lower)

mape.prophet1 <- mape(eval$y, eval$yhat)

cat("MAPE del modelo Prophet 1: ",`mape.prophet1`)

MAPE del modelo Prophet 1: 2.181772

Se grafica el pronóstico:

plot(m, forecast) +
  add_changepoints_to_plot(m) +
  geom_point(cleaned_test, mapping = aes(as.POSIXct(ds), y), col = "red") +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Pronóstico Modelo Prophet 1", y = "Salidas", x = "Año")

Vemos que este modelo también subestima el mes de febrero, por lo que, vamos a probar la técnica de validación cruzada para crear modelos aditivos y multiplicativos para ver cuáles generan el menor error de generalización:

Se creann los k-folds:

timeseries_k_folds <- sliding_period(cleaned_df %>% arrange(ds), 
                                     ds,
                                     period = "year",
                                     lookback = Inf,
                                     assess_stop = 1)

Con apoyo de una función se crearán los diferentes tipos de modelos:

tune_prophet <- function(splits){
  train_data <- analysis(splits)
  test_data <- assessment(splits)
  
  m1 <- prophet(df = train_data, seasonality.mode = "additive")
  m2 <- prophet(df = train_data, seasonality.mode = "multiplicative")
  
  future <- make_future_dataframe(m1, 
                                  periods = nrow(test_data), 
                                  freq = "month",
                                  include_history = FALSE)
  
  
  bind_rows(
    predict(m1, future) %>% select(ds, yhat) %>% mutate(type = "additive"),
    predict(m2, future) %>% select(ds, yhat) %>% mutate(type = "multiplicative")
  ) %>% left_join(test_data, by = "ds")
}

Se muestran los forecast, su tipo y métrica de error:

# se hacen los k-folds
ts_tune <- timeseries_k_folds %>% 
  slice_tail(n = 10) %>%
  mutate(res = map(splits, tune_prophet))

# errores
ts_tune %>% 
  select(id, res) %>% 
  unnest(res) %>% 
  group_by(id, type) %>% 
  arrange(ds) %>% 
  mutate(forecast = paste0("forecast_", row_number())) %>% 
  ungroup() %>% 
  select(forecast, type, ds, yhat, y) %>% 
  group_by(forecast, type) %>% 
  yardstick::mape(truth = y, estimate = yhat) %>% 
  ungroup() %>% 
  group_by(forecast) %>% 
  slice_min(.estimate) %>% 
  ungroup()

Veamos la mediana de los errores (MASE) en los modelos aditivos y multiplicativos, pero comparando los pronósticos con un modelo ingenuo (naive):

ts_tune %>% 
  select(id, res) %>% 
  unnest(res) %>% 
  group_by(id, type) %>% 
  arrange(ds) %>% 
  mutate(forecast = paste0("forecast_", row_number())) %>% 
  ungroup() %>% 
  select(forecast, type, ds, yhat, y) %>% 
  left_join(cleaned_df %>% 
              mutate(naive = lag(y, n = 1, order_by = ds)) %>% 
              drop_na() %>% 
              select(ds, naive),
            by =  "ds") %>% 
  group_by(forecast, type) %>% 
  summarise(mase = mean(abs(yhat - y)) / mean(abs(naive - y))) %>% 
  ungroup() %>% 
  group_by(type) %>% 
  summarise(median = median(mase))

ts_tune %>% 
  select(id, res) %>% 
  unnest(res) %>% 
  group_by(id, type) %>% 
  arrange(ds) %>% 
  mutate(forecast = paste0("forecast_", row_number())) %>% 
  ungroup() %>% 
  select(id, forecast, type, ds, yhat, y) %>% 
  group_by(id, type) %>% 
  mutate(naive_date = min(ds) - months(1)) %>% 
  ungroup() %>% 
  left_join(cleaned_df %>% 
              mutate(naive = lag(y, n = 1, order_by = ds)) %>% 
              drop_na() %>% 
              select(ds, naive),
            by =  c("naive_date" = "ds")) %>% 
  group_by(forecast, type) %>% 
  summarise(mase = mean(abs(yhat - y)) / mean(abs(naive - y))) %>% 
  ungroup() %>% 
  group_by(type) %>% 
  summarise(median = median(mase))

En ambos casos el menor error se presenta en el tipo multiplicativo, por lo que, crearemos un segundo modelo pero ahora con el tipo de estacionalidad multiplicativa:

m2 <- prophet(df = cleaned_train %>% arrange(ds),
              yearly.seasonality = T, 
              seasonality.mode = "multiplicative")

future2 <- make_future_dataframe(m2, periods = 12, freq = "month", include_history = TRUE)

forecast2 <- predict(m2, future2)

Se grafica la serie con el pronostico:

plot(m2, forecast2) + 
  add_changepoints_to_plot(m2) +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Modelo Prophet 2: Pronóstico a 12 meses", y = "Salidas", x = "Año")

Volvamos a ver cómo quedo el pronóstico, o sea, en cuántos meses quedó por debajo y por arriba del valor real:

# se valida el set de test
s2 <- forecast2 %>% 
  inner_join(cleaned_test, by = "ds") %>% 
  mutate(
    dif = (y-yhat), 
    esPositivo = dif >= 0)

s2 %>% 
  ggplot(aes(x = ds, y = dif, fill = esPositivo)) + 
  geom_col(position = "identity") +
  scale_fill_manual(values = c("firebrick4", "dodgerblue4")) +
  scale_y_continuous(labels = scales::comma) +
  theme_minimal() + 
  labs(title = "Diferencia Entre Salidas Esperadas y Actuales",
          subtitle = "Las Salidas pronosticadas quedaron por debajo durante la mayor parte del año", y = "Diferencia", x = "Año")

En la gráfica anterior, ya se aprecian los datos para los primeros meses más cercanos a las observaciones reales, sin embargo, este modelo subestimó el mes de julio. Veamos la gráfica de dispersión de las observaciones esperadas y reales:

qplot(s2$y, s2$yhat, log = "xy") +
  geom_abline() +
  scale_y_continuous(labels = scales::comma) +
  scale_x_continuous(labels = scales::comma) +
  labs(title = "Comparando las Observaciones Esperadas vs Reales", 
       y = "Salidas Esperadas",
       x = "Salidas Reales") +
  theme_minimal()

Se aprecia que las predicciones se ajustan mejor a los datos reales.

Se revisan los errores:

mape <- function(y, yhat) {
  return(mean(abs(y - yhat)/ y)*100)
}

eval2 <- cleaned_test %>% 
  left_join(forecast2, by = "ds") %>% 
  select(ds, y, yhat, yhat_upper, yhat_lower) 

mape.prophet.f <- mape(eval2$y, eval2$yhat)

cat("MAPE del modelo Prophet: ",`mape.prophet.f`)

## MAPE del modelo Prophet:  1.88846

Se grafica el pronóstico:

plot(m2, forecast2) +
  add_changepoints_to_plot(m2) +
  geom_point(cleaned_test, mapping = aes(as.POSIXct(ds), y), col = "red") +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Pronóstico Modelo Prophet Multiplicativo", y = "Salidas", x = "Año")

Este modelo se ve y se ajusta mejor, además, el error de test bajó a 1.88, por lo que, nos quedamos con este segundo modelo multiplicativo.

Ensamble

Son tres los modelos que generaron un bajo error MAPE en el conjunto de prueba, por lo que, los vamos a combinar esperando obtener un modelo más robusto. Para hacer lo anterior, crearemos una función que haga el ensamble de los siguientes modelos:

• ets con el alpha óptimo
• arima 1
• prophet multiplicativo

Se crea la función:

ensamble <- function(df, n_pronosticos) {
  
  N <- nrow(df)
  n <- n_pronosticos
  fin <- N - n
  
  inicio.train <- df$Fecha[fin]
  inicio.test <- df$Fecha[fin]+1
  
  df_train <- df[1:fin,]
  df_train <- ts(df_train[,2], start = inicio.train, frequency = 12)
  
  df_test <- df[(fin + 1):N,] %>% ts(start = inicio.test, frequency = 12)
  df_test <- ts(df_test[,2], start = inicio.test, frequency = 12)
  
  h <- n
  
  fit.ets <- ets(df_train,
                model = "ZZZ",
                damped = T,
                additive.only = T,
                lambda = "auto",
                alpha = alpha.min$alpha
                )
  prediccion.ets <- forecast(fit.ets, h = h)$mean %>% as.vector()
  
  
  arima <- auto.arima(df_train)
  prediccion.arima <- forecast(arima, h = h)$mean %>% as.vector()
  
  
  
  cleaned_train <- df[1:fin,] %>% 
    rename(ds = Fecha,
         y = Salidas) %>% 
  mutate(ds = as.Date(ds, "%d/%m/%Y"))

  cleaned_test <- df[(fin + 1):N,] %>% 
    rename(ds = Fecha,
         y = Salidas) %>% 
  mutate(ds = as.Date(ds, "%d/%m/%Y"))
  
  m.prophet <- prophet(df = cleaned_train %>% arrange(ds),
              yearly.seasonality = T, 
              seasonality.mode = "multiplicative")

  future <- make_future_dataframe(m.prophet, 
                                   periods = h,
                                   freq = "month",
                                   include_history = TRUE)

  fc.prophet <- predict(m.prophet, future)

  prediccion.prophet <- cleaned_test %>% 
  left_join(fc.prophet, by = "ds") %>% 
  select(yhat) %>% as.vector()
  
  
  
  prediccion.ensamble <- (prediccion.ets + prediccion.arima + prediccion.prophet)/3
  
  return(data_frame(Fecha = cleaned_test$ds,
                    Actual = df_test %>% as.vector(), 
                    Prediccion = prediccion.ensamble %>% as.vector() %>% round(0), 
                    Error = (Prediccion - Actual) %>% as.vector(),
                    Error_Porcentual = round((Error / Actual)*100, 2)))
} 

Veamos los resultados del ensamble:

# Usamos la función: 
df_ensamble <- ensamble(salidas, 12)

df_ensamble$Prediccion <- df_ensamble$Prediccion$yhat
df_ensamble$Error <- df_ensamble$Error$yhat


df_ensamble %>% kbl(align = "c",
                    digits = 2 , 
                    format.args = list(big.mark = ",")) %>%
  kable_paper("hover", full_width = F)

Fecha	Actual	Prediccion	Error	Error_Porcentual
2019-01-01	1,956,707	1,952,676	-4,031	-0.21
2019-02-01	1,755,724	1,763,365	7,641	0.44
2019-03-01	2,072,594	2,080,001	7,407	0.36
2019-04-01	2,091,546	2,034,959	-56,587	-2.71
2019-05-01	2,190,456	2,096,812	-93,644	-4.28
2019-06-01	2,121,607	2,066,573	-55,034	-2.59
2019-07-01	2,339,595	2,397,076	57,481	2.46
2019-08-01	2,207,163	2,231,714	24,551	1.11
2019-09-01	1,962,079	1,954,227	-7,852	-0.40
2019-10-01	2,102,570	2,102,308	-262	-0.01
2019-11-01	2,124,365	2,123,315	-1,050	-0.05
2019-12-01	2,254,497	2,275,157	20,660	0.92

Los pronósticos se acercan a los datos reales de test, que corresponden al año 2019, sin embargo, el mes de mayo quedó subestimado y, por el contrario, el mes de julio quedó sobreestimado esto es debido a la influencia del modelo Prophet. Ahora, comparemos los errores de los tres modelos con el del ensamble para ver si valio la pena aplicar esta técnica.

# Se crea una function que calcule algunas medidas de evaluación: 
eval_medidas <- function(df_modelo) {
  
  act <- df_modelo %>% pull(Actual)
  pred <- df_modelo %>% pull(Prediccion)
  err <- act - pred 
  per_err <- abs(err / act)
  
  # Mean Absolute Error (MAE): 
  mae <- err %>% abs() %>% mean()
  
  # Mean Squared Error (MSE): 
  mse <- mean(err^2)
  
  # Mean Absolute Percentage Error (MAPE): 
  mape <- mean(per_err)*100
  
  # Retorna los resultadoss: 
  return(data_frame(MAE = mae, MSE = mse, MAPE = mape, N = length(act)))
}

# Data Frame ets: 
df_ets <- data_frame(Actual = salidas.ts.test %>% as.vector(), 
                       Prediccion = fc.ets$mean %>% as.vector() %>% round(0), 
                       Error = Prediccion - Actual,
                       Error_Porcentual = round(Error / Actual, 2))

# Data Frame arima: 
df_arima <- data_frame(Actual = salidas.ts.test %>% as.vector(), 
                       Prediccion = fc.fitarima1$mean %>% as.vector() %>% round(0), 
                       Error = Prediccion - Actual,
                       Error_Porcentual = round(Error / Actual, 2))

# Data Frame prophet: 
df_prophet <- data_frame(Actual = salidas.ts.test %>% as.vector(), 
                       Prediccion = s2$yhat %>% as.vector() %>% round(0), 
                       Error = Prediccion - Actual,
                       Error_Porcentual = round(Error / Actual, 2))

# Comparando los tres modelos: 

bind_rows(df_ets %>% eval_medidas(), 
          df_arima %>% eval_medidas(), 
          df_prophet %>% eval_medidas(),
          df_ensamble %>% eval_medidas) %>%
  mutate(MODELO = c("ETS", "ARIMA", "PROPHET", "ENSAMBLE")) %>% 
  select(MODELO, everything()) %>% 
  mutate_if(is.numeric, function(x) {round(x, 2)}) %>%
  arrange(desc(-MAPE)) %>% 
  kbl(align = "c",
                    digits = 2 , 
                    format.args = list(big.mark = ",")) %>%
  kable_paper("hover", full_width = F)

MODELO	MAE	MSE	MAPE	N
ARIMA	23,194.50	1,452,653,203	1.09	12
ENSAMBLE	28,016.67	1,627,166,964	1.29	12
ETS	32,768.83	2,093,607,703	1.54	12
PROPHET	41,345.83	2,639,896,159	1.89	12

 

Después de hacer el ensamble de los tres mejores modelos a final de cuentas no se pudo mejorar el MAPE del mejor modelo, el arima, lo anterior, no es de sorprender, pues el ensamble se recomienda cuando los modelos que se juntarán presentan diferencias entre ellos, lo que en este caso no ocurrió.

Modelo final

 

Por lo anterior, nos quedamos con el modelo arima, que tuvo un MAPE de test de 1.09%.

Vamos a pronosticar con ese modelo el primer bimestre de 2020, ya que esos meses no están afectados por la pandemia en México, ya que el confinamiento empezó a finales de marzo:

# se cargan los datos
salidas.2020.C <- read.csv("Salidas_Nacionales_C.csv")
salidas.2020.C.ts <- ts(salidas.2020.C[,2], 
                      start = c(2008, 1),
                      end = c(2020, 2),
                      frequency = 12)

ts_info(salidas.2020.C.ts)

##  The salidas.2020.C.ts series is a ts object with 1 variable and 146 observations
##  Frequency: 12 
##  Start time: 2008 1 
##  End time: 2020 2

# Outliers
tsoutliers(salidas.2020.C.ts)

## $index
## [1] 14 17
## 
## $replacements
## [1]  874902.5 1039538.2

salidas.2020.C.ts <- tsclean(salidas.2020.C.ts)


# data real a pronosticar
salidas.2020 <- read.csv("Salidas_Nacionales_2020.csv")

salidas.2020.ts <- ts(salidas.2020[,2], 
                      start = c(2020, 1),
                      end = c(2020, 2),
                      frequency = 12)

ts_info(salidas.2020.ts)

##  The salidas.2020.ts series is a ts object with 1 variable and 2 observations
##  Frequency: 12 
##  Start time: 2020 1 
##  End time: 2020 2

Se crea el modelo:

fitarima.final <- auto.arima(salidas.ts)
fitarima.final

## Series: salidas.ts 
## ARIMA(0,1,1)(2,1,0)[12] 
## 
## Coefficients:
##           ma1     sar1     sar2
##       -0.3451  -0.7327  -0.4491
## s.e.   0.0917   0.0813   0.0844
## 
## sigma^2 estimated as 1.645e+09:  log likelihood=-1578.88
## AIC=3165.75   AICc=3166.07   BIC=3177.25

coeftest(fitarima.final)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1  -0.345105   0.091695 -3.7636 0.0001675 ***
## sar1 -0.732749   0.081267 -9.0165 < 2.2e-16 ***
## sar2 -0.449059   0.084408 -5.3201 1.037e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fitarima.final)

## Series: salidas.ts 
## ARIMA(0,1,1)(2,1,0)[12] 
## 
## Coefficients:
##           ma1     sar1     sar2
##       -0.3451  -0.7327  -0.4491
## s.e.   0.0917   0.0813   0.0844
## 
## sigma^2 estimated as 1.645e+09:  log likelihood=-1578.88
## AIC=3165.75   AICc=3166.07   BIC=3177.25
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 2409.473 38240.05 27829.13 0.1233916 2.058852 0.2475824 0.04054564

Análisis de los residuales:

checkresiduals(fitarima.final)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(2,1,0)[12]
## Q* = 14.857, df = 21, p-value = 0.83
## 
## Model df: 3.   Total lags used: 24

Análisis de la normalidad:

qqnorm(fitarima.final$residuals)
qqline(fitarima.final$residuals)

jarque.bera.test(fitarima.final$residuals)

## 
##  Jarque Bera Test
## 
## data:  fitarima.final$residuals
## X-squared = 4.6355, df = 2, p-value = 0.09849

ggtsdiag(fitarima.final)

Este modelo final no presenta autocorrelación en sus residuales, ya que su gráfica no muestra tendencia y sí pasa el test de Ljung-Boxy, además, éstos se distribuyen normalmente. En su correlograma los coeficientes no rebasan los límites. Por lo anterior, el modelo es apto para pronosticar.

Se hace la predicción de los próximos 2 meses:

fc.fitarima.final <- forecast(fitarima.final, h = 2)
forecast::accuracy(fc.fitarima.final, salidas.2020.ts)

##                    ME     RMSE       MAE       MPE     MAPE       MASE
## Training set 2409.473 38240.05 27829.128 0.1233916 2.058852 0.24758238
## Test set     8773.898 13153.31  9799.399 0.4636181 0.512745 0.08718054
##                     ACF1 Theil's U
## Training set  0.04054564        NA
## Test set     -0.50000000 0.1003116

autoplot(fc.fitarima.final) + 
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Pronóstico Arima Final: Salidas 2020", y = "Salidas")

test_forecast(forecast.obj = fc.fitarima.final, 
              actual = salidas.2020.C.ts,
              test = salidas.2020.ts) %>%
  layout(legend = list(x = 0.1, y = 0.95))

El error MAPE del modelo es: 0.512745. Hay que tener presente que los argumentos de este modelo son diferentes al modelo arima1 creado anteriormente, esto debido a que para su construcción se utilizó toda la serie histórica, desde enero de 2008 hasta diciembre de 2019, por lo tanto, ya no evaluamos el modelo con un conjunto de test, sino que directamente lo evaluamos sobre los datos reales de enero-febrero de 2020 (que son las observaciones que se pronosticaron).

# Data Frame arima final: 
data_frame(Fecha = salidas.2020$Fecha %>% as.vector(),
           Actual = salidas.2020.ts %>% as.vector(), 
           Prediccion = fc.fitarima.final$mean %>% as.vector() %>% round(0), 
           Error = abs(Prediccion - Actual),
           Error_Porcentual = round(Error / Actual, 2)) %>% 
  kbl(caption = "PRONÓSTICO ENERO-FEBRERO 2020", align = "c",
                    digits = 4 , 
                    format.args = list(big.mark = ",")) %>%
  kable_paper("hover", full_width = F)

PRONÓSTICO ENERO-FEBRERO 2020

Fecha	Actual	Prediccion	Error	Error_Porcentual
01/01/2020	2,087,450	2,088,476	1,026	0.00
01/02/2020	1,902,294	1,883,721	18,573	0.01

Conclusiones

La serie presentó cambio de tendencia en el año 2011 debido a diversos factores coyunturales que se presentaron a lo largo del periodo 2008-2011, como la crisis financiera de 2009, la pandemia de Influenza H1N1 y la quiebra de la empresa Mexinana de Aviación.

Se detectaron patrones estacionales en los meses de julio, agosto y diciembre (periodo de alta demanda) y en enero, febrero y septiembre (periodo de baja demanda).

Cuando se vaya a trabajar con series de tiempo es de suma importancia el hacer primero un análisis visual a la serie o series, ya que eso permitirá conocer si hay presencia de tendencia, su tipo y si existe sospecha de algún patrón estacional. Lo anterior, permitirá identificar el tipo de algoritmo a usar, ya que hay algunos que no están diseñados para modelar más de un patrón estacional como los de suavisado exponencial o arima.

Es altamente recomendable hacer un primer modelo que será la base a partir del cual se tratará de mejorarlo de acuerdo con la métrica elegida para su evaluación y, además, no olvidar que la evaluación se debe de hacer en la parte de prueba, es decir, sobre el conjunto de datos que el algoritmo no ha visto.

Cuando se opte por hacer un ensamble de modelos se debe tener en cuenta que esta estratégia funcionará mejor si los modelos a juntar presentan diferencias marcadas entre sí, porque de lo contrario se puede generar un peor resultado.