Sección 10 : Forecasting with Exponential Smoothing Models

En este capítulo, exploraremos más a fondo las funciones de suavizado y su uso en pronósticos. Los modelos que veremos pueden abordar diversas series temporales, con o sin tendencias y componentes estacionales. Iniciaremos con el modelo de media móvil y el suavizado exponencial simple, para luego incrementar la complejidad del modelo y su capacidad para tratar con series temporales más complejas.

The simple moving average

La función de media móvil simple (SMA), usada en el capítulo 5 para suavizar datos de series temporales, puede convertirse en un modelo de pronóstico con algunos ajustes.

Vamos a personalizar la función SMA para pronosticar los precios mensuales del café Robusta durante los próximos 12 meses. Los precios del café Robusta, que no siguen una tendencia o patrón estacional específico, presentan un componente cíclico con ciclos de magnitud y duración variables. Estos datos forman parte del conjunto Coffee_Prices disponible en el paquete TSstudio.

library (TSstudio) 
data(Coffee_Prices) 
ts_info(Coffee_Prices) 
##  The Coffee_Prices series is a mts object with 2 variables and 701 observations
##  Frequency: 12 
##  Start time: 1960 1 
##  End time: 2018 5

Este conjunto de datos es un objeto de tipo mt s y contiene los precios mensuales (USD por kg) tanto del café Robusta como del café Arábica entre 1960 y 2018. Extraeremos los precios mensuales del café Robusta del objeto Coffee_Prices.

robusta <- Coffee_Prices[,1] 

ts_plot (robusta, title = "The Robusta Coffee Monthly Prices", 
Ytitle = "Price in USD", 
Xtitle = "Year")

A continuación, crearemos una función básica de SMA (Media Móvil Simple) utilizando algunas de las funcionalidades del paquete tidyr.

library(tidyr)

sma_forecast <- function(df, h, m, w = NULL){
  
  # Error handling
  if(h > nrow(df)){
    stop("The length of the forecast horizon must be shorter than the length of the series")}
  
  if(m > nrow(df)){
    stop("The length of the rolling window must be shorter than the length of the series")}
  if(!is.null(w)){
    if(length(w) != m){
      stop("The weight argument is not aligned with the length of the rolling window")
    } else if(sum(w) !=1){
      stop("The sum of the average weight is different than 1")
    }
  }
  
  # Setting the average weigths
  if(is.null(w)){
    w <- rep(1/m, m)
  }
  
  # Setting the data frame 
  #-----------------------
  # Changing the Date object column name
  names(df)[1] <- "date" 
  # Setting the training and testing partition 
  # according to the forecast horizon
  df$type <- c(rep("train", nrow(df) - h), 
               rep("test", h)) 
  
  # Spreading the table by the partition type
  df1 <- df %>% spread(key = type, value = y)
  
  # Create the target variable
  df1$yhat <- df1$train
  
  
  # Simple moving average function
  for(i in (nrow(df1) - h + 1):nrow(df1)){
    r <- (i-m):(i-1) 
    df1$yhat[i] <- sum(df1$yhat[r] * w) 
  } 
  
  # dropping from the yhat variable the actual values
  # that were used for the rolling window
  df1$yhat <- ifelse(is.na(df1$test), NA, df1$yhat)
  
  df1$y <- ifelse(is.na(df1$test), df1$train, df1$test)
  
  return(df1)
}

Vamos a utilizar esta función para demostrar el rendimiento de la función SMA. Vamos a predecir los últimos 24 meses de la serie Robusta utilizando una ventana móvil de 3, 6, 12, 24 y 36 meses

robusta_df <- ts_to_prophet(robusta)

robusta_fc_m1 <-  sma_forecast(robusta_df, h = 24, m = 1)
robusta_fc_m6 <-  sma_forecast(robusta_df, h = 24, m = 6)
robusta_fc_m12 <- sma_forecast(robusta_df, h = 24, m = 12)
robusta_fc_m24 <- sma_forecast(robusta_df, h = 24, m = 24)
robusta_fc_m36 <- sma_forecast(robusta_df, h = 24, m = 36) 
library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
plot_ly(data = robusta_df[650:nrow(robusta_df),], x = ~ ds, y = ~ y,
        type = "scatter", mode = "lines", 
        name = "Actual") %>%
  add_lines(x = robusta_fc_m1$date, y = robusta_fc_m1$yhat, 
            name = "SMA - 1", line = list(dash = "dash")) %>%
  add_lines(x = robusta_fc_m6$date, y = robusta_fc_m6$yhat, 
            name = "SMA - 6", line = list(dash = "dash")) %>%
  add_lines(x = robusta_fc_m12$date, y = robusta_fc_m12$yhat, 
            name = "SMA - 12", line = list(dash = "dash")) %>%
  add_lines(x = robusta_fc_m24$date, y = robusta_fc_m24$yhat, 
            name = "SMA - 24", line = list(dash = "dash")) %>%
  add_lines(x = robusta_fc_m36$date, y = robusta_fc_m36$yhat, 
            name = "SMA - 36", line = list(dash = "dash")) %>%
  layout(title = "Forecasting the Robusta Coffee Monthly Prices",
         xaxis = list(title = ""),
         yaxis = list(title = "USD per Kg."))

Las principales observaciones del gráfico anterior son las siguientes: - Si la longitud de la ventana móvil es más corta: - El rango de la predicción está bastante cerca de las observaciones más recientes de la serie. - La predicción converge más rápido a algún valor constante. - Si la longitud de la ventana es más larga: - Toma más tiempo hasta que la predicción converge a algún valor constante. - Puede manejar mejor los shocks y los valores atípicos. - Un modelo de pronóstico SMA con una ventana móvil de una longitud de 1 es equivalente al modelo de pronóstico Naive.

Aunque la función SMA es bastante simple de usar y consume poco poder de cálculo, tiene algunas limitaciones: - El poder de pronóstico de la función SMA está limitado a un horizonte corto y puede tener un rendimiento deficiente a largo plazo. - Este método está limitado para datos de series temporales, sin tendencia o patrones estacionales. Esto afecta principalmente al promedio aritmético que suaviza el patrón estacional y se vuelve plano a largo plazo.

Weighted moving average

El promedio móvil ponderado (WMA) es una versión extendida de la función SMA, y se basa en el uso del promedio ponderado (en lugar del promedio aritmético). La principal ventaja de la función WMA, con respecto a la función SMA, es que te permite distribuir el peso de los retrasos en la ventana móvil. Esto puede ser útil cuando la serie tiene una alta correlación con algunos de sus retrasos

data(USgas)

USgas_df <- ts_to_prophet(USgas)
# -------- Code Chank 8 --------
USgas_fc_m12a <- sma_forecast(USgas_df,
                              h = 24,
                              m = 12,
                              w = c(1, rep(0,11)))
USgas_fc_m12b <- sma_forecast(USgas_df,
                              h = 24,
                              m = 12,
                              w = c(0.8, rep(0,10), 0.2))
plot_ly(data = USgas_df[190:nrow(USgas_df),], x = ~ ds, y = ~ y,
        type = "scatter", mode = "lines", 
        name = "Actual") %>%
  add_lines(x = USgas_fc_m12a$date, y = USgas_fc_m12a$yhat, 
            name = "WMA - Seasonal Lag", line = list(dash = "dash")) %>%
  add_lines(x = USgas_fc_m12b$date, y = USgas_fc_m12b$yhat, 
            name = "WMA - 12 (0.2/0.8)", line = list(dash = "dash")) %>%
  layout(title = "Forecasting the Monthly Consumption of Natural Gas in the US",
         xaxis = list(title = ""),
         yaxis = list(title = "Billion Cubic Feet"))

Como se puede ver en el gráfico anterior, ambos modelos capturaron hasta cierto punto la oscilación estacional de la serie. Establecer todo el peso en el retraso estacional es equivalente al modelo Naive estacional. Esta estrategia podría ser útil para una serie con un patrón estacional dominante, como Usgas. En el segundo ejemplo, ponderamos el promedio entre el retraso más reciente y el retraso estacional. Tendría sentido distribuir los pesos entre los diferentes retrasos cuando la serie tiene una alta correlación con esos retrasos.

Aunque WMA puede capturar el componente estacional de una serie, no puede capturar la tendencia de la serie (debido al efecto promedio). Por lo tanto, este método comenzará a perder su efectividad una vez que el horizonte de pronóstico cruce la longitud de la frecuencia de la serie (por ejemplo, más de un año para series mensuales).

Simple exponential smoothing model

Es el modelo de pronóstico más sencillo de la familia de alisado exponencial. La principal suposición de este modelo es que la serie se mantiene en el mismo nivel (es decir, la media local de la serie es constante) a lo largo del tiempo, y por lo tanto, este modelo es adecuado para series sin componentes de tendencia ni estacionales.

El modelo SES comparte algunos de los atributos del modelo WMA, ya que ambos modelos pronostican los valores futuros de la serie mediante un promedio ponderado de las observaciones pasadas de la serie. Por otro lado, la principal distinción entre los dos es que el modelo SES utiliza todas las observaciones anteriores, mientras que el modelo WMA solo utiliza las m observaciones más recientes (para un modelo con una ventana móvil de una longitud de #7).

El principal atributo del modelo SES es la técnica del promedio ponderado, que se basa en la disminución exponencial de los pesos de observación según su distancia cronológica (es decir, índice de serie o marca de tiempo) desde los primeros valores pronosticados. La tasa de disminución de los pesos de observación

alpha_df <- data.frame(index = seq(from = 1, to = 15, by = 1),
                       power = seq(from = 14, to = 0, by = -1))

alpha_df$alpha_0.01 <- 0.01 * (1 - 0.01) ^ alpha_df$power
alpha_df$alpha_0.2 <- 0.2 * (1 - 0.2) ^ alpha_df$power
alpha_df$alpha_0.4 <- 0.4 * (1 - 0.4) ^ alpha_df$power
alpha_df$alpha_0.6 <- 0.6 * (1 - 0.6) ^ alpha_df$power
alpha_df$alpha_0.8 <- 0.8 * (1 - 0.8) ^ alpha_df$power
alpha_df$alpha_1 <- 1 * (1 - 1) ^ alpha_df$power

plot_ly(data = alpha_df) %>%
  add_lines(x = ~ index, y = ~ alpha_0.01, name = "alpha = 0.01") %>%
  add_lines(x = ~ index, y = ~ alpha_0.2, name = "alpha = 0.2") %>%
  add_lines(x = ~ index, y = ~ alpha_0.4, name = "alpha = 0.4") %>%
  add_lines(x = ~ index, y = ~ alpha_0.6, name = "alpha = 0.6") %>%
  add_lines(x = ~ index, y = ~ alpha_0.8, name = "alpha = 0.8") %>%
  add_lines(x = ~ index, y = ~ alpha_1, name = "alpha = 1") %>%
  layout(title = "Decay Rate of the SES Weights",
         xaxis = list(title = "Index"),
         yaxis = list(title = "Weight"))

Forecasting with the ses function

El paquete de pronóstico proporciona un modelo SES personalizado con la función ses. Los principales argumentos de esta función son los siguientes:

initial: Define el método para inicializar el valor de ¥, que puede calcularse utilizando las primeras observaciones de la serie al establecer el argumento en simple, o estimándolo con el modelo et s (una versión avanzada del modelo Holt-Winters del paquete de pronóstico) al establecerlo en óptimo. alpha: Define el valor del parámetro de suavizado del modelo. Si se establece en NULL, la función lo estimará. h: Establece el horizonte de pronóstico. Vamos a utilizar la función ses para pronosticar nuevamente los precios mensuales del café Robusta. Dejaremos los últimos 12 meses de la serie como un conjunto de pruebas para evaluar el rendimiento del modelo.

robusta_par <- ts_split(robusta, sample.out = 12)
train <- robusta_par$train
test <- robusta_par$test

library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
fc_ses <- ses(train, h = 12, initial = "optimal")

fc_ses$model
## Simple exponential smoothing 
## 
## Call:
##  ses(y = train, h = 12, initial = "optimal") 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 0.6957 
## 
##   sigma:  0.161
## 
##      AIC     AICc      BIC 
## 1989.646 1989.681 2003.252
test_forecast(actual = robusta,
              forecast.obj = fc_ses,
              test = test) %>%
  layout(title = "Robusta Coffee Prices Forecast vs. Actual",
         xaxis = list(range = c(2010, max(time(robusta)))),
         yaxis = list(range = c(1, 3)))

Como podemos ver en el gráfico de pronóstico anterior, la función ses utiliza el conjunto de entrenamiento para identificar el nivel de la serie estimando el parámetro alfa y el nivel inicial del modelo (o Y1). El nivel del valor de pronóstico está bastante cerca del valor de la última observación de la serie, ya que el valor & en este caso, está cerca de 1. Dado que el objetivo del modelo SES es pronosticar el nivel de la serie, el modelo no capturará ninguna oscilación a corto plazo

plot_forecast(fc_ses) %>%
  add_lines(x = time(test) + deltat(test),
            y = as.numeric(test),
            name = "Testing Partition") %>%
  layout(title = "Robusta Coffee Prices Forecast vs. Actual",
         xaxis = list(range = c(2010, max(time(robusta)) +
                                  deltat(robusta))),
         yaxis = list(range = c(0, 4)))

En el caso de un pronóstico plano, los intervalos de confianza del modelo juegan un papel crítico, ya que el nivel de incertidumbre es mayor. Por lo tanto, será útil evaluar si los valores de pronóstico están dentro de los límites del intervalo de confianza del modelo.

Como puedes ver en el gráfico de pronóstico anterior, el conjunto de pruebas está dentro del rango del intervalo de confianza del 80%.